hpDIRC/time_imaging_reco/recoPdf_epic.cpp

0001 // Author: Nilanga Wickramaarachchi
0002 // Based on https://github.com/rdom/eicdirc - by Roman Dzhygadlo
0003
0004 #include "prttools.cpp"
0005 #include <TVirtualFitter.h>
0006 #include <TKey.h>
0007 #include <TRandom.h>
0008
0009 const int nch(24*256);
0010 TF1 *pdff[nch],*pdfs[nch];
0011 TH1F *hpdff[nch],*hpdfs[nch];
0012 double nhits_cut, nhits_cut_val;
0013
0014 void recoPdf_epic(TString in="eicrecon.root", TString pdf="eicrecon.pdf.root", int theta_ang = 30, double timeres=0.1, int pid =3, TString nameid="") // use pid = 0 for e+
0015 {
0016
0017   TGaxis::SetMaxDigits(4);
0018
0019   TCanvas *cc = new TCanvas("cc","cc");
0020
0021   TH1F *hl[5],*hll[5],*hnph[5];
0022   TH1F *Hist_leadtime = new TH1F("leadtime",";LE time [ns]; entries [#]", 2000, 0, 100);
0023   for(int i=0; i<5; i++){
0024     hl[i] = new TH1F(Form("hl_%d",i),";LE time [ns]; entries [#]", 2000,0,100);
0025     //hll[i]= new TH1F(Form("ll_i%d",i),";ln L("+prt_lname[pid]+") - ln L(#pi); entries [#]",240,-80,80);
0026     //hll[i]= new TH1F(Form("ll_i%d",i),";ln L("+prt_lname[pid]+") - ln L(#pi); entries [#]",300,-300,300);
0027     hll[i]= new TH1F(Form("ll_i%d",i),";ln L("+prt_lname[pid]+") - ln L(#pi); entries [#]",400,-200,200);
0028     hnph[i] = new TH1F(Form("hnph_%d",i),";multiplicity [#]; entries [#]", 220,0,220);
0029     hnph[i]->SetLineColor(prt_color[i]);
0030     hll[i]->SetLineColor(prt_color[i]);
0031   }
0032
0033   TFile f(pdf);
0034   TTree *tree_nph = (TTree*) f.Get("nph_pip");
0035   tree_nph->SetBranchAddress("nhits_cut",&nhits_cut);
0036   tree_nph->GetEntry(0);
0037   nhits_cut_val = nhits_cut;
0038   //std::cout << "number of hits cut at " << nhits_cut_val << std::endl;
0039
0040   int rebin=timeres/(100/2000.);
0041   std::cout<<"rebin "<<rebin <<std::endl;
0042
0043   int integ1(0), integ2(0);
0044   for(int i=0; i < nch; i++){
0045     hpdff[i] = (TH1F*)f.Get(Form("hf_%d",i));
0046     hpdfs[i] = (TH1F*)f.Get(Form("hs_%d",i));
0047     if(rebin >0) hpdff[i]->Rebin(rebin);
0048     if(rebin >0) hpdfs[i]->Rebin(rebin);
0049     integ1+= hpdff[i]->Integral();
0050     integ2+= hpdfs[i]->Integral();
0051     //hpdff[i]->Smooth(1);
0052     //hpdfs[i]->Smooth(1);
0053   }
0054
0055   prt_ch = new TChain("hpDIRCrawHits/dirctree");
0056   prt_ch->Add(in);
0057   int nEvents = prt_ch->GetEntries();
0058   std::cout << "Entries = " << nEvents << std::endl;
0059
0060   int hit_size = 0;
0061   int Particle_id = 0;
0062
0063   const int arr_size = 10000;
0064
0065   int mcp_num[arr_size], pixel_id[arr_size];
0066   Double_t lead_time[arr_size];
0067   double track_momentum_at_bar[arr_size][3];
0068   int track_id[arr_size];
0069
0070   prt_ch->SetBranchAddress("nhits", &hit_size);
0071   prt_ch->SetBranchAddress("mcp_id", &mcp_num);
0072   prt_ch->SetBranchAddress("pixel_id", &pixel_id);
0073   prt_ch->SetBranchAddress("hit_time", &lead_time);
0074
0075   int printstep=1000;
0076   double time = 0;
0077
0078   double nph[5]={0};
0079   double sigma_nph[5]={0};
0080   double nph_err[5] = {0};
0081   double sigma_nph_err[5] = {0};
0082   double sep_err;
0083
0084   int tnph(0),totalf(0),totals(0), ch(0);
0085
0086
0087   for (int ievent=0; ievent < nEvents; ievent++)
0088     {
0089       prt_ch->GetEntry(ievent);
0090
0091       if(ievent < nEvents/2) Particle_id = 211;
0092       else Particle_id = 321;
0093
0094       //if((Particle_id == 211 && ievent > 2499) || (Particle_id == 11 && ievent > 27499)) continue;
0095       if((Particle_id == 211 && ievent > 2499) || (Particle_id == 321 && ievent > 27499)) continue;
0096
0097       int nHits = hit_size;
0098
0099       //if(nHits < nhits_cut_val) continue;
0100       if(ievent%printstep==0 && ievent!=0) cout<< "Event # "<< ievent<< " # hits "<< nHits <<endl;
0101
0102
0103       int id = prt_get_pid(Particle_id);
0104       double aminf,amins, sum(0),sumf(0),sums(0);
0105       tnph = 0;
0106       if(hll[id]->GetEntries()>4000) continue;
0107
0108       for(int h=0; h < nHits; h++)
0109     {
0110       if(pixel_id[h] < 0) continue;
0111       ch = 256 * mcp_num[h] + pixel_id[h];
0112       time = lead_time[h] + gRandom->Gaus(0, timeres);
0113
0114       aminf = hpdff[ch]->GetBinContent(hpdff[ch]->FindBin(time));
0115       amins = hpdfs[ch]->GetBinContent(hpdfs[ch]->FindBin(time));
0116       tnph++;
0117
0118       //double noise = 1e-6;
0119       double noise = 0.5e-6;
0120
0121       sumf+=TMath::Log((aminf+noise));
0122       sums+=TMath::Log((amins+noise));
0123
0124       hl[id]->Fill(time);
0125       if(id == 2) Hist_leadtime->Fill(time);
0126     }
0127
0128       if(tnph>4) hnph[id]->Fill(tnph);
0129       sum = sumf-sums;
0130       if(fabs(sum)<0.1) continue;
0131
0132       hll[id]->Fill(sum);
0133     }
0134
0135   gStyle->SetOptStat(0);
0136
0137   prt_theta = theta_ang;
0138
0139   //TString name = Form("%1.0f_%1.2f_pik_6GeV_epic",prt_theta,timeres);
0140   TString name = "pik_6GeV_epic";
0141   //TString name = "pie_1.2GeV_epic";
0142
0143   prt_canvasAdd("nph_"+name,800,400);
0144   for(int i=0; i<5; i++)
0145     {
0146       if(hnph[i]->GetEntries()<20) continue;
0147       TFitResultPtr r = hnph[i]->Fit("gaus","SQ","",5,250);
0148       nph[i] = r->Parameter(1);
0149       nph_err[i] = r->ParError(1);
0150       sigma_nph[i] = r->Parameter(2);
0151       sigma_nph_err[i] = r->ParError(2);
0152     }
0153
0154   hnph[2]->Draw();
0155   hnph[pid]->Draw("sames hist");
0156
0157   prt_canvasGet("nph_"+name)->Update();
0158
0159   TPaveText *pt = new TPaveText();
0160   pt->AddText(Form("N(#pi^{+}) = %1.2f #pm %1.2f", nph[2], nph_err[2]));
0161   ((TText*)pt->GetListOfLines()->Last())->SetTextColor(prt_color[2]);
0162
0163   pt->AddText(Form("N(K^{+}) = %1.2f #pm %1.2f", nph[3], nph_err[3]));
0164   //pt->AddText(Form("N(e^{+}) = %1.2f #pm %1.2f", nph[0], nph_err[0]));
0165   ((TText*)pt->GetListOfLines()->Last())->SetTextColor(prt_color[3]);
0166   //((TText*)pt->GetListOfLines()->Last())->SetTextColor(prt_color[0]);
0167
0168   pt->SetX1NDC(0.12);
0169   pt->SetX2NDC(0.24);
0170   pt->SetY1NDC(0.65);
0171   pt->SetY2NDC(0.85);
0172   pt->SetShadowColor(0);
0173   pt->SetFillColor(0);
0174   pt->SetLineColor(0);
0175   pt->Draw();
0176
0177   prt_canvasAdd("sep_"+name,800,500);
0178   prt_normalize(hll,5);
0179
0180   TF1 *ff;
0181   double m1(0),m2(0),s1(100),s2(100);
0182   double dm1=0,dm2=0,ds1=0,ds2=0;
0183   if(hll[2]->GetEntries()>10){
0184     hll[2]->Fit("gaus","Sq");
0185     ff = hll[2]->GetFunction("gaus");
0186     ff->SetLineColor(1);
0187     m1=ff->GetParameter(1);
0188     s1=ff->GetParameter(2);
0189     dm1=ff->GetParError(1);
0190     ds1=ff->GetParError(2);
0191
0192     if(pid==0){ //handle tails
0193       hll[2]->Fit("gaus","S","",-300, m1+1.8*s1);
0194       ff = hll[2]->GetFunction("gaus");
0195       m1=ff->GetParameter(1);
0196       s1=ff->GetParameter(2);
0197       dm1=ff->GetParError(1);
0198       ds1=ff->GetParError(2);
0199     }
0200   }
0201
0202   if(hll[pid]->GetEntries()>10){
0203     hll[pid]->Fit("gaus","Sq");
0204     ff = hll[pid]->GetFunction("gaus");
0205     ff->SetLineColor(1);
0206     m2=ff->GetParameter(1);
0207     s2=ff->GetParameter(2);
0208     dm2=ff->GetParError(1);
0209     ds2=ff->GetParError(2);
0210
0211     if(pid==0){ ///handle tails
0212       hll[pid]->Fit("gaus","S","",m2-1.8*s2, 300);
0213       ff = hll[pid]->GetFunction("gaus");
0214       m2=ff->GetParameter(1);
0215       s2=ff->GetParameter(2);
0216       dm2=ff->GetParError(1);
0217       ds2=ff->GetParError(2);
0218     }
0219   }
0220
0221   double sep = (fabs(m1-m2))/(0.5*(s1+s2));
0222
0223   double e1,e2,e3,e4;
0224   e1=2/(s1+s2)*dm1;
0225   e2=-2/(s1+s2)*dm2;
0226   e3=-((2*(m1 - m2))/((s1 + s2)*(s1 + s2)))*ds1;
0227   e4=-((2*(m1 - m2))/((s1 + s2)*(s1 + s2)))*ds2;
0228   sep_err=sqrt(e1*e1+e2*e2+e3*e3+e4*e4);
0229
0230   /*std::cout << "mean1 = " << m1 << " +/- " << dm1 << std::endl;
0231   std::cout << "mean2 = " << m2 << " +/- " << dm2 <<std::endl;
0232   std::cout << "sigma1 = " << s1 << " +/- " << ds1 <<std::endl;
0233   std::cout << "sigma2 = " << s2 << " +/- " << ds2 <<std::endl;
0234   */
0235   std::cout<<in<<" separation "<< sep << "+/-" << sep_err << std::endl;
0236
0237   hll[2]->SetTitle(Form("#theta = %1.2f deg       separation = (%1.2f #pm %1.2f)#sigma",prt_theta, sep, sep_err));
0238   hll[2]->Draw();
0239   hll[pid]->Draw("same");
0240
0241   for(int i=0; i<5; i++){
0242     hl[i]->Scale(1/hl[i]->GetMaximum());
0243   }
0244
0245   prt_normalize(hl,3);
0246   prt_canvasAdd("tim_"+name,800,400);
0247   hl[2]->Draw();
0248   hl[pid]->SetLineColor(4);
0249   hl[2]->Draw("same");
0250   hl[pid]->SetLineColor(2);
0251   hl[pid]->Draw("same");
0252
0253   prt_canvasAdd("time_pip_"+name,800,400);
0254   gStyle->SetOptStat(1);
0255   Hist_leadtime->Draw();
0256
0257   prt_canvasSave("."+nameid,0);
0258
0259   TFile fc(in+"_r.root","recreate");
0260   TTree *tc = new TTree("reco","reco");
0261   tc->Branch("theta",&prt_theta,"theta/D");
0262   tc->Branch("sep",&sep,"sep/D");
0263   tc->Branch("sep_err",&sep_err,"sep_err/D");
0264   tc->Branch("timeres",&timeres,"timeres/D");
0265   tc->Branch("nph",&nph,"nph[5]/D");
0266   tc->Branch("nph_err",&nph_err,"nph_err[5]/D");
0267   tc->Branch("sigma_nph",&sigma_nph,"sigma_nph[5]/D");
0268   tc->Branch("sigma_nph_err",&sigma_nph_err,"sigma_nph_err[5]/D");
0269   tc->Fill();
0270   tc->Write();
0271 }