ATLAS Offline Software
Loading...
Searching...
No Matches
DumpObjects.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3*/
4
5#include "DumpObjects.h"
16
17#include "HepPDT/ParticleDataTable.hh"
18
24#include "TrkTrack/TrackInfo.h"
25
26#include "GaudiKernel/ITHistSvc.h"
27#include "TTree.h"
28
31
32#include <fstream>
33
34int InDet::compute_overlap_SP_flag(const int& eta_module_cl1,const int& phi_module_cl1,
35 const int& eta_module_cl2,const int& phi_module_cl2){
36 int flag=-999;
37
38 if( (eta_module_cl1==eta_module_cl2) && (phi_module_cl1==phi_module_cl2) ){
39 flag=0; // not an overlap Space Point
40 }
41 else if((eta_module_cl1!=eta_module_cl2) && (phi_module_cl1==phi_module_cl2) ){
42 flag=1; // overlap Space Point in eta only
43 }
44 else if((eta_module_cl1==eta_module_cl2) && (phi_module_cl1!=phi_module_cl2) ){
45 flag=2; // overlap Space Point in phi only
46 }
47 else{
48 flag=3; // "overlap" Space Point in eta and phi (not sure we can call it overlap)
49 }
50 return flag;
51}
52
53//-------------------------------------------------------------------------
54InDet::DumpObjects::DumpObjects(const std::string &name, ISvcLocator *pSvcLocator)
55 //-------------------------------------------------------------------------
56 : AthAlgorithm(name, pSvcLocator),
57 m_particlePropSvc("PartPropSvc", name) {
58 declareProperty("Offset", m_offset);
59 declareProperty("FileName", m_name = "");
60 //
61 declareProperty("NtupleFileName", m_ntupleFileName);
62 declareProperty("NtupleDirectoryName", m_ntupleDirName);
63 declareProperty("NtupleTreeName", m_ntupleTreeName);
64 declareProperty("maxCL", m_maxCL = 1500000);
65 declareProperty("maxPart", m_maxPart = 1500000);
66 declareProperty("maxSP", m_maxSP = 1500000);
67 declareProperty("maxTRK", m_maxTRK = 1500000);
68 declareProperty("maxDTT", m_maxDTT = 1500000);
69
70 declareProperty("rootFile", m_rootFile);
71}
72
73//-------------------------------------------
75//-------------------------------------------
77
78 // ReadHandle keys
79 ATH_CHECK(m_eventInfoKey.initialize());
81 ATH_CHECK(m_stripClusterKey.initialize());
82 ATH_CHECK(m_pixelClusterKey.initialize());
83 ATH_CHECK(m_pixelSDOKey.initialize());
84 ATH_CHECK(m_stripSDOKey.initialize());
88
89 ATH_CHECK(m_tracksKey.initialize());
90 ATH_CHECK(m_tracksTruthKey.initialize());
92
93
94
95
96 // Grab PixelID helper
97 ATH_CHECK (detStore()->retrieve(m_pixelID, "PixelID") );
98
100 detStore()->retrieve(m_pixelManager, "Pixel").isFailure()) {
101 // if Pixel retrieval fails, try ITkPixel
103 detStore()->retrieve(m_pixelManager, "ITkPixel").isFailure()) {
104 return StatusCode::FAILURE;
105 }
106 }
107
108 // Grab SCT_ID helper
109 ATH_CHECK (detStore()->retrieve(m_SCT_ID,"SCT_ID") );
110
112 detStore()->retrieve(m_SCT_Manager, "SCT").isFailure()) {
113 // if SCT retrieval fails, try ITkStrip
115 detStore()->retrieve(m_SCT_Manager, "ITkStrip").isFailure()) {
116 return StatusCode::FAILURE;
117 }
118 }
119
120 // particle property service
121 ATH_CHECK (m_particlePropSvc.retrieve());
122
123 // and the particle data table
125 if (m_particleDataTable == 0) {
126 ATH_MSG_ERROR("Could not get ParticleDataTable! Cannot associate pdg code with charge. Aborting. ");
127 return StatusCode::FAILURE;
128 }
129
130 // Define the TTree
131 //
132 SmartIF<ITHistSvc> tHistSvc{Gaudi::svcLocator()->service("THistSvc")};
133 ATH_CHECK(tHistSvc.isValid());
134 m_nt = new TTree(TString(m_ntupleTreeName), "Athena Dump for GNN4ITk");
135 // NB: we must not delete the tree, this is done by THistSvc
136 std::string fullNtupleName = m_ntupleFileName + m_ntupleDirName + m_ntupleTreeName;
137 StatusCode sc = tHistSvc->regTree(fullNtupleName, m_nt);
138 if (sc.isFailure()) {
139 ATH_MSG_ERROR("Unable to register TTree: " << fullNtupleName);
140 return sc;
141 }
142
143 if (m_rootFile) {
144 m_SEID = new int[m_maxCL];
145
146 m_CLindex = new int[m_maxCL];
147 m_CLhardware = new std::vector<std::string>;
148 m_CLx = new double[m_maxCL];
149 m_CLy = new double[m_maxCL];
150 m_CLz = new double[m_maxCL];
151 m_CLbarrel_endcap = new int[m_maxCL];
152 m_CLlayer_disk = new int[m_maxCL];
153 m_CLeta_module = new int[m_maxCL];
154 m_CLphi_module = new int[m_maxCL];
155 m_CLside = new int[m_maxCL];
156 m_CLmoduleID = new uint64_t[m_maxCL];
157 m_CLparticleLink_eventIndex = new std::vector<std::vector<int>>;
158 m_CLparticleLink_barcode = new std::vector<std::vector<int>>;
159 m_CLbarcodesLinked = new std::vector<std::vector<bool>>;
160 m_CLparticle_charge = new std::vector<std::vector<float>>;
161 m_CLphis = new std::vector<std::vector<int>>;
162 m_CLetas = new std::vector<std::vector<int>>;
163 m_CLtots = new std::vector<std::vector<int>>;
164 m_CLloc_direction1 = new double[m_maxCL];
165 m_CLloc_direction2 = new double[m_maxCL];
166 m_CLloc_direction3 = new double[m_maxCL];
167 m_CLJan_loc_direction1 = new double[m_maxCL];
168 m_CLJan_loc_direction2 = new double[m_maxCL];
169 m_CLJan_loc_direction3 = new double[m_maxCL];
170 m_CLpixel_count = new int[m_maxCL];
171 m_CLcharge_count = new float[m_maxCL];
172 m_CLloc_eta = new float[m_maxCL];
173 m_CLloc_phi = new float[m_maxCL];
174 m_CLglob_eta = new float[m_maxCL];
175 m_CLglob_phi = new float[m_maxCL];
176 m_CLeta_angle = new double[m_maxCL];
177 m_CLphi_angle = new double[m_maxCL];
178 m_CLnorm_x = new float[m_maxCL];
179 m_CLnorm_y = new float[m_maxCL];
180 m_CLnorm_z = new float[m_maxCL];
181 m_CLlocal_cov = new std::vector<std::vector<double>>;
182
184 m_Part_barcode = new int[m_maxPart];
185 m_Part_px = new float[m_maxPart];
186 m_Part_py = new float[m_maxPart];
187 m_Part_pz = new float[m_maxPart];
188 m_Part_pt = new float[m_maxPart];
189 m_Part_eta = new float[m_maxPart];
190 m_Part_vx = new float[m_maxPart];
191 m_Part_vy = new float[m_maxPart];
192 m_Part_vz = new float[m_maxPart];
193 m_Part_radius = new float[m_maxPart];
194 m_Part_status = new float[m_maxPart];
195 m_Part_charge = new float[m_maxPart];
196 m_Part_pdg_id = new int[m_maxPart];
197 m_Part_passed = new int[m_maxPart];
198 m_Part_vProdNin = new int[m_maxPart];
199 m_Part_vProdNout = new int[m_maxPart];
200 m_Part_vProdStatus = new int[m_maxPart];
202 m_Part_vParentID = new std::vector<std::vector<int>>;
203 m_Part_vParentBarcode = new std::vector<std::vector<int>>;
204
205 m_SPindex = new int[m_maxSP];
206 m_SPx = new double[m_maxSP];
207 m_SPy = new double[m_maxSP];
208 m_SPz = new double[m_maxSP];
209 m_SPCL1_index = new int[m_maxSP];
210 m_SPCL2_index = new int[m_maxSP];
211 m_SPisOverlap = new int[m_maxSP];
212
213 m_SPradius = new double[m_maxSP];
214 m_SPcovr = new double[m_maxSP];
215 m_SPcovz = new double[m_maxSP];
216 m_SPhl_topstrip = new float[m_maxSP];
217 m_SPhl_botstrip = new float[m_maxSP];
218 m_SPtopStripDirection = new std::vector<std::vector<float>>;
219 m_SPbottomStripDirection = new std::vector<std::vector<float>>;
220 m_SPstripCenterDistance = new std::vector<std::vector<float>>;
221 m_SPtopStripCenterPosition = new std::vector<std::vector<float>>;
222
223 m_TRKindex = new int[m_maxTRK];
224 m_TRKtrack_fitter = new int[m_maxTRK];
226 m_TRKproperties = new std::vector<std::vector<int>>;
227 m_TRKpattern = new std::vector<std::vector<int>>;
228 m_TRKndof = new int[m_maxTRK];
229 m_TRKmot = new int[m_maxTRK];
230 m_TRKoot = new int[m_maxTRK];
231 m_TRKchiSq = new float[m_maxTRK];
232 m_TRKmeasurementsOnTrack_pixcl_sctcl_index = new std::vector<std::vector<int>>;
233 m_TRKoutliersOnTrack_pixcl_sctcl_index = new std::vector<std::vector<int>>;
234 m_TRKcharge = new int[m_maxTRK];
235 m_TRKperigee_position = new std::vector<std::vector<double>>;
236 m_TRKperigee_momentum = new std::vector<std::vector<double>>;
237 m_TTCindex = new int[m_maxTRK];
238 m_TTCevent_index = new int[m_maxTRK];
239 m_TTCparticle_link = new int[m_maxTRK];
240 m_TTCprobability = new float[m_maxTRK];
241
242 m_DTTindex = new int[m_maxDTT];
243 m_DTTsize = new int[m_maxDTT];
244 m_DTTtrajectory_eventindex = new std::vector<std::vector<int>>;
245 m_DTTtrajectory_barcode = new std::vector<std::vector<int>>;
246 m_DTTstTruth_subDetType = new std::vector<std::vector<int>>;
247 m_DTTstTrack_subDetType = new std::vector<std::vector<int>>;
248 m_DTTstCommon_subDetType = new std::vector<std::vector<int>>;
249
250 m_nt->Branch("run_number", &m_run_number, "run_number/i");
251 m_nt->Branch("event_number", &m_event_number, "event_number/l");
252
253 m_nt->Branch("nSE", &m_nSE, "nSE/I");
254 m_nt->Branch("SEID", m_SEID, "SEID[nSE]/I");
255
256 m_nt->Branch("nCL", &m_nCL, "nCL/I");
257 m_nt->Branch("CLindex", m_CLindex, "CLindex[nCL]/I");
258 m_nt->Branch("CLhardware", &m_CLhardware);
259 m_nt->Branch("CLx", m_CLx, "CLx[nCL]/D");
260 m_nt->Branch("CLy", m_CLy, "CLy[nCL]/D");
261 m_nt->Branch("CLz", m_CLz, "CLz[nCL]/D");
262 m_nt->Branch("CLbarrel_endcap", m_CLbarrel_endcap, "CLbarrel_endcap[nCL]/I");
263 m_nt->Branch("CLlayer_disk", m_CLlayer_disk, "CLlayer_disk[nCL]/I");
264 m_nt->Branch("CLeta_module", m_CLeta_module, "CLeta_module[nCL]/I");
265 m_nt->Branch("CLphi_module", m_CLphi_module, "CLphi_module[nCL]/I");
266 m_nt->Branch("CLside", m_CLside, "CLside[nCL]/I");
267 m_nt->Branch("CLmoduleID", m_CLmoduleID, "CLmoduleID[nCL]/l");
268 m_nt->Branch("CLparticleLink_eventIndex", &m_CLparticleLink_eventIndex);
269 m_nt->Branch("CLparticleLink_barcode", &m_CLparticleLink_barcode);
270 m_nt->Branch("CLbarcodesLinked", &m_CLbarcodesLinked);
271 m_nt->Branch("CLparticle_charge", &m_CLparticle_charge);
272 m_nt->Branch("CLphis", &m_CLphis);
273 m_nt->Branch("CLetas", &m_CLetas);
274 m_nt->Branch("CLtots", &m_CLtots);
275 m_nt->Branch("CLloc_direction1", m_CLloc_direction1, "CLloc_direction1[nCL]/D");
276 m_nt->Branch("CLloc_direction2", m_CLloc_direction2, "CLloc_direction2[nCL]/D");
277 m_nt->Branch("CLloc_direction3", m_CLloc_direction3, "CLloc_direction3[nCL]/D");
278 m_nt->Branch("CLJan_loc_direction1", m_CLJan_loc_direction1, "CLJan_loc_direction1[nCL]/D");
279 m_nt->Branch("CLJan_loc_direction2", m_CLJan_loc_direction2, "CLJan_loc_direction2[nCL]/D");
280 m_nt->Branch("CLJan_loc_direction3", m_CLJan_loc_direction3, "CLJan_loc_direction3[nCL]/D");
281 m_nt->Branch("CLpixel_count", m_CLpixel_count, "CLpixel_count[nCL]/I");
282 m_nt->Branch("CLcharge_count", m_CLcharge_count, "CLcharge_count[nCL]/F");
283 m_nt->Branch("CLloc_eta", m_CLloc_eta, "CLloc_eta[nCL]/F");
284 m_nt->Branch("CLloc_phi", m_CLloc_phi, "CLloc_phi[nCL]/F");
285 m_nt->Branch("CLglob_eta", m_CLglob_eta, "CLglob_eta[nCL]/F");
286 m_nt->Branch("CLglob_phi", m_CLglob_phi, "CLglob_phi[nCL]/F");
287 m_nt->Branch("CLeta_angle", m_CLeta_angle, "CLeta_angle[nCL]/D");
288 m_nt->Branch("CLphi_angle", m_CLphi_angle, "CLphi_angle[nCL]/D");
289 m_nt->Branch("CLnorm_x", m_CLnorm_x, "CLnorm_x[nCL]/F");
290 m_nt->Branch("CLnorm_y", m_CLnorm_y, "CLnorm_y[nCL]/F");
291 m_nt->Branch("CLnorm_z", m_CLnorm_z, "CLnorm_z[nCL]/F");
292 m_nt->Branch("CLlocal_cov", &m_CLlocal_cov);
293
294 m_nt->Branch("nPartEVT", &m_nPartEVT, "nPartEVT/I");
295 m_nt->Branch("Part_event_number", m_Part_event_number, "Part_event_number[nPartEVT]/I");
296 m_nt->Branch("Part_barcode", m_Part_barcode, "Part_barcode[nPartEVT]/I");
297 m_nt->Branch("Part_px", m_Part_px, "Part_px[nPartEVT]/F");
298 m_nt->Branch("Part_py", m_Part_py, "Part_py[nPartEVT]/F");
299 m_nt->Branch("Part_pz", m_Part_pz, "Part_pz[nPartEVT]/F");
300 m_nt->Branch("Part_pt", m_Part_pt, "Part_pt[nPartEVT]/F");
301 m_nt->Branch("Part_eta", m_Part_eta, "Part_eta[nPartEVT]/F");
302 m_nt->Branch("Part_vx", m_Part_vx, "Part_vx[nPartEVT]/F");
303 m_nt->Branch("Part_vy", m_Part_vy, "Part_vy[nPartEVT]/F");
304 m_nt->Branch("Part_vz", m_Part_vz, "Part_vz[nPartEVT]/F");
305 m_nt->Branch("Part_radius", m_Part_radius, "Part_radius[nPartEVT]/F");
306 m_nt->Branch("Part_status", m_Part_status, "Part_status[nPartEVT]/F");
307 m_nt->Branch("Part_charge", m_Part_charge, "Part_charge[nPartEVT]/F");
308 m_nt->Branch("Part_pdg_id", m_Part_pdg_id, "Part_pdg_id[nPartEVT]/I");
309 m_nt->Branch("Part_passed", m_Part_passed, "Part_passed[nPartEVT]/I");
310 m_nt->Branch("Part_vProdNin", m_Part_vProdNin, "Part_vProdNin[nPartEVT]/I");
311 m_nt->Branch("Part_vProdNout", m_Part_vProdNout, "Part_vProdNout[nPartEVT]/I");
312 m_nt->Branch("Part_vProdStatus", m_Part_vProdStatus, "Part_vProdStatus[nPartEVT]/I");
313 m_nt->Branch("Part_vProdBarcode", m_Part_vProdBarcode, "Part_vProdBarcode[nPartEVT]/I");
314 m_nt->Branch("Part_vParentID", &m_Part_vParentID);
315 m_nt->Branch("Part_vParentBarcode", &m_Part_vParentBarcode);
316
317 m_nt->Branch("nSP", &m_nSP, "nSP/I");
318 m_nt->Branch("SPindex", m_SPindex, "SPindex[nSP]/I");
319 m_nt->Branch("SPx", m_SPx, "SPx[nSP]/D");
320 m_nt->Branch("SPy", m_SPy, "SPy[nSP]/D");
321 m_nt->Branch("SPz", m_SPz, "SPz[nSP]/D");
322 m_nt->Branch("SPCL1_index", m_SPCL1_index, "SPCL1_index[nSP]/I");
323 m_nt->Branch("SPCL2_index", m_SPCL2_index, "SPCL2_index[nSP]/I");
324 m_nt->Branch("SPisOverlap", m_SPisOverlap, "SPisOverlap[nSP]/I");
325 m_nt->Branch("SPradius",m_SPradius, "SPradius[nSP]/D");
326 m_nt->Branch("SPcovr",m_SPcovr, "SPradius[nSP]/D");
327 m_nt->Branch("SPcovz",m_SPcovz, "SPradius[nSP]/D");
328 m_nt->Branch("SPhl_topstrip",m_SPhl_topstrip, "SPhl_topstrip[nSP]/F");
329 m_nt->Branch("SPhl_botstrip",m_SPhl_botstrip, "SPhl_botstrip[nSP]/F");
330 m_nt->Branch("SPtopStripDirection",&m_SPtopStripDirection);
331 m_nt->Branch("SPbottomStripDirection",&m_SPbottomStripDirection);
332 m_nt->Branch("SPstripCenterDistance",&m_SPstripCenterDistance);
333 m_nt->Branch("SPtopStripCenterPosition",m_SPtopStripCenterPosition);
334
335 m_nt->Branch("nTRK", &m_nTRK, "nTRK/I");
336 m_nt->Branch("TRKindex", m_TRKindex, "TRKindex[nTRK]/I");
337 m_nt->Branch("TRKtrack_fitter", m_TRKtrack_fitter, "TRKtrack_fitter[nTRK]/I");
338 m_nt->Branch("TRKparticle_hypothesis", m_TRKparticle_hypothesis, "TRKparticle_hypothesis[nTRK]/I");
339 m_nt->Branch("TRKproperties", &m_TRKproperties);
340 m_nt->Branch("TRKpattern", &m_TRKpattern);
341 m_nt->Branch("TRKndof", m_TRKndof, "TRKndof[nTRK]/I");
342 m_nt->Branch("TRKmot", m_TRKmot, "TRKmot[nTRK]/I");
343 m_nt->Branch("TRKoot", m_TRKoot, "TRKoot[nTRK]/I");
344 m_nt->Branch("TRKchiSq", m_TRKchiSq, "TRKchiSq[nTRK]/F");
345 m_nt->Branch("TRKmeasurementsOnTrack_pixcl_sctcl_index", &m_TRKmeasurementsOnTrack_pixcl_sctcl_index);
346 m_nt->Branch("TRKoutliersOnTrack_pixcl_sctcl_index", &m_TRKoutliersOnTrack_pixcl_sctcl_index);
347 m_nt->Branch("TRKcharge", m_TRKcharge, "TRKcharge[nTRK]/I");
348 m_nt->Branch("TRKperigee_position", &m_TRKperigee_position);
349 m_nt->Branch("TRKperigee_momentum", &m_TRKperigee_momentum);
350 m_nt->Branch("TTCindex", m_TTCindex, "TTCindex[nTRK]/I");
351 m_nt->Branch("TTCevent_index", m_TTCevent_index, "TTCevent_index[nTRK]/I");
352 m_nt->Branch("TTCparticle_link", m_TTCparticle_link, "TTCparticle_link[nTRK]/I");
353 m_nt->Branch("TTCprobability", m_TTCprobability, "TTCprobability[nTRK]/F");
354
355 m_nt->Branch("nDTT", &m_nDTT, "nDTT/I");
356 m_nt->Branch("DTTindex", m_DTTindex, "DTTindex[nDTT]/I");
357 m_nt->Branch("DTTsize", m_DTTsize, "DTTsize[nDTT]/I");
358 m_nt->Branch("DTTtrajectory_eventindex", &m_DTTtrajectory_eventindex);
359 m_nt->Branch("DTTtrajectory_barcode", &m_DTTtrajectory_barcode);
360 m_nt->Branch("DTTstTruth_subDetType", &m_DTTstTruth_subDetType);
361 m_nt->Branch("DTTstTrack_subDetType", &m_DTTstTrack_subDetType);
362 m_nt->Branch("DTTstCommon_subDetType", &m_DTTstCommon_subDetType);
363 }
364
365 return StatusCode::SUCCESS;
366}
367
368//-------------------------------
370 //-------------------------------
371 //
372 const EventContext &ctx = Gaudi::Hive::currentContext();
373
374 m_event++;
375
376 // map cluster ID to an index
377 // in order to connect the cluster to spacepoints
378 std::map<Identifier, long int> clusterIDMapIdx;
379 m_selected = 0; // global indices for clusters
380
381 std::map<Identifier, long int> clusterIDMapSpacePointIdx; // key: cluster indentifier, value: spacepoint index
382
383 // create a container with HepMcParticleLink and list of clusters
384 // particle barcode --> is accepted and number of clusters
385 std::map<std::pair<int, int>, std::pair<bool, int>> allTruthParticles;
386
387 const McEventCollection *mcCollptr = nullptr;
388 SG::ReadHandle<McEventCollection> mcEventCollectionHandle{m_mcEventCollectionKey, ctx};
389 if (not mcEventCollectionHandle.isValid()) {
390 ATH_MSG_WARNING(" McEventCollection not found: " << m_mcEventCollectionKey.key());
391 return StatusCode::FAILURE;
392 }
393 mcCollptr = mcEventCollectionHandle.cptr();
394
395 // dump out event ID
396 const xAOD::EventInfo *eventInfo = nullptr;
398 if (not eventInfoHandle.isValid()) {
399 ATH_MSG_WARNING(" EventInfo not found: " << m_eventInfoKey.key());
400 return StatusCode::FAILURE;
401 }
402 eventInfo = eventInfoHandle.cptr();
403
404 m_run_number = eventInfo->runNumber();
405 m_event_number = eventInfo->eventNumber();
406
407 std::map<int, int> allSubEvents;
408
409 m_nSE = 0;
410
411 bool duplicateSubeventID = false;
412 for (unsigned int cntr = 0; cntr < mcCollptr->size(); ++cntr) {
413 int ID = mcCollptr->at(cntr)->event_number();
414 if (m_rootFile)
415 m_SEID[m_nSE++] = ID;
416
417 if (m_nSE == m_maxCL) {
418 ATH_MSG_WARNING("DUMP : hit max number of subevent ID");
419 break;
420 }
421 std::map<int, int>::iterator it = allSubEvents.find(ID);
422 if (it == allSubEvents.end())
423 allSubEvents.insert(std::make_pair(ID, 1));
424 else {
425 it->second++;
426 duplicateSubeventID = true;
427 }
428 }
429
430 if (duplicateSubeventID) {
431 ATH_MSG_WARNING("Duplicate subevent ID in event " << m_event);
432 }
433
434 m_nPartEVT = 0;
435
436 if (m_rootFile) {
437 (*m_Part_vParentID).clear();
438 (*m_Part_vParentBarcode).clear();
439 }
440
441 for (unsigned int cntr = 0; cntr < mcCollptr->size(); ++cntr) {
442 const HepMC::GenEvent *genEvt = (mcCollptr->at(cntr));
443
444 // for ( HepMC::GenEvent::particle_const_iterator p = genEvt->particles_begin(); p != genEvt->particles_end(); ++p )
445 // {
446
450
451 for (auto p : *genEvt) {
452 //*p is a GenParticle
453 float px, py, pz, pt, eta, vx, vy, vz, radius, status, charge = 0.;
454 std::vector<int> vParentID;
455 std::vector<int> vParentBarcode;
456
457 int vProdNin, vProdNout, vProdStatus, vProdBarcode;
458 bool passed = isPassed(p, px, py, pz, pt, eta, vx, vy, vz, radius, status, charge, vParentID, vParentBarcode,
459 vProdNin, vProdNout, vProdStatus, vProdBarcode);
460 allTruthParticles.insert(std::make_pair(std::make_pair(genEvt->event_number(), HepMC::barcode(p)),
461 std::make_pair(passed, 0))); // JB: HEPMC3 barcode() -> HepMC::barcode(p)
462 // subevent, barcode, px, py, pz, pt, eta, vx, vy, vz, radius, status, charge
463 if (m_rootFile) {
464 m_Part_event_number[m_nPartEVT] = genEvt->event_number();
465 m_Part_barcode[m_nPartEVT] = HepMC::barcode(p); // JB: HEPMC3 barcode() -> HepMC::barcode(p)
466 m_Part_px[m_nPartEVT] = px;
467 m_Part_py[m_nPartEVT] = py;
468 m_Part_pz[m_nPartEVT] = pz;
469 m_Part_pt[m_nPartEVT] = pt;
471 m_Part_vx[m_nPartEVT] = vx;
472 m_Part_vy[m_nPartEVT] = vy;
473 m_Part_vz[m_nPartEVT] = vz;
474 m_Part_radius[m_nPartEVT] = radius;
475 m_Part_status[m_nPartEVT] = status;
477 m_Part_pdg_id[m_nPartEVT] = p->pdg_id();
478 m_Part_passed[m_nPartEVT] = (passed ? true : false);
479 m_Part_vProdNin[m_nPartEVT] = vProdNin;
480 m_Part_vProdNout[m_nPartEVT] = vProdNout;
481 m_Part_vProdStatus[m_nPartEVT] = vProdStatus;
482 m_Part_vProdBarcode[m_nPartEVT] = vProdBarcode;
483 (*m_Part_vParentID).push_back(vParentID);
484 (*m_Part_vParentBarcode).push_back(vParentBarcode);
485 }
486
487 m_nPartEVT++;
488 if (m_nPartEVT == m_maxPart) {
489 ATH_MSG_WARNING("DUMP : hit max number of particle events");
490 break;
491 }
492 }
493 }
494
495 const InDet::PixelClusterContainer *PixelClusterContainer = 0;
496 SG::ReadHandle<InDet::PixelClusterContainer> pixelClusterContainerHandle{m_pixelClusterKey, ctx};
497 if (not pixelClusterContainerHandle.isValid()) {
498 ATH_MSG_WARNING(" PixelClusterContainer not found: " << m_pixelClusterKey.key());
499 return StatusCode::FAILURE;
500 }
501 PixelClusterContainer = pixelClusterContainerHandle.cptr();
502
503 const InDet::SCT_ClusterContainer *SCT_ClusterContainer = 0;
504 SG::ReadHandle<InDet::SCT_ClusterContainer> stripClusterContainerHandle{m_stripClusterKey, ctx};
505 if (not stripClusterContainerHandle.isValid()) {
506 ATH_MSG_WARNING(" SCT_ClusterContainer not found: " << m_stripClusterKey.key());
507 return StatusCode::FAILURE;
508 }
509 SCT_ClusterContainer = stripClusterContainerHandle.cptr();
510
511 auto cartesion_to_spherical = [](const Amg::Vector3D &xyzVec, float &eta_, float &phi_) {
512 float r3 = 0;
513 for (int idx = 0; idx < 3; ++idx) {
514 r3 += xyzVec[idx] * xyzVec[idx];
515 }
516 r3 = sqrt(r3);
517 phi_ = atan2(xyzVec[1], xyzVec[0]);
518 float theta_ = acos(xyzVec[2] / r3);
519 eta_ = log(tan(0.5 * theta_));
520 };
521
525
526 m_nCL = 0;
527 if (m_rootFile) {
528 (*m_CLhardware).clear();
529 (*m_CLparticleLink_eventIndex).clear();
530 (*m_CLparticleLink_barcode).clear();
531 (*m_CLbarcodesLinked).clear();
532 (*m_CLparticle_charge).clear();
533 (*m_CLphis).clear();
534 (*m_CLetas).clear();
535 (*m_CLtots).clear();
536 (*m_CLlocal_cov).clear();
537 }
538
539 if (PixelClusterContainer->size() > 0) {
540
541 const InDetSimDataCollection *sdoCollection = 0;
543 if (not sdoCollectionHandle.isValid()) {
544 ATH_MSG_WARNING(" InDetSimDataCollection not found: " << m_pixelSDOKey.key());
545 return StatusCode::FAILURE;
546 }
547 sdoCollection = sdoCollectionHandle.cptr();
548
549 for (const auto clusterCollection : *PixelClusterContainer) {
550 // skip empty collections
551 if (clusterCollection->empty())
552 continue;
553
554 int barrel_endcap = m_pixelID->barrel_ec(clusterCollection->identify());
555 int layer_disk = m_pixelID->layer_disk(clusterCollection->identify());
556 int eta_module = m_pixelID->eta_module(clusterCollection->identify());
557 int phi_module = m_pixelID->phi_module(clusterCollection->identify());
558
559 const InDetDD::SiDetectorElement *element = m_pixelManager->getDetectorElement(clusterCollection->identify());
560
561 Amg::Vector3D my_normal = element->normal();
562 float norm_x = fabs(my_normal.x()) > 1e-5 ? my_normal.x() : 0.;
563 float norm_y = fabs(my_normal.y()) > 1e-5 ? my_normal.y() : 0.;
564 float norm_z = fabs(my_normal.z()) > 1e-5 ? my_normal.z() : 0.;
565
566 const InDetDD::PixelModuleDesign *design(dynamic_cast<const InDetDD::PixelModuleDesign *>(&element->design()));
567
568 if (not design) {
569 ATH_MSG_ERROR("Dynamic cast failed at " << __LINE__ << " of MergedPixelsTool.cxx.");
570 return StatusCode::FAILURE;
571 }
572
573 // loop over collection
574 for (const auto cluster : *clusterCollection) {
575 Identifier clusterId = cluster->identify();
576 if (!clusterId.is_valid()) {
577 ATH_MSG_WARNING("Pixel cluster identifier is not valid");
578 }
579
580 const Amg::MatrixX &local_cov = cluster->localCovariance();
581
582 std::vector<std::pair<int, int>> barcodes = {};
583 std::vector<int> particleLink_eventIndex = {};
584 std::vector<int> particleLink_barcode = {};
585 std::vector<bool> barcodesLinked = {};
586 std::vector<float> charge = {};
587 std::vector<int> phis = {};
588 std::vector<int> etas = {};
589 std::vector<int> tots = {};
590 int min_eta = 999;
591 int min_phi = 999;
592 int max_eta = -999;
593 int max_phi = -999;
594
595 float charge_count = 0;
596 int pixel_count = 0;
597
598 for (unsigned int rdo = 0; rdo < cluster->rdoList().size(); rdo++) {
599 const auto &rdoID = cluster->rdoList().at(rdo);
600 int phi = m_pixelID->phi_index(rdoID);
601 int eta = m_pixelID->eta_index(rdoID);
602 if (min_eta > eta)
603 min_eta = eta;
604 if (min_phi > phi)
605 min_phi = phi;
606 if (max_eta < eta)
607 max_eta = eta;
608 if (max_phi < phi)
609 max_phi = phi;
610
611 ++pixel_count;
612 charge_count += cluster->totList().at(rdo);
613
614 phis.push_back(phi);
615 etas.push_back(eta);
616 tots.push_back(cluster->totList().at(rdo));
617
618 auto pos = sdoCollection->find(rdoID);
619 if (pos != sdoCollection->end()) {
620 for (auto deposit : pos->second.getdeposits()) {
621 const HepMcParticleLink &particleLink = deposit.first;
622 std::pair<int, int> barcode(particleLink.eventIndex(), particleLink.barcode());
623 // if (particleLink.isValid()) allTruthParticles.at(barcode).second++; // JB comment this out
624 if (std::find(barcodes.begin(), barcodes.end(), barcode) == barcodes.end()) {
625 barcodes.push_back(barcode);
626 particleLink_eventIndex.push_back(particleLink.eventIndex());
627 particleLink_barcode.push_back(particleLink.barcode());
628 charge.push_back(deposit.second);
629 barcodesLinked.push_back(particleLink.isValid());
630 }
631 }
632 }
633 }
634
635 InDetDD::SiLocalPosition localPos_entry = design->localPositionOfCell(InDetDD::SiCellId(min_phi, min_eta));
636 InDetDD::SiLocalPosition localPos_exit = design->localPositionOfCell(InDetDD::SiCellId(max_phi, max_eta));
637
638 Amg::Vector3D localStartPosition(localPos_entry.xEta() - 0.5 * element->etaPitch(),
639 localPos_entry.xPhi() - 0.5 * element->phiPitch(),
640 -0.5 * element->thickness());
641 Amg::Vector3D localEndPosition(localPos_exit.xEta() + 0.5 * element->etaPitch(),
642 localPos_exit.xPhi() + 0.5 * element->phiPitch(), 0.5 * element->thickness());
643
644 // local direction in local coordinates
645 // clusterShape: [lx, ly, lz]
646 Amg::Vector3D localDirection = localEndPosition - localStartPosition;
647
648 float loc_eta = 0, loc_phi = 0; // clusterShape: [leta, lphi]
649 cartesion_to_spherical(localDirection, loc_eta, loc_phi);
650
651 Amg::Vector3D globalStartPosition = element->globalPosition(localStartPosition);
652 Amg::Vector3D globalEndPosition = element->globalPosition(localEndPosition);
653
654 Amg::Vector3D direction = globalEndPosition - globalStartPosition;
655 float glob_eta = 0, glob_phi = 0; // clusterShape: [geta, gphi]
656 cartesion_to_spherical(direction, glob_eta, glob_phi);
657
658 Amg::Vector3D my_phiax = element->phiAxis();
659 Amg::Vector3D my_etaax = element->etaAxis();
660
661 float trkphicomp = direction.dot(my_phiax);
662 float trketacomp = direction.dot(my_etaax);
663 float trknormcomp = direction.dot(my_normal);
664 double phi_angle = atan2(trknormcomp, trkphicomp);
665 double eta_angle = atan2(trknormcomp, trketacomp);
666 // now dumping all the values now
667 clusterIDMapIdx[cluster->identify()] = m_selected;
668 std::vector<double> v_local_cov;
669 if (local_cov.size() > 0) {
670 for (size_t i = 0, nRows = local_cov.rows(), nCols = local_cov.cols(); i < nRows; i++) {
671 for (size_t j = 0; j < nCols; ++j) {
672 v_local_cov.push_back(local_cov(i, j));
673 }
674 }
675 }
676 if (m_rootFile) {
677 // fill TTree
679 (*m_CLhardware).push_back("PIXEL");
680 m_CLx[m_nCL] = cluster->globalPosition().x();
681 m_CLy[m_nCL] = cluster->globalPosition().y();
682 m_CLz[m_nCL] = cluster->globalPosition().z();
683 m_CLbarrel_endcap[m_nCL] = barrel_endcap;
684 m_CLlayer_disk[m_nCL] = layer_disk;
685 m_CLeta_module[m_nCL] = eta_module;
686 m_CLphi_module[m_nCL] = phi_module;
687 m_CLside[m_nCL] = 0;
688 m_CLmoduleID[m_nCL] = clusterCollection->identify().get_compact();
689 (*m_CLparticleLink_eventIndex).push_back(particleLink_eventIndex);
690 (*m_CLparticleLink_barcode).push_back(particleLink_barcode);
691 (*m_CLbarcodesLinked).push_back(barcodesLinked);
692 (*m_CLparticle_charge).push_back(charge);
693 (*m_CLetas).push_back(etas);
694 (*m_CLphis).push_back(phis);
695 (*m_CLtots).push_back(tots);
696 m_CLloc_direction1[m_nCL] = localDirection[0];
697 m_CLloc_direction2[m_nCL] = localDirection[1];
698 m_CLloc_direction3[m_nCL] = localDirection[2];
702 m_CLpixel_count[m_nCL] = pixel_count;
703 m_CLcharge_count[m_nCL] = charge_count;
704 m_CLloc_eta[m_nCL] = loc_eta;
705 m_CLloc_phi[m_nCL] = loc_phi;
706 m_CLglob_eta[m_nCL] = glob_eta;
707 m_CLglob_phi[m_nCL] = glob_phi;
708 m_CLeta_angle[m_nCL] = eta_angle;
709 m_CLphi_angle[m_nCL] = phi_angle;
710 m_CLnorm_x[m_nCL] = norm_x;
711 m_CLnorm_y[m_nCL] = norm_y;
712 m_CLnorm_z[m_nCL] = norm_z;
713 (*m_CLlocal_cov).push_back(v_local_cov);
714 }
715 m_nCL++;
716 m_selected++;
717 if (m_nCL == m_maxCL) {
718 ATH_MSG_WARNING("DUMP : hit max number of clusters");
719 break;
720 }
721 }
722 }
723 }
724
728
729 if (SCT_ClusterContainer->size() > 0) {
730 const InDetSimDataCollection *sdoCollection = 0;
732 if (not sdoCollectionHandle.isValid()) {
733 ATH_MSG_WARNING(" InDetSimDataCollection not found: " << m_stripSDOKey.key());
734 return StatusCode::FAILURE;
735 }
736 sdoCollection = sdoCollectionHandle.cptr();
737
738 for (const auto clusterCollection : *SCT_ClusterContainer) {
739 // skip empty collections
740 if (clusterCollection->empty())
741 continue;
742
743 int barrel_endcap = m_SCT_ID->barrel_ec(clusterCollection->identify());
744 int layer_disk = m_SCT_ID->layer_disk(clusterCollection->identify());
745 int eta_module = m_SCT_ID->eta_module(clusterCollection->identify());
746 int phi_module = m_SCT_ID->phi_module(clusterCollection->identify());
747 int side = m_SCT_ID->side(clusterCollection->identify());
748
749 const InDetDD::SiDetectorElement *element = m_SCT_Manager->getDetectorElement(clusterCollection->identify());
750
751 Amg::Vector3D my_normal = element->normal();
752 float norm_x = fabs(my_normal.x()) > 1e-5 ? my_normal.x() : 0.;
753 float norm_y = fabs(my_normal.y()) > 1e-5 ? my_normal.y() : 0.;
754 float norm_z = fabs(my_normal.z()) > 1e-5 ? my_normal.z() : 0.;
755
756 // loop over collection
757 for (const auto cluster : *clusterCollection) {
758 Identifier clusterId = cluster->identify();
759 if (!clusterId.is_valid()) {
760 ATH_MSG_WARNING("SCT cluster identifier is not valid");
761 }
762
763 const Amg::MatrixX &local_cov = cluster->localCovariance();
764
765 std::vector<std::pair<int, int>> barcodes = {};
766 std::vector<int> particleLink_eventIndex = {};
767 std::vector<int> particleLink_barcode = {};
768 std::vector<bool> barcodesLinked = {};
769 std::vector<float> charge = {};
770
771 std::vector<int> tots = {};
772 std::vector<int> strip_ids = {};
773 int min_strip = 999;
774 int max_strip = -999;
775
776 float charge_count = 0;
777 int pixel_count = 0;
778
779 for (unsigned int rdo = 0; rdo < cluster->rdoList().size(); rdo++) {
780 const auto &rdoID = cluster->rdoList().at(rdo);
781
782 int strip = m_SCT_ID->strip(rdoID);
783
784 if (min_strip > strip)
785 min_strip = strip;
786 if (max_strip < strip)
787 max_strip = strip;
788 strip_ids.push_back(strip);
789 // tots.push_back(cluster->totList().at(rdo));
790 tots.push_back(0); // FIXME
791 ++pixel_count;
792 // find barcodes of the truth particles
793 auto pos = sdoCollection->find(rdoID);
794 if (pos != sdoCollection->end()) {
795 for (auto deposit : pos->second.getdeposits()) {
796 const HepMcParticleLink &particleLink = deposit.first;
797 std::pair<int, int> barcode(particleLink.eventIndex(), particleLink.barcode());
798 // note that we are not filling the map allTruthParticles here - OK, we are not using this map for
799 // anything
800 if (std::find(barcodes.begin(), barcodes.end(), barcode) == barcodes.end()) {
801 barcodes.push_back(barcode);
802 particleLink_eventIndex.push_back(particleLink.eventIndex());
803 particleLink_barcode.push_back(particleLink.barcode());
804 charge.push_back(deposit.second);
805 barcodesLinked.push_back(particleLink.isValid());
806 }
807 }
808 }
809 }
810
811 // retrieve cluster shape
812 const InDetDD::SCT_ModuleSideDesign *design(
813 dynamic_cast<const InDetDD::SCT_ModuleSideDesign *>(&element->design()));
814 if (not design) {
815 ATH_MSG_ERROR("Failed at " << __LINE__ << " of accessing SCT ModuleSide Design");
816 return StatusCode::FAILURE;
817 }
818
819 Amg::Vector2D locpos = cluster->localPosition();
820 std::pair<Amg::Vector3D, Amg::Vector3D> ends(
821 element->endsOfStrip(InDetDD::SiLocalPosition(locpos.y(), locpos.x(), 0)));
822
823 Amg::Vector3D JanDirection = ends.second - ends.first;
824
825 InDetDD::SiLocalPosition localPos_entry = design->localPositionOfCell(InDetDD::SiCellId(min_strip));
826 InDetDD::SiLocalPosition localPos_exit = design->localPositionOfCell(InDetDD::SiCellId(max_strip));
827
828 Amg::Vector3D localStartPosition(localPos_entry.xEta() - 0.5 * element->etaPitch(),
829 localPos_entry.xPhi() - 0.5 * element->phiPitch(),
830 -0.5 * element->thickness());
831 Amg::Vector3D localEndPosition(localPos_exit.xEta() + 0.5 * element->etaPitch(),
832 localPos_exit.xPhi() + 0.5 * element->phiPitch(), 0.5 * element->thickness());
833
834 Amg::Vector3D localDirection = localEndPosition - localStartPosition;
835 float loc_eta = 0, loc_phi = 0; // clusterShape: [leta, lphi]
836 cartesion_to_spherical(localDirection, loc_eta, loc_phi);
837
838 Amg::Vector3D globalStartPosition = element->globalPosition(localStartPosition);
839 Amg::Vector3D globalEndPosition = element->globalPosition(localEndPosition);
840
841 Amg::Vector3D direction = globalEndPosition - globalStartPosition;
842 float glob_eta = 0, glob_phi = 0; // clusterShape: [geta, gphi]
843 cartesion_to_spherical(direction, glob_eta, glob_phi);
844
845 Amg::Vector3D my_phiax = element->phiAxis();
846 Amg::Vector3D my_etaax = element->etaAxis();
847
848 float trkphicomp = direction.dot(my_phiax);
849 float trketacomp = direction.dot(my_etaax);
850 float trknormcomp = direction.dot(my_normal);
851 double phi_angle = atan2(trknormcomp, trkphicomp);
852 double eta_angle = atan2(trknormcomp, trketacomp);
853
854 // now dumping all the values now
855 clusterIDMapIdx[cluster->identify()] = m_selected;
856 // cluster shape
857 std::vector<int> cst;
858 for (unsigned strip = 0; strip < strip_ids.size(); strip++) {
859 cst.push_back(-1);
860 }
861 std::vector<double> v_local_cov;
862 if (local_cov.size() > 0) {
863 for (size_t i = 0, nRows = local_cov.rows(), nCols = local_cov.cols(); i < nRows; i++) {
864 for (size_t j = 0; j < nCols; ++j) {
865 v_local_cov.push_back(local_cov(i, j));
866 }
867 }
868 }
869 if (m_rootFile) {
871 (*m_CLhardware).push_back("STRIP");
872 m_CLx[m_nCL] = cluster->globalPosition().x();
873 m_CLy[m_nCL] = cluster->globalPosition().y();
874 m_CLz[m_nCL] = cluster->globalPosition().z();
875 m_CLbarrel_endcap[m_nCL] = barrel_endcap;
876 m_CLlayer_disk[m_nCL] = layer_disk;
877 m_CLeta_module[m_nCL] = eta_module;
878 m_CLphi_module[m_nCL] = phi_module;
879 m_CLside[m_nCL] = side;
880 m_CLmoduleID[m_nCL] = clusterCollection->identify().get_compact();
881 (*m_CLparticleLink_eventIndex).push_back(particleLink_eventIndex);
882 (*m_CLparticleLink_barcode).push_back(particleLink_barcode);
883 (*m_CLbarcodesLinked).push_back(barcodesLinked);
884 (*m_CLparticle_charge).push_back(charge);
885 (*m_CLetas).push_back(strip_ids);
886 (*m_CLphis).push_back(cst);
887 (*m_CLtots).push_back(tots);
888 m_CLloc_direction1[m_nCL] = localDirection[0];
889 m_CLloc_direction2[m_nCL] = localDirection[1];
890 m_CLloc_direction3[m_nCL] = localDirection[2];
891 m_CLJan_loc_direction1[m_nCL] = JanDirection[0];
892 m_CLJan_loc_direction2[m_nCL] = JanDirection[1];
893 m_CLJan_loc_direction3[m_nCL] = JanDirection[2];
894 m_CLpixel_count[m_nCL] = pixel_count;
895 m_CLcharge_count[m_nCL] = charge_count;
896 m_CLloc_eta[m_nCL] = loc_eta;
897 m_CLloc_phi[m_nCL] = loc_phi;
898 m_CLglob_eta[m_nCL] = glob_eta;
899 m_CLglob_phi[m_nCL] = glob_phi;
900 m_CLeta_angle[m_nCL] = eta_angle;
901 m_CLphi_angle[m_nCL] = phi_angle;
902 m_CLnorm_x[m_nCL] = norm_x;
903 m_CLnorm_y[m_nCL] = norm_y;
904 m_CLnorm_z[m_nCL] = norm_z;
905 (*m_CLlocal_cov).push_back(v_local_cov);
906 }
907
908 m_nCL++;
909 m_selected++;
910 if (m_nCL == m_maxCL) {
911 ATH_MSG_WARNING("DUMP : hit max number of clusters");
912 break;
913 }
914 }
915 }
916 }
917
918
922
923 static const SG::Accessor< ElementLink<SpacePointCollection> > linkAcc("pixelSpacePointLink");
924 static const SG::Accessor< ElementLink< ::SpacePointCollection > > striplinkAcc("sctSpacePointLink");
925 static const SG::Accessor< ElementLink< ::SpacePointOverlapCollection > > stripOverlaplinkAcc("stripOverlapSpacePointLink");
926
927 // xAOD Containers
928 const xAOD::SpacePointContainer *xAODPixelSPContainer = nullptr;
929
931
932 if (not xAODPixelSpacePointContainerHandle.isValid()) {
933 ATH_MSG_ERROR(" SpacePointContainer not found: " << m_xaodPixelSpacePointContainerKey.key());
934 return StatusCode::FAILURE;
935 }
936
937 xAODPixelSPContainer = xAODPixelSpacePointContainerHandle.cptr();
938
939
940 const xAOD::SpacePointContainer *xAODStripSPContainer = 0;
942 if (not xAODStripSpacePointContainerHandle.isValid()) {
943 ATH_MSG_ERROR(" SpacePointContainer not found: " << m_xaodStripSpacePointContainerKey.key());
944 return StatusCode::FAILURE;
945 }
946 xAODStripSPContainer = xAODStripSpacePointContainerHandle.cptr();
947
948
949 const xAOD::SpacePointContainer *xAODStripSPOverlapContainer = 0;
950 SG::ReadHandle<xAOD::SpacePointContainer> xAODStripSpacePointOverlapContainerHandle{m_xaodStripSpacePointOverlapContainerKey, ctx};
951 if (not xAODStripSpacePointOverlapContainerHandle.isValid()) {
952 ATH_MSG_ERROR(" SpacePointContainer not found: " << m_xaodStripSpacePointOverlapContainerKey.key());
953 return StatusCode::FAILURE;
954 }
955 xAODStripSPOverlapContainer = xAODStripSpacePointOverlapContainerHandle.cptr();
956
957 int sp_index = 0;
958 m_nSP = 0;
959
960 if (xAODPixelSPContainer && xAODPixelSPContainer->size() > 0) {
961 for (const auto sp : *xAODPixelSPContainer) {
962
963 if (not linkAcc.isAvailable(*sp))
964 ATH_MSG_FATAL("no pixel SpacePoint link for xAOD::SpacePoint");
965
966
967 auto trk_sp = *linkAcc(*sp);
968 const InDet::SiCluster *cl = static_cast<const InDet::SiCluster*>(trk_sp->clusterList().first);
969
970 if (m_rootFile) {
971 m_SPindex[m_nSP] = sp_index;
972 m_SPx[m_nSP] = sp->globalPosition().x();
973 m_SPy[m_nSP] = sp->globalPosition().y();
974 m_SPz[m_nSP] = sp->globalPosition().z();
975 m_SPradius[m_nSP] = sp->radius();
976 m_SPcovr[m_nSP] = sp->varianceR();
977 m_SPcovz[m_nSP] = sp->varianceZ();
978 m_SPCL1_index[m_nSP] = clusterIDMapIdx[cl->identify()];
979 m_SPCL2_index[m_nSP] = -1;
980 m_SPisOverlap[m_nSP] = -1;
981 }
982
983 sp_index++;
984 m_nSP++;
985 if (m_nSP == m_maxSP) {
986 ATH_MSG_WARNING("DUMP : hit max number of space points");
987 break;
988 }
989 } // loop on container
990 } // container not empty
991
992 if (xAODStripSPContainer && xAODStripSPContainer->size() > 0) {
993
994 //loop over collection
995 for (const auto sp : *xAODStripSPContainer) {
996
997 ATH_CHECK(striplinkAcc.isAvailable(*sp));
998
999 auto trk_sp = *striplinkAcc(*sp);
1000 const InDet::SiCluster *cl_1 = static_cast<const InDet::SiCluster *>(trk_sp->clusterList().first);
1001 const InDet::SiCluster *cl_2 = static_cast<const InDet::SiCluster *>(trk_sp->clusterList().second);
1002
1003 if (m_rootFile) {
1004
1005 m_SPindex[m_nSP] = sp_index;
1006 m_SPx[m_nSP] = sp->globalPosition().x();
1007 m_SPy[m_nSP] = sp->globalPosition().y();
1008 m_SPz[m_nSP] = sp->globalPosition().z();
1009 m_SPradius[m_nSP] = sp->radius();
1010 m_SPcovr[m_nSP] = sp->varianceR();
1011 m_SPcovz[m_nSP] = sp->varianceZ();
1012 m_SPCL1_index[m_nSP] = clusterIDMapIdx[cl_1->identify()];
1013 m_SPCL2_index[m_nSP] = clusterIDMapIdx[cl_2->identify()];
1014 m_SPisOverlap[m_nSP] = 0;
1015 m_SPhl_topstrip[m_nSP] = sp->topHalfStripLength();
1016 m_SPhl_botstrip[m_nSP] = sp->bottomHalfStripLength();
1017
1018
1019 std::vector<float> topstripDir(sp->topStripDirection().data(),
1020 sp->topStripDirection().data() +
1021 sp->topStripDirection().size());
1022
1023 std::vector<float> botstripDir(sp->bottomStripDirection().data(),
1024 sp->bottomStripDirection().data() +
1025 sp->bottomStripDirection().size());
1026
1027 std::vector<float> DstripCnt(sp->stripCenterDistance().data(),
1028 sp->stripCenterDistance().data() +
1029 sp->stripCenterDistance().size());
1030
1031 std::vector<float> topstripCnt(sp->topStripCenter().data(),
1032 sp->topStripCenter().data() +
1033 sp->topStripCenter().size());
1034
1035 (*m_SPtopStripDirection).push_back(topstripDir);
1036 (*m_SPbottomStripDirection).push_back(botstripDir);
1037 (*m_SPstripCenterDistance).push_back(DstripCnt);
1038 (*m_SPtopStripCenterPosition).push_back(topstripCnt);
1039
1040 }
1041
1042 sp_index++;
1043 m_nSP++;
1044
1045 if (m_nSP == m_maxSP) {
1046 ATH_MSG_WARNING("DUMP : hit max number of space points");
1047 break;
1048 }
1049 }
1050 }
1051
1052
1053 if (xAODStripSPOverlapContainer && xAODStripSPOverlapContainer->size() > 0) {
1054
1055 //loop over collection
1056 for (const auto sp : *xAODStripSPOverlapContainer) {
1057
1058 ATH_CHECK(stripOverlaplinkAcc.isAvailable(*sp));
1059
1060 auto trk_sp = *stripOverlaplinkAcc(*sp);
1061 const InDet::SiCluster *cl_1 = static_cast<const InDet::SiCluster *>(trk_sp->clusterList().first);
1062 const InDet::SiCluster *cl_2 = static_cast<const InDet::SiCluster *>(trk_sp->clusterList().second);
1063
1064 if (m_rootFile) {
1065
1066 m_SPindex[m_nSP] = sp_index;
1067 m_SPx[m_nSP] = sp->globalPosition().x();
1068 m_SPy[m_nSP] = sp->globalPosition().y();
1069 m_SPz[m_nSP] = sp->globalPosition().z();
1070 m_SPradius[m_nSP] = sp->radius();
1071 m_SPcovr[m_nSP] = sp->varianceR();
1072 m_SPcovz[m_nSP] = sp->varianceZ();
1073 m_SPCL1_index[m_nSP] = clusterIDMapIdx[cl_1->identify()];
1074 m_SPCL2_index[m_nSP] = clusterIDMapIdx[cl_2->identify()];
1075
1076 int flag = compute_overlap_SP_flag(m_CLeta_module[clusterIDMapIdx[cl_1->identify()]],
1077 m_CLphi_module[clusterIDMapIdx[cl_1->identify()]],
1078 m_CLeta_module[clusterIDMapIdx[cl_2->identify()]],
1079 m_CLphi_module[clusterIDMapIdx[cl_2->identify()]]);
1080
1081 if ( flag<1 || flag > 3 )
1082 ATH_MSG_WARNING("Unexpected overlap SP flag: "<<flag);
1083
1084
1085 m_SPisOverlap[m_nSP] = flag;
1086 m_SPhl_topstrip[m_nSP] = sp->topHalfStripLength();
1087 m_SPhl_botstrip[m_nSP] = sp->bottomHalfStripLength();
1088
1089
1090 std::vector<float> topstripDir(sp->topStripDirection().data(),
1091 sp->topStripDirection().data() +
1092 sp->topStripDirection().size());
1093
1094 std::vector<float> botstripDir(sp->bottomStripDirection().data(),
1095 sp->bottomStripDirection().data() +
1096 sp->bottomStripDirection().size());
1097
1098 std::vector<float> DstripCnt(sp->stripCenterDistance().data(),
1099 sp->stripCenterDistance().data() +
1100 sp->stripCenterDistance().size());
1101
1102 std::vector<float> topstripCnt(sp->topStripCenter().data(),
1103 sp->topStripCenter().data() +
1104 sp->topStripCenter().size());
1105
1106 (*m_SPtopStripDirection).push_back(topstripDir);
1107 (*m_SPbottomStripDirection).push_back(botstripDir);
1108 (*m_SPstripCenterDistance).push_back(DstripCnt);
1109 (*m_SPtopStripCenterPosition).push_back(topstripCnt);
1110
1111 }
1112
1113 sp_index++;
1114 m_nSP++;
1115 if (m_nSP == m_maxSP) {
1116 ATH_MSG_WARNING("DUMP : hit max number of space points");
1117 break;
1118 }
1119 } // loop on container
1120 } // container not empty
1121
1122
1126
1127 const TrackCollection *trackCollection = 0;
1128 SG::ReadHandle<TrackCollection> trackCollectionHandle{m_tracksKey, ctx};
1129 if (not trackCollectionHandle.isValid()) {
1130 ATH_MSG_WARNING(" TrackCollection not found: " << m_tracksKey.key());
1131 return StatusCode::FAILURE;
1132 }
1133 trackCollection = trackCollectionHandle.cptr();
1134
1135 const TrackTruthCollection *trackTruthCollection = 0;
1136 SG::ReadHandle<TrackTruthCollection> trackTruthCollectionHandle{m_tracksTruthKey, ctx};
1137 if (not trackTruthCollectionHandle.isValid()) {
1138 ATH_MSG_WARNING(" TrackTruthCollection not found: " << m_tracksTruthKey.key());
1139 return StatusCode::FAILURE;
1140 }
1141 trackTruthCollection = trackTruthCollectionHandle.cptr();
1142
1143 int trk_index = 0;
1144
1145 // loop over tracks (and track truth) objects
1146 TrackCollection::const_iterator trackIterator = (*trackCollection).begin();
1147 m_nTRK = 0;
1148 if (m_rootFile) {
1149 (*m_TRKproperties).clear();
1150 (*m_TRKpattern).clear();
1151 (*m_TRKperigee_position).clear();
1152 (*m_TRKperigee_momentum).clear();
1153 (*m_TRKmeasurementsOnTrack_pixcl_sctcl_index).clear();
1154 (*m_TRKoutliersOnTrack_pixcl_sctcl_index).clear();
1155 }
1156
1157 for (; trackIterator < (*trackCollection).end(); ++trackIterator) {
1158 if (!((*trackIterator))) {
1159 ATH_MSG_WARNING("TrackCollection contains empty entries");
1160 continue;
1161 }
1162 const Trk::TrackInfo &info = (*trackIterator)->info();
1163 const Trk::FitQuality *fitQuality = (*trackIterator)->fitQuality();
1164 const Trk::Perigee *perigeeParameters = (*trackIterator)->perigeeParameters();
1165 const DataVector<const Trk::MeasurementBase> *measurementsOnTrack = (*trackIterator)->measurementsOnTrack();
1166 const DataVector<const Trk::MeasurementBase> *outliersOnTrack = (*trackIterator)->outliersOnTrack();
1167
1169 tracklink.setElement(const_cast<Trk::Track *>(*trackIterator));
1170 tracklink.setStorableObject(*trackCollection);
1171 const ElementLink<TrackCollection> tracklink2 = tracklink;
1172 TrackTruthCollection::const_iterator found = trackTruthCollection->find(tracklink2);
1173
1174 const std::bitset<Trk::TrackInfo::NumberOfTrackProperties> &properties = info.properties();
1175 std::vector<int> v_properties;
1176 for (std::size_t i = 0; i < properties.size(); i++) {
1177 if (properties[i]) {
1178 v_properties.push_back(i);
1179 }
1180 }
1181
1182 const std::bitset<Trk::TrackInfo::NumberOfTrackRecoInfo> &pattern = info.patternRecognition();
1183 std::vector<int> v_pattern;
1184 for (std::size_t i = 0; i < pattern.size(); i++) {
1185 if (pattern[i]) {
1186 v_pattern.push_back(i);
1187 }
1188 }
1189
1190 int ndof = -1;
1191 float chiSq = 0;
1192 if (fitQuality) {
1193 ndof = fitQuality->numberDoF();
1194 chiSq = fitQuality->chiSquared();
1195 }
1196 std::vector<double> position, momentum;
1197 int charge = 0;
1198 if (perigeeParameters) {
1199 position.push_back(perigeeParameters->position()[0]);
1200 position.push_back(perigeeParameters->position()[1]);
1201 position.push_back(perigeeParameters->position()[2]);
1202 momentum.push_back(perigeeParameters->momentum()[0]);
1203 momentum.push_back(perigeeParameters->momentum()[1]);
1204 momentum.push_back(perigeeParameters->momentum()[2]);
1205 charge = perigeeParameters->charge();
1206 } else {
1207 position.push_back(0);
1208 position.push_back(0);
1209 position.push_back(0);
1210 momentum.push_back(0);
1211 momentum.push_back(0);
1212 momentum.push_back(0);
1213 }
1214 int mot = 0;
1215 int oot = 0;
1216 if (measurementsOnTrack)
1217 mot = measurementsOnTrack->size();
1218 if (outliersOnTrack)
1219 oot = outliersOnTrack->size();
1220 std::vector<int> measurementsOnTrack_pixcl_sctcl_index, outliersOnTrack_pixcl_sctcl_index;
1221 int TTCindex, TTCevent_index, TTCparticle_link;
1222 float TTCprobability;
1223 if (measurementsOnTrack) {
1224 for (size_t i = 0; i < measurementsOnTrack->size(); i++) {
1225 const Trk::MeasurementBase *mb = (*measurementsOnTrack)[i];
1226 const InDet::PixelClusterOnTrack *pixcl = dynamic_cast<const InDet::PixelClusterOnTrack *>(mb);
1227 const InDet::SCT_ClusterOnTrack *sctcl = dynamic_cast<const InDet::SCT_ClusterOnTrack *>(mb);
1228 if (pixcl) {
1229 measurementsOnTrack_pixcl_sctcl_index.push_back(clusterIDMapIdx[pixcl->prepRawData()->identify()]);
1230 }
1231 else if (sctcl) {
1232 measurementsOnTrack_pixcl_sctcl_index.push_back(clusterIDMapIdx[sctcl->prepRawData()->identify()]);
1233 } else {
1234 measurementsOnTrack_pixcl_sctcl_index.push_back(-1);
1235 }
1236 }
1237 }
1238 if (outliersOnTrack) {
1239 for (size_t i = 0; i < outliersOnTrack->size(); i++) {
1240 const Trk::MeasurementBase *mb = (*outliersOnTrack)[i];
1241 const InDet::PixelClusterOnTrack *pixcl = dynamic_cast<const InDet::PixelClusterOnTrack *>(mb);
1242 const InDet::SCT_ClusterOnTrack *sctcl = dynamic_cast<const InDet::SCT_ClusterOnTrack *>(mb);
1243 if (pixcl) {
1244 outliersOnTrack_pixcl_sctcl_index.push_back(clusterIDMapIdx[pixcl->prepRawData()->identify()]);
1245 } else if (sctcl) {
1246 outliersOnTrack_pixcl_sctcl_index.push_back(clusterIDMapIdx[sctcl->prepRawData()->identify()]);
1247 } else {
1248 outliersOnTrack_pixcl_sctcl_index.push_back(-1);
1249 }
1250 }
1251 }
1252 if (found != trackTruthCollection->end()) {
1253 TTCindex = found->first.index();
1254 TTCevent_index = found->second.particleLink().eventIndex();
1255 TTCparticle_link = found->second.particleLink().barcode();
1256 TTCprobability = found->second.probability();
1257 } else {
1258 TTCindex = TTCevent_index = TTCparticle_link = -999;
1259 TTCprobability = -1;
1260 }
1261
1262 if (m_rootFile) {
1263 m_TRKindex[m_nTRK] = trk_index;
1264 m_TRKtrack_fitter[m_nTRK] = info.trackFitter();
1265 m_TRKndof[m_nTRK] = info.trackFitter();
1266 m_TRKparticle_hypothesis[m_nTRK] = info.particleHypothesis();
1267 (*m_TRKproperties).push_back(v_properties);
1268 (*m_TRKpattern).push_back(v_pattern);
1269 m_TRKndof[m_nTRK] = ndof;
1270 m_TRKchiSq[m_nTRK] = chiSq;
1271 (*m_TRKmeasurementsOnTrack_pixcl_sctcl_index).push_back(measurementsOnTrack_pixcl_sctcl_index);
1272 (*m_TRKoutliersOnTrack_pixcl_sctcl_index).push_back(outliersOnTrack_pixcl_sctcl_index);
1274 (*m_TRKperigee_position).push_back(position);
1275 (*m_TRKperigee_momentum).push_back(momentum);
1276 m_TRKmot[m_nTRK] = mot;
1277 m_TRKoot[m_nTRK] = oot;
1278 m_TTCindex[m_nTRK] = TTCindex;
1279 m_TTCevent_index[m_nTRK] = TTCevent_index;
1280 m_TTCparticle_link[m_nTRK] = TTCparticle_link;
1281 m_TTCprobability[m_nTRK] = TTCprobability;
1282 }
1283
1284 trk_index++;
1285 // index
1286 m_nTRK++;
1287 if (m_nTRK == m_maxTRK) {
1288 ATH_MSG_WARNING("DUMP : hit max number of track events");
1289 break;
1290 }
1291 }
1292
1293 const DetailedTrackTruthCollection *detailedTrackTruthCollection = 0;
1294 SG::ReadHandle<DetailedTrackTruthCollection> detailedTrackTruthCollectionHandle{m_detailedTracksTruthKey, ctx};
1295 if (not detailedTrackTruthCollectionHandle.isValid()) {
1296 ATH_MSG_WARNING(" DetailedTrackTruthCollection not found: " << m_detailedTracksTruthKey.key());
1297 return StatusCode::FAILURE;
1298 }
1299 detailedTrackTruthCollection = detailedTrackTruthCollectionHandle.cptr();
1300
1301 m_nDTT = 0;
1302 if (m_rootFile) {
1303 (*m_DTTtrajectory_eventindex).clear();
1304 (*m_DTTtrajectory_barcode).clear();
1305 (*m_DTTstTruth_subDetType).clear();
1306 (*m_DTTstTrack_subDetType).clear();
1307 (*m_DTTstCommon_subDetType).clear();
1308 }
1309
1310 // loop over DetailedTrackTruth objects
1311 DetailedTrackTruthCollection::const_iterator detailedTrackTruthIterator = (*detailedTrackTruthCollection).begin();
1312 for (; detailedTrackTruthIterator != (*detailedTrackTruthCollection).end(); ++detailedTrackTruthIterator) {
1313 std::vector<int> DTTtrajectory_eventindex, DTTtrajectory_barcode, DTTstTruth_subDetType, DTTstTrack_subDetType,
1314 DTTstCommon_subDetType;
1315 const TruthTrajectory &traj = detailedTrackTruthIterator->second.trajectory();
1316 for (size_t j = 0; j < traj.size(); j++) {
1317 DTTtrajectory_eventindex.push_back(traj[j].eventIndex());
1318 DTTtrajectory_barcode.push_back(traj[j].barcode());
1319 }
1320 const SubDetHitStatistics &stTruth = detailedTrackTruthIterator->second.statsTruth();
1321 const SubDetHitStatistics &stTrack = detailedTrackTruthIterator->second.statsTrack();
1322 const SubDetHitStatistics &stCommon = detailedTrackTruthIterator->second.statsCommon();
1323 for (unsigned j = 0; j < SubDetHitStatistics::NUM_SUBDETECTORS; j++) {
1324 DTTstTruth_subDetType.push_back(stTruth[SubDetHitStatistics::SubDetType(j)]);
1325 }
1326 for (unsigned j = 0; j < SubDetHitStatistics::NUM_SUBDETECTORS; j++) {
1327 DTTstTrack_subDetType.push_back(stTrack[SubDetHitStatistics::SubDetType(j)]);
1328 }
1329 for (unsigned j = 0; j < SubDetHitStatistics::NUM_SUBDETECTORS; j++) {
1330 DTTstCommon_subDetType.push_back(stCommon[SubDetHitStatistics::SubDetType(j)]);
1331 }
1332
1333 if (m_rootFile) {
1334 m_DTTindex[m_nDTT] = detailedTrackTruthIterator->first.index();
1335 m_DTTsize[m_nDTT] = traj.size();
1336 (*m_DTTtrajectory_eventindex).push_back(DTTtrajectory_eventindex);
1337 (*m_DTTtrajectory_barcode).push_back(DTTtrajectory_barcode);
1338 (*m_DTTstTruth_subDetType).push_back(DTTstTruth_subDetType);
1339 (*m_DTTstTrack_subDetType).push_back(DTTstTrack_subDetType);
1340 (*m_DTTstCommon_subDetType).push_back(DTTstCommon_subDetType);
1341 }
1342
1343 m_nDTT++;
1344 }
1345
1346 // Once all the information for this event has been filled in the arrays,
1347 // copy content of the arrays to the TTree
1348 if (m_rootFile)
1349 m_nt->Fill();
1350
1351 return StatusCode::SUCCESS;
1352}
1353
1354//--------------------------------
1356 //--------------------------------
1357 if (m_rootFile) {
1358 delete[] m_SEID;
1359
1360 delete[] m_CLindex;
1361 delete m_CLhardware;
1362 delete[] m_CLx;
1363 delete[] m_CLy;
1364 delete[] m_CLz;
1365 delete[] m_CLbarrel_endcap;
1366 delete[] m_CLlayer_disk;
1367 delete[] m_CLeta_module;
1368 delete[] m_CLphi_module;
1369 delete[] m_CLside;
1370 delete[] m_CLmoduleID;
1373 delete m_CLbarcodesLinked;
1374 delete m_CLparticle_charge;
1375 delete m_CLphis;
1376 delete m_CLetas;
1377 delete m_CLtots;
1378 delete[] m_CLloc_direction1;
1379 delete[] m_CLloc_direction2;
1380 delete[] m_CLloc_direction3;
1381 delete[] m_CLJan_loc_direction1;
1382 delete[] m_CLJan_loc_direction2;
1383 delete[] m_CLJan_loc_direction3;
1384 delete[] m_CLpixel_count;
1385 delete[] m_CLcharge_count;
1386 delete[] m_CLloc_eta;
1387 delete[] m_CLloc_phi;
1388 delete[] m_CLglob_eta;
1389 delete[] m_CLglob_phi;
1390 delete[] m_CLeta_angle;
1391 delete[] m_CLphi_angle;
1392 delete[] m_CLnorm_x;
1393 delete[] m_CLnorm_y;
1394 delete[] m_CLnorm_z;
1395 delete m_CLlocal_cov;
1396
1397 delete[] m_Part_event_number;
1398 delete[] m_Part_barcode;
1399 delete[] m_Part_px;
1400 delete[] m_Part_py;
1401 delete[] m_Part_pz;
1402 delete[] m_Part_pt;
1403 delete[] m_Part_eta;
1404 delete[] m_Part_vx;
1405 delete[] m_Part_vy;
1406 delete[] m_Part_vz;
1407 delete[] m_Part_radius;
1408 delete[] m_Part_status;
1409 delete[] m_Part_charge;
1410 delete[] m_Part_pdg_id;
1411 delete[] m_Part_passed;
1412
1413 delete[] m_Part_vProdNin;
1414 delete[] m_Part_vProdNout;
1415 delete[] m_Part_vProdStatus;
1416 delete[] m_Part_vProdBarcode;
1417 delete m_Part_vParentID;
1418 delete m_Part_vParentBarcode;
1419
1420 delete[] m_SPindex;
1421 delete[] m_SPx;
1422 delete[] m_SPy;
1423 delete[] m_SPz;
1424 delete[] m_SPCL1_index;
1425 delete[] m_SPCL2_index;
1426 delete[] m_SPisOverlap;
1427 delete[] m_SPradius;
1428 delete[] m_SPcovr;
1429 delete[] m_SPcovz;
1430 delete[] m_SPhl_topstrip;
1431 delete[] m_SPhl_botstrip;
1432 delete m_SPtopStripDirection;
1436
1437 delete[] m_TRKindex;
1438 delete[] m_TRKtrack_fitter;
1439 delete[] m_TRKparticle_hypothesis;
1440 delete m_TRKproperties;
1441 delete m_TRKpattern;
1442 delete[] m_TRKndof;
1443 delete[] m_TRKmot;
1444 delete[] m_TRKoot;
1445 delete[] m_TRKchiSq;
1448 delete[] m_TRKcharge;
1449 delete m_TRKperigee_position;
1450 delete m_TRKperigee_momentum;
1451 delete[] m_TTCindex;
1452 delete[] m_TTCevent_index;
1453 delete[] m_TTCparticle_link;
1454 delete[] m_TTCprobability;
1455
1456 delete[] m_DTTindex;
1457 delete[] m_DTTsize;
1463 }
1464
1465 return StatusCode::SUCCESS;
1466}
1467
1468//--------------------------------------------------------------------------------------------
1469bool InDet::DumpObjects::isPassed(HepMC::ConstGenParticlePtr particle, float &px, float &py, float &pz,
1470 float &pt, float &eta, float &vx, float &vy, float &vz, float &radius, float &status,
1471 float &charge, std::vector<int> &vParentID, std::vector<int> &vParentBarcode,
1472 int &vProdNin, int &vProdNout, int &vProdStatus, int &vProdBarcode) {
1473 //--------------------------------------------------------------------------------------------
1474
1475 px = particle->momentum().px();
1476 py = particle->momentum().py();
1477 pz = particle->momentum().pz();
1478
1479 pt = std::sqrt(px * px + py * py);
1480 eta = particle->momentum().eta();
1481
1482 int pdgCode = particle->pdg_id();
1483
1484 int absPdgCode = std::abs(pdgCode);
1485 // get the charge: ap->charge() is used later, DOES NOT WORK RIGHT NOW
1486 const HepPDT::ParticleData *ap = m_particleDataTable->particle(absPdgCode);
1487 charge = 1.;
1488 if (ap)
1489 charge = ap->charge();
1490 // since the PDT table only has abs(PID) values for the charge
1491 charge *= (pdgCode > 0.) ? 1. : -1.;
1492
1493 status = particle->status();
1494
1495 if (particle->production_vertex()) {
1496 vx = particle->production_vertex()->position().x();
1497 vy = particle->production_vertex()->position().y();
1498 vz = particle->production_vertex()->position().z();
1499 radius = particle->production_vertex()->position().perp();
1500 } else {
1501 vx = vy = vz = -1;
1502 radius = 999;
1503 if (status == 1)
1504 ATH_MSG_WARNING("no vertex for particle with status 1");
1505 }
1506
1507 if (particle->production_vertex()) {
1508 vProdNin = particle->production_vertex()->particles_in_size();
1509 vProdNout = particle->production_vertex()->particles_out_size();
1510 vProdStatus = particle->production_vertex()->id();
1511 vProdBarcode = HepMC::barcode(particle->production_vertex()); // JB: HEPMC3 barcode() -> HepMC::barcode(p)
1512#ifdef HEPMC3
1513 for (const auto &p : particle->production_vertex()->particles_in()) {
1514 vParentID.push_back(p->pdg_id());
1515 vParentBarcode.push_back(HepMC::barcode(p)); // JB: HEPMC3 barcode() -> HepMC::barcode(p)
1516 }
1517#else
1518 for (auto ip = particle->production_vertex()->particles_in_const_begin();
1519 ip != particle->production_vertex()->particles_in_const_end();
1520 ++ip)
1521 {
1522 vParentID.push_back((*ip)->pdg_id());
1523 vParentBarcode.push_back(HepMC::barcode(*ip)); // JB: HEPMC3 barcode() -> HepMC::barcode(p)
1524 }
1525#endif
1526 } else {
1527 vProdNin = 0;
1528 vProdNout = 0;
1529 vProdStatus = -999;
1530 vProdBarcode = 999;
1531 }
1532
1533 bool passEta = (pt > 0.1) ? (std::abs(eta) < m_max_eta) : false;
1534 if (not passEta)
1535 return false;
1536
1537 bool passPt = (pt > m_min_pt);
1538 if (not passPt)
1539 return false;
1540
1541 bool passBarcode = (HepMC::barcode(particle) < m_max_barcode); // JB: HEPMC3 barcode() -> HepMC::barcode(p)
1542 if (not passBarcode)
1543 return false;
1544
1545 bool passCharge = not(charge == 0.);
1546 if (not passCharge)
1547 return false;
1548
1549 bool passStatus = (status == 1);
1550 if (not passStatus)
1551 return false;
1552
1553 bool passProdRadius = (radius < m_maxProdVertex);
1554 if (not passProdRadius)
1555 return false;
1556
1557 return true;
1558}
Scalar eta() const
pseudorapidity method
Scalar phi() const
phi method
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_ERROR(x)
#define ATH_MSG_FATAL(x)
#define ATH_MSG_WARNING(x)
double charge(const T &p)
Definition AtlasPID.h:997
std::vector< Identifier > ID
bool passed(DecisionID id, const DecisionIDContainer &)
checks if required decision ID is in the set of IDs in the container
static Double_t sp
static Double_t sc
DataVector< Trk::Track > TrackCollection
This typedef represents a collection of Trk::Track objects.
AthAlgorithm(const std::string &name, ISvcLocator *pSvcLocator)
Constructor with parameters:
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
const ServiceHandle< StoreGateSvc > & detStore() const
Derived DataVector<T>.
Definition DataVector.h:795
DataModel_detail::const_iterator< DataVector > const_iterator
Definition DataVector.h:838
const T * at(size_type n) const
Access an element, as an rvalue.
size_type size() const noexcept
Returns the number of elements in the collection.
bool is_valid() const
Check if id is in a valid state.
Class used to describe the design of a module (diode segmentation and readout scheme)
virtual SiLocalPosition localPositionOfCell(const SiCellId &cellId) const
readout or diode id -> position.
Base class for the SCT module side design, extended by the Forward and Barrel module design.
virtual SiLocalPosition localPositionOfCell(const SiCellId &cellId) const override=0
id -> position
Identifier for the strip or pixel cell.
Definition SiCellId.h:29
Class to hold geometrical description of a silicon detector element.
virtual const SiDetectorDesign & design() const override final
access to the local description (inline):
double phiPitch() const
Pitch (inline methods)
std::pair< Amg::Vector3D, Amg::Vector3D > endsOfStrip(const Amg::Vector2D &position) const
Special method for SCT to retrieve the two ends of a "strip" Returned coordinates are in global frame...
Class to represent a position in the natural frame of a silicon sensor, for Pixel and SCT For Pixel: ...
double xPhi() const
position along phi direction:
double xEta() const
position along eta direction:
virtual const Amg::Vector3D & normal() const override final
Get reconstruction local normal axes in global frame.
HepGeom::Point3D< double > globalPosition(const HepGeom::Point3D< double > &localPos) const
transform a reconstruction local position into a global position (inline):
double etaPitch() const
Pitch (inline methods)
double * m_CLloc_direction3
virtual StatusCode execute() override final
std::vector< std::vector< int > > * m_CLtots
std::vector< std::vector< int > > * m_DTTtrajectory_eventindex
int m_maxCL
jobOption: maximum number of clusters
std::vector< std::vector< bool > > * m_CLbarcodesLinked
std::vector< std::vector< int > > * m_TRKpattern
std::vector< std::vector< int > > * m_DTTstTruth_subDetType
SG::ReadHandleKey< DetailedTrackTruthCollection > m_detailedTracksTruthKey
Definition DumpObjects.h:96
bool isPassed(HepMC::ConstGenParticlePtr particle, float &px, float &py, float &pz, float &pt, float &eta, float &vx, float &vy, float &vz, float &radius, float &status, float &charge, std::vector< int > &vParentID, std::vector< int > &vParentBarcode, int &vProdNin, int &vProdNout, int &vProdStatus, int &vProdBarcode)
std::vector< std::vector< int > > * m_CLparticleLink_barcode
std::vector< std::vector< int > > * m_TRKoutliersOnTrack_pixcl_sctcl_index
unsigned int m_run_number
double * m_CLloc_direction2
SG::ReadHandleKey< InDet::SCT_ClusterContainer > m_stripClusterKey
Definition DumpObjects.h:83
int m_maxSP
jobOption: maximum number of space points
const SCT_ID * m_SCT_ID
Definition DumpObjects.h:69
SG::ReadHandleKey< TrackTruthCollection > m_tracksTruthKey
Definition DumpObjects.h:95
std::vector< std::vector< double > > * m_TRKperigee_momentum
SG::ReadHandleKey< TrackCollection > m_tracksKey
Definition DumpObjects.h:94
unsigned long long m_event_number
double * m_CLJan_loc_direction3
SG::ReadHandleKey< xAOD::EventInfo > m_eventInfoKey
Definition DumpObjects.h:79
std::vector< std::vector< float > > * m_SPtopStripDirection
double * m_CLJan_loc_direction2
std::string m_name
Definition DumpObjects.h:76
uint64_t * m_CLmoduleID
int * m_TRKparticle_hypothesis
std::vector< std::vector< float > > * m_SPstripCenterDistance
const InDetDD::SCT_DetectorManager * m_SCT_Manager
Definition DumpObjects.h:71
std::vector< std::vector< double > > * m_TRKperigee_position
std::vector< std::vector< int > > * m_Part_vParentID
const PixelID * m_pixelID
Definition DumpObjects.h:68
std::string m_ntupleDirName
jobOption: Ntuple directory name
std::string m_ntupleFileName
jobOption: Ntuple file name
double * m_CLloc_direction1
std::vector< std::vector< int > > * m_DTTstTrack_subDetType
bool m_rootFile
jobOption: save data in root format
double * m_CLJan_loc_direction1
std::vector< std::vector< int > > * m_CLphis
ServiceHandle< IPartPropSvc > m_particlePropSvc
Definition DumpObjects.h:74
virtual StatusCode finalize() override final
SG::ReadHandleKey< McEventCollection > m_mcEventCollectionKey
Definition DumpObjects.h:80
std::vector< std::vector< int > > * m_DTTtrajectory_barcode
SG::ReadHandleKey< xAOD::SpacePointContainer > m_xaodPixelSpacePointContainerKey
Definition DumpObjects.h:88
SG::ReadHandleKey< InDetSimDataCollection > m_stripSDOKey
Definition DumpObjects.h:86
const InDetDD::PixelDetectorManager * m_pixelManager
Definition DumpObjects.h:70
DumpObjects(const std::string &name, ISvcLocator *pSvcLocator)
virtual StatusCode initialize() override final
std::string m_ntupleTreeName
jobOption: Ntuple tree name
std::vector< std::vector< int > > * m_TRKmeasurementsOnTrack_pixcl_sctcl_index
std::vector< std::vector< int > > * m_DTTstCommon_subDetType
std::vector< std::vector< int > > * m_TRKproperties
std::vector< std::string > * m_CLhardware
std::vector< std::vector< int > > * m_CLparticleLink_eventIndex
SG::ReadHandleKey< xAOD::SpacePointContainer > m_xaodStripSpacePointContainerKey
Definition DumpObjects.h:90
const HepPDT::ParticleDataTable * m_particleDataTable
Definition DumpObjects.h:75
SG::ReadHandleKey< InDetSimDataCollection > m_pixelSDOKey
Definition DumpObjects.h:85
SG::ReadHandleKey< InDet::PixelClusterContainer > m_pixelClusterKey
Definition DumpObjects.h:82
int m_maxPart
jobOption: maximum number of particles
std::vector< std::vector< int > > * m_CLetas
std::vector< std::vector< float > > * m_CLparticle_charge
std::vector< std::vector< float > > * m_SPtopStripCenterPosition
std::vector< std::vector< double > > * m_CLlocal_cov
SG::ReadHandleKey< xAOD::SpacePointContainer > m_xaodStripSpacePointOverlapContainerKey
Definition DumpObjects.h:92
std::vector< std::vector< int > > * m_Part_vParentBarcode
std::vector< std::vector< float > > * m_SPbottomStripDirection
Specific class to represent the pixel measurements.
virtual const PixelCluster * prepRawData() const override final
returns the PrepRawData - is a SiCluster in this scope
Specific class to represent the SCT measurements.
virtual const InDet::SCT_Cluster * prepRawData() const override final
returns the PrepRawData - is a SCT_Cluster in this scope
This defines the McEventCollection, which is really just an ObjectVector of McEvent objectsFile: Gene...
Helper class to provide type-safe access to aux data.
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
virtual bool isValid() override final
Can the handle be successfully dereferenced?
const_pointer_type cptr()
Dereference the pointer.
int numberDoF() const
returns the number of degrees of freedom of the overall track or vertex fit as integer
Definition FitQuality.h:60
double chiSquared() const
returns the of the overall track fit
Definition FitQuality.h:56
Class to represent and store fit qualities from track reconstruction in terms of and number of degre...
Definition FitQuality.h:97
This class is the pure abstract base class for all fittable tracking measurements.
Identifier identify() const
return the identifier
Contains information about the 'fitter' of this track.
A TruthTrajectory is a chain of charged MC particles connected through the mother-daughter relationsh...
uint32_t runNumber() const
The current event's run number.
uint64_t eventNumber() const
The current event's event number.
bool contains(const std::string &s, const std::string &regx)
does a string contain the substring
Definition hcg.cxx:114
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Eigen::Matrix< double, 2, 1 > Vector2D
Eigen::Matrix< double, 3, 1 > Vector3D
int barcode(const T *p)
Definition Barcode.h:16
const GenParticle * ConstGenParticlePtr
Definition GenParticle.h:38
int compute_overlap_SP_flag(const int &eta_module_cl1, const int &phi_module_cl1, const int &eta_module_cl2, const int &phi_module_cl2)
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
bool flag
Definition master.py:29
EventInfo_v1 EventInfo
Definition of the latest event info version.
SpacePointContainer_v1 SpacePointContainer
Define the version of the space point container.
perigeeParameters(double d0, double z0, double phi, double eta, double pt, double charge)