Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-03 08:14:09

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_m;
0039   float back_Vtx;
0040   float back_pt;
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}","",60,0,.18));
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   for (int i = 0; i < background->GetEntries(); ++i)
0242   {
0243     background->GetEvent(i);
0244     detaB_plot->Fill(detab);
0245   }
0246   detaB_plot->Scale(1/detaB_plot->Integral());
0247   detaS_plot->Scale(1/detaS_plot->Integral());
0248   detaS_plot->Write();
0249   detaB_plot->Write();
0250 }
0251 
0252 void testCuts(TChain* ttree,TFile* out_file){
0253   float dphi;
0254   float prob;
0255   float track_pT;
0256   float deta;
0257   float radius;
0258   int layer;
0259   int dlayer;
0260   ttree->SetBranchAddress("cluster_dphi",&dphi);
0261   ttree->SetBranchAddress("cluster_prob",&prob);
0262   ttree->SetBranchAddress("track_layer",&layer);
0263   ttree->SetBranchAddress("track_dlayer",&dlayer);
0264   ttree->SetBranchAddress("track_pT",&track_pT);
0265   ttree->SetBranchAddress("track_deta",&deta);
0266   ttree->SetBranchAddress("vtx_radius",&radius);
0267 
0268   TH1F *layerDist = new TH1F("layer","",16,-.5,15.5);
0269   TH1F *probDist = new TH1F("clust_prob","",30,-.5,1.);
0270   TH1F *deta_plot = new TH1F("deta","",30,-.001,.01);
0271   TH1F *dlayer_plot = new TH1F("dlayer","",11,-.5,10.5);
0272   TH1F *r_plot = new TH1F("signal_vtx_radius_dist","",26,-.5,25.5);
0273   layerDist->Sumw2();
0274   probDist->Sumw2();
0275   deta_plot->Sumw2();
0276   dlayer_plot->Sumw2();
0277   r_plot->Sumw2();
0278   unsigned badLayCount=0;
0279   unsigned badClusCount=0;
0280   unsigned bigDetaCount=0;
0281   unsigned shortRadiusCount=0;
0282 
0283   for (int event = 0; event < ttree->GetEntries(); ++event)
0284   {
0285     ttree->GetEvent(event);
0286     if(layer==0)badLayCount++;
0287     if(dphi<0&&track_pT>.6){
0288       badClusCount++;
0289     }
0290     if(track_pT>.6&&dphi>=0){
0291       layerDist->Fill(layer);
0292       probDist->Fill(prob);
0293       deta_plot->Fill(deta);
0294       dlayer_plot->Fill(TMath::Abs(dlayer));
0295       if(deta>.0082||TMath::Abs(dlayer)>9)bigDetaCount++;
0296       else{
0297         r_plot->Fill(radius);
0298         if(radius<1.84) shortRadiusCount++;
0299       }
0300     }
0301 
0302   }
0303   layerDist->Scale(1./ttree->GetEntries());
0304   deta_plot->Scale(1./ttree->GetEntries());
0305   dlayer_plot->Scale(1./ttree->GetEntries());
0306   r_plot->Scale(1./ttree->GetEntries());
0307   cout<<"Signal rejection through layer cut= "<<(float)badLayCount/ttree->GetEntries()<<endl;
0308   cout<<"error= "<<sqrt((float)badLayCount)/ttree->GetEntries()<<endl;
0309   cout<<"Signal rejection through clus cut= "<<(float)badClusCount/ttree->GetEntries()<<endl;
0310   cout<<"error= "<<sqrt((float)badClusCount)/ttree->GetEntries()<<endl;
0311   cout<<"Signal rejection through deta cut= "<<(float)bigDetaCount/ttree->GetEntries()<<endl;
0312   cout<<"error= "<<sqrt((float)bigDetaCount)/ttree->GetEntries()<<endl;
0313   cout<<"Signal rejection through radius cut= "<<(float)shortRadiusCount/ttree->GetEntries()<<endl;
0314   cout<<"error= "<<sqrt((float)shortRadiusCount)/ttree->GetEntries()<<endl;
0315   out_file->Write();
0316 }
0317 
0318 void makeRefitDist(TChain* ttree, TFile *out_file){
0319   float diffx;
0320   float diffy;
0321   float diffz;
0322   ttree->SetBranchAddress("refitdiffx",&diffx);
0323   ttree->SetBranchAddress("refitdiffy",&diffy);
0324   ttree->SetBranchAddress("refitdiffz",&diffz);
0325   TH1F *diffplotx= new TH1F("x_{0}-x_{refit}","",40,-10,10);
0326   TH1F *diffploty= new TH1F("y_{0}-y_{refit}","",40,-10,10);
0327   TH1F *diffplotz= new TH1F("z_{0}-z_{refit}","",40,-10,10);
0328 
0329   diffplotx->Sumw2();
0330   diffploty->Sumw2();
0331   diffplotz->Sumw2();
0332 
0333   for (int event = 0; event < ttree->GetEntries(); ++event)
0334   {
0335     ttree->GetEvent(event);
0336     diffplotx->Fill(diffx);
0337     diffploty->Fill(diffy);
0338     diffplotz->Fill(diffz);
0339   }
0340 
0341   diffplotx->Scale(1./ttree->GetEntries(),"width");
0342   diffploty->Scale(1./ttree->GetEntries(),"width");
0343   diffplotz->Scale(1./ttree->GetEntries(),"width");
0344 
0345   out_file->Write();
0346   ttree->ResetBranchAddresses();
0347 }
0348 
0349 void makepTCaloGraph(string filename,TFile* outfile){
0350   ifstream caloFile;
0351   caloFile.open(filename.c_str());
0352   double x,y;
0353   string s;
0354   vector<double> xData, yData;
0355   /*if(!(caloFile >>x>>y)){
0356     cout<<"file error"<<endl;
0357     if(!caloFile.is_open()) cout<<"file not opened"<<endl;
0358     }*/
0359   while(caloFile >>x>>s>>y){
0360     xData.push_back(x);
0361     yData.push_back(y);
0362   }
0363   double *xArray, *yArray;
0364   xArray=&xData[0];
0365   yArray=&yData[0];
0366   TGraph *pTResCaloGraph = new TGraph(xData.size(),xArray,yArray);
0367   pTResCaloGraph->SetNameTitle("calopTRes","calopTRes");
0368   pTResCaloGraph->Sort();
0369   pTResCaloGraph->Write();
0370   outfile->Write();
0371 }
0372 
0373 TH1F* makePythiaSpec(TChain* ttree,TFile* out_file,string type=""){
0374   vector<float>* tpT= NULL;
0375   ttree->SetBranchAddress("photon_pT",&tpT);
0376   string title;
0377   if(!type.empty()){
0378     title=type;
0379     title+="_photon_truth_pT";
0380   }  
0381   else title="photon_truth_pT";
0382   TH1F *tpTspec = new TH1F(title.c_str(),"",20,5,25);
0383   tpTspec->Sumw2();
0384   cout<<"pythia tree with: "<<ttree->GetEntries()<<" entries"<<endl;
0385   for (int event = 0; event < ttree->GetEntries(); ++event)
0386   {
0387     ttree->GetEvent(event);
0388     for(auto i: *tpT){
0389       if(i>5)tpTspec->Fill(i);
0390     }
0391   }
0392   out_file->Write();
0393   ttree->ResetBranchAddresses();
0394   return tpTspec;
0395 }
0396 
0397 TH1F removeFirstBin(TH1F h){
0398   if(h.GetNbinsX()>1){
0399     string rname = h.GetName();
0400     rname+=to_string(replot++);
0401     TH1F r = TH1F(rname.c_str(),h.GetTitle(),h.GetNbinsX()-1,h.GetBinLowEdge(2),h.GetBinLowEdge(h.GetNbinsX()+1));
0402     r.Sumw2();
0403     for(unsigned i=1; i<r.GetNbinsX();++i){
0404       r.SetBinContent(i,h.GetBinContent(i+1));
0405       r.SetBinError(i,h.GetBinError(i+1));
0406     }
0407     return r;
0408   }
0409   else {
0410     return h;
0411   }
0412 }
0413 
0414 double ADtoP(double ad){
0415   cout<<"AD:"<<ad<<endl;
0416   double r=-1;
0417   if(ad<.2) r= 1-TMath::Exp(-13.436+101.14*ad-223.74*ad*ad);
0418   else if (ad<.34) r= 1-TMath::Exp(-8.31+42.79*ad-59.938*ad*ad);
0419   else if (ad<.6) r= 1-TMath::Exp(.917-4.279*ad-1.38*ad*ad);
0420   else r= 1-TMath::Exp(1.29-5.709*ad+.0186*ad*ad);
0421   cout<<"P:"<<r<<endl;
0422   return r;
0423 }
0424 
0425 float smartChi2(TH1F hard, TH1F soft){
0426   hard.Scale(1/hard.Integral());
0427   soft.Scale(1/soft.Integral());
0428   return soft.Chi2Test(&hard,"WW CHI2/NDF");
0429 }
0430 
0431 void chopHard(TH1F hard,TH1F soft){
0432   unsigned bins = soft.GetNbinsX();
0433   while(ADtoP(hard.AndersonDarlingTest(&soft))>.05&&bins>2){
0434     cout<<"Hard/Soft p="<<ADtoP(hard.AndersonDarlingTest(&soft))<<" n bins="<<hard.GetNbinsX()<<endl;
0435     cout<<"chi2 p="<<smartChi2(hard,soft)<<" n bins="<<hard.GetNbinsX()<<endl;
0436     hard=removeFirstBin(hard);
0437     soft=removeFirstBin(soft);
0438     bins--;
0439   }
0440   TH1F* temp = (TH1F*) hard.Clone("hard_chopped");
0441   temp->Write();
0442 }
0443 
0444 void calculateConversionRate(TEfficiency* rate, TH1F *pythia,TFile* out_file){
0445   TH1F* conversion_rate = (TH1F*)  pythia->Clone("rate");
0446   TH1* uni_rate = (TH1F*)rate->GetPassedHistogram()->Clone("uni_rate");
0447   uni_rate->Divide(rate->GetTotalHistogram());
0448   conversion_rate->Multiply(uni_rate);
0449   conversion_rate->Scale(1/pythia->Integral());
0450   out_file->Write();
0451 }
0452 
0453 double hardWeightFactor(TH1F* hard,TH1F* soft, unsigned matchBin, double* error){
0454   double hardTotal=hard->Integral(matchBin,hard->GetNbinsX());
0455   double softTotal=soft->Integral(matchBin,soft->GetNbinsX());
0456   double r= softTotal/hardTotal;
0457   *error = (softTotal+TMath::Power(softTotal,1./2))/hardTotal;
0458   cout<<soft->Integral(matchBin,soft->GetNbinsX())<<endl;
0459   cout<<hard->Integral(matchBin,soft->GetNbinsX())<<endl;
0460   return r;
0461 }
0462 
0463 TH1F* addSpec(TH1F* soft,TH1F* hard,TFile* file){
0464   TH1F* combined = new TH1F("combinedpythia","",soft->GetNbinsX(),soft->GetBinLowEdge(1),hard->GetBinLowEdge(hard->GetNbinsX()+1));
0465   unsigned matchBin = 10;//getMatchingBin(hard,soft); //hard coded by eye
0466   for (int i = 1; i < matchBin; ++i)
0467   {
0468     combined->SetBinContent(i,soft->GetBinContent(i));
0469     combined->SetBinError(i,soft->GetBinError(i));
0470   }
0471   double systematic;
0472   hard->Scale(hardWeightFactor(hard,soft,matchBin,&systematic));
0473   for (int i = matchBin; i < combined->GetNbinsX()+1; ++i)
0474   {
0475     combined->SetBinContent(i,hard->GetBinContent(i));
0476     combined->SetBinError(i,sqrt(hard->GetBinError(i)*hard->GetBinError(i)+systematic*systematic));
0477   }
0478   file->Write();
0479   return combined;
0480 } 
0481 
0482 void photonEff()
0483 {
0484   TFile *out_file = new TFile("effplots.root","UPDATE");
0485   //string treePath = "/sphenix/user/vassalli/gammasample/truthconversionembededanalysis";
0486   string treeExtension = ".root";
0487   unsigned int nFiles=1000;
0488   string softPath = "/sphenix/user/vassalli/minBiasPythia/softana.root";
0489   //string hard0Path = "/sphenix/user/vassalli/minBiasPythia/hard0ana.root";
0490   string hard4Path = "/sphenix/user/vassalli/minBiasPythia/hard4ana.root";
0491   TChain *softTree = new TChain("photonTree");
0492   //TChain *hard0Tree = new TChain("photonTree");
0493   TChain *hard4Tree = new TChain("photonTree");
0494   softTree->Add(softPath.c_str());
0495   //hard0Tree->Add(hard0Path.c_str());
0496   hard4Tree->Add(hard4Path.c_str());
0497   /*TChain *ttree = handleFile(treePath,treeExtension,"cutTreeSignal",nFiles);
0498   TChain *pairBackTree = handleFile(treePath,treeExtension,"pairBackTree",nFiles);
0499   TChain *vtxBackTree = handleFile(treePath,treeExtension,"vtxBackTree",nFiles);
0500   TChain *vertexingTree = handleFile(treePath,treeExtension,"vtxingTree",nFiles);
0501   makephotonM(ttree,vtxBackTree,out_file);*/
0502   auto pythiaSpec=makePythiaSpec(softTree,out_file,"soft");
0503   //makePythiaSpec(hard0Tree,out_file,"hard0");
0504   auto hardSpec = makePythiaSpec(hard4Tree,out_file,"hard4");
0505   //chopHard(*hardSpec,*pythiaSpec);
0506   pythiaSpec = addSpec(pythiaSpec,hardSpec,out_file);
0507   //makepTDist(ttree,observations,out_file);
0508   /*makeVtxRes(ttree,out_file);
0509   makeVtxEff(ttree,out_file);
0510   makepTRes(ttree);
0511   compareDeta(ttree,pairBackTree);
0512   //testCuts(ttree,out_file);
0513   makepTCaloGraph("pTcalodata.csv",out_file);
0514   makeVtxR(ttree,vertexingTree,out_file);*/
0515   //makeRefitDist(ttree,out_file);
0516 }