TruthSvc.cxx 22.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
/*
  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/

///////////////////////////////////////////////////////////////////
// TruthSvc.cxx, (c) ATLAS Detector software
///////////////////////////////////////////////////////////////////

// class header
#include "TruthSvc.h"
// other ISF_HepMC includes
#include "ISF_HepMC_Interfaces/ITruthStrategy.h"
// ISF includes
#include "ISF_Event/ITruthIncident.h"
// Framework
#include "GaudiKernel/ISvcLocator.h"
#include "StoreGate/StoreGateSvc.h"
#include "GaudiKernel/SystemOfUnits.h"
// Barcode
#include "BarcodeInterfaces/IBarcodeSvc.h"
21
22
//
#include "TruthUtils/HepMCHelpers.h" // for MC::findChildren(...)
23
24
25
26
27
28
29
30
31
32
33
// HepMC includes
#include "HepMC/SimpleVector.h"
#include "HepMC/GenParticle.h"
#include "HepMC/GenEvent.h"
#include "HepMC/GenVertex.h"
// CLHEP includes
#include "CLHEP/Geometry/Point3D.h"

// DetectorDescription
#include "AtlasDetDescr/AtlasRegionHelper.h"

34
35
36
37
#include <sstream>

#undef DEBUG_TRUTHSVC

38
39
/** Constructor **/
ISF::TruthSvc::TruthSvc(const std::string& name,ISvcLocator* svc) :
40
  base_class(name,svc),
41
42
43
44
45
46
47
48
49
  m_barcodeSvc("BarcodeSvc",name),
  m_geoStrategies(),
  m_numStrategies(),
  m_skipIfNoChildren(true),
  m_skipIfNoParentBarcode(true),
  m_ignoreUndefinedBarcodes(false),
  m_passWholeVertex(true),
  m_forceEndVtxRegionsVec(),
  m_forceEndVtx(),
John Chapman's avatar
John Chapman committed
50
  m_quasiStableParticlesIncluded(false)
51
{
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
  // the barcode service (used to compute Vertex Barco des)
  declareProperty("BarcodeSvc",                        m_barcodeSvc              );
  // MCTruth writing strategies
  declareProperty("SkipIfNoChildren",                  m_skipIfNoChildren        );
  declareProperty("SkipIfNoParentBarcode",             m_skipIfNoParentBarcode   );
  declareProperty("IgnoreUndefinedBarcodes",           m_ignoreUndefinedBarcodes );
  declareProperty("PassWholeVertices",                 m_passWholeVertex         );
  // the truth strategies for the different AtlasDetDescr regions
  declareProperty("BeamPipeTruthStrategies",           m_geoStrategyHandles[AtlasDetDescr::fAtlasForward] );
  declareProperty("IDTruthStrategies",                 m_geoStrategyHandles[AtlasDetDescr::fAtlasID]      );
  declareProperty("CaloTruthStrategies",               m_geoStrategyHandles[AtlasDetDescr::fAtlasCalo]    );
  declareProperty("MSTruthStrategies",                 m_geoStrategyHandles[AtlasDetDescr::fAtlasMS]      );
  declareProperty("CavernTruthStrategies",             m_geoStrategyHandles[AtlasDetDescr::fAtlasCavern]  );
  // attach end-vertex if parent particle dies for the different AtlasDetDescr regions
  declareProperty("ForceEndVtxInRegions",              m_forceEndVtxRegionsVec );

  declareProperty("QuasiStableParticlesIncluded",      m_quasiStableParticlesIncluded);
69
70
71
72
73
74
75
76
77
}

ISF::TruthSvc::~TruthSvc()
{}


/** framework methods */
StatusCode ISF::TruthSvc::initialize()
{
78
79
80
81
  ATH_MSG_VERBOSE( "initialize()" );

  // Screen output
  ATH_MSG_DEBUG("--------------------------------------------------------");
82

83
84
85
86
87
  // retrieve BarcodeSvc
  if ( m_barcodeSvc.retrieve().isFailure() ) {
    ATH_MSG_FATAL("Could not retrieve BarcodeService. Abort.");
    return StatusCode::FAILURE;
  }
88

89
90
91
92
93
94
  // retrieve the strategies for each atlas region (Athena Alg Tools)
  // and setup whether we want to write end-vertices in this region whenever a truth particle dies
  for ( unsigned short geoID=AtlasDetDescr::fFirstAtlasRegion; geoID<AtlasDetDescr::fNumAtlasRegions; ++geoID) {
    if ( m_geoStrategyHandles[geoID].retrieve().isFailure() ) {
      ATH_MSG_FATAL("Could not retrieve TruthStrategy Tool Array for SimGeoID="
                    << AtlasDetDescr::AtlasRegionHelper::getName(geoID) << ". Abort.");
95
96
97
      return StatusCode::FAILURE;
    }

98
99
100
101
102
103
104
    // copy a pointer to the strategy instance to the local
    // array of pointers (for faster access)
    unsigned short curNumStrategies = m_geoStrategyHandles[geoID].size();
    m_numStrategies[geoID] = curNumStrategies;
    m_geoStrategies[geoID] = new ISF::ITruthStrategy*[curNumStrategies];
    for ( unsigned short i = 0; i < curNumStrategies; ++i) {
      m_geoStrategies[geoID][i] = &(*m_geoStrategyHandles[geoID][i]);
105
106
    }

107
108
109
110
111
112
113
114
115
    // create an end-vertex for all truth particles ending in the current AtlasRegion?
    bool forceEndVtx = std::find( m_forceEndVtxRegionsVec.begin(),
                                  m_forceEndVtxRegionsVec.end(),
                                  geoID ) != m_forceEndVtxRegionsVec.end();
    m_forceEndVtx[geoID] = forceEndVtx;
  }

  ATH_MSG_VERBOSE("initialize() successful");
  return StatusCode::SUCCESS;
116
117
118
119
120
121
}


/** framework methods */
StatusCode ISF::TruthSvc::finalize()
{
122
123
  ATH_MSG_VERBOSE("Finalizing ...");
  return StatusCode::SUCCESS;
124
125
126
127
128
129
130
131
132
}


/** Initialize the TruthSvc and the truthSvc */
StatusCode ISF::TruthSvc::initializeTruthCollection()
{
  return StatusCode::SUCCESS;
}

133
/** Delete child vertex */
John Chapman's avatar
John Chapman committed
134
135
136
137
138
139
140
141
142
143
144
void ISF::TruthSvc::deleteChildVertex(HepMC::GenVertex* vertex) const {
  std::vector<HepMC::GenVertex*> verticesToDelete;
  verticesToDelete.resize(0);
  verticesToDelete.push_back(vertex);
  for ( unsigned short i = 0; i<verticesToDelete.size(); ++i ) {
    HepMC::GenVertex* vtx = verticesToDelete.at(i);
    for (HepMC::GenVertex::particles_out_const_iterator iter = vtx->particles_out_const_begin();
         iter != vtx->particles_out_const_end(); ++iter) {
      if( (*iter) && (*iter)->end_vertex() ) {
        verticesToDelete.push_back( (*iter)->end_vertex() );
      }
145
    }
John Chapman's avatar
John Chapman committed
146
    vtx->parent_event()->remove_vertex(vtx);
147
148
149
150
  }
  return;
}

151
152
153
154
155
156
157

StatusCode ISF::TruthSvc::releaseEvent() {
  return StatusCode::SUCCESS;
}


/** Register a truth incident */
158
void ISF::TruthSvc::registerTruthIncident( ISF::ITruthIncident& ti, bool saveAllChildren) const {
159

160
  const bool passWholeVertex = m_passWholeVertex || saveAllChildren;
161
  // pass whole vertex or individual child particles
162
  ti.setPassWholeVertices(passWholeVertex);
163
164
165
166
167
168

  // the GeoID
  AtlasDetDescr::AtlasRegion geoID = ti.geoID();

  // check geoID assigned to the TruthIncident
  if ( !validAtlasRegion(geoID) ) {
169
170
171
    const auto& position = ti.position();
    ATH_MSG_ERROR("Unable to register truth incident with unknown SimGeoID="<< geoID
                  << " at position z=" << position.z() << " r=" << position.perp());
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
    return;
  }

  ATH_MSG_VERBOSE( "Registering TruthIncident for SimGeoID="
                   << AtlasDetDescr::AtlasRegionHelper::getName(geoID) );

  // number of child particles
  unsigned short numSec = ti.numberOfChildren();
  if ( m_skipIfNoChildren && (numSec==0) ) {
    ATH_MSG_VERBOSE( "No child particles present in the TruthIncident,"
                     << " will not record this TruthIncident.");
    return;
  }

  // the parent particle -> get its barcode
  Barcode::ParticleBarcode parentBC = ti.parentBarcode();
  if ( m_skipIfNoParentBarcode && (parentBC==Barcode::fUndefinedBarcode) ) {
    ATH_MSG_VERBOSE( "Parent particle in TruthIncident does not have a barcode,"
                     << " will not record this TruthIncident.");
    return;
  }

  // loop over registered truth strategies for given geoID
  bool pass = false;
  for ( unsigned short stratID=0; (!pass) && (stratID<m_numStrategies[geoID]); stratID++) {
    // (*) test if given TruthIncident passes current strategy
    pass = m_geoStrategies[geoID][stratID]->pass(ti);
  }

  if (pass) {
202
    ATH_MSG_VERBOSE("At least one TruthStrategy passed.");
203
    // at least one truth strategy returned true
204
    //  -> record incident
205
    recordIncidentToMCTruth(ti, passWholeVertex);
206
207
208
209
210
211

  } else {
    // none of the truth strategies returned true
    //  -> child particles will NOT be added to the TruthEvent collection
    // attach parent particle end vertex if it gets killed by this interaction
    if ( m_forceEndVtx[geoID] && !ti.parentSurvivesIncident() ) {
212
213
214
215
216
217
      ATH_MSG_VERBOSE("No TruthStrategies passed and parent destroyed - create end vertex.");
      const bool replaceVertex(true);
      this->createGenVertexFromTruthIncident( ti, replaceVertex );

#ifdef DEBUG_TRUTHSVC
      const std::string survival = (ti.parentSurvivesIncident()) ? "parent survives" : "parent destroyed";
218
219
220
221
222
      const std::string vtxType = (ti.interactionClassification()==ISF::STD_VTX) ? "Normal" : "Quasi-stable";
      ATH_MSG_INFO("TruthSvc: " << vtxType << " vertex + " << survival
                   << ", TI Class: " << ti.interactionClassification()
                   << ", ProcessType: " << ti.physicsProcessCategory()
                   << ", ProcessSubType: " << ti.physicsProcessCode());
223
224
#endif

225
226
227
228
229
230
231
232
233
234
    }

    //  -> assign shared barcode to all child particles (if barcode service supports it)
    setSharedChildParticleBarcode( ti);
  }

  return;
}

/** Record the given truth incident to the MC Truth */
235
void ISF::TruthSvc::recordIncidentToMCTruth( ISF::ITruthIncident& ti, bool passWholeVertex) const {
236
237
238
#ifdef  DEBUG_TRUTHSVC
  ATH_MSG_INFO("Starting recordIncidentToMCTruth(...)");
#endif
239
240
241
242
  Barcode::PhysicsProcessCode processCode = ti.physicsProcessCode();
  Barcode::ParticleBarcode       parentBC = ti.parentBarcode();

  // record the GenVertex
243
244
245
246
  const bool replaceVertex(false);
  HepMC::GenVertex *vtx = createGenVertexFromTruthIncident(ti, replaceVertex);
  const ISF::InteractionClass_t classification = ti.interactionClassification();
#ifdef DEBUG_TRUTHSVC
247
248
249
250
251
252
  const std::string survival = (ti.parentSurvivesIncident()) ? "parent survives" : "parent destroyed";
  const std::string vtxType = (ti.interactionClassification()==ISF::STD_VTX) ? "Normal" : "Quasi-stable";
  ATH_MSG_INFO("TruthSvc: " << vtxType << " vertex + " << survival
               << ", TI Class: " << ti.interactionClassification()
               << ", ProcessType: " << ti.physicsProcessCategory()
               << ", ProcessSubType: " << ti.physicsProcessCode());
253
#endif
254
255
256

  ATH_MSG_VERBOSE ( "Outgoing particles:" );
  // update parent barcode and add it to the vertex as outgoing particle
257
258
  Barcode::ParticleBarcode newPrimBC = Barcode::fUndefinedBarcode;
  if (classification == ISF::QS_SURV_VTX) {
259
260
261
262
    // Special case when a particle with a pre-defined decay interacts
    // and survives.
    // Set the barcode to the next available value below the simulation
    // barcode offset.
263
264
265
266
267
    newPrimBC = this->maxGeneratedParticleBarcode(ti.parentParticle()->parent_event())+1;
  }
  else {
    newPrimBC = m_barcodeSvc->incrementBarcode( parentBC, processCode);
  }
268
269
270
271
  if ( newPrimBC == Barcode::fUndefinedBarcode) {
    if (m_ignoreUndefinedBarcodes) {
      ATH_MSG_WARNING("Unable to generate new Particle Barcode. Continuing due to 'IgnoreUndefinedBarcodes'==True");
    } else {
272
      ATH_MSG_FATAL("Unable to generate new Particle Barcode. Aborting");
273
274
275
276
      abort();
    }
  }

277
278
  HepMC::GenParticle *parentBeforeIncident = ti.parentParticle();
  HepMC::GenParticle *parentAfterIncident = ti.parentParticleAfterIncident( newPrimBC ); // This call changes ti.parentParticle() output
279
280
  if(parentAfterIncident) {
    ATH_MSG_VERBOSE ( "Parent After Incident: " << *parentAfterIncident);
281
282
283
284
285
286
287
288
    if (classification==ISF::QS_SURV_VTX) {
      // Special case when a particle with a pre-defined decay
      // interacts and survives.
      // 1) As the parentParticleAfterIncident has a pre-defined decay
      // its status should be to 2.
      parentAfterIncident->set_status(2);
      // 2) A new GenVertex for the intermediate interaction should be
      // added.
289
290
291
292
293
      std::unique_ptr<HepMC::GenVertex> newVtx = std::make_unique<HepMC::GenVertex>( vtx->position(), vtx->id(), vtx->weights() );
#ifdef DEBUG_TRUTHSVC
      ATH_MSG_INFO("New GenVertex 1: " << *(newVtx.get()) );
      ATH_MSG_INFO("New QS GenVertex 1: " << *(newVtx.get()) );
#endif
294
      HepMC::GenEvent *mcEvent = parentBeforeIncident->parent_event();
295
296
297
298
299
300
301
302
303
304
305
306
      newVtx->suggest_barcode( this->maxGeneratedVertexBarcode(mcEvent)-1 );
#ifdef DEBUG_TRUTHSVC
      ATH_MSG_INFO("New QSGenVertex 2: " << *(newVtx.get()) );
#endif
      auto tmpVtx = newVtx.get();
#ifdef DEBUG_TRUTHSVC
      ATH_MSG_INFO("New QS GenVertex 3: " << (*tmpVtx) );
#endif
      if(!mcEvent->add_vertex( newVtx.release() )) {
        ATH_MSG_FATAL("Failed to add GenVertex to GenEvent.");
        abort();
      }
307
      tmpVtx->add_particle_in( parentBeforeIncident );
308
309
310
311
312
313
314
      tmpVtx->add_particle_out( parentAfterIncident );
      vtx->add_particle_in( parentAfterIncident );
      vtx = tmpVtx;
    }
    else {
      vtx->add_particle_out( parentAfterIncident );
    }
315
316
  }

317
  const bool isQuasiStableVertex = (classification == ISF::QS_PREDEF_VTX); // QS_DEST_VTX and QS_SURV_VTX should be treated as normal from now on.
318
319
  // add child particles to the vertex
  unsigned short numSec = ti.numberOfChildren();
320
  if (isQuasiStableVertex) {
321
322
323
324
325
326
327
328
329
330
331
332
333
334
    // Here we are checking if the existing GenVertex has the same
    // number of child particles as the truth incident.
    // FIXME should probably make this part a separate function and
    // also check if the pdgids of the child particles are the same
    // too.
    unsigned short nVertexChildren=vtx->particles_out_size();
    if(parentAfterIncident) { nVertexChildren-=1; }
    if(nVertexChildren!=numSec) {
      ATH_MSG_WARNING("Existing vertex has " << nVertexChildren << " children. " <<
                      "Number of secondaries in current truth incident = " << numSec);
    }
    ATH_MSG_VERBOSE("Existing vertex has " << nVertexChildren << " children. " <<
                 "Number of secondaries in current truth incident = " << numSec);
  }
335
  const std::vector<HepMC::GenParticle*> childParticleVector = (isQuasiStableVertex) ? MC::findChildren(ti.parentParticle()) : std::vector<HepMC::GenParticle*>();
336
  std::vector<HepMC::GenParticle*> matchedChildParticles;
337
338
  for ( unsigned short i=0; i<numSec; ++i) {

339
    bool writeOutChild = isQuasiStableVertex || passWholeVertex || ti.childPassedFilters(i);
340
341

    if (writeOutChild) {
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
      HepMC::GenParticle *p = nullptr;
      if(isQuasiStableVertex) {
        //Find matching GenParticle in GenVertex
        const int childPDGcode= ti.childPdgCode(i);
        bool noMatch(true);
        for(auto childParticle : childParticleVector) {
          if( (childParticle->pdg_id() == childPDGcode) && std::count(matchedChildParticles.begin(),matchedChildParticles.end(),childParticle)==0) {
            noMatch=false;
            ATH_MSG_VERBOSE("Found a matching Quasi-stable GenParticle with PDGcode " << childPDGcode << ":\n\t" << *childParticle );
            matchedChildParticles.push_back(childParticle);
            // FIXME There is a weakness in the code here for
            // vertices with multiple children with the same
            // pdgid. The code relies on the order of the children in
            // childParticleVector being the same as in the
            // truthIncident...
            p = ti.updateChildParticle( i, childParticle );
            break;
          }
        }
        if (noMatch) {
          std::ostringstream warnStream;
          warnStream << "Failed to find a Quasi-stable GenParticle with PDGID " << childPDGcode << ". Options are: \n";
          for(auto childParticle : childParticleVector) {
            warnStream << childParticle->pdg_id() <<"\n";
          }
          ATH_MSG_WARNING(warnStream.str());
368
369
        }
      }
370
371
372
373
374
375
376
377
378
379
      else {
        // generate a new barcode for the child particle
        Barcode::ParticleBarcode secBC = m_barcodeSvc->newSecondary( parentBC, processCode);
        if ( secBC == Barcode::fUndefinedBarcode) {
          if (m_ignoreUndefinedBarcodes)
            ATH_MSG_WARNING("Unable to generate new Secondary Particle Barcode. Continuing due to 'IgnoreUndefinedBarcodes'==True");
          else {
            ATH_MSG_ERROR("Unable to generate new Secondary Particle Barcode. Aborting");
            abort();
          }
380
        }
381
382
383
        p = ti.childParticle( i, secBC );
        // add particle to vertex
        vtx->add_particle_out( p);
384
      }
385
      ATH_MSG_VERBOSE ( "Writing out " << i << "th child particle: " << *p);
386
387
388
389
390
391
    } // <-- if write out child particle
    else {
      ATH_MSG_VERBOSE ( "Not writing out " << i << "th child particle." );
    }

  } // <-- loop over all child particles
392
  ATH_MSG_VERBOSE("--------------------------------------------------------");
393
394
395
}

/** Record the given truth incident to the MC Truth */
396
HepMC::GenVertex *ISF::TruthSvc::createGenVertexFromTruthIncident( ISF::ITruthIncident& ti,
John Chapman's avatar
John Chapman committed
397
                                                                   bool replaceExistingGenVertex) const {
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431

  Barcode::PhysicsProcessCode processCode = ti.physicsProcessCode();
  Barcode::ParticleBarcode       parentBC = ti.parentBarcode();

  std::vector<double> weights(1);
  Barcode::ParticleBarcode primaryBC = parentBC % m_barcodeSvc->particleGenerationIncrement();
  weights[0] = static_cast<double>( primaryBC );

  // Check for a previous end vertex on this particle.  If one existed, then we should put down next to this
  //  a new copy of the particle.  This is the agreed upon version of the quasi-stable particle truth, where
  //  the vertex at which we start Q-S simulation no longer conserves energy, but we keep both copies of the
  //  truth particles
  HepMC::GenParticle *parent = ti.parentParticle();
  if (!parent) {
    ATH_MSG_ERROR("Unable to write particle interaction to MC truth due to missing parent HepMC::GenParticle instance");
    abort();
  }
  HepMC::GenEvent *mcEvent = parent->parent_event();
  if (!mcEvent) {
    ATH_MSG_ERROR("Unable to write particle interaction to MC truth due to missing parent HepMC::GenEvent instance");
    abort();
  }

  // generate vertex
  Barcode::VertexBarcode vtxbcode = m_barcodeSvc->newVertex( parentBC, processCode );
  if ( vtxbcode == Barcode::fUndefinedBarcode) {
    if (m_ignoreUndefinedBarcodes) {
      ATH_MSG_WARNING("Unable to generate new Truth Vertex Barcode. Continuing due to 'IgnoreUndefinedBarcodes'==True");
    } else {
      ATH_MSG_ERROR("Unable to generate new Truth Vertex Barcode. Aborting");
      abort();
    }
  }
  int vtxID = 1000 + static_cast<int>(processCode);
432
  std::unique_ptr<HepMC::GenVertex> vtx = std::make_unique<HepMC::GenVertex>( ti.position(), vtxID, weights );
433
434
435
436
437
438
439
440
441
  vtx->suggest_barcode( vtxbcode );

  if (parent->end_vertex()){
    if(!m_quasiStableParticlesIncluded) {
      ATH_MSG_WARNING("Parent particle found with an end vertex attached.  This should only happen");
      ATH_MSG_WARNING("in the case of simulating quasi-stable particles.  That functionality");
      ATH_MSG_WARNING("is not yet validated in ISF, so you'd better know what you're doing.");
      ATH_MSG_WARNING("Will delete the old vertex and swap in the new one.");
    }
442
443
444
445
446
447
448
449
450
451
    auto* oldVertex = parent->end_vertex();
#ifdef DEBUG_TRUTHSVC
    ATH_MSG_VERBOSE("createGVfromTI Existing QS GenVertex 1: " << *oldVertex );
    ATH_MSG_VERBOSE("createGVfromTI QS Parent 1: " << *parent);
#endif
    if(replaceExistingGenVertex) {
      vtx->add_particle_in( parent );
      ATH_MSG_VERBOSE("createGVfromTI Replacement QS GenVertex: " << *(vtx.get()) );
      mcEvent->add_vertex( vtx.release() );
      // Delete oldVertex and children here
John Chapman's avatar
John Chapman committed
452
      this->deleteChildVertex(oldVertex);
453
    }
454
455
    else {
      //oldVertex->suggest_barcode( vtxbcode );
456
457
458
459
460
461
462
463
464
465
466
467
      const auto& old_pos=oldVertex->position();
      const auto& new_pos=ti.position();
      HepMC::ThreeVector diff(new_pos.x()-old_pos.x(),new_pos.y()-old_pos.y(),new_pos.z()-old_pos.z()); //complicated, but HepMC::ThreeVector and FourVector have no + or - operators
      
      if(diff.r()>1*Gaudi::Units::mm) { //Check for a change of the vertex position by more than 1mm
        ATH_MSG_WARNING("For particle: " << *parent);
        ATH_MSG_WARNING("  decay vertex before QS partice sim: " << *oldVertex );
        oldVertex->set_position( ti.position() );
        ATH_MSG_WARNING("  decay vertex after  QS partice sim:  " << *oldVertex );
      } else {
        oldVertex->set_position( ti.position() );
      }  
468
469
      oldVertex->set_id( vtxID );
      oldVertex->weights() = weights;
470
471
472
#ifdef DEBUG_TRUTHSVC
      ATH_MSG_VERBOSE("createGVfromTI Existing QS GenVertex 2: " << *oldVertex );
#endif
473
474
475
476
477
    }
#ifdef DEBUG_TRUTHSVC
    ATH_MSG_VERBOSE ( "createGVfromTI QS End Vertex representing process: " << processCode << ", for parent with barcode "<<parentBC<<". Creating." );
    ATH_MSG_VERBOSE ( "createGVfromTI QS Parent 2: " << *parent);
#endif
478
  } else { // Normal simulation
479
480
481
#ifdef DEBUG_TRUTHSVC
    ATH_MSG_VERBOSE ("createGVfromTI Parent 1: " << *parent);
#endif
482
483
    // add parent particle to vtx
    vtx->add_particle_in( parent );
484
485
486
487
488
#ifdef DEBUG_TRUTHSVC
    ATH_MSG_VERBOSE ( "createGVfromTI End Vertex representing process: " << processCode << ", for parent with barcode "<<parentBC<<". Creating." );
    ATH_MSG_VERBOSE ( "createGVfromTI Parent 2: " << *parent);
#endif
    mcEvent->add_vertex( vtx.release() );
489
490
  }

491
  return parent->end_vertex();
492
493
494
}

/** Set shared barcode for child particles particles */
John Chapman's avatar
John Chapman committed
495
void ISF::TruthSvc::setSharedChildParticleBarcode( ISF::ITruthIncident& ti) const {
496
497
498
499
500
501
502
503
504
505
506
507
508
509
  Barcode::PhysicsProcessCode processCode = ti.physicsProcessCode();
  Barcode::ParticleBarcode       parentBC = ti.parentBarcode();

  ATH_MSG_VERBOSE ( "End Vertex representing process: " << processCode << ". TruthIncident failed cuts. Skipping.");

  // generate one new barcode for all child particles
  Barcode::ParticleBarcode childBC = m_barcodeSvc->sharedChildBarcode( parentBC, processCode);

  // propagate this barcode into the TruthIncident only if
  // it is a proper barcode, ie !=fUndefinedBarcode
  if (childBC != Barcode::fUndefinedBarcode) {
    ti.setAllChildrenBarcodes( childBC );
  }
}
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535

int ISF::TruthSvc::maxGeneratedParticleBarcode(HepMC::GenEvent *genEvent) const {
  int maxBarcode=0;
  const int firstSecondaryParticleBarcode(m_barcodeSvc->secondaryParticleBcOffset());
  HepMC::GenEvent::particle_const_iterator currentGenParticleIter;
  for (currentGenParticleIter= genEvent->particles_begin();
       currentGenParticleIter!= genEvent->particles_end();
       ++currentGenParticleIter) {
    const int barcode((*currentGenParticleIter)->barcode());
    if(barcode > maxBarcode && barcode < firstSecondaryParticleBarcode) { maxBarcode=barcode; }
  }
  return maxBarcode;
}

int ISF::TruthSvc::maxGeneratedVertexBarcode(HepMC::GenEvent *genEvent) const {
  int maxBarcode=0;
  const int firstSecondaryVertexBarcode(m_barcodeSvc->secondaryVertexBcOffset());
  HepMC::GenEvent::vertex_const_iterator currentGenVertexIter;
  for (currentGenVertexIter= genEvent->vertices_begin();
       currentGenVertexIter!= genEvent->vertices_end();
       ++currentGenVertexIter) {
    const int barcode((*currentGenVertexIter)->barcode());
    if(barcode < maxBarcode && barcode > firstSecondaryVertexBarcode) { maxBarcode=barcode; }
  }
  return maxBarcode;
}