PhotonConversion/ana/minBiasRecoAna.C

0001 #include <iostream>
0002 using namespace std;
0003
0004 #include "TFile.h"
0005 #include "TTree.h"
0006 #include "TChain.h"
0007 #include "TLegend.h"
0008 #include "math.h"
0009 #include "TH1.h"
0010 #include "TH2.h"
0011 #include "TEfficiency.h"
0012
0013 TChain* handleFile(string name, string extension, string treename, unsigned int filecount){
0014   TChain *all = new TChain(treename.c_str());
0015   string temp;
0016   for (int i = 0; i < filecount; ++i)
0017   {
0018
0019     ostringstream s;
0020     s<<i;
0021     temp = name+string(s.str())+extension;
0022     all->Add(temp.c_str());
0023   }
0024   return all;
0025 }
0026
0027 int reportBackground(TChain* truthTree,TChain* recoTree,TFile* out_file){
0028   double signal = truthTree->GetEntries();
0029   signal*=(1-.012)*(1-.002)*(1-.011);
0030   int background = recoTree->GetEntries();
0031   background-=(int) signal;
0032   cout<<"For "<<truthTree->GetEntries()<<" truth conversions, I found "<<background<<" background events"<<endl;
0033   return (int) signal;
0034   /*float pT;
0035   float tpT;
0036   ttree->SetBranchAddress("photon_pT",&pT);
0037   ttree->SetBranchAddress("tphoton_pT",&tpT);
0038
0039   TH1F *pTeffPlot = new TH1F("#frac{#Delta#it{p}^{T}}{#it{p}_{#it{truth}}^{T}}","",40,-2,2);
0040   TH2F *pTefffuncPlot = new TH2F("pT_resolution_to_truthpt","",40,1,35,40,-1.5,1.5);
0041   pTeffPlot->Sumw2();
0042
0043   for (int event = 0; event < ttree->GetEntries(); ++event)
0044   {
0045     ttree->GetEvent(event);
0046     pTeffPlot->Fill((pT-tpT)/tpT);
0047     pTefffuncPlot->Fill(tpT,(pT-tpT)/tpT);
0048     }
0049     pTeffPlot->Scale(1./ttree->GetEntries(),"width");
0050     pTefffuncPlot->Scale(1./ttree->GetEntries(),"width");
0051     out_file->Write();*/
0052 }
0053
0054 void reportCuts(TChain* _treeBackground,int signal){
0055   unsigned totalTracks;
0056   unsigned passedpTEtaQ;
0057   unsigned passedCluster;
0058   unsigned passedPair;
0059   unsigned passedVtx;
0060   unsigned passedPhoton;
0061
0062   int sum_totalTracks=-1*signal;
0063   int sum_passedpTEtaQ=-1*signal;
0064   int sum_passedCluster=-1*signal;
0065   int sum_passedPair=-1*signal;
0066   int sum_passedVtx=-1*signal;
0067   int sum_passedPhoton=-1*signal;
0068
0069   _treeBackground->SetBranchAddress("total",   &totalTracks);
0070   _treeBackground->SetBranchAddress("tracksPEQ",  &passedpTEtaQ);
0071   _treeBackground->SetBranchAddress("tracks_clus", &passedCluster);
0072   _treeBackground->SetBranchAddress("pairs", &passedPair);
0073   _treeBackground->SetBranchAddress("vtx",    &passedVtx);
0074   _treeBackground->SetBranchAddress("photon",     &passedPhoton);
0075
0076   for(unsigned i=0; i<_treeBackground->GetEntries();i++){
0077     _treeBackground->GetEntry(i);
0078     sum_totalTracks+=totalTracks;
0079     sum_passedpTEtaQ+=passedpTEtaQ;
0080     sum_passedCluster+=passedCluster;
0081     sum_passedPair+=passedPair;
0082     sum_passedVtx+=passedVtx;
0083     sum_passedPhoton+=passedPhoton;
0084   }
0085   cout<<"Did RecoConversionEval with "<<sum_totalTracks<<" total tracks\n\t";
0086   cout<<1-(float)sum_passedpTEtaQ/sum_totalTracks<<"+/-"<<sqrt((float)sum_passedpTEtaQ)/sum_totalTracks<<" of tracks were rejected by pTEtaQ\n\t";
0087   cout<<1-(float)sum_passedCluster/sum_passedpTEtaQ<<"+/-"<<sqrt((float)sum_passedCluster)/sum_passedpTEtaQ<<" of remaining tracks were rejected by cluster\n\t";
0088   cout<<1-(float)sum_passedPair/sum_passedCluster<<"+/-"<<sqrt((float)sum_passedPair)/sum_passedCluster<<" of pairs were rejected by pair cuts\n\t";
0089   cout<<1-(float)sum_passedVtx/sum_passedPair<<"+/-"<<sqrt((float)sum_passedVtx)/sum_passedPair<<" of vtx were rejected by vtx cuts\n\t";
0090   cout<<1-(float)sum_passedPhoton/sum_passedVtx<<"+/-"<<sqrt((float)sum_passedPhoton)/sum_passedVtx<<" of photons were rejected by photon cuts\n";
0091   cout<<sum_passedPhoton<<" background events remain with "<<signal<<" signal events\n";
0092 }
0093
0094 int analyzeSignal(TChain* _signalCutTree,TFile *outFile){
0095   int _b_track_layer ;
0096   int _b_track_dlayer ;
0097   float _b_track_deta ;
0098   float _b_track_pT;
0099   float _b_ttrack_pT;
0100   float _b_vtx_radius ;
0101   float _b_tvtx_radius ;
0102   float _b_cluster_prob ;
0103   float _b_photon_pT ;
0104   float _b_tphoton_pT ;
0105   float _b_photon_m ;
0106   _signalCutTree->SetBranchAddress("track_deta", &_b_track_deta);
0107   _signalCutTree->SetBranchAddress("track_dlayer",&_b_track_dlayer);
0108   //    _signalCutTree->SetBranchAddress("track_layer", &_b_track_layer);
0109   _signalCutTree->SetBranchAddress("track_pT", &_b_track_pT);
0110   //_signalCutTree->SetBranchAddress("ttrack_pT", &_b_ttrack_pT);
0111    _signalCutTree->SetBranchAddress("vtx_radius", &_b_vtx_radius);
0112   //_signalCutTree->SetBranchAddress("tvtx_radius", &_b_tvtx_radius);
0113   // _signalCutTree->SetBranchAddress("vtx_chi2", &_b_vtx_chi2);
0114   _signalCutTree->SetBranchAddress("cluster_prob", &_b_cluster_prob);
0115   _signalCutTree->SetBranchAddress("photon_m", &_b_photon_m);
0116   _signalCutTree->SetBranchAddress("photon_pT", &_b_photon_pT);
0117   _signalCutTree->SetBranchAddress("tphoton_pT", &_b_tphoton_pT);
0118   unsigned pT=0, cluster=0, eta=0, photon=0,rsignal=0;
0119
0120   TH1F *rate_plot = new TH1F("rateplot","",60,0,15);
0121   rate_plot->Sumw2();
0122
0123   for(unsigned i=0; i<_signalCutTree->GetEntries();i++){
0124     _signalCutTree->GetEntry(i);
0125     if(_b_track_pT<.6) pT++;
0126     else if (_b_cluster_prob<0) cluster++;
0127     else if (_b_track_deta>.0082) eta++;
0128     else if (TMath::Abs(_b_track_dlayer)<=9&&_b_vtx_radius>1.84){
0129       if(_b_photon_m<.27||_b_photon_m>8.||_b_photon_pT<.039) photon++;
0130       else rsignal++;
0131     }
0132     rate_plot->Fill(_b_tphoton_pT);
0133   }
0134   cout<<"pT cut "<<pT<<" events\n";
0135   cout<<"cluster cut "<<cluster<<" events\n";
0136   cout<<"eta cut "<<eta<<" events\n";
0137   cout<<"photon cut "<<photon<<" events\n";
0138   rate_plot->Scale(1./6e5);
0139   outFile->Wrtie();
0140   return rsignal;
0141 }
0142
0143 void minBiasRecoAna()
0144 {
0145   string treePath = "/sphenix/user/vassalli/minBiasPythia/conversionembededminBiasanalysis";
0146   string truthTreePath = "/sphenix/user/vassalli/minBiasPythia/conversionembededminBiasTruthanalysis";
0147   string embedPath = "/sphenix/user/vassalli/RecoConversionTests/conversionembededonlineanalysis";
0148   string treeExtension = ".root";
0149   unsigned int nFiles=100;
0150   TFile *out_file = new TFile("minBiasplots.root","RECREATE");
0151   TChain *embed_truth_ttree = handleFile(embedPath,treeExtension,"cutTreeSignal",nFiles);
0152   TChain *embed_reco_ttree = handleFile(embedPath,treeExtension,"recoBackground",nFiles);
0153   TChain *cut_tree = new TChain("recoBackground");
0154   TChain *back_tree = new TChain("recoSignal");
0155   TChain *truth_ttree = new TChain("cutTreeSignal");
0156   string filename = "/sphenix/user/vassalli/minBiasConversion/conversionanaout.root";
0157   back_tree->Add(filename.c_str());
0158   cut_tree->Add(filename.c_str());
0159   filename = "/sphenix/user/vassalli/minBiasConversion/conversiontruthanaout.root";
0160   truth_ttree->Add(filename.c_str());
0161   reportCuts(cut_tree,analyzeSignal(truth_ttree));
0162 /*  cout<<"///////////////////////////////////////////////////////////\n";
0163   cout<<"EMBEDED ANALYSIS:\n";
0164   reportCuts(embed_reco_ttree,analyzeSignal(embed_truth_ttree));
0165 */
0166 }