PhotonConversion/ana/photonEff2.C

0001 #include <iostream>
0002 #include <fstream>
0003 using namespace std;
0004
0005 #include "TFile.h"
0006 #include "Scalar.h"
0007 #include "TTree.h"
0008 #include "TChain.h"
0009 #include "TLegend.h"
0010 #include "math.h"
0011 #include "TH1.h"
0012 #include "TH2.h"
0013 #include "TEfficiency.h"
0014 #include "TLine.h"
0015 #include "TGraphAsymmErrors.h"
0016
0017 namespace {
0018   unsigned replot=0;
0019 }
0020
0021 TChain* handleFile(string name, string extension, string treename, unsigned int filecount){
0022   TChain *all = new TChain(treename.c_str());
0023   string temp;
0024   for (int i = 0; i < filecount; ++i)
0025   {
0026
0027     ostringstream s;
0028     s<<i;
0029     temp = name+string(s.str())+extension;
0030     all->Add(temp.c_str());
0031   }
0032   return all;
0033 }
0034
0035 void makephotonM(TChain* ttree,TChain* back, TFile* out_file){
0036   float photon_m;
0037   float tphoton_m;
0038   float back_Vtx;
0039   float back_pt;
0040   float back_m;
0041   std::vector<TH1F*> plots;
0042   ttree->SetBranchAddress("photon_m",     &photon_m   );
0043   back->SetBranchAddress("photon_m",     &back_m   );
0044   back->SetBranchAddress("vtx_radius",     &back_Vtx   );
0045   back->SetBranchAddress("track_pT",     &back_pt   );
0046   //ttree->SetBranchAddress("rephoton_m",     &rephoton_m   );
0047   plots.push_back(new TH1F("m^{#gamma}_{reco}","",60,0,.18));
0048   plots.push_back(new TH1F("m^{bkgd}_{reco}_wide","",60,0,1));
0049   //plots.push_back(new TH1F("m^{#gamma}_{recoRefit}","",40,-2,10));
0050
0051   for (int i = 0; i < plots.size(); ++i)
0052   {
0053     plots[i]->Sumw2();
0054   }
0055   for (int event = 0; event < ttree->GetEntries(); ++event)
0056   {
0057     ttree->GetEvent(event);
0058     plots[0]->Fill(photon_m);
0059     // plots[1]->Fill(rephoton_m);
0060   }
0061   for (int event = 0; event < back->GetEntries(); ++event)
0062   {
0063     back->GetEvent(event);
0064     if(back_Vtx<21.&&back_pt>2.5)
0065     plots[1]->Fill(back_m);
0066     // plots[1]->Fill(rephoton_m);
0067   }
0068  // cout<<"Background mass in range count= "<<plots[1]->Integral()<<'\n';
0069   for (int i = 0; i < plots.size(); ++i)
0070   {
0071     plots[i]->Scale(1./ttree->GetEntries(),"width");
0072   }
0073   out_file->Write();
0074   ttree->ResetBranchAddresses();
0075 }
0076
0077 void makeVtxR(TChain* ttree, TTree* vtxTree,TFile* out_file){
0078   float vtxr;
0079   float vtxX;
0080   float vtxY;
0081   float tvtxr;
0082   ttree->SetBranchAddress("vtx_radius",&vtxr);
0083   ttree->SetBranchAddress("tvtx_radius",&tvtxr);
0084   vtxTree->SetBranchAddress("vtx_x",&vtxX);
0085   vtxTree->SetBranchAddress("vtx_y",&vtxY);
0086
0087   std::vector<TH1F*> plots;
0088   plots.push_back(new TH1F("vtx_reco","",30,0,40));
0089   plots.push_back(new TH1F("vtx_corrected","",30,0,40));
0090   plots.push_back(new TH1F("vtx_truth","",30,0,40));
0091
0092   plots[0]->Sumw2();
0093   plots[1]->Sumw2();
0094   plots[2]->Sumw2();
0095
0096   double calc=0;
0097   for (int event = 0; event < ttree->GetEntries(); ++event)
0098   {
0099     ttree->GetEvent(event);
0100     plots[1]->Fill(vtxr);
0101     plots[2]->Fill(tvtxr);
0102     calc+=TMath::Abs(vtxr-tvtxr);
0103   }
0104   calc/=ttree->GetEntries();
0105   for (int i = 0; i < vtxTree->GetEntries(); ++i)
0106   {
0107     vtxTree->GetEvent(i);
0108     plots[0]->Fill(sqrt(vtxX*vtxX+vtxY*vtxY));
0109   }
0110   for (int i = 0; i < 3; ++i)
0111   {
0112     plots[i]->Scale(1./plots[2]->GetEntries(),"width");
0113   }
0114   out_file->Write();
0115   std::cout<<"mean deviation="<<calc<<std::endl;
0116 }
0117
0118 void makeVtxRes(TChain* ttree,TFile* out_file){
0119   float r;
0120   float tr;
0121   ttree->SetBranchAddress("vtx_radius",&r);
0122   ttree->SetBranchAddress("tvtx_radius",&tr);
0123   TH1F *vtxeffPlot = new TH1F("#frac{#Deltar_{vtx}_^{#it{reco}}}{r_{vtx}^{#it{truth}}}","",40,-2,2);
0124   TH2F *vtxefffuncPlot = new TH2F("vtx_resolution_to_truthvtx","",20,0,21,40,-1.5,1.5);
0125   vtxeffPlot->Sumw2();
0126   vtxefffuncPlot->Sumw2();
0127   for (int event = 0; event < ttree->GetEntries(); ++event)
0128   {
0129     ttree->GetEvent(event);
0130     if(r<0) continue;
0131     vtxeffPlot->Fill((r-tr)/tr);
0132     vtxefffuncPlot->Fill(tr,(r-tr)/tr);
0133   }
0134   vtxeffPlot->Scale(1./ttree->GetEntries(),"width");
0135   vtxefffuncPlot->Scale(1./ttree->GetEntries(),"width");
0136   out_file->Write();
0137   ttree->ResetBranchAddresses();
0138 }
0139
0140 void makeVtxEff(TChain* ttree,TFile* out_file){
0141   float r;
0142   float tr;
0143   float pT;
0144   ttree->SetBranchAddress("vtx_radius",&r);
0145   ttree->SetBranchAddress("tvtx_radius",&tr);
0146   ttree->SetBranchAddress("track_pT",&pT);
0147   TEfficiency *vtxEff;
0148   TH1F *recoR= new TH1F("vtxrecoR","",30,0,40);
0149   TH1F *truthR= new TH1F("vtxtruthR","",30,0,40);
0150   recoR->Sumw2();
0151   truthR->Sumw2();
0152   for (int event = 0; event < ttree->GetEntries(); ++event)
0153   {
0154     ttree->GetEvent(event);
0155     if(r>0) recoR->Fill(tr);
0156     truthR->Fill(tr);
0157   }
0158   vtxEff = new TEfficiency(*recoR,*truthR);
0159   vtxEff->SetName("vtxEff");
0160   vtxEff->Write();
0161   out_file->Write();
0162   ttree->ResetBranchAddresses();
0163 }
0164
0165 void makepTRes(TChain* ttree){
0166   float pT;
0167   float tpT;
0168   ttree->SetBranchAddress("photon_pT",&pT);
0169   ttree->SetBranchAddress("tphoton_pT",&tpT);
0170   TH1F *pTRes = new TH1F("#frac{#it{p}^{T}_{reco}}{#it{p}_{#it{truth}}^{T}}","",40,0,2.5);
0171   pTRes->Sumw2();
0172   for (int event = 0; event < ttree->GetEntries(); ++event)
0173   {
0174     ttree->GetEvent(event);
0175     if(pT>0)pTRes->Fill(pT/tpT);
0176   }
0177   pTRes->Scale(1./pTRes->Integral());
0178   ttree->ResetBranchAddresses();
0179   pTRes->Write();
0180 }
0181
0182 TEfficiency* makepTDist(TChain* ttree,TTree* allTree,TFile* out_file){
0183   float pT;
0184   float tpT;
0185   float track_pT;
0186   ttree->SetBranchAddress("photon_pT",&pT);
0187   ttree->SetBranchAddress("tphoton_pT",&tpT);
0188
0189   vector<float> *allpT=NULL;
0190   allTree->SetBranchAddress("photon_pT",&allpT);
0191
0192   TH1F *pTeffPlot = new TH1F("#frac{#it{p}^{T}}{#it{p}_{#it{truth}}^{T}}","",40,-2,2);
0193   TH2F *pTefffuncPlot = new TH2F("pT_resolution_to_truthpt","",40,1,35,40,-1.5,1.5);
0194   TH1F *converted_pTspec = new TH1F("converted_photon_truth_pT","",25,5,30);
0195   TH1F *all_pTspec = new TH1F("all_photon_truth_pT","",25,5,30);
0196   //TH1F *trackpTDist = new TH1F("truthpt","",40,0,35);
0197   pTeffPlot->Sumw2();
0198   converted_pTspec->Sumw2();
0199   all_pTspec->Sumw2();
0200   pTefffuncPlot->Sumw2();
0201   for (int event = 0; event < ttree->GetEntries(); ++event)
0202   {
0203     ttree->GetEvent(event);
0204     if(pT>0)pTeffPlot->Fill(pT/tpT);
0205     converted_pTspec->Fill(tpT);
0206     if(pT>0)pTefffuncPlot->Fill(tpT,pT/tpT);
0207     //trackpTDist->Fill(track_pT);
0208   }
0209   for (int event = 0; event < allTree->GetEntries(); ++event)
0210   {
0211     allTree->GetEvent(event);
0212     for(auto i : *allpT){
0213       all_pTspec->Fill(i);
0214     }
0215   }
0216   TEfficiency* uni_rate = new TEfficiency(*converted_pTspec,*all_pTspec);
0217   pTeffPlot->Scale(1./ttree->GetEntries(),"width");
0218   pTefffuncPlot->Scale(1./ttree->GetEntries());
0219   TProfile* resProfile = pTefffuncPlot->ProfileX("func_prof",5,30);
0220   resProfile->Write();
0221   //trackpTDist->Scale(1./ttree->GetEntries(),"width");
0222   out_file->Write();
0223   ttree->ResetBranchAddresses();
0224   return uni_rate;
0225 }
0226
0227 void compareDeta(TTree* signalTree, TTree* background){
0228   float detas,detab;
0229   signalTree->SetBranchAddress("track_deta",&detas);
0230   background->SetBranchAddress("track_deta",&detab);
0231   TH1F *detaS_plot = new TH1F("detaS","",3000,-.001,1);
0232   TH1F *detaB_plot = new TH1F("detaB","",3000,-.001,1);
0233   detaS_plot->Sumw2();
0234   detaB_plot->Sumw2();
0235
0236   for (int i = 0; i < signalTree->GetEntries(); ++i)
0237   {
0238     signalTree->GetEvent(i);
0239     detaS_plot->Fill(detas);
0240   }
0241   cout<<"background has "<<background->GetEntries()<<" entries"<<endl;
0242   for (int i = 0; i < background->GetEntries(); ++i)
0243   {
0244     background->GetEvent(i);
0245     detaB_plot->Fill(detab);
0246   }
0247   detaB_plot->Scale(1/detaB_plot->Integral());
0248   detaS_plot->Scale(1/detaS_plot->Integral());
0249   detaS_plot->Write();
0250   detaB_plot->Write();
0251   signalTree->ResetBranchAddresses();
0252   background->ResetBranchAddresses();
0253 }
0254
0255 void testCuts(TChain* ttree,TFile* out_file){
0256   float dphi;
0257   float prob;
0258   float track_pT;
0259   float deta;
0260   float radius;
0261   int layer;
0262   int dlayer;
0263   ttree->SetBranchAddress("cluster_dphi",&dphi);
0264   ttree->SetBranchAddress("cluster_prob",&prob);
0265   ttree->SetBranchAddress("track_layer",&layer);
0266   ttree->SetBranchAddress("track_dlayer",&dlayer);
0267   ttree->SetBranchAddress("track_pT",&track_pT);
0268   ttree->SetBranchAddress("track_deta",&deta);
0269   ttree->SetBranchAddress("vtx_radius",&radius);
0270
0271   TH1F *layerDist = new TH1F("layer","",16,-.5,15.5);
0272   TH1F *probDist = new TH1F("clust_prob","",30,-.5,1.);
0273   TH1F *deta_plot = new TH1F("deta","",30,-.001,.01);
0274   TH1F *dlayer_plot = new TH1F("dlayer","",11,-.5,10.5);
0275   TH1F *r_plot = new TH1F("signal_vtx_radius_dist","",26,-.5,25.5);
0276   layerDist->Sumw2();
0277   probDist->Sumw2();
0278   deta_plot->Sumw2();
0279   dlayer_plot->Sumw2();
0280   r_plot->Sumw2();
0281   unsigned badLayCount=0;
0282   unsigned badClusCount=0;
0283   unsigned bigDetaCount=0;
0284   unsigned shortRadiusCount=0;
0285
0286   for (int event = 0; event < ttree->GetEntries(); ++event)
0287   {
0288     ttree->GetEvent(event);
0289     if(layer==0)badLayCount++;
0290     if(dphi<0&&track_pT>.6){
0291       badClusCount++;
0292     }
0293     if(track_pT>.6&&dphi>=0){
0294       layerDist->Fill(layer);
0295       probDist->Fill(prob);
0296       deta_plot->Fill(deta);
0297       dlayer_plot->Fill(TMath::Abs(dlayer));
0298       if(deta>.0082||TMath::Abs(dlayer)>9)bigDetaCount++;
0299       else{
0300         r_plot->Fill(radius);
0301         if(radius<1.84) shortRadiusCount++;
0302       }
0303     }
0304
0305   }
0306   layerDist->Scale(1./ttree->GetEntries());
0307   deta_plot->Scale(1./ttree->GetEntries());
0308   dlayer_plot->Scale(1./ttree->GetEntries());
0309   r_plot->Scale(1./ttree->GetEntries());
0310   cout<<"Signal rejection through layer cut= "<<(float)badLayCount/ttree->GetEntries()<<endl;
0311   cout<<"error= "<<sqrt((float)badLayCount)/ttree->GetEntries()<<endl;
0312   cout<<"Signal rejection through clus cut= "<<(float)badClusCount/ttree->GetEntries()<<endl;
0313   cout<<"error= "<<sqrt((float)badClusCount)/ttree->GetEntries()<<endl;
0314   cout<<"Signal rejection through deta cut= "<<(float)bigDetaCount/ttree->GetEntries()<<endl;
0315   cout<<"error= "<<sqrt((float)bigDetaCount)/ttree->GetEntries()<<endl;
0316   cout<<"Signal rejection through radius cut= "<<(float)shortRadiusCount/ttree->GetEntries()<<endl;
0317   cout<<"error= "<<sqrt((float)shortRadiusCount)/ttree->GetEntries()<<endl;
0318   out_file->Write();
0319 }
0320
0321 void makeRefitDist(TChain* ttree, TFile *out_file){
0322   float diffx;
0323   float diffy;
0324   float diffz;
0325   ttree->SetBranchAddress("refitdiffx",&diffx);
0326   ttree->SetBranchAddress("refitdiffy",&diffy);
0327   ttree->SetBranchAddress("refitdiffz",&diffz);
0328   TH1F *diffplotx= new TH1F("x_{0}-x_{refit}","",40,-10,10);
0329   TH1F *diffploty= new TH1F("y_{0}-y_{refit}","",40,-10,10);
0330   TH1F *diffplotz= new TH1F("z_{0}-z_{refit}","",40,-10,10);
0331
0332   diffplotx->Sumw2();
0333   diffploty->Sumw2();
0334   diffplotz->Sumw2();
0335
0336   for (int event = 0; event < ttree->GetEntries(); ++event)
0337   {
0338     ttree->GetEvent(event);
0339     diffplotx->Fill(diffx);
0340     diffploty->Fill(diffy);
0341     diffplotz->Fill(diffz);
0342   }
0343
0344   diffplotx->Scale(1./ttree->GetEntries(),"width");
0345   diffploty->Scale(1./ttree->GetEntries(),"width");
0346   diffplotz->Scale(1./ttree->GetEntries(),"width");
0347
0348   out_file->Write();
0349   ttree->ResetBranchAddresses();
0350 }
0351
0352 void makepTCaloGraph(string filename,TFile* outfile){
0353   ifstream caloFile;
0354   caloFile.open(filename.c_str());
0355   double x,y;
0356   string s;
0357   vector<double> xData, yData;
0358   /*if(!(caloFile >>x>>y)){
0359     cout<<"file error"<<endl;
0360     if(!caloFile.is_open()) cout<<"file not opened"<<endl;
0361     }*/
0362   while(caloFile >>x>>s>>y){
0363     xData.push_back(x);
0364     yData.push_back(y);
0365   }
0366   double *xArray, *yArray;
0367   xArray=&xData[0];
0368   yArray=&yData[0];
0369   TGraph *pTResCaloGraph = new TGraph(xData.size(),xArray,yArray);
0370   pTResCaloGraph->SetNameTitle("calopTRes","calopTRes");
0371   pTResCaloGraph->Sort();
0372   pTResCaloGraph->Write();
0373   outfile->Write();
0374 }
0375
0376 TH1F* makePythiaSpec(TChain* ttree,TFile* out_file,string type=""){
0377   vector<float>* tpT= NULL;
0378   ttree->SetBranchAddress("photon_pT",&tpT);
0379   string title;
0380   if(!type.empty()){
0381     title=type;
0382     title+="_photon_truth_pT";
0383   }
0384   else title="photon_truth_pT";
0385   TH1F *tpTspec = new TH1F(title.c_str(),"",20,5,25);
0386   tpTspec->Sumw2();
0387   cout<<"pythia tree with: "<<ttree->GetEntries()<<" entries"<<endl;
0388   for (int event = 0; event < ttree->GetEntries(); ++event)
0389   {
0390     ttree->GetEvent(event);
0391     for(auto i: *tpT){
0392       if(i>5)tpTspec->Fill(i);
0393     }
0394   }
0395   out_file->Write();
0396   ttree->ResetBranchAddresses();
0397   return tpTspec;
0398 }
0399
0400 TH1F removeFirstBin(TH1F h){
0401   if(h.GetNbinsX()>1){
0402     string rname = h.GetName();
0403     rname+=to_string(replot++);
0404     TH1F r = TH1F(rname.c_str(),h.GetTitle(),h.GetNbinsX()-1,h.GetBinLowEdge(2),h.GetBinLowEdge(h.GetNbinsX()+1));
0405     r.Sumw2();
0406     for(unsigned i=1; i<r.GetNbinsX();++i){
0407       r.SetBinContent(i,h.GetBinContent(i+1));
0408       r.SetBinError(i,h.GetBinError(i+1));
0409     }
0410     return r;
0411   }
0412   else {
0413     return h;
0414   }
0415 }
0416
0417 double ADtoP(double ad){
0418   cout<<"AD:"<<ad<<endl;
0419   double r=-1;
0420   if(ad<.2) r= 1-TMath::Exp(-13.436+101.14*ad-223.74*ad*ad);
0421   else if (ad<.34) r= 1-TMath::Exp(-8.31+42.79*ad-59.938*ad*ad);
0422   else if (ad<.6) r= 1-TMath::Exp(.917-4.279*ad-1.38*ad*ad);
0423   else r= 1-TMath::Exp(1.29-5.709*ad+.0186*ad*ad);
0424   cout<<"P:"<<r<<endl;
0425   return r;
0426 }
0427
0428 float smartChi2(TH1F hard, TH1F soft){
0429   hard.Scale(1/hard.Integral());
0430   soft.Scale(1/soft.Integral());
0431   return soft.Chi2Test(&hard,"WW CHI2/NDF");
0432 }
0433
0434 void chopHard(TH1F hard,TH1F soft){
0435   unsigned bins = soft.GetNbinsX();
0436   while(ADtoP(hard.AndersonDarlingTest(&soft))>.05&&bins>2){
0437     cout<<"Hard/Soft p="<<ADtoP(hard.AndersonDarlingTest(&soft))<<" n bins="<<hard.GetNbinsX()<<endl;
0438     cout<<"chi2 p="<<smartChi2(hard,soft)<<" n bins="<<hard.GetNbinsX()<<endl;
0439     hard=removeFirstBin(hard);
0440     soft=removeFirstBin(soft);
0441     bins--;
0442   }
0443   TH1F* temp = (TH1F*) hard.Clone("hard_chopped");
0444   temp->Write();
0445 }
0446
0447 void calculateConversionRate(TEfficiency* rate, TH1F *pythia,TFile* out_file){
0448   TH1F* conversion_rate = (TH1F*)  pythia->Clone("rate");
0449   TH1* uni_rate = (TH1F*)rate->GetPassedHistogram()->Clone("uni_rate");
0450   uni_rate->Divide(rate->GetTotalHistogram());
0451   conversion_rate->Multiply(uni_rate);
0452   conversion_rate->Scale(1/pythia->Integral());
0453   out_file->Write();
0454 }
0455
0456 double hardWeightFactor(TH1F* hard,TH1F* soft, unsigned matchBin){
0457   cout<<soft->Integral(matchBin,soft->GetNbinsX())<<endl;
0458   cout<<hard->Integral(matchBin,soft->GetNbinsX())<<endl;
0459   return soft->Integral(matchBin,soft->GetNbinsX())/hard->Integral(matchBin,hard->GetNbinsX());
0460 }
0461
0462 TH1F* addSpec(TH1F* soft,TH1F* hard,TFile* file){
0463   TH1F* combined = new TH1F("combinedpythia","",soft->GetNbinsX(),soft->GetBinLowEdge(1),hard->GetBinLowEdge(hard->GetNbinsX()+1));
0464   unsigned matchBin = 11;//getMatchingBin(hard,soft); //hard coded by eye
0465   for (int i = 1; i < matchBin; ++i)
0466   {
0467     combined->SetBinContent(i,soft->GetBinContent(i));
0468     combined->SetBinError(i,soft->GetBinError(i));
0469   }
0470   hard->Scale(hardWeightFactor(hard,soft,matchBin));
0471   for (int i = matchBin; i < combined->GetNbinsX()+1; ++i)
0472   {
0473     combined->SetBinContent(i,hard->GetBinContent(i));
0474     combined->SetBinError(i,hard->GetBinError(i));
0475   }
0476   file->Write();
0477   return combined;
0478 }
0479
0480 void photonEff2()
0481 {
0482   TFile *out_file = new TFile("effplots2.root","UPDATE");
0483   string treePath = "/sphenix/user/vassalli/gammasample/embeded/truthconversionembededanalysis";
0484   string treeExtension = ".root";
0485   unsigned int nFiles=1000;
0486   //string softPath = "/sphenix/user/vassalli/minBiasPythia/softana.root";
0487   //string hard0Path = "/sphenix/user/vassalli/minBiasPythia/hard0ana.root";
0488   //string hard4Path = "/sphenix/user/vassalli/minBiasPythia/hard4ana.root";
0489   //TChain *softTree = new TChain("photonTree");
0490   //TChain *hard0Tree = new TChain("photonTree");
0491   //TChain *hard4Tree = new TChain("photonTree");
0492   //softTree->Add(softPath.c_str());
0493   //hard0Tree->Add(hard0Path.c_str());
0494   //hard4Tree->Add(hard4Path.c_str());
0495   TChain *ttree = handleFile(treePath,treeExtension,"cutTreeSignal",nFiles);
0496   TChain *pairBackTree = handleFile(treePath,treeExtension,"pairBackTree",nFiles);
0497   TChain *vtxBackTree = handleFile(treePath,treeExtension,"vtxBackTree",nFiles);
0498   TChain *vertexingTree = handleFile(treePath,treeExtension,"vtxingTree",nFiles);
0499   makephotonM(ttree,vtxBackTree,out_file);
0500   //auto pythiaSpec=makePythiaSpec(softTree,out_file,"soft");
0501   //makePythiaSpec(hard0Tree,out_file,"hard0");
0502   //auto hardSpec = makePythiaSpec(hard4Tree,out_file,"hard4");
0503   //chopHard(*hardSpec,*pythiaSpec);
0504   //pythiaSpec = addSpec(pythiaSpec,hardSpec,out_file);
0505   //makepTDist(ttree,observations,out_file);
0506   makeVtxRes(ttree,out_file);
0507   makeVtxEff(ttree,out_file);
0508   makepTRes(ttree);
0509   compareDeta(ttree,pairBackTree);
0510   //testCuts(ttree,out_file);
0511   makepTCaloGraph("pTcalodata.csv",out_file);
0512   makeVtxR(ttree,vertexingTree,out_file);
0513   //makeRefitDist(ttree,out_file);
0514 }