ana/work/plotpirej.C

0001 void plotpirej()
0002 {
0003 gStyle->SetOptStat(0);
0004
0005 const int nbins = 5;
0006 double bins[nbins+1];
0007 for(int i=0; i<=nbins; i++) { bins[i] = 2. + 2.*double(i); }
0008 //cout << "bin edges:" << endl;
0009 //for(int i=0; i<=nbins; i++) { cout <<  bins[i] << endl; }
0010
0011 int neopbins = 100;
0012 char hname[999];
0013 TH1D* heop_pi[nbins+1];
0014 TH1D* heop_e[nbins+1];
0015 for(int i=0; i<=nbins; i++) {
0016   sprintf(hname,"heop_pi_%d",i);
0017   heop_pi[i] = new TH1D(hname,hname,neopbins,0.,2.);
0018   heop_pi[i]->Sumw2();
0019   sprintf(hname,"heop_e_%d",i);
0020   heop_e[i] = new TH1D(hname,hname,neopbins,0.,2.);
0021   heop_e[i]->Sumw2();
0022 }
0023   TH1D* hhoe_e = new TH1D("hhoe_e","hhoe_e",150,0.,1.5);
0024   TH1D* hhoe_pi = new TH1D("hhoe_pi","hhoe_pi",150,0.,1.5);
0025   TH2D* hhhoe_e = new TH2D("hhhoe_e","hhhoe_e",100,0.,2.0,100,0.,2.0);
0026   TH2D* hhhoe_pi = new TH2D("hhhoe_pi","hhhoe_pi",100,0.,2.0,100,0.,2.0);
0027
0028   TH1D* hprob_e = new TH1D("hprob_e","hprob_e",100,0.,1.);
0029   TH1D* hprob_pi = new TH1D("hprob_pi","hprob_pi",100,0.,1.);
0030   TH1D* hchi2_e = new TH1D("hchi2_e","hchi2_e",100,0.,10.);
0031   TH1D* hchi2_pi = new TH1D("hchi2_pi","hchi2_pi",100,0.,10.);
0032
0033 TChain* ntp_track_pi = new TChain("ntppid");
0034   ntp_track_pi->AddFile("pions0.root");
0035   ntp_track_pi->AddFile("pions1.root");
0036   ntp_track_pi->AddFile("pions2.root");
0037   ntp_track_pi->AddFile("pions4.root");
0038
0039 //string infname_piminus;
0040 //ifstream ifs_piminus;
0041 //ifs_piminus.open("piminuslist.txt");
0042 //if ( ifs_piminus.fail() ) { cout << "FAIL TO READ INPUT FILE 1" << endl; ifs_piminus.close(); return; }
0043 //int tmpcount=0;
0044 //while(!ifs_piminus.eof()) {
0045 //  ifs_piminus >> infname_piminus;
0046 //  ntp_track_pi->AddFile(infname_piminus.c_str());
0047 //  tmpcount++;
0048 //  //if(tmpcount>20) break;
0049 //}
0050 //ifs_piminus.close();
0051
0052 TChain* ntp_track_e = new TChain("ntppid");
0053   ntp_track_e->AddFile("electrons0.root");
0054   ntp_track_e->AddFile("electrons1.root");
0055   ntp_track_e->AddFile("electrons2.root");
0056   ntp_track_e->AddFile("electrons3.root");
0057
0058 //string infname_e;
0059 //ifstream ifs_e;
0060 //ifs_e.open("electronlist.txt");
0061 //if ( ifs_e.fail() ) { cout << "FAIL TO READ INPUT FILE 2" << endl; ifs_e.close(); return; }
0062 //tmpcount=0;
0063 //while(!ifs_e.eof()) {
0064 //  ifs_e >> infname_e;
0065 //  ntp_track_e->AddFile(infname_e.c_str());
0066 //  tmpcount++;
0067 //  //if(tmpcount>20) break;
0068 //}
0069 //ifs_e.close();
0070
0071 cout << "number of entries for pions = " << ntp_track_pi->GetEntries() << endl;
0072 cout << "number of entries for electrons = " << ntp_track_e->GetEntries() << endl;
0073
0074 float px,py,pz,pt,quality,nmvtx,ntpc,eta,chi2,ndf,mom;
0075 float cemc_e,cemc_e3x3,cemc_dphi,cemc_deta,cemc_ecore,cemc_prob,cemc_chi2;
0076 float dca2d,hcalin_e3x3;
0077
0078 ntp_track_pi->SetBranchAddress("nmvtx",&nmvtx);
0079 ntp_track_pi->SetBranchAddress("ntpc",&ntpc);
0080 ntp_track_pi->SetBranchAddress("pt",&pt);
0081 ntp_track_pi->SetBranchAddress("eta",&eta);
0082 ntp_track_pi->SetBranchAddress("mom",&mom);
0083 ntp_track_pi->SetBranchAddress("quality",&quality);
0084 ntp_track_pi->SetBranchAddress("dca2d",&dca2d);
0085 ntp_track_pi->SetBranchAddress("cemc_ecore",&cemc_ecore);
0086 ntp_track_pi->SetBranchAddress("cemc_prob",&cemc_prob);
0087 ntp_track_pi->SetBranchAddress("cemc_chi2",&cemc_chi2);
0088 ntp_track_pi->SetBranchAddress("cemc_e3x3",&cemc_e3x3);
0089 ntp_track_pi->SetBranchAddress("cemc_dphi",&cemc_dphi);
0090 ntp_track_pi->SetBranchAddress("cemc_deta",&cemc_deta);
0091 ntp_track_pi->SetBranchAddress("hcalin_e3x3",&hcalin_e3x3);
0092
0093   for(int j=0; j<ntp_track_pi->GetEntries(); j++) {
0094     ntp_track_pi->GetEvent(j);
0095     if(quality>5) continue;
0096     if(nmvtx<1) continue;
0097     if(ntpc<20) continue;
0098     if(mom<=0.) continue;
0099     //if(fabs(cemc_dphi)>0.015) continue;
0100     //if(fabs(cemc_deta)>0.015) continue;
0101     if(cemc_ecore==0.) continue;
0102     hhoe_pi->Fill(hcalin_e3x3/cemc_ecore);
0103     float eop = cemc_ecore/mom;
0104     hhhoe_pi->Fill(eop,hcalin_e3x3/cemc_ecore);
0105     if(eop>=2.0) eop = 1.999;
0106     heop_pi[nbins]->Fill(eop);
0107     for(int j=0; j<nbins; j++) { if(pt>bins[j] && pt<bins[j+1]) { heop_pi[j]->Fill(eop); } }
0108     hprob_pi->Fill(cemc_prob);
0109     hchi2_pi->Fill(cemc_chi2);
0110   }
0111
0112 //----------------------------------------------------------------------------------------------
0113
0114 ntp_track_e->SetBranchAddress("nmvtx",&nmvtx);
0115 ntp_track_e->SetBranchAddress("ntpc",&ntpc);
0116 ntp_track_e->SetBranchAddress("pt",&pt);
0117 ntp_track_e->SetBranchAddress("eta",&eta);
0118 ntp_track_e->SetBranchAddress("mom",&mom);
0119 ntp_track_e->SetBranchAddress("quality",&quality);
0120 ntp_track_e->SetBranchAddress("dca2d",&dca2d);
0121 ntp_track_e->SetBranchAddress("cemc_ecore",&cemc_ecore);
0122 ntp_track_e->SetBranchAddress("cemc_prob",&cemc_prob);
0123 ntp_track_e->SetBranchAddress("cemc_chi2",&cemc_chi2);
0124 ntp_track_e->SetBranchAddress("cemc_e3x3",&cemc_e3x3);
0125 ntp_track_e->SetBranchAddress("cemc_dphi",&cemc_dphi);
0126 ntp_track_e->SetBranchAddress("cemc_deta",&cemc_deta);
0127 ntp_track_e->SetBranchAddress("hcalin_e3x3",&hcalin_e3x3);
0128
0129   for(int j=0; j<ntp_track_e->GetEntries(); j++) {
0130     ntp_track_e->GetEvent(j);
0131     if(quality>5) continue;
0132     if(nmvtx<1) continue;
0133     if(ntpc<20) continue;
0134     if(mom<=0.) continue;
0135     //if(fabs(cemc_dphi)>0.015) continue;
0136     //if(fabs(cemc_deta)>0.015) continue;
0137     if(cemc_ecore==0.) continue;
0138     hhoe_e->Fill(hcalin_e3x3/cemc_ecore);
0139     float eop = cemc_ecore/mom;
0140     hhhoe_e->Fill(eop,hcalin_e3x3/cemc_ecore);
0141     if(eop>=2.0) eop = 1.999;
0142     heop_e[nbins]->Fill(eop);
0143     for(int j=0; j<nbins; j++) { if(pt>bins[j] && pt<bins[j+1]) { heop_e[j]->Fill(eop); } }
0144     hprob_e->Fill(cemc_prob);
0145     hchi2_e->Fill(cemc_chi2);
0146   }
0147
0148 //--------------------------------------------------------------------------
0149
0150 TCanvas* choe = new TCanvas("choe","hcalin/cemc",50,50,600,600);
0151 choe->SetLogy();
0152 hhoe_e->SetLineColor(kRed);
0153 hhoe_pi->Draw();
0154 hhoe_e->Draw("same");
0155
0156 TCanvas* choe2 = new TCanvas("choe2","hcalin/cemc vs e/p",100,100,600,600);
0157 hhhoe_e->SetLineColor(kRed);
0158 hhhoe_pi->Draw("box");
0159 hhhoe_e->Draw("boxsame");
0160
0161 TCanvas* cprob = new TCanvas("cprob","prob",150,150,600,500);
0162 hprob_e->SetLineColor(kRed);
0163 hprob_pi->Draw();
0164 hprob_e->Draw("same");
0165
0166 TCanvas* cchi2 = new TCanvas("cchi2","cemc chi2",200,200,600,500);
0167 hchi2_e->SetLineColor(kRed);
0168 hchi2_pi->Draw();
0169 hchi2_e->Draw("same");
0170
0171 TCanvas* ceop = new TCanvas("ceop","E/p",250,250,600,500);
0172 ceop->SetLogy();
0173 heop_pi[nbins]->Draw();
0174 heop_e[nbins]->SetLineColor(kRed);
0175 heop_e[nbins]->Draw("same");
0176
0177 double vsume[nbins];
0178 double vsumpi[nbins];
0179 for(int i=0; i<nbins; i++) {vsume[i]=0.; vsumpi[i]=0.;}
0180 double sume = 0.;
0181 double sumpi = 0.;
0182 int uplim = heop_e[nbins]->GetNbinsX();
0183 for(int i=1; i<=uplim; i++) {
0184   sume += heop_e[nbins]->GetBinContent(i);
0185   sumpi += heop_pi[nbins]->GetBinContent(i);
0186 }
0187 cout << "total number of pions = " << sumpi << endl;
0188 cout << "total number of electrons = " << sume << endl;
0189
0190   for(int j=0; j<nbins; j++) {
0191     for(int i=1; i<=uplim; i++) {
0192       vsume[j] += heop_e[j]->GetBinContent(i);
0193       vsumpi[j] += heop_pi[j]->GetBinContent(i);
0194     }
0195     cout << "bin " << j << " " << vsume[j] << " " << vsumpi[j] << endl;
0196   }
0197
0198 double eleff = 0.9;
0199
0200 double vsumecut[nbins];
0201 int vbincut[nbins];
0202 for(int i = 0; i<nbins; i++) { vsumecut[i] = 0.; vbincut[i] = 0; }
0203 double sumecut = 0.;
0204 int bincut = 0;
0205 for(int i=uplim; i>=1; i--) {
0206   sumecut += heop_e[nbins]->GetBinContent(i);
0207   if(sumecut>sume*eleff) {
0208     cout << "90% cut: " << i << " " << heop_e[nbins]->GetBinLowEdge(i) << endl;
0209     bincut = i;
0210     break;
0211   }
0212 }
0213   for(int j=0; j<nbins; j++) {
0214     for(int i=uplim; i>=1; i--) {
0215     vsumecut[j] += heop_e[j]->GetBinContent(i);
0216       if(vsumecut[j]>vsume[j]*eleff) {
0217         cout << "bin " << j << " 90% cut: " << i << " " << heop_e[j]->GetBinLowEdge(i) << endl;
0218         vbincut[j] = i;
0219         break;
0220       }
0221     }
0222   }
0223
0224 double vsumpicut[nbins];
0225 for(int i=0; i<nbins; i++) {vsumpicut[i]=0.;}
0226 double sumpicut = 0.;
0227 for(int i=bincut; i<=heop_pi[nbins]->GetNbinsX(); i++) {
0228   sumpicut += heop_pi[nbins]->GetBinContent(i);
0229 }
0230 cout << "number of misidentified pions = " << sumpicut << endl;
0231 cout << "rejection factor = " << sumpi/sumpicut << " +- " << (sumpi/sumpicut)*(1./sqrt(sumpicut)) << endl;
0232
0233   for(int j=0; j<nbins; j++) {
0234     for(int i=vbincut[j]; i<=heop_pi[j]->GetNbinsX(); i++) {
0235       vsumpicut[j] += heop_pi[j]->GetBinContent(i);
0236     }
0237     cout << "bin " << j << " number of misidentified pions = " << vsumpicut[j] << endl;
0238   }
0239
0240 double rejpi[nbins];
0241 double errejpi[nbins];
0242 double xx[nbins];
0243 for(int i=0; i<nbins; i++) {
0244  rejpi[i] = vsumpicut[i]/vsumpi[i];
0245  errejpi[i] = rejpi[i]/sqrt(vsumpicut[i]);
0246  xx[i] = (bins[i]+bins[i+1])/2.;
0247 }
0248
0249 TCanvas* c4 = new TCanvas("c4","pion rejection",120,120,600,500);
0250 c4->SetLogy();
0251 double lowlim = 0.001;
0252 if(eleff==0.7) lowlim = 0.002;
0253 TH2D* hh = new TH2D("hh","",10,1.,13.,10,lowlim,0.1);
0254 hh->SetTitle(";Transverse momentum, [GeV] ; Inverse pion rejection factor");
0255 hh->Draw();
0256
0257 TGraphErrors* gr1 = new TGraphErrors(nbins,xx,rejpi,0,errejpi);
0258 gr1->SetMarkerStyle(20);
0259 gr1->SetMarkerColor(kBlack);
0260 gr1->SetMarkerSize(1.5);
0261 gr1->SetLineColor(kBlack);
0262 gr1->SetLineWidth(2);
0263 gr1->Draw("p");
0264
0265 }
0266