ATLAS Offline Software
PhysicsValidationTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // PhysicsValidationTool.cxx, (c) ATLAS Detector software
8 
9 // class header
10 #include "PhysicsValidationTool.h"
11 
12 // Tracking
13 #include "TrkParameters/TrackParameters.h"
14 // ISF includes
15 #include "ISF_Event/ISFParticle.h"
16 #include "ISF_Event/ParticleUserInformation.h"
17 // CLHEP
18 #include "CLHEP/Units/SystemOfUnits.h"
19 #include "CLHEP/Units/PhysicalConstants.h"
20 //
21 // ROOT includes
22 #include "TTree.h"
23 
24 #include <cmath>
25 
26 #include <iostream>
27 
28 /*=========================================================================
29  * DESCRIPTION OF FUNCTION:
30  * ==> see headerfile
31  *=======================================================================*/
32 iFatras::PhysicsValidationTool::PhysicsValidationTool( const std::string& t,
33  const std::string& n,
34  const IInterface* p )
35  : base_class(t,n,p),
36  m_thistSvc("THistSvc",n),
37  m_validationStream("ISFFatras"),
38  m_eloss(0.),
39  m_ionloss(0.),
40  m_radloss(0.),
41  m_zOaTr(0.),
42  m_wZ(0.),
43  m_thIn(0.),
44  m_phIn(0.),
45  m_dIn(0.),
46  m_thEnd(0.),
47  m_phEnd(0.),
48  m_dEnd(0.),
49  m_X0(0.),
50  m_dt(0.),
51  m_interactions(nullptr),
52  m_process(0),
53  m_pdg_mother(0),
54  m_gen_mother(0),
55  m_nChild(0),
56  m_vtx_dist(0.),
57  m_vtx_theta(0.),
58  m_vtx_phi(0.),
59  m_vtx_e_diff(0.),
60  m_vtx_p_diff(0.),
61  m_vtx_plong_diff(0.),
62  m_vtx_pperp_diff(0.),
63  m_p_mother(0.)
64 {
65  // validation output section
66  declareProperty("ValidationStreamName",
67  m_validationStream = "ISFFatras",
68  "Name of the output stream" );
69  declareProperty("THistService",
70  m_thistSvc,
71  "The THistSvc" );
72 }
73 
74 /*=========================================================================
75  * DESCRIPTION OF FUNCTION:
76  * ==> see headerfile
77  *=======================================================================*/
78 iFatras::PhysicsValidationTool::~PhysicsValidationTool()
79 {
80  m_interactions = nullptr;
81  delete(m_interactions);
82 }
83 
84 /*=========================================================================
85  * DESCRIPTION OF FUNCTION:
86  * ==> see headerfile
87  *=======================================================================*/
88 StatusCode
89 iFatras::PhysicsValidationTool::initialize()
90 {
91 
92  ATH_MSG_VERBOSE( "initialize()" );
93 
94  // retrieve the histogram service
95  if ( m_thistSvc.retrieve().isSuccess() ) {
96  // Create the prefix of histogram names for the THistSvc
97  const char *treeName="particles";
98  const std::string prefix = "/" + m_validationStream + "/"+ treeName;
99  m_particles = new TTree( treeName, treeName );
100  m_particles->Branch("pdg" , &m_pdg , "pdg/I" );
101  m_particles->Branch("scIn" , &m_scIn , "scIn/I" );
102  m_particles->Branch("scOut" , &m_scEnd , "scOut/I" );
103  m_particles->Branch("gen" , &m_gen , "gen/I" );
104  m_particles->Branch("geoID" , &m_geoID , "geoID/I" );
105  m_particles->Branch("pth" , &m_pth , "pth/F" );
106  m_particles->Branch("pph" , &m_pph , "pph/F" );
107  m_particles->Branch("p" , &m_p , "p/F" );
108  m_particles->Branch("eloss" , &m_eloss , "eloss/F" );
109  m_particles->Branch("ionloss" , &m_ionloss , "ionloss/F" );
110  m_particles->Branch("radloss" , &m_radloss , "radloss/F" );
111  m_particles->Branch("zOaTr" , &m_zOaTr , "zOaTr/F" );
112  m_particles->Branch("wZ" , &m_wZ , "wZ/F" );
113  m_particles->Branch("X0" , &m_X0 , "X0/F" );
114  m_particles->Branch("dt" , &m_dt , "dt/F" );
115  m_particles->Branch("thIn" , &m_thIn , "thIn/F" );
116  m_particles->Branch("phIn" , &m_phIn , "phIn/F" );
117  m_particles->Branch("dIn" , &m_dIn , "dIn/F" );
118  m_particles->Branch("thEnd" , &m_thEnd , "thEnd/F" );
119  m_particles->Branch("phEnd" , &m_phEnd , "phEnd/F" );
120  m_particles->Branch("dEnd" , &m_dEnd , "dEnd/F" );
121 
122  // register the Tree to the THistSvc and return it's StatusCode
123  ATH_CHECK(m_thistSvc->regTree( prefix, m_particles) );
124 
125 
126  const char *treeNameInt="interactions";
127  const std::string prefixInt = "/" + m_validationStream + "/"+ treeNameInt;
128  m_interactions = new TTree( treeNameInt, treeNameInt );
129  m_interactions->Branch("process" , &m_process , "process/I" );
130  m_interactions->Branch("pdg_mother" , &m_pdg_mother , "pdg_mother/I" );
131  m_interactions->Branch("gen_mother" , &m_gen_mother , "gen_mother/I" );
132  m_interactions->Branch("nChild" , &m_nChild , "nch/I" );
133  m_interactions->Branch("vtx_dist" , &m_vtx_dist , "vtx_dist/F" );
134  m_interactions->Branch("vtx_theta" , &m_vtx_theta , "vtx_theta/F" );
135  m_interactions->Branch("vtx_phi" , &m_vtx_phi , "vtx_phi/F" );
136  m_interactions->Branch("vtx_e_diff" , &m_vtx_e_diff , "vtx_e_diff/F" );
137  m_interactions->Branch("vtx_p_diff" , &m_vtx_p_diff , "vtx_p_diff/F" );
138  m_interactions->Branch("vtx_plong_diff" , &m_vtx_plong_diff , "vtx_plong_diff/F" );
139  m_interactions->Branch("vtx_pperp_diff" , &m_vtx_pperp_diff , "vtx_pperp_diff/F" );
140  m_interactions->Branch("p_mother" , &m_p_mother , "p_mother/F" );
141  m_interactions->Branch("pdg_child" , m_pdg_child , "pdg_child[nch]/I" );
142  m_interactions->Branch("fp_child" , m_fp_child , "fp_child[nch]/F" );
143  m_interactions->Branch("oa_child" , m_oa_child , "oa_child[nch]/F" );
144  // register the Tree to the THistSvc and return it's StatusCode
145  ATH_CHECK(m_thistSvc->regTree( prefixInt, m_interactions) );
146 
147  }
148 
149  return StatusCode::SUCCESS;
150 }
151 
152 /*=========================================================================
153  * DESCRIPTION OF FUNCTION:
154  * ==> see headerfile
155  *=======================================================================*/
156 
158 void iFatras::PhysicsValidationTool::saveISFParticleInfo(const ISF::ISFParticle& isp,
159  const Trk::ExtrapolationCell<Trk::TrackParameters>& ec,
160  Trk::ExtrapolationCode ecode) const {
161 
162  m_p = isp.momentum().mag();
163  m_ionloss = ec.eLoss ? ec.eLoss->meanIoni() : 0.;
164  m_radloss = ec.eLoss ? ec.eLoss->meanRad() : 0.; // average estimate
165  m_zOaTr = ec.zOaTrX;
166  m_wZ = ec.zX;
167 
168  if (ecode == Trk::ExtrapolationCode::SuccessBoundaryReached) { // surviving particle
169  m_scEnd = 0;
170  m_eloss = ec.endParameters->momentum().mag()-m_p;
171  //m_dt = time - isp.timeStamp();
172  m_dt = ec.pathLength/CLHEP::c_light; // TODO get timing info from the cell
173 
174  m_thEnd = ec.endParameters->position().theta();
175  m_phEnd = ec.endParameters->position().phi();
176  m_dEnd = ec.endParameters->position().mag();
177 
178  } else {
179 
180  m_scEnd = -1; // undefined
181  if (ecode == Trk::ExtrapolationCode::SuccessPathLimit) m_scEnd = 201; // decay
182  else if (ecode == Trk::ExtrapolationCode::SuccessMaterialLimit ) m_scEnd = ec.materialProcess; // material interaction
183 
184  m_eloss = -m_p;
185  m_thEnd = -10.; // TODO : retrieve position of stopped particle
186  m_phEnd = -10.; // TODO : retrieve position of stopped particle
187  m_dEnd = -1.; // TODO : retrieve position of stopped particle
188  }
189  m_X0 = ec.materialX0;
190 
191  saveInfo(isp);
192 
193 }
194 
196 void iFatras::PhysicsValidationTool::saveISFParticleInfo(const ISF::ISFParticle& isp,
197  const Trk::ExtrapolationCell<Trk::NeutralParameters>& ec,
198  Trk::ExtrapolationCode ecode) const {
199 
200  m_p = isp.momentum().mag();
201  m_ionloss = 0.;
202  m_radloss = 0.;
203  m_zOaTr = ec.zOaTrX;
204  m_wZ = ec.zX;
205 
206  if (ecode == Trk::ExtrapolationCode::SuccessBoundaryReached) { // surviving particle
207  m_scEnd = 0;
208  m_eloss = ec.endParameters->momentum().mag()-m_p;
209  //m_dt = time - isp.timeStamp();
210  m_dt = ec.pathLength/CLHEP::c_light; // TODO get timing info from the cell
211 
212  m_thEnd = ec.endParameters->position().theta();
213  m_phEnd = ec.endParameters->position().phi();
214  m_dEnd = ec.endParameters->position().mag();
215 
216  } else {
217 
218  m_scEnd = -1; // undefined
219  if (ecode == Trk::ExtrapolationCode::SuccessPathLimit) m_scEnd = 201; // decay
220  else if (ecode == Trk::ExtrapolationCode::SuccessMaterialLimit ) m_scEnd = ec.materialProcess; // material interaction
221 
222  m_eloss = -m_p;
223  m_thEnd = -10.; // TODO : retrieve position of stopped particle
224  m_phEnd = -10.; // TODO : retrieve position of stopped particle
225  m_dEnd = -1.; // TODO : retrieve position of stopped particle
226  }
227  m_X0 = ec.materialX0;
228 
229  saveInfo(isp);
230 
231 }
232 
234 void iFatras::PhysicsValidationTool::saveISFParticleInfo(const ISF::ISFParticle& isp,
235  int endProcess,
236  const Trk::TrackParameters* ePar,
237  double time, double dX0 ) const {
238 
239  m_pdg = isp.pdgCode();
240  m_scIn = isp.getUserInformation()? isp.getUserInformation()->process() : 0; // assume primary track
241  m_gen = isp.getUserInformation()? isp.getUserInformation()->generation() : 0; // assume primary track
242  m_pth=isp.momentum().theta();
243  m_pph=isp.momentum().phi();
244  m_geoID = isp.nextGeoID();
245  m_dIn = isp.position().mag();
246  m_thIn = isp.position().theta();
247  m_phIn = isp.position().phi();
248  m_p = isp.momentum().mag();
249  m_scEnd = endProcess;
250  m_eloss = ePar? ePar->momentum().mag()-m_p : -m_p;
251  m_dt = time - isp.timeStamp();
252 
253  m_ionloss = 0.; // n.a.
254  m_radloss = 0.; // n.a.
255  m_zOaTr = 0.; // n.a.
256 
257  if (ePar) {
258 
259  //error propagation
260  //const Trk::MeasuredTrackParameters* mPar = dynamic_cast<const Trk::MeasuredTrackParameters*> (ePar);
261  //if (mPar) std::cout << "error propagation :"<< mPar->localErrorMatrix().error(Trk::qOverP) << std::endl;
262 
263  m_thEnd = ePar->position().theta();
264  m_phEnd = ePar->position().phi();
265  m_dEnd = ePar->position().mag();
266  } else {
267  m_thEnd = m_thIn;
268  m_phEnd = m_phIn;
269  m_dEnd = -1.; // particle stopped and position unknown
270  }
271  m_X0 = dX0;
272 
273  m_particles->Fill();
274 
275 }
276 
277 void iFatras::PhysicsValidationTool::saveInfo(const ISF::ISFParticle& isp) const {
278 
279  // ISF particle info ( common )
280 
281  m_pdg = isp.pdgCode();
282  m_scIn = isp.getUserInformation()? isp.getUserInformation()->process() : 0; // assume primary track
283  m_gen = isp.getUserInformation()? isp.getUserInformation()->generation() : 0; // assume primary track
284  m_pth=isp.momentum().theta();
285  m_pph=isp.momentum().phi();
286  m_geoID = isp.nextGeoID();
287  m_dIn = isp.position().mag();
288  m_thIn = isp.position().theta();
289  m_phIn = isp.position().phi();
290 
291  m_particles->Fill();
292 
293 }
294 
295 void iFatras::PhysicsValidationTool::saveISFVertexInfo(int process,Amg::Vector3D vertex,const ISF::ISFParticle& parent,Amg::Vector3D primIn,
296  Amg::Vector3D* primOut, const ISF::ISFParticleVector& children) const {
297 
298  m_process = process;
299  unsigned int nSec = children.size();
300 
301  int iPrimSurv = primOut ? 1 : 0;
302 
303  m_pdg_mother = parent.pdgCode();
304  m_gen_mother = parent.getUserInformation()? parent.getUserInformation()->generation() : 0; // assume primary track
305 
306  m_p_mother = primIn.mag();
307 
308  m_vtx_dist = vertex.mag();
309  m_vtx_theta = vertex.theta();
310  m_vtx_phi = vertex.phi();
311 
312  m_nChild = nSec+iPrimSurv;
313 
314  Amg::Vector3D pbal(primIn);
315 
316  if (iPrimSurv>0) {
317  m_pdg_child[0] = m_pdg_mother;
318  m_fp_child[0] = primOut->mag()/m_p_mother;
319  m_oa_child[0] = primOut->unit().dot( primIn.unit() );
320  pbal -= *primOut;
321  }
322 
323  unsigned int nSecMax = nSec+iPrimSurv> MAXCHILDREN ? MAXCHILDREN-iPrimSurv : nSec;
324  for (unsigned int isec=0; isec< nSec; isec++) {
325  if (isec<nSecMax) {
326  m_pdg_child[isec+iPrimSurv] = children[isec]->pdgCode();
327  m_fp_child[isec+iPrimSurv] = children[isec]->momentum().mag()/m_p_mother;
328  m_oa_child[isec+iPrimSurv] = children[isec]->momentum().unit().dot( primIn.unit() );
329  }
330  pbal -= children[isec]->momentum();
331  }
332 
333  m_vtx_p_diff = pbal.mag();
334  m_vtx_plong_diff = pbal.dot(primIn)/m_p_mother;
335  m_vtx_pperp_diff = std::sqrt(m_vtx_p_diff*m_vtx_p_diff-m_vtx_plong_diff*m_vtx_plong_diff);
336 
337  m_interactions->Fill();
338 }
iFatras::PhysicsValidationTool::initialize
virtual StatusCode initialize() override
AlgTool initialize method.
Definition: PhysicsValidationTool.cxx:89
TrackParameters.h
iFatras::PhysicsValidationTool::m_thistSvc
ServiceHandle< ITHistSvc > m_thistSvc
the histogram service
Definition: PhysicsValidationTool.h:94
Trk::ParametersBase::position
const Amg::Vector3D & position() const
Access method for the position.
iFatras::PhysicsValidationTool::saveISFVertexInfo
virtual void saveISFVertexInfo(int process, Amg::Vector3D vertex, const ISF::ISFParticle &isp, Amg::Vector3D primIn, Amg::Vector3D *primOut, const ISF::ISFParticleVector &children) const override
Definition: PhysicsValidationTool.cxx:295
ISF::ISFParticle
Definition: ISFParticle.h:42
iFatras::PhysicsValidationTool::saveISFParticleInfo
virtual void saveISFParticleInfo(const ISF::ISFParticle &isp, int endProcess, const Trk::TrackParameters *ePar, double time, double dX0) const override
ISFParticle info: old transport tool.
Definition: PhysicsValidationTool.cxx:234
ISF::ISFParticle::pdgCode
int pdgCode() const
PDG value.
Trk::ExtrapolationCode::SuccessBoundaryReached
@ SuccessBoundaryReached
Definition: ExtrapolationCell.h:113
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
SUSY_SimplifiedModel_PostInclude.process
string process
Definition: SUSY_SimplifiedModel_PostInclude.py:42
Trk::ExtrapolationCell::pathLength
double pathLength
the path length accumulated
Definition: ExtrapolationCell.h:269
Trk::ExtrapolationCell::zX
float zX
z*dInX0 (for average calculations)
Definition: ExtrapolationCell.h:300
Trk::ExtrapolationCell::eLoss
EnergyLoss * eLoss
cumulated energy loss
Definition: ExtrapolationCell.h:298
ISF::ISFParticle::position
const Amg::Vector3D & position() const
The current position of the ISFParticle.
iFatras::PhysicsValidationTool::saveInfo
void saveInfo(const ISF::ISFParticle &isp) const
Definition: PhysicsValidationTool.cxx:277
Trk::ExtrapolationCell::materialX0
double materialX0
collected material so far in units of X0
Definition: ExtrapolationCell.h:272
ISFParticle.h
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
ISF::ISFParticle::nextGeoID
AtlasDetDescr::AtlasRegion nextGeoID() const
next geoID the particle will be simulated in
iFatras::PhysicsValidationTool::PhysicsValidationTool
PhysicsValidationTool(const std::string &, const std::string &, const IInterface *)
Constructor.
Definition: PhysicsValidationTool.cxx:32
beamspotman.n
n
Definition: beamspotman.py:731
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
iFatras::PhysicsValidationTool::~PhysicsValidationTool
virtual ~PhysicsValidationTool()
Destructor.
Definition: PhysicsValidationTool.cxx:78
ISF::ISFParticleVector
std::vector< ISF::ISFParticle * > ISFParticleVector
ISFParticle vector.
Definition: ISFParticleContainer.h:26
Trk::ExtrapolationCode
Definition: ExtrapolationCell.h:105
Trk::ExtrapolationCell::endParameters
T * endParameters
by pointer - are newly created and can be optionally 0
Definition: ExtrapolationCell.h:238
checkCorrelInHIST.prefix
dictionary prefix
Definition: checkCorrelInHIST.py:391
test_pyathena.parent
parent
Definition: test_pyathena.py:15
dumpFileToPlots.treeName
string treeName
Definition: dumpFileToPlots.py:20
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Trk::ExtrapolationCode::SuccessPathLimit
@ SuccessPathLimit
Definition: ExtrapolationCell.h:115
Trk::ParametersBase
Definition: ParametersBase.h:55
Trk::ExtrapolationCell::materialProcess
int materialProcess
the material process to be generated
Definition: ExtrapolationCell.h:276
Trk::ExtrapolationCode::SuccessMaterialLimit
@ SuccessMaterialLimit
Definition: ExtrapolationCell.h:117
ParticleUserInformation.h
ISF::ISFParticle::timeStamp
double timeStamp() const
Timestamp of the ISFParticle.
PhysicsValidationTool.h
ISF::ISFParticle::momentum
const Amg::Vector3D & momentum() const
The current momentum vector of the ISFParticle.
python.PhysicalConstants.c_light
float c_light
Definition: PhysicalConstants.py:63
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
iFatras::PhysicsValidationTool::m_validationStream
std::string m_validationStream
validation THist stream name
Definition: PhysicsValidationTool.h:95
Trk::ExtrapolationCell::zOaTrX
float zOaTrX
z/A*rho*dInX0 (for average calculations)
Definition: ExtrapolationCell.h:299
Trk::ParametersBase::momentum
const Amg::Vector3D & momentum() const
Access method for the momentum.
Trk::vertex
@ vertex
Definition: MeasurementType.h:21
Trk::ExtrapolationCell
Definition: ExtrapolationCell.h:231
python.DecayParser.children
children
Definition: DecayParser.py:32
ISF::ParticleUserInformation::process
int process() const
Definition: ParticleUserInformation.h:84
ISF::ParticleUserInformation::generation
int generation() const
Definition: ParticleUserInformation.h:86
ISF::ISFParticle::getUserInformation
const ParticleUserInformation * getUserInformation() const
get/set ParticleUserInformation
MAXCHILDREN
#define MAXCHILDREN
Definition: PhysicsValidationTool.h:24
Generated on Thu Nov 7 2024 21:23:31 for ATLAS Offline Software by   doxygen 1.8.18