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