pfRICH/scripts/reco-ftbf.C

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