File indexing completed on 2024-09-10 02:59:06
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
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
0077 unsigned int total_num_apd_hits_seen[2];
0078 unsigned int num_apd_hits_no_parent[2];
0079
0080
0081
0082 unsigned int num_apd_hits_no_simtrack[2];
0083
0084
0085 unsigned int num_apd_hits_no_gen_particle[2];
0086
0087
0088
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
0175 const edm::Handle<edm::PCaloHitContainer>& caloHitEB = e.getHandle(tok_eb_);
0176
0177
0178 const edm::Handle<edm::PCaloHitContainer>& caloHitEE = e.getHandle(tok_ee_);
0179
0180
0181 const edm::Handle<edm::PCaloHitContainer>& caloHitES = e.getHandle(tok_es_);
0182
0183
0184 const edm::Handle<edm::PCaloHitContainer>& caloHitHC = e.getHandle(tok_hc_);
0185
0186
0187 const edm::Handle<edm::SimTrackContainer>& simTk = e.getHandle(tok_tk_);
0188
0189
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
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
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
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
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
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
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 }
0327 }
0328
0329 if (!found)
0330 ++noSimTrack[id];
0331 }
0332 if (histos)
0333 hitType[id]->Fill((double)(flag));
0334 }
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
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
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
0372 particle_type_count[apd_pid] = 1;
0373 else
0374 ++count_it->second;
0375
0376
0377
0378
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 }
0394 }
0395
0396 if (!found)
0397 ++num_apd_hits_no_simtrack[depth - 1];
0398 }
0399 }
0400 }
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);
0416
0417 edm::SimTrackContainer::const_iterator oldest_parent_track = parentSimTrack(simTk, simVtx, track);
0418
0419 int vertex_index = oldest_parent_track->vertIndex();
0420
0421
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
0478 if (thisTrkItr->trackId() == simTkId)
0479 return true;
0480
0481
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
0513 DEFINE_FWK_MODULE(HitParentTest);