FastCaloSimSvc.cxx 16.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
/*
  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/

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

// class header include
10
#include "FastCaloSimSvc.h"
11
12
13
14
15
16
17

// StoreGate
#include "StoreGate/StoreGateSvc.h"
#include "ISF_Interfaces/IParticleBroker.h"

// ISF includes
#include "ISF_Event/ISFParticle.h"
18
#include "ISF_Event/ISFParticleVector.h"
19
20
21
22
23
24
25
26
27
28
29
30

// HepMC include needed for FastCaloSim
#include "HepMC/GenParticle.h"
#include "HepMC/GenVertex.h"
#include "HepMC/SimpleVector.h"
#include "CLHEP/Units/SystemOfUnits.h"
// McEventCollection
#include "GeneratorObjects/McEventCollection.h"
// Calo Cell includes
#include "CaloEvent/CaloCell.h"
#include "CaloEvent/CaloCellContainer.h"
#include "NavFourMom/INavigable4MomentumCollection.h"
31
32
// Barcode
#include "BarcodeEvent/Barcode.h"
33
34
35
36
37
38
39
40
41
// use FastShowerCellBuilderTool for actual simulation
//#include "FastSimulationEvent/GenParticleEnergyDepositMap.h"
#include "FastCaloSim/FastShowerCellBuilderTool.h"
// PunchThrough Tool
#include "ISF_FastCaloSimInterfaces/IPunchThroughTool.h"

/** Constructor **/
ISF::FastCaloSimSvc::FastCaloSimSvc(const std::string& name,ISvcLocator* svc) :
  BaseSimulationSvc(name, svc),
42
  m_extrapolator(),
43
44
45
46
47
48
49
50
51
52
53
54
55
  m_ownPolicy(static_cast<int>(SG::VIEW_ELEMENTS)),
  m_batchProcessMcTruth(false),
  m_simulateUndefinedBCs(false),
  m_caloCellsOutputName("AllCalo"),
  m_caloCellHack(false),
  m_doPunchThrough(false),
  m_caloCellMakerTools_setup(),
  m_caloCellMakerTools_simulate(),
  m_caloCellMakerTools_release(),
  m_punchThroughTool(""),
  m_theContainer(0),
  m_particleBroker ("ISF_ParticleBroker",name)
{
56
  // where to go
57
58
59
60
61
62
63
64
  declareProperty("OwnPolicy",                         m_ownPolicy) ;
  declareProperty("CaloCellMakerTools_setup"   ,       m_caloCellMakerTools_setup) ;
  declareProperty("CaloCellMakerTools_simulate",       m_caloCellMakerTools_simulate) ;
  declareProperty("CaloCellMakerTools_release" ,       m_caloCellMakerTools_release) ;
  declareProperty("PunchThroughTool",                  m_punchThroughTool);
  declareProperty("CaloCellsOutputName",               m_caloCellsOutputName) ;
  declareProperty("CaloCellHack",                      m_caloCellHack) ;
  declareProperty("DoPunchThroughSimulation", 	       m_doPunchThrough) ;
65
  declareProperty("Extrapolator",                      m_extrapolator );
66
67
68
69
70
71
72
73
74
75
76
  declareProperty("SimulateUndefinedBarcodeParticles",
                  m_simulateUndefinedBCs,
                  "Whether or not to simulate paritcles with undefined barcode" );
  declareProperty("ParticleBroker",
                  m_particleBroker,
                  "ISF ParticleBroker Svc" );
  declareProperty("BatchProcessMcTruth",
                  m_batchProcessMcTruth=false,
                  "Run the FastShowerCellBuilders on the McTruth at the end of the event" );
}

77
ISF::FastCaloSimSvc::~FastCaloSimSvc()
78
79
80
81
82
83
84
85
{}

/** framework methods */
StatusCode ISF::FastCaloSimSvc::initialize()
{
   ATH_MSG_INFO ( m_screenOutputPrefix << "Initializing ...");

   // access tools and store them
86
   if ( retrieveTools<ICaloCellMakerTool>(m_caloCellMakerTools_setup).isFailure() )
87
        return StatusCode::FAILURE;
88
   if ( retrieveTools<ICaloCellMakerTool>(m_caloCellMakerTools_simulate).isFailure() )
89
        return StatusCode::FAILURE;
90
   if ( retrieveTools<ICaloCellMakerTool>(m_caloCellMakerTools_release).isFailure() )
91
92
        return StatusCode::FAILURE;

93
   if (m_doPunchThrough && m_punchThroughTool.retrieve().isFailure() )
94
95
96
   {
     ATH_MSG_ERROR (m_punchThroughTool.propertyName() << ": Failed to retrieve tool " << m_punchThroughTool.type());
     return StatusCode::FAILURE;
97
   }
98

99
   // Get TimedExtrapolator
100
101
   if (!m_extrapolator.empty() && m_extrapolator.retrieve().isFailure())
     return StatusCode::FAILURE;
102

103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
   ATH_MSG_DEBUG( m_screenOutputPrefix << " Output CaloCellContainer Name " << m_caloCellsOutputName );
   if (m_ownPolicy==SG::OWN_ELEMENTS){
       ATH_MSG_INFO( m_screenOutputPrefix << "...will OWN its cells." );
     }
     else
     {
       ATH_MSG_INFO( m_screenOutputPrefix << "...will VIEW its cells." );
     }

   if (m_caloCellHack){
     ATH_MSG_WARNING( m_screenOutputPrefix << " CaloCellContainer: " << m_caloCellsOutputName << "will be read in and modified !. To be used with care. " );
   }

  if ( m_particleBroker.retrieve().isFailure()) {
    ATH_MSG_FATAL ("Could not retrieve ISF ParticleBroker service " <<m_particleBroker);
    return StatusCode::FAILURE;
  }

   return StatusCode::SUCCESS;
}

/** framework methods */
StatusCode ISF::FastCaloSimSvc::finalize()
{
    ATH_MSG_INFO ( m_screenOutputPrefix << "Finalizing ...");
    return StatusCode::SUCCESS;
}

StatusCode ISF::FastCaloSimSvc::setupEvent()
132
{
133
  ATH_MSG_DEBUG ( m_screenOutputPrefix << "setup Event");
134

135
  if (!m_caloCellHack) {
136

137
138
139
140
    m_theContainer = new CaloCellContainer(static_cast<SG::OwnershipPolicy>(m_ownPolicy));

    StatusCode sc=StatusCode::SUCCESS;
    sc=evtStore()->record(m_theContainer,m_caloCellsOutputName);
141

142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
    if (sc.isFailure())  {
      ATH_MSG_FATAL( m_screenOutputPrefix << "cannot record CaloCellContainer " << m_caloCellsOutputName );
      return StatusCode::FAILURE;
    }

    // also symLink as INavigable4MomentumCollection!
    INavigable4MomentumCollection* theNav4Coll = 0;
    sc = evtStore()->symLink(m_theContainer,theNav4Coll);

    if (sc.isFailure()) {
      ATH_MSG_WARNING( m_screenOutputPrefix << "Error symlinking CaloCellContainer to INavigable4MomentumCollection " );
      return StatusCode::FAILURE;
    }
  }
  else {
    // take CaloCellContainer from input and cast away constness
    const CaloCellContainer * theConstContainer ;
159

160
161
162
163
164
165
166
167
168
169
170
171
172
173
    StatusCode sc=StatusCode::SUCCESS;
    sc=evtStore()->retrieve(theConstContainer,m_caloCellsOutputName);
    if (sc.isFailure() || theConstContainer==0)
    {
      ATH_MSG_FATAL( m_screenOutputPrefix << "Could not retrieve CaloCellContainer " << m_caloCellsOutputName );
      return StatusCode::FAILURE;
    }
    m_theContainer = const_cast<CaloCellContainer *> (theConstContainer);
  }

  // loop on setup tools
  ToolHandleArray<ICaloCellMakerTool>::iterator itrTool=m_caloCellMakerTools_setup.begin();
  ToolHandleArray<ICaloCellMakerTool>::iterator endTool=m_caloCellMakerTools_setup.end();
  for (;itrTool!=endTool;++itrTool){
174
    ATH_MSG_DEBUG( m_screenOutputPrefix << "Calling tool " << itrTool->name() );
175
    std::string chronoName=this->name()+"_"+ itrTool->name();
176

177
178
179
180
181
182
183
184
185
    if (m_chrono) m_chrono -> chronoStart( chronoName);
    StatusCode sc = (*itrTool)->process(m_theContainer);
    if (m_chrono) {
      m_chrono -> chronoStop( chronoName );
      ATH_MSG_DEBUG( m_screenOutputPrefix << "Chrono stop : delta " << m_chrono->chronoDelta (chronoName,IChronoStatSvc::USER ) * CLHEP::microsecond / CLHEP::second << " second " );
    }

    if (sc.isFailure()) {
      ATH_MSG_ERROR( m_screenOutputPrefix << "Error executing tool " << itrTool->name() );
186
    }
187
188
189
190
191
192
193
194
195
196
197
  }

  // loop on simulate tools
  itrTool=m_caloCellMakerTools_simulate.begin();
  endTool=m_caloCellMakerTools_simulate.end();
  for (;itrTool!=endTool;++itrTool) {
    FastShowerCellBuilderTool* fcs=dynamic_cast< FastShowerCellBuilderTool* >(&(*(*itrTool)));
    if(fcs) {
      if(fcs->setupEvent().isFailure()) {
        ATH_MSG_ERROR( m_screenOutputPrefix << "Error executing tool " << itrTool->name() << " in setupEvent");
        return StatusCode::FAILURE;
198
      }
199
200
201
202
203
204
205
    }
  }

  return StatusCode::SUCCESS;
}

StatusCode ISF::FastCaloSimSvc::releaseEvent()
206
{
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
  ATH_MSG_DEBUG ( m_screenOutputPrefix << "release Event");

  // the return value
  StatusCode sc = StatusCode::SUCCESS;

  // (1.) batch processing of all particles from McTruth
  //       (for MC12 parametrization)
  //
  if ( m_batchProcessMcTruth) {
    // -> run the FastShowerCellBuilder tools
    //        (in Python they should be configured to pick up the modified truth collection)
    ToolHandleArray<ICaloCellMakerTool>::iterator itrTool=m_caloCellMakerTools_simulate.begin();
    ToolHandleArray<ICaloCellMakerTool>::iterator endTool=m_caloCellMakerTools_simulate.end();
    for (;itrTool!=endTool;++itrTool) {
      FastShowerCellBuilderTool* fcs=dynamic_cast< FastShowerCellBuilderTool* >(&(*(*itrTool)));
      if(!fcs) {
223
        ATH_MSG_WARNING( m_screenOutputPrefix << "tool " << itrTool->name()<< "is not a FastShowerCellBuilderTool" );
224
225
        continue;
      }
226
227

      ATH_MSG_VERBOSE( m_screenOutputPrefix << "Calling tool " << itrTool->name() );
228
229

      if( fcs->process(m_theContainer).isFailure()) {
230
        ATH_MSG_WARNING( m_screenOutputPrefix << "batch simulation of FastCaloSim particles failed" );
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
        sc = StatusCode::FAILURE;
      }
    }
  } // <-- end of batch particle processing

  // (2.) finalize simulation tools in a loop
  //
  ToolHandleArray<ICaloCellMakerTool>::iterator itrTool=m_caloCellMakerTools_simulate.begin();
  ToolHandleArray<ICaloCellMakerTool>::iterator endTool=m_caloCellMakerTools_simulate.end();
  for (;itrTool!=endTool;++itrTool) {
    FastShowerCellBuilderTool* fcs=dynamic_cast< FastShowerCellBuilderTool* >(&(*(*itrTool)));
    if(fcs) {
      if(fcs->releaseEvent().isFailure()) {
        ATH_MSG_ERROR( m_screenOutputPrefix << "Error executing tool " << itrTool->name() << " in releaseEvent");
        return StatusCode::FAILURE;
      }
    }
  }
249

250
251
252
253
254
  // (3.) run release tools in a loop
  //
  itrTool=m_caloCellMakerTools_release.begin();
  endTool=m_caloCellMakerTools_release.end();
  for (;itrTool!=endTool;++itrTool){
255
    ATH_MSG_DEBUG( m_screenOutputPrefix << "Calling tool " << itrTool->name() );
256
    std::string chronoName=this->name()+"_"+ itrTool->name();
257

258
259
260
261
262
263
264
265
266
267
268
269
    if (m_chrono) m_chrono -> chronoStart( chronoName);
    sc = (*itrTool)->process(m_theContainer);
    if (m_chrono) {
      m_chrono -> chronoStop( chronoName );
      ATH_MSG_DEBUG( m_screenOutputPrefix << "Chrono stop : delta " << m_chrono->chronoDelta (chronoName,IChronoStatSvc::USER ) * CLHEP::microsecond / CLHEP::second << " second " );
    }

    if (sc.isFailure()) {
      ATH_MSG_ERROR( m_screenOutputPrefix << "Error executing tool " << itrTool->name() );
    }
  }

270
  return StatusCode::SUCCESS;
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
}


/** Simulation Call */
StatusCode ISF::FastCaloSimSvc::simulate(const ISF::ISFParticle& isfp)
{
  // read the particle's barcode
  Barcode::ParticleBarcode bc = isfp.barcode();
//lets do punch-through here
//----------------------------------------------------------

  // punch-through simulation

  if (m_doPunchThrough) {
    // call punch-through simulation
286
    const ISF::ISFParticleVector* isfpVec = m_punchThroughTool->computePunchThroughParticles(isfp);
287
288
289

    // add punch-through particles to the ISF particle broker
    if (isfpVec) {
290
291
      ISF::ISFParticleVector::const_iterator partIt    = isfpVec->begin();
      ISF::ISFParticleVector::const_iterator partItEnd = isfpVec->end();
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
      for ( ; partIt!=partItEnd; ++partIt) {
        m_particleBroker->push( *partIt, &isfp);
      }
    }
  }

  // (a.) batch process mode, ignore the incoming particle for now
  if ( m_batchProcessMcTruth) {
    ATH_MSG_DEBUG("particle is ignored now, will run Calo simulation using ID McTruth at the end of the event");
    return StatusCode::SUCCESS;
  }
  // (b.) throw away particles with undefined Barcode if m_simulateUndefinedBCs==False
  else if ( (!m_simulateUndefinedBCs) && (bc == Barcode::fUndefinedBarcode)) {
    ATH_MSG_DEBUG("particle has undefined barcode, will not simulate it");
    return StatusCode::SUCCESS;
  }
  // (c.) individual particle processing
309
310
  ATH_MSG_DEBUG("particle is simulated individually");
  return processOneParticle( isfp);
311
312
313
314
315
316
317
318
319
}


StatusCode ISF::FastCaloSimSvc::processOneParticle( const ISF::ISFParticle& isfp) {
  ATH_MSG_VERBOSE ( m_screenOutputPrefix << "Simulating pdgid = "<< isfp.pdgCode());

  ToolHandleArray<ICaloCellMakerTool>::iterator itrTool=m_caloCellMakerTools_simulate.begin();
  ToolHandleArray<ICaloCellMakerTool>::iterator endTool=m_caloCellMakerTools_simulate.end();

320
321
322
323
324
325
326
  std::vector<Trk::HitInfo>* hitVector= caloHits(isfp);

  if (!hitVector || !hitVector->size()) {
    ATH_MSG_WARNING ( "ISF_FastCaloSim: no hits in calo");
    return StatusCode::FAILURE;
  }

327
328
329
330
  // loop on tools
  for (;itrTool!=endTool;++itrTool) {
    FastShowerCellBuilderTool* fcs=dynamic_cast< FastShowerCellBuilderTool* >(&(*(*itrTool)));
    if(!fcs) {
331
      ATH_MSG_WARNING( m_screenOutputPrefix << "tool " << itrTool->name()<< "is not a FastShowerCellBuilderTool" );
332
333
      continue;
    }
334
335

    ATH_MSG_VERBOSE( m_screenOutputPrefix << "Calling tool " << itrTool->name() );
336
    std::string chronoName=this->name()+"_"+ itrTool->name();
337

338
339
340
    if (m_chrono) m_chrono->chronoStart( chronoName);

    //sc = (*itrTool)->process(m_theContainer);
341
342
    if(fcs->process_particle(m_theContainer,hitVector,
			     isfp.momentum(),isfp.mass(),isfp.pdgCode(),isfp.barcode()).isFailure()) {
343
      ATH_MSG_WARNING( m_screenOutputPrefix << "simulation of particle pdgid=" << isfp.pdgCode()<< " failed" );
344
      return StatusCode::FAILURE;
345
346
347
348
349
350
    }

    if (m_chrono) m_chrono->chronoStop( chronoName );
    //ATH_MSG_DEBUG( m_screenOutputPrefix << "Chrono stop : delta " << m_chrono->chronoDelta (chronoName,IChronoStatSvc::USER ) * CLHEP::microsecond / CLHEP::second << " second " );
  } //end of for-loop

351
352
353
354
355
  if(hitVector) {
    for(std::vector<Trk::HitInfo>::iterator it = hitVector->begin();it < hitVector->end();++it)  {
      if((*it).trackParms) {
        delete (*it).trackParms;
        (*it).trackParms=0;
356
      }
357
358
    }
    delete hitVector;
359
  }
360

361
  //  ATH_MSG_VERBOSE ( m_screenOutputPrefix << "kill the particle in the end");
362
  return StatusCode::SUCCESS;
363
364
}

365
366
367
368
369
370

std::vector<Trk::HitInfo>* ISF::FastCaloSimSvc::caloHits(const ISF::ISFParticle& isp) const
{
  // Start calo extrapolation
  ATH_MSG_VERBOSE ("[ fastCaloSim transport ] processing particle "<<isp.pdgCode() );

371
  std::vector<Trk::HitInfo>*     hitVector =  new std::vector<Trk::HitInfo>;
372
373
374
375
376
377
378
379
380
381
382
383

  int  absPdg         = abs(isp.pdgCode());
  bool charged        = isp.charge()*isp.charge() > 0 ;

  // particle Hypothesis for the extrapolation

  Trk::ParticleHypothesis pHypothesis = m_pdgToParticleHypothesis.convert(isp.pdgCode(),isp.charge());

  // geantinos not handled by PdgToParticleHypothesis - fix there
  if ( absPdg == 999 ) pHypothesis = Trk::geantino;

  // choose the extrapolator
384
  //const Trk::ITimedExtrapolator* processor = &(*m_extrapolator);
385
386
387
388

  // input parameters : curvilinear parameters
  Trk::CurvilinearParameters inputPar(isp.position(),isp.momentum(),isp.charge());

389
390
  // stable vs. unstable check : ADAPT for FASTCALOSIM
  //double freepath = ( !m_particleDecayHelper.empty()) ? m_particleDecayHelper->freePath(isp) : - 1.;
391
392
393
394
395
396
397
398
399
400
  double freepath = -1.;
  ATH_MSG_VERBOSE( "[ fatras transport ] Particle free path : " << freepath);
  // path limit -> time limit  ( TODO : extract life-time directly from decay helper )
  double tDec = freepath > 0. ? freepath : -1.;
  int decayProc = 0;

  /* stable particles only for the moment
  // beta calculated here for further use in validation
  double mass = m_particleMasses.mass[pHypothesis];
  double mom = isp.momentum().mag();
401
  double beta = mom/sqrt(mom*mom+mass*mass);
402
403
404

  if ( tDec>0.) {
    tDec = tDec/beta/CLHEP::c_light + isp.timeStamp();
405
    decayProc = 201;
406
407
408
409
  }
  */

  Trk::TimeLimit timeLim(tDec,isp.timeStamp(),decayProc);
410

411
412
413
414
415
416
417
418
419
420
421
422
  // prompt decay
  //if ( freepath>0. && freepath<0.01 ) {
  //  if (!m_particleDecayHelper.empty()) {
  //    ATH_MSG_VERBOSE( "[ fatras transport ] Decay is triggered for input particle.");
  //    m_particleDecayHelper->decay(isp);
  //  }
  //  return 0;
  //}

  // presample interactions - ADAPT FOR FASTCALOSIM ( non-interacting )
  Trk::PathLimit pathLim(-1.,0);
  //if (absPdg!=999 && pHypothesis<99) pathLim = m_samplingTool->sampleProcess(mom,isp.charge(),pHypothesis);
423

424
  Trk::GeometrySignature nextGeoID = static_cast<Trk::GeometrySignature>(isp.nextGeoID());
425
426

  // save Calo entry hit (fallback info)
427
428
  hitVector->push_back(Trk::HitInfo(inputPar.clone(),isp.timeStamp(),nextGeoID,0.));

429
430
431
432
  const Trk::TrackParameters* eParameters = 0;

  if ( !charged ) {

433
434
    eParameters = m_extrapolator->transportNeutralsWithPathLimit(inputPar,pathLim,timeLim,Trk::alongMomentum,pHypothesis,hitVector,nextGeoID);

435
  } else {
436

437
438
439
440
    eParameters = m_extrapolator->extrapolateWithPathLimit(inputPar,pathLim,timeLim,Trk::alongMomentum,pHypothesis,hitVector,nextGeoID);

  }
  // save Calo exit hit (fallback info)
441
  if (eParameters) hitVector->push_back(Trk::HitInfo(eParameters,timeLim.time,nextGeoID,0.));
442
443
444
445
446
447

  ATH_MSG_VERBOSE( "[ fastCaloSim transport ] number of intersections "<< hitVector->size());

  return hitVector;

}