13 #include "TEfficiency.h"
15 #include "TGraphAsymmErrors.h"
22 TChain *all =
new TChain(treename.c_str());
24 for (
int i = 0;
i < filecount; ++
i)
29 temp = name+
string(s.str())+extension;
30 all->Add(temp.c_str());
41 std::vector<TH1F*> plots;
42 ttree->SetBranchAddress(
"photon_m", &photon_m );
43 back->SetBranchAddress(
"photon_m", &back_m );
44 back->SetBranchAddress(
"vtx_radius", &back_Vtx );
45 back->SetBranchAddress(
"track_pT", &back_pt );
47 plots.push_back(
new TH1F(
"m^{#gamma}_{reco}",
"",60,0,.18));
48 plots.push_back(
new TH1F(
"m^{bkgd}_{reco}",
"",60,0,.18));
51 for (
int i = 0;
i < plots.size(); ++
i)
55 for (
int event = 0;
event < ttree->GetEntries(); ++
event)
57 ttree->GetEvent(
event);
58 plots[0]->Fill(photon_m);
61 for (
int event = 0;
event < back->GetEntries(); ++
event)
63 back->GetEvent(
event);
64 if(back_Vtx<21.&&back_pt>2.5)
65 plots[1]->Fill(back_m);
68 cout<<
"Background mass in range count= "<<plots[1]->Integral()<<
'\n';
69 for (
int i = 0;
i < plots.size(); ++
i)
71 plots[
i]->Scale(1./ttree->GetEntries(),
"width");
74 ttree->ResetBranchAddresses();
77 void makeVtxR(TChain* ttree, TTree* vtxTree,TFile* out_file){
82 ttree->SetBranchAddress(
"vtx_radius",&vtxr);
83 ttree->SetBranchAddress(
"tvtx_radius",&tvtxr);
84 vtxTree->SetBranchAddress(
"vtx_x",&vtxX);
85 vtxTree->SetBranchAddress(
"vtx_y",&vtxY);
87 std::vector<TH1F*> plots;
88 plots.push_back(
new TH1F(
"vtx_reco",
"",30,0,40));
89 plots.push_back(
new TH1F(
"vtx_corrected",
"",30,0,40));
90 plots.push_back(
new TH1F(
"vtx_truth",
"",30,0,40));
97 for (
int event = 0;
event < ttree->GetEntries(); ++
event)
99 ttree->GetEvent(
event);
100 plots[1]->Fill(vtxr);
101 plots[2]->Fill(tvtxr);
102 calc+=TMath::Abs(vtxr-tvtxr);
104 calc/=ttree->GetEntries();
105 for (
int i = 0;
i < vtxTree->GetEntries(); ++
i)
107 vtxTree->GetEvent(
i);
108 plots[0]->Fill(sqrt(vtxX*vtxX+vtxY*vtxY));
110 for (
int i = 0;
i < 3; ++
i)
112 plots[
i]->Scale(1./plots[2]->GetEntries(),
"width");
115 std::cout<<
"mean deviation="<<calc<<std::endl;
121 ttree->SetBranchAddress(
"vtx_radius",&r);
122 ttree->SetBranchAddress(
"tvtx_radius",&tr);
123 TH1F *vtxeffPlot =
new TH1F(
"#frac{#Deltar_{vtx}_^{#it{reco}}}{r_{vtx}^{#it{truth}}}",
"",40,-2,2);
124 TH2F *vtxefffuncPlot =
new TH2F(
"vtx_resolution_to_truthvtx",
"",20,0,21,40,-1.5,1.5);
126 vtxefffuncPlot->Sumw2();
127 for (
int event = 0;
event < ttree->GetEntries(); ++
event)
129 ttree->GetEvent(
event);
131 vtxeffPlot->Fill((r-tr)/tr);
132 vtxefffuncPlot->Fill(tr,(r-tr)/tr);
134 vtxeffPlot->Scale(1./ttree->GetEntries(),
"width");
135 vtxefffuncPlot->Scale(1./ttree->GetEntries(),
"width");
137 ttree->ResetBranchAddresses();
144 ttree->SetBranchAddress(
"vtx_radius",&r);
145 ttree->SetBranchAddress(
"tvtx_radius",&tr);
146 ttree->SetBranchAddress(
"track_pT",&pT);
148 TH1F *recoR=
new TH1F(
"vtxrecoR",
"",30,0,40);
149 TH1F *truthR=
new TH1F(
"vtxtruthR",
"",30,0,40);
152 for (
int event = 0;
event < ttree->GetEntries(); ++
event)
154 ttree->GetEvent(
event);
155 if(r>0) recoR->Fill(tr);
158 vtxEff =
new TEfficiency(*recoR,*truthR);
159 vtxEff->SetName(
"vtxEff");
162 ttree->ResetBranchAddresses();
168 ttree->SetBranchAddress(
"photon_pT",&pT);
169 ttree->SetBranchAddress(
"tphoton_pT",&tpT);
170 TH1F *pTRes =
new TH1F(
"#frac{#it{p}^{T}_{reco}}{#it{p}_{#it{truth}}^{T}}",
"",40,0,2.5);
172 for (
int event = 0;
event < ttree->GetEntries(); ++
event)
174 ttree->GetEvent(
event);
175 if(pT>0)pTRes->Fill(pT/tpT);
177 pTRes->Scale(1./pTRes->Integral());
178 ttree->ResetBranchAddresses();
182 TEfficiency*
makepTDist(TChain* ttree,TTree* allTree,TFile* out_file){
186 ttree->SetBranchAddress(
"photon_pT",&pT);
187 ttree->SetBranchAddress(
"tphoton_pT",&tpT);
189 vector<float> *allpT=NULL;
190 allTree->SetBranchAddress(
"photon_pT",&allpT);
192 TH1F *pTeffPlot =
new TH1F(
"#frac{#it{p}^{T}}{#it{p}_{#it{truth}}^{T}}",
"",40,-2,2);
193 TH2F *pTefffuncPlot =
new TH2F(
"pT_resolution_to_truthpt",
"",40,1,35,40,-1.5,1.5);
194 TH1F *converted_pTspec =
new TH1F(
"converted_photon_truth_pT",
"",25,5,30);
195 TH1F *all_pTspec =
new TH1F(
"all_photon_truth_pT",
"",25,5,30);
198 converted_pTspec->Sumw2();
200 pTefffuncPlot->Sumw2();
201 for (
int event = 0;
event < ttree->GetEntries(); ++
event)
203 ttree->GetEvent(
event);
204 if(pT>0)pTeffPlot->Fill(pT/tpT);
205 converted_pTspec->Fill(tpT);
206 if(pT>0)pTefffuncPlot->Fill(tpT,pT/tpT);
209 for (
int event = 0;
event < allTree->GetEntries(); ++
event)
211 allTree->GetEvent(
event);
212 for(
auto i : *allpT){
216 TEfficiency* uni_rate =
new TEfficiency(*converted_pTspec,*all_pTspec);
217 pTeffPlot->Scale(1./ttree->GetEntries(),
"width");
218 pTefffuncPlot->Scale(1./ttree->GetEntries());
219 TProfile* resProfile = pTefffuncPlot->ProfileX(
"func_prof",5,30);
223 ttree->ResetBranchAddresses();
229 signalTree->SetBranchAddress(
"track_deta",&detas);
230 background->SetBranchAddress(
"track_deta",&detab);
231 TH1F *detaS_plot =
new TH1F(
"detaS",
"",3000,-.001,1);
232 TH1F *detaB_plot =
new TH1F(
"detaB",
"",3000,-.001,1);
236 for (
int i = 0;
i < signalTree->GetEntries(); ++
i)
238 signalTree->GetEvent(
i);
239 detaS_plot->Fill(detas);
241 for (
int i = 0;
i < background->GetEntries(); ++
i)
243 background->GetEvent(
i);
244 detaB_plot->Fill(detab);
246 detaB_plot->Scale(1/detaB_plot->Integral());
247 detaS_plot->Scale(1/detaS_plot->Integral());
260 ttree->SetBranchAddress(
"cluster_dphi",&dphi);
261 ttree->SetBranchAddress(
"cluster_prob",&prob);
262 ttree->SetBranchAddress(
"track_layer",&layer);
263 ttree->SetBranchAddress(
"track_dlayer",&dlayer);
264 ttree->SetBranchAddress(
"track_pT",&track_pT);
265 ttree->SetBranchAddress(
"track_deta",&deta);
266 ttree->SetBranchAddress(
"vtx_radius",&radius);
268 TH1F *layerDist =
new TH1F(
"layer",
"",16,-.5,15.5);
269 TH1F *probDist =
new TH1F(
"clust_prob",
"",30,-.5,1.);
270 TH1F *deta_plot =
new TH1F(
"deta",
"",30,-.001,.01);
271 TH1F *dlayer_plot =
new TH1F(
"dlayer",
"",11,-.5,10.5);
272 TH1F *r_plot =
new TH1F(
"signal_vtx_radius_dist",
"",26,-.5,25.5);
276 dlayer_plot->Sumw2();
278 unsigned badLayCount=0;
279 unsigned badClusCount=0;
280 unsigned bigDetaCount=0;
281 unsigned shortRadiusCount=0;
283 for (
int event = 0;
event < ttree->GetEntries(); ++
event)
285 ttree->GetEvent(
event);
286 if(layer==0)badLayCount++;
287 if(dphi<0&&track_pT>.6){
290 if(track_pT>.6&&dphi>=0){
291 layerDist->Fill(layer);
292 probDist->Fill(prob);
293 deta_plot->Fill(deta);
294 dlayer_plot->Fill(TMath::Abs(dlayer));
295 if(deta>.0082||TMath::Abs(dlayer)>9)bigDetaCount++;
297 r_plot->Fill(radius);
298 if(radius<1.84) shortRadiusCount++;
303 layerDist->Scale(1./ttree->GetEntries());
304 deta_plot->Scale(1./ttree->GetEntries());
305 dlayer_plot->Scale(1./ttree->GetEntries());
306 r_plot->Scale(1./ttree->GetEntries());
307 cout<<
"Signal rejection through layer cut= "<<(float)badLayCount/ttree->GetEntries()<<endl;
308 cout<<
"error= "<<sqrt((
float)badLayCount)/ttree->GetEntries()<<endl;
309 cout<<
"Signal rejection through clus cut= "<<(float)badClusCount/ttree->GetEntries()<<endl;
310 cout<<
"error= "<<sqrt((
float)badClusCount)/ttree->GetEntries()<<endl;
311 cout<<
"Signal rejection through deta cut= "<<(float)bigDetaCount/ttree->GetEntries()<<endl;
312 cout<<
"error= "<<sqrt((
float)bigDetaCount)/ttree->GetEntries()<<endl;
313 cout<<
"Signal rejection through radius cut= "<<(float)shortRadiusCount/ttree->GetEntries()<<endl;
314 cout<<
"error= "<<sqrt((
float)shortRadiusCount)/ttree->GetEntries()<<endl;
322 ttree->SetBranchAddress(
"refitdiffx",&diffx);
323 ttree->SetBranchAddress(
"refitdiffy",&diffy);
324 ttree->SetBranchAddress(
"refitdiffz",&diffz);
325 TH1F *diffplotx=
new TH1F(
"x_{0}-x_{refit}",
"",40,-10,10);
326 TH1F *diffploty=
new TH1F(
"y_{0}-y_{refit}",
"",40,-10,10);
327 TH1F *diffplotz=
new TH1F(
"z_{0}-z_{refit}",
"",40,-10,10);
333 for (
int event = 0;
event < ttree->GetEntries(); ++
event)
335 ttree->GetEvent(
event);
336 diffplotx->Fill(diffx);
337 diffploty->Fill(diffy);
338 diffplotz->Fill(diffz);
341 diffplotx->Scale(1./ttree->GetEntries(),
"width");
342 diffploty->Scale(1./ttree->GetEntries(),
"width");
343 diffplotz->Scale(1./ttree->GetEntries(),
"width");
346 ttree->ResetBranchAddresses();
351 caloFile.open(filename.c_str());
354 vector<double> xData, yData;
359 while(caloFile >>x>>s>>y){
363 double *xArray, *yArray;
366 TGraph *pTResCaloGraph =
new TGraph(xData.size(),xArray,yArray);
367 pTResCaloGraph->SetNameTitle(
"calopTRes",
"calopTRes");
368 pTResCaloGraph->Sort();
369 pTResCaloGraph->Write();
374 vector<float>* tpT= NULL;
375 ttree->SetBranchAddress(
"photon_pT",&tpT);
379 title+=
"_photon_truth_pT";
381 else title=
"photon_truth_pT";
382 TH1F *tpTspec =
new TH1F(title.c_str(),
"",20,5,25);
384 cout<<
"pythia tree with: "<<ttree->GetEntries()<<
" entries"<<endl;
385 for (
int event = 0;
event < ttree->GetEntries(); ++
event)
387 ttree->GetEvent(
event);
389 if(
i>5)tpTspec->Fill(
i);
393 ttree->ResetBranchAddresses();
399 string rname = h.GetName();
401 TH1F
r = TH1F(rname.c_str(),h.GetTitle(),h.GetNbinsX()-1,h.GetBinLowEdge(2),h.GetBinLowEdge(h.GetNbinsX()+1));
403 for(
unsigned i=1;
i<r.GetNbinsX();++
i){
404 r.SetBinContent(
i,h.GetBinContent(
i+1));
405 r.SetBinError(
i,h.GetBinError(
i+1));
415 cout<<
"AD:"<<ad<<endl;
417 if(ad<.2) r= 1-TMath::Exp(-13.436+101.14*ad-223.74*ad*ad);
418 else if (ad<.34) r= 1-TMath::Exp(-8.31+42.79*ad-59.938*ad*ad);
419 else if (ad<.6) r= 1-TMath::Exp(.917-4.279*ad-1.38*ad*ad);
420 else r= 1-TMath::Exp(1.29-5.709*ad+.0186*ad*ad);
426 hard.Scale(1/hard.Integral());
427 soft.Scale(1/soft.Integral());
428 return soft.Chi2Test(&hard,
"WW CHI2/NDF");
432 unsigned bins = soft.GetNbinsX();
433 while(
ADtoP(hard.AndersonDarlingTest(&soft))>.05&&bins>2){
434 cout<<
"Hard/Soft p="<<
ADtoP(hard.AndersonDarlingTest(&soft))<<
" n bins="<<hard.GetNbinsX()<<endl;
435 cout<<
"chi2 p="<<
smartChi2(hard,soft)<<
" n bins="<<hard.GetNbinsX()<<endl;
440 TH1F* temp = (TH1F*) hard.Clone(
"hard_chopped");
445 TH1F* conversion_rate = (TH1F*) pythia->Clone(
"rate");
446 TH1* uni_rate = (TH1F*)rate->GetPassedHistogram()->Clone(
"uni_rate");
447 uni_rate->Divide(rate->GetTotalHistogram());
448 conversion_rate->Multiply(uni_rate);
449 conversion_rate->Scale(1/pythia->Integral());
454 double hardTotal=hard->Integral(matchBin,hard->GetNbinsX());
455 double softTotal=soft->Integral(matchBin,soft->GetNbinsX());
456 double r= softTotal/hardTotal;
457 *error = (softTotal+TMath::Power(softTotal,1./2))/hardTotal;
458 cout<<soft->Integral(matchBin,soft->GetNbinsX())<<endl;
459 cout<<hard->Integral(matchBin,soft->GetNbinsX())<<endl;
464 TH1F* combined =
new TH1F(
"combinedpythia",
"",soft->GetNbinsX(),soft->GetBinLowEdge(1),hard->GetBinLowEdge(hard->GetNbinsX()+1));
465 unsigned matchBin = 10;
466 for (
int i = 1;
i < matchBin; ++
i)
468 combined->SetBinContent(
i,soft->GetBinContent(
i));
469 combined->SetBinError(
i,soft->GetBinError(
i));
473 for (
int i = matchBin;
i < combined->GetNbinsX()+1; ++
i)
475 combined->SetBinContent(
i,hard->GetBinContent(
i));
476 combined->SetBinError(
i,sqrt(hard->GetBinError(
i)*hard->GetBinError(
i)+systematic*systematic));
484 TFile *out_file =
new TFile(
"effplots.root",
"UPDATE");
486 string treeExtension =
".root";
487 unsigned int nFiles=1000;
488 string softPath =
"/sphenix/user/vassalli/minBiasPythia/softana.root";
490 string hard4Path =
"/sphenix/user/vassalli/minBiasPythia/hard4ana.root";
491 TChain *softTree =
new TChain(
"photonTree");
493 TChain *hard4Tree =
new TChain(
"photonTree");
494 softTree->Add(softPath.c_str());
496 hard4Tree->Add(hard4Path.c_str());
506 pythiaSpec =
addSpec(pythiaSpec,hardSpec,out_file);