CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/FastSimulation/Event/src/FBaseSimEvent.cc

Go to the documentation of this file.
00001 //HepMC Headers
00002 #include "HepMC/GenEvent.h"
00003 #include "HepMC/GenVertex.h"
00004 #include "HepMC/GenParticle.h"
00005 
00006 //Framework Headers
00007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00008 
00009 //CMSSW Data Formats
00010 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00011 #include "DataFormats/Candidate/interface/Candidate.h"
00012 
00013 //FAMOS Headers
00014 #include "FastSimulation/Event/interface/FBaseSimEvent.h"
00015 #include "FastSimulation/Event/interface/FSimTrack.h"
00016 #include "FastSimulation/Event/interface/FSimVertex.h"
00017 #include "FastSimulation/Event/interface/KineParticleFilter.h"
00018 #include "FastSimulation/BaseParticlePropagator/interface/BaseParticlePropagator.h"
00019 #include "FastSimulation/Event/interface/BetaFuncPrimaryVertexGenerator.h"
00020 #include "FastSimulation/Event/interface/GaussianPrimaryVertexGenerator.h"
00021 #include "FastSimulation/Event/interface/FlatPrimaryVertexGenerator.h"
00022 #include "FastSimulation/Event/interface/NoPrimaryVertexGenerator.h"
00023 
00024 #include "FastSimDataFormats/NuclearInteractions/interface/FSimVertexType.h"
00025 //#include "FastSimulation/Utilities/interface/Histos.h"
00026 
00027 using namespace HepPDT;
00028 
00029 // system include
00030 #include <iostream>
00031 #include <iomanip>
00032 #include <map>
00033 #include <string>
00034 
00035 FBaseSimEvent::FBaseSimEvent(const edm::ParameterSet& kine) 
00036   :
00037   nSimTracks(0),
00038   nSimVertices(0),
00039   nGenParticles(0),
00040   nChargedParticleTracks(0),
00041   initialSize(5000),
00042   random(0)
00043 {
00044 
00045   theVertexGenerator = new NoPrimaryVertexGenerator();
00046   theBeamSpot = math::XYZPoint(0.0,0.0,0.0);
00047 
00048   // Initialize the vectors of particles and vertices
00049   theGenParticles = new std::vector<HepMC::GenParticle*>(); 
00050   theSimTracks = new std::vector<FSimTrack>;
00051   theSimVertices = new std::vector<FSimVertex>;
00052   theChargedTracks = new std::vector<unsigned>();
00053   theFSimVerticesType = new FSimVertexTypeCollection();
00054 
00055   // Reserve some size to avoid mutiple copies
00056   /* */
00057   theSimTracks->resize(initialSize);
00058   theSimVertices->resize(initialSize);
00059   theGenParticles->resize(initialSize);
00060   theChargedTracks->resize(initialSize);
00061   theFSimVerticesType->resize(initialSize);
00062   theTrackSize = initialSize;
00063   theVertexSize = initialSize;
00064   theGenSize = initialSize;
00065   theChargedSize = initialSize;
00066   /* */
00067 
00068   // Initialize the Particle filter
00069   myFilter = new KineParticleFilter(kine);
00070 
00071 }
00072 
00073 FBaseSimEvent::FBaseSimEvent(const edm::ParameterSet& vtx,
00074                              const edm::ParameterSet& kine,
00075                              const RandomEngine* engine) 
00076   :
00077   nSimTracks(0),
00078   nSimVertices(0),
00079   nGenParticles(0),
00080   nChargedParticleTracks(0), 
00081   initialSize(5000),
00082   theVertexGenerator(0), 
00083   random(engine)
00084 {
00085 
00086   // Initialize the vertex generator
00087   std::string vtxType = vtx.getParameter<std::string>("type");
00088   if ( vtxType == "Gaussian" ) 
00089     theVertexGenerator = new GaussianPrimaryVertexGenerator(vtx,random);
00090   else if ( vtxType == "Flat" ) 
00091     theVertexGenerator = new FlatPrimaryVertexGenerator(vtx,random);
00092   else if ( vtxType == "BetaFunc" )
00093     theVertexGenerator = new BetaFuncPrimaryVertexGenerator(vtx,random);
00094   else
00095     theVertexGenerator = new NoPrimaryVertexGenerator();
00096   // Initialize the beam spot, if not read from the DataBase
00097   theBeamSpot = math::XYZPoint(0.0,0.0,0.0);
00098 
00099   // Initialize the distance from (0,0,0) after which *generated* particles are 
00100   // no longer considered - because the mother could have interacted before.
00101   // unit : cm x cm
00102   lateVertexPosition = 2.5*2.5;
00103 
00104   // Initialize the vectors of particles and vertices
00105   theGenParticles = new std::vector<HepMC::GenParticle*>(); 
00106   theSimTracks = new std::vector<FSimTrack>;
00107   theSimVertices = new std::vector<FSimVertex>;
00108   theChargedTracks = new std::vector<unsigned>();
00109   theFSimVerticesType = new FSimVertexTypeCollection();
00110 
00111   // Reserve some size to avoid mutiple copies
00112   /* */
00113   theSimTracks->resize(initialSize);
00114   theSimVertices->resize(initialSize);
00115   theGenParticles->resize(initialSize);
00116   theChargedTracks->resize(initialSize);
00117   theFSimVerticesType->resize(initialSize);
00118   theTrackSize = initialSize;
00119   theVertexSize = initialSize;
00120   theGenSize = initialSize;
00121   theChargedSize = initialSize;
00122   /* */
00123 
00124   // Initialize the Particle filter
00125   myFilter = new KineParticleFilter(kine);
00126 
00127   // Get the Famos Histos pointer
00128   //  myHistos = Histos::instance();
00129 
00130   // Initialize a few histograms
00131   /*
00132   myHistos->book("hvtx",100,-0.1,0.1);
00133   myHistos->book("hvty",100,-0.1,0.1);
00134   myHistos->book("hvtz",100,-500.,500.);
00135   */
00136 }
00137  
00138 FBaseSimEvent::~FBaseSimEvent(){
00139 
00140   // Clear the vectors
00141   theGenParticles->clear();
00142   theSimTracks->clear();
00143   theSimVertices->clear();
00144   theChargedTracks->clear();
00145   theFSimVerticesType->clear();
00146 
00147   // Delete 
00148   delete theGenParticles;
00149   delete theSimTracks;
00150   delete theSimVertices;
00151   delete theChargedTracks;
00152   delete theFSimVerticesType;
00153   delete myFilter;
00154 
00155   //Write the histograms
00156   //  myHistos->put("histos.root");
00157   //  delete myHistos;
00158 }
00159 
00160 void 
00161 FBaseSimEvent::initializePdt(const HepPDT::ParticleDataTable* aPdt) { 
00162 
00163   pdt = aPdt; 
00164 
00165 }
00166 
00167 /*
00168 const HepPDT::ParticleDataTable*
00169 FBaseSimEvent::theTable() const {
00170   return pdt;
00171 }
00172 */
00173 
00174 void
00175 FBaseSimEvent::fill(const HepMC::GenEvent& myGenEvent) {
00176   
00177   // Clear old vectors
00178   clear();
00179 
00180   // Add the particles in the FSimEvent
00181   addParticles(myGenEvent);
00182 
00183   /*
00184   std::cout << "The MC truth! " << std::endl;
00185   printMCTruth(myGenEvent);
00186 
00187   std::cout << std::endl  << "The FAMOS event! " << std::endl;
00188   print();
00189   */
00190 
00191 }
00192 
00193 void
00194 FBaseSimEvent::fill(const reco::GenParticleCollection& myGenParticles) {
00195   
00196   // Clear old vectors
00197   clear();
00198 
00199   // Add the particles in the FSimEvent
00200   addParticles(myGenParticles);
00201 
00202 }
00203 
00204 void
00205 FBaseSimEvent::fill(const std::vector<SimTrack>& simTracks, 
00206                     const std::vector<SimVertex>& simVertices) {
00207 
00208   // Watch out there ! A SimVertex is in mm (stupid), 
00209   //            while a FSimVertex is in cm (clever).
00210   
00211   clear();
00212 
00213   unsigned nVtx = simVertices.size();
00214   unsigned nTks = simTracks.size();
00215 
00216   // Empty event, do nothin'
00217   if ( nVtx == 0 ) return;
00218 
00219   // Two arrays for internal use.
00220   std::vector<int> myVertices(nVtx,-1);
00221   std::vector<int> myTracks(nTks,-1);
00222 
00223   // create a map associating geant particle id and position in the 
00224   // event SimTrack vector
00225   
00226   std::map<unsigned, unsigned> geantToIndex;
00227   for( unsigned it=0; it<simTracks.size(); ++it ) {
00228     geantToIndex[ simTracks[it].trackId() ] = it;
00229   }  
00230 
00231   // Create also a map associating a SimTrack with its endVertex
00232   /*
00233   std::map<unsigned, unsigned> endVertex;
00234   for ( unsigned iv=0; iv<simVertices.size(); ++iv ) { 
00235     endVertex[ simVertices[iv].parentIndex() ] = iv;
00236   }
00237   */
00238 
00239   // Set the main vertex for the kine particle filter
00240   // SimVertices were in mm until 110_pre2
00241   //  HepLorentzVector primaryVertex = simVertices[0].position()/10.;
00242   // SImVertices are now in cm
00243   // Also : position is copied until SimVertex switches to Mathcore.
00244   //  XYZTLorentzVector primaryVertex = simVertices[0].position();
00245   // The next 5 lines to be then replaced by the previous line
00246   XYZTLorentzVector primaryVertex(simVertices[0].position().x(),
00247                                   simVertices[0].position().y(),
00248                                   simVertices[0].position().z(),
00249                                   simVertices[0].position().t());
00250   //
00251   myFilter->setMainVertex(primaryVertex);
00252   // Add the main vertex to the list.
00253   addSimVertex(myFilter->vertex(), -1, FSimVertexType::PRIMARY_VERTEX);
00254   myVertices[0] = 0;
00255 
00256   for( unsigned trackId=0; trackId<nTks; ++trackId ) {
00257 
00258     // The track
00259     const SimTrack& track = simTracks[trackId];
00260     //    std::cout << std::endl << "SimTrack " << trackId << " " << track << std::endl;
00261 
00262     // The origin vertex
00263     int vertexId = track.vertIndex();
00264     const SimVertex& vertex = simVertices[vertexId];
00265     //std::cout << "Origin vertex " << vertexId << " " << vertex << std::endl;
00266 
00267     // The mother track 
00268     int motherId = -1;
00269     if( !vertex.noParent() ) { // there is a parent to this vertex
00270       // geant id of the mother
00271       unsigned motherGeantId =   vertex.parentIndex(); 
00272       std::map<unsigned, unsigned >::iterator association  
00273         = geantToIndex.find( motherGeantId );
00274       if(association != geantToIndex.end() )
00275         motherId = association->second;
00276     }
00277     int originId = motherId == - 1 ? -1 : myTracks[motherId];
00278     //std::cout << "Origin id " << originId << std::endl;
00279 
00280     /*
00281     if ( endVertex.find(trackId) != endVertex.end() ) 
00282       std::cout << "End vertex id = " << endVertex[trackId] << std::endl;
00283     else
00284       std::cout << "No endVertex !!! " << std::endl;
00285     std::cout << "Tracker surface position " << track.trackerSurfacePosition() << std::endl;
00286     */
00287 
00288     // Add the vertex (if it does not already exist!)
00289     XYZTLorentzVector position(vertex.position().px(),vertex.position().py(),
00290                                vertex.position().pz(),vertex.position().e());
00291     if ( myVertices[vertexId] == -1 ) 
00292       // Momentum and position are copied until SimTrack and SimVertex
00293       // switch to Mathcore.
00294       //      myVertices[vertexId] = addSimVertex(vertex.position(),originId); 
00295       // The next line to be then replaced by the previous line
00296       myVertices[vertexId] = addSimVertex(position,originId); 
00297 
00298     // Add the track (with protection for brem'ing electrons and muons)
00299     int motherType = motherId == -1 ? 0 : simTracks[motherId].type();
00300 
00301     bool notBremInDetector =
00302       (abs(motherType) != 11 && abs(motherType) != 13) ||
00303       motherType != track.type() ||
00304       position.Perp2() < lateVertexPosition ;
00305 
00306     if ( notBremInDetector ) {
00307       // Momentum and position are copied until SimTrack and SimVertex
00308       // switch to Mathcore.
00309       //      RawParticle part(track.momentum(), vertex.position());
00310       // The next 3 lines to be then replaced by the previous line
00311       XYZTLorentzVector momentum(track.momentum().px(),track.momentum().py(),
00312                                  track.momentum().pz(),track.momentum().e());
00313       RawParticle part(momentum,position);
00314       //
00315       part.setID(track.type()); 
00316       //std::cout << "Ctau  = " << part.PDGcTau() << std::endl;
00317       // Don't save tracks that have decayed immediately but for which no daughters
00318       // were saved (probably due to cuts on E, pT and eta)
00319       //  if ( part.PDGcTau() > 0.1 || endVertex.find(trackId) != endVertex.end() ) 
00320       myTracks[trackId] = addSimTrack(&part,myVertices[vertexId],track.genpartIndex());
00321       if ( myTracks[trackId] >= 0 ) { 
00322         (*theSimTracks)[ myTracks[trackId] ].setTkPosition(track.trackerSurfacePosition());
00323         (*theSimTracks)[ myTracks[trackId] ].setTkMomentum(track.trackerSurfaceMomentum());
00324       }
00325     } else {
00326 
00327       myTracks[trackId] = myTracks[motherId];
00328       if ( myTracks[trackId] >= 0 ) { 
00329         (*theSimTracks)[ myTracks[trackId] ].setTkPosition(track.trackerSurfacePosition());
00330         (*theSimTracks)[ myTracks[trackId] ].setTkMomentum(track.trackerSurfaceMomentum());
00331       }
00332     }
00333     
00334   }
00335 
00336   // Now loop over the remaining end vertices !
00337   for( unsigned vertexId=0; vertexId<nVtx; ++vertexId ) {
00338 
00339     // if the vertex is already saved, just ignore.
00340     if ( myVertices[vertexId] != -1 ) continue;
00341 
00342     // The yet unused vertex
00343     const SimVertex& vertex = simVertices[vertexId];
00344 
00345     // The mother track 
00346     int motherId = -1;
00347     if( !vertex.noParent() ) { // there is a parent to this vertex
00348 
00349       // geant id of the mother
00350       unsigned motherGeantId =   vertex.parentIndex(); 
00351       std::map<unsigned, unsigned >::iterator association  
00352         = geantToIndex.find( motherGeantId );
00353       if(association != geantToIndex.end() )
00354         motherId = association->second;
00355     }
00356     int originId = motherId == - 1 ? -1 : myTracks[motherId];
00357 
00358     // Add the vertex
00359     // Momentum and position are copied until SimTrack and SimVertex
00360     // switch to Mathcore.
00361     //    myVertices[vertexId] = addSimVertex(vertex.position(),originId);
00362     // The next 3 lines to be then replaced by the previous line
00363     XYZTLorentzVector position(vertex.position().px(),vertex.position().py(),
00364                                vertex.position().pz(),vertex.position().e());
00365     myVertices[vertexId] = addSimVertex(position,originId); 
00366   }
00367 
00368   // Finally, propagate all particles to the calorimeters
00369   BaseParticlePropagator myPart;
00370   XYZTLorentzVector mom;
00371   XYZTLorentzVector pos;
00372 
00373   // Loop over the tracks
00374   for( int fsimi=0; fsimi < (int)nTracks() ; ++fsimi) {
00375 
00376     
00377     FSimTrack& myTrack = track(fsimi);
00378     double trackerSurfaceTime = myTrack.vertex().position().t() 
00379                               + myTrack.momentum().e()/myTrack.momentum().pz()
00380                               * ( myTrack.trackerSurfacePosition().z()
00381                                 - myTrack.vertex().position().z() );
00382     pos = XYZTLorentzVector(myTrack.trackerSurfacePosition().x(),
00383                             myTrack.trackerSurfacePosition().y(),
00384                             myTrack.trackerSurfacePosition().z(),
00385                                     trackerSurfaceTime);
00386     mom = XYZTLorentzVector(myTrack.trackerSurfaceMomentum().x(),
00387                             myTrack.trackerSurfaceMomentum().y(),
00388                             myTrack.trackerSurfaceMomentum().z(),
00389                             myTrack.trackerSurfaceMomentum().t());
00390 
00391     if ( mom.T() >  0. ) {  
00392       // The particle to be propagated
00393       myPart = BaseParticlePropagator(RawParticle(mom,pos),0.,0.,4.);
00394       myPart.setCharge(myTrack.charge());
00395       
00396       // Propagate to Preshower layer 1
00397       myPart.propagateToPreshowerLayer1(false);
00398       if ( myTrack.notYetToEndVertex(myPart.vertex()) && myPart.getSuccess()>0 )
00399         myTrack.setLayer1(myPart,myPart.getSuccess());
00400       
00401       // Propagate to Preshower Layer 2 
00402       myPart.propagateToPreshowerLayer2(false);
00403       if ( myTrack.notYetToEndVertex(myPart.vertex()) && myPart.getSuccess()>0 )
00404         myTrack.setLayer2(myPart,myPart.getSuccess());
00405       
00406       // Propagate to Ecal Endcap
00407       myPart.propagateToEcalEntrance(false);
00408       if ( myTrack.notYetToEndVertex(myPart.vertex()) )
00409         myTrack.setEcal(myPart,myPart.getSuccess());
00410       
00411       // Propagate to HCAL entrance
00412       myPart.propagateToHcalEntrance(false);
00413       if ( myTrack.notYetToEndVertex(myPart.vertex()) )
00414         myTrack.setHcal(myPart,myPart.getSuccess());
00415       
00416       // Propagate to VFCAL entrance
00417       myPart.propagateToVFcalEntrance(false);
00418       if ( myTrack.notYetToEndVertex(myPart.vertex()) )
00419         myTrack.setVFcal(myPart,myPart.getSuccess());
00420     }
00421  
00422   }
00423 
00424 }
00425 
00426 
00427 void
00428 FBaseSimEvent::addParticles(const HepMC::GenEvent& myGenEvent) {
00429 
00431   int genEventSize = myGenEvent.particles_size();
00432   std::vector<int> myGenVertices(genEventSize, static_cast<int>(0));
00433 
00434   // If no particles, no work to be done !
00435   if ( myGenEvent.particles_empty() ) return;
00436 
00437   // Are there particles in the FSimEvent already ? 
00438   int offset = nGenParts();
00439 
00440   // Primary vertex (already smeared by the SmearedVtx module)
00441   HepMC::GenVertex* primaryVertex = *(myGenEvent.vertices_begin());
00442 
00443   // Beginning of workaround a bug in pythia particle gun
00444   unsigned primaryMother = primaryVertex->particles_in_size();
00445   if ( primaryMother ) {
00446     unsigned partId = (*(primaryVertex->particles_in_const_begin()))->pdg_id();
00447     if ( abs(partId) == 2212 ) primaryMother = 0;
00448   }
00449   // End of workaround a bug in pythia particle gun
00450 
00451   XYZTLorentzVector primaryVertexPosition(primaryVertex->position().x()/10.,
00452                                           primaryVertex->position().y()/10.,
00453                                           primaryVertex->position().z()/10.,
00454                                           primaryVertex->position().t()/10.);
00455   // Actually this is the true end of the workaround
00456   primaryVertexPosition *= (1-primaryMother);
00457   // THE END.
00458 
00459   // Smear the main vertex if needed
00460   // Now takes the origin from the database
00461   XYZTLorentzVector smearedVertex; 
00462   if ( primaryVertexPosition.Vect().Mag2() < 1E-16 ) {
00463     theVertexGenerator->generate();
00464     smearedVertex = XYZTLorentzVector(
00465       theVertexGenerator->X()-theVertexGenerator->beamSpot().X()+theBeamSpot.X(),
00466       theVertexGenerator->Y()-theVertexGenerator->beamSpot().Y()+theBeamSpot.Y(),
00467       theVertexGenerator->Z()-theVertexGenerator->beamSpot().Z()+theBeamSpot.Z(),
00468       0.);
00469   }
00470 
00471   // Set the main vertex
00472   myFilter->setMainVertex(primaryVertexPosition+smearedVertex);
00473 
00474   // This is the smeared main vertex
00475   int mainVertex = addSimVertex(myFilter->vertex(), -1, FSimVertexType::PRIMARY_VERTEX);
00476 
00477   HepMC::GenEvent::particle_const_iterator piter;
00478   HepMC::GenEvent::particle_const_iterator pbegin = myGenEvent.particles_begin();
00479   HepMC::GenEvent::particle_const_iterator pend = myGenEvent.particles_end();
00480 
00481   int initialBarcode = 0; 
00482   if ( pbegin != pend ) initialBarcode = (*pbegin)->barcode();
00483   // Loop on the particles of the generated event
00484   for ( piter = pbegin; piter != pend; ++piter ) {
00485 
00486     // This is the generated particle pointer - for the signal event only
00487     HepMC::GenParticle* p = *piter;
00488 
00489     if  ( !offset ) {
00490       (*theGenParticles)[nGenParticles++] = p;
00491       if ( nGenParticles/theGenSize*theGenSize == nGenParticles ) { 
00492         theGenSize *= 2;
00493         theGenParticles->resize(theGenSize);
00494       }
00495 
00496     }
00497 
00498     // Reject particles with late origin vertex (i.e., coming from late decays)
00499     // This should not happen, but one never knows what users may be up to!
00500     // For example exotic particles might decay late - keep the decay products in the case.
00501     XYZTLorentzVector productionVertexPosition(0.,0.,0.,0.);
00502     HepMC::GenVertex* productionVertex = p->production_vertex();
00503     if ( productionVertex ) { 
00504       unsigned productionMother = productionVertex->particles_in_size();
00505       if ( productionMother ) {
00506         unsigned motherId = (*(productionVertex->particles_in_const_begin()))->pdg_id();
00507         if ( abs(motherId) < 1000000 ) 
00508           productionVertexPosition = 
00509             XYZTLorentzVector(productionVertex->position().x()/10.,
00510                               productionVertex->position().y()/10.,
00511                               productionVertex->position().z()/10.,
00512                               productionVertex->position().t()/10.) + smearedVertex;
00513       }
00514     }
00515     if ( !myFilter->accept(productionVertexPosition) ) continue;
00516 
00517     int abspdgId = abs(p->pdg_id());
00518     HepMC::GenVertex* endVertex = p->end_vertex();
00519 
00520     // Keep only: 
00521     // 1) Stable particles (watch out! New status code = 1001!)
00522     bool testStable = p->status()%1000==1;
00523     // Declare stable standard particles that decay after a macroscopic path length
00524     // (except if exotic)
00525     if ( p->status() == 2 && abspdgId < 1000000) {
00526       if ( endVertex ) { 
00527         XYZTLorentzVector decayPosition = 
00528           XYZTLorentzVector(endVertex->position().x()/10.,
00529                             endVertex->position().y()/10.,
00530                             endVertex->position().z()/10.,
00531                             endVertex->position().t()/10.) + smearedVertex;
00532         // If the particle flew enough to be beyond the beam pipe enveloppe, just declare it stable
00533         if ( decayPosition.Perp2() > lateVertexPosition ) testStable = true;
00534       }
00535     }      
00536 
00537     // 2) or particles with stable daughters (watch out! New status code = 1001!)
00538     bool testDaugh = false;
00539     if ( !testStable && 
00540          p->status() == 2 &&
00541          endVertex && 
00542          endVertex->particles_out_size() ) { 
00543       HepMC::GenVertex::particles_out_const_iterator firstDaughterIt = 
00544         endVertex->particles_out_const_begin();
00545       HepMC::GenVertex::particles_out_const_iterator lastDaughterIt = 
00546         endVertex->particles_out_const_end();
00547       for ( ; firstDaughterIt != lastDaughterIt ; ++firstDaughterIt ) {
00548         HepMC::GenParticle* daugh = *firstDaughterIt;
00549         if ( daugh->status()%1000==1 ) {
00550           // Check that it is not a "prompt electron or muon brem":
00551           if (abspdgId == 11 || abspdgId == 13) {
00552             if ( endVertex ) { 
00553               XYZTLorentzVector endVertexPosition = XYZTLorentzVector(endVertex->position().x()/10.,
00554                                                                       endVertex->position().y()/10.,
00555                                                                       endVertex->position().z()/10.,
00556                                                                       endVertex->position().t()/10.);
00557               // If the particle flew enough to be beyond the beam pipe enveloppe, just declare it stable
00558               if ( endVertexPosition.Perp2() < lateVertexPosition ) {
00559                 break;
00560               }
00561             }
00562           }
00563           testDaugh=true;
00564           break;
00565         }
00566       }
00567     }
00568 
00569     // 3) or particles that fly more than one micron.
00570     double dist = 0.;
00571     if ( !testStable && !testDaugh && p->production_vertex() ) {
00572       XYZTLorentzVector 
00573         productionVertexPosition(p->production_vertex()->position().x()/10.,
00574                                  p->production_vertex()->position().y()/10.,
00575                                  p->production_vertex()->position().z()/10.,
00576                                  p->production_vertex()->position().t()/10.);
00577       dist = (primaryVertexPosition-productionVertexPosition).Vect().Mag2();
00578     }
00579     bool testDecay = ( dist > 1e-8 ) ? true : false; 
00580 
00581     // Save the corresponding particle and vertices
00582     if ( testStable || testDaugh || testDecay ) {
00583       
00584       /*
00585       const HepMC::GenParticle* mother = p->production_vertex() ?
00586         *(p->production_vertex()->particles_in_const_begin()) : 0;
00587       */
00588 
00589       int motherBarcode = p->production_vertex() && 
00590         p->production_vertex()->particles_in_const_begin() !=
00591         p->production_vertex()->particles_in_const_end() ?
00592         (*(p->production_vertex()->particles_in_const_begin()))->barcode() : 0;
00593 
00594       int originVertex = 
00595         motherBarcode && myGenVertices[motherBarcode-initialBarcode] ?
00596         myGenVertices[motherBarcode-initialBarcode] : mainVertex;
00597 
00598       XYZTLorentzVector momentum(p->momentum().px(),
00599                                  p->momentum().py(),
00600                                  p->momentum().pz(),
00601                                  p->momentum().e());
00602       RawParticle part(momentum, vertex(originVertex).position());
00603       part.setID(p->pdg_id());
00604 
00605       // Add the particle to the event and to the various lists
00606       
00607       int theTrack = testStable && p->end_vertex() ? 
00608         // The particle is scheduled to decay
00609         addSimTrack(&part,originVertex, nGenParts()-offset,p->end_vertex()) :
00610         // The particle is not scheduled to decay 
00611         addSimTrack(&part,originVertex, nGenParts()-offset);
00612 
00613       if ( 
00614           // This one deals with particles with no end vertex
00615           !p->end_vertex() ||
00616           // This one deals with particles that have a pre-defined
00617           // decay proper time, but have not decayed yet
00618           ( testStable && p->end_vertex() &&  !p->end_vertex()->particles_out_size() ) 
00619           // In both case, just don't add a end vertex in the FSimEvent 
00620           ) continue; 
00621       
00622       // Add the vertex to the event and to the various lists
00623       XYZTLorentzVector decayVertex = 
00624         XYZTLorentzVector(p->end_vertex()->position().x()/10.,
00625                           p->end_vertex()->position().y()/10.,
00626                           p->end_vertex()->position().z()/10.,
00627                           p->end_vertex()->position().t()/10.) +
00628         smearedVertex;
00629       //        vertex(mainVertex).position();
00630       int theVertex = addSimVertex(decayVertex,theTrack, FSimVertexType::DECAY_VERTEX);
00631 
00632       if ( theVertex != -1 ) myGenVertices[p->barcode()-initialBarcode] = theVertex;
00633 
00634       // There we are !
00635     }
00636   }
00637 
00638 }
00639 
00640 void
00641 FBaseSimEvent::addParticles(const reco::GenParticleCollection& myGenParticles) {
00642 
00643   // If no particles, no work to be done !
00644   unsigned int nParticles = myGenParticles.size();
00645   nGenParticles = nParticles;
00646 
00647   if ( !nParticles ) return;
00648 
00650   std::map<const reco::Candidate*,int> myGenVertices;
00651 
00652   // Are there particles in the FSimEvent already ? 
00653   int offset = nTracks();
00654 
00655   // Skip the incoming protons
00656   nGenParticles = 0;
00657   unsigned int ip = 0;
00658   if ( nParticles > 1 && 
00659        myGenParticles[0].pdgId() == 2212 &&
00660        myGenParticles[1].pdgId() == 2212 ) { 
00661     ip = 2;
00662     nGenParticles = 2;
00663   }
00664 
00665   // Primary vertex (already smeared by the SmearedVtx module)
00666   XYZTLorentzVector primaryVertex (myGenParticles[ip].vx(),
00667                                    myGenParticles[ip].vy(),
00668                                    myGenParticles[ip].vz(),
00669                                    0.);
00670 
00671   // Smear the main vertex if needed
00672   XYZTLorentzVector smearedVertex;
00673   if ( primaryVertex.mag() < 1E-8 ) {
00674     theVertexGenerator->generate();
00675     smearedVertex = XYZTLorentzVector(
00676       theVertexGenerator->X()-theVertexGenerator->beamSpot().X()+theBeamSpot.X(),
00677       theVertexGenerator->Y()-theVertexGenerator->beamSpot().Y()+theBeamSpot.Y(),
00678       theVertexGenerator->Z()-theVertexGenerator->beamSpot().Z()+theBeamSpot.Z(),
00679       0.);
00680   }
00681 
00682   // Set the main vertex
00683   myFilter->setMainVertex(primaryVertex+smearedVertex);
00684 
00685   // This is the smeared main vertex
00686   int mainVertex = addSimVertex(myFilter->vertex(), -1, FSimVertexType::PRIMARY_VERTEX);
00687 
00688   // Loop on the particles of the generated event
00689   for ( ; ip<nParticles; ++ip ) { 
00690     
00691     // nGenParticles = ip;
00692     
00693     nGenParticles++;
00694     const reco::GenParticle& p = myGenParticles[ip];
00695 
00696     // Reject particles with late origin vertex (i.e., coming from late decays)
00697     // This should not happen, but one never knows what users may be up to!
00698     // For example exotic particles might decay late - keep the decay products in the case.
00699     XYZTLorentzVector productionVertexPosition(0.,0.,0.,0.);
00700     const reco::Candidate* productionMother = p.numberOfMothers() ? p.mother(0) : 0;
00701     if ( productionMother ) {
00702       unsigned motherId = productionMother->pdgId();
00703       if ( abs(motherId) < 1000000 )
00704         productionVertexPosition = XYZTLorentzVector(p.vx(), p.vy(), p.vz(), 0.) + smearedVertex;
00705     }
00706     if ( !myFilter->accept(productionVertexPosition) ) continue;
00707 
00708     // Keep only: 
00709     // 1) Stable particles
00710     bool testStable = p.status()%1000==1;
00711     // Declare stable standard particles that decay after a macroscopic path length 
00712     // (except if exotic particle)
00713     if ( p.status() == 2 && abs(p.pdgId()) < 1000000 ) {
00714       unsigned int nDaughters = p.numberOfDaughters();
00715       if ( nDaughters ) { 
00716         const reco::Candidate* daughter = p.daughter(0);
00717         XYZTLorentzVector decayPosition = 
00718           XYZTLorentzVector(daughter->vx(), daughter->vy(), daughter->vz(), 0.) + smearedVertex;
00719         // If the particle flew enough to be beyond the beam pipe enveloppe, just declare it stable
00720         if ( decayPosition.Perp2() > lateVertexPosition ) testStable = true;
00721       }
00722     }
00723 
00724     // 2) or particles with stable daughters
00725     bool testDaugh = false;
00726     unsigned int nDaughters = p.numberOfDaughters();
00727     if ( !testStable  && 
00728          //      p.status() == 2 && 
00729          nDaughters ) {  
00730       for ( unsigned iDaughter=0; iDaughter<nDaughters; ++iDaughter ) {
00731         const reco::Candidate* daughter = p.daughter(iDaughter);
00732         if ( daughter->status()%1000==1 ) {
00733           testDaugh=true;
00734           break;
00735         }
00736       }
00737     }
00738 
00739     // 3) or particles that fly more than one micron.
00740     double dist = 0.;
00741     if ( !testStable && !testDaugh ) {
00742       XYZTLorentzVector productionVertex(p.vx(),p.vy(),p.vz(),0.);
00743       dist = (primaryVertex-productionVertex).Vect().Mag2();
00744     }
00745     bool testDecay = ( dist > 1e-8 ) ? true : false; 
00746 
00747     // Save the corresponding particle and vertices
00748     if ( testStable || testDaugh || testDecay ) {
00749       
00750       const reco::Candidate* mother = p.numberOfMothers() ? p.mother(0) : 0;
00751 
00752       int originVertex = 
00753         mother &&  
00754         myGenVertices.find(mother) != myGenVertices.end() ? 
00755         myGenVertices[mother] : mainVertex;
00756       
00757       XYZTLorentzVector momentum(p.px(),p.py(),p.pz(),p.energy());
00758       RawParticle part(momentum, vertex(originVertex).position());
00759       part.setID(p.pdgId());
00760 
00761       // Add the particle to the event and to the various lists
00762       int theTrack = addSimTrack(&part,originVertex, nGenParts()-offset);
00763 
00764       // It there an end vertex ?
00765       if ( !nDaughters ) continue; 
00766       const reco::Candidate* daughter = p.daughter(0);
00767 
00768       // Add the vertex to the event and to the various lists
00769       XYZTLorentzVector decayVertex = 
00770         XYZTLorentzVector(daughter->vx(), daughter->vy(),
00771                           daughter->vz(), 0.) + smearedVertex;
00772       int theVertex = addSimVertex(decayVertex,theTrack, FSimVertexType::DECAY_VERTEX);
00773 
00774       if ( theVertex != -1 ) myGenVertices[&p] = theVertex;
00775 
00776       // There we are !
00777     }
00778   }
00779 
00780   // There is no GenParticle's in that case...
00781   // nGenParticles=0;
00782 
00783 }
00784 
00785 int 
00786 FBaseSimEvent::addSimTrack(const RawParticle* p, int iv, int ig, 
00787                            const HepMC::GenVertex* ev) { 
00788   
00789   // Check that the particle is in the Famos "acceptance"
00790   // Keep all primaries of pile-up events, though
00791   if ( !myFilter->accept(p) && ig >= -1 ) return -1;
00792 
00793   // The new track index
00794   int trackId = nSimTracks++;
00795   if ( nSimTracks/theTrackSize*theTrackSize == nSimTracks ) {
00796     theTrackSize *= 2;
00797     theSimTracks->resize(theTrackSize);
00798   }
00799 
00800   // Attach the particle to the origin vertex, and to the mother
00801   vertex(iv).addDaughter(trackId);
00802   if ( !vertex(iv).noParent() ) {
00803     track(vertex(iv).parent().id()).addDaughter(trackId);
00804 
00805     if ( ig == -1 ) {
00806       int motherId = track(vertex(iv).parent().id()).genpartIndex();
00807       if ( motherId < -1 ) ig = motherId;
00808     }
00809   }
00810     
00811   // Some transient information for FAMOS internal use
00812   (*theSimTracks)[trackId] = ev ? 
00813     // A proper decay time is scheduled
00814     FSimTrack(p,iv,ig,trackId,this,
00815               ev->position().t()/10.
00816               * p->PDGmass()
00817               / std::sqrt(p->momentum().Vect().Mag2())) : 
00818     // No proper decay time is scheduled
00819     FSimTrack(p,iv,ig,trackId,this);
00820 
00821   return trackId;
00822 
00823 }
00824 
00825 int
00826 FBaseSimEvent::addSimVertex(const XYZTLorentzVector& v, int im, FSimVertexType::VertexType type) {
00827   
00828   // Check that the vertex is in the Famos "acceptance"
00829   if ( !myFilter->accept(v) ) return -1;
00830 
00831   // The number of vertices
00832   int vertexId = nSimVertices++;
00833   if ( nSimVertices/theVertexSize*theVertexSize == nSimVertices ) {
00834     theVertexSize *= 2;
00835     theSimVertices->resize(theVertexSize);
00836     theFSimVerticesType->resize(theVertexSize);
00837   }
00838 
00839   // Attach the end vertex to the particle (if accepted)
00840   if ( im !=-1 ) track(im).setEndVertex(vertexId);
00841 
00842   // Some transient information for FAMOS internal use
00843   (*theSimVertices)[vertexId] = FSimVertex(v,im,vertexId,this);
00844 
00845   (*theFSimVerticesType)[vertexId] = FSimVertexType(type);
00846 
00847   return vertexId;
00848 
00849 }
00850 
00851 void
00852 FBaseSimEvent::printMCTruth(const HepMC::GenEvent& myGenEvent) {
00853   
00854   std::cout << "Id  Gen Name       eta    phi     pT     E    Vtx1   " 
00855             << " x      y      z   " 
00856             << "Moth  Vtx2  eta   phi     R      Z   Da1  Da2 Ecal?" << std::endl;
00857 
00858   for ( HepMC::GenEvent::particle_const_iterator 
00859           piter  = myGenEvent.particles_begin();
00860           piter != myGenEvent.particles_end(); 
00861         ++piter ) {
00862     
00863     HepMC::GenParticle* p = *piter;
00864      /* */
00865     int partId = p->pdg_id();
00866     std::string name;
00867 
00868     if ( pdt->particle(ParticleID(partId)) !=0 ) {
00869       name = (pdt->particle(ParticleID(partId)))->name();
00870     } else {
00871       name = "none";
00872     }
00873   
00874     XYZTLorentzVector momentum1(p->momentum().px(),
00875                                 p->momentum().py(),
00876                                 p->momentum().pz(),
00877                                 p->momentum().e());
00878 
00879     int vertexId1 = 0;
00880 
00881     if ( !p->production_vertex() ) continue;
00882 
00883     XYZVector vertex1 (p->production_vertex()->position().x()/10.,
00884                        p->production_vertex()->position().y()/10.,
00885                        p->production_vertex()->position().z()/10.);
00886     vertexId1 = p->production_vertex()->barcode();
00887     
00888     std::cout.setf(std::ios::fixed, std::ios::floatfield);
00889     std::cout.setf(std::ios::right, std::ios::adjustfield);
00890     
00891     std::cout << std::setw(4) << p->barcode() << " " 
00892          << name;
00893     
00894     for(unsigned int k=0;k<11-name.length() && k<12; k++) std::cout << " ";  
00895     
00896     double eta = momentum1.eta();
00897     if ( eta > +10. ) eta = +10.;
00898     if ( eta < -10. ) eta = -10.;
00899     std::cout << std::setw(6) << std::setprecision(2) << eta << " " 
00900               << std::setw(6) << std::setprecision(2) << momentum1.phi() << " " 
00901               << std::setw(7) << std::setprecision(2) << momentum1.pt() << " " 
00902               << std::setw(7) << std::setprecision(2) << momentum1.e() << " " 
00903               << std::setw(4) << vertexId1 << " " 
00904               << std::setw(6) << std::setprecision(1) << vertex1.x() << " " 
00905               << std::setw(6) << std::setprecision(1) << vertex1.y() << " " 
00906               << std::setw(6) << std::setprecision(1) << vertex1.z() << " ";
00907 
00908     const HepMC::GenParticle* mother = 
00909       *(p->production_vertex()->particles_in_const_begin());
00910 
00911     if ( mother )
00912       std::cout << std::setw(4) << mother->barcode() << " ";
00913     else 
00914       std::cout << "     " ;
00915     
00916     if ( p->end_vertex() ) {  
00917       XYZTLorentzVector vertex2(p->end_vertex()->position().x()/10.,
00918                                 p->end_vertex()->position().y()/10.,
00919                                 p->end_vertex()->position().z()/10.,
00920                                 p->end_vertex()->position().t()/10.);
00921       int vertexId2 = p->end_vertex()->barcode();
00922 
00923       std::vector<const HepMC::GenParticle*> children;
00924       HepMC::GenVertex::particles_out_const_iterator firstDaughterIt = 
00925         p->end_vertex()->particles_out_const_begin();
00926       HepMC::GenVertex::particles_out_const_iterator lastDaughterIt = 
00927         p->end_vertex()->particles_out_const_end();
00928       for ( ; firstDaughterIt != lastDaughterIt ; ++firstDaughterIt ) {
00929         children.push_back(*firstDaughterIt);
00930       }      
00931 
00932       std::cout << std::setw(4) << vertexId2 << " "
00933                 << std::setw(6) << std::setprecision(2) << vertex2.eta() << " " 
00934                 << std::setw(6) << std::setprecision(2) << vertex2.phi() << " " 
00935                 << std::setw(5) << std::setprecision(1) << vertex2.pt() << " " 
00936                 << std::setw(6) << std::setprecision(1) << vertex2.z() << " ";
00937       for ( unsigned id=0; id<children.size(); ++id )
00938         std::cout << std::setw(4) << children[id]->barcode() << " ";
00939     }
00940     std::cout << std::endl;
00941 
00942   }
00943 
00944 }
00945 
00946 void
00947 FBaseSimEvent::print() const {
00948 
00949   std::cout << "  Id  Gen Name       eta    phi     pT     E    Vtx1   " 
00950             << " x      y      z   " 
00951             << "Moth  Vtx2  eta   phi     R      Z   Daughters Ecal?" << std::endl;
00952 
00953   for( int i=0; i<(int)nTracks(); i++ ) 
00954     std::cout << track(i) << std::endl;
00955 
00956   for( int i=0; i<(int)nVertices(); i++ ) 
00957     std::cout << "i = " << i << "  " << vertexType(i) << std::endl;
00958 
00959   
00960 
00961 }
00962 
00963 void 
00964 FBaseSimEvent::clear() {
00965 
00966   nSimTracks = 0;
00967   nSimVertices = 0;
00968   nGenParticles = 0;
00969   nChargedParticleTracks = 0;
00970 
00971 }
00972 
00973 void 
00974 FBaseSimEvent::addChargedTrack(int id) { 
00975   (*theChargedTracks)[nChargedParticleTracks++] = id;
00976   if ( nChargedParticleTracks/theChargedSize*theChargedSize 
00977        == nChargedParticleTracks ) {
00978     theChargedSize *= 2;
00979     theChargedTracks->resize(theChargedSize);
00980   }
00981 }
00982 
00983 int
00984 FBaseSimEvent::chargedTrack(int id) const {
00985   if (id>=0 && id<(int)nChargedParticleTracks) 
00986     return (*theChargedTracks)[id]; 
00987   else 
00988     return -1;
00989 }
00990 
00991 /* 
00992 const SimTrack & 
00993 FBaseSimEvent::embdTrack(int i) const {  
00994   return (*theSimTracks)[i].simTrack();
00995 }
00996 
00997 const SimVertex & 
00998 FBaseSimEvent::embdVertex(int i) const { 
00999   return (*theSimVertices)[i].simVertex();
01000 }
01001 */
01002 
01003 const HepMC::GenParticle* 
01004 FBaseSimEvent::embdGenpart(int i) const {
01005   return (*theGenParticles)[i]; 
01006 }
01007 
01008 /*
01009 FSimTrack&  
01010 FBaseSimEvent::track(int id) const { 
01011   return (*theSimTracks)[id];
01012 }
01013 
01014 
01015 FSimVertex&  
01016 FBaseSimEvent::vertex(int id) const { 
01017   return (*theSimVertices)[id];
01018 }
01019 */