McEventCollectionCnv_p2.cxx 11.5 KB
Newer Older
1
2
3
///////////////////////// -*- C++ -*- /////////////////////////////

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

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


// STL includes
#include <utility>

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


21
///////////////////////////////////////////////////////////////////
22
// Constructors
23
///////////////////////////////////////////////////////////////////
24
25

McEventCollectionCnv_p2::McEventCollectionCnv_p2() :
26
  Base_t( )
27
28
29
{}

McEventCollectionCnv_p2::McEventCollectionCnv_p2( const McEventCollectionCnv_p2& rhs ) :
30
  Base_t( rhs )
31
32
{}

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

42
///////////////////////////////////////////////////////////////////
43
// Destructor
44
///////////////////////////////////////////////////////////////////
45
46
47
48
49
50
51

McEventCollectionCnv_p2::~McEventCollectionCnv_p2()
{
}


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

  // elements are managed by DataPool
  transObj->clear(SG::VIEW_ELEMENTS);
60
61
62
63
  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 );
64
  }
65
66
67
  const unsigned int nParts = persObj->m_genParticles.size();
  if ( datapools.part.capacity() - datapools.part.allocated() < nParts ) {
    datapools.part.reserve( datapools.part.allocated() + nParts );
68
69
  }
  const unsigned int nEvts = persObj->m_genEvents.size();
70
71
  if ( datapools.evt.capacity() - datapools.evt.allocated() < nEvts ) {
    datapools.evt.reserve( datapools.evt.allocated() + nEvts );
72
  }
73

74
75
76
  transObj->reserve( nEvts );
  const std::vector<GenEvent_p2>::const_iterator itrEnd = persObj->m_genEvents.end();
  for ( std::vector<GenEvent_p2>::const_iterator itr = persObj->m_genEvents.begin();
77
78
        itr != itrEnd;
        ++itr ) {
79
    const GenEvent_p2& persEvt = *itr;
80
    HepMC::GenEvent * genEvt        = datapools.getGenEvent();
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
#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));

    transObj->push_back( genEvt );

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

    // create the vertices
    const unsigned int endVtx = persEvt.m_verticesEnd;
    for ( unsigned int iVtx= persEvt.m_verticesBegin; iVtx != endVtx; ++iVtx ) {
      genEvt->add_vertex( createGenVertex( *persObj,
                                           persObj->m_genVertices[iVtx],
                                           partToEndVtx,
100
                                           datapools ) );
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
    } //> end loop over vertices

    // set the signal process vertex
    const int sigProcVtx = persEvt.m_signalProcessVtx;
    if ( sigProcVtx ) HepMC::set_signal_process_vertex(genEvt, HepMC::barcode_to_vertex(genEvt, sigProcVtx ) );

    // connect particles to their end vertices
    for ( auto p:  partToEndVtx) {
      auto decayVtx = HepMC::barcode_to_vertex(genEvt, p.second );
      if ( decayVtx ) {
        decayVtx->add_particle_in( p.first );
      } else {
        msg << MSG::ERROR<< "GenParticle points to null end vertex !!"<< endmsg;
      }
    }
#else
117
118
119
120
121
122
123
124
125
126
127
128
129
    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 );
130
131

    // create a temporary map associating the barcode of an end-vtx to its
132
133
134
135
    // 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-
136
137
                                  persEvt.m_particlesBegin)/2 );

138
139
140
    // create the vertices
    const unsigned int endVtx = persEvt.m_verticesEnd;
    for ( unsigned int iVtx= persEvt.m_verticesBegin; iVtx != endVtx; ++iVtx ) {
141
142
143
      genEvt->add_vertex( createGenVertex( *persObj,
                                           persObj->m_genVertices[iVtx],
                                           partToEndVtx,
144
                                           datapools ) );
145
    } //> end loop over vertices
146

147
148
149
150
151
152
153
154
155
    // set the signal process vertex
    const int sigProcVtx = persEvt.m_signalProcessVtx;
    if ( sigProcVtx != 0 ) {
      genEvt->set_signal_process_vertex( genEvt->barcode_to_vertex( sigProcVtx ) );
    }

    // connect particles to their end vertices
    const ParticlesMap_t::iterator endItr = partToEndVtx.end();
    for ( ParticlesMap_t::iterator p = partToEndVtx.begin();
156
157
          p != endItr;
          ++p ) {
Andrii Verbytskyi's avatar
Andrii Verbytskyi committed
158
      auto decayVtx = genEvt->barcode_to_vertex( p->second );
159
      if ( decayVtx ) {
160
        decayVtx->add_particle_in( p->first );
161
      } else {
162
163
        msg << MSG::ERROR
            << "GenParticle points to null end vertex !!"
164
            << endmsg;
165
166
      }
    }
167
#endif    
168
169
170
  } //> end loop over m_genEvents

  msg << MSG::DEBUG << "Loaded McEventCollection from persistent state [OK]"
171
      << endmsg;
172
173
174
175

  return;
}

176
177
178
void McEventCollectionCnv_p2::transToPers( const McEventCollection*,
                                           McEventCollection_p2*,
                                           MsgStream& msg )
179
180
{
  msg << MSG::DEBUG << "Creating persistent state of McEventCollection..."
181
      << endmsg;
182
183
184

  msg << MSG::ERROR
      << "This transient-to-persistent converter method has been RETIRED !!"
185
      << endmsg
186
      << "You are not supposed to end-up here ! Go away !"
187
      << endmsg;
188
189
190
191
192
193

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


194
HepMC::GenVertexPtr
195
McEventCollectionCnv_p2::createGenVertex( const McEventCollection_p2& persEvt,
196
                                          const GenVertex_p2& persVtx,
197
                                          ParticlesMap_t& partToEndVtx, HepMC::DataPool& datapools, HepMC::GenEvent* parent ) const
198
{
199
  HepMC::GenVertexPtr vtx = datapools.getGenVertex();
200
  if (parent) parent->add_vertex(vtx);
201
202
#ifdef HEPMC3
  vtx->set_position( HepMC::FourVector(persVtx.m_x,persVtx.m_y, persVtx.m_z, persVtx.m_t) );
203
  vtx->set_status(persVtx.m_id);
204
205
206
207
208
209
  vtx->add_attribute("weights",std::make_shared<HepMC3::VectorFloatAttribute>(persVtx.m_weights));
  vtx->add_attribute("barcode",std::make_shared<HepMC3::IntAttribute>(persVtx.m_barcode));

  // handle the in-going (orphans) particles
  const unsigned int nPartsIn = persVtx.m_particlesIn.size();
  for ( unsigned int i = 0; i != nPartsIn; ++i ) {
210
    createGenParticle( persEvt.m_genParticles[persVtx.m_particlesIn[i]], partToEndVtx, datapools );
211
212
213
214
215
  }

  // now handle the out-going particles
  const unsigned int nPartsOut = persVtx.m_particlesOut.size();
  for ( unsigned int i = 0; i != nPartsOut; ++i ) {
216
     createGenParticle( persEvt.m_genParticles[persVtx.m_particlesOut[i]], partToEndVtx, datapools, vtx );
217
218
  }
#else
219
220
221
222
223
224
225
226
227
  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(),
228
                                    persVtx.m_weights.end() );
229
230
  vtx->m_event   = 0;
  vtx->m_barcode = persVtx.m_barcode;
231

232
233
234
  // handle the in-going (orphans) particles
  const unsigned int nPartsIn = persVtx.m_particlesIn.size();
  for ( unsigned int i = 0; i != nPartsIn; ++i ) {
235
236
237
    createGenParticle( persEvt.m_genParticles[persVtx.m_particlesIn[i]],
                       partToEndVtx,
                       datapools );
238
  }
239

240
241
242
  // now handle the out-going particles
  const unsigned int nPartsOut = persVtx.m_particlesOut.size();
  for ( unsigned int i = 0; i != nPartsOut; ++i ) {
243
244
245
    vtx->add_particle_out( createGenParticle( persEvt.m_genParticles[persVtx.m_particlesOut[i]],
                                              partToEndVtx,
                                              datapools ) );
246
  }
247
#endif
248
249
250
251

  return vtx;
}

252
HepMC::GenParticlePtr
253
McEventCollectionCnv_p2::createGenParticle( const GenParticle_p2& persPart,
254
                                            ParticlesMap_t& partToEndVtx, HepMC::DataPool& datapools, HepMC::GenVertexPtr parent ) const
255
{
256
257
  HepMC::GenParticlePtr p    = datapools.getGenParticle();

258
  if (parent) parent->add_particle_out(p);
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
#ifdef HEPMC3
  p->set_momentum( HepMC::FourVector(persPart.m_px,persPart.m_py,persPart.m_pz, persPart.m_ene ));
  p->set_pdg_id(               persPart.m_pdgId);
  p->set_status(               persPart.m_status);
  p->add_attribute("phi",std::make_shared<HepMC3::DoubleAttribute>(persPart.m_phiPolarization));
  p->add_attribute("theta",std::make_shared<HepMC3::DoubleAttribute>(persPart.m_thetaPolarization));
  p->add_attribute("barcode",std::make_shared<HepMC3::IntAttribute>(persPart.m_barcode));

  // 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
277
278
279
280
281
282
283
284
285
286
287
288
289
290
  p->m_momentum.setPx( persPart.m_px  );
  p->m_momentum.setPy( persPart.m_py  );
  p->m_momentum.setPz( persPart.m_pz  );
  p->m_momentum.setE ( persPart.m_ene );
  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;

  // fillin' the flow
  const unsigned int nFlow = persPart.m_flow.size();
291
  p->m_flow.clear();
292
  for ( unsigned int iFlow= 0; iFlow != nFlow; ++iFlow ) {
293
294
    p->m_flow.set_icode( persPart.m_flow[iFlow].first,
                         persPart.m_flow[iFlow].second );
295
  }
296
#endif
297
298
299
300
301
302
303

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

  return p;
}