ATLAS Offline Software
Loading...
Searching...
No Matches
SingleTrackValidation.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
3*/
4
5#include "SingleTrackValidation.h"
6#include "StoreGate/StoreGateSvc.h"
7#include "StoreGate/ReadCondHandle.h"
8
9// In order to be able to Ntuple & cet:
10#include "GaudiKernel/NTuple.h"
11#include "GaudiKernel/INTupleSvc.h"
12#include "GaudiKernel/SmartDataPtr.h"
13
14// To get the Magnetic Field:
15#include "CLHEP/GenericFunctions/FixedConstant.hh"
16
17// TruthUtils
18#include "TruthUtils/HepMCHelpers.h"
19
20
21// To extrapolate:
22#include "AtlasBComponent.h"
23#include "GeoXPEngine.h"
24
25// To interpret LAr Geometry information:
26#include "CaloDetDescr/CaloDetDescrElement.h"
27#include "LArSimEvent/LArHitContainer.h"
28#include "CaloIdentifier/CaloCell_ID.h"
29#include "CaloIdentifier/LArID_Exception.h"
30
31// pi etc
32#include "CLHEP/Units/PhysicalConstants.h"
33#include "AthenaKernel/Units.h"
34
35// System libraries
36#include <signal.h>
37#include <unistd.h>
38#include <sys/types.h>
39#include <sstream>
40#include <fstream>
41#include <iostream>
42#include <sys/times.h>
43#include <string>
44
45namespace Units = Athena::Units;
46using HepGeom::Point3D;
47using CLHEP::HepLorentzVector;
48using CLHEP::mm;
49using CLHEP::pi;
50using CLHEP::twopi;
51
52inline int iabs (int a) { return (a>=0)? a : -a;}
53inline int getCpu() {
54 int cpuTime = 0;
55 {
56 static const int ticksPerJiffy = sysconf(_SC_CLK_TCK)/100;
57 struct tms buff{};
58 times(&buff);
59 cpuTime=(buff.tms_utime + buff.tms_stime + buff.tms_cutime + buff.tms_cstime)/ticksPerJiffy;
60 }
61 return cpuTime;
62}
63
64// We will define some actual data, later.
65class SingleTrackValidation::Clockwork {
66
67public:
68 Clockwork() {}
69 ~Clockwork() {}
70
71 IPartPropSvc* partPropSvc{nullptr};
72 ITHistSvc* histSvc{nullptr};
73 const CaloCell_ID* cellId{nullptr};
74 NTuple::Tuple* nt = nullptr;
75
76 NTuple::Item<double> eta;
77 NTuple::Item<double> pt;
78 NTuple::Item<double> phi;
79
80 // Impact point
81 NTuple::Item<double> x;
82 NTuple::Item<double> y;
83 NTuple::Item<double> z;
84
85 // Indices run from 0 to 11: 4 in barrel, 4 in end cap, 3 in FCAL
86 // Energy within each sampling layer (directly struck cell):
87 NTuple::Item<double> s_c00[15];
88 NTuple::Item<double> s_hits[15];
89 NTuple::Item<double> s_sumE[15];
90 NTuple::Item<double> s_deltaPhi[15];
91 NTuple::Item<double> s_sigmaPhi[15];
92 NTuple::Item<double> s_deltaEta[15];
93 NTuple::Item<double> s_sigmaEta[15];
94 NTuple::Item<double> s_t00[15];
95 NTuple::Item<double> s_widthX[15];
96 NTuple::Item<double> s_widthY[15];
97
98 NTuple::Item<double> cpuTime;
99 NTuple::Item<double> Energy;
100 NTuple::Item<double> PDG;
101 NTuple::Item<double> RunNo;
102 NTuple::Item<double> EventNo;
103 NTuple::Item<double> E_Deposit;
104
105 bool coolField{true};
106};
107
108SingleTrackValidation::SingleTrackValidation (const std::string & name, ISvcLocator * pSvcLocator) :
109 AthAlgorithm(name,pSvcLocator),m_c(new Clockwork())
110{
111 for (unsigned int i=0;i<162;++i) m_histos[i] = nullptr;
112}
113
114SingleTrackValidation::~SingleTrackValidation () {
115 if (m_c!=nullptr){ delete m_c; m_c=nullptr; }
116}
117
118StatusCode SingleTrackValidation::initialize() {
119 std::string names[162] = { "eta", "pt", "phi", "pos_x", "pos_y", "pos_z",
120 "emb0_cell", "emb1_cell", "emb2_cell", "emb3_cell", "emec0_cell", "emec1_cell", "emec2_cell", "emec3_cell",
121 "hec0_cell", "hec1_cell", "hec2_cell", "hec3_cell", "fc1_cell", "fc2_cell", "fc3_cell",
122 "emb0_hits", "emb1_hits", "emb2_hits", "emb3_hits", "emec0_hits", "emec1_hits", "emec2_hits", "emec3_hits",
123 "hec0_hits", "hec1_hits", "hec2_hits", "hec3_hits", "fc1_hits", "fc2_hits", "fc3_hits",
124 "emb0_sumE", "emb1_sumE", "emb2_sumE", "emb3_sumE", "emec0_sumE", "emec1_sumE", "emec2_sumE", "emec3_sumE",
125 "hec0_sumE", "hec1_sumE", "hec2_sumE", "hec3_sumE", "fc1_sumE", "fc2_sumE", "fc3_sumE",
126 "emb0_dPhi", "emb1_dPhi", "emb2_dPhi", "emb3_dPhi", "emec0_dPhi", "emec1_dPhi", "emec2_dPhi", "emec3_dPhi",
127 "hec0_dPhi", "hec1_dPhi", "hec2_dPhi", "hec3_dPhi", "fc1_dX", "fc2_dX", "fc3_dX",
128 "emb0_sPhi", "emb1_sPhi", "emb2_sPhi", "emb3_sPhi", "emec0_sPhi", "emec1_sPhi", "emec2_sPhi", "emec3_sPhi",
129 "hec0_sPhi", "hec1_sPhi", "hec2_sPhi", "hec3_sPhi", "fc1_sX", "fc2_sX", "fc3_sX",
130 "emb0_dEta", "emb1_dEta", "emb2_dEta", "emb3_dEta", "emec0_dEta", "emec1_dEta", "emec2_dEta", "emec3_dEta",
131 "hec0_dEta", "hec1_dEta", "hec2_dEta", "hec3_dEta", "fc1_dY", "fc2_dY", "fc3_dY",
132 "emb0_sEta", "emb1_sEta", "emb2_sEta", "emb3_sEta", "emec0_sEta", "emec1_sEta", "emec2_sEta", "emec3_sEta",
133 "hec0_sEta", "hec1_sEta", "hec2_sEta", "hec3_sEta", "fc1_sY", "fc2_sY", "fc3_sY",
134 "emb0_time", "emb1_time", "emb2_time", "emb3_time", "emec0_time", "emec1_time", "emec2_time", "emec3_time",
135 "hec0_time", "hec1_time", "hec2_time", "hec3_time", "fc1_time", "fc2_time", "fc3_time",
136 "emb0_widthX", "emb1_widthX", "emb2_widthX", "emb3_widthX", "emec0_widthX", "emec1_widthX", "emec2_widthX", "emec3_widthX",
137 "hec0_widthX", "hec1_widthX", "hec2_widthX", "hec3_widthX", "fc1_widthX", "fc2_widthX", "fc3_widthX",
138 "emb0_widthY", "emb1_widthY", "emb2_widthY", "emb3_widthY", "emec0_widthY", "emec1_widthY", "emec2_widthY", "emec3_widthY",
139 "hec0_widthY", "hec1_widthY", "hec2_widthY", "hec3_widthY", "fc1_widthY", "fc2_widthY", "fc3_widthY",
140 "cpuTime", "Energy", "PDG_ID", "RunNo", "EventNo", "E_Dep" };
141
142 double lim[162][2] = { {0,5}, {0,100}, {-4,4}, {-1600,1600}, {-1600,1600}, {-4000,4000},
143 {0,0.25}, {0,3}, {0,6}, {0,0.1}, {0,0.25}, {0,2}, {0,3}, {0,0.1}, {0,1}, {0,10}, {0,10}, {0,10}, {0,0.1}, {0,0.1}, {0,0.1}, //cell
144 {0,100}, {0,600}, {0,600}, {0,50}, {0,50}, {0,400}, {0,500}, {0,200}, {0,100}, {0,1000}, {0,1000}, {0,100}, {0,150}, {0,50}, {0,10}, //hits
145 {0,0.4}, {0,5}, {0,10}, {0,0.2}, {0,0.1}, {0,3}, {0,5}, {0,0.2}, {0,1}, {0,10}, {0,10}, {0,0.1}, {0,1}, {0,0.1}, {0,0.1}, //sumE
146 {-500,500}, {-100,100}, {-15,15}, {-200,200}, {-3000,3000}, {-60,60}, {-25,25}, {-200,200}, {-50,50}, {-50,50}, {-50,50}, {-50,50}, {-60,60}, {-500,500}, {-200,200}, //dPhi
147 {0,1000}, {0,500}, {0,200}, {0,500}, {0,2000}, {0,250}, {0,300}, {0,500}, {0,200}, {0,200}, {0,200}, {0,200}, {0,100}, {0,100}, {0,50}, //sPhi
148 {-150,150}, {-15,15}, {-20,20}, {-200,200}, {-0,2500}, {-50,20}, {-15,15}, {-150,150}, {-50,50}, {-50,50}, {-50,50}, {-50,50}, {-60,60}, {-500,500}, {-200,200}, //dEta
149 {0,500}, {0,100}, {0,100}, {0,400}, {0,1000}, {0,150}, {0,100}, {0,400}, {0,200}, {0,200}, {0,200}, {0,200}, {0,60}, {0,100}, {0,50}, //sPhi
150 {0,750}, {0,25}, {0,20}, {0,1000}, {0,10000}, {0,40}, {0,30}, {0,1000}, {0,1000}, {0,100}, {0,100}, {0,200}, {0,10}, {0,500}, {0,100}, //time
151 {-150,150}, {-15,15}, {-20,20}, {-200,200}, {-0,2500}, {-50,20}, {-15,15}, {-150,150}, {-50,50}, {-50,50}, {-50,50}, {-50,50}, {-60,60}, {-500,500}, {-200,200}, //widthX
152 {-150,150}, {-15,15}, {-20,20}, {-200,200}, {-0,2500}, {-50,20}, {-15,15}, {-150,150}, {-50,50}, {-50,50}, {-50,50}, {-50,50}, {-60,60}, {-500,500}, {-200,200}, //widthY
153 {0,50}, {0,100}, {-25,25}, {0,10}, {0,1000}, {0,10}};
154
155 //-------------------------------------------------------------------------//
156 // //
157 // Initialize the Particle Property Service. This is necessary in order //
158 // to obtain charge & type & other properties of the primary particle and //
159 // other particles that may turn up in the debris. //
160 // //
161 ATH_CHECK(m_ppSvc.retrieve());
162 m_c->partPropSvc = m_ppSvc.get();
163 ATH_CHECK(m_histSvc.retrieve());
164 m_c->histSvc = m_histSvc.get();
165 ATH_CHECK(detStore()->retrieve(m_c->cellId, "CaloCell_ID"));
166 ATH_CHECK(m_truthKey.initialize());
167 ATH_CHECK(m_caloMgrKey.initialize());
168 ATH_CHECK(m_fieldCacheCondObjInputKey.initialize());
169
170 //----------------Now initialize the ntuples ----------------------//
171 // //
172 //==~ ~ ~==//
173 std::string filename="";
174 for (int i=0;i<162;i++){
175 m_histos[i] = new TH1F( names[i].data(), names[i].data(),50,lim[i][0],lim[i][1]);
176 filename = "/file1/Electron/";
177 filename.append(names[i]);
178 if (m_c->histSvc->regHist( filename.data() , m_histos[i] ).isFailure()){
179 ATH_MSG_WARNING( "Failed to register historam " << names[i] << ". Not sure what will happen now..." );
180 }
181 }
182
183
184 NTupleFilePtr file(ntupleSvc(),"/NTUPLES/FILE");
185 if (!file) throw std::runtime_error ("Ntuple MGR not open");
186 NTuple::Directory *col=ntupleSvc()->createDirectory("/NTUPLES/FILE/COL");
187 NTuplePtr nt(ntupleSvc(),"/NTUPLES/FILE/COL/SingleTrackValidation");
188 if (!nt) nt=ntupleSvc()->book(col, 1, CLID_ColumnWiseTuple, "SingleTrackValidation");
189
190 if (nt->addItem("Eta", m_c->eta ).isFailure() ||
191 nt->addItem("Pt", m_c->pt ).isFailure() ||
192 nt->addItem("BarrelX", m_c->x ).isFailure() ||
193 nt->addItem("BarrelY", m_c->y ).isFailure() ||
194 nt->addItem("BarrelZ", m_c->z ).isFailure() ||
195 nt->addItem("Phi", m_c->phi ).isFailure() ){
196 ATH_MSG_WARNING( "Registration of some of the ntuple branches failed. No idea what will happen next..." );
197 }
198
199
200 char title[80];
201
202 // Indices from 0 to 11: 4 EMB, 4 EMEC, and 4 FCAL sampling layers
203 // Handling FCAL layers separately for different titles (could do
204 // something more complicated if we wanted)
205 for (int i=0;i<15;i++){
206 if (i<12) sprintf(title,"S%i_C00",i);
207 else sprintf(title,"FC%i_C00",i-11);
208 if (nt->addItem(title,m_c->s_c00[i]).isFailure()) ATH_MSG_INFO( "Registration of a branch failed in the ntupler..." );
209 if (i<12) sprintf(title,"S%i_SumE",i);
210 else sprintf(title,"FC%i_SumE",i-11);
211 if (nt->addItem(title,m_c->s_sumE[i]).isFailure()) ATH_MSG_INFO( "Registration of a branch failed in the ntupler..." );
212 if (i<12) sprintf(title,"S%i_Hits",i);
213 else sprintf(title,"FC%i_Hits",i-11);
214 if (nt->addItem(title,m_c->s_hits[i]).isFailure()) ATH_MSG_INFO( "Registration of a branch failed in the ntupler..." );
215 if (i<12) sprintf(title,"S%i_DeltaPhi",i);
216 else sprintf(title,"FC%i_DeltaX",i-11);
217 if (nt->addItem(title,m_c->s_deltaPhi[i]).isFailure()) ATH_MSG_INFO( "Registration of a branch failed in the ntupler..." );
218 if (i<12) sprintf(title,"S%i_SigmaPhi",i);
219 else sprintf(title,"FC%i_SigmaX",i-11);
220 if (nt->addItem(title,m_c->s_sigmaPhi[i]).isFailure()) ATH_MSG_INFO( "Registration of a branch failed in the ntupler..." );
221 if (i<12) sprintf(title,"S%i_DeltaEta",i);
222 else sprintf(title,"FC%i_DeltaY",i-11);
223 if (nt->addItem(title,m_c->s_deltaEta[i]).isFailure()) ATH_MSG_INFO( "Registration of a branch failed in the ntupler..." );
224 if (i<12) sprintf(title,"S%i_SigmaEta",i);
225 else sprintf(title,"FC%i_SigmaY",i-11);
226 if (nt->addItem(title,m_c->s_sigmaEta[i]).isFailure()) ATH_MSG_INFO( "Registration of a branch failed in the ntupler..." );
227 if (i<12) sprintf(title,"S%i_Time",i);
228 else sprintf(title,"FC%i_Time",i-11);
229 if (nt->addItem(title,m_c->s_t00[i]).isFailure()) ATH_MSG_INFO( "Registration of a branch failed in the ntupler..." );
230 if (i<12) sprintf(title,"S%i_WidthX",i);
231 else sprintf(title,"FC%i_WidthX",i-11);
232 if (nt->addItem(title,m_c->s_widthX[i]).isFailure()) ATH_MSG_INFO( "Registration of a branch failed in the ntupler..." );
233 if (i<12) sprintf(title,"S%i_WidthY",i);
234 else sprintf(title,"FC%i_WidthY",i-11);
235 if (nt->addItem(title,m_c->s_widthY[i]).isFailure()) ATH_MSG_INFO( "Registration of a branch failed in the ntupler..." );
236 }
237
238 if (nt->addItem("CPU" , m_c->cpuTime ).isFailure() ||
239 nt->addItem("TrackEnergy" , m_c->Energy ).isFailure() ||
240 nt->addItem("ParticleID" , m_c->PDG ).isFailure() ||
241 nt->addItem("Run#" , m_c->RunNo ).isFailure() ||
242 nt->addItem("Event#" , m_c->EventNo ).isFailure() ||
243 nt->addItem("DepositedEnergy", m_c->E_Deposit ).isFailure() ){
244 ATH_MSG_WARNING( "Registration of some of the ntuple branches failed. No idea what will happen next..." );
245 }
246
247 m_c->cpuTime=0.0;
248 m_c->nt = nt;
249
250 //==~ ~ ~==//
251 // //
252 //------------------------Done with initializations------------------------//
253
254 return StatusCode::SUCCESS;
255}
256
257StatusCode SingleTrackValidation::execute(const EventContext& ctx) {
258
259 if (m_c->cpuTime==0) {
260 m_c->cpuTime=getCpu();
261 m_c->cpuTime+=1; // Fill the histogram with -1 for the first event
262 }
263 m_c->cpuTime= getCpu()-m_c->cpuTime;
264 m_histos[156]->Fill( m_c->cpuTime/100. , 1. );
265
266 int RunNum=ctx.eventID().run_number();
267 int EvtNum=ctx.eventID().event_number();
268 double RunStr=double(RunNum);
269 double EvtStr=double(EvtNum);
270 m_c->EventNo=EvtStr;
271 m_c->RunNo=RunStr;
272 m_histos[160]->Fill(EvtStr);
273 m_histos[159]->Fill(RunNum);
274
275 MagField::AtlasFieldCache fieldCache;
276 // Get field cache object
277 SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{m_fieldCacheCondObjInputKey, ctx};
278 const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
279 if (fieldCondObj == nullptr) {
280 ATH_MSG_ERROR("Failed to retrieve AtlasFieldCacheCondObj with key " << m_fieldCacheCondObjInputKey.key());
281 return StatusCode::FAILURE;
282 }
283 fieldCondObj->getInitializedCache (fieldCache);
284
285 // Get the MC Truth Information
286 SG::ReadHandle<McEventCollection> mcEvent{m_truthKey, ctx};
287 for (const HepMC::GenEvent* e : *mcEvent) {
288
289 // Get just the primary, call it "theParticle"
290 auto theParticle = *HepMC::begin(*e);
291
292 // Get the kinematic variables:
293 HepLorentzVector momentum(theParticle->momentum().px(),
294 theParticle->momentum().py(),
295 theParticle->momentum().pz(),
296 theParticle->momentum().e());
297 Point3D<double> origin(theParticle->production_vertex()->position().x(),
298 theParticle->production_vertex()->position().y(),
299 theParticle->production_vertex()->position().z());
300 double charge = MC::charge(theParticle->pdg_id());
301 // Put Eta and Phi into the Ntuple
302 m_c->phi = theParticle->momentum().phi();
303 m_c->eta = -log(tan(theParticle->momentum().theta()/2));
304 if (!finite(m_c->eta)) m_c->eta=0;
305 m_c->pt = theParticle->momentum().perp();
306 int partID = theParticle->pdg_id();
307 double pID = double(partID);
308 m_c->PDG = pID;
309 m_c->Energy = theParticle->momentum().e();
310 m_histos[2]->Fill( theParticle->momentum().phi() );
311 double myEta = -log(tan(theParticle->momentum().theta()/2));
312 if (!finite(myEta)) m_histos[0]->Fill(0);
313 else m_histos[0]->Fill( myEta );
314 m_histos[158]->Fill( pID );
315 m_histos[1]->Fill( theParticle->momentum().perp()/Units::GeV );
316 m_histos[157]->Fill( m_c->Energy = theParticle->momentum().e()/Units::GeV );
317
318 // Make an extrapolator:
319 const Genfun::AtlasBComponent Bx(0,&fieldCache);
320 const Genfun::AtlasBComponent By(1,&fieldCache);
321 const Genfun::AtlasBComponent Bz(2,&fieldCache);
322 GeoXPEngine extrapolator(Bx, By, Bz, origin, momentum, charge);
323
324 // Extrapolate to the first layer in the barrel:
325 m_c->x = 0;
326 m_c->y = 0;
327 m_c->z = 0;
328 double x=0,y=0,z=0;
329 bool hitsBarrel=false;
330 //bool hitsEndcap=false;
331 for (double t = 0; t< 50; t += 0.1) {
332 x = extrapolator.x()(t);
333 y = extrapolator.y()(t);
334 z = extrapolator.z()(t);
335 double magicZ=3640.0*mm;
336 double magicR=1500.0*mm;
337 if (x*x+y*y > magicR*magicR) {
338 m_c->x = x;
339 m_c->y = y;
340 m_c->z = z;
341 hitsBarrel=true;
342 break;
343 }
344 else if (z*z > magicZ*magicZ) {
345 m_c->x = x;
346 m_c->y = y;
347 m_c->z = z;
348 //hitsEndcap=true;
349 break;
350 }
351 }
352
353 m_histos[3]->Fill( x );
354 m_histos[4]->Fill( y );
355 m_histos[5]->Fill( z );
356
357 // You have an x,y, and z position. Now go and get the Element corresponding to
358 // that hit position. There are four, one for each sampling layer:
359 double radImpact = std::sqrt(x*x+y*y+z*z);
360 double phiImpact = std::atan2(y,x);
361 double thetaImpact = std::acos(z/radImpact);
362 double etaImpact = -std::log(std::tan(thetaImpact/2));
363
364 SG::ReadCondHandle<CaloDetDescrManager> caloMgrHandle{m_caloMgrKey, ctx};
365 ATH_CHECK(caloMgrHandle.isValid());
366 const CaloDetDescrManager* caloMgr = *caloMgrHandle;
367 const CaloDetDescrElement *element[15]={nullptr};
368
369 for (int i=0;i<4;i++) {
370 try {
371 element[i] = caloMgr->get_element(CaloCell_ID::LAREM,i, hitsBarrel, etaImpact, phiImpact);
372 }
373 catch (const LArID_Exception & e) {
374 std::cerr << "SingleTrackValidation EXCEPTION (LAREM)" << e.message() << std::endl;
375 }
376 }
377 for (int i=0;i<4;i++) {
378 try {
379 element[i+4] = caloMgr->get_element(CaloCell_ID::LAREM,i, hitsBarrel, etaImpact, phiImpact);
380 }
381 catch (const LArID_Exception & e) {
382 std::cerr << "SingleTrackValidation EXCEPTION (LAREM)" << e.message() << std::endl;
383 }
384 }
385 for (int i=0;i<4;i++) {
386 try {
387 element[i+8] = caloMgr->get_element(CaloCell_ID::LARHEC,i, hitsBarrel, etaImpact, phiImpact);
388 }
389 catch (const LArID_Exception & e) {
390 std::cerr << "SingleTrackValidation EXCEPTION in (LARHEC)" << e.message() << std::endl;
391 }
392 }
393 for (int i=1;i<4;i++) {
394 try {
395 element[i+11] = caloMgr->get_element(CaloCell_ID::LARFCAL,i, hitsBarrel, etaImpact, phiImpact);
396 }
397 catch (const LArID_Exception & e) {
398 std::cerr << "SingleTrackValidation EXCEPTIONin LARFCAL" << e.message() << std::endl;
399 }
400 }
401
402
403 // Now go and find out how much energy is there:
404 for (int z=0;z<15;z++){
405 m_c->s_c00[z]=0;
406 m_c->s_t00[z]=0;
407 }
408
409 std::string lArKey = hitsBarrel ? "LArHitEMB" : "LArHitEMEC" ;
410
411 double eSum [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
412 double eEta [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
413 double eEtaEta [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
414 double ePhi [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
415 double ePhiPhi [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
416 int hit_count [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
417 double eX [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
418 double eXX [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
419 double eY [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
420 double eYY [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
421 double c00 [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
422 double t00 [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
423 // double width [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
424 double e_dep = 0;
425
426 for (int i=0;i<4;i++) { // Loop over the four LAr collections
427 if (i==0) {
428 lArKey="LArHitEMB";
429 } else if (i==1) {
430 lArKey="LArHitEMEC";
431 } else if (i==2) {
432 lArKey="LArHitHEC";
433 } else if (i==3) {
434 lArKey="LArHitFCAL";
435 }
436
437 SG::ReadHandle<LArHitContainer> larHitContainer(lArKey, "StoreGateSvc");
438 if (larHitContainer.isValid()) {
439 for (const LArHit* larHit : *larHitContainer) {
440 const CaloDetDescrElement *hitElement = caloMgr->get_element(larHit->cellID());
441 int samplingLayer = m_c->cellId->sampling(larHit->cellID());
442 double energy = larHit->energy();
443
444 for (int j=0;j<15;j++) {
445 if (hitElement==element[j]) {
446 c00[j] += energy;
447 }
448 }
449
450 if (lArKey=="LArHitEMEC") samplingLayer += 4;
451 else if (lArKey=="LArHitHEC") samplingLayer += 8;
452 else if (lArKey=="LArHitFCAL") samplingLayer += 11;
453
454 // Calculate phi of hit w.r.t phiImpact
455 // to avoid problems due to the 2pi modulus of phi
456 // we are thus calculating deltaphi directly, so don't subtract phiImpact again later
457 double hitPhi = hitElement->phi() - phiImpact;
458 if (hitPhi < -pi) hitPhi += twopi;
459 if (hitPhi > pi) hitPhi -= twopi;
460
461 eSum [samplingLayer]+=energy;
462 eEta [samplingLayer]+=energy*hitElement->eta();
463 eEtaEta [samplingLayer]+=energy*hitElement->eta()*hitElement->eta();
464 ePhi [samplingLayer]+=energy*hitPhi;
465 ePhiPhi [samplingLayer]+=energy*hitPhi*hitPhi;
466 eX [samplingLayer]+=energy*hitElement->x();
467 eXX [samplingLayer]+=energy*hitElement->x()*hitElement->x();
468 eY [samplingLayer]+=energy*hitElement->y();
469 eYY [samplingLayer]+=energy*hitElement->y()*hitElement->y();
470 t00 [samplingLayer]+=energy*larHit->time();
471 hit_count[samplingLayer]+=1;
472 }
473 }
474 }
475
476 for (int i=0;i<15;i++) {
477 if (eSum[i]!=0) eEta[i] /= eSum[i];
478 if (eSum[i]!=0) eEtaEta[i] /= eSum[i];
479 if (eSum[i]!=0) ePhi[i] /= eSum[i];
480 if (eSum[i]!=0) ePhiPhi[i] /= eSum[i];
481 if (eSum[i]!=0) eY[i] /= eSum[i];
482 if (eSum[i]!=0) eYY[i] /= eSum[i];
483 if (eSum[i]!=0) eX[i] /= eSum[i];
484 if (eSum[i]!=0) eXX[i] /= eSum[i];
485 if (eSum[i]!=0) t00[i] /= eSum[i];
486 e_dep+=eSum[i];
487 }
488
489 double dThetaDEta = -1.0/cosh(etaImpact);
490
491 //Fill the E Sum, center cell E, hit count fields:
492 m_c->E_Deposit=e_dep;
493 for (int z=0;z<15;z++){
494 m_c->s_sumE[z]=eSum[z];
495 m_c->s_c00[z]=c00[z];
496 m_c->s_t00[z]=t00[z];
497 m_c->s_hits[z]=hit_count[z];
498 if (z<12){
499 m_c->s_deltaPhi[z]=radImpact*std::sin(thetaImpact)*(ePhi[z]);
500 m_c->s_sigmaPhi[z]=radImpact*std::sin(thetaImpact)*std::sqrt(ePhiPhi[z]- ePhi[z]*ePhi[z]);
501 m_c->s_deltaEta[z]=radImpact*dThetaDEta*(eEta[z]-etaImpact);
502 m_c->s_sigmaEta[z]=radImpact*std::fabs(dThetaDEta)*std::sqrt(eEtaEta[z]- eEta[z]*eEta[z]);
503 } else {
504 m_c->s_deltaPhi[z]=(eX[z]-x);
505 m_c->s_sigmaPhi[z]=std::sqrt(eXX[z]- eX[z]*eX[z]);
506 m_c->s_deltaEta[z]=(eY[z]-y);
507 m_c->s_sigmaEta[z]=std::sqrt(eYY[z]-eY[z]*eY[z]);
508 }
509 m_c->s_widthX[z]=std::sqrt(eXX[z]-eX[z]*eX[z]);
510 m_c->s_widthY[z]=std::sqrt(eYY[z]-eY[z]*eY[z]);
511 }
512
513 m_histos[161]->Fill(e_dep/Units::GeV);
514 for (int i=0;i<15;i++){
515 m_histos[6+i]->Fill( c00[i]/Units::GeV );
516 m_histos[21+i]->Fill( hit_count[i] );
517 m_histos[36+i]->Fill( eSum[i]/Units::GeV );
518 m_histos[111+i]->Fill( t00[i] );
519 m_histos[126+i]->Fill( sqrt(eXX[i]-eX[i]*eX[i]) );
520 m_histos[141+i]->Fill( sqrt(eYY[i]-eY[i]*eY[i]) );
521 if (i<8){
522 m_histos[51+i]->Fill( radImpact*std::sin(thetaImpact)*ePhi[i] );
523 m_histos[66+i]->Fill( radImpact*std::sin(thetaImpact)*std::sqrt(ePhiPhi[i]-ePhi[i]*ePhi[i]) );
524 m_histos[81+i]->Fill( radImpact*dThetaDEta*(eEta[i]-etaImpact) );
525 m_histos[96+i]->Fill( radImpact*std::fabs(dThetaDEta)*std::sqrt(eEtaEta[i]-eEta[i]*eEta[i]) );
526 } else {
527 m_histos[51+i]->Fill( eX[i]-x );
528 m_histos[66+i]->Fill( std::sqrt(eXX[i]-eX[i]*eX[i]) );
529 m_histos[81+i]->Fill( eY[i]-y );
530 m_histos[96+i]->Fill( std::sqrt(eYY[i]-eY[i]*eY[i]) );
531 }
532 }
533
534 ATH_CHECK(ntupleSvc()->writeRecord(m_c->nt));
535
536 }
537
538 m_c->cpuTime= getCpu();
539
540 return StatusCode::SUCCESS;
541
542}
543
544StatusCode SingleTrackValidation::finalize() {
545 return StatusCode::SUCCESS;
546}
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_INFO
#define ATH_MSG_INFO(x)
Definition AthMsgStreamMacros.h:31
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition AthMsgStreamMacros.h:32
charge
double charge(const T &p)
Definition AtlasPID.h:997
CaloCell_ID.h
HepMCHelpers.h
ATLAS-specific HepMC functions.
LArID_Exception.h
a
static Double_t a
Definition LArPhysWaveHECTool.cxx:38
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:14
iabs
int iabs(int a)
Definition SingleTrackValidation.cxx:52
getCpu
int getCpu()
Definition SingleTrackValidation.cxx:53
pi
#define pi
Definition TileMuonFitter.cxx:65
Units.h
Wrapper to avoid constant divisions when using units.
y
#define y
x
#define x
z
#define z
twopi
constexpr double twopi
Definition VertexPointEstimator.cxx:16
AthAlgorithm::AthAlgorithm
AthAlgorithm(const std::string &name, ISvcLocator *pSvcLocator)
Constructor.
Definition AthAlgorithm.cxx:10
AthCommonAlgorithm< Gaudi::Algorithm >::detStore
const ServiceHandle< StoreGateSvc > & detStore() const
Definition AthCommonDataStore.h:95
AtlasFieldCacheCondObj::getInitializedCache
void getInitializedCache(MagField::AtlasFieldCache &cache) const
get B field cache for evaluation as a function of 2-d or 3-d position.
Definition AtlasFieldCacheCondObj.h:32
CaloCell_Base_ID::LAREM
@ LAREM
Definition CaloCell_Base_ID.h:45
CaloCell_Base_ID::LARFCAL
@ LARFCAL
Definition CaloCell_Base_ID.h:45
CaloCell_Base_ID::LARHEC
@ LARHEC
Definition CaloCell_Base_ID.h:45
CaloCell_ID
Helper class for offline cell identifiers.
Definition CaloCell_ID.h:34
CaloDetDescrElement
This class groups all DetDescr information related to a CaloCell.
Definition Calorimeter/CaloDetDescr/CaloDetDescr/CaloDetDescrElement.h:66
CaloDetDescrManager_Base::get_element
const CaloDetDescrElement * get_element(const Identifier &cellId) const
get element by its identifier
Definition CaloDetDescrManager.cxx:159
CaloDetDescrManager
This class provides the client interface for accessing the detector description information common to...
Definition CaloDetDescrManager.h:469
GeoXPEngine::y
const Genfun::AbsFunction & y() const
Definition GeoXPEngine.cxx:80
GeoXPEngine::z
const Genfun::AbsFunction & z() const
Definition GeoXPEngine.cxx:83
GeoXPEngine::x
const Genfun::AbsFunction & x() const
Definition GeoXPEngine.cxx:77
LArHit
Class to store hit energy and time in LAr cell from G4 simulation.
Definition LArHit.h:25
LArID_Exception
Exception class for LAr Identifiers.
Definition LArID_Exception.h:20
MagField::AtlasFieldCache
Local cache for magnetic field (based on MagFieldServices/AtlasFieldSvcTLS.h).
Definition AtlasFieldCache.h:43
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
SingleTrackValidation::Clockwork::s_deltaEta
NTuple::Item< double > s_deltaEta[15]
Definition SingleTrackValidation.cxx:92
SingleTrackValidation::Clockwork::s_sigmaPhi
NTuple::Item< double > s_sigmaPhi[15]
Definition SingleTrackValidation.cxx:91
SingleTrackValidation::Clockwork::s_sigmaEta
NTuple::Item< double > s_sigmaEta[15]
Definition SingleTrackValidation.cxx:93
SingleTrackValidation::Clockwork::s_deltaPhi
NTuple::Item< double > s_deltaPhi[15]
Definition SingleTrackValidation.cxx:90
SingleTrackValidation::m_caloMgrKey
SG::ReadCondHandleKey< CaloDetDescrManager > m_caloMgrKey
Definition SingleTrackValidation.h:34
SingleTrackValidation::m_ppSvc
ServiceHandle< IPartPropSvc > m_ppSvc
Definition SingleTrackValidation.h:45
SingleTrackValidation::m_histSvc
ServiceHandle< ITHistSvc > m_histSvc
Definition SingleTrackValidation.h:46
SingleTrackValidation::m_truthKey
SG::ReadHandleKey< McEventCollection > m_truthKey
Definition SingleTrackValidation.h:32
SingleTrackValidation::finalize
StatusCode finalize() override
Definition SingleTrackValidation.cxx:544
SingleTrackValidation::initialize
StatusCode initialize() override
Definition SingleTrackValidation.cxx:118
SingleTrackValidation::execute
StatusCode execute(const EventContext &ctx) override
Execute method.
Definition SingleTrackValidation.cxx:257
SingleTrackValidation::SingleTrackValidation
SingleTrackValidation(const std::string &name, ISvcLocator *pSvcLocator)
Definition SingleTrackValidation.cxx:108
SingleTrackValidation::m_fieldCacheCondObjInputKey
SG::ReadCondHandleKey< AtlasFieldCacheCondObj > m_fieldCacheCondObjInputKey
Definition SingleTrackValidation.h:40
HepMC::begin
std::vector< HepMC3::GenParticlePtr >::const_iterator begin(HepMC3::GenEvent &e)
Definition GenEvent.h:355
HepMC::GenEvent
HepMC3::GenEvent GenEvent
Definition GenEvent.h:39
MC::charge
double charge(const T &p)
Definition HepMCHelpers.h:998
file
TFile * file
Definition tile_monitor.h:29
Generated on for ATLAS Offline Software by doxygen 1.16.1