Event/src/FBaseSimEvent.cc

0001 //HepMC Headers
0002 #include "HepMC/GenEvent.h"
0003 #include "HepMC/GenVertex.h"
0004 #include "HepMC/GenParticle.h"
0005
0006 //Framework Headers
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008
0009 //CMSSW Data Formats
0010 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0011 #include "DataFormats/Candidate/interface/Candidate.h"
0012
0013 //FAMOS Headers
0014 #include "FastSimulation/Particle/interface/pdg_functions.h"
0015 #include "FastSimulation/Particle/interface/makeParticle.h"
0016 #include "FastSimulation/Event/interface/FBaseSimEvent.h"
0017 #include "FastSimulation/Event/interface/FSimTrack.h"
0018 #include "FastSimulation/Event/interface/FSimVertex.h"
0019 #include "FastSimulation/Event/interface/KineParticleFilter.h"
0020 #include "CommonTools/BaseParticlePropagator/interface/BaseParticlePropagator.h"
0021
0022 #include "FastSimDataFormats/NuclearInteractions/interface/FSimVertexType.h"
0023
0024 using namespace HepPDT;
0025
0026 // system include
0027 #include <iostream>
0028 #include <iomanip>
0029 #include <map>
0030 #include <string>
0031
0032 FBaseSimEvent::FBaseSimEvent(const edm::ParameterSet& kine)
0033     : nSimTracks(0), nSimVertices(0), nGenParticles(0), nChargedParticleTracks(0), initialSize(5000) {
0034   // Initialize the vectors of particles and vertices
0035   theGenParticles = new std::vector<HepMC::GenParticle*>();
0036   theSimTracks = new std::vector<FSimTrack>;
0037   theSimVertices = new std::vector<FSimVertex>;
0038   theChargedTracks = new std::vector<unsigned>();
0039   theFSimVerticesType = new FSimVertexTypeCollection();
0040
0041   // Reserve some size to avoid mutiple copies
0042   /* */
0043   theSimTracks->resize(initialSize);
0044   theSimVertices->resize(initialSize);
0045   theGenParticles->resize(initialSize);
0046   theChargedTracks->resize(initialSize);
0047   theFSimVerticesType->resize(initialSize);
0048   theTrackSize = initialSize;
0049   theVertexSize = initialSize;
0050   theGenSize = initialSize;
0051   theChargedSize = initialSize;
0052   /* */
0053
0054   // Initialize the Particle filter
0055   myFilter = new KineParticleFilter(kine);
0056
0057   // Initialize the distance from (0,0,0) after which *generated* particles are
0058   // no longer considered - because the mother could have interacted before.
0059   // unit : cm x cm
0060   lateVertexPosition = 2.5 * 2.5;
0061 }
0062
0063 FBaseSimEvent::~FBaseSimEvent() {
0064   // Clear the vectors
0065   theGenParticles->clear();
0066   theSimTracks->clear();
0067   theSimVertices->clear();
0068   theChargedTracks->clear();
0069   theFSimVerticesType->clear();
0070
0071   // Delete
0072   delete theGenParticles;
0073   delete theSimTracks;
0074   delete theSimVertices;
0075   delete theChargedTracks;
0076   delete theFSimVerticesType;
0077   delete myFilter;
0078 }
0079
0080 void FBaseSimEvent::initializePdt(const HepPDT::ParticleDataTable* aPdt) { pdt = aPdt; }
0081
0082 void FBaseSimEvent::fill(const HepMC::GenEvent& myGenEvent) {
0083   // Clear old vectors
0084   clear();
0085
0086   // Add the particles in the FSimEvent
0087   addParticles(myGenEvent);
0088 }
0089
0090 void FBaseSimEvent::fill(const std::vector<SimTrack>& simTracks, const std::vector<SimVertex>& simVertices) {
0091   // Watch out there ! A SimVertex is in mm (stupid),
0092   //            while a FSimVertex is in cm (clever).
0093
0094   clear();
0095
0096   unsigned nVtx = simVertices.size();
0097   unsigned nTks = simTracks.size();
0098
0099   // Empty event, do nothin'
0100   if (nVtx == 0)
0101     return;
0102
0103   // Two arrays for internal use.
0104   std::vector<int> myVertices(nVtx, -1);
0105   std::vector<int> myTracks(nTks, -1);
0106
0107   // create a map associating geant particle id and position in the
0108   // event SimTrack vector
0109
0110   std::map<unsigned, unsigned> geantToIndex;
0111   for (unsigned it = 0; it < simTracks.size(); ++it) {
0112     geantToIndex[simTracks[it].trackId()] = it;
0113   }
0114
0115   // Create also a map associating a SimTrack with its endVertex
0116   /*
0117   std::map<unsigned, unsigned> endVertex;
0118   for ( unsigned iv=0; iv<simVertices.size(); ++iv ) {
0119     endVertex[ simVertices[iv].parentIndex() ] = iv;
0120   }
0121   */
0122
0123   // Set the main vertex for the kine particle filter
0124   // SimVertices were in mm until 110_pre2
0125   //  HepLorentzVector primaryVertex = simVertices[0].position()/10.;
0126   // SImVertices are now in cm
0127   // Also : position is copied until SimVertex switches to Mathcore.
0128   //  XYZTLorentzVector primaryVertex = simVertices[0].position();
0129   // The next 5 lines to be then replaced by the previous line
0130   XYZTLorentzVector primaryVertex(simVertices[0].position().x(),
0131                                   simVertices[0].position().y(),
0132                                   simVertices[0].position().z(),
0133                                   simVertices[0].position().t());
0134   //
0135   //myFilter->setMainVertex(primaryVertex);
0136   // Add the main vertex to the list.
0137   addSimVertex(/*myFilter->vertex()*/ primaryVertex, -1, FSimVertexType::PRIMARY_VERTEX);
0138   myVertices[0] = 0;
0139
0140   for (unsigned trackId = 0; trackId < nTks; ++trackId) {
0141     // The track
0142     const SimTrack& track = simTracks[trackId];
0143     //    std::cout << std::endl << "SimTrack " << trackId << " " << track << std::endl;
0144
0145     // The origin vertex
0146     int vertexId = track.vertIndex();
0147     const SimVertex& vertex = simVertices[vertexId];
0148     //std::cout << "Origin vertex " << vertexId << " " << vertex << std::endl;
0149
0150     // The mother track
0151     int motherId = -1;
0152     if (!vertex.noParent()) {  // there is a parent to this vertex
0153       // geant id of the mother
0154       unsigned motherGeantId = vertex.parentIndex();
0155       std::map<unsigned, unsigned>::iterator association = geantToIndex.find(motherGeantId);
0156       if (association != geantToIndex.end())
0157         motherId = association->second;
0158     }
0159     int originId = motherId == -1 ? -1 : myTracks[motherId];
0160     //std::cout << "Origin id " << originId << std::endl;
0161
0162     /*
0163     if ( endVertex.find(trackId) != endVertex.end() )
0164       std::cout << "End vertex id = " << endVertex[trackId] << std::endl;
0165     else
0166       std::cout << "No endVertex !!! " << std::endl;
0167     std::cout << "Tracker surface position " << track.trackerSurfacePosition() << std::endl;
0168     */
0169
0170     // Add the vertex (if it does not already exist!)
0171     XYZTLorentzVector position(
0172         vertex.position().px(), vertex.position().py(), vertex.position().pz(), vertex.position().e());
0173     if (myVertices[vertexId] == -1)
0174       // Momentum and position are copied until SimTrack and SimVertex
0175       // switch to Mathcore.
0176       //      myVertices[vertexId] = addSimVertex(vertex.position(),originId);
0177       // The next line to be then replaced by the previous line
0178       myVertices[vertexId] = addSimVertex(position, originId);
0179
0180     // Add the track (with protection for brem'ing electrons and muons)
0181     int motherType = motherId == -1 ? 0 : simTracks[motherId].type();
0182
0183     bool notBremInDetector = (abs(motherType) != 11 && std::abs(motherType) != 13) || motherType != track.type() ||
0184                              position.Perp2() < lateVertexPosition;
0185
0186     if (notBremInDetector) {
0187       // Momentum and position are copied until SimTrack and SimVertex
0188       // switch to Mathcore.
0189       //      RawParticle part(track.momentum(), vertex.position());
0190       // The next 3 lines to be then replaced by the previous line
0191       XYZTLorentzVector momentum(
0192           track.momentum().px(), track.momentum().py(), track.momentum().pz(), track.momentum().e());
0193       RawParticle part = makeParticle(theTable(), track.type(), momentum, position);
0194       //
0195       //std::cout << "Ctau  = " << part.PDGcTau() << std::endl;
0196       // Don't save tracks that have decayed immediately but for which no daughters
0197       // were saved (probably due to cuts on E, pT and eta)
0198       //  if ( part.PDGcTau() > 0.1 || endVertex.find(trackId) != endVertex.end() )
0199       myTracks[trackId] = addSimTrack(&part, myVertices[vertexId], track.genpartIndex());
0200       if (myTracks[trackId] >= 0) {
0201         (*theSimTracks)[myTracks[trackId]].setTkPosition(track.trackerSurfacePosition());
0202         (*theSimTracks)[myTracks[trackId]].setTkMomentum(track.trackerSurfaceMomentum());
0203       }
0204     } else {
0205       myTracks[trackId] = myTracks[motherId];
0206       if (myTracks[trackId] >= 0) {
0207         (*theSimTracks)[myTracks[trackId]].setTkPosition(track.trackerSurfacePosition());
0208         (*theSimTracks)[myTracks[trackId]].setTkMomentum(track.trackerSurfaceMomentum());
0209       }
0210     }
0211   }
0212
0213   // Now loop over the remaining end vertices !
0214   for (unsigned vertexId = 0; vertexId < nVtx; ++vertexId) {
0215     // if the vertex is already saved, just ignore.
0216     if (myVertices[vertexId] != -1)
0217       continue;
0218
0219     // The yet unused vertex
0220     const SimVertex& vertex = simVertices[vertexId];
0221
0222     // The mother track
0223     int motherId = -1;
0224     if (!vertex.noParent()) {  // there is a parent to this vertex
0225
0226       // geant id of the mother
0227       unsigned motherGeantId = vertex.parentIndex();
0228       std::map<unsigned, unsigned>::iterator association = geantToIndex.find(motherGeantId);
0229       if (association != geantToIndex.end())
0230         motherId = association->second;
0231     }
0232     int originId = motherId == -1 ? -1 : myTracks[motherId];
0233
0234     // Add the vertex
0235     // Momentum and position are copied until SimTrack and SimVertex
0236     // switch to Mathcore.
0237     //    myVertices[vertexId] = addSimVertex(vertex.position(),originId);
0238     // The next 3 lines to be then replaced by the previous line
0239     XYZTLorentzVector position(
0240         vertex.position().px(), vertex.position().py(), vertex.position().pz(), vertex.position().e());
0241     myVertices[vertexId] = addSimVertex(position, originId);
0242   }
0243
0244   // Finally, propagate all particles to the calorimeters
0245   BaseParticlePropagator myPart;
0246   XYZTLorentzVector mom;
0247   XYZTLorentzVector pos;
0248
0249   // Loop over the tracks
0250   for (int fsimi = 0; fsimi < (int)nTracks(); ++fsimi) {
0251     FSimTrack& myTrack = track(fsimi);
0252     double trackerSurfaceTime =
0253         myTrack.vertex().position().t() + myTrack.momentum().e() / myTrack.momentum().pz() *
0254                                               (myTrack.trackerSurfacePosition().z() - myTrack.vertex().position().z());
0255     pos = XYZTLorentzVector(myTrack.trackerSurfacePosition().x(),
0256                             myTrack.trackerSurfacePosition().y(),
0257                             myTrack.trackerSurfacePosition().z(),
0258                             trackerSurfaceTime);
0259     mom = XYZTLorentzVector(myTrack.trackerSurfaceMomentum().x(),
0260                             myTrack.trackerSurfaceMomentum().y(),
0261                             myTrack.trackerSurfaceMomentum().z(),
0262                             myTrack.trackerSurfaceMomentum().t());
0263
0264     if (mom.T() > 0.) {
0265       // The particle to be propagated
0266       myPart = BaseParticlePropagator(RawParticle(mom, pos, myTrack.charge()), 0., 0., 4.);
0267
0268       // Propagate to Preshower layer 1
0269       myPart.propagateToPreshowerLayer1(false);
0270       if (myTrack.notYetToEndVertex(myPart.particle().vertex()) && myPart.getSuccess() > 0)
0271         myTrack.setLayer1(myPart.particle(), myPart.getSuccess());
0272
0273       // Propagate to Preshower Layer 2
0274       myPart.propagateToPreshowerLayer2(false);
0275       if (myTrack.notYetToEndVertex(myPart.particle().vertex()) && myPart.getSuccess() > 0)
0276         myTrack.setLayer2(myPart.particle(), myPart.getSuccess());
0277
0278       // Propagate to Ecal Endcap
0279       myPart.propagateToEcalEntrance(false);
0280       if (myTrack.notYetToEndVertex(myPart.particle().vertex()))
0281         myTrack.setEcal(myPart.particle(), myPart.getSuccess());
0282
0283       // Propagate to HCAL entrance
0284       myPart.propagateToHcalEntrance(false);
0285       if (myTrack.notYetToEndVertex(myPart.particle().vertex()))
0286         myTrack.setHcal(myPart.particle(), myPart.getSuccess());
0287
0288       // Attempt propagation to HF for low pt and high eta
0289       if (myPart.particle().cos2ThetaV() > 0.8 || mom.T() < 3.) {
0290         // Propagate to VFCAL entrance
0291         myPart.propagateToVFcalEntrance(false);
0292         if (myTrack.notYetToEndVertex(myPart.particle().vertex()))
0293           myTrack.setVFcal(myPart.particle(), myPart.getSuccess());
0294
0295         // Otherwise propagate to the HCAL exit and HO.
0296       } else {
0297         // Propagate to HCAL exit
0298         myPart.propagateToHcalExit(false);
0299         if (myTrack.notYetToEndVertex(myPart.particle().vertex()))
0300           myTrack.setHcalExit(myPart.particle(), myPart.getSuccess());
0301         // Propagate to HOLayer entrance
0302         myPart.setMagneticField(0);
0303         myPart.propagateToHOLayer(false);
0304         if (myTrack.notYetToEndVertex(myPart.particle().vertex()))
0305           myTrack.setHO(myPart.particle(), myPart.getSuccess());
0306       }
0307     }
0308   }
0309 }
0310
0311 void FBaseSimEvent::addParticles(const HepMC::GenEvent& myGenEvent) {
0312   /// Some internal array to work with.
0313   int genEventSize = myGenEvent.particles_size();
0314   std::vector<int> myGenVertices(genEventSize, static_cast<int>(0));
0315
0316   // If no particles, no work to be done !
0317   if (myGenEvent.particles_empty())
0318     return;
0319
0320   // Are there particles in the FSimEvent already ?
0321   int offset = nGenParts();
0322
0323   // Primary vertex
0324   HepMC::GenVertex* primaryVertex = *(myGenEvent.vertices_begin());
0325
0326   // unit transformation (needs review)
0327   XYZTLorentzVector primaryVertexPosition(primaryVertex->position().x() / 10.,
0328                                           primaryVertex->position().y() / 10.,
0329                                           primaryVertex->position().z() / 10.,
0330                                           primaryVertex->position().t() / 10.);
0331
0332   // This is the main vertex index
0333   int mainVertex = addSimVertex(primaryVertexPosition, -1, FSimVertexType::PRIMARY_VERTEX);
0334
0335   int initialBarcode = 0;
0336   if (myGenEvent.particles_begin() != myGenEvent.particles_end()) {
0337     initialBarcode = (*myGenEvent.particles_begin())->barcode();
0338   }
0339
0340   // Loop on the particles of the generated event
0341   for (auto piter = myGenEvent.particles_begin(); piter != myGenEvent.particles_end(); ++piter) {
0342     // This is the generated particle pointer - for the signal event only
0343     HepMC::GenParticle* p = *piter;
0344
0345     if (!offset) {
0346       (*theGenParticles)[nGenParticles++] = p;
0347       if (nGenParticles / theGenSize * theGenSize == nGenParticles) {
0348         theGenSize *= 2;
0349         theGenParticles->resize(theGenSize);
0350       }
0351     }
0352
0353     // Reject particles with late origin vertex (i.e., coming from late decays)
0354     // This should not happen, but one never knows what users may be up to!
0355     // For example exotic particles might decay late - keep the decay products in the case.
0356     XYZTLorentzVector productionVertexPosition(0., 0., 0., 0.);
0357     HepMC::GenVertex* productionVertex = p->production_vertex();
0358     if (productionVertex) {
0359       unsigned productionMother = productionVertex->particles_in_size();
0360       if (productionMother) {
0361         unsigned motherId = (*(productionVertex->particles_in_const_begin()))->pdg_id();
0362         if (motherId < 1000000)
0363           productionVertexPosition = XYZTLorentzVector(productionVertex->position().x() / 10.,
0364                                                        productionVertex->position().y() / 10.,
0365                                                        productionVertex->position().z() / 10.,
0366                                                        productionVertex->position().t() / 10.);
0367       }
0368     }
0369     if (!myFilter->acceptVertex(productionVertexPosition))
0370       continue;
0371
0372     int abspdgId = std::abs(p->pdg_id());
0373     HepMC::GenVertex* endVertex = p->end_vertex();
0374
0375     // Keep only:
0376     // 1) Stable particles (watch out! New status code = 1001!)
0377     bool testStable = p->status() % 1000 == 1;
0378     // Declare stable standard particles that decay after a macroscopic path length
0379     // (except if exotic)
0380     if (p->status() == 2 && abspdgId < 1000000) {
0381       if (endVertex) {
0382         XYZTLorentzVector decayPosition = XYZTLorentzVector(endVertex->position().x() / 10.,
0383                                                             endVertex->position().y() / 10.,
0384                                                             endVertex->position().z() / 10.,
0385                                                             endVertex->position().t() / 10.);
0386         // If the particle flew enough to be beyond the beam pipe enveloppe, just declare it stable
0387         if (decayPosition.Perp2() > lateVertexPosition)
0388           testStable = true;
0389       }
0390     }
0391
0392     // 2) or particles with stable daughters (watch out! New status code = 1001!)
0393     bool testDaugh = false;
0394     if (!testStable && p->status() == 2 && endVertex && endVertex->particles_out_size()) {
0395       HepMC::GenVertex::particles_out_const_iterator firstDaughterIt = endVertex->particles_out_const_begin();
0396       HepMC::GenVertex::particles_out_const_iterator lastDaughterIt = endVertex->particles_out_const_end();
0397       for (; firstDaughterIt != lastDaughterIt; ++firstDaughterIt) {
0398         HepMC::GenParticle* daugh = *firstDaughterIt;
0399         if (daugh->status() % 1000 == 1) {
0400           // Check that it is not a "prompt electron or muon brem":
0401           if (abspdgId == 11 || abspdgId == 13) {
0402             if (endVertex) {
0403               XYZTLorentzVector endVertexPosition = XYZTLorentzVector(endVertex->position().x() / 10.,
0404                                                                       endVertex->position().y() / 10.,
0405                                                                       endVertex->position().z() / 10.,
0406                                                                       endVertex->position().t() / 10.);
0407               // If the particle flew enough to be beyond the beam pipe enveloppe, just declare it stable
0408               if (endVertexPosition.Perp2() < lateVertexPosition) {
0409                 break;
0410               }
0411             }
0412           }
0413           testDaugh = true;
0414           break;
0415         }
0416       }
0417     }
0418
0419     // 3) or particles that fly more than one micron.
0420     double dist = 0.;
0421     if (!testStable && !testDaugh && p->production_vertex()) {
0422       XYZTLorentzVector productionVertexPosition(p->production_vertex()->position().x() / 10.,
0423                                                  p->production_vertex()->position().y() / 10.,
0424                                                  p->production_vertex()->position().z() / 10.,
0425                                                  p->production_vertex()->position().t() / 10.);
0426       dist = (primaryVertexPosition - productionVertexPosition).Vect().Mag2();
0427     }
0428     bool testDecay = (dist > 1e-8) ? true : false;
0429
0430     // Save the corresponding particle and vertices
0431     if (testStable || testDaugh || testDecay) {
0432       /*
0433       const HepMC::GenParticle* mother = p->production_vertex() ?
0434     *(p->production_vertex()->particles_in_const_begin()) : 0;
0435       */
0436
0437       int motherBarcode = p->production_vertex() && p->production_vertex()->particles_in_const_begin() !=
0438                                                         p->production_vertex()->particles_in_const_end()
0439                               ? (*(p->production_vertex()->particles_in_const_begin()))->barcode()
0440                               : 0;
0441
0442       int originVertex = motherBarcode && myGenVertices[motherBarcode - initialBarcode]
0443                              ? myGenVertices[motherBarcode - initialBarcode]
0444                              : mainVertex;
0445
0446       XYZTLorentzVector momentum(p->momentum().px(), p->momentum().py(), p->momentum().pz(), p->momentum().e());
0447       RawParticle part = makeParticle(theTable(), p->pdg_id(), momentum, vertex(originVertex).position());
0448
0449       // Add the particle to the event and to the various lists
0450
0451       int theTrack = testStable && p->end_vertex() ?
0452                                                    // The particle is scheduled to decay
0453                          addSimTrack(&part, originVertex, nGenParts() - offset, p->end_vertex())
0454                                                    :
0455                                                    // The particle is not scheduled to decay
0456                          addSimTrack(&part, originVertex, nGenParts() - offset);
0457
0458       if (
0459           // This one deals with particles with no end vertex
0460           !p->end_vertex() ||
0461           // This one deals with particles that have a pre-defined
0462           // decay proper time, but have not decayed yet
0463           (testStable && p->end_vertex() && !p->end_vertex()->particles_out_size())
0464           // In both case, just don't add a end vertex in the FSimEvent
0465       )
0466         continue;
0467
0468       // Add the vertex to the event and to the various lists
0469       XYZTLorentzVector decayVertex = XYZTLorentzVector(p->end_vertex()->position().x() / 10.,
0470                                                         p->end_vertex()->position().y() / 10.,
0471                                                         p->end_vertex()->position().z() / 10.,
0472                                                         p->end_vertex()->position().t() / 10.);
0473       //    vertex(mainVertex).position();
0474       int theVertex = addSimVertex(decayVertex, theTrack, FSimVertexType::DECAY_VERTEX);
0475
0476       if (theVertex != -1)
0477         myGenVertices[p->barcode() - initialBarcode] = theVertex;
0478
0479       // There we are !
0480     }
0481   }
0482 }
0483
0484 int FBaseSimEvent::addSimTrack(const RawParticle* p, int iv, int ig, const HepMC::GenVertex* ev) {
0485   // Check that the particle is in the Famos "acceptance"
0486   // Keep all primaries of pile-up events, though
0487   if (!myFilter->acceptParticle(*p) && ig >= -1)
0488     return -1;
0489
0490   // The new track index
0491   int trackId = nSimTracks++;
0492   if (nSimTracks / theTrackSize * theTrackSize == nSimTracks) {
0493     theTrackSize *= 2;
0494     theSimTracks->resize(theTrackSize);
0495   }
0496
0497   // Attach the particle to the origin vertex, and to the mother
0498   vertex(iv).addDaughter(trackId);
0499   if (!vertex(iv).noParent()) {
0500     track(vertex(iv).parent().id()).addDaughter(trackId);
0501
0502     if (ig == -1) {
0503       int motherId = track(vertex(iv).parent().id()).genpartIndex();
0504       if (motherId < -1)
0505         ig = motherId;
0506     }
0507   }
0508
0509   // Some transient information for FAMOS internal use
0510   (*theSimTracks)[trackId] =
0511       ev ?
0512          // A proper decay time is scheduled
0513           FSimTrack(p,
0514                     iv,
0515                     ig,
0516                     trackId,
0517                     this,
0518                     ev->position().t() / 10. * pdg::mass(p->pid(), theTable()) / std::sqrt(p->momentum().Vect().Mag2()))
0519          :
0520          // No proper decay time is scheduled
0521           FSimTrack(p, iv, ig, trackId, this);
0522
0523   return trackId;
0524 }
0525
0526 int FBaseSimEvent::addSimVertex(const XYZTLorentzVector& v, int im, FSimVertexType::VertexType type) {
0527   // Check that the vertex is in the Famos "acceptance"
0528   if (!myFilter->acceptVertex(v))
0529     return -1;
0530
0531   // The number of vertices
0532   int vertexId = nSimVertices++;
0533   if (nSimVertices / theVertexSize * theVertexSize == nSimVertices) {
0534     theVertexSize *= 2;
0535     theSimVertices->resize(theVertexSize);
0536     theFSimVerticesType->resize(theVertexSize);
0537   }
0538
0539   // Attach the end vertex to the particle (if accepted)
0540   if (im != -1)
0541     track(im).setEndVertex(vertexId);
0542
0543   // Some transient information for FAMOS internal use
0544   (*theSimVertices)[vertexId] = FSimVertex(v, im, vertexId, this);
0545
0546   (*theFSimVerticesType)[vertexId] = FSimVertexType(type);
0547
0548   return vertexId;
0549 }
0550
0551 void FBaseSimEvent::printMCTruth(const HepMC::GenEvent& myGenEvent) {
0552   std::cout << "Id  Gen Name       eta    phi     pT     E    Vtx1   "
0553             << " x      y      z   "
0554             << "Moth  Vtx2  eta   phi     R      Z   Da1  Da2 Ecal?" << std::endl;
0555
0556   for (HepMC::GenEvent::particle_const_iterator piter = myGenEvent.particles_begin();
0557        piter != myGenEvent.particles_end();
0558        ++piter) {
0559     HepMC::GenParticle* p = *piter;
0560     /* */
0561     int partId = p->pdg_id();
0562     std::string name;
0563
0564     if (pdt->particle(ParticleID(partId)) != nullptr) {
0565       name = (pdt->particle(ParticleID(partId)))->name();
0566     } else {
0567       name = "none";
0568     }
0569
0570     XYZTLorentzVector momentum1(p->momentum().px(), p->momentum().py(), p->momentum().pz(), p->momentum().e());
0571
0572     int vertexId1 = 0;
0573
0574     if (!p->production_vertex())
0575       continue;
0576
0577     XYZVector vertex1(p->production_vertex()->position().x() / 10.,
0578                       p->production_vertex()->position().y() / 10.,
0579                       p->production_vertex()->position().z() / 10.);
0580     vertexId1 = p->production_vertex()->barcode();
0581
0582     std::cout.setf(std::ios::fixed, std::ios::floatfield);
0583     std::cout.setf(std::ios::right, std::ios::adjustfield);
0584
0585     std::cout << std::setw(4) << p->barcode() << " " << name;
0586
0587     for (unsigned int k = 0; k < 11 - name.length() && k < 12; k++)
0588       std::cout << " ";
0589
0590     double eta = momentum1.eta();
0591     if (eta > +10.)
0592       eta = +10.;
0593     if (eta < -10.)
0594       eta = -10.;
0595     std::cout << std::setw(6) << std::setprecision(2) << eta << " " << std::setw(6) << std::setprecision(2)
0596               << momentum1.phi() << " " << std::setw(7) << std::setprecision(2) << momentum1.pt() << " " << std::setw(7)
0597               << std::setprecision(2) << momentum1.e() << " " << std::setw(4) << vertexId1 << " " << std::setw(6)
0598               << std::setprecision(1) << vertex1.x() << " " << std::setw(6) << std::setprecision(1) << vertex1.y()
0599               << " " << std::setw(6) << std::setprecision(1) << vertex1.z() << " ";
0600
0601     const HepMC::GenParticle* mother = *(p->production_vertex()->particles_in_const_begin());
0602
0603     if (mother)
0604       std::cout << std::setw(4) << mother->barcode() << " ";
0605     else
0606       std::cout << "     ";
0607
0608     if (p->end_vertex()) {
0609       XYZTLorentzVector vertex2(p->end_vertex()->position().x() / 10.,
0610                                 p->end_vertex()->position().y() / 10.,
0611                                 p->end_vertex()->position().z() / 10.,
0612                                 p->end_vertex()->position().t() / 10.);
0613       int vertexId2 = p->end_vertex()->barcode();
0614
0615       std::vector<const HepMC::GenParticle*> children;
0616       HepMC::GenVertex::particles_out_const_iterator firstDaughterIt = p->end_vertex()->particles_out_const_begin();
0617       HepMC::GenVertex::particles_out_const_iterator lastDaughterIt = p->end_vertex()->particles_out_const_end();
0618       for (; firstDaughterIt != lastDaughterIt; ++firstDaughterIt) {
0619         children.push_back(*firstDaughterIt);
0620       }
0621
0622       std::cout << std::setw(4) << vertexId2 << " " << std::setw(6) << std::setprecision(2) << vertex2.eta() << " "
0623                 << std::setw(6) << std::setprecision(2) << vertex2.phi() << " " << std::setw(5) << std::setprecision(1)
0624                 << vertex2.pt() << " " << std::setw(6) << std::setprecision(1) << vertex2.z() << " ";
0625       for (unsigned id = 0; id < children.size(); ++id)
0626         std::cout << std::setw(4) << children[id]->barcode() << " ";
0627     }
0628     std::cout << std::endl;
0629   }
0630 }
0631
0632 void FBaseSimEvent::print() const {
0633   std::cout << "  Id  Gen Name       eta    phi     pT     E    Vtx1   "
0634             << " x      y      z   "
0635             << "Moth  Vtx2  eta   phi     R      Z   Daughters Ecal?" << std::endl;
0636
0637   for (int i = 0; i < (int)nTracks(); i++)
0638     std::cout << track(i) << std::endl;
0639
0640   for (int i = 0; i < (int)nVertices(); i++)
0641     std::cout << "i = " << i << "  " << vertexType(i) << std::endl;
0642 }
0643
0644 void FBaseSimEvent::clear() {
0645   nSimTracks = 0;
0646   nSimVertices = 0;
0647   nGenParticles = 0;
0648   nChargedParticleTracks = 0;
0649 }
0650
0651 void FBaseSimEvent::addChargedTrack(int id) {
0652   (*theChargedTracks)[nChargedParticleTracks++] = id;
0653   if (nChargedParticleTracks / theChargedSize * theChargedSize == nChargedParticleTracks) {
0654     theChargedSize *= 2;
0655     theChargedTracks->resize(theChargedSize);
0656   }
0657 }
0658
0659 int FBaseSimEvent::chargedTrack(int id) const {
0660   if (id >= 0 && id < (int)nChargedParticleTracks)
0661     return (*theChargedTracks)[id];
0662   else
0663     return -1;
0664 }
0665
0666 const HepMC::GenParticle* FBaseSimEvent::embdGenpart(int i) const { return (*theGenParticles)[i]; }