5 #include "GaudiKernel/IPartPropSvc.h"
6 #include "GaudiKernel/ServiceHandle.h"
13 #include "HepPDT/ParticleDataTable.hh"
24 (
const std::string&
name,ISvcLocator* pSvcLocator)
26 m_spacepointsSCTname(
"SCT_SpacePoints"),
27 m_spacepointsPixelname(
"PixelSpacePoints"),
28 m_spacepointsOverlapname(
"OverlapSpacePoints"),
29 m_clustersSCTname(
"SCT_Clusters"),
30 m_clustersPixelname(
"PixelClusters"),
31 m_clustersTRTname(
"TRT_DriftCircles"),
32 m_truth_locationPixel(
"PRD_MultiTruthPixel" ),
33 m_truth_locationSCT(
"PRD_MultiTruthSCT" ),
34 m_truth_locationTRT(
"PRD_MultiTruthTRT" )
107 IPartPropSvc* partPropSvc =
nullptr;
108 sc = service(
"PartPropSvc", partPropSvc,
true);
109 if (
sc.isFailure()) {
111 return StatusCode::FAILURE;
119 return StatusCode::FAILURE;
170 std::vector<SG::ReadHandle<PRD_MultiTruthCollection> > read_handle;
171 read_handle.reserve(3);
174 if (not read_handle.back().isValid()) {
176 return StatusCode::FAILURE;
178 event_data.
m_truthPIX = &(*read_handle.back());
183 if (not read_handle.back().isValid()) {
185 return StatusCode::FAILURE;
187 event_data.
m_truthSCT = &(*read_handle.back());
192 if (not read_handle.back().isValid()) {
194 return StatusCode::FAILURE;
196 event_data.
m_truthTRT = &(*read_handle.back());
213 for (
unsigned int i=0;
i< m_trackCollectionStat.size(); ++
i ) {
223 return StatusCode::SUCCESS;
231 if(m_eventStat.m_events<=0)
return StatusCode::SUCCESS;
232 const auto & w13 = std::setw(13);
233 const auto & p5 = std::setprecision(5);
234 const auto topNtail=[](
const std::string &
str){
return "|" +
str +
"|";};
235 const std::string lineSeparator(83,
'-');
236 const std::string spaceSeparator(83,
' ');
237 std::stringstream
out;
239 out<<topNtail(lineSeparator)<<
"\n";
240 out<<
"| TrackClusterAssValidation statistic for charge truth particles with |\n";
241 out<<topNtail(spaceSeparator)<<
"\n";
242 out<<
"| pT >="<<w13<<p5<<
m_ptcut<<
" MeV"<<
" |\n";
247 out<<
"| max vertex radius <="<<w13<<p5<<
m_rmax<<
" mm"<<
" |\n";
248 out<<
"| min vertex radius >="<<w13<<p5<<
m_rmin<<
" mm"<<
" |\n";
249 out<<
"| particles pdg ="<<std::setw(8)<<
m_pdg<<
" |\n";
251 auto yesNo=[](
const bool predicate){
return predicate ?
"yes" :
"no ";};
252 out<<
"| use Pixels information "<<yesNo(
m_usePIX)<<
" |\n";
253 out<<
"| use SCT information "<<yesNo(
m_useSCT)<<
" |\n";
254 out<<
"| use TRT information "<<yesNo(
m_useTRT)<<
" |\n";
256 out<<topNtail(spaceSeparator)<<
"\n";
258 enum Regions{Barrel, Transition, Endcap, Forward, NRegions};
259 auto incrementArray=[](
auto &
array,
const int idx){
for (
int j{};j!=NRegions;++j)
array[
idx][j] +=
array[
idx+1][j];};
260 for(
int i=48;
i>=0; --
i) {
261 m_eventStat.m_particleClusters [
i] +=m_eventStat.m_particleClusters [
i+1];
262 incrementArray(m_eventStat.m_particleClustersBTE,
i);
263 m_eventStat.m_particleSpacePoints [
i] +=m_eventStat.m_particleSpacePoints [
i+1];
264 incrementArray(m_eventStat.m_particleSpacePointsBTE,
i);
266 auto coerceToOne=[](
const double & initialVal)->
double {
return (initialVal<1.) ? 1. : initialVal; };
273 double pa = coerceToOne(m_eventStat.m_particleClusters[0]);
274 std::array<double, 10> pc2ff{};
275 size_t clusterIdx = 2;
276 for (
auto & thisCluster: pc2ff){
277 thisCluster =
double(m_eventStat.m_particleClusters[ clusterIdx++ ])/ pa;
280 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Barrel]);
281 std::array<double, 10> pcBarrel2ff{};
282 size_t clusterBarrelIdx = 2;
283 for (
auto & thisClusterB: pcBarrel2ff){
284 thisClusterB =
double(m_eventStat.m_particleClustersBTE[ clusterBarrelIdx++ ][Barrel])/ pa;
287 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Transition]);
288 std::array<double, 10> pcTransition2ff{};
289 size_t clusterTransitionIdx = 2;
290 for (
auto & thisClusterT: pcTransition2ff){
291 thisClusterT =
double(m_eventStat.m_particleClustersBTE[ clusterTransitionIdx++ ][Transition])/ pa;
294 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Endcap]);
295 std::array<double, 10> pcEndcap2ff{};
296 size_t clusterEndcapIdx = 2;
297 for (
auto & thisClusterE: pcEndcap2ff){
298 thisClusterE =
double(m_eventStat.m_particleClustersBTE[ clusterEndcapIdx++ ][Endcap])/ pa;
301 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][3]);
302 std::array<double, 10> pcFwd2ff{};
303 size_t clusterFwdIdx = 2;
304 for (
auto & thisClusterD: pcFwd2ff){
305 thisClusterD =
double(m_eventStat.m_particleClustersBTE[ clusterFwdIdx++ ][Forward])/ pa;
311 pa = coerceToOne(m_eventStat.m_particleSpacePoints[0]);
312 std::array<double, 10> sp2ff{};
313 size_t spacepointIdx = 2;
314 for (
auto & thisSpacepoint: sp2ff){
315 thisSpacepoint =
double(m_eventStat.m_particleSpacePoints[ spacepointIdx++ ])/ pa;
318 pa = coerceToOne(m_eventStat.m_particleSpacePointsBTE[0][Barrel]);
319 std::array<double, 10> spBarrel2ff{};
320 size_t spacepointBarrelIdx = 2;
321 for (
auto & thisSpacepoint: spBarrel2ff){
322 thisSpacepoint =
double(m_eventStat.m_particleSpacePointsBTE[ spacepointBarrelIdx++ ][Barrel])/ pa;
325 pa = coerceToOne(m_eventStat.m_particleSpacePointsBTE[0][Transition]);
326 std::array<double, 10> spTransition2ff{};
327 size_t spacepointTransitionIdx = 2;
328 for (
auto & thisSpacepoint: spTransition2ff){
329 thisSpacepoint =
double(m_eventStat.m_particleSpacePointsBTE[ spacepointTransitionIdx++ ][Transition])/ pa;
332 pa = coerceToOne(m_eventStat.m_particleSpacePointsBTE[0][Endcap]);
333 std::array<double, 10> spEndcap2ff{};
334 size_t spacepointEndcapIdx = 2;
335 for (
auto & thisSpacepoint: spEndcap2ff){
336 thisSpacepoint =
double(m_eventStat.m_particleSpacePointsBTE[ spacepointEndcapIdx++ ][Endcap])/ pa;
339 pa = coerceToOne(m_eventStat.m_particleSpacePointsBTE[0][Forward]);
340 std::array<double, 10> spFwd2ff{};
341 size_t spacepointFwdIdx = 2;
342 for (
auto & thisSpacepoint: spFwd2ff){
343 thisSpacepoint =
double(m_eventStat.m_particleSpacePointsBTE[ spacepointFwdIdx++ ][Forward])/ pa;
345 auto w8=std::setw(8);
346 out<<
"| Probability for such charge particles to have some number silicon |\n";
347 out<<
"| clusters | space points |\n";
348 out<<
"| Total Barrel Transi Endcap Forward | Total Barrel Transi Endcap Forward |\n";
351 out<<
"| >= "<<
idx+2<< std::string((
idx<8)?
" ":
" ")
353 <<w8<<p5<<pcBarrel2ff[
idx]
354 <<w8<<p5<<pcTransition2ff[
idx]
355 <<w8<<p5<<pcEndcap2ff[
idx]
356 <<w8<<p5<<pcFwd2ff[
idx]<<
" | "
359 <<w8<<p5<<spBarrel2ff[
idx]
360 <<w8<<p5<<spTransition2ff[
idx]
361 <<w8<<p5<<spEndcap2ff[
idx]
362 <<w8<<p5<<spFwd2ff[
idx]
366 out<<topNtail(spaceSeparator)<<
"\n";
367 out<<
"| Additional cuts for truth particles are |\n";
368 out<<
"| number silicon clusters >="<<w13<<
m_clcut<<
" |\n";
370 out<<
"| number space points >="<<w13<<
m_spcut<<
" |\n";
372 pa = coerceToOne(m_eventStat.m_particleClusters[0]);
373 out<<
"| Probability find truth particles with this cuts is "<<w8<<p5<<
double(m_eventStat.m_events)/pa<<
" |\n";
374 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Barrel]);
375 out<<
"| For barrel region "<<w8<<p5<<
double(m_eventStat.m_eventsBTE[Barrel])/pa<<
" |\n";
376 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Transition]);
377 out<<
"| For transition region "<<w8<<p5<<
double(m_eventStat.m_eventsBTE[Transition])/pa<<
" |\n";
378 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Endcap]);
379 out<<
"| For endcap region "<<w8<<p5<<
double(m_eventStat.m_eventsBTE[Endcap])/pa<<
" |\n";
380 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Forward]);
381 out<<
"| For forward region "<<w8<<p5<<
double(m_eventStat.m_eventsBTE[Forward])/pa<<
" |\n";
383 pa = coerceToOne(m_eventStat.m_nclustersNegBP);
384 double ratio =
double(m_eventStat.m_nclustersPosBP)/pa;
386 out<<
"| Ratio barrel pixels clusters for +/- particles ="<<w8<<p5<<
ratio<<
" +-"<<w8<<p5<<
eratio<<
" |\n";
387 pa = coerceToOne(m_eventStat.m_nclustersNegEP);
390 out<<
"| Ratio endcap pixels clusters for +/- particles ="<<w8<<p5<<
ratio<<
" +-"<<w8<<p5<<
eratio<<
" |\n";
391 pa = coerceToOne(m_eventStat.m_nclustersNegBS);
394 out<<
"| Ratio barrel SCT clusters for +/- particles ="<<w8<<p5<<
ratio<<
" +-"<<w8<<p5<<
eratio<<
" |\n";
395 pa = coerceToOne(m_eventStat.m_nclustersNegES);
398 out<<
"| Ratio endcap SCT clusters for +/- particles ="<<w8<<p5<<
ratio<<
" +-"<<w8<<p5<<
eratio<<
" |\n";
399 pa =
double(m_eventStat.m_eventsNEG);
if(pa < 1.) pa = 1.;
402 out<<
"| Number truth particles and +/- ratio ="<<std::setw(10)<<m_eventStat.m_events<<w8<<p5<<
ratio<<
" +-"<<w8<<p5<<
eratio<<
" |\n";
404 if(m_eventStat.m_nclustersPTOT!=0)
ratio =
double(m_eventStat.m_nclustersPTOTt)/
double(m_eventStat.m_nclustersPTOT);
406 out<<
"| Number pix clusters, truth clusters and ratio = "
407 <<std::setw(10)<<m_eventStat.m_nclustersPTOT
408 <<std::setw(10)<<m_eventStat.m_nclustersPTOTt
409 <<std::setw(12)<<std::setprecision(5)<<
ratio<<
" |\n";
411 if(m_eventStat.m_nclustersSTOT!=0)
ratio =
double(m_eventStat.m_nclustersSTOTt)/
double(m_eventStat.m_nclustersSTOT);
412 out<<
"| Number sct clusters, truth clusters and ratio = "
413 <<std::setw(10)<<m_eventStat.m_nclustersSTOT
414 <<std::setw(10)<<m_eventStat.m_nclustersSTOTt
415 <<std::setw(12)<<std::setprecision(5)<<
ratio<<
" |\n";
417 out<<
"|-----------------------------------------------------------------------------------|\n\n";
422 int n = 47-(
t->key().size());
423 std::string s1;
for(
int i=0;
i<
n; ++
i) s1.append(
" "); s1.append(
"|");
426 out<<
"|-----------------------------------------------------------------------------------|\n";
427 out<<
"| Statistic for "<<(
t->key())<<s1<<
"\n";
429 double ne =
double(m_eventStat.m_events);
if(ne < 1.) ne = 1.;
430 double ef [6];
for(
int i=0;
i!=6; ++
i)
ef [
i] =
double(m_trackCollectionStat[
nc].m_efficiency [
i]) /ne;
431 double ef0[6];
for(
int i=0;
i!=6; ++
i) ef0[
i] =
double(m_trackCollectionStat[
nc].m_efficiencyN [
i][0])/ne;
432 double ef1[6];
for(
int i=0;
i!=6; ++
i) ef1[
i] =
double(m_trackCollectionStat[
nc].m_efficiencyN [
i][1])/ne;
433 double ef2[6];
for(
int i=0;
i!=6; ++
i) ef2[
i] =
double(m_trackCollectionStat[
nc].m_efficiencyN [
i][2])/ne;
436 using EffArray_t = std::array<double, 6>;
438 auto makeEffArray = [](
const auto & threeDimArray,
const size_t secondIdx,
const size_t thirdIdx,
const double denom){
441 auto invDenom = 1./
denom;
443 entry = threeDimArray[
idx++][secondIdx][thirdIdx]*invDenom;
448 const auto & efficiencyArrayInput = m_trackCollectionStat[
nc].m_efficiencyBTE;
450 double neBTE = coerceToOne(m_eventStat.m_eventsBTE[Barrel]);
451 const EffArray_t efB0 = makeEffArray(efficiencyArrayInput,0,Barrel,neBTE);
452 const EffArray_t efB1 = makeEffArray(efficiencyArrayInput,1,Barrel,neBTE);
453 const EffArray_t efB2 = makeEffArray(efficiencyArrayInput,2,Barrel,neBTE);
454 const EffArray_t efB3 = makeEffArray(efficiencyArrayInput,3,Barrel,neBTE);
455 const EffArray_t efB4 = makeEffArray(efficiencyArrayInput,4,Barrel,neBTE);
457 neBTE = coerceToOne(m_eventStat.m_eventsBTE[Transition]);
458 const EffArray_t efT0 = makeEffArray(efficiencyArrayInput,0,Transition,neBTE);
459 const EffArray_t efT1 = makeEffArray(efficiencyArrayInput,1,Transition,neBTE);
460 const EffArray_t efT2 = makeEffArray(efficiencyArrayInput,2,Transition,neBTE);
461 const EffArray_t efT3 = makeEffArray(efficiencyArrayInput,3,Transition,neBTE);
462 const EffArray_t efT4 = makeEffArray(efficiencyArrayInput,4,Transition,neBTE);
464 neBTE = coerceToOne(m_eventStat.m_eventsBTE[Endcap]);
465 const EffArray_t efE0 = makeEffArray(efficiencyArrayInput,0,Endcap,neBTE);
466 const EffArray_t efE1 = makeEffArray(efficiencyArrayInput,1,Endcap,neBTE);
467 const EffArray_t efE2 = makeEffArray(efficiencyArrayInput,2,Endcap,neBTE);
468 const EffArray_t efE3 = makeEffArray(efficiencyArrayInput,3,Endcap,neBTE);
469 const EffArray_t efE4 = makeEffArray(efficiencyArrayInput,4,Endcap,neBTE);
471 neBTE = coerceToOne(m_eventStat.m_eventsBTE[Forward]);
472 const EffArray_t efD0 = makeEffArray(efficiencyArrayInput,0,Forward,neBTE);
473 const EffArray_t efD1 = makeEffArray(efficiencyArrayInput,1,Forward,neBTE);
474 const EffArray_t efD2 = makeEffArray(efficiencyArrayInput,2,Forward,neBTE);
475 const EffArray_t efD3 = makeEffArray(efficiencyArrayInput,3,Forward,neBTE);
476 const EffArray_t efD4 = makeEffArray(efficiencyArrayInput,4,Forward,neBTE);
479 double efrec = ef0[0]+ef0[1]+ef0[2]+ef1[0]+ef1[1]+ef2[0];
480 double efrecB = efB0[0]+efB0[1]+efB0[2]+efB1[0]+efB1[1]+efB2[0];
481 double efrecT = efT0[0]+efT0[1]+efT0[2]+efT1[0]+efT1[1]+efT2[0];
482 double efrecE = efE0[0]+efE0[1]+efE0[2]+efE1[0]+efE1[1]+efE2[0];
483 double efrecD = efD0[0]+efD0[1]+efD0[2]+efD1[0]+efD1[1]+efD2[0];
485 ne = coerceToOne(m_eventStat.m_eventsPOS);
486 double efP[6];
for(
int i=0;
i!=6; ++
i) efP[
i] =
double(m_trackCollectionStat[
nc].m_efficiencyPOS[
i])/ne;
487 ne = coerceToOne(m_eventStat.m_eventsNEG);
488 double efN[6];
for(
int i=0;
i!=6; ++
i) efN[
i] =
double(m_trackCollectionStat[
nc].m_efficiencyNEG[
i])/ne;
490 out<<
"|-----------------------------------------------------------------------------------|\n";
491 out<<
"| Probability to lose 0 1 2 3 4 >=5 clusters |\n";
492 out<<
"|-----------------------------------------------------------------------------------|\n";
494 auto formattedOutput=[&
out](
auto & effArray){
495 for (
size_t i{};
i!=6;++
i){
496 out<<std::setw(9)<<std::setprecision(4)<<effArray[
i];
501 out<<
"| For all particles ";
503 out<<
"| For + particles ";
504 formattedOutput(efP);
505 out<<
"| For - particles ";
506 formattedOutput(efN);
507 out<<
"|-----------------------------------------------------------------------------------|\n";
508 out<<
"| Barrel region |\n";
509 out<<
"| 0 wrong clusters ";
510 formattedOutput(efB0);
511 out<<
"| 1 wrong clusters ";
512 formattedOutput(efB1);
513 out<<
"| 2 wrong clusters ";
514 formattedOutput(efB2);
515 out<<
"| 3 wrong clusters ";
516 formattedOutput(efB3);
517 out<<
"| >=4 wrong clusters ";
518 formattedOutput(efB4);
519 out<<
"|-----------------------------------------------------------------------------------|\n";
520 out<<
"| Transition region |\n";
521 out<<
"| 0 wrong clusters ";
522 formattedOutput(efT0);
523 out<<
"| 1 wrong clusters ";
524 formattedOutput(efT1);
525 out<<
"| 2 wrong clusters ";
526 formattedOutput(efT2);
527 out<<
"| 3 wrong clusters ";
528 formattedOutput(efT3);
529 out<<
"| >=4 wrong clusters ";
530 formattedOutput(efT4);
531 out<<
"|-----------------------------------------------------------------------------------|\n";
532 out<<
"| Endcap region |\n";
533 out<<
"| 0 wrong clusters ";
534 formattedOutput(efE0);
535 out<<
"| 1 wrong clusters ";
536 formattedOutput(efE1);
537 out<<
"| 2 wrong clusters ";
538 formattedOutput(efE2);
539 out<<
"| 3 wrong clusters ";
540 formattedOutput(efE3);
541 out<<
"| >=4 wrong clusters ";
542 formattedOutput(efE4);
543 out<<
"|-----------------------------------------------------------------------------------|\n";
544 out<<
"| Forward region |\n";
545 out<<
"| 0 wrong clusters ";
546 formattedOutput(efD0);
547 out<<
"| 1 wrong clusters ";
548 formattedOutput(efD1);
549 out<<
"| 2 wrong clusters ";
550 formattedOutput(efD2);
551 out<<
"| 3 wrong clusters ";
552 formattedOutput(efD3);
553 out<<
"| >=4 wrong clusters ";
554 formattedOutput(efD4);
556 out<<
"|-----------------------------------------------------------------------------------|\n";
557 pa = coerceToOne(m_eventStat.m_particleClusters[0]);
558 out<<
"| Efficiency reconstruction (number lose+wrong < 3) = "
559 <<std::setw(9)<<std::setprecision(5)<<efrec
561 <<std::setw(9)<<std::setprecision(5)<<efrec*
double(m_eventStat.m_events)/pa
564 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Barrel]);
565 out<<
"| For barrel region = "
566 <<std::setw(9)<<std::setprecision(5)<<efrecB
568 <<std::setw(9)<<std::setprecision(5)<<efrecB*
double(m_eventStat.m_eventsBTE[Barrel])/pa
571 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Transition]);
572 out<<
"| For transition region = "
573 <<std::setw(9)<<std::setprecision(5)<<efrecT
575 <<std::setw(9)<<std::setprecision(5)<<efrecT*
double(m_eventStat.m_eventsBTE[Transition])/pa
578 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Endcap]);
579 out<<
"| For endcap region = "
580 <<std::setw(9)<<std::setprecision(5)<<efrecE
582 <<std::setw(9)<<std::setprecision(5)<<efrecE*
double(m_eventStat.m_eventsBTE[Endcap])/pa
585 pa = coerceToOne(m_eventStat.m_particleClustersBTE[0][Forward]);
586 out<<
"| For forward region = "
587 <<std::setw(9)<<std::setprecision(5)<<efrecD
589 <<std::setw(9)<<std::setprecision(5)<<efrecD*
double(m_eventStat.m_eventsBTE[Forward])/pa
593 out<<
"|-----------------------------------------------------------------------------------|\n";
594 out<<
"| Reconstructed tracks + - +/-ratio error |\n";
595 out<<
"|-----------------------------------------------------------------------------------|\n";
597 pa = coerceToOne(m_trackCollectionStat[
nc].m_ntracksNEGB);
602 <<std::setw(10)<<m_trackCollectionStat[
nc].m_ntracksPOSB
603 <<std::setw(11)<<m_trackCollectionStat[
nc].m_ntracksNEGB
604 <<std::setw(11)<<std::setprecision(5)<<
ratio
605 <<std::setw(11)<<std::setprecision(5)<<
eratio<<
" |\n";
606 pa = coerceToOne(m_trackCollectionStat[
nc].m_ntracksNEGE);
611 <<std::setw(10)<<m_trackCollectionStat[
nc].m_ntracksPOSE
612 <<std::setw(11)<<m_trackCollectionStat[
nc].m_ntracksNEGE
613 <<std::setw(11)<<std::setprecision(5)<<
ratio
614 <<std::setw(11)<<std::setprecision(5)<<
eratio<<
" |\n";
615 pa = coerceToOne(m_trackCollectionStat[
nc].m_ntracksNEGFWD);
616 ratio =
double(m_trackCollectionStat[
nc].m_ntracksPOSFWD)/pa;
620 <<std::setw(10)<<m_trackCollectionStat[
nc].m_ntracksPOSFWD
621 <<std::setw(11)<<m_trackCollectionStat[
nc].m_ntracksNEGFWD
622 <<std::setw(11)<<std::setprecision(5)<<
ratio
623 <<std::setw(11)<<std::setprecision(5)<<
eratio<<
" |\n";
630 for(
int k = 0;
k!=50; ++
k) {
631 nt+=m_trackCollectionStat[
nc].m_total[
k];
632 ft+=m_trackCollectionStat[
nc].m_fake [
k];
633 if(!kf &&
nt) kf =
k;
638 out<<
"|-----------------------------------------------------------------------------------|\n";
639 out<<
"| Fake tracks rate for different number of clusters on track |\n";
640 out<<
"|-----------------------------------------------------------------------------------|\n";
642 for(
int k = kf;
k!=kf+6; ++
k) {
643 out<<
"| >= "<<std::setw(2)<<
k<<
" ";
647 for(
int k = kf;
k!=kf+6; ++
k) {
649 out<<
"|"<<std::setw(12)<<std::setprecision(5)<<
eff<<
" ";
650 nt-=m_trackCollectionStat[
nc].m_total[
k];
651 ft-=m_trackCollectionStat[
nc].m_fake [
k];
654 out<<
"|-----------------------------------------------------------------------------------|\n";
659 return StatusCode::SUCCESS;
672 out<<
"|--------------------------------------------------------------------------------------------------------------------|\n";
674 n = 65-
t->key().size();
675 std::string s1;
for(
int i=0;
i<
n; ++
i) s1.append(
" "); s1.append(
"|");
676 out<<
"| Location of input tracks | "<<
t->key()<<s1<<
"\n";
678 auto padString = [](
const std::string &
s){
679 const int n = 65 -
s.size();
680 return s + std::string(
n,
' ');
695 out<<
"| Pixel space points | "<<
s2<<
"|\n";
696 out<<
"| SCT space points | "<<
s3<<
"|\n";
697 out<<
"| Overlap space points | "<<
s4<<
"|\n";
698 out<<
"| Pixel clusters | "<<s5<<
"|\n";
699 out<<
"| SCT clusters | "<<s6<<
"|\n";
700 out<<
"| TRT clusters | "<<s9<<
"|\n";
701 out<<
"| Truth location for pixels | "<<s7<<
"|\n";
702 out<<
"| Truth location for sct | "<<s8<<
"|\n";
703 out<<
"| Truth location for trt | "<<s10<<
"|\n";
705 <<std::setw(14)<<std::setprecision(5)<<
m_ptcut
707 out<<
"| max pT cut | "
708 <<std::setw(14)<<std::setprecision(5)<<
m_ptcutmax
710 out<<
"| rapidity cut | "
711 <<std::setw(14)<<std::setprecision(5)<<
m_rapcut
713 out<<
"| min Radius | "
714 <<std::setw(14)<<std::setprecision(5)<<
m_rmin
716 out<<
"| max Radius | "
717 <<std::setw(14)<<std::setprecision(5)<<
m_rmax
719 out<<
"| Min. number clusters for generated track | "
720 <<std::setw(14)<<std::setprecision(5)<<
m_clcut
722 out<<
"| Min. number sp.points for generated track | "
723 <<std::setw(14)<<std::setprecision(5)<<
m_spcut
725 out<<
"| Min. number TRT clusters for generated track | "
726 <<std::setw(14)<<std::setprecision(5)<<
m_clcutTRT
728 out<<
"|--------------------------------------------------------------------------------------------------------------------|\n";
739 auto formatOutput = [&
out](
const auto val){
743 out<<
"|---------------------------------------------------------------------|\n";
744 out<<
"| m_nspacepoints | ";
746 out<<
"| m_nclusters | ";
748 out<<
"| Kine-Clusters size | ";
750 out<<
"| Kine-TRTClusters size | ";
752 out<<
"| Kine-SpacePoints size | ";
754 out<<
"| Number good kine tracks | ";
756 out<<
"|---------------------------------------------------------------------|\n";
771 std::unique_ptr<SG::ReadHandle<SiClusterContainer> > pixelcontainer;
772 std::unique_ptr<SG::ReadHandle<SiClusterContainer> > sctcontainer;
773 std::unique_ptr<SG::ReadHandle<TRT_DriftCircleContainer> > trtcontainer;
776 pixelcontainer = std::make_unique<SG::ReadHandle<SiClusterContainer> >(
m_clustersPixelname,ctx);
783 sctcontainer = std::make_unique<SG::ReadHandle<SiClusterContainer> >(
m_clustersSCTname,ctx);
790 trtcontainer = std::make_unique<SG::ReadHandle<TRT_DriftCircleContainer> >(
m_clustersTRTname,ctx);
799 if(pixelcontainer && pixelcontainer->isValid()) {
800 InDet::SiClusterContainer::const_iterator
w = (*pixelcontainer)->begin();
801 InDet::SiClusterContainer::const_iterator we = (*pixelcontainer)->end ();
805 InDet::SiClusterCollection::const_iterator
c = (*w)->begin();
806 InDet::SiClusterCollection::const_iterator ce = (*w)->end ();
811 ++m_eventStat.m_nclustersPTOT;
812 if(
isTruth(event_data,(*
c))) ++m_eventStat.m_nclustersPTOTt;
815 int nk =
kine(event_data,(*
c),Kine,999);
816 for(
int i=0;
i!=nk; ++
i) {
829 if(sctcontainer && sctcontainer->isValid()) {
830 InDet::SiClusterContainer::const_iterator
w = (*sctcontainer)->begin();
831 InDet::SiClusterContainer::const_iterator we = (*sctcontainer)->end ();
835 InDet::SiClusterCollection::const_iterator
c = (*w)->begin();
836 InDet::SiClusterCollection::const_iterator ce = (*w)->end ();
841 ++m_eventStat.m_nclustersSTOT;
842 if(
isTruth(event_data,(*
c))) ++m_eventStat.m_nclustersSTOTt;
844 int nk =
kine(event_data,(*
c),Kine,999);
845 for(
int i=0;
i!=nk; ++
i) {
853 if(trtcontainer && trtcontainer->isValid()) {
856 InDet::TRT_DriftCircleContainer::const_iterator
w = (*trtcontainer)->begin();
857 InDet::TRT_DriftCircleContainer::const_iterator we = (*trtcontainer)->end ();
861 InDet::TRT_DriftCircleCollection::const_iterator
c = (*w)->begin();
862 InDet::TRT_DriftCircleCollection::const_iterator ce = (*w)->end ();
867 int nk =
kine(event_data,(*
c),Kine,999);
892 for(; spc != spce; ++spc) {
896 for(; sp != spe; ++sp) {
899 int nk =
kine(event_data,(*sp)->clusterList().first,Kine,999);
900 for(
int i=0;
i!=nk; ++
i) {
921 for(; spc != spce; ++spc) {
926 for(; sp != spe; ++sp) {
930 int nk =
kine(event_data,(*sp)->clusterList().first,(*sp)->clusterList().second,Kine,999);
931 for(
int i=0;
i!=nk; ++
i) {
951 for (; sp!=spe; ++sp) {
954 int nk =
kine(event_data,(*sp)->clusterList().first,(*sp)->clusterList().second,Kine,999);
955 for(
int i=0;
i!=nk; ++
i) {
986 std::list<int> worskine;
992 unsigned int nc = 1 ;
994 auto coerceTo49 = [] (
const size_t idx){
995 return (
idx<50) ?
idx : 49;
1000 if((*c).first==
k0) {++
nc;
continue;}
1003 const size_t clusterIdx =coerceTo49(
nc);
1008 const size_t spacepointIdx =coerceTo49(
ns);
1036 for(
auto & pThisCluster: worskine) {
1047 if (not de)
continue;
1051 if(de->isBarrel()) {
1060 if(de->isBarrel()) {
1089 if(++
nc >= 100)
return;
1101 int KINE[200],NKINE[200];
1106 se = (*t)->trackStateOnSurfaces()->end();
1115 const AmgVector(5)& Vpf = tpf ->parameters ();
1116 double pTf = std::abs(
std::sin(Vpf[3])/Vpf[4]);
1160 int Kine[1000], nk=
kine0(event_data,rd,Kine,999); ++
NC;
if(!nk) ++N0;
1162 for(
int k = 0;
k!=nk; ++
k) {
1165 for(;
n!=NK; ++
n) {
if(Kine[
k]==KINE[
n]) {++NKINE[
n];
break;}}
1166 if(
n==NK) {KINE[NK] = Kine[
k]; NKINE[NK] = 1;
if (NK < 200) ++NK;}
1168 for(
int n=0;
n!=NK; ++
n) {
if(NKINE[
n]>nkm) nkm = NKINE[
n];}
1171 for(
int n=0;
n!=NK; ++
n) {
1173 int NQ = 1000*NKINE[
n]+(
NC-NKINE[
n]);
1175 event_data.
m_tracks[
nc].insert(std::make_pair(KINE[
n],NQ));
1199 std::multimap<int,int>::const_iterator
t, te = event_data.
m_tracks[
nc].end();
1203 int k = (*p).barcode();
1209 if((*t).first!=
k)
break;
1210 int ts = (*t).second/1000;
1211 int ws = (*t).second%1000;
1212 if (
ts >
m ) {
m =
ts;
w = ws;}
1213 else if(
ts==
m &&
w > ws) {
w = ws;}
1215 int d =
n-
m;
if(
d<0)
d = 0;
else if(
d > 5)
d=5;
if(
w>4)
w = 4;
1220 int ch = (*p).charge();
1243 int Kine1[1000],Kine2[1000];
1245 int n2 =
kine(event_data,
d2,Kine2,
nmax);
if(!n2)
return nkine;
1247 for(
int i = 0;
i!=
n1; ++
i) {
1248 for(
int j = 0; j!=n2; ++j) {
1249 if(Kine1[
i]==Kine2[j]) {Kine[nkine++] = Kine1[
i];
break;}
1263 PRD_MultiTruthCollection::const_iterator mce;
1264 PRD_MultiTruthCollection::const_iterator
mc =
findTruth(event_data,
d,mce);
1269 for(;
mc!=mce; ++
mc) {
1271 if( (*mc).first !=
ID )
return nkine;
1273 int k = (*mc).second.barcode();
if(
k<=0)
continue;
1276 if(!pa or !pa->production_vertex())
continue;
1278 int pdg = std::abs(pa->pdg_id());
if(
m_pdg &&
m_pdg != pdg )
continue;
1281 if(!pd or std::abs(pd->charge()) < .5)
continue;
1285 double px = pa->momentum().px();
1286 double py = pa->momentum().py();
1287 double pz = pa->momentum().pz();
1289 if( pt < m_ptcut || pt >
m_ptcutmax)
continue;
1293 double t = std::abs(
pz)/
pt;
1298 double vx = pa->production_vertex()->position().x();
1299 double vy = pa->production_vertex()->position().y();
1300 double r = std::sqrt(vx*vx+vy*vy);
1301 if( r < m_rmin || r >
m_rmax)
continue;
1303 Kine[nkine] =
k;
if(++nkine >=
nmax)
break;
1316 PRD_MultiTruthCollection::const_iterator mce;
1317 PRD_MultiTruthCollection::const_iterator
mc =
findTruth(event_data,
d,mce);
1322 for(;
mc!=mce; ++
mc) {
1324 if( (*mc).first !=
ID )
return nkine;
1326 int k = (*mc).second.barcode();
if(
k<=0)
continue;
1327 Kine[nkine] =
k;
if(++nkine >=
nmax)
break;
1339 PRD_MultiTruthCollection::const_iterator mce;
1340 PRD_MultiTruthCollection::const_iterator
mc =
findTruth(event_data,
d,mce);
1351 std::multimap<int,const Trk::PrepRawData*>::const_iterator
k = event_data.
m_kinecluster.find(K);
1354 if((*k).first!= K)
return false;
1355 if(
d->detectorElement()==(*k).second->detectorElement())
return true;
1370 std::multimap<int,const Trk::SpacePoint*>::const_iterator
k = event_data.
m_kinespacepoint.find(K);
1375 if((*k).first!= K)
return false;
1388 if((*k).first!= K)
return false;
1414 std::list<int>::const_iterator dif = event_data.
m_difference[
nc].begin();
1416 std::multimap<int,const Trk::PrepRawData*>::const_iterator
c,ce = event_data.
m_kinecluster.end();
1419 std::string s1;
for(
int i=0;
i<
n; ++
i) s1.append(
" "); s1.append(
"|");
1420 std::stringstream
out;
1421 out<<
"|----------------------------------------------------------------------------------------|\n";
1423 out<<
"|----------------------------------------------------------------------------------------|\n";
1424 out<<
"| # pdg kine Ncl Ntr Nsp Lose pT(MeV) rapidity radius z |\n";
1425 out<<
"|----------------------------------------------------------------------------------------|\n";
1430 int k = (*p).barcode();
1435 PRD_MultiTruthCollection::const_iterator mce;
1436 PRD_MultiTruthCollection::const_iterator
mc =
findTruth(event_data,
d,mce);
1440 for(;
mc!=mce; ++
mc) {
1441 if((*mc).first !=
ID)
break;
1442 if((*mc).second.barcode()==
k) {Q=
true;
break;}
1449 double px = pa->momentum().px();
1450 double py = pa->momentum().py();
1451 double pz = pa->momentum().pz();
1452 double vx = pa->production_vertex()->position().x();
1453 double vy = pa->production_vertex()->position().y();
1454 double vz = pa->production_vertex()->position().z();
1456 double t = std::atan2(
pt,
pz);
1458 double r = std::sqrt(vx*vx+vy*vy);
1462 <<std::setw(6)<<pa->pdg_id()
1467 <<std::setw(4)<<(*dif)
1468 <<std::setw(12)<<std::setprecision(5)<<
pt
1469 <<std::setw(12)<<std::setprecision(5)<<ra
1470 <<std::setw(12)<<std::setprecision(5)<<
r
1471 <<std::setw(12)<<std::setprecision(5)<<vz
1476 out<<
"|----------------------------------------------------------------------------------------|\n";
1486 PRD_MultiTruthCollection::const_iterator
1490 PRD_MultiTruthCollection::const_iterator& mce)
1496 PRD_MultiTruthCollection::const_iterator
mc;
1503 for (
int i=0;
i<3;
i++) {
1505 mce=truth[
i]->end();
1506 return truth[
i]->end();
1509 throw std::runtime_error(
"Neither Pixel, SCT nor TRT truth.");
1522 PRD_MultiTruthCollection::const_iterator mce;
1523 PRD_MultiTruthCollection::const_iterator
mc =
findTruth(event_data,
d,mce);
1525 for(;
mc!=mce; ++
mc) {
1526 if((*mc).second.barcode()==
k) {
1531 double px =
pat->momentum().px();
1532 double py =
pat->momentum().py();
1533 double pz =
pat->momentum().pz();
1535 double t = std::atan2(
pt,
pz) ;
1542 ra > 1.6 ? rap = 2 : ra > .8 ? rap = 1 : rap = 0;
1544 int pdg =
pat->pdg_id();
1547 double ch = pd->charge();
if(pdg < 0)
ch = -
ch;
1548 if(
ch > .5)
return 1;
1549 if(
ch < -.5)
return -1;