Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-16 00:33:31

0001 #include "Validation/TrackingMCTruth/interface/TrackingTruthValid.h"
0002 #include "DataFormats/Common/interface/Handle.h"
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/Framework/interface/EventSetup.h"
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006 #include "FWCore/PluginManager/interface/ModuleDef.h"
0007 #include "FWCore/ServiceRegistry/interface/Service.h"
0008 #include "DQMServices/Core/interface/DQMStore.h"
0009 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
0010 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
0011 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0012 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0013 
0014 #include <cmath>
0015 
0016 typedef edm::RefVector<std::vector<TrackingParticle>> TrackingParticleContainer;
0017 
0018 typedef TrackingParticleRefVector::iterator tp_iterator;
0019 typedef TrackingParticle::g4t_iterator g4t_iterator;
0020 typedef TrackingParticle::genp_iterator genp_iterator;
0021 typedef TrackingVertex::genv_iterator genv_iterator;
0022 typedef TrackingVertex::g4v_iterator g4v_iterator;
0023 
0024 TrackingTruthValid::TrackingTruthValid(const edm::ParameterSet &conf)
0025     : runStandalone(conf.getParameter<bool>("runStandalone")),
0026       outputFile(conf.getParameter<std::string>("outputFile")),
0027       dbe_(nullptr),
0028       vec_TrackingParticle_Token_(consumes<TrackingParticleCollection>(conf.getParameter<edm::InputTag>("src"))) {}
0029 
0030 void TrackingTruthValid::bookHistograms(DQMStore::IBooker &ibooker, const edm::Run &run, const edm::EventSetup &es) {
0031   dbe_ = edm::Service<DQMStore>().operator->();
0032   ibooker.setCurrentFolder("Tracking/TrackingMCTruth/TrackingParticle");
0033 
0034   meTPMass = ibooker.book1D("TPMass", "Tracking Particle Mass", 100, -1, +5.);
0035   meTPCharge = ibooker.book1D("TPCharge", "Tracking Particle Charge", 10, -5, 5);
0036   meTPId = ibooker.book1D("TPId", "Tracking Particle Id", 500, -5000, 5000);
0037   meTPProc = ibooker.book1D("TPProc", "Tracking Particle Proc", 20, -0.5, 19.5);
0038   meTPAllHits = ibooker.book1D("TPAllHits", "Tracking Particle All Hits", 200, -0.5, 199.5);
0039   meTPMatchedHits = ibooker.book1D("TPMatchedHits", "Tracking Particle Matched Hits", 100, -0.5, 99.5);
0040   meTPPt = ibooker.book1D("TPPt", "Tracking Particle Pt", 100, 0, 100.);
0041   meTPEta = ibooker.book1D("TPEta", "Tracking Particle Eta", 100, -7., 7.);
0042   meTPPhi = ibooker.book1D("TPPhi", "Tracking Particle Phi", 100, -4., 4);
0043   meTPVtxX = ibooker.book1D("TPVtxX", "Tracking Particle VtxX", 100, -100, 100.);
0044   meTPVtxY = ibooker.book1D("TPVtxY", "Tracking Particle VtxY", 100, -100, 100.);
0045   meTPVtxZ = ibooker.book1D("TPVtxZ", "Tracking Particle VtxZ", 100, -100, 100.);
0046   meTPtip = ibooker.book1D("TPtip", "Tracking Particle tip", 100, 0, 1000.);
0047   meTPlip = ibooker.book1D("TPlip", "Tracking Particle lip", 100, 0, 100.);
0048 
0049   // Prepare Axes Labels for Processes
0050   meTPProc->setBinLabel(1, "Undefined");              // value =  0
0051   meTPProc->setBinLabel(2, "Unknown");                // value =  1
0052   meTPProc->setBinLabel(3, "Primary");                // value =  2
0053   meTPProc->setBinLabel(4, "Hadronic");               // value =  3
0054   meTPProc->setBinLabel(5, "Decay");                  // value =  4
0055   meTPProc->setBinLabel(6, "Compton");                // value =  5
0056   meTPProc->setBinLabel(7, "Annihilation");           // value =  6
0057   meTPProc->setBinLabel(8, "EIoni");                  // value =  7
0058   meTPProc->setBinLabel(9, "HIoni");                  // value =  8
0059   meTPProc->setBinLabel(10, "MuIoni");                // value =  9
0060   meTPProc->setBinLabel(11, "Photon");                // value = 10
0061   meTPProc->setBinLabel(12, "MuPairProd");            // value = 11
0062   meTPProc->setBinLabel(13, "Conversions");           // value = 12
0063   meTPProc->setBinLabel(14, "EBrem");                 // value = 13
0064   meTPProc->setBinLabel(15, "SynchrotronRadiation");  // value = 14
0065   meTPProc->setBinLabel(16, "MuBrem");                // value = 15
0066   meTPProc->setBinLabel(17, "MuNucl");                // value = 16
0067   meTPProc->setBinLabel(18, "");
0068   meTPProc->setBinLabel(19, "");
0069   meTPProc->setBinLabel(20, "");
0070 }
0071 
0072 void TrackingTruthValid::analyze(const edm::Event &event, const edm::EventSetup &c) {
0073   edm::Handle<TrackingParticleCollection> TruthTrackContainer;
0074   //  edm::Handle<TrackingVertexCollection>    TruthVertexContainer;
0075 
0076   event.getByToken(vec_TrackingParticle_Token_, TruthTrackContainer);
0077 
0078   const TrackingParticleCollection *tPC = TruthTrackContainer.product();
0079 
0080   // Loop over TrackingParticle's
0081   for (TrackingParticleCollection::const_iterator t = tPC->begin(); t != tPC->end(); ++t) {
0082     // if(t -> trackerPSimHit().size() ==0) cout << " Track with 0 SimHit " <<
0083     // endl;
0084 
0085     meTPMass->Fill(t->mass());
0086     meTPCharge->Fill(t->charge());
0087     meTPId->Fill(t->pdgId());
0088     meTPPt->Fill(sqrt(t->momentum().perp2()));
0089     meTPEta->Fill(t->momentum().eta());
0090     meTPPhi->Fill(t->momentum().Phi());
0091     //#warning "This file has been modified just to get it to compile without
0092     // any regard as to whether it still functions as intended" #ifdef
0093     // REMOVED_JUST_TO_GET_IT_TO_COMPILE__THIS_CODE_NEEDS_TO_BE_CHECKED
0094     //    std::vector<PSimHit> trackerPSimHit( t->trackPSimHit(DetId::Tracker)
0095     //    );
0096     //#endif
0097     meTPAllHits->Fill(t->numberOfTrackerHits());
0098     // get the process of the first hit
0099     //#warning "This file has been modified just to get it to compile without
0100     // any regard as to whether it still functions as intended" #ifdef
0101     // REMOVED_JUST_TO_GET_IT_TO_COMPILE__THIS_CODE_NEEDS_TO_BE_CHECKED
0102     //    if(trackerPSimHit.size() !=0) meTPProc->Fill(
0103     //    trackerPSimHit.front().processType());
0104     //#endif
0105 
0106     // there is no more the PSimHits collection !!! how to deal w/ the
0107     // processType ?
0108     //    if(t->numberOfTrackerHits() !=0) meTPProc->Fill(
0109     //    trackerPSimHit.front().processType());
0110 
0111     meTPMatchedHits->Fill(t->numberOfTrackerLayers());
0112     meTPVtxX->Fill(t->vx());
0113     meTPVtxY->Fill(t->vy());
0114     meTPVtxZ->Fill(t->vz());
0115     meTPtip->Fill(sqrt(t->vertex().perp2()));
0116     meTPlip->Fill(t->vz());
0117 
0118     /*
0119  // Compare momenta from sources
0120   cout << "T.P.   Track mass, Momentum, q , ID, & Event # "
0121         << t -> mass()  << " "
0122         << t -> p4()    << " " << t -> charge() << " "
0123         << t -> pdgId() << " "
0124         << t -> eventId().bunchCrossing() << "." << t -> eventId().event() <<
0125  endl;
0126 
0127   if(t->mass() < 0) cout << "======= WARNING, this particle has negative mass: "
0128  << t->mass()
0129                          << " and pdgId: " << t->pdgId() << endl;
0130   if(t->pdgId() == 0) cout << "======= WARNING, this particle has pdgId = 0: "
0131  << t->pdgId() << endl; cout << " Hits for this track: " << t ->
0132  trackerPSimHit().size() << endl;
0133   */
0134 
0135     /*
0136       std::cout << std::endl << "### Tracking Particle ###" << std::endl;
0137       std::cout << (*t) << std::endl;
0138       std::cout << "\t Tracker: " << t->trackerPSimHit().size() << std::endl;
0139       std::cout << "\t Muon: "    << t->muonPSimHit().size()    << std::endl;
0140       std::cout << (*t) << std::endl;
0141     */
0142   }  // End loop over TrackingParticle
0143 
0144   // Loop over TrackingVertex's
0145   /*
0146   cout << "Dumping sample vertex info" << endl;
0147   for (TrackingVertexCollection::const_iterator v = tVC -> begin(); v != tVC ->
0148   end(); ++v) { cout << " Vertex Position & Event #" << v -> position() << " "
0149   << v -> eventId().bunchCrossing() << "." << v -> eventId().event() << endl;
0150     cout << "  Associated with " << v -> daughterTracks().size() << " tracks" <<
0151   endl;
0152     // Get Geant and HepMC positions
0153     for (genv_iterator genV = v -> genVertices_begin(); genV != v ->
0154   genVertices_end(); ++genV) { cout << "  HepMC vertex position " <<
0155   (*(*genV)).position() << endl;
0156     }
0157     for (g4v_iterator g4V = v -> g4Vertices_begin(); g4V != v ->
0158   g4Vertices_end(); ++g4V) { cout << "  Geant vertex position " <<
0159   (*g4V).position() << endl;
0160       // Probably empty all the time, currently
0161     }
0162 
0163     // Loop over daughter track(s)
0164     for (tp_iterator iTP = v -> daughterTracks_begin(); iTP != v ->
0165   daughterTracks_end(); ++iTP) { cout << "  Daughter starts:      " <<
0166   (*(*iTP)).vertex(); for (g4t_iterator g4T  = (*(*iTP)).g4Track_begin(); g4T !=
0167   (*(*iTP)).g4Track_end(); ++g4T) { cout << " p " << g4T->momentum();
0168       }
0169       cout << endl;
0170     }
0171 
0172     // Loop over source track(s) (can be multiple since vertices are collapsed)
0173     for (tp_iterator iTP = v -> sourceTracks_begin(); iTP != v ->
0174   sourceTracks_end(); ++iTP) { cout << "  Source   starts: " <<
0175   (*(*iTP)).vertex(); for (g4t_iterator g4T  = (*iTP)->g4Track_begin(); g4T !=
0176   (*iTP)->g4Track_end(); ++g4T) { cout << ", p " <<  g4T ->momentum();
0177       }
0178       cout << endl;
0179     }
0180   }  // End loop over TrackingVertex
0181   */
0182 }