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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
#include "DumpMCParticleAlg.h"
#include "GaudiKernel/DataObject.h"
#include "GaudiKernel/IHistogramSvc.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "DetInterface/IGeomSvc.h"
#include "DataHelper/HelixClass.h"
#include "DD4hep/Detector.h"
#include "DD4hep/DD4hepUnits.h"
#include "CLHEP/Units/SystemOfUnits.h"
#include <math.h>
DECLARE_COMPONENT( DumpMCParticleAlg )
//------------------------------------------------------------------------------
DumpMCParticleAlg::DumpMCParticleAlg( const std::string& name, ISvcLocator* pSvcLocator )
: Algorithm( name, pSvcLocator ) {
declareProperty("MCParticleCollection", _inMCColHdl, "Handle of the Input MCParticle collection");
m_thisName = name;
}
//------------------------------------------------------------------------------
StatusCode DumpMCParticleAlg::initialize(){
info() << "Booking Ntuple" << endmsg;
NTuplePtr nt1(ntupleSvc(), "MyTuples/MC");
if ( !nt1 ) {
m_tuple = ntupleSvc()->book("MyTuples/MC",CLID_ColumnWiseTuple,"MC truth");
if ( 0 != m_tuple ) {
m_tuple->addItem ("nmc", m_nParticles, 0, 1000 ).ignore();
m_tuple->addIndexedItem ("pdg", m_nParticles, m_pdgID ).ignore();
m_tuple->addIndexedItem ("genStatus", m_nParticles, m_genStatus ).ignore();
m_tuple->addIndexedItem ("simStatus", m_nParticles, m_simStatus ).ignore();
m_tuple->addIndexedItem ("charge", m_nParticles, m_charge ).ignore();
m_tuple->addIndexedItem ("time", m_nParticles, m_time ).ignore();
m_tuple->addIndexedItem ("mass", m_nParticles, m_mass ).ignore();
m_tuple->addIndexedItem ("vx", m_nParticles, m_vx ).ignore();
m_tuple->addIndexedItem ("vy", m_nParticles, m_vy ).ignore();
m_tuple->addIndexedItem ("vz", m_nParticles, m_vz ).ignore();
m_tuple->addIndexedItem ("px", m_nParticles, m_px ).ignore();
m_tuple->addIndexedItem ("py", m_nParticles, m_py ).ignore();
m_tuple->addIndexedItem ("pz", m_nParticles, m_pz ).ignore();
m_tuple->addIndexedItem ("d0", m_nParticles, m_d0 ).ignore();
m_tuple->addIndexedItem ("phi0", m_nParticles, m_phi0 ).ignore();
m_tuple->addIndexedItem ("omega", m_nParticles, m_omega ).ignore();
m_tuple->addIndexedItem ("z0", m_nParticles, m_z0 ).ignore();
m_tuple->addIndexedItem ("tanLambda", m_nParticles, m_tanLambda ).ignore();
}
else { // did not manage to book the N tuple....
fatal() << "Cannot bool MyTuples/MC " << endmsg;
return StatusCode::FAILURE;
}
}
else{
m_tuple = nt1;
}
auto geomSvc = service<IGeomSvc>("GeomSvc");
if(geomSvc){
const dd4hep::Direction& field = geomSvc->lcdd()->field().magneticField(dd4hep::Position(0,0,0));
m_field = field.z()/dd4hep::tesla;
info() << "Magnetic field will obtain from GeomSvc = " << m_field << " tesla" << endmsg;
}
else{
info() << "Failed to find GeomSvc ..." << endmsg;
info() << "Magnetic field will use what input through python option for this algorithm namse as Field, now " << m_field << " tesla" << endmsg;
}
_nEvt = 0;
return StatusCode::SUCCESS;
}
//------------------------------------------------------------------------------
StatusCode DumpMCParticleAlg::execute(){
const edm4hep::MCParticleCollection* mcCols = nullptr;
try {
mcCols = _inMCColHdl.get();
}
catch ( GaudiException &e ) {
debug() << "Collection " << _inMCColHdl.fullKey() << " is unavailable in event " << _nEvt << endmsg;
}
if(mcCols){
m_nParticles = 0;
for(auto particle : *mcCols){
m_pdgID[m_nParticles] = particle.getPDG();
m_genStatus[m_nParticles] = particle.getGeneratorStatus();
m_simStatus[m_nParticles] = particle.getSimulatorStatus();
m_charge[m_nParticles] = particle.getCharge();
m_time[m_nParticles] = particle.getTime();
m_mass[m_nParticles] = particle.getMass();
const auto& vertex = particle.getVertex();
m_vx[m_nParticles] = vertex.x;
m_vy[m_nParticles] = vertex.y;
m_vz[m_nParticles] = vertex.z;
const auto& momentum = particle.getMomentum();
101
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
m_px[m_nParticles] = momentum.x;
m_py[m_nParticles] = momentum.y;
m_pz[m_nParticles] = momentum.z;
HelixClass helix;
float posV[3] = {vertex.x,vertex.y,vertex.z};
float momV[3] = {momentum.x,momentum.y,momentum.z};
helix.Initialize_VP(posV,momV,particle.getCharge(),m_field);
float phiMC = helix.getPhi0();
if(phiMC>CLHEP::pi) phiMC = phiMC - CLHEP::twopi;
m_phi0[m_nParticles] = phiMC;
m_d0[m_nParticles] = helix.getD0();
m_omega[m_nParticles] = helix.getOmega();
m_z0[m_nParticles] = helix.getZ0();
m_tanLambda[m_nParticles] = helix.getTanLambda();
m_nParticles++;
}
debug() << "MCParticle: " << m_nParticles <<endmsg;
}
m_tuple->write();
_nEvt++;
return StatusCode::SUCCESS;
}
//------------------------------------------------------------------------------
StatusCode DumpMCParticleAlg::finalize(){
debug() << "Finalizing..." << endmsg;
return StatusCode::SUCCESS;
}