McEventCollectionCnv_p4.cxx 28.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_p4.cxx
8
9
// Implementation file for class McEventCollectionCnv_p4
// Author: S.Binet<binet@cern.ch>
10
///////////////////////////////////////////////////////////////////
11
12
13
14
15
16
17
18
19
20
21


// STL includes
#include <utility>
#include <cmath>
#include <cfloat> // for DBL_EPSILON

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

22
#include "GenInterfaces/IHepMCWeightSvc.h"
23

24
#include "McEventCollectionCnv_utils.h"
25

26
///////////////////////////////////////////////////////////////////
27
// Constructors
28
///////////////////////////////////////////////////////////////////
29
30
31
32


McEventCollectionCnv_p4::McEventCollectionCnv_p4() :
  Base_t( ),
33
  m_isPileup(false),m_hepMCWeightSvc("HepMCWeightSvc","McEventCollectionCnv_p4")
34
35
36
37
{}

McEventCollectionCnv_p4::McEventCollectionCnv_p4( const McEventCollectionCnv_p4& rhs ) :
  Base_t( rhs ),
38
  m_isPileup(false),m_hepMCWeightSvc("HepMCWeightSvc","McEventCollectionCnv_p4")
39
40
{}

41
McEventCollectionCnv_p4&
42
43
44
45
46
McEventCollectionCnv_p4::operator=( const McEventCollectionCnv_p4& rhs )
{
  if ( this != &rhs ) {
    Base_t::operator=( rhs );
    m_isPileup=rhs.m_isPileup;
47
48
    m_hepMCWeightSvc = rhs.m_hepMCWeightSvc;
  }
49
50
51
52
53
54
55
56
57
58
  return *this;
}

// Destructor
///////////////

McEventCollectionCnv_p4::~McEventCollectionCnv_p4()
{
}

59
60
///////////////////////////////////////////////////////////////////
// Const methods:
61
62
63
///////////////////////////////////////////////////////////////////

void McEventCollectionCnv_p4::persToTrans( const McEventCollection_p4* persObj,
64
65
                                           McEventCollection* transObj,
                                           MsgStream& msg )
66
67
{
  msg << MSG::DEBUG << "Loading McEventCollection from persistent state..."
68
      << endmsg;
69
70

  // elements are managed by DataPool
71
72
73
  if (!m_isPileup)
    {
      transObj->clear(SG::VIEW_ELEMENTS);
74
    }
75
76
77
78
79
  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 );
80
    }
81
82
83
84
  const unsigned int nParts = persObj->m_genParticles.size();
  if ( datapools.part.capacity() - datapools.part.allocated() < nParts )
    {
      datapools.part.reserve( datapools.part.allocated() + nParts );
85
    }
86
87
88
89
  const unsigned int nEvts = persObj->m_genEvents.size();
  if ( datapools.evt.capacity() - datapools.evt.allocated() < nEvts )
    {
      datapools.evt.reserve( datapools.evt.allocated() + nEvts );
90
91
    }

92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
  transObj->reserve( nEvts );
  for ( std::vector<GenEvent_p4>::const_iterator
          itr = persObj->m_genEvents.begin(),
          itrEnd = persObj->m_genEvents.end();
        itr != itrEnd;
        ++itr )
    {
      const GenEvent_p4& persEvt = *itr;
      HepMC::GenEvent * genEvt(0);
      if(m_isPileup)
        {
          genEvt = new HepMC::GenEvent();
        }
      else
        {
107
          genEvt        =  datapools.getGenEvent();
108
        }
109
110
111
112
113
114
115
116
117
118
119
#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->add_attribute("signal_process_vertex",std::make_shared<HepMC3::IntAttribute>(0));
      genEvt->weights()= persEvt.m_weights;
      genEvt->add_attribute("random_states",std::make_shared<HepMC3::VectorLongIntAttribute>(persEvt.m_randomStates));
      //restore weight names from the dedicated svc (which was keeping them in metadata for efficiency)
120
      if(!genEvt->run_info()) genEvt->set_run_info(std::make_shared<HepMC3::GenRunInfo>());
121
122
123
124
125
126
127
      if(genEvt->run_info()) genEvt->run_info()->set_weight_names(name_index_map_to_names(m_hepMCWeightSvc->weightNames()));


       // pdfinfo restore
      if (!persEvt.m_pdfinfo.empty())
        {
          const std::vector<double>& pdf = persEvt.m_pdfinfo;
128
              HepMC3::GenPdfInfoPtr pi = std::make_shared<HepMC3::GenPdfInfo>();
129
130
131
132
133
134
135
136
137
138
              pi->set(
              static_cast<int>(pdf[6]), // id1
              static_cast<int>(pdf[5]), // id2
              pdf[4],                   // x1
              pdf[3],                   // x2
              pdf[2],                   // scalePDF
              pdf[1],                   // pdf1
              pdf[0] );                 // pdf2
              genEvt->set_pdf_info(pi);
        }
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153

      transObj->push_back( genEvt );

      // create a temporary map associating the barcode of an end-vtx to its
      // 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-persEvt.m_particlesBegin)/2 );
      // 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 );
        } //> end loop over vertices

154
#else
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
      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();
      //restore weight names from the dedicated svc (which was keeping them in metadata for efficiency)
      genEvt->m_weights.m_names = m_hepMCWeightSvc->weightNames();

      // pdfinfo restore
      delete genEvt->m_pdf_info; genEvt->m_pdf_info = 0;
      if (!persEvt.m_pdfinfo.empty())
        {
          const std::vector<double>& pdf = persEvt.m_pdfinfo;
          genEvt->m_pdf_info = new HepMC::PdfInfo
            ( static_cast<int>(pdf[6]), // id1
              static_cast<int>(pdf[5]), // id2
              pdf[4],                   // x1
              pdf[3],                   // x2
              pdf[2],                   // scalePDF
              pdf[1],                   // pdf1
              pdf[0] );                 // pdf2
        }
182

183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199

      transObj->push_back( genEvt );

      // create a temporary map associating the barcode of an end-vtx to its
      // 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-
                                    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,
200
                                               datapools ) );
201
        } //> end loop over vertices
202
#endif
203
204
205
206
      // set the signal process vertex
      const int sigProcVtx = persEvt.m_signalProcessVtx;
      if ( sigProcVtx != 0 )
        {
207
          HepMC::set_signal_process_vertex(genEvt,HepMC::barcode_to_vertex(genEvt, sigProcVtx ) );
208
209
210
211
212
213
214
215
216
        }

      // connect particles to their end vertices
      for ( ParticlesMap_t::iterator
              p = partToEndVtx.begin(),
              endItr = partToEndVtx.end();
            p != endItr;
            ++p )
        {
217
          auto decayVtx= HepMC::barcode_to_vertex(genEvt, p->second );
218
219
220
221
222
223
224
225
          if ( decayVtx )
            {
              decayVtx->add_particle_in( p->first );
            }
          else
            {
              msg << MSG::ERROR
                  << "GenParticle points to null end vertex !!"
226
                  << endmsg;
227
228
229
            }
        }
    } //> end loop over m_genEvents
230
231

  msg << MSG::DEBUG << "Loaded McEventCollection from persistent state [OK]"
232
      << endmsg;
233
234
235
236

  return;
}

237
238
239
void McEventCollectionCnv_p4::transToPers( const McEventCollection* transObj,
                                           McEventCollection_p4* persObj,
                                           MsgStream& msg )
240
241
{
  msg << MSG::DEBUG << "Creating persistent state of McEventCollection..."
242
      << endmsg;
243
244
245
246
247
248
249
250
  persObj->m_genEvents.reserve( transObj->size() );

  const std::pair<unsigned int,unsigned int> stats = nbrParticlesAndVertices( transObj );
  persObj->m_genParticles.reserve( stats.first  );
  persObj->m_genVertices.reserve ( stats.second );

  const McEventCollection::const_iterator itrEnd = transObj->end();
  for ( McEventCollection::const_iterator itr = transObj->begin();
251
252
253
        itr != itrEnd;
        ++itr )
    {
254
255
256
257
258
#ifdef HEPMC3 
      const unsigned int nPersVtx   = persObj->m_genVertices.size();
      const unsigned int nPersParts = persObj->m_genParticles.size();
      const HepMC::GenEvent* genEvt = *itr;
      //save the weight names to metadata via the HepMCWeightSvc
259
      if (genEvt->run_info()) {
260
261
262
263
264
265
266
        if (!genEvt->run_info()->weight_names().empty()) {
          m_hepMCWeightSvc->setWeightNames(  names_to_name_index_map(genEvt->weight_names()) ).ignore();
        } else {
          //AV : This to be decided if one would like to have default names.
          //std::vector<std::string> names{"0"};
          //m_hepMCWeightSvc->setWeightNames( names_to_name_index_map(names) );
        }
267
      }
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
      auto A_signal_process_id=genEvt->attribute<HepMC3::IntAttribute>("signal_process_id");    
      auto A_event_scale=genEvt->attribute<HepMC3::DoubleAttribute>("event_scale");    
      auto A_alphaQCD=genEvt->attribute<HepMC3::DoubleAttribute>("alphaQCD");    
      auto A_alphaQED=genEvt->attribute<HepMC3::DoubleAttribute>("alphaQED");    
      auto A_signal_process_vertex=genEvt->attribute<HepMC3::IntAttribute>("signal_process_vertex");    
      auto A_random_states=genEvt->attribute<HepMC3::VectorLongIntAttribute>("random_states");    

      persObj->m_genEvents.
      push_back( GenEvent_p4( A_signal_process_id?(A_signal_process_id->value()):0,
                                genEvt->event_number(),
                                A_event_scale?(A_event_scale->value()):0.0, 
                                A_alphaQCD?(A_alphaQCD->value()):0.0, 
                                A_alphaQED?(A_alphaQED->value()):0.0, 
                                A_signal_process_vertex?(A_signal_process_vertex->value()):0, 
                                genEvt->weights(), 
                                std::vector<double>(),//No idea why it is empty 
                                A_random_states?(A_random_states->value()):std::vector<long>(), 
                                nPersVtx,
                                nPersVtx + genEvt->vertices().size(),
                                nPersParts,
                                nPersParts + genEvt->particles().size() ) );

      //PdfInfo encoding
   if (genEvt->pdf_info())
        {
          auto pi=genEvt->pdf_info();
          GenEvent_p4& persEvt = persObj->m_genEvents.back();
          std::vector<double>& pdfinfo = persEvt.m_pdfinfo;
          pdfinfo.resize(7);
          pdfinfo[6] = static_cast<double>(pi->parton_id[0]);
          pdfinfo[5] = static_cast<double>(pi->parton_id[1]);
          pdfinfo[4] = pi->x[0];
          pdfinfo[3] = pi->x[1];
          pdfinfo[2] = pi->scale;
          pdfinfo[1] = pi->xf[0];
          pdfinfo[0] = pi->xf[1];
        }
      // create vertices
      for ( auto v: genEvt->vertices())
        {
          writeGenVertex( v, *persObj );
        }

#else
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
      const unsigned int nPersVtx   = persObj->m_genVertices.size();
      const unsigned int nPersParts = persObj->m_genParticles.size();
      const HepMC::GenEvent* genEvt = *itr;
      const int signalProcessVtx = genEvt->m_signal_process_vertex
        ? genEvt->m_signal_process_vertex->barcode()
        : 0;
      //save the weight names to metadata via the HepMCWeightSvc
      m_hepMCWeightSvc->setWeightNames( genEvt->m_weights.m_names ).ignore(); 
      persObj->m_genEvents.
        push_back( GenEvent_p4( genEvt->m_signal_process_id,
                                genEvt->m_event_number,
                                genEvt->m_event_scale,
                                genEvt->m_alphaQCD,
                                genEvt->m_alphaQED,
                                signalProcessVtx,
                                genEvt->m_weights.m_weights,
                                std::vector<double>(),
                                genEvt->m_random_states,
                                nPersVtx,
                                nPersVtx + genEvt->vertices_size(),
                                nPersParts,
                                nPersParts + genEvt->particles_size() ) );
      //PdfInfo encoding
      if (genEvt->m_pdf_info)
        {
          GenEvent_p4& persEvt = persObj->m_genEvents.back();
          std::vector<double>& pdfinfo = persEvt.m_pdfinfo;
          pdfinfo.resize(7);
          pdfinfo[6] = static_cast<double>(genEvt->m_pdf_info->m_id1);
          pdfinfo[5] = static_cast<double>(genEvt->m_pdf_info->m_id2);
          pdfinfo[4] = genEvt->m_pdf_info->m_x1;
          pdfinfo[3] = genEvt->m_pdf_info->m_x2;
          pdfinfo[2] = genEvt->m_pdf_info->m_scalePDF;
          pdfinfo[1] = genEvt->m_pdf_info->m_pdf1;
          pdfinfo[0] = genEvt->m_pdf_info->m_pdf2;
        }

      // create vertices
      const HepMC::GenEvent::vertex_const_iterator endVtx=genEvt->vertices_end();
      for ( HepMC::GenEvent::vertex_const_iterator i = genEvt->vertices_begin();
            i != endVtx;
            ++i )
        {
          writeGenVertex( **i, *persObj );
        }
357
#endif
358
359

    } //> end loop over GenEvents
360
361

  msg << MSG::DEBUG << "Created persistent state of HepMC::GenEvent [OK]"
362
      << endmsg;
363
364
365
366
  return;
}


367
HepMC::GenVertexPtr
368
McEventCollectionCnv_p4::createGenVertex( const McEventCollection_p4& persEvt,
369
370
                                          const GenVertex_p4& persVtx,
                                          ParticlesMap_t& partToEndVtx,
371
                                          HepMC::DataPool& datapools, HepMC::GenEvent* parent ) const
372
{
373
374
375
376
377
378
379
  HepMC::GenVertexPtr vtx(0);
  if(m_isPileup)
    {
      vtx=HepMC::newGenVertexPtr();
    }
  else
    {
380
      vtx = datapools.getGenVertex();
381
    }
382
  if (parent) parent->add_vertex(vtx);
383
#ifdef HEPMC3
384
  vtx->set_position(HepMC::FourVector( persVtx.m_x , persVtx.m_y , persVtx.m_z ,persVtx.m_t ));
385
  vtx->set_status(persVtx.m_id);
386
387
388
  vtx->add_attribute("weights",std::make_shared<HepMC3::VectorFloatAttribute>(persVtx.m_weights));
  vtx->add_attribute("barcode",std::make_shared<HepMC3::IntAttribute>(persVtx.m_barcode));

389
390
391
392
393
394
395
396
397
398
399
400
401
402
  // handle the in-going (orphans) particles
  //Is this needed in 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 );
    }
403
#else
404
405
406
407
408
409
410
411
412
  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(),
413
                                    persVtx.m_weights.end() );
414
415
  vtx->m_event   = 0;
  vtx->m_barcode = persVtx.m_barcode;
416

417
418
  // handle the in-going (orphans) particles
  const unsigned int nPartsIn = persVtx.m_particlesIn.size();
419
420
421
422
423
424
425
  for ( unsigned int i = 0; i != nPartsIn; ++i )
    {
      createGenParticle( persEvt.m_genParticles[persVtx.m_particlesIn[i]],
                         partToEndVtx,
                         datapools );
    }

426
427
  // now handle the out-going particles
  const unsigned int nPartsOut = persVtx.m_particlesOut.size();
428
429
430
431
432
433
  for ( unsigned int i = 0; i != nPartsOut; ++i )
    {
      vtx->add_particle_out( createGenParticle( persEvt.m_genParticles[persVtx.m_particlesOut[i]],
                                                partToEndVtx,
                                                datapools ) );
    }
434
#endif
435
436
437
438

  return vtx;
}

439
HepMC::GenParticlePtr
440
McEventCollectionCnv_p4::createGenParticle( const GenParticle_p4& persPart,
441
                                            ParticlesMap_t& partToEndVtx,
442
                                            HepMC::DataPool& datapools, HepMC::GenVertexPtr parent ) const
443
{
444
445
446
447
448
449
450
  HepMC::GenParticlePtr p(0);
  if (m_isPileup)
    {
      p = HepMC::newGenParticlePtr();
    }
  else
    {
451
      p    = datapools.getGenParticle();
452
    }
453
  if (parent) parent->add_particle_out(p);
454
#ifdef HEPMC3
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
  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));  

  // 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 )
    {
      double 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 );
472
      p->set_momentum( HepMC::FourVector(persPart.m_px,persPart.m_py,persPart.m_pz,temp_e));
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
    }
  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 );
     p->set_momentum( HepMC::FourVector( persPart.m_px,
                         persPart.m_py,
                         persPart.m_pz,
                         signEne * persPart_ene ));
    }

  // setup 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
498
499
500
501
502
503
504
505
506
507
508
509
510
  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.)
511
512
513
514
515
516
  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);
517
      double temp_e = std::sqrt( (long double)(persPart.m_px)*persPart.m_px +
518
519
520
521
522
523
524
525
526
                            (long double)(persPart.m_py)*persPart.m_py +
                            (long double)(persPart.m_pz)*persPart.m_pz +
                            (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 =
527
        std::sqrt( std::abs((long double)(persPart.m_px)*persPart.m_px +
528
529
530
531
532
533
534
535
536
                  (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 );
      p->m_momentum.set( persPart.m_px,
                         persPart.m_py,
                         persPart.m_pz,
                         signEne * persPart_ene );
    }
537
538
539

  // setup flow
  const unsigned int nFlow = persPart.m_flow.size();
540
  p->m_flow.clear();
541
542
543
544
545
  for ( unsigned int iFlow= 0; iFlow != nFlow; ++iFlow )
    {
      p->m_flow.set_icode( persPart.m_flow[iFlow].first,
                           persPart.m_flow[iFlow].second );
    }
546
#endif
547

548
549
550
551
  if ( persPart.m_endVtx != 0 )
    {
      partToEndVtx[p] = persPart.m_endVtx;
    }
552
553
554
555

  return p;
}

556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
#ifdef HEPMC3
void McEventCollectionCnv_p4::writeGenVertex( HepMC::ConstGenVertexPtr vtx,
                                              McEventCollection_p4& persEvt ) const
{                                              
  const HepMC::FourVector& position = vtx->position();
  auto A_weights=vtx->attribute<HepMC3::VectorDoubleAttribute>("weights"); 
  auto A_barcode=vtx->attribute<HepMC3::IntAttribute>("barcode"); 
  std::vector<double> weights=A_weights?(A_weights->value()):std::vector<double>();
  persEvt.m_genVertices.push_back(
                                  GenVertex_p4( position.x(),
                                                position.y(),
                                                position.z(),
                                                position.t(),
                                                vtx->id(),
                                                weights.begin(),
                                                weights.end(),
                                                A_barcode?(A_barcode->value()):vtx->id()
                                                ) );
  GenVertex_p4& persVtx = persEvt.m_genVertices.back();
  // we write only the orphans in-coming particles
  persVtx.m_particlesIn.reserve(vtx->particles_in().size());
  for ( auto p: vtx->particles_in())
    {
      if ( !p->production_vertex() )
        {
          persVtx.m_particlesIn.push_back( writeGenParticle(p, persEvt ));
        }
    }
  persVtx.m_particlesOut.reserve(vtx->particles_out().size());
  for ( auto p: vtx->particles_out())
    {
      persVtx.m_particlesOut.push_back( writeGenParticle(p, persEvt ) );
    }
  return;
}
#else
592
void McEventCollectionCnv_p4::writeGenVertex( const HepMC::GenVertex& vtx,
593
                                              McEventCollection_p4& persEvt ) const
594
595
{
  const HepMC::FourVector& position = vtx.m_position;
596
597
598
599
600
601
602
603
604
  persEvt.m_genVertices.push_back(
                                  GenVertex_p4( position.x(),
                                                position.y(),
                                                position.z(),
                                                position.t(),
                                                vtx.m_id,
                                                vtx.m_weights.m_weights.begin(),
                                                vtx.m_weights.m_weights.end(),
                                                vtx.m_barcode ) );
605
606
607
  GenVertex_p4& persVtx = persEvt.m_genVertices.back();

  // we write only the orphans in-coming particles
608
  const std::vector<HepMC::GenParticlePtr>::const_iterator endInVtx = vtx.m_particles_in.end();
609
  persVtx.m_particlesIn.reserve(vtx.m_particles_in.size());
610
  for ( std::vector<HepMC::GenParticlePtr>::const_iterator p = vtx.m_particles_in.begin();
611
612
613
614
615
616
617
        p != endInVtx;
        ++p )
    {
      if ( 0 == (*p)->production_vertex() )
        {
          persVtx.m_particlesIn.push_back( writeGenParticle( **p, persEvt ) );
        }
618
619
    }

620
  const std::vector<HepMC::GenParticlePtr>::const_iterator endOutVtx = vtx.m_particles_out.end();
621
  persVtx.m_particlesOut.reserve(vtx.m_particles_out.size());
622
  for ( std::vector<HepMC::GenParticlePtr>::const_iterator p = vtx.m_particles_out.begin();
623
624
625
626
627
628
        p != endOutVtx;
        ++p )
    {
      persVtx.m_particlesOut.push_back( writeGenParticle( **p, persEvt ) );
    }

629
630
  return;
}
631
#endif
632

633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
#ifdef HEPMC3
int McEventCollectionCnv_p4::writeGenParticle( HepMC::ConstGenParticlePtr p,
                                               McEventCollection_p4& persEvt ) const
{                                               
  const HepMC::FourVector& mom = p->momentum();
  const double ene = mom.e();
  const double m2  = mom.m2();

  // Definitions of Bool isTimeLilike, isSpacelike and isLightlike according to HepLorentzVector definition
  const bool useP2M2 = !(m2 > 0) &&   // !isTimelike
    (m2 < 0) &&   //  isSpacelike
    !(std::abs(m2) < 2.0*DBL_EPSILON*ene*ene); // !isLightlike

    const short recoMethod = ( !useP2M2 ? 0: ( ene >= 0. ? 1 : 2 ) );
    auto A_theta=p->attribute<HepMC3::DoubleAttribute>("theta"); 
    auto A_phi=p->attribute<HepMC3::DoubleAttribute>("phi"); 
    auto A_flows=p->attribute<HepMC3::VectorIntAttribute>("flows"); 
    
    
    persEvt.m_genParticles.push_back( GenParticle_p4( mom.px(),
                               mom.py(),
                               mom.pz(),
                               mom.m(),
                               p->pdg_id(),
                               p->status(),
                               A_flows?(A_flows->value().size()):0,
                               A_theta?(A_theta->value()):0.0,
                               A_phi?(A_phi->value()):0.0,
                               p->production_vertex()?(HepMC::barcode(p->production_vertex())):0,
                               p->end_vertex()?(HepMC::barcode(p->end_vertex())):0,
                               HepMC::barcode(p),
                               recoMethod ) );
  
  std::vector< std::pair<int,int> > flow_hepmc2;
  if(A_flows) flow_hepmc2=vector_to_vector_int_int(A_flows->value());
  persEvt.m_genParticles.back().m_flow.assign( flow_hepmc2.begin(),flow_hepmc2.end() );
  // we return the index of the particle in the big vector of particles
  // (contained by the persistent GenEvent)
  return (persEvt.m_genParticles.size() - 1);
}                                               
#else
674
int McEventCollectionCnv_p4::writeGenParticle( const HepMC::GenParticle& p,
675
                                               McEventCollection_p4& persEvt ) const
676
677
678
679
680
681
682
{
  const HepMC::FourVector& mom = p.m_momentum;
  const double ene = mom.e();
  const double m2  = mom.m2();

  // Definitions of Bool isTimeLilike, isSpacelike and isLightlike according to HepLorentzVector definition
  const bool useP2M2 = !(m2 > 0) &&   // !isTimelike
683
    (m2 < 0) &&   //  isSpacelike
684
    !(std::abs(m2) < 2.0*DBL_EPSILON*ene*ene); // !isLightlike
685

686
687
688
689
690
691
  const short recoMethod = ( !useP2M2
                             ? 0
                             : ( ene >= 0. //*GeV
                                 ? 1
                                 : 2 ) );

692
693
  persEvt.m_genParticles.
    push_back( GenParticle_p4( mom.px(),
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
                               mom.py(),
                               mom.pz(),
                               mom.m(),
                               p.m_pdg_id,
                               p.m_status,
                               p.m_flow.size(),
                               p.m_polarization.theta(),
                               p.m_polarization.phi(),
                               p.m_production_vertex
                               ? p.m_production_vertex->barcode()
                               : 0,
                               p.m_end_vertex
                               ? p.m_end_vertex->barcode()
                               : 0,
                               p.m_barcode,
                               recoMethod ) );
  persEvt.m_genParticles.back().m_flow.assign( p.m_flow.begin(),
                                               p.m_flow.end() );
712
713
714
715
716

  // we return the index of the particle in the big vector of particles
  // (contained by the persistent GenEvent)
  return (persEvt.m_genParticles.size() - 1);
}
717
#endif
718

719
720
721
void McEventCollectionCnv_p4::setPileup()
{
  m_isPileup = true;
722
}