10 #include "GaudiKernel/NTuple.h"
11 #include "GaudiKernel/INTupleSvc.h"
12 #include "GaudiKernel/SmartDataPtr.h"
15 #include "CLHEP/GenericFunctions/FixedConstant.hh"
18 #include "HepPDT/ParticleDataTable.hh"
19 #include "HepPDT/ParticleData.hh"
26 #include "CaloDetDescr/CaloDetDescrElement.h"
32 #include "CLHEP/Units/PhysicalConstants.h"
38 #include <sys/types.h>
42 #include <sys/times.h>
46 using HepGeom::Point3D;
47 using CLHEP::HepLorentzVector;
52 inline int iabs (
int a) {
return (
a>=0)?
a : -
a;}
56 static const int ticksPerJiffy = sysconf(_SC_CLK_TCK)/100;
59 cpuTime=(
buff.tms_utime +
buff.tms_stime +
buff.tms_cutime +
buff.tms_cstime)/ticksPerJiffy;
74 NTuple::Tuple*
nt =
nullptr;
76 NTuple::Item<double>
eta;
77 NTuple::Item<double>
pt;
78 NTuple::Item<double>
phi;
81 NTuple::Item<double>
x;
82 NTuple::Item<double>
y;
83 NTuple::Item<double>
z;
115 if (
m_c!=
nullptr){
delete m_c;
m_c=
nullptr; }
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" };
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},
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},
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},
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},
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},
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},
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},
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},
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},
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},
153 {0,50}, {0,100}, {-25,25}, {0,10}, {0,1000}, {0,10}};
174 for (
int i=0;
i<162;
i++){
179 ATH_MSG_WARNING(
"Failed to register historam " <<
names[
i] <<
". Not sure what will happen now..." );
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");
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..." );
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);
209 if (
i<12) sprintf(
title,
"S%i_SumE",
i);
210 else sprintf(
title,
"FC%i_SumE",
i-11);
212 if (
i<12) sprintf(
title,
"S%i_Hits",
i);
213 else sprintf(
title,
"FC%i_Hits",
i-11);
215 if (
i<12) sprintf(
title,
"S%i_DeltaPhi",
i);
216 else sprintf(
title,
"FC%i_DeltaX",
i-11);
218 if (
i<12) sprintf(
title,
"S%i_SigmaPhi",
i);
219 else sprintf(
title,
"FC%i_SigmaX",
i-11);
221 if (
i<12) sprintf(
title,
"S%i_DeltaEta",
i);
222 else sprintf(
title,
"FC%i_DeltaY",
i-11);
224 if (
i<12) sprintf(
title,
"S%i_SigmaEta",
i);
225 else sprintf(
title,
"FC%i_SigmaY",
i-11);
227 if (
i<12) sprintf(
title,
"S%i_Time",
i);
228 else sprintf(
title,
"FC%i_Time",
i-11);
230 if (
i<12) sprintf(
title,
"S%i_WidthX",
i);
231 else sprintf(
title,
"FC%i_WidthX",
i-11);
233 if (
i<12) sprintf(
title,
"S%i_WidthY",
i);
234 else sprintf(
title,
"FC%i_WidthY",
i-11);
239 nt->addItem(
"TrackEnergy" ,
m_c->
Energy ).isFailure() ||
240 nt->addItem(
"ParticleID" ,
m_c->
PDG ).isFailure() ||
241 nt->addItem(
"Run#" ,
m_c->
RunNo ).isFailure() ||
244 ATH_MSG_WARNING(
"Registration of some of the ntuple branches failed. No idea what will happen next..." );
254 return StatusCode::SUCCESS;
266 const EventContext& context = getContext();
267 int RunNum=context.eventID().run_number();
268 int EvtNum=context.eventID().event_number();
269 double RunStr=
double(RunNum);
270 double EvtStr=
double(EvtNum);
280 if (fieldCondObj ==
nullptr) {
282 return StatusCode::FAILURE;
284 fieldCondObj->getInitializedCache (fieldCache);
288 for (
const HepMC::GenEvent*
e : *mcEvent) {
294 const HepPDT::ParticleDataTable * dataTable =
m_c->
partPropSvc->PDT();
295 const HepPDT::ParticleData * particleData = dataTable->particle(
iabs(theParticle->pdg_id()));
298 HepLorentzVector
momentum(theParticle->momentum().px(),
299 theParticle->momentum().py(),
300 theParticle->momentum().pz(),
301 theParticle->momentum().e());
302 Point3D<double> origin(theParticle->production_vertex()->position().x(),
303 theParticle->production_vertex()->position().y(),
304 theParticle->production_vertex()->position().z());
305 double charge = theParticle->pdg_id() > 0 ? particleData->charge() : - particleData->charge();
307 m_c->
phi = theParticle->momentum().phi();
308 m_c->
eta = -
log(
tan(theParticle->momentum().theta()/2));
310 m_c->
pt = theParticle->momentum().perp();
311 int partID = theParticle->pdg_id();
314 m_c->
Energy = theParticle->momentum().e();
315 m_histos[2]->Fill( theParticle->momentum().phi() );
316 double myEta = -
log(
tan(theParticle->momentum().theta()/2));
317 if (!finite(myEta))
m_histos[0]->Fill(0);
334 bool hitsBarrel=
false;
336 for (
double t = 0;
t< 50;
t += 0.1) {
337 x = extrapolator.
x()(
t);
338 y = extrapolator.
y()(
t);
339 z = extrapolator.
z()(
t);
340 double magicZ=3640.0*
mm;
341 double magicR=1500.0*
mm;
342 if (
x*
x+
y*
y > magicR*magicR) {
349 else if (
z*
z > magicZ*magicZ) {
364 double radImpact = std::sqrt(
x*
x+
y*
y+
z*
z);
365 double phiImpact = std::atan2(
y,
x);
366 double thetaImpact = std::acos(
z/radImpact);
374 for (
int i=0;
i<4;
i++) {
379 std::cerr <<
"SingleTrackValidation EXCEPTION (LAREM)" <<
e.message() << std::endl;
382 for (
int i=0;
i<4;
i++) {
387 std::cerr <<
"SingleTrackValidation EXCEPTION (LAREM)" <<
e.message() << std::endl;
390 for (
int i=0;
i<4;
i++) {
395 std::cerr <<
"SingleTrackValidation EXCEPTION in (LARHEC)" <<
e.message() << std::endl;
398 for (
int i=1;
i<4;
i++) {
403 std::cerr <<
"SingleTrackValidation EXCEPTIONin LARFCAL" <<
e.message() << std::endl;
409 for (
int z=0;
z<15;
z++){
414 std::string lArKey = hitsBarrel ?
"LArHitEMB" :
"LArHitEMEC" ;
416 double eSum [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
417 double eEta [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
418 double eEtaEta [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
419 double ePhi [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
420 double ePhiPhi [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
421 int hit_count [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
422 double eX [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
423 double eXX [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
424 double eY [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
425 double eYY [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
426 double c00 [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
427 double t00 [15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
431 for (
int i=0;
i<4;
i++) {
443 if (larHitContainer.
isValid()) {
444 for (
const LArHit* larHit : *larHitContainer) {
447 double energy = larHit->energy();
449 for (
int j=0;j<15;j++) {
450 if (hitElement==element[j]) {
455 if (lArKey==
"LArHitEMEC") samplingLayer += 4;
456 else if (lArKey==
"LArHitHEC") samplingLayer += 8;
457 else if (lArKey==
"LArHitFCAL") samplingLayer += 11;
462 double hitPhi = hitElement->
phi() - phiImpact;
463 if (hitPhi < -
pi) hitPhi +=
twopi;
464 if (hitPhi >
pi) hitPhi -=
twopi;
466 eSum [samplingLayer]+=
energy;
467 eEta [samplingLayer]+=
energy*hitElement->
eta();
468 eEtaEta [samplingLayer]+=
energy*hitElement->
eta()*hitElement->
eta();
469 ePhi [samplingLayer]+=
energy*hitPhi;
470 ePhiPhi [samplingLayer]+=
energy*hitPhi*hitPhi;
471 eX [samplingLayer]+=
energy*hitElement->
x();
472 eXX [samplingLayer]+=
energy*hitElement->
x()*hitElement->
x();
473 eY [samplingLayer]+=
energy*hitElement->
y();
474 eYY [samplingLayer]+=
energy*hitElement->
y()*hitElement->
y();
475 t00 [samplingLayer]+=
energy*larHit->time();
476 hit_count[samplingLayer]+=1;
481 for (
int i=0;
i<15;
i++) {
482 if (eSum[
i]!=0) eEta[
i] /= eSum[
i];
483 if (eSum[
i]!=0) eEtaEta[
i] /= eSum[
i];
484 if (eSum[
i]!=0) ePhi[
i] /= eSum[
i];
485 if (eSum[
i]!=0) ePhiPhi[
i] /= eSum[
i];
486 if (eSum[
i]!=0) eY[
i] /= eSum[
i];
487 if (eSum[
i]!=0) eYY[
i] /= eSum[
i];
488 if (eSum[
i]!=0) eX[
i] /= eSum[
i];
489 if (eSum[
i]!=0) eXX[
i] /= eSum[
i];
490 if (eSum[
i]!=0) t00[
i] /= eSum[
i];
494 double dThetaDEta = -1.0/cosh(etaImpact);
498 for (
int z=0;
z<15;
z++){
507 m_c->
s_sigmaEta[
z]=radImpact*std::fabs(dThetaDEta)*std::sqrt(eEtaEta[
z]- eEta[
z]*eEta[
z]);
519 for (
int i=0;
i<15;
i++){
529 m_histos[81+
i]->Fill( radImpact*dThetaDEta*(eEta[
i]-etaImpact) );
530 m_histos[96+
i]->Fill( radImpact*std::fabs(dThetaDEta)*std::sqrt(eEtaEta[
i]-eEta[
i]*eEta[
i]) );
545 return StatusCode::SUCCESS;
550 return StatusCode::SUCCESS;