obsolete/cc/multi-eval.C

0001
0002
0003 void multi_eval()
0004 {
0005   const char *dfname = "pfrich.root";
0006   const char *cfname =0;
0007   auto *reco = new ReconstructionFactory(dfname, cfname, "pfRICH");
0008
0009   //
0010   // Factory configuration part;
0011   //
0012   //reco->IgnoreTimingInChiSquare();
0013   //reco->IgnorePoissonTermInChiSquare();
0014   //reco->SetSingleHitCCDFcut(0.005);
0015   //reco->RemoveAmbiguousHits();
0016   // This only affects the calibration stage;
0017   //reco->SetDefaultSinglePhotonThetaResolution(0.0040);
0018   // Sensor active area pixelated will be pixellated NxN in digitization;
0019   //reco->SetSensorActiveAreaPixellation(24);
0020   // [rad] (should match SPE sigma) & [ns];
0021   //auto *a1 = reco->UseRadiator("Aerogel225",      0.0040);
0022   auto *a1 = reco->UseRadiator("BelleIIAerogel1");
0023   //reco->SetSinglePhotonTimingResolution(0.030);
0024   //reco->SetQuietMode();
0025   reco->AddHypothesis("pi+");
0026   reco->AddHypothesis(321);
0027
0028   // Mark all pads hit by "calibration" (above the QE curve) photons as "useless";
0029   reco->AddBlackoutRadiator("QuartzWindow");
0030   reco->AddBlackoutRadiator("Acrylic");
0031   reco->AddBlackoutRadiator("GasVolume");
0032   // Carelessly remove (0x1 << n)x(0x1 << n) square area "around" these hits;
0033   reco->SetBlackoutBlowupValue(3);
0034   TH2D* hthphi[5]; for(int i=0; i<5; i++){TString hName; hName.Form("theta_phi_%d",i);  hthphi[i]= new TH2D(hName,hName, 360,-180,180,200,270,310);}
0035   auto hmatch = new TH1D("hmatch", "PID evaluation correctness",       2,    0,      2);
0036   auto hthtr1 = new TH1D("thtr1",  "Cherenkov angle (track)",        400,  -100,    100);
0037   // For a dual aerogel configuration;
0038   //auto hthtr2  = new TH1D("thtr2",   "Cherenkov angle (track)",        200,  220,    320);
0039   auto htt   = new TH2D("htt","Cherenkov photon residual for particle thetas",5,151,175,100,-30,30);
0040   reco->UseAutomaticCalibration();
0041   // This call is mandatory; second argument: statistics (default: all events);
0042   reco->PerformCalibration(200);
0043   {
0044     CherenkovEvent *event;
0045
0046     // Event loop;
0047     while((event = reco->GetNextEvent())) {
0048       for(auto mcparticle: event->ChargedParticles()) {
0049     if (mcparticle->GetPDG() == mcparticle->GetRecoPdgCode()) {
0050       hmatch->Fill(0.5);
0051     } else {
0052       hmatch->Fill(1.5);
0053     } //if
0054         auto pp = mcparticle->GetVertexMomentum();
0055     double pmag = pp.Mag();
0056     double thhypo = 1000* (TMath::ACos(TMath::Sqrt(pmag*pmag+(0.139*0.139))/(1.04556*pmag)));    //costheta= 1/nbeta --> theta =arccos((p^2+m^2)/n.p)
0057     double theta = pp.Theta()*(180./M_PI);
0058         double eta = pp.PseudoRapidity();
0059     for(unsigned id =0; id<mcparticle->GetRecoCherenkovHitCount(); id++){
0060       double phth  = 1000*mcparticle->GetRecoCherenkovPhotonTheta(id);
0061           double phphi = (180./M_PI)*mcparticle->GetRecoCherenkovPhotonPhi(id);
0062        if(theta>151 && theta<=156)hthphi[0]->Fill(phphi,phth);
0063        if(theta>156 && theta<=161)hthphi[1]->Fill(phphi,phth);
0064        if(theta>161 && theta<=166)hthphi[2]->Fill(phphi,phth);
0065        if(theta>166 && theta<=171)hthphi[3]->Fill(phphi,phth);
0066        if(theta>171 && theta<=176)hthphi[4]->Fill(phphi,phth);
0067           //for(int ii=0; ii<100; ii++)htt->Fill(theta,reco->hthph()->GetBinCenter(ii),reco->hthph()->GetBinContent(ii));}
0068       //htt->Fill(theta,reco->GetThetaDiff());
0069       //if(theta>151. && theta<=156.) hthphi->Fill(phphi,phth);
0070     }// Cherenkov Hit Loop
0071     for(int ii=1; ii<=100; ii++){htt->Fill(theta,reco->hthph()->GetBinCenter(ii),reco->hthph()->GetBinContent(ii));}
0072     hthtr1->Fill(1000*mcparticle->GetRecoCherenkovAverageTheta(a1) - thhypo);
0073     //hthtr2->Fill(1000*mcparticle->GetRecoCherenkovAverageTheta(a2));
0074       } //for mcparticle
0075     } //while
0076   }
0077   vector<double> sigma; vector<double>mean;
0078   vector<double> thetabin; vector<double> binsize;
0079   vector<TH1D *> hh;
0080   /*
0081   for(int bin =1; bin<=htt->GetNbinsX()-1; bin++) {
0082     TString hName;
0083         hName.Form("h_%lf_%lf_py",((TAxis*)htt->GetXaxis())->GetBinCenter(bin),((TAxis*) htt->GetXaxis())->GetBinCenter(bin+1));
0084         thetabin.push_back(((TAxis*)htt->GetXaxis())->GetBinCenter(bin));
0085         binsize.push_back(((TAxis*)htt->GetXaxis())->GetBinCenter(bin+1)-((TAxis*)htt->GetXaxis())->GetBinCenter(bin));
0086     auto htmp= htt->ProjectionY(hName, bin, bin+1, "d");
0087     htmp->Fit("gaus","R","",-15,15);
0088     auto f1 =  htmp->GetFunction("gaus");
0089     sigma.push_back(f1->GetParameter(2));
0090     mean.push_back(f1->GetParameter(1));
0091     hh.push_back(htmp);
0092   }
0093   int size = sigma.size();
0094   auto GG = new TGraphErrors(size,&thetabin[0],&mean[0],&binsize[0],&sigma[0]);
0095   */
0096   auto cc = new TCanvas("cc","",1600,1000);
0097   //cc->Divide(3,1);
0098   //cc->cd(1);
0099   //hthphi->Draw("colz");
0100   //cc->cd(2);
0101   htt->Draw("COLZ");
0102   //cc->cd(3);
0103   //GG->Draw("AP *");
0104
0105   auto cc1 = new TCanvas("cc1","",1600,1000);
0106   cc1->Divide(3,2);
0107   for(int i=0; i<5; i++){cc1->cd(i+1); hthphi[i]->Draw("colz");}
0108
0109
0110   auto cv = new TCanvas("cv", "", 1600, 1000);
0111   cv->Divide(4, 3);
0112   cv->cd(1); reco->hthph()->Fit("gaus");
0113   cv->cd(2); reco->hccdfph()->SetMinimum(0); reco->hccdfph()->Draw();
0114   cv->cd(3); reco->hccdftr()->SetMinimum(0); reco->hccdftr()->Draw();
0115   cv->cd(4); reco->hccdfev()->SetMinimum(0); reco->hccdfev()->Draw();
0116   cv->cd(5); reco->hnpetr()->Draw();
0117   cv->cd(6); reco->hmatch()->SetMinimum(0); reco->hmatch()->Draw();
0118   cv->cd(7);       hmatch  ->SetMinimum(0);       hmatch  ->Draw();
0119   cv->cd(8); reco->hdtph()->Fit("gaus");
0120   cv->cd(9); hthtr1->Draw();
0121   //cv->cd(10); hthtr2->Fit("gaus");
0122   cv->cd(10); reco->hwl()->Draw();
0123   cv->cd(11); reco->hvtx()->Draw();
0124   cv->cd(12); reco->hri()->Draw();
0125   auto f = new TFile("out.root","RECREATE");
0126   reco->hthph()->Write();
0127   hthtr1->Write();
0128   reco->hnpetr()->Write();
0129   reco->hri()->Write();
0130   f->Close();
0131 } // multi_eval()