irt/src/RICHIRT.cpp

0001 #include "RICHIRT.h"
0002 #include <stdexcept>
0003 ///////////////////////////
0004 RICHIRT::RICHIRT(){
0005   std::cout<<"RICH IRT Called\n"<<std::endl;
0006 }
0007 ///////////////////////////
0008 void RICHIRT::SetCollections(){
0009   m_IrtRadiatorInfo = get_irtRICH()+"IrtRadiatorInfo";
0010   m_IrtParticleCollection = get_irtRICH()+"IrtParticles";
0011   m_Tracks =  get_irtRICH()+"Tracks";
0012   m_MCParticles = "MCParticles";
0013   m_CentralCKFTracks="CentralCKFTracks";
0014 }
0015 /////////////////////
0016 void RICHIRT::Process(){
0017   auto reader = podio::makeReader(get_inFile());
0018   std::cout << "Events in file: " << reader.getEvents() << std::endl;
0019
0020   std::vector<int> mc_truth;
0021   for (size_t evt = 0; evt <reader.getEvents(); ++evt) {
0022     auto frame = reader.readEvent(evt);
0023     //std::cout << "\n=== Event " << evt << " ===" << std::endl;
0024
0025     const auto& mcParts = frame.get<edm4hep::MCParticleCollection>(m_MCParticles);
0026     // IRT Radiators, Particles and Track Segment
0027     const auto& irtRads = frame.get<edm4eic::IrtRadiatorInfoCollection>(m_IrtRadiatorInfo);
0028     const auto& irtParts = frame.get<edm4eic::IrtParticleCollection>(m_IrtParticleCollection);
0029     const auto& tracks = frame.get<edm4eic::TrackSegmentCollection>(m_Tracks);
0030     //Central track
0031     const auto& cCKFTracks=frame.get<edm4eic::TrackCollection>(m_CentralCKFTracks);
0032     //Dump MC Part PDG
0033     for(auto mcpart: mcParts){
0034       if(mcpart.getGeneratorStatus()!=1) continue;
0035       mc_truth.push_back(PDGToBin(mcpart.getPDG()));
0036     }
0037     //Get IRT-Tracks
0038     TVector3 p3;
0039     for(auto trs: tracks){
0040       //std::cout<<"DEBUG TRACK SEGMENT:: "<<trs.getLength()<<std::endl;
0041       auto tr = trs.getTrack();
0042       if(tr.isAvailable()){
0043         auto p = tr.getMomentum();
0044         p3 = TVector3 (p.x,p.y,p.z);
0045     m_RICH_Mom->Fill(p3.Mag());
0046     m_RICH_Eta->Fill(p3.PseudoRapidity());
0047         //std::cout<<"DEBUG: "<<p3.Mag()<<" "<<p3.PseudoRapidity()<<std::endl;
0048       }else std::cerr<<"DEBUG No Track "<<std::endl;
0049     }//IRT-Tracks
0050     int idx =0;
0051     //IRT Particles
0052     for(auto part: irtParts){
0053       int recPDG = part.getPDG();
0054       int recNp = part.getNpe();
0055       //if(recNp <=2 )std::cout<<recPDG<<" "<<recNp<<std::endl;
0056       //if(recPDG!= 211 ) std::cout<<"Info: "<<recPDG<<" "<<p3.Mag()<<" "<<p3.PseudoRapidity()<<std::endl;
0057       int rec = PDGToBin(part.getPDG());
0058       int truth = mc_truth.at(idx); idx++;
0059       FillConfusion(rec,truth);
0060
0061       ///// OneToMany : Radiator
0062       //std::cout<<m_Rad_Npe.size()<<"; "<<part.getRadiators().size()<<std::endl;
0063       if(m_Rad_Npe.size()!= part.getRadiators().size()){
0064     throw std::runtime_error("Radiator size (" + std::to_string(part.getRadiators().size())
0065                  + ") and histogram assignment (" + std::to_string(m_Rad_Npe.size())
0066                  + ") mismatch");
0067       }
0068       for (const auto& rad : part.getRadiators()) {
0069         auto id = rad.id();
0070         //if(part.getPDG() == 0) continue;
0071         m_Rad_Theta[id.index]->Fill(rad.getAngle());
0072         m_Rad_Npe[id.index]->Fill(rad.getNpe());
0073         m_Rad_NHits[id.index]->Fill(rad.getNhits());
0074     m_Theta_Mom[id.index]->Fill(p3.Mag(),rad.getAngle());
0075         //CKF_Angle->Fill(momCKF.Mag(),rad.getAngle());
0076         //MC_Angle->Fill(momMCP.Mag(),rad.getAngle());
0077       }//Radiators
0078     }//Irt-Particles
0079   }//events
0080   BookHistos();
0081 }
0082 void RICHIRT::BookHistos(){
0083   m_OutRootFile->cd();
0084   m_RICH_Mom->Write();
0085   m_RICH_Eta->Write();
0086   for(auto& h: m_Rad_Theta) h->Write();
0087   for(auto& h: m_Rad_Npe) h->Write();
0088   for(auto& h: m_Rad_NHits) h->Write();
0089   for(auto& h: m_Theta_Mom) h->Write();
0090   BookConfusionMatrix();
0091   m_OutRootFile->Write();
0092   m_OutRootFile->Close();
0093
0094 }
0095 RICHIRT::~RICHIRT(){
0096    std::cout<<"RICH IRT Destroyed\n"<<std::endl;
0097 }