Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
photonEff2.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file photonEff2.C
1 #include <iostream>
2 #include <fstream>
3 using namespace std;
4 
5 #include "TFile.h"
6 #include "Scalar.h"
7 #include "TTree.h"
8 #include "TChain.h"
9 #include "TLegend.h"
10 #include "math.h"
11 #include "TH1.h"
12 #include "TH2.h"
13 #include "TEfficiency.h"
14 #include "TLine.h"
15 #include "TGraphAsymmErrors.h"
16 
17 namespace {
18  unsigned replot=0;
19 }
20 
21 TChain* handleFile(string name, string extension, string treename, unsigned int filecount){
22  TChain *all = new TChain(treename.c_str());
23  string temp;
24  for (int i = 0; i < filecount; ++i)
25  {
26 
27  ostringstream s;
28  s<<i;
29  temp = name+string(s.str())+extension;
30  all->Add(temp.c_str());
31  }
32  return all;
33 }
34 
35 void makephotonM(TChain* ttree,TChain* back, TFile* out_file){
36  float photon_m;
37  float tphoton_m;
38  float back_Vtx;
39  float back_pt;
40  float back_m;
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 );
46  //ttree->SetBranchAddress("rephoton_m", &rephoton_m );
47  plots.push_back(new TH1F("m^{#gamma}_{reco}","",60,0,.18));
48  plots.push_back(new TH1F("m^{bkgd}_{reco}_wide","",60,0,1));
49  //plots.push_back(new TH1F("m^{#gamma}_{recoRefit}","",40,-2,10));
50 
51  for (int i = 0; i < plots.size(); ++i)
52  {
53  plots[i]->Sumw2();
54  }
55  for (int event = 0; event < ttree->GetEntries(); ++event)
56  {
57  ttree->GetEvent(event);
58  plots[0]->Fill(photon_m);
59  // plots[1]->Fill(rephoton_m);
60  }
61  for (int event = 0; event < back->GetEntries(); ++event)
62  {
63  back->GetEvent(event);
64  if(back_Vtx<21.&&back_pt>2.5)
65  plots[1]->Fill(back_m);
66  // plots[1]->Fill(rephoton_m);
67  }
68  // cout<<"Background mass in range count= "<<plots[1]->Integral()<<'\n';
69  for (int i = 0; i < plots.size(); ++i)
70  {
71  plots[i]->Scale(1./ttree->GetEntries(),"width");
72  }
73  out_file->Write();
74  ttree->ResetBranchAddresses();
75 }
76 
77 void makeVtxR(TChain* ttree, TTree* vtxTree,TFile* out_file){
78  float vtxr;
79  float vtxX;
80  float vtxY;
81  float tvtxr;
82  ttree->SetBranchAddress("vtx_radius",&vtxr);
83  ttree->SetBranchAddress("tvtx_radius",&tvtxr);
84  vtxTree->SetBranchAddress("vtx_x",&vtxX);
85  vtxTree->SetBranchAddress("vtx_y",&vtxY);
86 
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));
91 
92  plots[0]->Sumw2();
93  plots[1]->Sumw2();
94  plots[2]->Sumw2();
95 
96  double calc=0;
97  for (int event = 0; event < ttree->GetEntries(); ++event)
98  {
99  ttree->GetEvent(event);
100  plots[1]->Fill(vtxr);
101  plots[2]->Fill(tvtxr);
102  calc+=TMath::Abs(vtxr-tvtxr);
103  }
104  calc/=ttree->GetEntries();
105  for (int i = 0; i < vtxTree->GetEntries(); ++i)
106  {
107  vtxTree->GetEvent(i);
108  plots[0]->Fill(sqrt(vtxX*vtxX+vtxY*vtxY));
109  }
110  for (int i = 0; i < 3; ++i)
111  {
112  plots[i]->Scale(1./plots[2]->GetEntries(),"width");
113  }
114  out_file->Write();
115  std::cout<<"mean deviation="<<calc<<std::endl;
116 }
117 
118 void makeVtxRes(TChain* ttree,TFile* out_file){
119  float r;
120  float tr;
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);
125  vtxeffPlot->Sumw2();
126  vtxefffuncPlot->Sumw2();
127  for (int event = 0; event < ttree->GetEntries(); ++event)
128  {
129  ttree->GetEvent(event);
130  if(r<0) continue;
131  vtxeffPlot->Fill((r-tr)/tr);
132  vtxefffuncPlot->Fill(tr,(r-tr)/tr);
133  }
134  vtxeffPlot->Scale(1./ttree->GetEntries(),"width");
135  vtxefffuncPlot->Scale(1./ttree->GetEntries(),"width");
136  out_file->Write();
137  ttree->ResetBranchAddresses();
138 }
139 
140 void makeVtxEff(TChain* ttree,TFile* out_file){
141  float r;
142  float tr;
143  float pT;
144  ttree->SetBranchAddress("vtx_radius",&r);
145  ttree->SetBranchAddress("tvtx_radius",&tr);
146  ttree->SetBranchAddress("track_pT",&pT);
147  TEfficiency *vtxEff;
148  TH1F *recoR= new TH1F("vtxrecoR","",30,0,40);
149  TH1F *truthR= new TH1F("vtxtruthR","",30,0,40);
150  recoR->Sumw2();
151  truthR->Sumw2();
152  for (int event = 0; event < ttree->GetEntries(); ++event)
153  {
154  ttree->GetEvent(event);
155  if(r>0) recoR->Fill(tr);
156  truthR->Fill(tr);
157  }
158  vtxEff = new TEfficiency(*recoR,*truthR);
159  vtxEff->SetName("vtxEff");
160  vtxEff->Write();
161  out_file->Write();
162  ttree->ResetBranchAddresses();
163 }
164 
165 void makepTRes(TChain* ttree){
166  float pT;
167  float tpT;
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);
171  pTRes->Sumw2();
172  for (int event = 0; event < ttree->GetEntries(); ++event)
173  {
174  ttree->GetEvent(event);
175  if(pT>0)pTRes->Fill(pT/tpT);
176  }
177  pTRes->Scale(1./pTRes->Integral());
178  ttree->ResetBranchAddresses();
179  pTRes->Write();
180 }
181 
182 TEfficiency* makepTDist(TChain* ttree,TTree* allTree,TFile* out_file){
183  float pT;
184  float tpT;
185  float track_pT;
186  ttree->SetBranchAddress("photon_pT",&pT);
187  ttree->SetBranchAddress("tphoton_pT",&tpT);
188 
189  vector<float> *allpT=NULL;
190  allTree->SetBranchAddress("photon_pT",&allpT);
191 
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);
196  //TH1F *trackpTDist = new TH1F("truthpt","",40,0,35);
197  pTeffPlot->Sumw2();
198  converted_pTspec->Sumw2();
199  all_pTspec->Sumw2();
200  pTefffuncPlot->Sumw2();
201  for (int event = 0; event < ttree->GetEntries(); ++event)
202  {
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);
207  //trackpTDist->Fill(track_pT);
208  }
209  for (int event = 0; event < allTree->GetEntries(); ++event)
210  {
211  allTree->GetEvent(event);
212  for(auto i : *allpT){
213  all_pTspec->Fill(i);
214  }
215  }
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);
220  resProfile->Write();
221  //trackpTDist->Scale(1./ttree->GetEntries(),"width");
222  out_file->Write();
223  ttree->ResetBranchAddresses();
224  return uni_rate;
225 }
226 
227 void compareDeta(TTree* signalTree, TTree* background){
228  float detas,detab;
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);
233  detaS_plot->Sumw2();
234  detaB_plot->Sumw2();
235 
236  for (int i = 0; i < signalTree->GetEntries(); ++i)
237  {
238  signalTree->GetEvent(i);
239  detaS_plot->Fill(detas);
240  }
241  cout<<"background has "<<background->GetEntries()<<" entries"<<endl;
242  for (int i = 0; i < background->GetEntries(); ++i)
243  {
244  background->GetEvent(i);
245  detaB_plot->Fill(detab);
246  }
247  detaB_plot->Scale(1/detaB_plot->Integral());
248  detaS_plot->Scale(1/detaS_plot->Integral());
249  detaS_plot->Write();
250  detaB_plot->Write();
251  signalTree->ResetBranchAddresses();
252  background->ResetBranchAddresses();
253 }
254 
255 void testCuts(TChain* ttree,TFile* out_file){
256  float dphi;
257  float prob;
258  float track_pT;
259  float deta;
260  float radius;
261  int layer;
262  int dlayer;
263  ttree->SetBranchAddress("cluster_dphi",&dphi);
264  ttree->SetBranchAddress("cluster_prob",&prob);
265  ttree->SetBranchAddress("track_layer",&layer);
266  ttree->SetBranchAddress("track_dlayer",&dlayer);
267  ttree->SetBranchAddress("track_pT",&track_pT);
268  ttree->SetBranchAddress("track_deta",&deta);
269  ttree->SetBranchAddress("vtx_radius",&radius);
270 
271  TH1F *layerDist = new TH1F("layer","",16,-.5,15.5);
272  TH1F *probDist = new TH1F("clust_prob","",30,-.5,1.);
273  TH1F *deta_plot = new TH1F("deta","",30,-.001,.01);
274  TH1F *dlayer_plot = new TH1F("dlayer","",11,-.5,10.5);
275  TH1F *r_plot = new TH1F("signal_vtx_radius_dist","",26,-.5,25.5);
276  layerDist->Sumw2();
277  probDist->Sumw2();
278  deta_plot->Sumw2();
279  dlayer_plot->Sumw2();
280  r_plot->Sumw2();
281  unsigned badLayCount=0;
282  unsigned badClusCount=0;
283  unsigned bigDetaCount=0;
284  unsigned shortRadiusCount=0;
285 
286  for (int event = 0; event < ttree->GetEntries(); ++event)
287  {
288  ttree->GetEvent(event);
289  if(layer==0)badLayCount++;
290  if(dphi<0&&track_pT>.6){
291  badClusCount++;
292  }
293  if(track_pT>.6&&dphi>=0){
294  layerDist->Fill(layer);
295  probDist->Fill(prob);
296  deta_plot->Fill(deta);
297  dlayer_plot->Fill(TMath::Abs(dlayer));
298  if(deta>.0082||TMath::Abs(dlayer)>9)bigDetaCount++;
299  else{
300  r_plot->Fill(radius);
301  if(radius<1.84) shortRadiusCount++;
302  }
303  }
304 
305  }
306  layerDist->Scale(1./ttree->GetEntries());
307  deta_plot->Scale(1./ttree->GetEntries());
308  dlayer_plot->Scale(1./ttree->GetEntries());
309  r_plot->Scale(1./ttree->GetEntries());
310  cout<<"Signal rejection through layer cut= "<<(float)badLayCount/ttree->GetEntries()<<endl;
311  cout<<"error= "<<sqrt((float)badLayCount)/ttree->GetEntries()<<endl;
312  cout<<"Signal rejection through clus cut= "<<(float)badClusCount/ttree->GetEntries()<<endl;
313  cout<<"error= "<<sqrt((float)badClusCount)/ttree->GetEntries()<<endl;
314  cout<<"Signal rejection through deta cut= "<<(float)bigDetaCount/ttree->GetEntries()<<endl;
315  cout<<"error= "<<sqrt((float)bigDetaCount)/ttree->GetEntries()<<endl;
316  cout<<"Signal rejection through radius cut= "<<(float)shortRadiusCount/ttree->GetEntries()<<endl;
317  cout<<"error= "<<sqrt((float)shortRadiusCount)/ttree->GetEntries()<<endl;
318  out_file->Write();
319 }
320 
321 void makeRefitDist(TChain* ttree, TFile *out_file){
322  float diffx;
323  float diffy;
324  float diffz;
325  ttree->SetBranchAddress("refitdiffx",&diffx);
326  ttree->SetBranchAddress("refitdiffy",&diffy);
327  ttree->SetBranchAddress("refitdiffz",&diffz);
328  TH1F *diffplotx= new TH1F("x_{0}-x_{refit}","",40,-10,10);
329  TH1F *diffploty= new TH1F("y_{0}-y_{refit}","",40,-10,10);
330  TH1F *diffplotz= new TH1F("z_{0}-z_{refit}","",40,-10,10);
331 
332  diffplotx->Sumw2();
333  diffploty->Sumw2();
334  diffplotz->Sumw2();
335 
336  for (int event = 0; event < ttree->GetEntries(); ++event)
337  {
338  ttree->GetEvent(event);
339  diffplotx->Fill(diffx);
340  diffploty->Fill(diffy);
341  diffplotz->Fill(diffz);
342  }
343 
344  diffplotx->Scale(1./ttree->GetEntries(),"width");
345  diffploty->Scale(1./ttree->GetEntries(),"width");
346  diffplotz->Scale(1./ttree->GetEntries(),"width");
347 
348  out_file->Write();
349  ttree->ResetBranchAddresses();
350 }
351 
352 void makepTCaloGraph(string filename,TFile* outfile){
353  ifstream caloFile;
354  caloFile.open(filename.c_str());
355  double x,y;
356  string s;
357  vector<double> xData, yData;
358  /*if(!(caloFile >>x>>y)){
359  cout<<"file error"<<endl;
360  if(!caloFile.is_open()) cout<<"file not opened"<<endl;
361  }*/
362  while(caloFile >>x>>s>>y){
363  xData.push_back(x);
364  yData.push_back(y);
365  }
366  double *xArray, *yArray;
367  xArray=&xData[0];
368  yArray=&yData[0];
369  TGraph *pTResCaloGraph = new TGraph(xData.size(),xArray,yArray);
370  pTResCaloGraph->SetNameTitle("calopTRes","calopTRes");
371  pTResCaloGraph->Sort();
372  pTResCaloGraph->Write();
373  outfile->Write();
374 }
375 
376 TH1F* makePythiaSpec(TChain* ttree,TFile* out_file,string type=""){
377  vector<float>* tpT= NULL;
378  ttree->SetBranchAddress("photon_pT",&tpT);
379  string title;
380  if(!type.empty()){
381  title=type;
382  title+="_photon_truth_pT";
383  }
384  else title="photon_truth_pT";
385  TH1F *tpTspec = new TH1F(title.c_str(),"",20,5,25);
386  tpTspec->Sumw2();
387  cout<<"pythia tree with: "<<ttree->GetEntries()<<" entries"<<endl;
388  for (int event = 0; event < ttree->GetEntries(); ++event)
389  {
390  ttree->GetEvent(event);
391  for(auto i: *tpT){
392  if(i>5)tpTspec->Fill(i);
393  }
394  }
395  out_file->Write();
396  ttree->ResetBranchAddresses();
397  return tpTspec;
398 }
399 
400 TH1F removeFirstBin(TH1F h){
401  if(h.GetNbinsX()>1){
402  string rname = h.GetName();
403  rname+=to_string(replot++);
404  TH1F r = TH1F(rname.c_str(),h.GetTitle(),h.GetNbinsX()-1,h.GetBinLowEdge(2),h.GetBinLowEdge(h.GetNbinsX()+1));
405  r.Sumw2();
406  for(unsigned i=1; i<r.GetNbinsX();++i){
407  r.SetBinContent(i,h.GetBinContent(i+1));
408  r.SetBinError(i,h.GetBinError(i+1));
409  }
410  return r;
411  }
412  else {
413  return h;
414  }
415 }
416 
417 double ADtoP(double ad){
418  cout<<"AD:"<<ad<<endl;
419  double r=-1;
420  if(ad<.2) r= 1-TMath::Exp(-13.436+101.14*ad-223.74*ad*ad);
421  else if (ad<.34) r= 1-TMath::Exp(-8.31+42.79*ad-59.938*ad*ad);
422  else if (ad<.6) r= 1-TMath::Exp(.917-4.279*ad-1.38*ad*ad);
423  else r= 1-TMath::Exp(1.29-5.709*ad+.0186*ad*ad);
424  cout<<"P:"<<r<<endl;
425  return r;
426 }
427 
428 float smartChi2(TH1F hard, TH1F soft){
429  hard.Scale(1/hard.Integral());
430  soft.Scale(1/soft.Integral());
431  return soft.Chi2Test(&hard,"WW CHI2/NDF");
432 }
433 
434 void chopHard(TH1F hard,TH1F soft){
435  unsigned bins = soft.GetNbinsX();
436  while(ADtoP(hard.AndersonDarlingTest(&soft))>.05&&bins>2){
437  cout<<"Hard/Soft p="<<ADtoP(hard.AndersonDarlingTest(&soft))<<" n bins="<<hard.GetNbinsX()<<endl;
438  cout<<"chi2 p="<<smartChi2(hard,soft)<<" n bins="<<hard.GetNbinsX()<<endl;
439  hard=removeFirstBin(hard);
440  soft=removeFirstBin(soft);
441  bins--;
442  }
443  TH1F* temp = (TH1F*) hard.Clone("hard_chopped");
444  temp->Write();
445 }
446 
447 void calculateConversionRate(TEfficiency* rate, TH1F *pythia,TFile* out_file){
448  TH1F* conversion_rate = (TH1F*) pythia->Clone("rate");
449  TH1* uni_rate = (TH1F*)rate->GetPassedHistogram()->Clone("uni_rate");
450  uni_rate->Divide(rate->GetTotalHistogram());
451  conversion_rate->Multiply(uni_rate);
452  conversion_rate->Scale(1/pythia->Integral());
453  out_file->Write();
454 }
455 
456 double hardWeightFactor(TH1F* hard,TH1F* soft, unsigned matchBin){
457  cout<<soft->Integral(matchBin,soft->GetNbinsX())<<endl;
458  cout<<hard->Integral(matchBin,soft->GetNbinsX())<<endl;
459  return soft->Integral(matchBin,soft->GetNbinsX())/hard->Integral(matchBin,hard->GetNbinsX());
460 }
461 
462 TH1F* addSpec(TH1F* soft,TH1F* hard,TFile* file){
463  TH1F* combined = new TH1F("combinedpythia","",soft->GetNbinsX(),soft->GetBinLowEdge(1),hard->GetBinLowEdge(hard->GetNbinsX()+1));
464  unsigned matchBin = 11;//getMatchingBin(hard,soft); //hard coded by eye
465  for (int i = 1; i < matchBin; ++i)
466  {
467  combined->SetBinContent(i,soft->GetBinContent(i));
468  combined->SetBinError(i,soft->GetBinError(i));
469  }
470  hard->Scale(hardWeightFactor(hard,soft,matchBin));
471  for (int i = matchBin; i < combined->GetNbinsX()+1; ++i)
472  {
473  combined->SetBinContent(i,hard->GetBinContent(i));
474  combined->SetBinError(i,hard->GetBinError(i));
475  }
476  file->Write();
477  return combined;
478 }
479 
480 void photonEff2()
481 {
482  TFile *out_file = new TFile("effplots2.root","UPDATE");
483  string treePath = "/sphenix/user/vassalli/gammasample/embeded/truthconversionembededanalysis";
484  string treeExtension = ".root";
485  unsigned int nFiles=1000;
486  //string softPath = "/sphenix/user/vassalli/minBiasPythia/softana.root";
487  //string hard0Path = "/sphenix/user/vassalli/minBiasPythia/hard0ana.root";
488  //string hard4Path = "/sphenix/user/vassalli/minBiasPythia/hard4ana.root";
489  //TChain *softTree = new TChain("photonTree");
490  //TChain *hard0Tree = new TChain("photonTree");
491  //TChain *hard4Tree = new TChain("photonTree");
492  //softTree->Add(softPath.c_str());
493  //hard0Tree->Add(hard0Path.c_str());
494  //hard4Tree->Add(hard4Path.c_str());
495  TChain *ttree = handleFile(treePath,treeExtension,"cutTreeSignal",nFiles);
496  TChain *pairBackTree = handleFile(treePath,treeExtension,"pairBackTree",nFiles);
497  TChain *vtxBackTree = handleFile(treePath,treeExtension,"vtxBackTree",nFiles);
498  TChain *vertexingTree = handleFile(treePath,treeExtension,"vtxingTree",nFiles);
499  makephotonM(ttree,vtxBackTree,out_file);
500  //auto pythiaSpec=makePythiaSpec(softTree,out_file,"soft");
501  //makePythiaSpec(hard0Tree,out_file,"hard0");
502  //auto hardSpec = makePythiaSpec(hard4Tree,out_file,"hard4");
503  //chopHard(*hardSpec,*pythiaSpec);
504  //pythiaSpec = addSpec(pythiaSpec,hardSpec,out_file);
505  //makepTDist(ttree,observations,out_file);
506  makeVtxRes(ttree,out_file);
507  makeVtxEff(ttree,out_file);
508  makepTRes(ttree);
509  compareDeta(ttree,pairBackTree);
510  //testCuts(ttree,out_file);
511  makepTCaloGraph("pTcalodata.csv",out_file);
512  makeVtxR(ttree,vertexingTree,out_file);
513  //makeRefitDist(ttree,out_file);
514 }