obsolete/ak/display.C

0001
0002 void display(const char *dfname, const char *cfname = 0)
0003 {
0004   auto *reco = new ReconstructionFactory(dfname, cfname, "pfRICH");
0005
0006   // [rad] (should match SPE sigma) & [ns];
0007   auto *a1 = reco->UseRadiator("Aerogel225", 0.0041);
0008   //reco->AddHypothesis("mu+");
0009   reco->AddHypothesis("pi+");
0010   //reco->AddHypothesis(321);
0011
0012   // Mark all pads hit by "calibration" (above the QE curve) photons as "useless";
0013   reco->AddBlackoutRadiator("QuartzWindow");
0014   reco->AddBlackoutRadiator("GasVolume");
0015   // Carelessly remove (0x1 << n)x(0x1 << n) square area "around" these hits;
0016   reco->SetBlackoutBlowupValue(3);
0017
0018   auto hxy = new TH2D("hxy", "", 650, -650., 650., 650, -650.0, 650.);
0019
0020   // This call is mandatory;
0021   reco->PerformCalibration();
0022   {
0023     CherenkovEvent *event;
0024
0025     // Event loop;
0026     while((event = reco->GetNextEvent())) {
0027       for(auto mcparticle: event->ChargedParticles()) {
0028     for(auto rhistory: mcparticle->GetRadiatorHistory()) {
0029       auto history  = mcparticle->GetHistory (rhistory);
0030
0031       for(auto photon: history->Photons()) {
0032         //printf("%d %d\n", photon->WasDetected(), photon->IsMarkedAsBlackout());
0033         if (!photon->WasDetected() ) continue;
0034         //if (!photon->IsUsefulForCalibration()) continue;
0035
0036         if (photon->IsMarkedAsBlackout()) continue;
0037
0038         TVector3 phx = photon->GetDetectionPosition();
0039         hxy->Fill(phx.X(), phx.Y());
0040       } //for photon
0041     } //for rhistory
0042       } //for mcparticle
0043     } //while
0044   }
0045
0046   gStyle->SetOptStat(0);
0047   auto cv = new TCanvas("cv", "", 1000, 1000);
0048   hxy->GetXaxis()->SetTitle("Sensor plane X, [mm]");
0049   hxy->GetYaxis()->SetTitle("Sensor plane Y, [mm]");
0050   hxy->Draw("COL");
0051 } // display()