Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
photonEff.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file photonEff.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_m;
39  float back_Vtx;
40  float back_pt;
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}","",60,0,.18));
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  for (int i = 0; i < background->GetEntries(); ++i)
242  {
243  background->GetEvent(i);
244  detaB_plot->Fill(detab);
245  }
246  detaB_plot->Scale(1/detaB_plot->Integral());
247  detaS_plot->Scale(1/detaS_plot->Integral());
248  detaS_plot->Write();
249  detaB_plot->Write();
250 }
251 
252 void testCuts(TChain* ttree,TFile* out_file){
253  float dphi;
254  float prob;
255  float track_pT;
256  float deta;
257  float radius;
258  int layer;
259  int dlayer;
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);
267 
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);
273  layerDist->Sumw2();
274  probDist->Sumw2();
275  deta_plot->Sumw2();
276  dlayer_plot->Sumw2();
277  r_plot->Sumw2();
278  unsigned badLayCount=0;
279  unsigned badClusCount=0;
280  unsigned bigDetaCount=0;
281  unsigned shortRadiusCount=0;
282 
283  for (int event = 0; event < ttree->GetEntries(); ++event)
284  {
285  ttree->GetEvent(event);
286  if(layer==0)badLayCount++;
287  if(dphi<0&&track_pT>.6){
288  badClusCount++;
289  }
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++;
296  else{
297  r_plot->Fill(radius);
298  if(radius<1.84) shortRadiusCount++;
299  }
300  }
301 
302  }
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;
315  out_file->Write();
316 }
317 
318 void makeRefitDist(TChain* ttree, TFile *out_file){
319  float diffx;
320  float diffy;
321  float diffz;
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);
328 
329  diffplotx->Sumw2();
330  diffploty->Sumw2();
331  diffplotz->Sumw2();
332 
333  for (int event = 0; event < ttree->GetEntries(); ++event)
334  {
335  ttree->GetEvent(event);
336  diffplotx->Fill(diffx);
337  diffploty->Fill(diffy);
338  diffplotz->Fill(diffz);
339  }
340 
341  diffplotx->Scale(1./ttree->GetEntries(),"width");
342  diffploty->Scale(1./ttree->GetEntries(),"width");
343  diffplotz->Scale(1./ttree->GetEntries(),"width");
344 
345  out_file->Write();
346  ttree->ResetBranchAddresses();
347 }
348 
349 void makepTCaloGraph(string filename,TFile* outfile){
350  ifstream caloFile;
351  caloFile.open(filename.c_str());
352  double x,y;
353  string s;
354  vector<double> xData, yData;
355  /*if(!(caloFile >>x>>y)){
356  cout<<"file error"<<endl;
357  if(!caloFile.is_open()) cout<<"file not opened"<<endl;
358  }*/
359  while(caloFile >>x>>s>>y){
360  xData.push_back(x);
361  yData.push_back(y);
362  }
363  double *xArray, *yArray;
364  xArray=&xData[0];
365  yArray=&yData[0];
366  TGraph *pTResCaloGraph = new TGraph(xData.size(),xArray,yArray);
367  pTResCaloGraph->SetNameTitle("calopTRes","calopTRes");
368  pTResCaloGraph->Sort();
369  pTResCaloGraph->Write();
370  outfile->Write();
371 }
372 
373 TH1F* makePythiaSpec(TChain* ttree,TFile* out_file,string type=""){
374  vector<float>* tpT= NULL;
375  ttree->SetBranchAddress("photon_pT",&tpT);
376  string title;
377  if(!type.empty()){
378  title=type;
379  title+="_photon_truth_pT";
380  }
381  else title="photon_truth_pT";
382  TH1F *tpTspec = new TH1F(title.c_str(),"",20,5,25);
383  tpTspec->Sumw2();
384  cout<<"pythia tree with: "<<ttree->GetEntries()<<" entries"<<endl;
385  for (int event = 0; event < ttree->GetEntries(); ++event)
386  {
387  ttree->GetEvent(event);
388  for(auto i: *tpT){
389  if(i>5)tpTspec->Fill(i);
390  }
391  }
392  out_file->Write();
393  ttree->ResetBranchAddresses();
394  return tpTspec;
395 }
396 
397 TH1F removeFirstBin(TH1F h){
398  if(h.GetNbinsX()>1){
399  string rname = h.GetName();
400  rname+=to_string(replot++);
401  TH1F r = TH1F(rname.c_str(),h.GetTitle(),h.GetNbinsX()-1,h.GetBinLowEdge(2),h.GetBinLowEdge(h.GetNbinsX()+1));
402  r.Sumw2();
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));
406  }
407  return r;
408  }
409  else {
410  return h;
411  }
412 }
413 
414 double ADtoP(double ad){
415  cout<<"AD:"<<ad<<endl;
416  double r=-1;
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);
421  cout<<"P:"<<r<<endl;
422  return r;
423 }
424 
425 float smartChi2(TH1F hard, TH1F soft){
426  hard.Scale(1/hard.Integral());
427  soft.Scale(1/soft.Integral());
428  return soft.Chi2Test(&hard,"WW CHI2/NDF");
429 }
430 
431 void chopHard(TH1F hard,TH1F soft){
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;
436  hard=removeFirstBin(hard);
437  soft=removeFirstBin(soft);
438  bins--;
439  }
440  TH1F* temp = (TH1F*) hard.Clone("hard_chopped");
441  temp->Write();
442 }
443 
444 void calculateConversionRate(TEfficiency* rate, TH1F *pythia,TFile* out_file){
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());
450  out_file->Write();
451 }
452 
453 double hardWeightFactor(TH1F* hard,TH1F* soft, unsigned matchBin, double* error){
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;
460  return r;
461 }
462 
463 TH1F* addSpec(TH1F* soft,TH1F* hard,TFile* file){
464  TH1F* combined = new TH1F("combinedpythia","",soft->GetNbinsX(),soft->GetBinLowEdge(1),hard->GetBinLowEdge(hard->GetNbinsX()+1));
465  unsigned matchBin = 10;//getMatchingBin(hard,soft); //hard coded by eye
466  for (int i = 1; i < matchBin; ++i)
467  {
468  combined->SetBinContent(i,soft->GetBinContent(i));
469  combined->SetBinError(i,soft->GetBinError(i));
470  }
471  double systematic;
472  hard->Scale(hardWeightFactor(hard,soft,matchBin,&systematic));
473  for (int i = matchBin; i < combined->GetNbinsX()+1; ++i)
474  {
475  combined->SetBinContent(i,hard->GetBinContent(i));
476  combined->SetBinError(i,sqrt(hard->GetBinError(i)*hard->GetBinError(i)+systematic*systematic));
477  }
478  file->Write();
479  return combined;
480 }
481 
482 void photonEff()
483 {
484  TFile *out_file = new TFile("effplots.root","UPDATE");
485  //string treePath = "/sphenix/user/vassalli/gammasample/truthconversionembededanalysis";
486  string treeExtension = ".root";
487  unsigned int nFiles=1000;
488  string softPath = "/sphenix/user/vassalli/minBiasPythia/softana.root";
489  //string hard0Path = "/sphenix/user/vassalli/minBiasPythia/hard0ana.root";
490  string hard4Path = "/sphenix/user/vassalli/minBiasPythia/hard4ana.root";
491  TChain *softTree = new TChain("photonTree");
492  //TChain *hard0Tree = new TChain("photonTree");
493  TChain *hard4Tree = new TChain("photonTree");
494  softTree->Add(softPath.c_str());
495  //hard0Tree->Add(hard0Path.c_str());
496  hard4Tree->Add(hard4Path.c_str());
497  /*TChain *ttree = handleFile(treePath,treeExtension,"cutTreeSignal",nFiles);
498  TChain *pairBackTree = handleFile(treePath,treeExtension,"pairBackTree",nFiles);
499  TChain *vtxBackTree = handleFile(treePath,treeExtension,"vtxBackTree",nFiles);
500  TChain *vertexingTree = handleFile(treePath,treeExtension,"vtxingTree",nFiles);
501  makephotonM(ttree,vtxBackTree,out_file);*/
502  auto pythiaSpec=makePythiaSpec(softTree,out_file,"soft");
503  //makePythiaSpec(hard0Tree,out_file,"hard0");
504  auto hardSpec = makePythiaSpec(hard4Tree,out_file,"hard4");
505  //chopHard(*hardSpec,*pythiaSpec);
506  pythiaSpec = addSpec(pythiaSpec,hardSpec,out_file);
507  //makepTDist(ttree,observations,out_file);
508  /*makeVtxRes(ttree,out_file);
509  makeVtxEff(ttree,out_file);
510  makepTRes(ttree);
511  compareDeta(ttree,pairBackTree);
512  //testCuts(ttree,out_file);
513  makepTCaloGraph("pTcalodata.csv",out_file);
514  makeVtxR(ttree,vertexingTree,out_file);*/
515  //makeRefitDist(ttree,out_file);
516 }