Actually write out the variables into the map.
< number of clusters with Elab > m_EminLab outside of high background endcap region
< number of clusters with Elab > m_EminLab4Cluster outside of high background endcap region
< number of clusters with Elab > m_EminLab3Cluster outside of high background endcap region
< number of clusters with Ecms > m_Ehigh outside of high background endcap region
85{
86 calculationResult["nTrkLoose"] = 0;
87 calculationResult["nTrkTight"] = 0;
88 calculationResult["ee2leg"] = 0;
89 calculationResult["nEmedium"] = 0;
90 calculationResult["nElow"] = 0;
91 calculationResult["nEhigh"] = 0;
92 calculationResult["nE180Lab"] = 0;
93 calculationResult["nE300Lab"] = 0;
94 calculationResult["nE500Lab"] = 0;
95 calculationResult["nE2000CMS"] = 0;
96 calculationResult["nE4000CMS"] = 0;
97 calculationResult["nE250Lab"] = 0;
98 calculationResult["nMaxEPhotonAcc"] = 0;
99 calculationResult["dphiCmsClust"] = NAN;
101 calculationResult["netChargeLoose"] = 0;
102 calculationResult["maximumPCMS"] = NAN;
103 calculationResult["maximumPLab"] = NAN;
104 calculationResult["eexx"] = 0;
105 calculationResult["ee1leg1trk"] = 0;
106 calculationResult["nEhighLowAng"] = 0;
107 calculationResult["nEsingleClust"] = 0;
108 calculationResult["nEsinglePhotonBarrel"] = 0;
109 calculationResult["nEsinglePhotonExtendedBarrel"] = 0;
110 calculationResult["nEsinglePhotonEndcap"] = 0;
111 calculationResult["nEsingleElectronBarrel"] = 0;
112 calculationResult["nEsingleElectronExtendedBarrel"] = 0;
113 calculationResult["nReducedEsinglePhotonReducedBarrel"] = 0;
114 calculationResult["nVetoClust"] = 0;
115 calculationResult["chrgClust2GeV"] = 0;
116 calculationResult["neutClust045GeVAcc"] = 0;
117 calculationResult["neutClust045GeVBarrel"] = 0;
118 calculationResult["singleTagLowMass"] = 0;
119 calculationResult["singleTagHighMass"] = 0;
120 calculationResult["n2GeVNeutBarrel"] = 0;
121 calculationResult["n2GeVNeutEndcap"] = 0;
122 calculationResult["n2GeVChrg"] = 0;
123 calculationResult["n2GeVPhotonBarrel"] = 0;
124 calculationResult["n2GeVPhotonEndcap"] = 0;
125 calculationResult["ee1leg"] = 0;
126 calculationResult["ee1leg1clst"] = 0;
127 calculationResult["ee1leg1e"] = 0;
128 calculationResult["ee2clst"] = 0;
129 calculationResult["gg2clst"] = 0;
130 calculationResult["eeee"] = 0;
131 calculationResult["eemm"] = 0;
132 calculationResult["eexxSelect"] = 0;
133 calculationResult["radBhabha"] = 0;
134 calculationResult["eeBrem"] = 0;
135 calculationResult["isrRadBhabha"] = 0;
136 calculationResult["muonPairECL"] = 0;
137 calculationResult["selectee1leg1trk"] = 0;
138 calculationResult["selectee1leg1clst"] = 0;
139 calculationResult["selectee"] = 0;
140 calculationResult["ggBarrelVL"] = 0;
141 calculationResult["ggBarrelLoose"] = 0;
142 calculationResult["ggBarrelTight"] = 0;
143 calculationResult["ggEndcapVL"] = 0;
144 calculationResult["ggEndcapLoose"] = 0;
145 calculationResult["ggEndcapTight"] = 0;
146 calculationResult["muonPairV"] = 0;
147 calculationResult["selectmumu"] = 0;
148 calculationResult["singleMuon"] = 0;
149 calculationResult["cosmic"] = 0;
150 calculationResult["displacedVertex"] = 0;
151 calculationResult["eeFlat0"] = 0;
152 calculationResult["eeFlat1"] = 0;
153 calculationResult["eeFlat2"] = 0;
154 calculationResult["eeFlat3"] = 0;
155 calculationResult["eeFlat4"] = 0;
156 calculationResult["eeFlat5"] = 0;
157 calculationResult["eeFlat6"] = 0;
158 calculationResult["eeFlat7"] = 0;
159 calculationResult["eeFlat8"] = 0;
160 calculationResult["eeOneClust"] = 0;
161
162
163 calculationResult["nTrkLooseB"] = 0;
164 calculationResult["nTrkTightB"] = 0;
165 calculationResult["maximumPCMSB"] = NAN;
166 calculationResult["netChargeLooseB"] = 0;
167 calculationResult["muonPairVB"] = 0;
168 calculationResult["eeBremB"] = 0;
169 calculationResult["singleTagLowMassB"] = 0;
170 calculationResult["singleTagHighMassB"] = 0;
171 calculationResult["radBhabhaB"] = 0;
174
179 bool bha3d;
180 try {
182 } catch (const std::exception&) {
183 bha3d = false;
184 }
185 bool bhapurPsnm;
186 try {
188 } catch (const std::exception&) {
189 bhapurPsnm = false;
190 }
191 bool bhapurFtdl;
192 try {
194 } catch (const std::exception&) {
195 bhapurFtdl = false;
196 }
197 bool lml1;
198 try {
200 } catch (const std::exception&) {
201 lml1 = false;
202 }
203 bool l1_bit_f;
204 try {
206 } catch (const std::exception&) {
207 l1_bit_f = false;
208 }
209 calculationResult["bha3d"] = bha3d;
210 calculationResult["bhapur"] = bhapurPsnm;
211 calculationResult["bhapur_lml1"] = lml1 and bhapurFtdl;
212 calculationResult["l1_bit_f"] = l1_bit_f;
213 } else {
214 calculationResult["l1_trigger_random"] = 1;
215 calculationResult["l1_trigger_delayed_bhabha"] = 0;
216 calculationResult["l1_trigger_poisson"] = 0;
217 calculationResult["bha3d"] = 0;
218 calculationResult["bhapur"] = 0;
219 calculationResult["bhapur_lml1"] = 0;
220 calculationResult["l1_bit_f"] = 0;
221 }
222
223
224 calculationResult["l1_trg_NN_info"] = 0;
225 if (
m_bitsNN.isValid() and
m_bitsNN.getEntries() > 0) {calculationResult[
"l1_trg_NN_info"] = 1;}
226
227 calculationResult["true"] = 1;
228 calculationResult["false"] = 0;
229
230
231
233 const ROOT::Math::XYZVector clustervertex = cUtil.
GetIPPosition();
235 ROOT::Math::PxPyPzEVector p4ofCOM;
236 p4ofCOM.SetPxPyPzE(0, 0, 0, boostrotate.
getCMSEnergy());
237
238
239 std::map<short, std::optional<MaximumPtTrack>> maximumPtTracks = {
240 { -1, {}}, {1, {}}
241 };
242
243
244 std::map<short, std::optional<MaximumPtTrack>> maximumPtTracksWithoutZCut = {
245 { -1, {}}, {1, {}}
246 };
247
248
251 if (not trackFitResult) {
252
253 continue;
254 }
255
256
257 if (trackFitResult->getHitPatternCDC().getNHits() == 0) {
258 continue;
259 }
260
261
262 const double z0 = trackFitResult->getZ0();
264 calculationResult["nTrkTight"] += 1;
265 }
266
267
268 const short charge = trackFitResult->getChargeSign();
269 if (charge == 0) {
270 continue;
271 }
272
273 const ROOT::Math::PxPyPzEVector& momentumLab = trackFitResult->get4Momentum();
274 const ROOT::Math::PxPyPzEVector momentumCMS = boostrotate.
rotateLabToCms() * momentumLab;
275 double pCMS = momentumCMS.P();
276 double pLab = momentumLab.P();
277
278
279 const double pT = trackFitResult->getTransverseMomentum();
280 const auto& currentMaximum = maximumPtTracksWithoutZCut.at(charge);
281 if (not currentMaximum or pT > currentMaximum->pT) {
284 newMaximum.
track = &track;
285 newMaximum.
pCMS = pCMS;
286 newMaximum.
pLab = pLab;
287 newMaximum.
p4CMS = momentumCMS;
288 newMaximum.
p4Lab = momentumLab;
289 maximumPtTracksWithoutZCut[
charge] = newMaximum;
290 }
291
292
294 calculationResult["nTrkLoose"] += 1;
295 calculationResult[
"netChargeLoose"] +=
charge;
296
297 if (std::isnan(calculationResult["maximumPCMS"]) or pCMS > calculationResult["maximumPCMS"]) {
298 calculationResult["maximumPCMS"] = pCMS;
299 }
300
301 if (std::isnan(calculationResult["maximumPLab"]) or pLab > calculationResult["maximumPLab"]) {
302 calculationResult["maximumPLab"] = pLab;
303 }
304
305
306 const double pTLoose = trackFitResult->getTransverseMomentum();
307 const auto& currentMaximumLoose = maximumPtTracks.at(charge);
308 if (not currentMaximumLoose or pTLoose > currentMaximumLoose->pT) {
310 newMaximum.
pT = pTLoose;
311 newMaximum.
track = &track;
312 newMaximum.
pCMS = pCMS;
313 newMaximum.
pLab = momentumLab.P();
314 newMaximum.
p4CMS = momentumCMS;
315 newMaximum.
p4Lab = momentumLab;
316 maximumPtTracks[
charge] = newMaximum;
317 }
318 }
319
320
321
322 if (trackFitResult->getHitPatternCDC().getNHits() >= 5) {
323 if (std::abs(z0) < 1.) {calculationResult["nTrkTightB"] += 1;}
324 if (std::abs(z0) < 5.) {
325 calculationResult["nTrkLooseB"] += 1;
326 calculationResult[
"netChargeLooseB"] +=
charge;
327 if (std::isnan(calculationResult["maximumPCMSB"]) or pCMS > calculationResult["maximumPCMSB"]) {
328 calculationResult["maximumPCMSB"] = pCMS;
329
330 }
331 }
332 }
333
334 }
335
336
337
338 std::vector<SelectedECLCluster> selectedClusters;
340
341 const double time = cluster.getTime();
343 std::abs(time) > 200 or
345 continue;
346 }
347
348 const double dt99 = cluster.getDeltaTime99();
349
350
355 selectedCluster.
cluster = &cluster;
359 selectedCluster.
isTrack = cluster.isTrack();
360
361 selectedClusters.push_back(selectedCluster);
362
364 calculationResult["nElow"] += 1;
365 }
367 calculationResult["nEmedium"] += 1;
368 }
370 calculationResult["nEhigh"] += 1;
371 }
372
374 calculationResult["nVetoClust"] += 1;
375 }
376
377
378
379 const double thetaLab = selectedCluster.
p4Lab.Theta() * TMath::RadToDeg();
380 const double zmva = cluster.getZernikeMVA();
381 const bool photon = zmva > 0.5 and not selectedCluster.
isTrack;
382 const bool electron = zmva > 0.5 and selectedCluster.
isTrack;
383
384
386 calculationResult["chrgClust2GeV"] += 1;
387 }
389 const bool isInAcceptance = 17. < thetaLab and thetaLab < 150.;
390 if (isInAcceptance) {calculationResult["neutClust045GeVAcc"] += 1;}
391 const bool isInBarrel = 30. < thetaLab and thetaLab < 130.;
392 if (isInBarrel) {calculationResult["neutClust045GeVBarrel"] += 1;}
393 }
394
395
396
397 const bool notInHighBackgroundEndcapRegion = 18.5 < thetaLab and thetaLab < 139.3;
399 calculationResult["nE180Lab"] += 1;
400 }
401
403 calculationResult["nE300Lab"] += 1;
404 }
405
407 calculationResult["nE500Lab"] += 1;
408 }
409
410 if (selectedCluster.
energyCMS >
m_Ehigh and notInHighBackgroundEndcapRegion) {
411 calculationResult["nE2000CMS"] += 1;
412 }
413
414
416 calculationResult["nE250Lab"] += 1;
417 }
419 calculationResult["nE4000CMS"] += 1;
420 }
421
422
424 calculationResult["nEsingleClust"] += 1;
425
426 const bool barrelRegion = thetaLab > 45 and thetaLab < 115;
427 const bool extendedBarrelRegion = thetaLab > 30 and thetaLab < 130;
428 const bool endcapRegion = (thetaLab > 22 and thetaLab < 45) or (thetaLab > 115 and thetaLab < 145);
429
430 if (photon and barrelRegion) {
431 calculationResult["nEsinglePhotonBarrel"] += 1;
432 }
433
434 if (photon and extendedBarrelRegion) {
435 calculationResult["nEsinglePhotonExtendedBarrel"] += 1;
436 }
437
438 if (electron and barrelRegion) {
439 calculationResult["nEsingleElectronBarrel"] += 1;
440 }
441
442 if (electron and extendedBarrelRegion) {
443 calculationResult["nEsingleElectronExtendedBarrel"] += 1;
444 }
445
446 if (photon and endcapRegion) {
447 calculationResult["nEsinglePhotonEndcap"] += 1;
448 }
449 }
450
452 const bool reducedBarrelRegion = thetaLab > 44 and thetaLab < 98;
453
454 if (photon and reducedBarrelRegion) {
455 calculationResult["nReducedEsinglePhotonReducedBarrel"] += 1;
456 }
457 }
458
460
461 const bool barrelRegion = thetaLab > 32 and thetaLab < 130;
462 const bool endcapRegion = (thetaLab > 22 and thetaLab < 32) or (thetaLab > 130 and thetaLab < 145);
463 const bool lowAngleRegion = thetaLab < 22 or thetaLab > 145;
464
465 if (not selectedCluster.
isTrack and barrelRegion) {
466 calculationResult["n2GeVNeutBarrel"] += 1;
467 }
468 if (not selectedCluster.
isTrack and endcapRegion) {
469 calculationResult["n2GeVNeutEndcap"] += 1;
470 }
471 if (selectedCluster.
isTrack and not lowAngleRegion) {
472 calculationResult["n2GeVChrg"] += 1;
473 }
474 if (lowAngleRegion) {
475 calculationResult["nEhighLowAng"] += 1;
476 }
477 if (photon and barrelRegion) {
478 calculationResult["n2GeVPhotonBarrel"] += 1;
479 }
480 if (photon and endcapRegion) {
481 calculationResult["n2GeVPhotonEndcap"] += 1;
482 }
483 }
484 }
485 std::sort(selectedClusters.begin(), selectedClusters.end(), [](const auto & lhs, const auto & rhs) {
486 return lhs.energyCMS > rhs.energyCMS;
487 });
488
489
490 for (short charge : { -1, 1}) {
491 auto& maximumPtTrack = maximumPtTracks.at(charge);
492 if (not maximumPtTrack) {
493 continue;
494 }
495
496
497
498 for (
auto& cluster : maximumPtTrack->track->getRelationsTo<
ECLCluster>()) {
501 }
502 }
503 maximumPtTrack->clusterEnergySumCMS =
504 maximumPtTrack->clusterEnergySumLab * maximumPtTrack->pCMS / maximumPtTrack->pLab;
505
506
507 if (maximumPtTrack->clusterEnergySumCMS > 4.5) {
508 calculationResult["ee1leg"] = 1;
509 }
510
511
512 if (maximumPtTrack->pCMS > 3 and maximumPtTrack->clusterEnergySumLab > 0. and maximumPtTrack->clusterEnergySumLab < 1.) {
513 calculationResult["singleMuon"] = 1;
514 }
515 }
516
517
518 if (maximumPtTracks.at(-1) and maximumPtTracks.at(1)) {
521
522 double dphi = std::abs(negativeTrack.
p4CMS.Phi() - positiveTrack.
p4CMS.Phi()) * TMath::RadToDeg();
523 if (dphi > 180) {
524 dphi = 360 - dphi;
525 }
526
529 const double negativeP = negativeTrack.
pCMS;
532 const double positiveP = positiveTrack.
pCMS;
533
534
535 const double thetaSum = (negativeTrack.
p4CMS.Theta() + positiveTrack.
p4CMS.Theta()) * TMath::RadToDeg();
536 const double dthetaSum = std::abs(thetaSum - 180);
537 const auto back2back = dphi > 175 and dthetaSum < 15;
538 if (back2back and negativeClusterSum > 3 and positiveClusterSum > 3 and
539 (negativeClusterSum > 4.5 or positiveClusterSum > 4.5)) {
540 calculationResult["ee2leg"] = 1;
541 }
542
543
544 if (back2back and ((negativeClusterSum > 4.5 and positiveP > 3) or (positiveClusterSum > 4.5 and negativeP > 3))) {
545 calculationResult["ee1leg1trk"] = 1;
546 }
547
548
549
550 if ((negativeClusterSum > 4.5 and positiveClusterSum > 0.8 * positiveP) or
551 (positiveClusterSum > 4.5 and negativeClusterSum > 0.8 * negativeP)) {
552 calculationResult["ee1leg1e"] = 1;
553 }
554
555
556 const ROOT::Math::PxPyPzEVector p4Miss = p4ofCOM - negativeTrack.
p4CMS - positiveTrack.
p4CMS;
557 const double pmissTheta = p4Miss.Theta() * TMath::RadToDeg();
558 const double pmissp = p4Miss.P();
561
562 const bool electronEP = positiveClusterSum > 0.8 * positiveP or negativeClusterSum > 0.8 * negativeP;
563 const bool notMuonPair = negativeClusterSumLab > 1 or positiveClusterSumLab > 1;
564 const double highp = std::max(negativeP, positiveP);
565 const double lowp = std::min(negativeP, positiveP);
566 const bool lowEdep = negativeClusterSumLab < 0.5 and positiveClusterSumLab < 0.5;
567
568
569 if (calculationResult["maximumPCMS"] < 2 and dphi > 160 and (pmissTheta < 25. or pmissTheta > 155.)) {
570 calculationResult["eexxSelect"] = 1;
571 if (electronEP) {
572 calculationResult["eeee"] = 1;
573 } else {
574 calculationResult["eemm"] = 1;
575 }
576 }
577
578
579 if ((pmissTheta < 20. or pmissTheta > 160.) and
580 ((calculationResult["maximumPCMS"] < 1.2 and dphi > 150.) or
581 (calculationResult["maximumPCMS"] < 2. and 175. < dphi))) {
582 calculationResult["eexx"] = 1;
583 }
584
585
586 if (negativeP > 1. and pmissTheta > 10. and pmissTheta < 170. and positiveP > 1. and dphi < 170. and pmissp > 1. and electronEP) {
587 if (calculationResult["nTrkLoose"] == 2 and calculationResult["nTrkTight"] >= 1) {
588 calculationResult["radBhabha"] = 1;
589 }
590 if (calculationResult["nTrkLooseB"] == 2 and calculationResult["nTrkTightB"] >= 1) {
591 calculationResult["radBhabhaB"] = 1;
592 }
593 }
594
595
596 if (negativeP > 2. and positiveP > 2. and 2 == calculationResult["nTrkLoose"] and
597 calculationResult["nTrkTight"] >= 1 and dphi > 175. and
598 (pmissTheta < 5. or pmissTheta > 175.) and electronEP) {
599 calculationResult["isrRadBhabha"] = 1;
600 }
601
602
603 if (highp > 4.5 and notMuonPair and pmissp > 1. and (relMissAngle0 < 5. or relMissAngle1 < 5.)) {
604 if (calculationResult["nTrkLoose"] == 2) { calculationResult["eeBrem"] = 1;}
605 if (calculationResult["nTrkLooseB"] >= 1) { calculationResult["eeBremB"] = 1;}
606 }
607
608
609 if (highp > 4.5 and lowEdep and thetaSum > 175. and thetaSum < 185. and dphi > 175.) {
610 if (calculationResult["nTrkLoose"] == 2) {calculationResult["muonPairV"] = 1;}
611 if (calculationResult["nTrkLooseB"] == 2) {calculationResult["muonPairVB"] = 1;}
612 }
613
614
615 if (highp > 3. and lowp > 2.5 and dphi > 165. and
616 ((negativeClusterSumLab > 0. and negativeClusterSumLab < 1.) or
617 (positiveClusterSumLab > 0. and positiveClusterSumLab < 1.))) {
618 calculationResult["selectmumu"] = 1;
619 }
620 }
621
622
623 if (selectedClusters.size() >= 2) {
626
627 double dphi = std::abs(firstCluster.
p4CMS.Phi() - secondCluster.
p4CMS.Phi()) * TMath::RadToDeg();
628 if (dphi > 180) {
629 dphi = 360 - dphi;
630 }
631 double thetaSum = (firstCluster.
p4CMS.Theta() + secondCluster.
p4CMS.Theta()) * TMath::RadToDeg();
632 double dthetaSum = std::abs(thetaSum - 180);
633
634
635 calculationResult["dphiCmsClust"] = dphi;
636 for (int ic = 0; ic < 2; ic++) {
637 const double thetaLab = selectedClusters[ic].p4Lab.Theta() * TMath::RadToDeg();
638 const bool isInAcceptance = 17. < thetaLab and thetaLab < 150.;
639 const ECLCluster* cluster = selectedClusters[ic].cluster;
640 const double zmva = cluster->getZernikeMVA();
641 const bool photon = zmva > 0.5 and not selectedClusters[ic].isTrack;
642 if (isInAcceptance and photon) {calculationResult["nMaxEPhotonAcc"] += 1;}
643 }
644
645 const double firstEnergy = firstCluster.
p4CMS.E();
646 const double secondEnergy = secondCluster.
p4CMS.E();
647
648 const bool highEnergetic = firstEnergy > 3 and secondEnergy > 3 and (firstEnergy > 4.5 or secondEnergy > 4.5);
649
650 if (dphi > 160 and dphi < 178 and dthetaSum < 15 and highEnergetic) {
651 calculationResult["ee2clst"] = 1;
652 }
653
654 if (dphi > 178 and dthetaSum < 15 and highEnergetic) {
655 calculationResult["gg2clst"] = 1;
656 }
657
658 if ((calculationResult["ee2clst"] == 1 or calculationResult["gg2clst"] == 1) and
659 calculationResult["ee1leg"] == 1) {
660 calculationResult["ee1leg1clst"] = 1;
661 }
662
663 const double Elab0 = firstCluster.
p4Lab.E();
664 const double Elab1 = secondCluster.
p4Lab.E();
665
666
667 if (firstEnergy > 2 and secondEnergy > 2) {
668 const double thetaLab0 = firstCluster.
p4Lab.Theta() * TMath::RadToDeg();
669 const double thetaLab1 = secondCluster.
p4Lab.Theta() * TMath::RadToDeg();
670
671 const bool barrel0 = 32. < thetaLab0 and thetaLab0 < 130.;
672 const bool barrel1 = 32. < thetaLab1 and thetaLab1 < 130.;
673 const bool oneClustersAbove4 = firstEnergy > 4 or secondEnergy > 4;
674 const bool oneIsNeutral = not firstCluster.
isTrack or not secondCluster.
isTrack;
675 const bool bothAreNeutral = not firstCluster.
isTrack and not secondCluster.
isTrack;
676 const bool oneIsBarrel = barrel0 or barrel1;
677 const bool dphiCutExtraLoose = dphi > 175;
678 const bool dphiCutLoose = dphi > 177;
679 const bool dphiCutTight = dphi > 177.5;
680
681 if (dphiCutExtraLoose and oneIsNeutral and oneIsBarrel) {
682 calculationResult["ggBarrelVL"] = 1;
683 }
684 if (oneClustersAbove4 and dphiCutLoose and oneIsNeutral and oneIsBarrel) {
685 calculationResult["ggBarrelLoose"] = 1;
686 }
687 if (oneClustersAbove4 and dphiCutTight and bothAreNeutral and oneIsBarrel) {
688 calculationResult["ggBarrelTight"] = 1;
689 }
690 if (dphiCutExtraLoose and oneIsNeutral and not oneIsBarrel) {
691 calculationResult["ggEndcapVL"] = 1;
692 }
693 if (oneClustersAbove4 and dphiCutLoose and oneIsNeutral and not oneIsBarrel) {
694 calculationResult["ggEndcapLoose"] = 1;
695 }
696 if (oneClustersAbove4 and dphiCutTight and bothAreNeutral and not oneIsBarrel) {
697 calculationResult["ggEndcapTight"] = 1;
698 }
699 }
700
701 const double minEnergy = std::min(Elab0, Elab1);
702 const double maxEnergy = std::max(Elab0, Elab1);
703 if (dphi > 155 and thetaSum > 165 and thetaSum < 195 and minEnergy > 0.15 and minEnergy < 0.5 and
704 maxEnergy > 0.15 and maxEnergy < 0.5) {
705 calculationResult["muonPairECL"] = 1;
706 }
707
708 }
709
710
711
712 double thetaFlatten = 0;
713
714
715 for (short charge : { -1, 1}) {
716 const auto& maximumPtTrack = maximumPtTracks.at(charge);
717 if (not maximumPtTrack) {
718 continue;
719 }
720
721 if (maximumPtTrack->clusterEnergySumCMS > 1.5) {
722 double invMass = 0.;
723 double tempFlatten = 0.;
725 double tempInvMass = (maximumPtTrack->p4Lab + cluster.p4Lab).M();
726 if (tempInvMass > invMass) {
727 invMass = tempInvMass;
728 if (charge == 1) {
729 tempFlatten = cluster.p4Lab.Theta() * TMath::RadToDeg();
730 }
731 }
732 }
733 if (charge == -1) {
734 tempFlatten = maximumPtTrack->p4Lab.Theta() * TMath::RadToDeg();
735 }
736 if (invMass > 5.29) {
737 calculationResult["selectee1leg1clst"] = 1;
738 thetaFlatten = tempFlatten;
739 }
740 }
741 }
742
743
744 if (maximumPtTracks.at(-1) and maximumPtTracks.at(1)) {
747 const double invMass = (negativeTrack.
p4Lab + positiveTrack.
p4Lab).M();
749 calculationResult["selectee1leg1trk"] = 1;
750 }
751
752 thetaFlatten = negativeTrack.
p4Lab.Theta() * TMath::RadToDeg();
753
754
757 if ((invMass > 9.) and (missNegClust or missPosClust)) {
758 calculationResult["eeOneClust"] = 1;
759 }
760 }
761
762 if (calculationResult["selectee1leg1trk"] == 1 or calculationResult["selectee1leg1clst"] == 1) {
763 for (int iflat = 0; iflat < 9; iflat++) {
764 const std::string& eeFlatName = "eeFlat" + std::to_string(iflat);
765 calculationResult[eeFlatName] =
766 thetaFlatten > flatBoundaries[iflat] and thetaFlatten < flatBoundaries[iflat + 1];
767 if (calculationResult[eeFlatName]) {
768 calculationResult["selectee"] = 1;
769 }
770 }
771 }
772
773
774 if (calculationResult["nTrkLoose"] == 1 and calculationResult["maximumPCMS"] > 0.8 and selectedClusters.size() >= 2) {
775
776 decltype(selectedClusters) selectedSingleTagClusters(selectedClusters.size());
777 auto lastItem = std::copy_if(selectedClusters.begin(), selectedClusters.end(), selectedSingleTagClusters.begin(),
778 [](auto & cluster) {
779 const bool isNeutralCluster = not cluster.isTrack;
780 const bool hasEnoughEnergy = cluster.energyLab > 0.1;
781 const double clusterThetaLab = cluster.p4Lab.Theta() * TMath::RadToDeg();
782 const bool isInAcceptance = 17 < clusterThetaLab and clusterThetaLab < 150.;
783 return isNeutralCluster and hasEnoughEnergy and isInAcceptance;
784 });
785 selectedSingleTagClusters.resize(std::distance(selectedSingleTagClusters.begin(), lastItem));
786
787 if (selectedSingleTagClusters.size() >= 2) {
788
789 const auto& track = maximumPtTracks.at(-1) ? *maximumPtTracks.at(-1) : *maximumPtTracks.at(1);
790
791 const auto& firstCluster = selectedSingleTagClusters[0];
792 const auto& secondCluster = selectedSingleTagClusters[1];
793
794 const ROOT::Math::PxPyPzEVector trackP4CMS = track.p4CMS;
795 const ROOT::Math::PxPyPzEVector pi0P4CMS = firstCluster.
p4CMS + secondCluster.
p4CMS;
796
797 const bool passPi0ECMS = pi0P4CMS.E() > 1. and pi0P4CMS.E() < 0.525 * p4ofCOM.M();
798 const double thetaSumCMS = (pi0P4CMS.Theta() + trackP4CMS.Theta()) * TMath::RadToDeg();
799 const bool passThetaSum = thetaSumCMS < 170. or thetaSumCMS > 190.;
800
801 double dphiCMS = std::abs(trackP4CMS.Phi() - pi0P4CMS.Phi()) * TMath::RadToDeg();
802 if (dphiCMS > 180) {
803 dphiCMS = 360 - dphiCMS;
804 }
805 const bool passdPhi = dphiCMS > 160.;
806
807 if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() < 0.7) {
808 calculationResult["singleTagLowMass"] = 1;
809 } else if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() > 0.7) {
810 calculationResult["singleTagHighMass"] = 1;
811 }
812 }
813 }
814
815
816 if (calculationResult["nTrkLooseB"] == 1 and calculationResult["maximumPCMSB"] > 0.8 and selectedClusters.size() >= 2) {
817
818 decltype(selectedClusters) selectedSingleTagClusters(selectedClusters.size());
819 auto lastItem = std::copy_if(selectedClusters.begin(), selectedClusters.end(), selectedSingleTagClusters.begin(),
820 [](auto & cluster) {
821 const bool isNeutralCluster = not cluster.isTrack;
822 const bool hasEnoughEnergy = cluster.energyLab > 0.1;
823 const double clusterThetaLab = cluster.p4Lab.Theta() * TMath::RadToDeg();
824 const bool isInAcceptance = 17 < clusterThetaLab and clusterThetaLab < 150.;
825 return isNeutralCluster and hasEnoughEnergy and isInAcceptance;
826 });
827 selectedSingleTagClusters.resize(std::distance(selectedSingleTagClusters.begin(), lastItem));
828
829 if (selectedSingleTagClusters.size() >= 2) {
830
831 const auto& track = maximumPtTracks.at(-1) ? *maximumPtTracks.at(-1) : *maximumPtTracks.at(1);
832
833 const auto& firstCluster = selectedSingleTagClusters[0];
834 const auto& secondCluster = selectedSingleTagClusters[1];
835
836 const ROOT::Math::PxPyPzEVector trackP4CMS = track.p4CMS;
837 const ROOT::Math::PxPyPzEVector pi0P4CMS = firstCluster.
p4CMS + secondCluster.
p4CMS;
838
839 const bool passPi0ECMS = pi0P4CMS.E() > 1. and pi0P4CMS.E() < 0.525 * p4ofCOM.M();
840 const double thetaSumCMS = (pi0P4CMS.Theta() + trackP4CMS.Theta()) * TMath::RadToDeg();
841 const bool passThetaSum = thetaSumCMS < 170. or thetaSumCMS > 190.;
842
843 double dphiCMS = std::abs(trackP4CMS.Phi() - pi0P4CMS.Phi()) * TMath::RadToDeg();
844 if (dphiCMS > 180) {
845 dphiCMS = 360 - dphiCMS;
846 }
847 const bool passdPhi = dphiCMS > 160.;
848
849 if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() < 0.7) {
850 calculationResult["singleTagLowMassB"] = 1;
851 } else if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() > 0.7) {
852 calculationResult["singleTagHighMassB"] = 1;
853 }
854 }
855 }
856
857
859
860 const auto negTrack = maximumPtTracksWithoutZCut.at(-1);
861 const auto posTrack = maximumPtTracksWithoutZCut.at(1);
862
863 if (negTrack and posTrack) {
864
865 const double maxNegpT = negTrack->pT;
866 const double maxPospT = posTrack->pT;
867
868 auto accumulatePhotonEnergy = [](double result, const auto & cluster) {
871 };
872
873 const auto& clustersOfNegTrack = negTrack->track->getRelationsTo<
ECLCluster>();
875
876 const double maxClusterENeg = std::accumulate(clustersOfNegTrack.begin(), clustersOfNegTrack.end(), 0.0, accumulatePhotonEnergy);
877 const double maxClusterEPos = std::accumulate(clustersOfPosTrack.begin(), clustersOfPosTrack.end(), 0.0, accumulatePhotonEnergy);
878
879 const ROOT::Math::PxPyPzEVector& momentumLabNeg(negTrack->p4Lab);
880 const ROOT::Math::PxPyPzEVector& momentumLabPos(posTrack->p4Lab);
881
882 const double& z0Neg = negTrack->track->getTrackFitResultWithClosestMass(
Const::pion)->getZ0();
883 const double& d0Neg = negTrack->track->getTrackFitResultWithClosestMass(
Const::pion)->getD0();
884 const double& z0Pos = posTrack->track->getTrackFitResultWithClosestMass(
Const::pion)->getZ0();
885 const double& d0Pos = posTrack->track->getTrackFitResultWithClosestMass(
Const::pion)->getD0();
886
887
892 double dphiLab = std::abs(momentumLabNeg.Phi() - momentumLabPos.Phi()) * TMath::RadToDeg();
893 if (dphiLab > 180) {
894 dphiLab = 360 - dphiLab;
895 }
896
897 const double thetaSumLab = (momentumLabNeg.Theta() + momentumLabPos.Theta()) * TMath::RadToDeg();
898
899 constexpr double phiBackToBackTolerance = 2.;
900 constexpr double thetaBackToBackTolerance = 2.;
901 if ((180 - dphiLab) < phiBackToBackTolerance and std::abs(180 - thetaSumLab) < thetaBackToBackTolerance) {
902 calculationResult["cosmic"] = 1;
903 }
904 }
905 }
906 }
907
908
909
910 if (maximumPtTracksWithoutZCut.at(-1) and maximumPtTracksWithoutZCut.at(1)) {
911
912
913 const auto nTrack = maximumPtTracksWithoutZCut.at(-1)->track;
915
916 const auto pTrack = maximumPtTracksWithoutZCut.at(1)->track;
918
919
924
925
926 const double chisq = vertexFit->
getCHIsq();
927 const int ndf = vertexFit->
getNDF();
928 const double vertexProb = TMath::Prob(chisq, ndf);
929 const auto vertexLocation = vertexFit->
getVertex();
930 const double vertexXY = vertexLocation.perp();
931 const double vertexTheta = vertexLocation.theta() * TMath::RadToDeg();
932
933
934
935
936
937 const ROOT::Math::PxPyPzEVector& momentumLabNeg(maximumPtTracksWithoutZCut.at(-1)->p4Lab);
938 const ROOT::Math::PxPyPzEVector& momentumLabPos(maximumPtTracksWithoutZCut.at(1)->p4Lab);
939 const double thetaSumLab = (momentumLabNeg.Theta() + momentumLabPos.Theta()) * TMath::RadToDeg();
940 double dPhiLab = std::abs(momentumLabNeg.Phi() - momentumLabPos.Phi()) * TMath::RadToDeg();
941 if (dPhiLab > 180) {
942 dPhiLab = 360 - dPhiLab;
943 }
944 const double backToBackTolerance = 10.;
945 const bool backToBackLab = std::abs(thetaSumLab - 180.) < backToBackTolerance and std::abs(dPhiLab - 180.) < backToBackTolerance;
946
947
948 const double minProbChiVertex = 0.005;
949 const double minXYVertex = 3.;
950 const double maxXYVertex = 60.;
951 const double minThetaVertex = 30.;
952 const double maxThetaVertex = 120.;
953 if (vertexProb > minProbChiVertex and vertexXY > minXYVertex and vertexXY < maxXYVertex and vertexTheta > minThetaVertex
954 and vertexTheta < maxThetaVertex
955 and not backToBackLab) {calculationResult["displacedVertex"] = 1;}
956
957
958 delete nParticle;
959 delete pParticle;
960 delete vertexFit;
961
962 }
963
964}
DataType Angle(const B2Vector3< DataType > &q) const
The angle w.r.t.
Class to provide momentum-related information from ECLClusters.
const ROOT::Math::PxPyPzEVector Get4MomentumFromCluster(const ECLCluster *cluster, ECLCluster::EHypothesisBit hypo)
Returns four momentum vector.
const ROOT::Math::XYZVector GetIPPosition()
Returns default IP position from beam parameters.
static const ChargedStable pion
charged pion particle
@ c_nPhotons
CR is split into n photons (N1)
Class to store reconstructed particles.
RelationVector< TO > getRelationsTo(const std::string &name="", const std::string &namedRelation="") const
Get the relations that point from this object to another store array.
double m_EminLab
which lab energy defines nE180Lab
double m_reducedEsinglePhoton
which CMS energy defines nReducedEsingle clusters
StoreObjPtr< TRGSummary > m_l1Trigger
Store Object with the trigger result.
double m_E0min
which CMS energy defines nEmedium
double m_goodMagneticRegionZ0
maximum z0 for well understood magnetic field (cm)
double m_cosmicMinPt
which LAB pt defines a cosmic
double m_EsinglePhoton
which CMS energy defines nEsingleClust
double m_E2min
which CMS energy defines nElow
double m_tightTrkZ0
which Z0 defines a tight track
StoreArray< Track > m_tracks
Store Array of the tracks to be used.
double m_EminLab3Cluster
which lab energy defines nE500Lab
double m_Ehigh
which CMS energy defines nEhigh
StoreArray< ECLCluster > m_eclClusters
Store Array of the ecl clusters to be used.
double m_EminLab4Cluster
which lab energy defines nE300Lab
double m_cosmicMaxClusterEnergy
which LAB cluster energy vetoes a cosmic candidate
double m_goodMagneticRegionD0
minimum d0 for well understood magnetic field, if z0 is large (cm)
double m_looseTrkZ0
which Z0 defines a loose track
@ TTYP_DPHY
delayed physics events for background
@ TTYP_POIS
poisson random trigger
@ TTYP_RAND
random trigger events
Values of the result of a track fit with a given particle hypothesis.
Class that bundles various TrackFitResults.
virtual int getNDF(void) const
Get an NDF of the fit.
enum KFitError::ECode addParticle(const Particle *particle)
Add a particle to the fitter.
VertexFitKFit is a derived class from KFitBase to perform vertex-constraint kinematical fit.
double getCHIsq(void) const override
Get a chi-square of the fit.
enum KFitError::ECode doFit(void)
Perform a vertex-constraint fit.
const HepPoint3D getVertex(const int flag=KFitConst::kAfterFit) const
Get a vertex position.
B2Vector3< double > B2Vector3D
typedef for common usage with double
double charge(int pdgCode)
Returns electric charge of a particle with given pdg code.
Temporary data structure holding the track(s) with the maximum pT.
double pLab
the momentum magnitude in lab system
double pT
the pT of the track
ROOT::Math::PxPyPzEVector p4CMS
the 4 momentum in CMS system
double clusterEnergySumLab
the sum of related cluster energies in lab system
ROOT::Math::PxPyPzEVector p4Lab
the 4 momentum in lab system
double clusterEnergySumCMS
the sum of related cluster energies in CMS system
double pCMS
the momentum magnitude in CMS system
const Track * track
the track
Temporary data structure holding the ECL clusters used for this analysis.
double energyLab
the energy in Lab system
double clusterTime
the time of the cluster
ROOT::Math::PxPyPzEVector p4CMS
the 4 momentum in CMS system
ROOT::Math::PxPyPzEVector p4Lab
the 4 momentum in lab system
const ECLCluster * cluster
The ECL cluster.
bool isTrack
is this ECL cluster likely from a track (or a photon) = is it charged?
double energyCMS
the energy in CMS system