00001
00002 #include "HepMC/GenEvent.h"
00003 #include "HepMC/GenVertex.h"
00004 #include "HepMC/GenParticle.h"
00005
00006
00007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00008
00009
00010 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
00011 #include "DataFormats/Candidate/interface/Candidate.h"
00012
00013
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
00026
00027 using namespace HepPDT;
00028
00029
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
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
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
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
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
00097 theBeamSpot = math::XYZPoint(0.0,0.0,0.0);
00098
00099
00100
00101
00102 lateVertexPosition = 2.5*2.5;
00103
00104
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
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
00125 myFilter = new KineParticleFilter(kine);
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136 }
00137
00138 FBaseSimEvent::~FBaseSimEvent(){
00139
00140
00141 theGenParticles->clear();
00142 theSimTracks->clear();
00143 theSimVertices->clear();
00144 theChargedTracks->clear();
00145 theFSimVerticesType->clear();
00146
00147
00148 delete theGenParticles;
00149 delete theSimTracks;
00150 delete theSimVertices;
00151 delete theChargedTracks;
00152 delete theFSimVerticesType;
00153 delete myFilter;
00154
00155
00156
00157
00158 }
00159
00160 void
00161 FBaseSimEvent::initializePdt(const HepPDT::ParticleDataTable* aPdt) {
00162
00163 pdt = aPdt;
00164
00165 }
00166
00167
00168
00169
00170
00171
00172
00173
00174 void
00175 FBaseSimEvent::fill(const HepMC::GenEvent& myGenEvent) {
00176
00177
00178 clear();
00179
00180
00181 addParticles(myGenEvent);
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191 }
00192
00193 void
00194 FBaseSimEvent::fill(const reco::GenParticleCollection& myGenParticles) {
00195
00196
00197 clear();
00198
00199
00200 addParticles(myGenParticles);
00201
00202 }
00203
00204 void
00205 FBaseSimEvent::fill(const std::vector<SimTrack>& simTracks,
00206 const std::vector<SimVertex>& simVertices) {
00207
00208
00209
00210
00211 clear();
00212
00213 unsigned nVtx = simVertices.size();
00214 unsigned nTks = simTracks.size();
00215
00216
00217 if ( nVtx == 0 ) return;
00218
00219
00220 std::vector<int> myVertices(nVtx,-1);
00221 std::vector<int> myTracks(nTks,-1);
00222
00223
00224
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
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
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
00253 addSimVertex(myFilter->vertex(), -1, FSimVertexType::PRIMARY_VERTEX);
00254 myVertices[0] = 0;
00255
00256 for( unsigned trackId=0; trackId<nTks; ++trackId ) {
00257
00258
00259 const SimTrack& track = simTracks[trackId];
00260
00261
00262
00263 int vertexId = track.vertIndex();
00264 const SimVertex& vertex = simVertices[vertexId];
00265
00266
00267
00268 int motherId = -1;
00269 if( !vertex.noParent() ) {
00270
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
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289 XYZTLorentzVector position(vertex.position().px(),vertex.position().py(),
00290 vertex.position().pz(),vertex.position().e());
00291 if ( myVertices[vertexId] == -1 )
00292
00293
00294
00295
00296 myVertices[vertexId] = addSimVertex(position,originId);
00297
00298
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
00308
00309
00310
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
00317
00318
00319
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
00337 for( unsigned vertexId=0; vertexId<nVtx; ++vertexId ) {
00338
00339
00340 if ( myVertices[vertexId] != -1 ) continue;
00341
00342
00343 const SimVertex& vertex = simVertices[vertexId];
00344
00345
00346 int motherId = -1;
00347 if( !vertex.noParent() ) {
00348
00349
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
00359
00360
00361
00362
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
00369 BaseParticlePropagator myPart;
00370 XYZTLorentzVector mom;
00371 XYZTLorentzVector pos;
00372
00373
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
00393 myPart = BaseParticlePropagator(RawParticle(mom,pos),0.,0.,4.);
00394 myPart.setCharge(myTrack.charge());
00395
00396
00397 myPart.propagateToPreshowerLayer1(false);
00398 if ( myTrack.notYetToEndVertex(myPart.vertex()) && myPart.getSuccess()>0 )
00399 myTrack.setLayer1(myPart,myPart.getSuccess());
00400
00401
00402 myPart.propagateToPreshowerLayer2(false);
00403 if ( myTrack.notYetToEndVertex(myPart.vertex()) && myPart.getSuccess()>0 )
00404 myTrack.setLayer2(myPart,myPart.getSuccess());
00405
00406
00407 myPart.propagateToEcalEntrance(false);
00408 if ( myTrack.notYetToEndVertex(myPart.vertex()) )
00409 myTrack.setEcal(myPart,myPart.getSuccess());
00410
00411
00412 myPart.propagateToHcalEntrance(false);
00413 if ( myTrack.notYetToEndVertex(myPart.vertex()) )
00414 myTrack.setHcal(myPart,myPart.getSuccess());
00415
00416
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
00435 if ( myGenEvent.particles_empty() ) return;
00436
00437
00438 int offset = nGenParts();
00439
00440
00441 HepMC::GenVertex* primaryVertex = *(myGenEvent.vertices_begin());
00442
00443
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
00450
00451 XYZTLorentzVector primaryVertexPosition(primaryVertex->position().x()/10.,
00452 primaryVertex->position().y()/10.,
00453 primaryVertex->position().z()/10.,
00454 primaryVertex->position().t()/10.);
00455
00456 primaryVertexPosition *= (1-primaryMother);
00457
00458
00459
00460
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
00472 myFilter->setMainVertex(primaryVertexPosition+smearedVertex);
00473
00474
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
00484 for ( piter = pbegin; piter != pend; ++piter ) {
00485
00486
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
00499
00500
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
00521
00522 bool testStable = p->status()%1000==1;
00523
00524
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
00533 if ( decayPosition.Perp2() > lateVertexPosition ) testStable = true;
00534 }
00535 }
00536
00537
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
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
00558 if ( endVertexPosition.Perp2() < lateVertexPosition ) {
00559 break;
00560 }
00561 }
00562 }
00563 testDaugh=true;
00564 break;
00565 }
00566 }
00567 }
00568
00569
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
00582 if ( testStable || testDaugh || testDecay ) {
00583
00584
00585
00586
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
00606
00607 int theTrack = testStable && p->end_vertex() ?
00608
00609 addSimTrack(&part,originVertex, nGenParts()-offset,p->end_vertex()) :
00610
00611 addSimTrack(&part,originVertex, nGenParts()-offset);
00612
00613 if (
00614
00615 !p->end_vertex() ||
00616
00617
00618 ( testStable && p->end_vertex() && !p->end_vertex()->particles_out_size() )
00619
00620 ) continue;
00621
00622
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
00630 int theVertex = addSimVertex(decayVertex,theTrack, FSimVertexType::DECAY_VERTEX);
00631
00632 if ( theVertex != -1 ) myGenVertices[p->barcode()-initialBarcode] = theVertex;
00633
00634
00635 }
00636 }
00637
00638 }
00639
00640 void
00641 FBaseSimEvent::addParticles(const reco::GenParticleCollection& myGenParticles) {
00642
00643
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
00653 int offset = nTracks();
00654
00655
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
00666 XYZTLorentzVector primaryVertex (myGenParticles[ip].vx(),
00667 myGenParticles[ip].vy(),
00668 myGenParticles[ip].vz(),
00669 0.);
00670
00671
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
00683 myFilter->setMainVertex(primaryVertex+smearedVertex);
00684
00685
00686 int mainVertex = addSimVertex(myFilter->vertex(), -1, FSimVertexType::PRIMARY_VERTEX);
00687
00688
00689 for ( ; ip<nParticles; ++ip ) {
00690
00691
00692
00693 nGenParticles++;
00694 const reco::GenParticle& p = myGenParticles[ip];
00695
00696
00697
00698
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
00709
00710 bool testStable = p.status()%1000==1;
00711
00712
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
00720 if ( decayPosition.Perp2() > lateVertexPosition ) testStable = true;
00721 }
00722 }
00723
00724
00725 bool testDaugh = false;
00726 unsigned int nDaughters = p.numberOfDaughters();
00727 if ( !testStable &&
00728
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
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
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
00762 int theTrack = addSimTrack(&part,originVertex, nGenParts()-offset);
00763
00764
00765 if ( !nDaughters ) continue;
00766 const reco::Candidate* daughter = p.daughter(0);
00767
00768
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
00777 }
00778 }
00779
00780
00781
00782
00783 }
00784
00785 int
00786 FBaseSimEvent::addSimTrack(const RawParticle* p, int iv, int ig,
00787 const HepMC::GenVertex* ev) {
00788
00789
00790
00791 if ( !myFilter->accept(p) && ig >= -1 ) return -1;
00792
00793
00794 int trackId = nSimTracks++;
00795 if ( nSimTracks/theTrackSize*theTrackSize == nSimTracks ) {
00796 theTrackSize *= 2;
00797 theSimTracks->resize(theTrackSize);
00798 }
00799
00800
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
00812 (*theSimTracks)[trackId] = ev ?
00813
00814 FSimTrack(p,iv,ig,trackId,this,
00815 ev->position().t()/10.
00816 * p->PDGmass()
00817 / std::sqrt(p->momentum().Vect().Mag2())) :
00818
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
00829 if ( !myFilter->accept(v) ) return -1;
00830
00831
00832 int vertexId = nSimVertices++;
00833 if ( nSimVertices/theVertexSize*theVertexSize == nSimVertices ) {
00834 theVertexSize *= 2;
00835 theSimVertices->resize(theVertexSize);
00836 theFSimVerticesType->resize(theVertexSize);
00837 }
00838
00839
00840 if ( im !=-1 ) track(im).setEndVertex(vertexId);
00841
00842
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
00993
00994
00995
00996
00997
00998
00999
01000
01001
01002
01003 const HepMC::GenParticle*
01004 FBaseSimEvent::embdGenpart(int i) const {
01005 return (*theGenParticles)[i];
01006 }
01007
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019