pfRICH/examples/multi-eval.C

0001
0002
0003 void multi_eval(const char *dfname, const char *cfname = 0)
0004 {
0005   auto *reco = new ReconstructionFactory(dfname, cfname, "pfRICH");
0006
0007   //
0008   // Factory configuration part;
0009   //
0010   //reco->IgnoreTimingInChiSquare();
0011   //reco->IgnorePoissonTermInChiSquare();
0012   //reco->SetSingleHitCCDFcut(0.005);
0013   //reco->RemoveAmbiguousHits();
0014   // Sensor active area pixelated will be pixellated NxN in digitization;
0015   //reco->SetSensorActiveAreaPixellation(24);
0016   auto *a1 = reco->UseRadiator("BelleIIAerogel1");
0017   //reco->SetSinglePhotonTimingResolution(0.030);
0018   //reco->SetQuietMode();
0019   reco->AddHypothesis("pi+");
0020   reco->AddHypothesis(321);
0021
0022   // Mark all pads hit by "calibration" (above the QE curve) photons as "useless";
0023   reco->AddBlackoutRadiator("QuartzWindow");
0024   reco->AddBlackoutRadiator("Acrylic");
0025   reco->AddBlackoutRadiator("GasVolume");
0026   // Carelessly remove (0x1 << n)x(0x1 << n) square area "around" these hits;
0027   reco->SetBlackoutBlowupValue(3);
0028
0029   auto hmatch = new TH1D("hmatch", "PID evaluation correctness",       2,    0,      2);
0030   auto hthtr1 = new TH1D("thtr1",  "Cherenkov angle (track)",        200,  220,    320);
0031   // For a dual aerogel configuration;
0032   //auto hthtr2  = new TH1D("thtr2",   "Cherenkov angle (track)",        200,  220,    320);
0033
0034   reco->UseAutomaticCalibration();
0035   // This call is mandatory; second argument: statistics (default: all events);
0036   reco->PerformCalibration(200);
0037   {
0038     CherenkovEvent *event;
0039
0040     // Event loop;
0041     while((event = reco->GetNextEvent())) {
0042       for(auto mcparticle: event->ChargedParticles()) {
0043     if (!mcparticle->IsPrimary()) continue;
0044
0045     if (mcparticle->GetPDG() == mcparticle->GetRecoPdgCode()) {
0046       hmatch->Fill(0.5);
0047     } else {
0048       hmatch->Fill(1.5);
0049     } //if
0050
0051     hthtr1->Fill(1000*mcparticle->GetRecoCherenkovAverageTheta(a1));
0052     //hthtr2->Fill(1000*mcparticle->GetRecoCherenkovAverageTheta(a2));
0053       } //for mcparticle
0054     } //while
0055   }
0056
0057   auto cv = new TCanvas("cv", "", 1600, 1000);
0058   cv->Divide(4, 3);
0059   cv->cd(1); reco->hthph()->Fit("gaus");
0060   cv->cd(2); reco->hccdfph()->SetMinimum(0); reco->hccdfph()->Draw();
0061   cv->cd(3); reco->hccdftr()->SetMinimum(0); reco->hccdftr()->Draw();
0062   cv->cd(4); reco->hccdfev()->SetMinimum(0); reco->hccdfev()->Draw();
0063   cv->cd(5); reco->hnpetr()->Draw();
0064   cv->cd(6); reco->hmatch()->SetMinimum(0); reco->hmatch()->Draw();
0065   cv->cd(7);       hmatch  ->SetMinimum(0);       hmatch  ->Draw();
0066   cv->cd(8); reco->hdtph()->Fit("gaus");
0067   cv->cd(9); hthtr1->Fit("gaus");
0068   //cv->cd(10); hthtr2->Fit("gaus");
0069   cv->cd(10); reco->hwl()->Draw();
0070   cv->cd(11); reco->hvtx()->Draw();
0071   cv->cd(12); reco->hri()->Draw();
0072 } // multi_eval()