PhotonConversion/ana/observHister.C

0001 #include <iostream>
0002 using namespace std;
0003
0004 #include "TFile.h"
0005 #include "Scalar.h"
0006 #include "TTree.h"
0007 #include "TChain.h"
0008 #include "TLegend.h"
0009 #include "math.h"
0010 #include "TH1.h"
0011 #include "TH2.h"
0012 #include "TEfficiency.h"
0013 #include "TLine.h"
0014 #include "TGraphAsymmErrors.h"
0015
0016 TChain* handleFile(string name, string extension, string treename, unsigned int filecount){
0017   TChain *all = new TChain(treename.c_str());
0018   string temp;
0019   for (int i = 0; i < filecount; ++i)
0020   {
0021
0022     ostringstream s;
0023     s<<i;
0024     temp = name+string(s.str())+extension;
0025     all->Add(temp.c_str());
0026   }
0027   return all;
0028 }
0029
0030 void matchPlot(TChain* ttree,TFile* out_file){
0031   int b_match;
0032   int b_unmatch;
0033   ttree->SetBranchAddress("nMatched",     & b_match   );
0034   ttree->SetBranchAddress("nUnmatched",     &b_unmatch   );
0035   unsigned totalmatch=0;
0036   unsigned totaltrack=0;
0037   for (int event = 0; event < ttree->GetEntries(); ++event)
0038   {
0039     ttree->GetEvent(event);
0040
0041     totalmatch+=b_match;
0042     totaltrack+=b_match;
0043     totaltrack+=b_unmatch;
0044   }
0045   cout<<"Match "<<totalmatch<<" tracks\n";
0046   cout<<"Matched "<<(float)totalmatch/totaltrack*100.<<"\% of tracks"<<endl;
0047 }
0048
0049 void makeVtxR(TChain* ttree,TFile* out_file){
0050   float vtxr;
0051   float tvtxr;
0052   ttree->SetBranchAddress("vtx_radius",&vtxr);
0053   ttree->SetBranchAddress("tvtx_radius",&tvtxr);
0054
0055   std::vector<TH1F*> plots;
0056   plots.push_back(new TH1F("vtx_reco","",40,0,30));
0057   plots.push_back(new TH1F("vtx_truth","",40,0,30));
0058
0059   plots[0]->Sumw2();
0060   plots[1]->Sumw2();
0061
0062   double calc=0;
0063   for (int event = 0; event < ttree->GetEntries(); ++event)
0064   {
0065     ttree->GetEvent(event);
0066     plots[0]->Fill(vtxr);
0067     plots[1]->Fill(tvtxr);
0068     calc+=TMath::Abs(vtxr-tvtxr);
0069   }
0070   calc/=ttree->GetEntries();
0071   for (int i = 0; i < 2; ++i)
0072   {
0073     plots[i]->Scale(1./ttree->GetEntries(),"width");
0074   }
0075   out_file->Write();
0076   std::cout<<"mean deviation="<<calc<<std::endl;
0077 }
0078
0079 void makepTEff(TChain* ttree,TFile* out_file){
0080   float pT;
0081   float tpT;
0082   float track_pT;
0083   ttree->SetBranchAddress("photon_pT",&pT);
0084   ttree->SetBranchAddress("tphoton_pT",&tpT);
0085   ttree->SetBranchAddress("track_pT",&track_pT);
0086
0087   TH1F *pTeffPlot = new TH1F("#frac{#Delta#it{p}^{T}}{#it{p}_{#it{truth}}^{T}}","",40,-2,2);
0088   TH2F *pTefffuncPlot = new TH2F("pT_resolution_to_truthpt","",40,1,35,40,-1.5,1.5);
0089   TH1F *trackpTDist = new TH1F("truthpt","",40,0,35);
0090   pTeffPlot->Sumw2();
0091   pTefffuncPlot->Sumw2();
0092   trackpTDist->Sumw2();
0093   unsigned lowpTCount=0;
0094
0095   for (int event = 0; event < ttree->GetEntries(); ++event)
0096   {
0097     ttree->GetEvent(event);
0098     pTeffPlot->Fill((pT-tpT)/tpT);
0099     pTefffuncPlot->Fill(tpT,(pT-tpT)/tpT);
0100     trackpTDist->Fill(track_pT);
0101     if(tpT<.6)lowpTCount++;
0102   }
0103   pTeffPlot->Scale(1./ttree->GetEntries(),"width");
0104   pTefffuncPlot->Scale(1./ttree->GetEntries(),"width");
0105   trackpTDist->Scale(1./ttree->GetEntries(),"width");
0106   cout<<"Signal rejection through pT cut= "<<(float)lowpTCount/ttree->GetEntries()<<endl;
0107   cout<<"Signal rejection through pT cut= "<<sqrt((float)lowpTCount)/ttree->GetEntries()<<endl;
0108   out_file->Write();
0109 }
0110 void testCuts(TChain* ttree,TFile* out_file){
0111   float dphi;
0112   float prob;
0113   float track_pT;
0114   float deta;
0115   float radius;
0116   int layer;
0117   int dlayer;
0118   ttree->SetBranchAddress("cluster_dphi",&dphi);
0119   ttree->SetBranchAddress("cluster_prob",&prob);
0120   ttree->SetBranchAddress("track_layer",&layer);
0121   ttree->SetBranchAddress("track_dlayer",&dlayer);
0122   ttree->SetBranchAddress("track_pT",&track_pT);
0123   ttree->SetBranchAddress("track_deta",&deta);
0124   ttree->SetBranchAddress("vtx_radius",&radius);
0125
0126   TH1F *layerDist = new TH1F("layer","",16,-.5,15.5);
0127   TH1F *probDist = new TH1F("clust_prob","",30,-.5,1.);
0128   TH1F *deta_plot = new TH1F("deta","",30,-.001,.01);
0129   TH1F *dlayer_plot = new TH1F("dlayer","",11,-.5,10.5);
0130   TH1F *r_plot = new TH1F("signal_vtx_radius_dist","",26,-.5,25.5);
0131   layerDist->Sumw2();
0132   probDist->Sumw2();
0133   deta_plot->Sumw2();
0134   dlayer_plot->Sumw2();
0135   r_plot->Sumw2();
0136   unsigned badLayCount=0;
0137   unsigned badClusCount=0;
0138   unsigned bigDetaCount=0;
0139   unsigned shortRadiusCount=0;
0140
0141   for (int event = 0; event < ttree->GetEntries(); ++event)
0142   {
0143     ttree->GetEvent(event);
0144     if(layer==0)badLayCount++;
0145     if(dphi<0&&track_pT>.6){
0146       badClusCount++;
0147     }
0148     if(track_pT>.6&&dphi>=0){
0149       layerDist->Fill(layer);
0150       probDist->Fill(prob);
0151       deta_plot->Fill(deta);
0152       dlayer_plot->Fill(TMath::Abs(dlayer));
0153       if(deta>.0082||TMath::Abs(dlayer)>9)bigDetaCount++;
0154       else{
0155         r_plot->Fill(radius);
0156         if(radius<1.84) shortRadiusCount++;
0157       }
0158     }
0159
0160   }
0161   layerDist->Scale(1./ttree->GetEntries());
0162   deta_plot->Scale(1./ttree->GetEntries());
0163   dlayer_plot->Scale(1./ttree->GetEntries());
0164   r_plot->Scale(1./ttree->GetEntries());
0165   cout<<"Signal rejection through layer cut= "<<(float)badLayCount/ttree->GetEntries()<<endl;
0166   cout<<"error= "<<sqrt((float)badLayCount)/ttree->GetEntries()<<endl;
0167   cout<<"Signal rejection through clus cut= "<<(float)badClusCount/ttree->GetEntries()<<endl;
0168   cout<<"error= "<<sqrt((float)badClusCount)/ttree->GetEntries()<<endl;
0169   cout<<"Signal rejection through deta cut= "<<(float)bigDetaCount/ttree->GetEntries()<<endl;
0170   cout<<"error= "<<sqrt((float)bigDetaCount)/ttree->GetEntries()<<endl;
0171   cout<<"Signal rejection through radius cut= "<<(float)shortRadiusCount/ttree->GetEntries()<<endl;
0172   cout<<"error= "<<sqrt((float)shortRadiusCount)/ttree->GetEntries()<<endl;
0173   out_file->Write();
0174 }
0175 void makeRefitDist(TChain* ttree, TFile *out_file){
0176   float diffx;
0177   float diffy;
0178   float diffz;
0179   ttree->SetBranchAddress("refitdiffx",&diffx);
0180   ttree->SetBranchAddress("refitdiffy",&diffy);
0181   ttree->SetBranchAddress("refitdiffz",&diffz);
0182   TH1F *diffplotx= new TH1F("x_{0}-x_{refit}","",40,-10,10);
0183   TH1F *diffploty= new TH1F("y_{0}-y_{refit}","",40,-10,10);
0184   TH1F *diffplotz= new TH1F("z_{0}-z_{refit}","",40,-10,10);
0185
0186   diffplotx->Sumw2();
0187   diffploty->Sumw2();
0188   diffplotz->Sumw2();
0189
0190   for (int event = 0; event < ttree->GetEntries(); ++event)
0191   {
0192     ttree->GetEvent(event);
0193     diffplotx->Fill(diffx);
0194     diffploty->Fill(diffy);
0195     diffplotz->Fill(diffz);
0196   }
0197
0198   diffplotx->Scale(1./ttree->GetEntries(),"width");
0199   diffploty->Scale(1./ttree->GetEntries(),"width");
0200   diffplotz->Scale(1./ttree->GetEntries(),"width");
0201
0202   out_file->Write();
0203 }
0204
0205 void observHister()
0206 {
0207   string treePath = "/sphenix/user/vassalli/minBiasConversion/conversiontruthanaout.root";
0208   string treeExtension = ".root";
0209   unsigned int nFiles=100;
0210   TFile *out_file = new TFile("obsrvplots.root","RECREATE");
0211   TChain *ttree = new TChain("observTree");
0212   ttree->Add(treePath.c_str());
0213   cout<<"Total events= "<<ttree->GetEntries()<<'\n';
0214   matchPlot(ttree,out_file);
0215 }