Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:29:51 

0001 #include <algorithm>
0002 #include <fstream>
0003 #include <iomanip>
0004 #include <iostream>
0005 #include <map>
0006 #include <memory>
0007 #include <string>
0008 #include <vector>
0009 
0010 #include <TH1F.h>
0011 
0012 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0013 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/Framework/interface/Frameworkfwd.h"
0017 #include "FWCore/Framework/interface/MakerMacros.h"
0018 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/ServiceRegistry/interface/Service.h"
0022 #include "FWCore/Utilities/interface/Exception.h"
0023 #include "FWCore/Utilities/interface/InputTag.h"
0024 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0025 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0026 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0027 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0028 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0029 #include "SimG4CMS/Calo/interface/CaloHitID.h"
0030 
0031 class HitParentTest : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0032 public:
0033   HitParentTest(const edm::ParameterSet& ps);
0034   ~HitParentTest() override {}
0035   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0036 
0037 protected:
0038   void beginJob() override;
0039   void beginRun(edm::Run const&, edm::EventSetup const&) override {}
0040   void analyze(const edm::Event& e, const edm::EventSetup& c) override;
0041   void endRun(edm::Run const&, edm::EventSetup const&) override {}
0042   void endJob() override;
0043 
0044 private:
0045   /** performs some checks on hits */
0046   void analyzeHits(const std::vector<PCaloHit>&,
0047                    const edm::Handle<edm::SimTrackContainer>&,
0048                    const edm::Handle<edm::SimVertexContainer>&,
0049                    int type);
0050   void analyzeAPDHits(const std::vector<PCaloHit>&,
0051                       const edm::Handle<edm::SimTrackContainer>&,
0052                       const edm::Handle<edm::SimVertexContainer>&,
0053                       int depth);
0054 
0055   bool simTrackPresent(const edm::Handle<edm::SimTrackContainer>&, int id);
0056   math::XYZTLorentzVectorD getOldestParentVertex(const edm::Handle<edm::SimTrackContainer>&,
0057                                                  const edm::Handle<edm::SimVertexContainer>&,
0058                                                  edm::SimTrackContainer::const_iterator track);
0059   edm::SimTrackContainer::const_iterator parentSimTrack(const edm::Handle<edm::SimTrackContainer>&,
0060                                                         const edm::Handle<edm::SimVertexContainer>&,
0061                                                         edm::SimTrackContainer::const_iterator thisTrkItr);
0062   bool validSimTrack(const edm::Handle<edm::SimTrackContainer>&,
0063                      const edm::Handle<edm::SimVertexContainer>&,
0064                      unsigned int simTkId,
0065                      edm::SimTrackContainer::const_iterator thisTrkItr);
0066 
0067 private:
0068   const std::string sourceLabel, g4Label_, hitLabEB_, hitLabEE_, hitLabES_, hitLabHC_;
0069   const edm::EDGetTokenT<edm::PCaloHitContainer> tok_eb_;
0070   const edm::EDGetTokenT<edm::PCaloHitContainer> tok_ee_;
0071   const edm::EDGetTokenT<edm::PCaloHitContainer> tok_es_;
0072   const edm::EDGetTokenT<edm::PCaloHitContainer> tok_hc_;
0073   const edm::EDGetTokenT<edm::SimTrackContainer> tok_tk_;
0074   const edm::EDGetTokenT<edm::SimVertexContainer> tok_vtx_;
0075 
0076   /** error and other counters */
0077   unsigned int total_num_apd_hits_seen[2];
0078   unsigned int num_apd_hits_no_parent[2];
0079 
0080   /** number of apd hits for which the parent sim track id was not found in
0081       the simtrack collection. */
0082   unsigned int num_apd_hits_no_simtrack[2];
0083 
0084   /** number of APD hits for which no generator particle was found */
0085   unsigned int num_apd_hits_no_gen_particle[2];
0086 
0087   /** 'histogram' of types of particles going through the APD. Maps from numeric particle code
0088       to the number of occurences. */
0089   std::map<int, unsigned> particle_type_count;
0090 
0091   std::string detector[7];
0092   bool histos;
0093   unsigned int totalHits[7], noParent[7];
0094   unsigned int noSimTrack[7], noGenParticle[7];
0095   TH1F *hitType[7], *hitRho[7], *hitZ[7];
0096 };
0097 
0098 HitParentTest::HitParentTest(const edm::ParameterSet& ps)
0099     : sourceLabel(ps.getUntrackedParameter<std::string>("SourceLabel", "VtxSmeared")),
0100       g4Label_(ps.getUntrackedParameter<std::string>("ModuleLabel", "g4SimHits")),
0101       hitLabEB_(ps.getUntrackedParameter<std::string>("EBCollection", "EcalHitsEB")),
0102       hitLabEE_(ps.getUntrackedParameter<std::string>("EECollection", "EcalHitsEE")),
0103       hitLabES_(ps.getUntrackedParameter<std::string>("ESCollection", "EcalHitsES")),
0104       hitLabHC_(ps.getUntrackedParameter<std::string>("HCCollection", "HcalHits")),
0105       tok_eb_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, hitLabEB_))),
0106       tok_ee_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, hitLabEE_))),
0107       tok_es_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, hitLabES_))),
0108       tok_hc_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label_, hitLabHC_))),
0109       tok_tk_(consumes<edm::SimTrackContainer>(edm::InputTag(g4Label_))),
0110       tok_vtx_(consumes<edm::SimVertexContainer>(edm::InputTag(g4Label_))) {
0111   usesResource(TFileService::kSharedResource);
0112 
0113   edm::LogVerbatim("HitParentTest") << "Module Label: " << g4Label_ << "   Hits: " << hitLabEB_ << ", " << hitLabEE_
0114                                     << ", " << hitLabES_ << ", " << hitLabHC_;
0115 
0116   for (unsigned int i = 0; i < 2; ++i) {
0117     total_num_apd_hits_seen[i] = 0;
0118     num_apd_hits_no_parent[i] = 0;
0119     num_apd_hits_no_simtrack[i] = 0;
0120     num_apd_hits_no_gen_particle[i] = 0;
0121   }
0122   std::string dets[7] = {"EB", "EE", "ES", "HB", "HE", "HO", "HF"};
0123   for (unsigned int i = 0; i < 7; ++i) {
0124     detector[i] = dets[i];
0125     totalHits[i] = 0;
0126     noParent[i] = 0;
0127     noSimTrack[i] = 0;
0128     noGenParticle[i] = 0;
0129   }
0130 }
0131 
0132 void HitParentTest::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0133   edm::ParameterSetDescription desc;
0134   desc.addUntracked<std::string>("SourceLabel", "generatorSmeared");
0135   desc.addUntracked<std::string>("ModuleLabel", "g4SimHits");
0136   desc.addUntracked<std::string>("EBCollection", "EcalHitsEB");
0137   desc.addUntracked<std::string>("EECollection", "EcalHitsEE");
0138   desc.addUntracked<std::string>("ESCollection", "EcalHitsES");
0139   desc.addUntracked<std::string>("HCCollection", "HcalHits");
0140   descriptions.add("HitParentTest", desc);
0141 }
0142 
0143 void HitParentTest::beginJob() {
0144   edm::Service<TFileService> tfile;
0145   if (!tfile.isAvailable()) {
0146     edm::LogVerbatim("HitParentTest") << "TFileService unavailable: no histograms";
0147     histos = false;
0148   } else {
0149     char name[20], title[100];
0150     histos = true;
0151     for (unsigned int i = 0; i < 7; ++i) {
0152       sprintf(name, "HitType%d", i);
0153       sprintf(title, "Hit types for %s", detector[i].c_str());
0154       hitType[i] = tfile->make<TH1F>(name, title, 10, 0., 10.);
0155       hitType[i]->GetXaxis()->SetTitle(title);
0156       hitType[i]->GetYaxis()->SetTitle("Hits");
0157       sprintf(name, "RhoVertex%d", i);
0158       sprintf(title, "#rho of the vertex for %s Hits", detector[i].c_str());
0159       hitRho[i] = tfile->make<TH1F>(name, title, 10000, 0., 100.);
0160       hitRho[i]->GetXaxis()->SetTitle(title);
0161       hitRho[i]->GetYaxis()->SetTitle("Hits");
0162       sprintf(name, "ZVertex%d", i);
0163       sprintf(title, "z of the vertex for %s Hits", detector[i].c_str());
0164       hitZ[i] = tfile->make<TH1F>(name, title, 2000, -100., 100.);
0165       hitZ[i]->GetXaxis()->SetTitle(title);
0166       hitZ[i]->GetYaxis()->SetTitle("Hits");
0167     }
0168   }
0169 }
0170 
0171 void HitParentTest::analyze(const edm::Event& e, const edm::EventSetup&) {
0172   edm::LogVerbatim("HitParentTest") << "HitParentTes::Run = " << e.id().run() << " Event = " << e.id().event();
0173 
0174   // get PCaloHits for ecal barrel
0175   const edm::Handle<edm::PCaloHitContainer>& caloHitEB = e.getHandle(tok_eb_);
0176 
0177   // get PCaloHits for ecal endcap
0178   const edm::Handle<edm::PCaloHitContainer>& caloHitEE = e.getHandle(tok_ee_);
0179 
0180   // get PCaloHits for preshower
0181   const edm::Handle<edm::PCaloHitContainer>& caloHitES = e.getHandle(tok_es_);
0182 
0183   // get PCaloHits for hcal
0184   const edm::Handle<edm::PCaloHitContainer>& caloHitHC = e.getHandle(tok_hc_);
0185 
0186   // get sim tracks
0187   const edm::Handle<edm::SimTrackContainer>& simTk = e.getHandle(tok_tk_);
0188 
0189   // get sim vertices
0190   const edm::Handle<edm::SimVertexContainer>& simVtx = e.getHandle(tok_vtx_);
0191 
0192   edm::LogVerbatim("HitParentTest") << "HitParentTest: hits valid[EB]: " << caloHitEB.isValid()
0193                                     << " valid[EE]: " << caloHitEE.isValid() << " valid[ES]: " << caloHitES.isValid()
0194                                     << " valid[HC]: " << caloHitHC.isValid();
0195 
0196   if (caloHitEB.isValid()) {
0197     for (int depth = 1; depth <= 2; ++depth)
0198       analyzeAPDHits(*caloHitEB, simTk, simVtx, depth);
0199     analyzeHits(*caloHitEB, simTk, simVtx, 0);
0200   }
0201   if (caloHitEE.isValid())
0202     analyzeHits(*caloHitEE, simTk, simVtx, 1);
0203   if (caloHitES.isValid())
0204     analyzeHits(*caloHitES, simTk, simVtx, 2);
0205   if (caloHitHC.isValid())
0206     analyzeHits(*caloHitHC, simTk, simVtx, 3);
0207 }
0208 
0209 /** define the comparison for sorting the particle ids counting map */
0210 class HitParentTestComparison {
0211 protected:
0212   const std::map<int, unsigned>& particle_type_count;
0213 
0214 public:
0215   HitParentTestComparison(const std::map<int, unsigned>& _particle_type_count)
0216       : particle_type_count(_particle_type_count) {}
0217 
0218   bool operator()(int pid1, int pid2) const { return particle_type_count.at(pid1) > particle_type_count.at(pid2); }
0219 };
0220 
0221 void HitParentTest::endJob() {
0222   edm::LogVerbatim("HitParentTest") << "HitParentTest::Total number of APD hits seen: "
0223                                     << total_num_apd_hits_seen[0] + total_num_apd_hits_seen[1];
0224   for (unsigned int depth = 0; depth < 2; ++depth) {
0225     edm::LogVerbatim("HitParentTest") << "APD Hits in depth = " << depth + 1
0226                                       << " Total = " << total_num_apd_hits_seen[depth];
0227     edm::LogVerbatim("HitParentTest") << "summary of errors:\n"
0228                                       << " number of APD hits with zero geant track id: "
0229                                       << num_apd_hits_no_parent[depth] << "\n"
0230                                       << "number of APD hits for which the parent simtrack was not found "
0231                                       << "in the simtrack collection: " << num_apd_hits_no_simtrack[depth] << "\n"
0232                                       << "number of APD hits for which no generator particle was "
0233                                       << "found: " << num_apd_hits_no_gen_particle[depth] << "\n";
0234   }
0235 
0236   for (unsigned int det = 0; det < 7; ++det) {
0237     edm::LogVerbatim("HitParentTest") << "Total number of hits seen in " << detector[det] << ": " << totalHits[det];
0238     edm::LogVerbatim("HitParentTest") << "summary of errors:\n"
0239                                       << "number of hits with zero geant track id: " << noParent[det] << "\n"
0240                                       << "number of hits for which the parent simtrack was not found in "
0241                                       << "the simtrack collection: " << noSimTrack[det] << "\n"
0242                                       << "number of hits for which no generator particle was found: "
0243                                       << noGenParticle[det] << "\n";
0244   }
0245 
0246   // sort in decreasing order of occurence
0247   std::vector<int> sorted_pids;
0248   for (std::map<int, unsigned>::const_iterator it = particle_type_count.begin(); it != particle_type_count.end(); ++it)
0249     sorted_pids.push_back(it->first);
0250 
0251   // now sort the pids
0252 
0253   std::sort(sorted_pids.begin(), sorted_pids.end(), HitParentTestComparison(particle_type_count));
0254 
0255   edm::LogVerbatim("HitParentTest") << "HitParentTest::Frequency particle types through the APD "
0256                                     << "(pid/frequency):";
0257   for (unsigned i = 0; i < sorted_pids.size(); ++i) {
0258     int pid = sorted_pids[i];
0259     edm::LogVerbatim("HitParentTest") << "  pid " << std::setw(6) << pid << ": count " << std::setw(6)
0260                                       << particle_type_count[pid];
0261   }
0262 }
0263 
0264 void HitParentTest::analyzeHits(const std::vector<PCaloHit>& hits,
0265                                 const edm::Handle<edm::SimTrackContainer>& simTk,
0266                                 const edm::Handle<edm::SimVertexContainer>& simVtx,
0267                                 int type) {
0268   for (std::vector<PCaloHit>::const_iterator hit_it = hits.begin(); hit_it != hits.end(); ++hit_it) {
0269     int id(type), flag(0);
0270     if (type == 3) {
0271       HcalDetId id_ = HcalDetId(hit_it->id());
0272       int subdet = id_.subdet();
0273       if (subdet == static_cast<int>(HcalEndcap))
0274         id = type + 1;
0275       else if (subdet == static_cast<int>(HcalOuter))
0276         id = type + 2;
0277       else if (subdet == static_cast<int>(HcalForward))
0278         id = type + 3;
0279     }
0280     ++totalHits[id];
0281 
0282     // get the geant track id
0283     int hit_geant_track_id = hit_it->geantTrackId();
0284 
0285     if (hit_geant_track_id <= 0) {
0286       ++noParent[id];
0287       flag = 1;
0288     } else {
0289       bool found = false;
0290       flag = 2;
0291       // check whether this id is actually there in the list of simtracks
0292       for (edm::SimTrackContainer::const_iterator simTrkItr = simTk->begin(); simTrkItr != simTk->end() && !found;
0293            ++simTrkItr) {
0294         if (hit_geant_track_id == (int)(simTrkItr->trackId())) {
0295           found = true;
0296           flag = 3;
0297           bool match = validSimTrack(simTk, simVtx, hit_geant_track_id, simTrkItr);
0298 
0299           edm::LogVerbatim("HitParentTest")
0300               << "[" << detector[type] << "] Match = " << match << " hit_geant_track_id=" << hit_geant_track_id
0301               << " particle id=" << simTrkItr->type();
0302 
0303           if (!match) {
0304             edm::LogVerbatim("HitParentTest") << "NO MATCH FOUND !";
0305             ++noGenParticle[id];
0306           }
0307 
0308           // beam pipe...
0309           int pid = simTrkItr->type();
0310           math::XYZTLorentzVectorD oldest_parent_vertex = getOldestParentVertex(simTk, simVtx, simTrkItr);
0311 
0312           edm::SimTrackContainer::const_iterator oldest_parent_track = parentSimTrack(simTk, simVtx, simTrkItr);
0313 
0314           edm::LogVerbatim("HitParentTest")
0315               << "[" << detector[type] << "] Hit pid = " << pid << "  hit track id = " << hit_geant_track_id
0316               << " Oldest Parent's Vertex: " << oldest_parent_vertex << " rho = " << oldest_parent_vertex.Rho()
0317               << " Oldest Parent's pid: " << oldest_parent_track->type()
0318               << " Oldest Parent's track id: " << oldest_parent_track->trackId()
0319               << "\nHit vertex index: " << simTrkItr->vertIndex() << " (tot #vertices: " << simVtx->size() << ")"
0320               << "\nHit vertex parent track: " << (*simVtx)[simTrkItr->vertIndex()].parentIndex()
0321               << " present=" << simTrackPresent(simTk, (*simVtx)[simTrkItr->vertIndex()].parentIndex());
0322           if (histos) {
0323             hitRho[id]->Fill(oldest_parent_vertex.Rho());
0324             hitZ[id]->Fill(oldest_parent_vertex.Z());
0325           }
0326         }  // a match was found
0327       }    // geant track id found in simtracks
0328 
0329       if (!found)
0330         ++noSimTrack[id];
0331     }  // hits with a valid SimTrack Id
0332     if (histos)
0333       hitType[id]->Fill((double)(flag));
0334   }  // loop over all calohits (of the given depth)
0335 }
0336 
0337 void HitParentTest::analyzeAPDHits(const std::vector<PCaloHit>& hits,
0338                                    const edm::Handle<edm::SimTrackContainer>& simTk,
0339                                    const edm::Handle<edm::SimVertexContainer>& simVtx,
0340                                    int depth) {
0341   for (std::vector<PCaloHit>::const_iterator hit_it = hits.begin(); hit_it != hits.end(); ++hit_it) {
0342     if (hit_it->depth() == depth) {
0343       ++total_num_apd_hits_seen[depth - 1];
0344 
0345       // get the geant track id
0346       int hit_geant_track_id = hit_it->geantTrackId();
0347 
0348       if (hit_geant_track_id <= 0) {
0349         ++num_apd_hits_no_parent[depth - 1];
0350       } else {
0351         bool found = false;
0352         // check whether this id is actually there in the list of simtracks
0353         for (edm::SimTrackContainer::const_iterator simTrkItr = simTk->begin(); simTrkItr != simTk->end() && !found;
0354              ++simTrkItr) {
0355           if (hit_geant_track_id == (int)(simTrkItr->trackId())) {
0356             found = true;
0357             bool match = validSimTrack(simTk, simVtx, hit_geant_track_id, simTrkItr);
0358 
0359             edm::LogVerbatim("HitParentTest")
0360                 << "APDHIT Match = " << match << " hit_geant_track_id = " << hit_geant_track_id
0361                 << " particle id=" << simTrkItr->type();
0362 
0363             if (!match) {
0364               edm::LogVerbatim("HitParentTest") << "NO MATCH FOUND !";
0365               ++num_apd_hits_no_gen_particle[depth - 1];
0366             }
0367 
0368             int apd_pid = simTrkItr->type();
0369             std::map<int, unsigned>::iterator count_it = particle_type_count.find(apd_pid);
0370             if (count_it == particle_type_count.end())
0371               // first occurence of this particle pid
0372               particle_type_count[apd_pid] = 1;
0373             else
0374               ++count_it->second;
0375 
0376             //--------------------
0377             // check where the oldest parent has its vertex. Should be close to the
0378             // beam pipe...
0379             math::XYZTLorentzVectorD oldest_parent_vertex = getOldestParentVertex(simTk, simVtx, simTrkItr);
0380 
0381             edm::SimTrackContainer::const_iterator oldest_parent_track = parentSimTrack(simTk, simVtx, simTrkItr);
0382 
0383             edm::LogVerbatim("HitParentTest")
0384                 << "APD hit pid = " << apd_pid << " APD hit track id = " << hit_geant_track_id
0385                 << " depth = " << hit_it->depth() << " OLDEST PARENT'S VERTEX: " << oldest_parent_vertex
0386                 << " rho = " << oldest_parent_vertex.Rho() << " OLDEST PARENT'S PID: " << oldest_parent_track->type()
0387                 << " OLDEST PARENT'S track id: " << oldest_parent_track->trackId() << "\n"
0388                 << "APD hit vertex index: " << simTrkItr->vertIndex() << " (tot #vertices: " << simVtx->size() << ")"
0389                 << "\n"
0390                 << "APD hit vertex parent track: " << (*simVtx)[simTrkItr->vertIndex()].parentIndex()
0391                 << " present=" << simTrackPresent(simTk, (*simVtx)[simTrkItr->vertIndex()].parentIndex());
0392 
0393           }  // a match was found
0394         }    // geant track id found in simtracks
0395 
0396         if (!found)
0397           ++num_apd_hits_no_simtrack[depth - 1];
0398       }  // hits with a valid SimTrack Id
0399     }    // right depth index
0400   }      // loop over all calohits (of the given depth)
0401 }
0402 
0403 bool HitParentTest::simTrackPresent(const edm::Handle<edm::SimTrackContainer>& simTk, int id) {
0404   for (edm::SimTrackContainer::const_iterator simTrkItr = simTk->begin(); simTrkItr != simTk->end(); ++simTrkItr) {
0405     if ((int)(simTrkItr->trackId()) == id)
0406       return true;
0407   }
0408   return false;
0409 }
0410 
0411 math::XYZTLorentzVectorD HitParentTest::getOldestParentVertex(const edm::Handle<edm::SimTrackContainer>& simTk,
0412                                                               const edm::Handle<edm::SimVertexContainer>& simVtx,
0413                                                               edm::SimTrackContainer::const_iterator track) {
0414   static const math::XYZTLorentzVectorD invalid_vertex(
0415       10000, 10000, 10000, 10000);  // default value if no valid vertex found
0416 
0417   edm::SimTrackContainer::const_iterator oldest_parent_track = parentSimTrack(simTk, simVtx, track);
0418 
0419   int vertex_index = oldest_parent_track->vertIndex();
0420 
0421   // sanity checks
0422   if (vertex_index < 0 || vertex_index >= (int)(simVtx->size()))
0423     return invalid_vertex;
0424 
0425   return (*simVtx)[vertex_index].position();
0426 }
0427 
0428 edm::SimTrackContainer::const_iterator HitParentTest::parentSimTrack(const edm::Handle<edm::SimTrackContainer>& simTk,
0429                                                                      const edm::Handle<edm::SimVertexContainer>& simVtx,
0430                                                                      edm::SimTrackContainer::const_iterator thisTrkItr) {
0431   edm::SimTrackContainer::const_iterator itr = simTk->end();
0432 
0433   int vertIndex = thisTrkItr->vertIndex();
0434   int type = thisTrkItr->type();
0435   int charge = (int)thisTrkItr->charge();
0436   edm::LogVerbatim("HitParentTest") << "SimTrackParent:: " << thisTrkItr->trackId() << " Vertex " << vertIndex
0437                                     << " Type " << type << " Charge " << charge;
0438 
0439   if (vertIndex == -1)
0440     return thisTrkItr;
0441   else if (vertIndex >= (int)simVtx->size())
0442     return itr;
0443 
0444   edm::SimVertexContainer::const_iterator simVtxItr = simVtx->begin();
0445   for (int iv = 0; iv < vertIndex; iv++)
0446     simVtxItr++;
0447   int parent = simVtxItr->parentIndex();
0448 
0449   if (parent < 0 && simVtxItr != simVtx->begin()) {
0450     const math::XYZTLorentzVectorD pos1 = simVtxItr->position();
0451     for (simVtxItr = simVtx->begin(); simVtxItr != simVtx->end(); ++simVtxItr) {
0452       if (simVtxItr->parentIndex() > 0) {
0453         const math::XYZTLorentzVectorD pos2 = pos1 - simVtxItr->position();
0454         double dist = pos2.P();
0455         if (dist < 0.001) {
0456           parent = simVtxItr->parentIndex();
0457           break;
0458         }
0459       }
0460     }
0461   }
0462   for (edm::SimTrackContainer::const_iterator simTrkItr = simTk->begin(); simTrkItr != simTk->end(); simTrkItr++) {
0463     if ((int)simTrkItr->trackId() == parent && simTrkItr != thisTrkItr)
0464       return parentSimTrack(simTk, simVtx, simTrkItr);
0465   }
0466 
0467   return thisTrkItr;
0468 }
0469 
0470 bool HitParentTest::validSimTrack(const edm::Handle<edm::SimTrackContainer>& simTk,
0471                                   const edm::Handle<edm::SimVertexContainer>& simVtx,
0472                                   unsigned int simTkId,
0473                                   edm::SimTrackContainer::const_iterator thisTrkItr) {
0474   edm::LogVerbatim("HitParentTest") << "Inside validSimTrack: trackId " << thisTrkItr->trackId() << " vtxIndex "
0475                                     << thisTrkItr->vertIndex() << " to be matched to " << simTkId;
0476 
0477   //This track originates from simTkId
0478   if (thisTrkItr->trackId() == simTkId)
0479     return true;
0480 
0481   //Otherwise trace back the history using SimTracks and SimVertices
0482   int vertIndex = thisTrkItr->vertIndex();
0483   if (vertIndex == -1 || vertIndex >= (int)simVtx->size())
0484     return false;
0485   edm::SimVertexContainer::const_iterator simVtxItr = simVtx->begin();
0486   for (int iv = 0; iv < vertIndex; iv++)
0487     simVtxItr++;
0488   int parent = simVtxItr->parentIndex();
0489   edm::LogVerbatim("HitParentTest") << "validSimTrack:: parent index " << parent << " ";
0490   if (parent < 0 && simVtxItr != simVtx->begin()) {
0491     const math::XYZTLorentzVectorD pos1 = simVtxItr->position();
0492     for (simVtxItr = simVtx->begin(); simVtxItr != simVtx->end(); ++simVtxItr) {
0493       if (simVtxItr->parentIndex() > 0) {
0494         const math::XYZTLorentzVectorD pos2 = pos1 - simVtxItr->position();
0495         double dist = pos2.P();
0496         if (dist < 0.001) {
0497           parent = simVtxItr->parentIndex();
0498           break;
0499         }
0500       }
0501     }
0502   }
0503   edm::LogVerbatim("HitParentTest") << "HitParentTes::final index " << parent;
0504   for (edm::SimTrackContainer::const_iterator simTrkItr = simTk->begin(); simTrkItr != simTk->end(); simTrkItr++) {
0505     if ((int)simTrkItr->trackId() == parent && simTrkItr != thisTrkItr)
0506       return validSimTrack(simTk, simVtx, simTkId, simTrkItr);
0507   }
0508 
0509   return false;
0510 }
0511 
0512 //define this as a plug-in
0513 DEFINE_FWK_MODULE(HitParentTest);
[ 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!