NativeFastCaloSimSvc.cxx 13.1 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
/*
  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/

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

// class header include
#include "ISF_FastCaloSimParametrization/NativeFastCaloSimSvc.h"

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

// ISF includes
#include "ISF_Event/ISFParticle.h"
#include "ISF_Event/ISFParticleContainer.h"

// 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"
// PunchThrough Tool
#include "ISF_FastCaloSimInterfaces/IPunchThroughTool.h"

/** Constructor **/
ISF::NativeFastCaloSimSvc::NativeFastCaloSimSvc(const std::string& name,ISvcLocator* svc) :
  BaseSimulationSvc(name, svc),
  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)
{
50
  // where to go
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
  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) ;
  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" );
}

70
ISF::NativeFastCaloSimSvc::~NativeFastCaloSimSvc()
71
72
73
74
75
76
77
78
{}

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

   // access tools and store them
79
   if ( retrieveTools<ICaloCellMakerTool>(m_caloCellMakerTools_setup).isFailure() )
80
        return StatusCode::FAILURE;
81
   if ( retrieveTools<ICaloCellMakerTool>(m_caloCellMakerTools_simulate).isFailure() )
82
        return StatusCode::FAILURE;
83
   if ( retrieveTools<ICaloCellMakerTool>(m_caloCellMakerTools_release).isFailure() )
84
85
        return StatusCode::FAILURE;

86
   if (m_doPunchThrough && m_punchThroughTool.retrieve().isFailure() )
87
88
89
   {
     ATH_MSG_ERROR (m_punchThroughTool.propertyName() << ": Failed to retrieve tool " << m_punchThroughTool.type());
     return StatusCode::FAILURE;
90
91
   }

92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
   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::NativeFastCaloSimSvc::finalize()
{
    ATH_MSG_INFO ( m_screenOutputPrefix << "Finalizing ...");
    return StatusCode::SUCCESS;
}

StatusCode ISF::NativeFastCaloSimSvc::setupEvent()
121
{
122
  ATH_MSG_DEBUG ( m_screenOutputPrefix << "setup Event");
123

124
  if (!m_caloCellHack) {
125

126
127
128
129
    m_theContainer = new CaloCellContainer(static_cast<SG::OwnershipPolicy>(m_ownPolicy));

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

131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
    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 ;
148

149
150
151
152
153
154
155
156
157
158
159
160
161
162
    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){
163
    ATH_MSG_DEBUG( m_screenOutputPrefix << "Calling tool " << itrTool->name() );
164
    std::string chronoName=this->name()+"_"+ itrTool->name();
165

166
167
168
169
170
171
172
173
174
    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() );
175
    }
176
177
178
179
180
181
182
183
184
185
186
187
  }

  /*
  // 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;
188
      }
189
190
191
192
193
194
195
196
    }
  }
  */

  return StatusCode::SUCCESS;
}

StatusCode ISF::NativeFastCaloSimSvc::releaseEvent()
197
{
198
199
200
201
202
203
204
205
206
207
208
  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)
209

210
211
212
213
214
215
216
    //ZH commented out these to avoid warnings at the moment
    //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) {
217
        ATH_MSG_WARNING( m_screenOutputPrefix << "tool " << itrTool->name()<< "is not a FastShowerCellBuilderTool" );
218
219
        continue;
      }
220
221

      ATH_MSG_VERBOSE( m_screenOutputPrefix << "Calling tool " << itrTool->name() );
222
223

      if( fcs->process(m_theContainer).isFailure()) {
224
        ATH_MSG_WARNING( m_screenOutputPrefix << "batch simulation of FastCaloSim particles failed" );
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
        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;
      }
    }
  }
  */
246

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

255
256
257
258
259
260
261
262
263
264
265
266
    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() );
    }
  }

267
  return StatusCode::SUCCESS;
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
}


/** Simulation Call */
StatusCode ISF::NativeFastCaloSimSvc::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
283
    const ISF::ISFParticleVector* isfpVec = m_punchThroughTool->computePunchThroughParticles(isfp);
284
285
286

    // add punch-through particles to the ISF particle broker
    if (isfpVec) {
287
288
      for (ISF::ISFParticle *particle : *isfpVec) {
        m_particleBroker->push( particle, &isfp);
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
      }
    }
  }

  // (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
304
305
  ATH_MSG_DEBUG("particle is simulated individually");
  return processOneParticle( isfp);
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
}


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

  StatusCode sc(StatusCode::SUCCESS);
  //ZH commented out these to avoid warnings at the moment
  //ToolHandleArray<ICaloCellMakerTool>::iterator itrTool=m_caloCellMakerTools_simulate.begin();
  //ToolHandleArray<ICaloCellMakerTool>::iterator endTool=m_caloCellMakerTools_simulate.end();

  /*
  // loop on tools
  for (;itrTool!=endTool;++itrTool) {
    FastShowerCellBuilderTool* fcs=dynamic_cast< FastShowerCellBuilderTool* >(&(*(*itrTool)));
    if(!fcs) {
322
      ATH_MSG_WARNING( m_screenOutputPrefix << "tool " << itrTool->name()<< "is not a FastShowerCellBuilderTool" );
323
324
      continue;
    }
325
326

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

329
330
331
332
333
334
335
336
337
338
    if (m_chrono) m_chrono->chronoStart( chronoName);

    HepMC::FourVector momentum(isfp.momentum().x(),isfp.momentum().y(),isfp.momentum().z(), sqrt( isfp.mass()*isfp.mass() + isfp.momentum().mag2()) );
    HepMC::GenParticle* part=new HepMC::GenParticle(momentum,isfp.pdgCode());
    HepMC::FourVector position(isfp.position().x(),isfp.position().y(),isfp.position().z(),isfp.timeStamp());
    HepMC::GenVertex vertex(position);
    vertex.add_particle_out( part ); // HepMC::GenVertex destructor will delete the particle

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

    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
  */

  //  ATH_MSG_VERBOSE ( m_screenOutputPrefix << "kill the particle in the end");
  return sc;
}