57 gROOT->LoadMacro(
"$VMCWORKDIR/gconfig/rootlogon.C");
61 gROOT->SetStyle(
"Plain");
62 gStyle->SetOptStat(10);
68 const Int_t nBins = 90;
73 TString inTreeFile =
"out_test.root";
74 TFile *
inFile = TFile::Open(inTreeFile,
"READ");
75 TTree *
tree=(TTree *) inFile->Get(
"nt") ;
78 Int_t
nEvents = nt->GetEntriesFast();
80 Int_t skipcount = 0, moreThanTwentyPercentCount = 0, goodEvents = 0;
86 TString histotitle =
"p_{z} for ";
88 histotitle +=
" #mu^{-} with 0.5 GeV/c < p < 5 GeV/c ";
89 TString histotitleinvpz = histotitle +
"(from 1/p_{z} HS)";
90 TString histotitleinvpzWithB = histotitle +
"(from 1/p_{z} HS with mean B field)";
91 TString histotitlepz = histotitle +
"(from p_{z} HS)";
92 TString histotitlepzWithB = histotitle +
"(from p_{z} HS with mean B field)";
98 TProfile *hpzFromInvpz =
new TProfile(
"pzFrominvpz",histotitleinvpz, nBins, 0.5, 5, 0, 8,
"s");
100 hpzFromInvpz->GetYaxis()->SetTitle(
"reconstructed p_{z} [GeV/c]");
101 hpzFromInvpz->GetXaxis()->SetTitle(
"p_{z} from MC Truth [GeV/c]");
103 TProfile* hpzFromInvpzWithBfield =
new TProfile(
"pzFrominvpzWithB", histotitleinvpzWithB, nBins, 0.5, 5, 0, 8,
"s");
104 hpzFromInvpzWithBfield->GetYaxis()->SetTitle(
"reconstructed p_{z} [GeV/c]");
105 hpzFromInvpzWithBfield->GetXaxis()->SetTitle(
"p_{z} from MC Truth [GeV/c]");
107 TProfile* hpzFrompz =
new TProfile(
"pzFrompz", histotitlepz, nBins, 0.5, 5, 0, 8,
"s");
108 hpzFrompz->GetYaxis()->SetTitle(
"reconstructed p_{z} [GeV/c]");
109 hpzFrompz->GetXaxis()->SetTitle(
"p_{z} from MC Truth [GeV/c]");
111 TProfile* hpzFrompzWithBfield =
new TProfile(
"pzFrompzWithB", histotitlepzWithB, nBins, 0.5, 5, 0, 8,
"s");
112 hpzFrompzWithBfield->GetYaxis()->SetTitle(
"reconstructed p_{z} [GeV/c]");
113 hpzFrompzWithBfield->GetXaxis()->SetTitle(
"p_{z} from MC Truth [GeV/c]");
118 TH1F* hhitshiftinx =
new TH1F(
"hitshiftinx",
"hitshiftinx for single track #mu^{-} p = 3 GeV/c", 200, -50, 50);
119 hhitshiftinx->GetXaxis()->SetTitle(
"hitshiftinx [cm]");
120 hhitshiftinx->GetYaxis()->SetTitle(
"Counts");
165 Float_t mc_px = 0., mc_py = 0., mc_pz = 0., mc_theta = 0., mc_theta_xz = 0., mc_phi = 0.;
167 Float_t nHitsForHoughSpaceLine = 0, nHitsForHoughSpaceParabola = 0;
168 Float_t pzFromInvpz = 0., pzFrompz = 0., pzFromInvpzWithBField = 0., pzFrompzWithBField = 0., hitshiftinx =0.;
169 Float_t thetaPeakLine = 0., thetaPeakParabolaFrominvpz = 0., thetaPeakParabolaFrompz = 0.;
170 Float_t mc_theta = 0., mc_theta_xz = 0.;
172 nt->SetBranchAddress(
"mc_pz", &mc_pz);
173 nt->SetBranchAddress(
"mc_theta", &mc_theta);
174 nt->SetBranchAddress(
"mc_theta_xz", &mc_theta_xz);
175 nt->SetBranchAddress(
"pzFromInvpz", &pzFromInvpz);
176 nt->SetBranchAddress(
"pzFromInvpzWithBField", &pzFromInvpzWithBField);
177 nt->SetBranchAddress(
"pzFrompz", &pzFrompz);
178 nt->SetBranchAddress(
"pzFrompzWithBField", &pzFrompzWithBField);
179 nt->SetBranchAddress(
"hitshiftinx", &hitshiftinx);
180 nt->SetBranchAddress(
"nHitsForHoughSpaceLine",&nHitsForHoughSpaceLine);
181 nt->SetBranchAddress(
"nHitsForHoughSpaceParabola",&nHitsForHoughSpaceParabola);
182 nt->SetBranchAddress(
"thetaPeakLine", &thetaPeakLine);
183 nt->SetBranchAddress(
"thetaPeakParabolaFrompz", &thetaPeakParabolaFrompz);
184 nt->SetBranchAddress(
"thetaPeakParabolaFrominvpz", &thetaPeakParabolaFrominvpz);
187 const Int_t maxMomResPics = 90;
188 const Int_t picsPerCanvas = 12;
189 const Int_t maxCanvas = maxMomResPics/picsPerCanvas+1;
191 Int_t ResolutionX = 800, ResolutionY = 600;
192 TCanvas *cresinvpz[maxCanvas];
199 const Int_t bins = 200;
200 TH1F *hresolution1invpz[maxMomResPics];
201 TF1 *func1[maxMomResPics];
207 for (Int_t iMomRes=0; iMomRes<maxMomResPics; ++iMomRes)
209 const Double_t momSpacingForResPlot = 0.05;
213 hresolution1invpz[iMomRes] =
new TH1F(
"Relative resolution for muons with low*iMomRes < p_{z} < high GeV/c",
"Relative resolution for muons with low*iMomRes < p_{z} < high GeV/c", bins, -1.0, 1.0);
214 hresolution1invpz[iMomRes]->GetYaxis()->SetTitle(
"counts");
215 hresolution1invpz[iMomRes]->GetXaxis()->SetTitle(
"(reconstructed p_{z} - p_{z} from MC Truth)/ p_{z} from MC Truth");
217 Int_t iCanvas = iMomRes/picsPerCanvas;
218 if (0==iMomRes%picsPerCanvas){
219 cresinvpz[iCanvas]=
new TCanvas(
"cresinvpz"+iCanvas,
"cresinvpz",ResolutionX,ResolutionY);
220 cresinvpz[iCanvas]->Divide(4,3);
225 for (Int_t iEvent=0; iEvent<
nEvents; ++iEvent)
227 nt->GetEntry(iEvent);
228 if (0 == iEvent%1000) { cout <<
"processing event " << iEvent <<
"\n"; }
230 if (2 >= nHitsForHoughSpaceLine) {
continue; }
231 if (4 >= nHitsForHoughSpaceParabola) {
continue; }
236 pzFromInvpz = pzFromInvpz*
cos(thetaPeakParabolaFrominvpz/180.*3.14159265);
241 pzFromInvpz = (pzFromInvpz-p0invpz)/p1invpz;
254 if (maxMomResPics==iMomRes+1){
256 hpzFromInvpz->Fill(mc_pz,pzFromInvpz);
257 hpzFromInvpzWithBfield->Fill(mc_pz,pzFromInvpzWithBField);
258 hpzFrompz->Fill(mc_pz,pzFrompz);
259 hpzFrompzWithBfield->Fill(mc_pz,pzFrompzWithBField);
260 hhitshiftinx->Fill(hitshiftinx);
265 Double_t relativeDeviation_pz = (pzFromInvpz-mc_pz)/mc_pz;
267 if (maxMomResPics==iMomRes+1){
269 if (0.2 <= relativeDeviation_pz)
271 ++moreThanTwentyPercentCount;
272 std::cout <<
"Event " << iEvent <<
" MC p_z = " << mc_pz <<
" FTS p_z = " << pzFromInvpz <<
" is off by " << relativeDeviation_pz*100 <<
" per cent." << std::endl;
274 totalDeviation +=
fabs(relativeDeviation_pz);
283 if ((low < mc_pz) && (mc_pz <= high))
287 hresolution1invpz[iMomRes]->Fill(relativeDeviation_pz);
304 cresinvpz[iCanvas]->cd(iMomRes-iCanvas*picsPerCanvas+1);
305 func1[iMomRes] =
new TF1(
"func1",
"gaus",-0.4,0.4);
306 TF1 *funcptr = func1[iMomRes];
307 int maxbin = hresolution1invpz[iMomRes]->GetMaximumBin();
308 funcptr->SetParameters(hresolution1invpz[iMomRes]->GetBinContent(maxbin), hresolution1invpz[iMomRes]->GetBinCenter(maxbin), 0.1);
312 hresolution1invpz[iMomRes]->Fit(funcptr,
"R",
"",-0.4,0.4);
313 hresolution1invpz[iMomRes]->Draw(
"ep");
317 Double_t sigmaErr = funcptr->GetParError(2);
320 xVal[iMomRes]=(low+
high)/2.0;
322 xErr[iMomRes]=momSpacingForResPlot/2.0;
323 yErr[iMomRes]=sigmaErr;
325 cout <<
"sigma = " << sigma << endl;
326 cout <<
"sigmaerror = " << sigmaErr << endl;
329 cout <<
"End of momentum loop " << iMomRes <<
" with events from "<< low <<
" to " << high <<
" GeV/c"<< endl;
337 TCanvas *cMomRes =
new TCanvas(
"c1",
"A Simple Graph with error bars",800,600);
339 TGraphErrors *momResPlot =
new TGraphErrors(maxMomResPics,xVal,yVal,xErr,yErr);
340 momResPlot->SetTitle(
"Momentum Resolution");
341 momResPlot->SetMarkerColor(4);
342 momResPlot->SetMarkerStyle(21);
343 momResPlot->Draw(
"AP");
345 momResPlot->GetYaxis()->SetTitle(
"relative p_{z} resolution");
346 momResPlot->SetMinimum(0.0);
347 momResPlot->SetMaximum(0.09);
348 momResPlot->GetXaxis()->SetTitle(
"p_{z, MC} [GeV/c]");
353 Int_t ResolutionX = 400, ResolutionY = 1200;
355 TCanvas *cinvpz =
new TCanvas(
"cpzMCvsHoughFromInvpz",
"cpzMCvsHoughFromInvpz", ResolutionX, ResolutionY);
356 TF1 *p1fit =
new TF1(
"P1",
"pol1",0.5,4.95);
357 hpzFromInvpz->Fit(p1fit,
"RS");
358 hpzFromInvpz->Draw();
373 Int_t ResolutionX = 800, ResolutionY = 600;
375 TCanvas *chitshiftinx=
new TCanvas(
"chitshiftinx",
"chitshiftinx", ResolutionX, ResolutionY);
376 hhitshiftinx->Draw();
383 std::cout << moreThanTwentyPercentCount <<
" from " << goodEvents <<
" events deviate more than 20 per cent from MC truth value! Mean absolute relative deviation for pz is " << totalDeviation/goodEvents*100 <<
" per cent." << std::endl;
friend F32vec4 cos(const F32vec4 &a)
friend F32vec4 fabs(const F32vec4 &a)