5 #include "GaudiKernel/IPartPropSvc.h"
6 #include "GaudiKernel/ServiceHandle.h"
13 #include "HepPDT/ParticleDataTable.hh"
25 (
const std::string&
name,ISvcLocator* pSvcLocator)
58 IPartPropSvc* partPropSvc =
nullptr;
59 sc = service(
"PartPropSvc", partPropSvc,
true);
62 return StatusCode::FAILURE;
67 m_particleDataTable = partPropSvc->PDT();
68 if(!m_particleDataTable) {
70 return StatusCode::FAILURE;
75 m_pdg = std::abs(m_pdg) ;
77 m_trackCollectionStat.resize(m_tracklocation.size());
81 ATH_CHECK(m_truth_locationStrip.initialize(m_useStrip));
82 ATH_CHECK(m_clustersStripname.initialize(m_useStrip));
83 ATH_CHECK(m_spacepointsStripname.initialize(m_useStrip));
85 ATH_CHECK( m_clustersPixelname.initialize(m_usePix));
86 ATH_CHECK( m_spacepointsPixelname.initialize(m_usePix));
87 ATH_CHECK( m_truth_locationPixel.initialize(m_usePix));
89 ATH_CHECK( m_spacepointsOverlapname.initialize(m_useStrip));
94 ATH_CHECK(m_pixelDetEleCollKey.initialize());
95 ATH_CHECK(m_StripDetEleCollKey.initialize());
112 if(!m_usePix && !m_useStrip)
return StatusCode::SUCCESS;
115 std::vector<SG::ReadHandle<PRD_MultiTruthCollection> > read_handle;
116 read_handle.reserve(3);
118 read_handle.emplace_back(m_truth_locationPixel,ctx);
119 if (not read_handle.back().isValid()) {
121 return StatusCode::FAILURE;
123 event_data.
m_truthPix = &(*read_handle.back());
127 read_handle.emplace_back(m_truth_locationStrip,ctx);
128 if (not read_handle.back().isValid()) {
130 return StatusCode::FAILURE;
135 newClustersEvent (ctx,event_data);
136 newSpacePointsEvent (ctx,event_data);
137 event_data.
m_nqtracks = qualityTracksSelection(event_data);
138 tracksComparison (ctx,event_data);
141 efficiencyReconstruction(event_data);
142 if(msgLvl(
MSG::DEBUG)) noReconstructedParticles(event_data);
147 std::lock_guard<std::mutex> lock(m_statMutex);
149 for (
unsigned int i=0;
i< m_trackCollectionStat.size(); ++
i ) {
156 dumpevent(
msg(),event_data);
159 return StatusCode::SUCCESS;
167 if(m_eventStat.m_events<=0)
return StatusCode::SUCCESS;
168 const auto & w13 = std::setw(13);
169 const auto & p5 = std::setprecision(5);
170 const auto topNtail=[](
const std::string &
str){
return "|" +
str +
"|";};
171 const std::string lineSeparator(83,
'-');
172 const std::string spaceSeparator(83,
' ');
173 std::stringstream
out;
175 out<<topNtail(lineSeparator)<<
"\n";
176 out<<
"| TrackClusterAssValidation statistic for charge truth particles with |\n";
177 out<<topNtail(spaceSeparator)<<
"\n";
179 out<<
"| eta bins for eta dependent variables = [0.0, ";
180 for (
unsigned int etabin = 1; etabin<(m_etabins.size()-1); etabin++)
181 out << std::setw(2) << std::setprecision(1) << m_etabins.value().at(etabin) <<
", ";
182 out << std::setw(2) << m_etabins.value().back() <<
"] |\n";
183 out<<
"| eta dependent pT [MeV] >= [";
184 for (
unsigned int ptbin = 0; ptbin<(m_ptcuts.size()-1); ptbin++)
185 out<<std::setw(6)<<std::setprecision(2)<<m_ptcuts.value().at(ptbin)<<
", ";
186 out<<std::setw(6)<<std::setprecision(2)<<m_ptcuts.value().back()<<
"] |\n";
187 if(m_ptcutmax < 1000000.) {
188 out<<
"| pT <="<<w13<<p5<<m_ptcutmax<<
" MeV"<<
" |\n";
190 out<<
"| |rapidity| <="<<w13<<p5<<m_rapcut<<
" |\n";
191 out<<
"| max vertex radius <="<<w13<<p5<<m_rmax<<
" mm"<<
" |\n";
192 out<<
"| min vertex radius >="<<w13<<p5<<m_rmin<<
" mm"<<
" |\n";
193 out<<
"| particles pdg ="<<std::setw(8)<<m_pdg<<
" |\n";
195 auto yesNo=[](
const bool predicate){
return predicate ?
"yes" :
"no ";};
196 out<<
"| use Pixels information "<<yesNo(m_usePix)<<
" |\n";
197 out<<
"| use Strip information "<<yesNo(m_useStrip)<<
" |\n";
198 out<<
"| take into account outliers "<<yesNo(m_useOutliers)<<
" |\n";
199 out<<topNtail(spaceSeparator)<<
"\n";
200 if(!m_usePix && !m_useStrip)
return StatusCode::SUCCESS;
201 enum Regions{Barrel, Transition, Endcap, Forward, NRegions};
202 auto incrementArray=[](
auto &
array,
const int idx){
for (
int j{};j!=NRegions;++j)
array[
idx][j] +=
array[
idx+1][j];};
203 for(
int i=48;
i>=0; --
i) {
204 m_eventStat.m_particleClusters [
i] +=m_eventStat.m_particleClusters [
i+1];
205 incrementArray(m_eventStat.m_particleClustersBTE,
i);
206 m_eventStat.m_particleSpacePoints [
i] +=m_eventStat.m_particleSpacePoints [
i+1];
207 incrementArray(m_eventStat.m_particleSpacePointsBTE,
i);
209 auto coerceToOne=[](
const double & initialVal)->
double {
return (initialVal<1.) ? 1. : initialVal; };
216 double pa = coerceToOne(m_eventStat.m_particleClusters[0]);
217 std::array<double, 10> pc2ff{};
218 size_t clusterIdx = 2;
219 for (
auto & thisCluster: pc2ff){
220 thisCluster =
double(m_eventStat.m_particleClusters[ clusterIdx++ ])/ pa;
223 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Barrel]);
224 std::array<double, 10> pcBarrel2ff{};
225 size_t clusterBarrelIdx = 2;
226 for (
auto & thisClusterB: pcBarrel2ff){
227 thisClusterB =
double(m_eventStat.m_particleClustersBTE[ clusterBarrelIdx++ ][Barrel])/ pa;
230 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Transition]);
231 std::array<double, 10> pcTransition2ff{};
232 size_t clusterTransitionIdx = 2;
233 for (
auto & thisClusterT: pcTransition2ff){
234 thisClusterT =
double(m_eventStat.m_particleClustersBTE[ clusterTransitionIdx++ ][Transition])/ pa;
237 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Endcap]);
238 std::array<double, 10> pcEndcap2ff{};
239 size_t clusterEndcapIdx = 2;
240 for (
auto & thisClusterE: pcEndcap2ff){
241 thisClusterE =
double(m_eventStat.m_particleClustersBTE[ clusterEndcapIdx++ ][Endcap])/ pa;
244 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][3]);
245 std::array<double, 10> pcFwd2ff{};
246 size_t clusterFwdIdx = 2;
247 for (
auto & thisClusterD: pcFwd2ff){
248 thisClusterD =
double(m_eventStat.m_particleClustersBTE[ clusterFwdIdx++ ][Forward])/ pa;
254 pa = coerceToOne(m_eventStat.m_particleSpacePoints[0]);
255 std::array<double, 10> sp2ff{};
256 size_t spacepointIdx = 2;
257 for (
auto & thisSpacepoint: sp2ff){
258 thisSpacepoint =
double(m_eventStat.m_particleSpacePoints[ spacepointIdx++ ])/ pa;
261 pa = coerceToOne(m_eventStat.m_particleSpacePointsBTE[0][Barrel]);
262 std::array<double, 10> spBarrel2ff{};
263 size_t spacepointBarrelIdx = 2;
264 for (
auto & thisSpacepoint: spBarrel2ff){
265 thisSpacepoint =
double(m_eventStat.m_particleSpacePointsBTE[ spacepointBarrelIdx++ ][Barrel])/ pa;
268 pa = coerceToOne(m_eventStat.m_particleSpacePointsBTE[0][Transition]);
269 std::array<double, 10> spTransition2ff{};
270 size_t spacepointTransitionIdx = 2;
271 for (
auto & thisSpacepoint: spTransition2ff){
272 thisSpacepoint =
double(m_eventStat.m_particleSpacePointsBTE[ spacepointTransitionIdx++ ][Transition])/ pa;
275 pa = coerceToOne(m_eventStat.m_particleSpacePointsBTE[0][Endcap]);
276 std::array<double, 10> spEndcap2ff{};
277 size_t spacepointEndcapIdx = 2;
278 for (
auto & thisSpacepoint: spEndcap2ff){
279 thisSpacepoint =
double(m_eventStat.m_particleSpacePointsBTE[ spacepointEndcapIdx++ ][Endcap])/ pa;
282 pa = coerceToOne(m_eventStat.m_particleSpacePointsBTE[0][Forward]);
283 std::array<double, 10> spFwd2ff{};
284 size_t spacepointFwdIdx = 2;
285 for (
auto & thisSpacepoint: spFwd2ff){
286 thisSpacepoint =
double(m_eventStat.m_particleSpacePointsBTE[ spacepointFwdIdx++ ][Forward])/ pa;
288 auto w8=std::setw(8);
289 out<<
"| Probability for such charge particles to have some number silicon |\n";
290 out<<
"| clusters | space points |\n";
291 out<<
"| Total Barrel Transi Endcap Forward | Total Barrel Transi Endcap Forward |\n";
294 out<<
"| >= "<<
idx+2<< std::string((
idx<8)?
" ":
" ")
296 <<w8<<p5<<pcBarrel2ff[
idx]
297 <<w8<<p5<<pcTransition2ff[
idx]
298 <<w8<<p5<<pcEndcap2ff[
idx]
299 <<w8<<p5<<pcFwd2ff[
idx]<<
" | "
302 <<w8<<p5<<spBarrel2ff[
idx]
303 <<w8<<p5<<spTransition2ff[
idx]
304 <<w8<<p5<<spEndcap2ff[
idx]
305 <<w8<<p5<<spFwd2ff[
idx]
309 out<<topNtail(spaceSeparator)<<
"\n";
310 out<<
"| Additional cuts for truth particles are |\n";
311 out<<
"| eta dependent number of silicon clusters >= [";
312 for (
unsigned int clbin = 0; clbin<(m_clcuts.size()-1); clbin++)
313 out<<std::setw(2)<<m_clcuts.value().at(clbin)<<
", ";
314 out<<std::setw(2)<<m_clcuts.value().back()<<
"] |\n";
315 out<<
"| number space points >="<<w13<<m_spcut.value()<<
" |\n";
317 pa = coerceToOne(m_eventStat.m_particleClusters[0]);
318 out<<
"| Probability find truth particles with this cuts is "<<w8<<p5<<
double(m_eventStat.m_events)/pa<<
" |\n";
319 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Barrel]);
320 out<<
"| For barrel region "<<w8<<p5<<
double(m_eventStat.m_eventsBTE[Barrel])/pa<<
" |\n";
321 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Transition]);
322 out<<
"| For transition region "<<w8<<p5<<
double(m_eventStat.m_eventsBTE[Transition])/pa<<
" |\n";
323 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Endcap]);
324 out<<
"| For endcap region "<<w8<<p5<<
double(m_eventStat.m_eventsBTE[Endcap])/pa<<
" |\n";
325 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Forward]);
326 out<<
"| For forward region "<<w8<<p5<<
double(m_eventStat.m_eventsBTE[Forward])/pa<<
" |\n";
328 pa = coerceToOne(m_eventStat.m_nclustersNegBP);
329 double ratio =
double(m_eventStat.m_nclustersPosBP)/pa;
331 out<<
"| Ratio barrel pixels clusters for +/- particles ="<<w8<<p5<<
ratio<<
" +-"<<w8<<p5<<
eratio<<
" |\n";
332 pa = coerceToOne(m_eventStat.m_nclustersNegEP);
335 out<<
"| Ratio endcap pixels clusters for +/- particles ="<<w8<<p5<<
ratio<<
" +-"<<w8<<p5<<
eratio<<
" |\n";
336 pa = coerceToOne(m_eventStat.m_nclustersNegBS);
339 out<<
"| Ratio barrel Strip clusters for +/- particles ="<<w8<<p5<<
ratio<<
" +-"<<w8<<p5<<
eratio<<
" |\n";
340 pa = coerceToOne(m_eventStat.m_nclustersNegES);
343 out<<
"| Ratio endcap Strip clusters for +/- particles ="<<w8<<p5<<
ratio<<
" +-"<<w8<<p5<<
eratio<<
" |\n";
344 pa =
double(m_eventStat.m_eventsNEG);
if(pa < 1.) pa = 1.;
347 out<<
"| Number truth particles and +/- ratio ="<<std::setw(10)<<m_eventStat.m_events<<w8<<p5<<
ratio<<
" +-"<<w8<<p5<<
eratio<<
" |\n";
349 if(m_eventStat.m_nclustersPTOT!=0)
ratio =
double(m_eventStat.m_nclustersPTOTt)/
double(m_eventStat.m_nclustersPTOT);
351 out<<
"| Number pix clusters, truth clusters and ratio = "
352 <<std::setw(10)<<m_eventStat.m_nclustersPTOT
353 <<std::setw(10)<<m_eventStat.m_nclustersPTOTt
354 <<std::setw(12)<<std::setprecision(5)<<
ratio<<
" |\n";
356 if(m_eventStat.m_nclustersSTOT!=0)
ratio =
double(m_eventStat.m_nclustersSTOTt)/
double(m_eventStat.m_nclustersSTOT);
357 out<<
"| Number strip clusters, truth clusters and ratio = "
358 <<std::setw(10)<<m_eventStat.m_nclustersSTOT
359 <<std::setw(10)<<m_eventStat.m_nclustersSTOTt
360 <<std::setw(12)<<std::setprecision(5)<<
ratio<<
" |\n";
362 out<<
"|-----------------------------------------------------------------------------------|\n\n";
367 int n = 47-(
t->key().size());
368 std::string s1;
for(
int i=0;
i<
n; ++
i) s1.append(
" "); s1.append(
"|");
371 out<<
"|-----------------------------------------------------------------------------------|\n";
372 out<<
"| Statistic for "<<(
t->key())<<s1<<
"\n";
374 double ne =
double(m_eventStat.m_events);
if(ne < 1.) ne = 1.;
375 double ef [6];
for(
int i=0;
i!=6; ++
i)
ef [
i] =
double(m_trackCollectionStat[
nc].m_efficiency [
i]) /ne;
376 double ef0[6];
for(
int i=0;
i!=6; ++
i) ef0[
i] =
double(m_trackCollectionStat[
nc].m_efficiencyN [
i][0])/ne;
377 double ef1[6];
for(
int i=0;
i!=6; ++
i) ef1[
i] =
double(m_trackCollectionStat[
nc].m_efficiencyN [
i][1])/ne;
378 double ef2[6];
for(
int i=0;
i!=6; ++
i) ef2[
i] =
double(m_trackCollectionStat[
nc].m_efficiencyN [
i][2])/ne;
381 using EffArray_t = std::array<double, 6>;
383 auto makeEffArray = [](
const auto & threeDimArray,
const size_t secondIdx,
const size_t thirdIdx,
const double denom){
386 auto invDenom = 1./
denom;
388 entry = threeDimArray[
idx++][secondIdx][thirdIdx]*invDenom;
393 const auto & efficiencyArrayInput = m_trackCollectionStat[
nc].m_efficiencyBTE;
395 double neBTE = coerceToOne(m_eventStat.m_eventsBTE[Barrel]);
396 const EffArray_t efB0 = makeEffArray(efficiencyArrayInput,0,Barrel,neBTE);
397 const EffArray_t efB1 = makeEffArray(efficiencyArrayInput,1,Barrel,neBTE);
398 const EffArray_t efB2 = makeEffArray(efficiencyArrayInput,2,Barrel,neBTE);
399 const EffArray_t efB3 = makeEffArray(efficiencyArrayInput,3,Barrel,neBTE);
400 const EffArray_t efB4 = makeEffArray(efficiencyArrayInput,4,Barrel,neBTE);
402 neBTE = coerceToOne(m_eventStat.m_eventsBTE[Transition]);
403 const EffArray_t efT0 = makeEffArray(efficiencyArrayInput,0,Transition,neBTE);
404 const EffArray_t efT1 = makeEffArray(efficiencyArrayInput,1,Transition,neBTE);
405 const EffArray_t efT2 = makeEffArray(efficiencyArrayInput,2,Transition,neBTE);
406 const EffArray_t efT3 = makeEffArray(efficiencyArrayInput,3,Transition,neBTE);
407 const EffArray_t efT4 = makeEffArray(efficiencyArrayInput,4,Transition,neBTE);
409 neBTE = coerceToOne(m_eventStat.m_eventsBTE[Endcap]);
410 const EffArray_t efE0 = makeEffArray(efficiencyArrayInput,0,Endcap,neBTE);
411 const EffArray_t efE1 = makeEffArray(efficiencyArrayInput,1,Endcap,neBTE);
412 const EffArray_t efE2 = makeEffArray(efficiencyArrayInput,2,Endcap,neBTE);
413 const EffArray_t efE3 = makeEffArray(efficiencyArrayInput,3,Endcap,neBTE);
414 const EffArray_t efE4 = makeEffArray(efficiencyArrayInput,4,Endcap,neBTE);
416 neBTE = coerceToOne(m_eventStat.m_eventsBTE[Forward]);
417 const EffArray_t efD0 = makeEffArray(efficiencyArrayInput,0,Forward,neBTE);
418 const EffArray_t efD1 = makeEffArray(efficiencyArrayInput,1,Forward,neBTE);
419 const EffArray_t efD2 = makeEffArray(efficiencyArrayInput,2,Forward,neBTE);
420 const EffArray_t efD3 = makeEffArray(efficiencyArrayInput,3,Forward,neBTE);
421 const EffArray_t efD4 = makeEffArray(efficiencyArrayInput,4,Forward,neBTE);
424 double efrec = ef0[0]+ef0[1]+ef0[2]+ef1[0]+ef1[1]+ef2[0];
425 double efrecB = efB0[0]+efB0[1]+efB0[2]+efB1[0]+efB1[1]+efB2[0];
426 double efrecT = efT0[0]+efT0[1]+efT0[2]+efT1[0]+efT1[1]+efT2[0];
427 double efrecE = efE0[0]+efE0[1]+efE0[2]+efE1[0]+efE1[1]+efE2[0];
428 double efrecD = efD0[0]+efD0[1]+efD0[2]+efD1[0]+efD1[1]+efD2[0];
430 ne = coerceToOne(m_eventStat.m_eventsPOS);
431 double efP[6];
for(
int i=0;
i!=6; ++
i) efP[
i] =
double(m_trackCollectionStat[
nc].m_efficiencyPOS[
i])/ne;
432 ne = coerceToOne(m_eventStat.m_eventsNEG);
433 double efN[6];
for(
int i=0;
i!=6; ++
i) efN[
i] =
double(m_trackCollectionStat[
nc].m_efficiencyNEG[
i])/ne;
435 out<<
"|-----------------------------------------------------------------------------------|\n";
436 out<<
"| Probability to lose 0 1 2 3 4 >=5 clusters |\n";
437 out<<
"|-----------------------------------------------------------------------------------|\n";
439 auto formattedOutput=[&
out](
auto & effArray){
440 for (
size_t i{};
i!=6;++
i){
441 out<<std::setw(9)<<std::setprecision(4)<<effArray[
i];
446 out<<
"| For all particles ";
448 out<<
"| For + particles ";
449 formattedOutput(efP);
450 out<<
"| For - particles ";
451 formattedOutput(efN);
452 out<<
"|-----------------------------------------------------------------------------------|\n";
453 out<<
"| Barrel region |\n";
454 out<<
"| 0 wrong clusters ";
455 formattedOutput(efB0);
456 out<<
"| 1 wrong clusters ";
457 formattedOutput(efB1);
458 out<<
"| 2 wrong clusters ";
459 formattedOutput(efB2);
460 out<<
"| 3 wrong clusters ";
461 formattedOutput(efB3);
462 out<<
"| >=4 wrong clusters ";
463 formattedOutput(efB4);
464 out<<
"|-----------------------------------------------------------------------------------|\n";
465 out<<
"| Transition region |\n";
466 out<<
"| 0 wrong clusters ";
467 formattedOutput(efT0);
468 out<<
"| 1 wrong clusters ";
469 formattedOutput(efT1);
470 out<<
"| 2 wrong clusters ";
471 formattedOutput(efT2);
472 out<<
"| 3 wrong clusters ";
473 formattedOutput(efT3);
474 out<<
"| >=4 wrong clusters ";
475 formattedOutput(efT4);
476 out<<
"|-----------------------------------------------------------------------------------|\n";
477 out<<
"| Endcap region |\n";
478 out<<
"| 0 wrong clusters ";
479 formattedOutput(efE0);
480 out<<
"| 1 wrong clusters ";
481 formattedOutput(efE1);
482 out<<
"| 2 wrong clusters ";
483 formattedOutput(efE2);
484 out<<
"| 3 wrong clusters ";
485 formattedOutput(efE3);
486 out<<
"| >=4 wrong clusters ";
487 formattedOutput(efE4);
488 out<<
"|-----------------------------------------------------------------------------------|\n";
489 out<<
"| Forward region |\n";
490 out<<
"| 0 wrong clusters ";
491 formattedOutput(efD0);
492 out<<
"| 1 wrong clusters ";
493 formattedOutput(efD1);
494 out<<
"| 2 wrong clusters ";
495 formattedOutput(efD2);
496 out<<
"| 3 wrong clusters ";
497 formattedOutput(efD3);
498 out<<
"| >=4 wrong clusters ";
499 formattedOutput(efD4);
501 out<<
"|-----------------------------------------------------------------------------------|\n";
502 pa = coerceToOne(m_eventStat.m_particleClusters[0]);
503 out<<
"| Efficiency reconstruction (number lose+wrong < 3) = "
504 <<std::setw(9)<<std::setprecision(5)<<efrec
506 <<std::setw(9)<<std::setprecision(5)<<efrec*
double(m_eventStat.m_events)/pa
509 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Barrel]);
510 out<<
"| For barrel region = "
511 <<std::setw(9)<<std::setprecision(5)<<efrecB
513 <<std::setw(9)<<std::setprecision(5)<<efrecB*
double(m_eventStat.m_eventsBTE[Barrel])/pa
516 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Transition]);
517 out<<
"| For transition region = "
518 <<std::setw(9)<<std::setprecision(5)<<efrecT
520 <<std::setw(9)<<std::setprecision(5)<<efrecT*
double(m_eventStat.m_eventsBTE[Transition])/pa
523 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Endcap]);
524 out<<
"| For endcap region = "
525 <<std::setw(9)<<std::setprecision(5)<<efrecE
527 <<std::setw(9)<<std::setprecision(5)<<efrecE*
double(m_eventStat.m_eventsBTE[Endcap])/pa
530 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Forward]);
531 out<<
"| For forward region = "
532 <<std::setw(9)<<std::setprecision(5)<<efrecD
534 <<std::setw(9)<<std::setprecision(5)<<efrecD*
double(m_eventStat.m_eventsBTE[Forward])/pa
538 out<<
"|-----------------------------------------------------------------------------------|\n";
539 out<<
"| Reconstructed tracks + - +/-ratio error |\n";
540 out<<
"|-----------------------------------------------------------------------------------|\n";
542 pa = coerceToOne(m_trackCollectionStat[
nc].m_ntracksNEGB);
547 <<std::setw(10)<<m_trackCollectionStat[
nc].m_ntracksPOSB
548 <<std::setw(11)<<m_trackCollectionStat[
nc].m_ntracksNEGB
549 <<std::setw(11)<<std::setprecision(5)<<
ratio
550 <<std::setw(11)<<std::setprecision(5)<<
eratio<<
" |\n";
551 pa = coerceToOne(m_trackCollectionStat[
nc].m_ntracksNEGE);
556 <<std::setw(10)<<m_trackCollectionStat[
nc].m_ntracksPOSE
557 <<std::setw(11)<<m_trackCollectionStat[
nc].m_ntracksNEGE
558 <<std::setw(11)<<std::setprecision(5)<<
ratio
559 <<std::setw(11)<<std::setprecision(5)<<
eratio<<
" |\n";
560 pa = coerceToOne(m_trackCollectionStat[
nc].m_ntracksNEGFWD);
561 ratio =
double(m_trackCollectionStat[
nc].m_ntracksPOSFWD)/pa;
565 <<std::setw(10)<<m_trackCollectionStat[
nc].m_ntracksPOSFWD
566 <<std::setw(11)<<m_trackCollectionStat[
nc].m_ntracksNEGFWD
567 <<std::setw(11)<<std::setprecision(5)<<
ratio
568 <<std::setw(11)<<std::setprecision(5)<<
eratio<<
" |\n";
575 for(
int k = 0;
k!=50; ++
k) {
576 nt+=m_trackCollectionStat[
nc].m_total[
k];
577 ft+=m_trackCollectionStat[
nc].m_fake [
k];
578 if(!kf &&
nt) kf =
k;
583 out<<
"|-----------------------------------------------------------------------------------|\n";
584 out<<
"| Fake tracks rate for different number of clusters on track |\n";
585 out<<
"|-----------------------------------------------------------------------------------|\n";
587 for(
int k = kf;
k!=kf+6; ++
k) {
588 out<<
"| >= "<<std::setw(2)<<
k<<
" ";
592 for(
int k = kf;
k!=kf+6; ++
k) {
594 out<<
"|"<<std::setw(12)<<std::setprecision(5)<<
eff<<
" ";
595 nt-=m_trackCollectionStat[
nc].m_total[
k];
596 ft-=m_trackCollectionStat[
nc].m_fake [
k];
599 out<<
"|-----------------------------------------------------------------------------------|\n";
604 return StatusCode::SUCCESS;
617 out<<
"|--------------------------------------------------------------------------------------------------------------------|\n";
619 n = 65-
t->key().size();
620 std::string s1;
for(
int i=0;
i<
n; ++
i) s1.append(
" "); s1.append(
"|");
621 out<<
"| Location of input tracks | "<<
t->key()<<s1<<
"\n";
623 auto padString = [](
const std::string &
s){
624 const int n = 65 -
s.size();
625 return s + std::string(
n,
' ');
627 std::string
s2 = padString(m_spacepointsPixelname.key());
628 std::string
s3 = padString(m_spacepointsStripname.key());
629 std::string
s4 = padString(m_spacepointsOverlapname.key());
630 std::string s5 = padString(m_clustersPixelname.key());
631 std::string s6 = padString(m_clustersStripname.key());
633 std::string s7 = padString(m_truth_locationPixel.key());
634 std::string s8 = padString(m_truth_locationStrip.key());
636 out<<
"| Pixel space points | "<<
s2<<
"|\n";
637 out<<
"| Strip space points | "<<
s3<<
"|\n";
638 out<<
"| Overlap space points | "<<
s4<<
"|\n";
639 out<<
"| Pixel clusters | "<<s5<<
"|\n";
640 out<<
"| Strip clusters | "<<s6<<
"|\n";
641 out<<
"| Truth location for pixels | "<<s7<<
"|\n";
642 out<<
"| Truth location for strips | "<<s8<<
"|\n";
643 out<<
"| max pT cut | "
644 <<std::setw(14)<<std::setprecision(5)<<m_ptcutmax
646 out<<
"| rapidity cut | "
647 <<std::setw(14)<<std::setprecision(5)<<m_rapcut
649 out<<
"| min Radius | "
650 <<std::setw(14)<<std::setprecision(5)<<m_rmin
652 out<<
"| max Radius | "
653 <<std::setw(14)<<std::setprecision(5)<<m_rmax
655 out<<
"| Min. number sp.points for generated track | "
656 <<std::setw(14)<<std::setprecision(5)<<m_spcut
658 out<<
"|--------------------------------------------------------------------------------------------------------------------|\n";
669 auto formatOutput = [&
out](
const auto val){
673 out<<
"|---------------------------------------------------------------------|\n";
674 out<<
"| m_nspacepoints | ";
676 out<<
"| m_nclusters | ";
678 out<<
"| Kine-Clusters size | ";
680 out<<
"| Kine-SpacePoints size | ";
682 out<<
"| Number good kine tracks | ";
684 out<<
"|---------------------------------------------------------------------|\n";
695 std::lock_guard<std::mutex> lock(m_statMutex);
699 std::unique_ptr<SG::ReadHandle<InDet::SiClusterContainer> > pixelcontainer;
700 std::unique_ptr<SG::ReadHandle<InDet::SiClusterContainer> > stripcontainer;
703 pixelcontainer = std::make_unique<SG::ReadHandle<InDet::SiClusterContainer> >(m_clustersPixelname,ctx);
704 if (!pixelcontainer->isValid())
ATH_MSG_DEBUG(
"Failed to create Pixel clusters container read handle with key " << m_clustersPixelname.key());
710 stripcontainer = std::make_unique<SG::ReadHandle<InDet::SiClusterContainer> >(m_clustersStripname,ctx);
711 if (!stripcontainer->isValid())
ATH_MSG_DEBUG(
"Failed to create Strip clusters container read handle with key " << m_clustersStripname.key());
719 if(pixelcontainer && pixelcontainer->isValid()) {
720 InDet::SiClusterContainer::const_iterator
w = (*pixelcontainer)->begin();
721 InDet::SiClusterContainer::const_iterator we = (*pixelcontainer)->end ();
725 InDet::SiClusterCollection::const_iterator
c = (*w)->begin();
726 InDet::SiClusterCollection::const_iterator ce = (*w)->end ();
731 ++m_eventStat.m_nclustersPTOT;
732 if(isTruth(event_data,(*
c))) ++m_eventStat.m_nclustersPTOTt;
735 int nk = kine(event_data,(*
c),Kine,999);
736 for(
int i=0;
i!=nk; ++
i) {
737 if(!isTheSameDetElement(event_data,Kine[
i],(*
c))) {
749 if(stripcontainer && stripcontainer->isValid()) {
750 InDet::SiClusterContainer::const_iterator
w = (*stripcontainer)->begin();
751 InDet::SiClusterContainer::const_iterator we = (*stripcontainer)->end ();
755 InDet::SiClusterCollection::const_iterator
c = (*w)->begin();
756 InDet::SiClusterCollection::const_iterator ce = (*w)->end ();
761 ++m_eventStat.m_nclustersSTOT;
762 if(isTruth(event_data,(*
c))) ++m_eventStat.m_nclustersSTOTt;
764 int nk = kine(event_data,(*
c),Kine,999);
765 for(
int i=0;
i!=nk; ++
i) {
766 if(!isTheSameDetElement(event_data,Kine[
i],(*
c))) event_data.
m_kinecluster.insert(std::make_pair(Kine[
i],(*
c)));
784 if(m_usePix && !m_spacepointsPixelname.key().empty()) {
787 ATH_MSG_DEBUG(
"Invalid Pixels space points container read handle for key " << m_spacepointsPixelname.key() );
792 for(; spc != spce; ++spc) {
796 for(; sp != spe; ++sp) {
799 int nk = kine(event_data,(*sp)->clusterList().first,Kine,999);
800 for(
int i=0;
i!=nk; ++
i) {
802 if(!isTheSameDetElement(event_data,Kine[
i],(*sp))) {
813 if(m_useStrip && !m_spacepointsStripname.key().empty()) {
816 ATH_MSG_DEBUG(
"Invalid Strip space points container read handle for key " << m_spacepointsStripname.key() );
822 for(; spc != spce; ++spc) {
827 for(; sp != spe; ++sp) {
831 int nk = kine(event_data,(*sp)->clusterList().first,(*sp)->clusterList().second,Kine,999);
832 for(
int i=0;
i!=nk; ++
i) {
833 if(!isTheSameDetElement(event_data,Kine[
i],(*sp))) {
844 if(m_useStrip && !m_spacepointsOverlapname.key().empty()) {
845 event_data.
m_spacepointsOverlap=std::make_unique< SG::ReadHandle<SpacePointOverlapCollection> >(m_spacepointsOverlapname,ctx);
847 ATH_MSG_DEBUG(
"Invalid overlap space points container read handle for key " << m_spacepointsOverlapname.key() );
853 for (; sp!=spe; ++sp) {
856 int nk = kine(event_data,(*sp)->clusterList().first,(*sp)->clusterList().second,Kine,999);
857 for(
int i=0;
i!=nk; ++
i) {
858 if(!isTheSameDetElement(event_data,Kine[
i],(*sp))) {
883 std::list<int> worskine;
890 unsigned int nc = 1 ;
892 auto coerceTo49 = [] (
const size_t idx){
893 return (
idx<50) ?
idx : 49;
898 if((*c).first==
k0) {++
nc;
continue;}
901 const size_t clusterIdx =coerceTo49(
nc);
906 const size_t spacepointIdx =coerceTo49(
ns);
910 if (
nc < minclusters(
eta) ) worskine.push_back(
k0);
927 if (
nc < minclusters(
eta) ) worskine.push_back(
k0);
932 for(
auto & pThisCluster: worskine) {
942 if (not de)
continue;
984 if(++
nc >= 100)
return;
996 int KINE[200],NKINE[200];
1001 se = (*t)->trackStateOnSurfaces()->end();
1010 const AmgVector(5)& Vpf = tpf ->parameters ();
1011 double pTf = std::abs(
std::sin(Vpf[3])/Vpf[4]);
1012 double etaf = std::abs(
log(
tan(.5*Vpf[3])));
1013 bool qTf = pTf > minpT(etaf);
1025 double minpt = minpT(rap);
1026 if (
pT > minpt &&
pT < m_ptcutmax) {
1031 else if(pT < -minpt && pT > -m_ptcutmax) {
1057 int Kine[1000], nk=kine0(event_data,rd,Kine,999); ++
NC;
if(!nk) ++N0;
1059 for(
int k = 0;
k!=nk; ++
k) {
1062 for(;
n!=NK; ++
n) {
if(Kine[
k]==KINE[
n]) {++NKINE[
n];
break;}}
1063 if(
n==NK) {KINE[NK] = Kine[
k]; NKINE[NK] = 1;
if (NK < 200) ++NK;}
1065 for(
int n=0;
n!=NK; ++
n) {
if(NKINE[
n]>nkm) nkm = NKINE[
n];}
1068 for(
int n=0;
n!=NK; ++
n) {
1070 int NQ = 1000*NKINE[
n]+(
NC-NKINE[
n]);
1072 event_data.
m_tracks[
nc].insert(std::make_pair(KINE[
n],NQ));
1096 std::multimap<int,int>::const_iterator
t, te = event_data.
m_tracks[
nc].end();
1100 int k = (*p).barcode();
1106 if((*t).first!=
k)
break;
1107 int ts = (*t).second/1000;
1108 int ws = (*t).second%1000;
1109 if (
ts >
m ) {
m =
ts;
w = ws;}
1110 else if(
ts==
m &&
w > ws) {
w = ws;}
1112 int d =
n-
m;
if(
d<0)
d = 0;
else if(
d > 5)
d=5;
if(
w>4)
w = 4;
1117 int ch = (*p).charge();
1140 int Kine1[1000],Kine2[1000];
1141 int n1 = kine(event_data,
d1,Kine1,
nmax);
if(!
n1)
return nkine;
1142 int n2 = kine(event_data,
d2,Kine2,
nmax);
if(!n2)
return nkine;
1144 for(
int i = 0;
i!=
n1; ++
i) {
1145 for(
int j = 0; j!=n2; ++j) {
1146 if(Kine1[
i]==Kine2[j]) {Kine[nkine++] = Kine1[
i];
break;}
1160 PRD_MultiTruthCollection::const_iterator mce;
1161 PRD_MultiTruthCollection::const_iterator
mc = findTruth(event_data,
d,mce);
1166 for(;
mc!=mce; ++
mc) {
1168 if( (*mc).first !=
ID )
return nkine;
1170 int k = (*mc).second.barcode();
if(
k<=0)
continue;
1173 if(!pa or !pa->production_vertex())
continue;
1175 int pdg = std::abs(pa->pdg_id());
if(m_pdg && m_pdg != pdg )
continue;
1177 const HepPDT::ParticleData* pd = m_particleDataTable->particle(pdg);
1178 if(!pd or std::abs(pd->charge()) < .5)
continue;
1182 double px = pa->momentum().px();
1183 double py = pa->momentum().py();
1184 double pz = pa->momentum().pz();
1186 if( pt < m_ptcut || pt > m_ptcutmax)
continue;
1190 double t = std::abs(
pz)/
pt;
1191 if(
t > m_tcut )
continue;
1195 double vx = pa->production_vertex()->position().x();
1196 double vy = pa->production_vertex()->position().y();
1197 double r = std::sqrt(vx*vx+vy*vy);
1198 if( r < m_rmin || r > m_rmax)
continue;
1200 Kine[nkine] =
k;
if(++nkine >=
nmax)
break;
1213 PRD_MultiTruthCollection::const_iterator mce;
1214 PRD_MultiTruthCollection::const_iterator
mc = findTruth(event_data,
d,mce);
1219 for(;
mc!=mce; ++
mc) {
1221 if( (*mc).first !=
ID )
return nkine;
1223 int k = (*mc).second.barcode();
if(
k<=0)
continue;
1224 Kine[nkine] =
k;
if(++nkine >=
nmax)
break;
1236 PRD_MultiTruthCollection::const_iterator mce;
1237 PRD_MultiTruthCollection::const_iterator
mc = findTruth(event_data,
d,mce);
1248 std::multimap<int,const Trk::PrepRawData*>::const_iterator
k = event_data.
m_kinecluster.find(K);
1251 if((*k).first!= K)
return false;
1252 if(
d->detectorElement()==(*k).second->detectorElement())
return true;
1267 std::multimap<int,const Trk::SpacePoint*>::const_iterator
k = event_data.
m_kinespacepoint.find(K);
1272 if((*k).first!= K)
return false;
1285 if((*k).first!= K)
return false;
1311 std::list<int>::const_iterator dif = event_data.
m_difference[
nc].begin();
1313 std::multimap<int,const Trk::PrepRawData*>::const_iterator
c,ce = event_data.
m_kinecluster.end();
1315 int n = 69-m_tracklocation[
nc].key().size();
1316 std::string s1;
for(
int i=0;
i<
n; ++
i) s1.append(
" "); s1.append(
"|");
1317 std::stringstream
out;
1318 out<<
"|----------------------------------------------------------------------------------------|\n";
1319 out<<
"| "<<m_tracklocation[
nc]<<s1<<
"\n";
1320 out<<
"|----------------------------------------------------------------------------------------|\n";
1321 out<<
"| # pdg kine Ncl Ntr Nsp Lose pT(MeV) rapidity radius z |\n";
1322 out<<
"|----------------------------------------------------------------------------------------|\n";
1327 int k = (*p).barcode();
1332 PRD_MultiTruthCollection::const_iterator mce;
1333 PRD_MultiTruthCollection::const_iterator
mc = findTruth(event_data,
d,mce);
1337 for(;
mc!=mce; ++
mc) {
1338 if((*mc).first !=
ID)
break;
1339 if((*mc).second.barcode()==
k) {Q=
true;
break;}
1346 double px = pa->momentum().px();
1347 double py = pa->momentum().py();
1348 double pz = pa->momentum().pz();
1349 double vx = pa->production_vertex()->position().x();
1350 double vy = pa->production_vertex()->position().y();
1351 double vz = pa->production_vertex()->position().z();
1353 double t = std::atan2(
pt,
pz);
1355 double r = std::sqrt(vx*vx+vy*vy);
1359 <<std::setw(6)<<pa->pdg_id()
1363 <<std::setw(4)<<(*dif)
1364 <<std::setw(12)<<std::setprecision(5)<<
pt
1365 <<std::setw(12)<<std::setprecision(5)<<ra
1366 <<std::setw(12)<<std::setprecision(5)<<
r
1367 <<std::setw(12)<<std::setprecision(5)<<vz
1372 out<<
"|----------------------------------------------------------------------------------------|\n";
1382 PRD_MultiTruthCollection::const_iterator
1386 PRD_MultiTruthCollection::const_iterator& mce)
1391 PRD_MultiTruthCollection::const_iterator
mc;
1397 for (
int i=0;
i<3;
i++) {
1399 mce=truth[
i]->end();
1400 return truth[
i]->end();
1403 throw std::runtime_error(
"Neither Pixel nor Strip truth.");
1416 PRD_MultiTruthCollection::const_iterator mce;
1417 PRD_MultiTruthCollection::const_iterator
mc = findTruth(event_data,
d,mce);
1419 for(;
mc!=mce; ++
mc) {
1420 if((*mc).second.barcode()==
k) {
1425 double px =
pat->momentum().px();
1426 double py =
pat->momentum().py();
1427 double pz =
pat->momentum().pz();
1429 double t = std::atan2(
pt,
pz) ;
1436 eta > 1.6 ? rap = 2 :
eta > .8 ? rap = 1 : rap = 0;
1438 int pdg =
pat->pdg_id();
1439 const HepPDT::ParticleData* pd = m_particleDataTable->particle(abs(pdg));
1441 double ch = pd->charge();
if(pdg < 0)
ch = -
ch;
1442 if(
ch > .5)
return 1;
1443 if(
ch < -.5)
return -1;