hpDIRC/time_imaging_reco/draw_hp.C

0001 // Author: Nilanga Wickramaarachchi
0002 // Based on https://github.com/rdom/eicdirc - by Roman Dzhygadlo
0003
0004 #include <iostream>
0005 #include "prttools.cpp"
0006
0007 void draw_hp(TString infile)
0008 {
0009   TChain *fChain = new TChain("hpDIRCrawHits/dirctree");
0010   fChain->Add(infile);
0011
0012   int nEvents = fChain->GetEntries();
0013
0014   TGaxis::SetMaxDigits(3);
0015   prt_setRootPalette(1);
0016   prt_initDigi(2);
0017
0018   TH1F* h_pip_nph = new TH1F("pip_nph",";multiplicity [#]; entries [#]", 220,0,220);
0019
0020   const int arr_size = 5000;
0021
0022   int hit_size = 0;
0023
0024   fChain->SetBranchAddress("nhits", &hit_size);
0025
0026   int mcp_id[arr_size], pixel_id[arr_size];
0027
0028   fChain->SetBranchAddress("mcp_id", &mcp_id);
0029   fChain->SetBranchAddress("pixel_id", &pixel_id);
0030
0031   double mean_nph_pip, mean_nph_pip_err, sigma_nph_pip, sigma_nph_pip_err;
0032   double nhits_cut;
0033
0034   for (int ievent=0; ievent < nEvents/2; ievent++)
0035     {
0036       fChain->GetEntry(ievent);
0037
0038       int numHits = hit_size;
0039       int nph_pip = 0;
0040
0041       for(int i=0; i < numHits; i++)
0042     {
0043       nph_pip++;
0044     }
0045
0046       h_pip_nph->Fill(nph_pip);
0047     }
0048
0049   TFitResultPtr ptr = h_pip_nph->Fit("gaus","SQ","",0,220);
0050   mean_nph_pip = ptr->Parameter(1);
0051   mean_nph_pip_err = ptr->ParError(1);
0052   sigma_nph_pip = ptr->Parameter(2);
0053   sigma_nph_pip_err = ptr->ParError(2);
0054
0055   nhits_cut = mean_nph_pip - 3*sigma_nph_pip;
0056   double nhits_cut_val = nhits_cut;
0057
0058   //std::cout << "number of hits cut at " << nhits_cut << std::endl;
0059
0060   for (int ievent=0; ievent < 2500; ievent++){
0061     fChain->GetEntry(ievent);
0062     int nHits = hit_size;
0063     //if(nHits < nhits_cut) continue;
0064     if(ievent%1000==0) std::cout<<"Event # "<< ievent << " has "<< nHits <<" hits"<<std::endl;
0065
0066     for(int h=0; h < nHits; h++)
0067       {
0068     if(pixel_id[h] < 0) continue;
0069     int mcp = mcp_id[h];
0070     int pix = pixel_id[h];
0071
0072     prt_hdigi[mcp]->Fill(pix%16, pix/16);
0073       }
0074   }
0075   //------------------
0076
0077   { // hp
0078     auto cdigi = prt_drawDigi(2030);
0079     cdigi->SetName("hp");
0080     prt_canvasAdd(cdigi);
0081   }
0082
0083   prt_canvasSave(".", 0,0,0);
0084
0085 }