{ // === CHOOSE ========================== int ifile = 4; double mean_inf = 0.39; double mean_sup = 1.23; double sigma_inf = 0.39; double sigma_sup = 1.23; // y upper limit Double_t yup = 30.; // steps for gaussians (e canvas) double nsteps = 50; int cx = 5; int cy = 10; double limit = 100/nsteps; double flimit = floor(100/nsteps); if( limit != flimit){ nsteps = 100/flimit; limit = flimit; cout << "changing nsteps to " << nsteps << endl; } TString fundedx = "[0] * (1./x)**2 + [1] * TMath::Log(x) + [2]"; // TString fundedx = "pol6(0)"; int npardedx = 3; // TString funsig = "pol4(0)"; // TString funsig_forsum = "pol4(3)"; // int nparsig = 5; TString funsig = "[0] * (1./x)**2 + [1] * TMath::Log(x) + [2]"; TString funsig_forsum = "([3] * (1./x)**2 + [4] * TMath::Log(x) + [5])"; int nparsig = 3; int npartot = npardedx + nparsig; // ==================================== TCanvas *c = new TCanvas("c", "dedx_p"); TFile infile("dedx_out.root", "READ"); TH2F *h; switch(ifile) { case 0: h = (TH2F*) infile.Get("hdedx_p_e"); break; case 1: h = (TH2F*) infile.Get("hdedx_p_mu"); break; case 2: h = (TH2F*) infile.Get("hdedx_p_pi"); break; case 3: h = (TH2F*) infile.Get("hdedx_p_k"); break; case 4: h = (TH2F*) infile.Get("hdedx_p_p"); } c->cd(); h->Draw(); // procedure ------------------------------------------- TH1D *htemp[nsteps]; TF1 *fgaus[nsteps]; Double_t x[nsteps]; // array of mean/sigma Double_t mean[nsteps], sigma[nsteps], meanerr[nsteps], sigmaerr[nsteps]; // gauss TCanvas *c2 = new TCanvas("c2", "slices"); c2->Divide(cx, cy); for(int i = 0; i < nsteps; i++) { double stepsize = 1.5/nsteps; x[i] = (i + 0.5) * stepsize; TString name = "htemp"; name += i; TString title = ""; TString fname = name; TString fgname = name; switch(ifile) { case 0: fname.Prepend("e_"); fgname.Prepend("fe_"); break; case 1: fname.Prepend("mu_"); fgname.Prepend("fmu_"); break; case 2: fname.Prepend("pi_"); fgname.Prepend("fpi_"); break; case 3: fname.Prepend("k_"); fgname.Prepend("fk_"); break; case 4: fname.Prepend("p_"); fgname.Prepend("fp_"); } title = fname; title += "("; title += (x[i] - 0.5 * stepsize); title += ", "; title += (x[i] + 0.5 * stepsize); title += ") GeV"; title += ((i + 1) * limit - 1); htemp[i] = new TH1D(fname, title, 100, 0, yup); h->ProjectionY(fname, limit * i, (i + 1) * limit - 1,"o"); c2->cd(i + 1); htemp[i].Draw(); double tmpmu = htemp[i]->GetMean() ; double tmpsig = htemp[i]->GetRMS(); double gaus_inf = tmpmu - 3 * tmpsig; double gaus_sup = tmpmu + 3 * tmpsig; fgaus[i] = new TF1("fgname", "gaus(0)", gaus_inf, gaus_sup); fgaus[i]->SetParameters(1, tmpmu, tmpsig); htemp[i]->Fit("fgname", "R"); mean[i] = fgaus[i]->GetParameter(1); sigma[i] = fgaus[i]->GetParameter(2); meanerr[i] = fgaus[i]->GetParError(1); sigmaerr[i] = fgaus[i]->GetParError(2); } // GRAPHS ======================================================== TCanvas *graphs = new TCanvas("graphs", "graphs"); graphs->Divide(2, 1); Double_t masses[5] = {0.000510999, 0.1056584, 0.139570, 0.49368, 0.9382723 }; Double_t xerr[nsteps]; for(int istep = 0; istep < nsteps; istep++) xerr[istep] = 0; // ================= MEAN ================= TGraphErrors *g = new TGraphErrors(nsteps, x, mean, xerr, meanerr); TF1 *fdedx = new TF1("fdedx", fundedx, mean_inf, mean_sup); fdedx->SetParameters(1, 1, 1); graphs->cd(1); g->Draw("A*P"); g->GetYaxis()->SetRangeUser(0, 100); g->Fit("fdedx", "R"); Double_t par[npardedx]; fdedx->GetParameters(&par[0]); // ========================================= // ================= SIGMA ================= for(int i = 0; i < nsteps; i++) sigma[i] = fabs(sigma[i]); TGraphErrors *gs = new TGraphErrors(nsteps, x, sigma, xerr, sigmaerr); TF1 *sfdedx = new TF1("sfdedx", funsig, sigma_inf, sigma_sup); graphs->cd(2); gs->Draw("A*P"); gs->GetYaxis()->SetRangeUser(0, 5); gs->Fit("sfdedx", "R"); Double_t spar[nparsig]; sfdedx->GetParameters(&spar[0]); // ========================================= // ==== DRAW ALL ON SCATTER PLOT =========== TF1 *fdedx2 = new TF1("fdedx2", fundedx, 0., 1.5); fdedx2->SetParameters(par); c->cd(); fdedx2->SetLineColor(kOrange); fdedx2->Draw("SAME"); double parcomp[npartot]; for(int i = 0; i < npartot; i++) { if(i < npardedx) parcomp[i] = par[i]; else parcomp[i] = spar[i - npardedx]; } TString funsum = fundedx + " + " + funsig_forsum; TF1 *fsum = new TF1("fsum", funsum, 0., 1.5); fsum->SetParameters(parcomp); TString fundiff = fundedx + " - " + funsig_forsum; TF1 *fdiff = new TF1("fdiff", fundiff, 0., 1.5); fdiff->SetParameters(parcomp); c->cd(); fsum->Draw("SAME"); fdiff->Draw("SAME"); // ========================================= // ==== PRINT OUTS ========================= cout << "ifile " << ifile << endl; cout << " momentum sampled over " << nsteps << " with step width " << 1.5/nsteps << endl; cout << endl; cout << "MEAN DEDX PARAMETRIZATION" << endl; cout << "mom limits " << mean_inf << " " << mean_sup << endl; cout << "mu: "; for(int param = 0; param < npardedx; param++) cout << fdedx->GetParameter(param) << ", "; cout << endl; cout << endl; cout << "SIGMA DEDX PARAMETRIZATION" << endl; cout << "mom limits " << sigma_inf << " " << sigma_sup << endl; cout << "sig: "; for(int param = 0; param < nparsig; param++) cout << sfdedx->GetParameter(param) << ", " ; cout << endl; // ========================================= }