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
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
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
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
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
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
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
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
0356
0357
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;
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
0486 string treeExtension = ".root";
0487 unsigned int nFiles=1000;
0488 string softPath = "/sphenix/user/vassalli/minBiasPythia/softana.root";
0489
0490 string hard4Path = "/sphenix/user/vassalli/minBiasPythia/hard4ana.root";
0491 TChain *softTree = new TChain("photonTree");
0492
0493 TChain *hard4Tree = new TChain("photonTree");
0494 softTree->Add(softPath.c_str());
0495
0496 hard4Tree->Add(hard4Path.c_str());
0497
0498
0499
0500
0501
0502 auto pythiaSpec=makePythiaSpec(softTree,out_file,"soft");
0503
0504 auto hardSpec = makePythiaSpec(hard4Tree,out_file,"hard4");
0505
0506 pythiaSpec = addSpec(pythiaSpec,hardSpec,out_file);
0507
0508
0509
0510
0511
0512
0513
0514
0515
0516 }