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