hpDIRC/time_imaging_reco/createPdf_epic.cpp

0001 // Author: Nilanga Wickramaarachchi
0002 // Based on https://github.com/rdom/eicdirc - by Roman Dzhygadlo
0003
0004 #include "prttools.cpp"
0005
0006 #include <iostream>
0007
0008 const int nch(24*256);
0009 TH1F *hlef[nch], *hles[nch];
0010
0011 void createPdf_epic(TString in="eicrecon.root", int pid=321) // use pid=11 for e+
0012 {
0013   for(int i=0; i<nch; i++){
0014     hlef[i] = new TH1F(Form("lef_%d",i),";LE time [ns]; entries/N [#]", 2000,0,100);
0015     hles[i] = new TH1F(Form("les_%d",i),";LE time [ns]; entries/N [#]", 2000,0,100);
0016   }
0017
0018   TH1F* h_pip_nph = new TH1F("pip_nph",";multiplicity [#]; entries [#]", 220,0,220);
0019
0020   prt_ch = new TChain("hpDIRCrawHits/dirctree");
0021   prt_ch->Add(in);
0022   int nEvents = prt_ch->GetEntries();
0023   std::cout << "Entries = " << nEvents << std::endl;
0024
0025   const int arr_size = 10000;
0026
0027   int hit_size = 0;
0028   int Particle_id = 0;
0029
0030   prt_ch->SetBranchAddress("nhits", &hit_size);
0031
0032   int mcp_num[arr_size], pixel_id[arr_size];
0033   Double_t lead_time[arr_size];
0034
0035   prt_ch->SetBranchAddress("mcp_id", &mcp_num);
0036   prt_ch->SetBranchAddress("pixel_id", &pixel_id);
0037   prt_ch->SetBranchAddress("hit_time", &lead_time);
0038
0039   int printstep = 2000;
0040   double time = 0;
0041
0042   int pdg(0), totalf(0),totals(0), ch(0);
0043   double mean_nph_pip, mean_nph_pip_err, sigma_nph_pip, sigma_nph_pip_err;
0044   double nhits_cut;
0045
0046   for (int ievent=0; ievent < nEvents; ievent++)
0047     {
0048       prt_ch->GetEntry(ievent);
0049
0050       if(ievent < nEvents/2) Particle_id = 211;
0051       else Particle_id = 321;
0052
0053       if(Particle_id != 211) continue;
0054       int numHits = hit_size;
0055       int nph_pip = 0;
0056
0057       for(int i=0; i < numHits; i++)
0058     {
0059       nph_pip++;
0060     }
0061
0062       h_pip_nph->Fill(nph_pip);
0063     }
0064
0065   prt_canvasAdd("nph_fit",800,400);
0066
0067   TFitResultPtr ptr = h_pip_nph->Fit("gaus","SQ","",0,220);
0068   mean_nph_pip = ptr->Parameter(1);
0069   mean_nph_pip_err = ptr->ParError(1);
0070   sigma_nph_pip = ptr->Parameter(2);
0071   sigma_nph_pip_err = ptr->ParError(2);
0072
0073   nhits_cut = mean_nph_pip - 3*sigma_nph_pip;
0074   double nhits_cut_val = nhits_cut; // can be used to reject events with low no. of hits (tail on the left side in photon yield)
0075
0076   //std::cout << "number of hits cut at " << nhits_cut << std::endl;
0077
0078   prt_canvasGet("nph_fit")->Update();
0079
0080   TPaveText *pt = new TPaveText();
0081   pt->AddText(Form("mean = %1.2f #pm %1.2f", mean_nph_pip, mean_nph_pip_err));
0082   pt->AddText(Form("#sigma = %1.2f #pm %1.2f", sigma_nph_pip, sigma_nph_pip_err));
0083   pt->AddText(Form("mean - 3#sigma = %1.2f", nhits_cut));
0084
0085   pt->SetX1NDC(0.12);
0086   pt->SetX2NDC(0.24);
0087   pt->SetY1NDC(0.65);
0088   pt->SetY2NDC(0.85);
0089   pt->SetShadowColor(0);
0090   pt->SetFillColor(0);
0091   pt->SetLineColor(0);
0092   pt->Draw();
0093
0094   prt_canvasSave(".",0);
0095
0096  for (int ievent=0; ievent < nEvents; ievent++)
0097    {
0098      prt_ch->GetEntry(ievent);
0099
0100      if(ievent < nEvents/2) Particle_id = 211;
0101      else Particle_id = 321;
0102      pdg = Particle_id;
0103
0104      //if((pdg == 211 && ievent < 2500) || (pdg == 11 && ievent < 27500)) continue;
0105      if((pdg == 211 && ievent < 2500) || (pdg == 321 && ievent < 27500)) continue;
0106
0107      int nHits = hit_size;
0108      if(ievent%printstep==0 && ievent!=0) std::cout << "Event # " << ievent << " # hits "<< nHits << std::endl;
0109
0110      if(nHits < 5) continue;
0111      //if(nHits < nhits_cut) continue;
0112      for(int i=0; i < nHits; i++)
0113       {
0114     if(pixel_id[i] < 0) continue;
0115     ch = mcp_num[i]*256 + pixel_id[i];
0116     time = lead_time[i] + gRandom->Gaus(0,0.1);
0117
0118     if(pdg==pid){
0119       totalf++;
0120       hlef[ch]->Fill(time);
0121     }
0122     if(pdg==211){
0123       totals++;
0124       hles[ch]->Fill(time);
0125     }
0126       }
0127    }
0128
0129  std::cout<<"#1 "<< totalf <<"  #2 "<<totals <<std::endl;
0130
0131  if(totalf>=0 && totals>0)
0132    {
0133     in.ReplaceAll(".root",".pdf.root");
0134     TFile efile(in,"RECREATE");
0135
0136     for(int i=0; i<nch; i++){
0137       hlef[i]->Scale(1/(double)totalf);
0138       hles[i]->Scale(1/(double)totals);
0139
0140       hlef[i]->SetName(Form("hf_%d",i));
0141       hles[i]->SetName(Form("hs_%d",i));
0142       hlef[i]->Write();
0143       hles[i]->Write();
0144
0145       if(0){
0146     int nrow=4, ncol=6, p = i/256;
0147     int np =p%ncol*nrow + p/ncol;
0148     hles[i]->SetName(Form("mcp%dpix%d",np,i%256));
0149     prt_canvasAdd(Form("pdf_mcp%dpix%d",np,i%256),800,400);
0150         prt_normalize(hlef[i],hles[i]);
0151         hlef[i]->SetLineColor(2);
0152     hles[i]->SetLineColor(4);
0153     hles[i]->GetXaxis()->SetRangeUser(5, 20);
0154         hles[i]->Draw("hist");
0155         hlef[i]->Draw("hist same");
0156     prt_canvasSave("data/pdfs_7",0,0,1);
0157       }
0158     }
0159
0160     TTree *tree_nph = new TTree("nph_pip","nph_pip");
0161     tree_nph->Branch("nhits_cut",&nhits_cut_val,"nhits_cut/D");
0162     tree_nph->Fill();
0163     tree_nph->Write();
0164
0165     efile.Write();
0166     efile.Close();
0167    }
0168
0169 }