Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:32:02 

0001 /*\class to compute not containment parameter */
0002 
0003 #include "FWCore/Framework/interface/MakerMacros.h"
0004 #include "Validation/EcalClusters/interface/ContainmentCorrectionAnalyzer.h"
0005 
0006 using namespace cms;
0007 using namespace edm;
0008 using namespace std;
0009 using namespace reco;
0010 
0011 ContainmentCorrectionAnalyzer::ContainmentCorrectionAnalyzer(const ParameterSet &ps) : pTopologyToken(esConsumes()) {
0012   BarrelSuperClusterCollection_ =
0013       consumes<reco::SuperClusterCollection>(ps.getParameter<InputTag>("BarrelSuperClusterCollection"));
0014   EndcapSuperClusterCollection_ =
0015       consumes<reco::SuperClusterCollection>(ps.getParameter<InputTag>("EndcapSuperClusterCollection"));
0016   reducedBarrelRecHitCollection_ =
0017       consumes<EcalRecHitCollection>(ps.getParameter<InputTag>("reducedBarrelRecHitCollection"));
0018   reducedEndcapRecHitCollection_ =
0019       consumes<EcalRecHitCollection>(ps.getParameter<InputTag>("reducedEndcapRecHitCollection"));
0020   SimTrackCollection_ = consumes<SimTrackContainer>(ps.getParameter<InputTag>("simTrackCollection"));
0021   SimVertexCollection_ = consumes<SimVertexContainer>(ps.getParameter<InputTag>("simVertexCollection"));
0022 }
0023 
0024 ContainmentCorrectionAnalyzer::~ContainmentCorrectionAnalyzer() {}
0025 
0026 void ContainmentCorrectionAnalyzer::beginJob() {
0027   Service<TFileService> fs;
0028 
0029   // Define reference histograms
0030   h_EB_eRecoEtrueReference = fs->make<TH1F>("EB_eRecoEtrueReference", "EB_eRecoEtrueReference", 440, 0., 1.1);
0031   h_EB_e9EtrueReference = fs->make<TH1F>("EB_e9EtrueReference", "EB_e9EtrueReference", 440, 0., 1.1);
0032   h_EB_e25EtrueReference = fs->make<TH1F>("EB_e25EtrueReference", "EB_e25EtrueReference", 440, 0., 1.1);
0033   h_EE_eRecoEtrueReference = fs->make<TH1F>("EE_eRecoEtrueReference", "EE_eRecoEtrueReference", 440, 0., 1.1);
0034   h_EE_e9EtrueReference = fs->make<TH1F>("EE_e9EtrueReference", "EE_e9EtrueReference", 440, 0., 1.1);
0035   h_EE_e25EtrueReference = fs->make<TH1F>("EE_e25EtrueReference", "EE_e25EtrueReference", 440, 0., 1.1);
0036   h_EB_eTrue = fs->make<TH1F>("EB_eTrue", "EB_eTrue", 41, 40., 60.);
0037   h_EE_eTrue = fs->make<TH1F>("EE_eTrue", "EE_eTrue", 41, 40., 60.);
0038   h_EB_converted = fs->make<TH1F>("EB_converted", "EB_converted", 2, 0., 2.);
0039   h_EE_converted = fs->make<TH1F>("EE_converted", "EE_converted", 2, 0., 2.);
0040 }
0041 
0042 void ContainmentCorrectionAnalyzer::analyze(const Event &evt, const EventSetup &es) {
0043   LogInfo("ContainmentCorrectionAnalyzer") << "Analyzing event " << evt.id() << "\n";
0044 
0045   // taking the needed collections
0046   std::vector<SimTrack> theSimTracks;
0047   Handle<SimTrackContainer> SimTk;
0048   evt.getByToken(SimTrackCollection_, SimTk);
0049   Labels l;
0050   labelsForToken(SimTrackCollection_, l);
0051 
0052   if (SimTk.isValid())
0053     theSimTracks.insert(theSimTracks.end(), SimTk->begin(), SimTk->end());
0054   else {
0055     LogError("ContainmentCorrectionAnalyzer") << "Error! can't get collection with label " << l.module;
0056   }
0057 
0058   std::vector<SimVertex> theSimVertexes;
0059   Handle<SimVertexContainer> SimVtx;
0060   evt.getByToken(SimVertexCollection_, SimVtx);
0061   labelsForToken(SimVertexCollection_, l);
0062 
0063   if (SimVtx.isValid())
0064     theSimVertexes.insert(theSimVertexes.end(), SimVtx->begin(), SimVtx->end());
0065   else {
0066     LogError("ContainmentCorrectionAnalyzer") << "Error! can't get collection with label " << l.module;
0067   }
0068 
0069   const reco::SuperClusterCollection *BarrelSuperClusters = nullptr;
0070   Handle<reco::SuperClusterCollection> pHybridBarrelSuperClusters;
0071   evt.getByToken(BarrelSuperClusterCollection_, pHybridBarrelSuperClusters);
0072   labelsForToken(BarrelSuperClusterCollection_, l);
0073 
0074   if (pHybridBarrelSuperClusters.isValid()) {
0075     BarrelSuperClusters = pHybridBarrelSuperClusters.product();
0076   } else {
0077     LogError("ContainmentCorrectionAnalyzer") << "Error! can't get collection with label " << l.module;
0078   }
0079 
0080   const reco::SuperClusterCollection *EndcapSuperClusters = nullptr;
0081   Handle<reco::SuperClusterCollection> pMulti5x5EndcapSuperClusters;
0082   evt.getByToken(EndcapSuperClusterCollection_, pMulti5x5EndcapSuperClusters);
0083   labelsForToken(EndcapSuperClusterCollection_, l);
0084 
0085   if (pMulti5x5EndcapSuperClusters.isValid())
0086     EndcapSuperClusters = pMulti5x5EndcapSuperClusters.product();
0087   else {
0088     LogError("ContainmentCorrectionAnalyzer") << "Error! can't get collection with label " << l.module;
0089   }
0090 
0091   const EcalRecHitCollection *ebRecHits = nullptr;
0092   Handle<EcalRecHitCollection> pEBRecHits;
0093   evt.getByToken(reducedBarrelRecHitCollection_, pEBRecHits);
0094   labelsForToken(reducedBarrelRecHitCollection_, l);
0095 
0096   if (pEBRecHits.isValid())
0097     ebRecHits = pEBRecHits.product();
0098   else {
0099     LogError("ContainmentCorrectionAnalyzer") << "Error! can't get collection with label " << l.module;
0100   }
0101 
0102   const EcalRecHitCollection *eeRecHits = nullptr;
0103   Handle<EcalRecHitCollection> pEERecHits;
0104   evt.getByToken(reducedEndcapRecHitCollection_, pEERecHits);
0105   labelsForToken(reducedEndcapRecHitCollection_, l);
0106 
0107   if (pEERecHits.isValid())
0108     eeRecHits = pEERecHits.product();
0109   else {
0110     LogError("ContainmentCorrectionAnalyzer") << "Error! can't get collection with label " << l.module;
0111   }
0112 
0113   const CaloTopology *topology = nullptr;
0114   auto pTopology = es.getHandle(pTopologyToken);
0115   if (pTopology.isValid())
0116     topology = &es.getData(pTopologyToken);
0117 
0118   std::vector<EcalSimPhotonMCTruth> photons = findMcTruth(theSimTracks, theSimVertexes);
0119 
0120   nMCphotons = 0;
0121   nRECOphotons = 0;
0122 
0123   int mc_size = photons.size();
0124   mcEnergy.resize(mc_size);
0125   mcEta.resize(mc_size);
0126   mcPhi.resize(mc_size);
0127   mcPt.resize(mc_size);
0128   isConverted.resize(mc_size);
0129   x_vtx.resize(mc_size);
0130   y_vtx.resize(mc_size);
0131   z_vtx.resize(mc_size);
0132 
0133   superClusterEnergy.resize(mc_size);
0134   superClusterEta.resize(mc_size);
0135   superClusterPhi.resize(mc_size);
0136   superClusterEt.resize(mc_size);
0137   e1.resize(mc_size);
0138   e9.resize(mc_size);
0139   e25.resize(mc_size);
0140   seedXtal.resize(mc_size);
0141   iMC.resize(mc_size);
0142 
0143   // loop over MC truth photons
0144   for (unsigned int ipho = 0; ipho < photons.size(); ipho++) {
0145     math::XYZTLorentzVectorD vtx = photons[ipho].primaryVertex();
0146     double phiTrue = photons[ipho].fourMomentum().phi();
0147     double vtxPerp = sqrt(vtx.x() * vtx.x() + vtx.y() * vtx.y());
0148     double etaTrue = ecalEta(photons[ipho].fourMomentum().eta(), vtx.z(), vtxPerp);
0149     double etTrue = photons[ipho].fourMomentum().e() / cosh(etaTrue);
0150     nMCphotons++;
0151     mcEnergy[nMCphotons - 1] = photons[ipho].fourMomentum().e();
0152     mcEta[nMCphotons - 1] = etaTrue;
0153     mcPhi[nMCphotons - 1] = phiTrue;
0154     mcPt[nMCphotons - 1] = etTrue;
0155     isConverted[nMCphotons - 1] = photons[ipho].isAConversion();
0156     x_vtx[nMCphotons - 1] = vtx.x();
0157     y_vtx[nMCphotons - 1] = vtx.y();
0158     z_vtx[nMCphotons - 1] = vtx.z();
0159 
0160     // check histos for MC truth
0161     if (std::fabs(etaTrue) < 1.479) {
0162       h_EB_eTrue->Fill(photons[ipho].fourMomentum().e());
0163       h_EB_converted->Fill(photons[ipho].isAConversion());
0164     }
0165     if (std::fabs(etaTrue) >= 1.6) {
0166       h_EE_eTrue->Fill(photons[ipho].fourMomentum().e());
0167       h_EE_converted->Fill(photons[ipho].isAConversion());
0168     }
0169 
0170     // barrel
0171     if (std::fabs(etaTrue) < 1.479) {
0172       double etaCurrent;  // , etaFound = 0; // UNUSED
0173       double phiCurrent;
0174       // double etCurrent,  etFound  = 0; // UNUSED
0175       const reco::SuperCluster *nearSC = nullptr;
0176 
0177       double closestParticleDistance = 999;
0178       for (reco::SuperClusterCollection::const_iterator aClus = BarrelSuperClusters->begin();
0179            aClus != BarrelSuperClusters->end();
0180            aClus++) {
0181         etaCurrent = aClus->position().eta();
0182         phiCurrent = aClus->position().phi();
0183         // etCurrent  = aClus->energy()/std::cosh(etaCurrent); // UNUSED
0184         double deltaR = std::sqrt(std::pow(etaCurrent - etaTrue, 2) + std::pow(phiCurrent - phiTrue, 2));
0185         if (deltaR < closestParticleDistance) {
0186           // etFound  = etCurrent; // UNUSED
0187           // etaFound = etaCurrent; // UNUSED
0188           closestParticleDistance = deltaR;
0189           nearSC = &(*aClus);
0190         }
0191       }
0192 
0193       if (closestParticleDistance < 0.3) {
0194         // Is a matched particle dumping informations
0195         nRECOphotons++;
0196         superClusterEnergy[nRECOphotons - 1] = nearSC->rawEnergy();
0197         superClusterEta[nRECOphotons - 1] = nearSC->position().eta();
0198         superClusterPhi[nRECOphotons - 1] = nearSC->position().phi();
0199         superClusterEt[nRECOphotons - 1] = nearSC->rawEnergy() / std::cosh(nearSC->position().eta());
0200         iMC[nRECOphotons - 1] = nMCphotons - 1;
0201 
0202         const reco::CaloClusterPtr &theSeed = nearSC->seed();
0203         seedXtal[nRECOphotons - 1] = EcalClusterTools::getMaximum(*theSeed, ebRecHits).first;
0204         e1[nRECOphotons - 1] = EcalClusterTools::eMax(*theSeed, ebRecHits);
0205         e9[nRECOphotons - 1] = EcalClusterTools::e3x3(*theSeed, ebRecHits, topology);
0206         e25[nRECOphotons - 1] = EcalClusterTools::e5x5(*theSeed, ebRecHits, topology);
0207       }
0208     }
0209 
0210     // endcap
0211     if (std::fabs(etaTrue) >= 1.6) {
0212       double etaCurrent;  // , etaFound = 0; // UNUSED
0213       double phiCurrent;
0214       // double etCurrent,  etFound  = 0; // UNUSED
0215       const reco::SuperCluster *nearSC = nullptr;
0216 
0217       double closestParticleDistance = 999;
0218       for (reco::SuperClusterCollection::const_iterator aClus = EndcapSuperClusters->begin();
0219            aClus != EndcapSuperClusters->end();
0220            aClus++) {
0221         etaCurrent = aClus->position().eta();
0222         phiCurrent = aClus->position().phi();
0223         // etCurrent  =  aClus->energy()/std::cosh(etaCurrent);
0224         double deltaR = std::sqrt(std::pow(etaCurrent - etaTrue, 2) + std::pow(phiCurrent - phiTrue, 2));
0225         if (deltaR < closestParticleDistance) {
0226           // etFound  = etCurrent; // UNUSED
0227           // etaFound = etaCurrent; // UNUSED
0228           closestParticleDistance = deltaR;
0229           nearSC = &(*aClus);
0230         }
0231       }
0232 
0233       if (closestParticleDistance < 0.3) {
0234         // Is a matched particle dumping informations
0235         nRECOphotons++;
0236         float psEnergy = nearSC->preshowerEnergy();
0237         superClusterEnergy[nRECOphotons - 1] = nearSC->rawEnergy() + psEnergy;
0238         superClusterEta[nRECOphotons - 1] = nearSC->position().eta();
0239         superClusterPhi[nRECOphotons - 1] = nearSC->position().phi();
0240         superClusterEt[nRECOphotons - 1] = (nearSC->rawEnergy() + psEnergy) / std::cosh(nearSC->position().eta());
0241         iMC[nRECOphotons - 1] = nMCphotons - 1;
0242 
0243         const reco::CaloClusterPtr &theSeed = nearSC->seed();
0244         seedXtal[nRECOphotons - 1] = EcalClusterTools::getMaximum(*theSeed, eeRecHits).first;
0245         e1[nRECOphotons - 1] = EcalClusterTools::eMax(*theSeed, eeRecHits) + psEnergy;
0246         e9[nRECOphotons - 1] = EcalClusterTools::e3x3(*theSeed, eeRecHits, topology) + psEnergy;
0247         e25[nRECOphotons - 1] = EcalClusterTools::e5x5(*theSeed, eeRecHits, topology) + psEnergy;
0248       }
0249     }
0250   }
0251 
0252   // containment analysis for unconverted photons in the reference region only
0253   for (int i = 0; i < nRECOphotons; i++) {
0254     // barrel
0255     if (fabs(superClusterEta[i]) < 1.479) {
0256       if (isConverted[iMC[i]] != 1) {
0257         int ietaAbs = (seedXtal[i] >> 9) & 0x7F;
0258         int iphi = seedXtal[i] & 0x1FF;
0259         if (ietaAbs > 5 && ietaAbs < 21 && ((iphi % 20) > 5) && ((iphi % 20) < 16)) {
0260           h_EB_eRecoEtrueReference->Fill(superClusterEnergy[i] / mcEnergy[iMC[i]]);
0261           h_EB_e9EtrueReference->Fill(e9[i] / mcEnergy[iMC[i]]);
0262           h_EB_e25EtrueReference->Fill(e25[i] / mcEnergy[iMC[i]]);
0263         }
0264       }
0265     }
0266 
0267     // endcap
0268     if (fabs(superClusterEta[i]) > 1.6) {
0269       if (isConverted[iMC[i]] != 1) {
0270         if (fabs(superClusterEta[i]) > 1.7 && fabs(superClusterEta[i] < 2.3) &&
0271             ((superClusterPhi[i] > -CLHEP::pi / 2. + 0.1 && superClusterPhi[i] < CLHEP::pi / 2. - 0.1) ||
0272              (superClusterPhi[i] > CLHEP::pi / 2. + 0.1) || (superClusterPhi[i] < -CLHEP::pi / 2. - 0.1))) {
0273           h_EE_eRecoEtrueReference->Fill(superClusterEnergy[i] / mcEnergy[iMC[i]]);
0274           h_EE_e9EtrueReference->Fill(e9[i] / mcEnergy[iMC[i]]);
0275           h_EE_e25EtrueReference->Fill(e25[i] / mcEnergy[iMC[i]]);
0276         }
0277       }
0278     }
0279 
0280   }  // loop over reco photons
0281 }
0282 
0283 void ContainmentCorrectionAnalyzer::endJob() {}
0284 
0285 float ContainmentCorrectionAnalyzer::ecalEta(float EtaParticle, float Zvertex, float plane_Radius) {
0286   const float R_ECAL = 136.5;
0287   const float Z_Endcap = 328.0;
0288   const float etaBarrelEndcap = 1.479;
0289 
0290   if (EtaParticle != 0.) {
0291     float Theta = 0.0;
0292     float ZEcal = (R_ECAL - plane_Radius) * sinh(EtaParticle) + Zvertex;
0293 
0294     if (ZEcal != 0.0)
0295       Theta = atan(R_ECAL / ZEcal);
0296     if (Theta < 0.0)
0297       Theta = Theta + Geom::pi();
0298 
0299     float ETA = -log(tan(0.5 * Theta));
0300 
0301     if (fabs(ETA) > etaBarrelEndcap) {
0302       float Zend = Z_Endcap;
0303       if (EtaParticle < 0.0)
0304         Zend = -Zend;
0305       float Zlen = Zend - Zvertex;
0306       float RR = Zlen / sinh(EtaParticle);
0307       Theta = atan((RR + plane_Radius) / Zend);
0308       if (Theta < 0.0)
0309         Theta = Theta + Geom::pi();
0310       ETA = -log(tan(0.5 * Theta));
0311     }
0312 
0313     return ETA;
0314   } else {
0315     LogWarning("") << "[ContainmentCorrectionAnalyzer::ecalEta] Warning: Eta "
0316                       "equals to zero, not correcting";
0317     return EtaParticle;
0318   }
0319 }
0320 
0321 // taken from an old version of RecoEgamma/EgammaMCTools/src/PhotonMCTruthFinder
0322 std::vector<EcalSimPhotonMCTruth> ContainmentCorrectionAnalyzer::findMcTruth(std::vector<SimTrack> &theSimTracks,
0323                                                                              std::vector<SimVertex> &theSimVertices) {
0324   std::vector<EcalSimPhotonMCTruth> result;
0325 
0326   geantToIndex_.clear();
0327   // int   idTrk1_[10]; // UNUSED
0328   // int   idTrk2_[10]; // UNUSED
0329 
0330   // Local variables
0331   // const int SINGLE=1; // UNUSED
0332   // const int DOUBLE=2; // UNUSED
0333   // const int PYTHIA=3; // UNUSED
0334   const int ELECTRON_FLAV = 1;
0335   const int PIZERO_FLAV = 2;
0336   const int PHOTON_FLAV = 3;
0337 
0338   // int ievtype=0; // UNUSED
0339   int ievflav = 0;
0340   std::vector<SimTrack *> photonTracks;
0341   std::vector<SimTrack *> pizeroTracks;
0342   std::vector<const SimTrack *> trkFromConversion;
0343   SimVertex primVtx;
0344   std::vector<int> convInd;
0345 
0346   fillMcTruth(theSimTracks, theSimVertices);
0347   int iPV = -1;
0348   int partType1 = 0;
0349   int partType2 = 0;
0350   std::vector<SimTrack>::iterator iFirstSimTk = theSimTracks.begin();
0351   if (!(*iFirstSimTk).noVertex()) {
0352     iPV = (*iFirstSimTk).vertIndex();
0353     int vtxId = (*iFirstSimTk).vertIndex();
0354     primVtx = theSimVertices[vtxId];
0355     partType1 = (*iFirstSimTk).type();
0356   }
0357 
0358   // Look at a second track
0359   iFirstSimTk++;
0360   if (iFirstSimTk != theSimTracks.end()) {
0361     if ((*iFirstSimTk).vertIndex() == iPV) {
0362       partType2 = (*iFirstSimTk).type();
0363     }
0364   }
0365   int npv = 0;
0366   int iPho = 0;
0367   for (std::vector<SimTrack>::iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk) {
0368     if ((*iSimTk).noVertex())
0369       continue;
0370     // int vertexId = (*iSimTk).vertIndex(); // UNUSED
0371     // SimVertex vertex = theSimVertices[vertexId]; // UNUSED
0372     if ((*iSimTk).vertIndex() == iPV) {
0373       npv++;
0374       if ((*iSimTk).type() == 22) {
0375         convInd.push_back(0);
0376         photonTracks.push_back(&(*iSimTk));
0377         // math::XYZTLorentzVectorD momentum = (*iSimTk).momentum(); // UNUSED
0378       }
0379     }
0380   }
0381 
0382   if (npv > 4) {  // ievtype = PYTHIA; // UNUSED
0383   } else if (npv == 1) {
0384     if (abs(partType1) == 11) { /* ievtype = SINGLE; ==UNUSED== */
0385       ievflav = ELECTRON_FLAV;
0386     } else if (partType1 == 111) { /* ievtype = SINGLE; ==UNUSED== */
0387       ievflav = PIZERO_FLAV;
0388     } else if (partType1 == 22) { /* ievtype = SINGLE; ==UNUSED== */
0389       ievflav = PHOTON_FLAV;
0390     }
0391   } else if (npv == 2) {
0392     if (abs(partType1) == 11 && abs(partType2) == 11) { /* ievtype = DOUBLE; ==UNUSED== */
0393       ievflav = ELECTRON_FLAV;
0394     } else if (partType1 == 111 && partType2 == 111) { /* ievtype = DOUBLE; ==UNUSED== */
0395       ievflav = PIZERO_FLAV;
0396     } else if (partType1 == 22 && partType2 == 22) { /* ievtype = DOUBLE; ==UNUSED== */
0397       ievflav = PHOTON_FLAV;
0398     }
0399   }
0400 
0401   //  Look into converted photons
0402   int isAconversion = 0;
0403   if (ievflav == PHOTON_FLAV) {
0404     int nConv = 0;
0405     iPho = 0;
0406     for (std::vector<SimTrack *>::iterator iPhoTk = photonTracks.begin(); iPhoTk != photonTracks.end(); ++iPhoTk) {
0407       trkFromConversion.clear();
0408       for (std::vector<SimTrack>::iterator iSimTk = theSimTracks.begin(); iSimTk != theSimTracks.end(); ++iSimTk) {
0409         if ((*iSimTk).noVertex())
0410           continue;
0411         if ((*iSimTk).vertIndex() == iPV)
0412           continue;
0413         if (abs((*iSimTk).type()) != 11)
0414           continue;
0415         int vertexId = (*iSimTk).vertIndex();
0416         SimVertex vertex = theSimVertices[vertexId];
0417         int motherId = -1;
0418         if (vertex.parentIndex()) {
0419           unsigned motherGeantId = vertex.parentIndex();
0420           std::map<unsigned, unsigned>::iterator association = geantToIndex_.find(motherGeantId);
0421           if (association != geantToIndex_.end())
0422             motherId = association->second;
0423           // int motherType = motherId == -1 ? 0 :
0424           // theSimTracks[motherId].type();
0425 
0426           if (theSimTracks[motherId].trackId() == (*iPhoTk)->trackId()) {
0427             /// store this electron since it's from a converted photon
0428             trkFromConversion.push_back(&(*iSimTk));
0429           }
0430         }
0431       }  // loop over the SimTracks
0432 
0433       if (!trkFromConversion.empty()) {
0434         isAconversion = 1;
0435         nConv++;
0436         convInd[iPho] = nConv;
0437         int convVtxId = trkFromConversion[0]->vertIndex();
0438         SimVertex convVtx = theSimVertices[convVtxId];
0439         const math::XYZTLorentzVectorD &vtxPosition = convVtx.position();
0440         // math::XYZTLorentzVectorD momentum = (*iPhoTk)->momentum(); // UNUSED
0441 
0442         result.push_back(EcalSimPhotonMCTruth(isAconversion,
0443                                               (*iPhoTk)->momentum(),
0444                                               vtxPosition.pt(),
0445                                               vtxPosition.z(),
0446                                               vtxPosition,
0447                                               primVtx.position(),
0448                                               trkFromConversion));
0449       } else {
0450         isAconversion = 0;
0451         math::XYZTLorentzVectorD vtxPosition(0., 0., 0., 0.);
0452         result.push_back(EcalSimPhotonMCTruth(isAconversion,
0453                                               (*iPhoTk)->momentum(),
0454                                               vtxPosition.pt(),
0455                                               vtxPosition.z(),
0456                                               vtxPosition,
0457                                               primVtx.position(),
0458                                               trkFromConversion));
0459       }
0460       iPho++;
0461     }  // loop over the primary photons
0462   }    // Event with one or two photons
0463 
0464   return result;
0465 }
0466 
0467 void ContainmentCorrectionAnalyzer::fillMcTruth(std::vector<SimTrack> &simTracks, std::vector<SimVertex> &simVertices) {
0468   unsigned nVtx = simVertices.size();
0469   unsigned nTks = simTracks.size();
0470   if (nVtx == 0)
0471     return;
0472   // create a map associating geant particle id and position in the event
0473   // SimTrack vector
0474   for (unsigned it = 0; it < nTks; ++it) {
0475     geantToIndex_[simTracks[it].trackId()] = it;
0476   }
0477 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!