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