McEventCollectionCnv_p3.cxx 12.9 KB
Newer Older
1
2
3
///////////////////////// -*- C++ -*- /////////////////////////////

/*
4
  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
5
6
*/

7
// McEventCollectionCnv_p3.cxx
8
9
// Implementation file for class McEventCollectionCnv_p3
// Author: S.Binet<binet@cern.ch>
10
///////////////////////////////////////////////////////////////////
11
12
13
14
15
16
17
18
19
20
21


// STL includes
#include <utility>
#include <cmath>

// GeneratorObjectsTPCnv includes
#include "GeneratorObjectsTPCnv/McEventCollectionCnv_p3.h"
#include "HepMcDataPool.h"


22
///////////////////////////////////////////////////////////////////
23
// Constructors
24
///////////////////////////////////////////////////////////////////
25
26
27


McEventCollectionCnv_p3::McEventCollectionCnv_p3() :
28
  Base_t( )
29
30
31
{}

McEventCollectionCnv_p3::McEventCollectionCnv_p3( const McEventCollectionCnv_p3& rhs ) :
32
  Base_t( rhs )
33
34
{}

35
McEventCollectionCnv_p3&
36
37
38
39
McEventCollectionCnv_p3::operator=( const McEventCollectionCnv_p3& rhs )
{
  if ( this != &rhs ) {
    Base_t::operator=( rhs );
40
  }
41
42
43
  return *this;
}

44
///////////////////////////////////////////////////////////////////
45
// Destructor
46
///////////////////////////////////////////////////////////////////
47
48
49
50
51
52
53

McEventCollectionCnv_p3::~McEventCollectionCnv_p3()
{
}


void McEventCollectionCnv_p3::persToTrans( const McEventCollection_p3* persObj,
54
55
                                           McEventCollection* transObj,
                                           MsgStream& msg )
56
57
{
  msg << MSG::DEBUG << "Loading McEventCollection from persistent state..."
58
      << endmsg;
59
60
61

  // elements are managed by DataPool
  transObj->clear(SG::VIEW_ELEMENTS);
62
63
64
65
  HepMC::DataPool datapools;
  const unsigned int nVertices = persObj->m_genVertices.size();
  if ( datapools.vtx.capacity() - datapools.vtx.allocated() < nVertices ) {
    datapools.vtx.reserve( datapools.vtx.allocated() + nVertices );
66
  }
67
68
69
  const unsigned int nParts = persObj->m_genParticles.size();
  if ( datapools.part.capacity() - datapools.part.allocated() < nParts ) {
    datapools.part.reserve( datapools.part.allocated() + nParts );
70
71
  }
  const unsigned int nEvts = persObj->m_genEvents.size();
72
73
  if ( datapools.evt.capacity() - datapools.evt.allocated() < nEvts ) {
    datapools.evt.reserve( datapools.evt.allocated() + nEvts );
74
  }
75

76
  transObj->reserve( nEvts );
77
78
79
80
81
  for ( std::vector<GenEvent_p3>::const_iterator
          itr = persObj->m_genEvents.begin(),
          itrEnd = persObj->m_genEvents.end();
        itr != itrEnd;
        ++itr ) {
82
83
    const GenEvent_p3& persEvt = *itr;

84
    HepMC::GenEvent * genEvt        = datapools.getGenEvent();
85
86
87
88
89
90
91
92
#ifdef HEPMC3
    genEvt->add_attribute("signal_process_id",std::make_shared<HepMC3::IntAttribute>(persEvt.m_signalProcessId));
    genEvt->set_event_number(persEvt.m_eventNbr);
    genEvt->add_attribute("event_scale",std::make_shared<HepMC3::DoubleAttribute>(persEvt.m_eventScale));
    genEvt->add_attribute("alphaQCD",std::make_shared<HepMC3::DoubleAttribute>(persEvt.m_alphaQCD));
    genEvt->add_attribute("alphaQED",std::make_shared<HepMC3::DoubleAttribute>(persEvt.m_alphaQED));
    genEvt->weights()= persEvt.m_weights;
    genEvt->add_attribute("random_states",std::make_shared<HepMC3::VectorLongIntAttribute>(persEvt.m_randomStates));
93
94
95
    transObj->push_back( genEvt );

    ParticlesMap_t partToEndVtx( (persEvt.m_particlesEnd- persEvt.m_particlesBegin)/2 );
96

97
98
99
100
101
    // create the vertices
    const unsigned int endVtx = persEvt.m_verticesEnd;
    for ( unsigned int iVtx= persEvt.m_verticesBegin; iVtx != endVtx; ++iVtx ) {
      createGenVertex( *persObj, persObj->m_genVertices[iVtx],partToEndVtx, datapools, genEvt );
    } 
102
#else
103
104
105
106
107
108
109
110
111
112
113
114
115
    genEvt->m_signal_process_id     = persEvt.m_signalProcessId;
    genEvt->m_event_number          = persEvt.m_eventNbr;
    genEvt->m_event_scale           = persEvt.m_eventScale;
    genEvt->m_alphaQCD              = persEvt.m_alphaQCD;
    genEvt->m_alphaQED              = persEvt.m_alphaQED;
    genEvt->m_signal_process_vertex = 0;
    genEvt->m_weights               = persEvt.m_weights;
    genEvt->m_random_states         = persEvt.m_randomStates;
    genEvt->m_vertex_barcodes.clear();
    genEvt->m_particle_barcodes.clear();
    genEvt->m_pdf_info = 0;         //> not available at that time...

    transObj->push_back( genEvt );
116
117

    // create a temporary map associating the barcode of an end-vtx to its
118
119
120
121
    // particle.
    // As not all particles are stable (d'oh!) we take 50% of the number of
    // particles as an initial size of the hash-map (to prevent re-hash)
    ParticlesMap_t partToEndVtx( (persEvt.m_particlesEnd-
122
123
                                  persEvt.m_particlesBegin)/2 );

124
125
126
    // create the vertices
    const unsigned int endVtx = persEvt.m_verticesEnd;
    for ( unsigned int iVtx= persEvt.m_verticesBegin; iVtx != endVtx; ++iVtx ) {
127
128
129
      genEvt->add_vertex( createGenVertex( *persObj,
                                           persObj->m_genVertices[iVtx],
                                           partToEndVtx,
130
                                           datapools ) );
131
    } //> end loop over vertices
132
#endif
133

134
135
136
    // set the signal process vertex
    const int sigProcVtx = persEvt.m_signalProcessVtx;
    if ( sigProcVtx != 0 ) {
137
      HepMC::set_signal_process_vertex(genEvt, HepMC::barcode_to_vertex(genEvt, sigProcVtx ) );
138
139
140
    }

    // connect particles to their end vertices
141
142
143
144
145
    for ( ParticlesMap_t::iterator
            p = partToEndVtx.begin(),
            endItr = partToEndVtx.end();
          p != endItr;
          ++p ) {
146
      auto decayVtx=HepMC::barcode_to_vertex(genEvt, p->second );
147
      if ( decayVtx ) {
148
        decayVtx->add_particle_in( p->first );
149
      } else {
150
151
        msg << MSG::ERROR
            << "GenParticle points to null end vertex !!"
152
            << endmsg;
153
154
155
156
157
      }
    }
  } //> end loop over m_genEvents

  msg << MSG::DEBUG << "Loaded McEventCollection from persistent state [OK]"
158
      << endmsg;
159
160
161
162

  return;
}

163
164
165
void McEventCollectionCnv_p3::transToPers( const McEventCollection* /*transObj*/,
                                           McEventCollection_p3* /*persObj*/,
                                           MsgStream& msg )
166
167
168
{
  msg << MSG::ERROR
      << "This transient-to-persistent converter method has been RETIRED !!"
169
      << endmsg
170
      << "You are not supposed to end-up here ! Go away !"
171
      << endmsg;
172
173
174
175
176
177

  throw std::runtime_error( "Retired McEventCollectionCnv_p3::transToPers() !!" );
  return;
}


178
HepMC::GenVertexPtr
179
McEventCollectionCnv_p3::createGenVertex( const McEventCollection_p3& persEvt,
180
181
                                          const GenVertex_p3& persVtx,
                                          ParticlesMap_t& partToEndVtx,
182
                                          HepMC::DataPool& datapools, HepMC::GenEvent* parent ) const
183
{
184
  HepMC::GenVertexPtr vtx = datapools.getGenVertex();
185
  if (parent) parent->add_vertex(vtx);
186
187
188

#ifdef HEPMC3
  vtx->set_position( HepMC::FourVector(persVtx.m_x,persVtx.m_y,persVtx.m_z,persVtx.m_t) );
189
  vtx->set_status(persVtx.m_id);
190
191
  vtx->add_attribute("weights",std::make_shared<HepMC3::VectorFloatAttribute>(persVtx.m_weights));
  vtx->add_attribute("barcode",std::make_shared<HepMC3::IntAttribute>(persVtx.m_barcode));
192
193
194
195
196
197
198
199
200
201
202
  // handle the in-going (orphans) particles
  //Is this needed for HEPMC3?
  const unsigned int nPartsIn = persVtx.m_particlesIn.size();
  for ( unsigned int i = 0; i != nPartsIn; ++i ) {
    createGenParticle( persEvt.m_genParticles[persVtx.m_particlesIn[i]], partToEndVtx, datapools );
  }
  // now handle the out-going particles
  const unsigned int nPartsOut = persVtx.m_particlesOut.size();
  for ( unsigned int i = 0; i != nPartsOut; ++i ) {
     createGenParticle( persEvt.m_genParticles[persVtx.m_particlesOut[i]], partToEndVtx, datapools, vtx );
  }
203
#else
204
205
206
207
208
209
210
211
212
  vtx->m_position.setX( persVtx.m_x );
  vtx->m_position.setY( persVtx.m_y );
  vtx->m_position.setZ( persVtx.m_z );
  vtx->m_position.setT( persVtx.m_t );
  vtx->m_particles_in.clear();
  vtx->m_particles_out.clear();
  vtx->m_id      = persVtx.m_id;
  vtx->m_weights.m_weights.reserve( persVtx.m_weights.size() );
  vtx->m_weights.m_weights.assign ( persVtx.m_weights.begin(),
213
                                    persVtx.m_weights.end() );
214
215
  vtx->m_event   = 0;
  vtx->m_barcode = persVtx.m_barcode;
216

217
218
219
  // handle the in-going (orphans) particles
  const unsigned int nPartsIn = persVtx.m_particlesIn.size();
  for ( unsigned int i = 0; i != nPartsIn; ++i ) {
220
221
222
    createGenParticle( persEvt.m_genParticles[persVtx.m_particlesIn[i]],
                       partToEndVtx,
                       datapools );
223
  }
224

225
226
227
  // now handle the out-going particles
  const unsigned int nPartsOut = persVtx.m_particlesOut.size();
  for ( unsigned int i = 0; i != nPartsOut; ++i ) {
228
229
230
    vtx->add_particle_out( createGenParticle( persEvt.m_genParticles[persVtx.m_particlesOut[i]],
                                              partToEndVtx,
                                              datapools ) );
231
  }
232
#endif
233
234
235
236

  return vtx;
}

237
HepMC::GenParticlePtr
238
McEventCollectionCnv_p3::createGenParticle( const GenParticle_p3& persPart,
239
                                            ParticlesMap_t& partToEndVtx,
240
                                            HepMC::DataPool& datapools, HepMC::GenVertexPtr parent ) const
241
{
242
  HepMC::GenParticlePtr p    = datapools.getGenParticle();
243
  if (parent) parent->add_particle_out(p);
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279

#ifdef HEPMC3
  p->set_pdg_id(persPart.m_pdgId);
  p->set_status(persPart.m_status);
  // Note: do the E calculation in extended (long double) precision.
  // That happens implicitly on x86 with optimization on; saying it
  // explicitly ensures that we get the same results with and without
  // optimization.  (If this is a performance issue for platforms
  // other than x86, one could change to double for those platforms.)
  double temp_e=0.0;
  if ( 0 == persPart.m_recoMethod ) {
    temp_e = std::sqrt( (long double)(persPart.m_px)*persPart.m_px +
                          (long double)(persPart.m_py)*persPart.m_py +
                          (long double)(persPart.m_pz)*persPart.m_pz +
                          (long double)(persPart.m_m) *persPart.m_m );
  } else {
    const int signM2 = ( persPart.m_m >= 0. ? 1 : -1 );
    const double persPart_ene =
      std::sqrt( std::abs((long double)(persPart.m_px)*persPart.m_px +
                (long double)(persPart.m_py)*persPart.m_py +
                (long double)(persPart.m_pz)*persPart.m_pz +
                signM2* (long double)(persPart.m_m)* persPart.m_m));
    const int signEne = ( persPart.m_recoMethod == 1 ? 1 : -1 );
    temp_e=signEne * persPart_ene;
  }
  p->set_momentum(HepMC::FourVector(persPart.m_px,persPart.m_py,persPart.m_pz,temp_e));
  // setup flow
  // fillin' the flow
  std::vector<int> flows;
  const unsigned int nFlow = persPart.m_flow.size();
  for ( unsigned int iFlow= 0; iFlow != nFlow; ++iFlow ) {
  flows.push_back(persPart.m_flow[iFlow].second );
  }
  //We construct it here as vector w/o gaps.
  p->add_attribute("flows", std::make_shared<HepMC3::VectorIntAttribute>(flows));
#else
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
  p->m_pdg_id              = persPart.m_pdgId;
  p->m_status              = persPart.m_status;
  p->m_polarization.m_theta= static_cast<double>(persPart.m_thetaPolarization);
  p->m_polarization.m_phi  = static_cast<double>(persPart.m_phiPolarization  );
  p->m_production_vertex   = 0;
  p->m_end_vertex          = 0;
  p->m_barcode             = persPart.m_barcode;

  // Note: do the E calculation in extended (long double) precision.
  // That happens implicitly on x86 with optimization on; saying it
  // explicitly ensures that we get the same results with and without
  // optimization.  (If this is a performance issue for platforms
  // other than x86, one could change to double for those platforms.)
  if ( 0 == persPart.m_recoMethod ) {
    p->m_momentum.setPx( persPart.m_px );
    p->m_momentum.setPy( persPart.m_py );
    p->m_momentum.setPz( persPart.m_pz );
297
    double temp_e = std::sqrt( (long double)(persPart.m_px)*persPart.m_px +
298
299
                          (long double)(persPart.m_py)*persPart.m_py +
                          (long double)(persPart.m_pz)*persPart.m_pz +
300
301
302
303
304
                          (long double)(persPart.m_m) *persPart.m_m );
    p->m_momentum.setE( temp_e );
  } else {
    const int signM2 = ( persPart.m_m >= 0. ? 1 : -1 );
    const double persPart_ene =
305
      std::sqrt( std::abs((long double)(persPart.m_px)*persPart.m_px +
306
307
                (long double)(persPart.m_py)*persPart.m_py +
                (long double)(persPart.m_pz)*persPart.m_pz +
308
                signM2* (long double)(persPart.m_m)* persPart.m_m));
309
310
    const int signEne = ( persPart.m_recoMethod == 1 ? 1 : -1 );
    p->m_momentum.set( persPart.m_px,
311
312
313
                       persPart.m_py,
                       persPart.m_pz,
                       signEne * persPart_ene );
314
315
316
317
  }

  // setup flow
  const unsigned int nFlow = persPart.m_flow.size();
318
  p->m_flow.clear();
319
  for ( unsigned int iFlow= 0; iFlow != nFlow; ++iFlow ) {
320
321
    p->m_flow.set_icode( persPart.m_flow[iFlow].first,
                         persPart.m_flow[iFlow].second );
322
  }
323
#endif
324
325
326
327
328
329
330

  if ( persPart.m_endVtx != 0 ) {
    partToEndVtx[p] = persPart.m_endVtx;
  }

  return p;
}