obsolete/ak/tof-eval.C

0001
0002
0003 void tof_eval(const char *dfname, const char *cfname = 0)
0004 {
0005   auto *reco = new ReconstructionFactory(dfname, cfname, "pfRICH");
0006
0007   auto hnhit1d = new TH1D("nhit1d", "Digitized hit count",                       50,     0,     50);
0008   auto hnhit2d = new TH2D("nhit2d", "Digitized hit count",  50,    -2.5, -2.3,   50,     0,     50);
0009   auto htime1d = new TH1D("ntime1d","Timing",  50, 5.2, 6.2);
0010   auto htime2d = new TH2D("ntime2d","Timing", 100,    -2.5, -2.3, 50, 5.5, 6.5);
0011   auto hnpe    = new TH1D("npe",    "Photon count per hit",                      50,     0,     50);
0012
0013   reco->SetSinglePhotonTimingResolution(0.050);
0014   reco->SetExperimentalMode();
0015
0016   {
0017     CherenkovEvent *event;
0018
0019     std::map<std::pair<unsigned, unsigned>, unsigned> stat;
0020     std::multimap<unsigned, std::pair<unsigned, unsigned>> istat;
0021
0022     // Event loop;
0023     while((event = reco->GetNextEvent())) {
0024       //printf("%3d\n", (int)reco->GetDigitizedHits().size());
0025
0026       auto mcparticle = *event->ChargedParticles().begin();
0027       double eta = mcparticle->GetVertexMomentum().Eta();
0028       //printf("%f\n", eta);
0029
0030       unsigned tstat = 0;
0031       double tsum = 0.0;
0032       for(unsigned iq=0; iq<reco->GetDigitizedHits().size(); iq++) {
0033     auto &hit = reco->GetDigitizedHits()[iq];
0034     //printf("%2d %2d\n", hit.m_iX, hit.m_iY);
0035
0036     stat[std::make_pair(hit.m_iX, hit.m_iY)]++;
0037
0038     //printf("%f\n", hit.GetAverageDetectionTime());
0039     //if (/*eta > -2.35 &&*/ hit.m_iX == 20 && hit.m_iY == 3) {
0040     if ((hit.m_iY >= 5 && hit.m_iY <= 6) && (hit.m_iX >= 19 && hit.m_iX <= 23) /*&&
0041                         hit.GetPhotonCount() >= 5 && hit.GetPhotonCount() <= 15*/) {
0042       //for(unsigned ip=0; ip<hit.GetPhotonCount(); ip++)
0043       //printf("%2d %f\n", ip, hit.GetDetectionTime(ip));
0044
0045       hnpe ->Fill(hit.GetPhotonCount());
0046       if (hit.GetPhotonCount() >= 1) {
0047         //printf("%f\n", hit.GetAverageDetectionTime());
0048         //htime1d->Fill(hit.GetAverageDetectionTime());
0049         htime1d->Fill(hit.GetDetectionTime(0));
0050         //htime2d->Fill(eta, hit.GetAverageDetectionTime());
0051       }
0052
0053       tstat++;
0054       tsum += hit.GetDetectionTime(0);
0055       //tsum += hit.GetAverageDetectionTime();
0056     } //if
0057
0058       } //for iq
0059
0060       if (tstat) {
0061     tsum /= tstat;//reco->GetDigitizedHits().size();
0062     //printf("%f\n", tsum);
0063     htime1d->Fill(tsum);//hit.GetDetectionTime(0));
0064       }
0065
0066       hnhit1d->Fill(reco->GetDigitizedHits().size());
0067       hnhit2d->Fill(eta, reco->GetDigitizedHits().size());
0068     } //while
0069
0070     for(auto &entry: stat)
0071       istat.insert(std::make_pair(entry.second, entry.first));
0072     for(auto &entry: istat)
0073       printf("%3d -> %2d %2d\n", entry.first, entry.second.first, entry.second.second);
0074   }
0075
0076
0077   auto cv = new TCanvas("cv", "", 1600, 750);
0078   cv->Divide(4, 2);
0079   cv->cd(1); hnhit1d->Draw();
0080   cv->cd(2); hnhit2d->Draw("COLZ");
0081   cv->cd(3); htime1d->Fit("gaus");
0082   cv->cd(4); hnpe->Draw();
0083   cv->cd(5); htime2d->Draw("COLZ");
0084 } // tof_eval()