88 calculationResult[
"nTrkLoose"] = 0;
89 calculationResult[
"nTrkTight"] = 0;
90 calculationResult[
"ee2leg"] = 0;
91 calculationResult[
"nEmedium"] = 0;
92 calculationResult[
"nElow"] = 0;
93 calculationResult[
"nEhigh"] = 0;
94 calculationResult[
"nE180Lab"] = 0;
95 calculationResult[
"nE300Lab"] = 0;
96 calculationResult[
"nE500Lab"] = 0;
97 calculationResult[
"nE2000CMS"] = 0;
98 calculationResult[
"nE4000CMS"] = 0;
99 calculationResult[
"nE250Lab"] = 0;
100 calculationResult[
"nMaxEPhotonAcc"] = 0;
101 calculationResult[
"dphiCmsClust"] = NAN;
103 calculationResult[
"netChargeLoose"] = 0;
104 calculationResult[
"maximumPCMS"] = NAN;
105 calculationResult[
"maximumPLab"] = NAN;
106 calculationResult[
"eexx"] = 0;
107 calculationResult[
"ee1leg1trk"] = 0;
108 calculationResult[
"nEhighLowAng"] = 0;
109 calculationResult[
"nEsingleClust"] = 0;
110 calculationResult[
"nEsinglePhotonBarrel"] = 0;
111 calculationResult[
"nEsinglePhotonExtendedBarrel"] = 0;
112 calculationResult[
"nEsinglePhotonEndcap"] = 0;
113 calculationResult[
"nEsingleElectronBarrel"] = 0;
114 calculationResult[
"nEsingleElectronExtendedBarrel"] = 0;
115 calculationResult[
"nReducedEsinglePhotonReducedBarrel"] = 0;
116 calculationResult[
"nVetoClust"] = 0;
117 calculationResult[
"chrgClust2GeV"] = 0;
118 calculationResult[
"neutClust045GeVAcc"] = 0;
119 calculationResult[
"neutClust045GeVBarrel"] = 0;
120 calculationResult[
"singleTagLowMass"] = 0;
121 calculationResult[
"singleTagHighMass"] = 0;
122 calculationResult[
"n2GeVNeutBarrel"] = 0;
123 calculationResult[
"n2GeVNeutEndcap"] = 0;
124 calculationResult[
"n2GeVChrg"] = 0;
125 calculationResult[
"n2GeVPhotonBarrel"] = 0;
126 calculationResult[
"n2GeVPhotonEndcap"] = 0;
127 calculationResult[
"ee1leg"] = 0;
128 calculationResult[
"ee1leg1clst"] = 0;
129 calculationResult[
"ee1leg1e"] = 0;
130 calculationResult[
"ee2clst"] = 0;
131 calculationResult[
"gg2clst"] = 0;
132 calculationResult[
"eeee"] = 0;
133 calculationResult[
"eemm"] = 0;
134 calculationResult[
"eexxSelect"] = 0;
135 calculationResult[
"radBhabha"] = 0;
136 calculationResult[
"eeBrem"] = 0;
137 calculationResult[
"isrRadBhabha"] = 0;
138 calculationResult[
"muonPairECL"] = 0;
139 calculationResult[
"ggHighPt"] = 0;
140 calculationResult[
"selectee1leg1trk"] = 0;
141 calculationResult[
"selectee1leg1clst"] = 0;
142 calculationResult[
"selectee"] = 0;
143 calculationResult[
"ggBarrelVL"] = 0;
144 calculationResult[
"ggBarrelLoose"] = 0;
145 calculationResult[
"ggBarrelTight"] = 0;
146 calculationResult[
"ggEndcapVL"] = 0;
147 calculationResult[
"ggEndcapLoose"] = 0;
148 calculationResult[
"ggEndcapTight"] = 0;
149 calculationResult[
"muonPairV"] = 0;
150 calculationResult[
"selectmumu"] = 0;
151 calculationResult[
"singleMuon"] = 0;
152 calculationResult[
"cosmic"] = 0;
153 calculationResult[
"displacedVertex"] = 0;
154 calculationResult[
"eeFlat0"] = 0;
155 calculationResult[
"eeFlat1"] = 0;
156 calculationResult[
"eeFlat2"] = 0;
157 calculationResult[
"eeFlat3"] = 0;
158 calculationResult[
"eeFlat4"] = 0;
159 calculationResult[
"eeFlat5"] = 0;
160 calculationResult[
"eeFlat6"] = 0;
161 calculationResult[
"eeFlat7"] = 0;
162 calculationResult[
"eeFlat8"] = 0;
163 calculationResult[
"eeOneClust"] = 0;
166 calculationResult[
"nTrkLooseB"] = 0;
167 calculationResult[
"nTrkTightB"] = 0;
168 calculationResult[
"maximumPCMSB"] = NAN;
169 calculationResult[
"netChargeLooseB"] = 0;
170 calculationResult[
"muonPairVB"] = 0;
171 calculationResult[
"eeBremB"] = 0;
172 calculationResult[
"singleTagLowMassB"] = 0;
173 calculationResult[
"singleTagHighMassB"] = 0;
174 calculationResult[
"radBhabhaB"] = 0;
177 calculationResult[
"nTrackC"] = 0;
178 calculationResult[
"maximumPCMSC"] = NAN;
179 calculationResult[
"pCmsNegC"] = NAN;
180 calculationResult[
"clusterENegC"] = NAN;
181 calculationResult[
"pCmsPosC"] = NAN;
182 calculationResult[
"clusterEPosC"] = NAN;
183 calculationResult[
"dPhiCmsC"] = NAN;
194 }
catch (
const std::exception&) {
200 }
catch (
const std::exception&) {
206 }
catch (
const std::exception&) {
212 }
catch (
const std::exception&) {
218 }
catch (
const std::exception&) {
221 calculationResult[
"bha3d"] = bha3d;
222 calculationResult[
"bhapur"] = bhapurPsnm;
223 calculationResult[
"bhapur_lml1"] = lml1 and bhapurFtdl;
224 calculationResult[
"l1_bit_f"] = l1_bit_f;
226 calculationResult[
"l1_trigger_random"] = 1;
227 calculationResult[
"l1_trigger_delayed_bhabha"] = 0;
228 calculationResult[
"l1_trigger_poisson"] = 0;
229 calculationResult[
"bha3d"] = 0;
230 calculationResult[
"bhapur"] = 0;
231 calculationResult[
"bhapur_lml1"] = 0;
232 calculationResult[
"l1_bit_f"] = 0;
236 calculationResult[
"l1_trg_NN_info"] = 0;
237 if (
m_bitsNN.isValid() and
m_bitsNN.getEntries() > 0) {calculationResult[
"l1_trg_NN_info"] = 1;}
239 calculationResult[
"true"] = 1;
240 calculationResult[
"false"] = 0;
245 const ROOT::Math::XYZVector clustervertex = cUtil.
GetIPPosition();
247 ROOT::Math::PxPyPzEVector p4ofCOM;
248 p4ofCOM.SetPxPyPzE(0, 0, 0, boostrotate.
getCMSEnergy());
251 std::map<short, std::optional<MaximumPtTrack>> maximumPtTracks = {
256 std::map<short, std::optional<MaximumPtTrack>> maximumPtTracksWithoutZCut = {
261 std::map<short, std::optional<MaximumPtTrack>> maximumPCmsTracksC = {
269 if (not trackFitResult) {
276 if (charge == 0) {
continue;}
277 const double z0 = trackFitResult->
getZ0();
278 const ROOT::Math::PxPyPzEVector& momentumLab = trackFitResult->
get4Momentum();
279 const ROOT::Math::PxPyPzEVector momentumCMS = boostrotate.
rotateLabToCms() * momentumLab;
280 const double pCMS = momentumCMS.P();
281 const double pLab = momentumLab.P();
282 const double trackTime = track.getTrackTime();
285 double clusterELab = 0.;
286 for (
auto& cluster : track.getRelationsTo<
ECLCluster>()) {
291 const double clusterECMS = clusterELab * pCMS / pLab;
296 bool goodTrackCTime =
true;
297 if (std::abs(trackTime) > 15.) {goodTrackCTime =
false;}
298 if (std::abs(z0) < 1. and goodTrackCTime and pCMS > 0.2) {
299 calculationResult[
"nTrackC"] += 1;
300 if (std::isnan(calculationResult[
"maximumPCMSC"]) or pCMS > calculationResult[
"maximumPCMSC"]) {
301 calculationResult[
"maximumPCMSC"] = pCMS;
306 const auto& currentMaximumC = maximumPCmsTracksC.at(charge);
307 if (not currentMaximumC or pCMS > currentMaximumC->pCMS) {
310 newMaximum.
track = &track;
311 newMaximum.
pCMS = pCMS;
312 newMaximum.
pLab = pLab;
313 newMaximum.
p4CMS = momentumCMS;
314 newMaximum.
p4Lab = momentumLab;
317 maximumPCmsTracksC[charge] = newMaximum;
329 calculationResult[
"nTrkTight"] += 1;
335 const auto& currentMaximum = maximumPtTracksWithoutZCut.at(charge);
336 if (not currentMaximum or pT > currentMaximum->pT) {
339 newMaximum.
track = &track;
340 newMaximum.
pCMS = pCMS;
341 newMaximum.
pLab = pLab;
342 newMaximum.
p4CMS = momentumCMS;
343 newMaximum.
p4Lab = momentumLab;
346 maximumPtTracksWithoutZCut[charge] = newMaximum;
353 calculationResult[
"nTrkLoose"] += 1;
354 calculationResult[
"netChargeLoose"] += charge;
356 if (std::isnan(calculationResult[
"maximumPCMS"]) or pCMS > calculationResult[
"maximumPCMS"]) {
357 calculationResult[
"maximumPCMS"] = pCMS;
360 if (std::isnan(calculationResult[
"maximumPLab"]) or pLab > calculationResult[
"maximumPLab"]) {
361 calculationResult[
"maximumPLab"] = pLab;
366 const auto& currentMaximumLoose = maximumPtTracks.at(charge);
367 if (not currentMaximumLoose or pTLoose > currentMaximumLoose->pT) {
369 newMaximum.
pT = pTLoose;
370 newMaximum.
track = &track;
371 newMaximum.
pCMS = pCMS;
372 newMaximum.
pLab = momentumLab.P();
373 newMaximum.
p4CMS = momentumCMS;
374 newMaximum.
p4Lab = momentumLab;
377 maximumPtTracks[charge] = newMaximum;
384 if (std::abs(z0) < 1.) {calculationResult[
"nTrkTightB"] += 1;}
385 if (std::abs(z0) < 5.) {
386 calculationResult[
"nTrkLooseB"] += 1;
387 calculationResult[
"netChargeLooseB"] += charge;
388 if (std::isnan(calculationResult[
"maximumPCMSB"]) or pCMS > calculationResult[
"maximumPCMSB"]) {
389 calculationResult[
"maximumPCMSB"] = pCMS;
400 for (
short charge : { -1, 1}) {
401 auto& maximumPcmsTrackC = maximumPCmsTracksC.at(charge);
402 if (not maximumPcmsTrackC) {
408 calculationResult[
"pCmsNegC"] = maximumPcmsTrackC->pCMS;
409 calculationResult[
"clusterENegC"] = maximumPcmsTrackC->clusterEnergySumLab;
411 calculationResult[
"pCmsPosC"] = maximumPcmsTrackC->pCMS;
412 calculationResult[
"clusterEPosC"] = maximumPcmsTrackC->clusterEnergySumLab;
417 if (maximumPCmsTracksC.at(-1) and maximumPCmsTracksC.at(1)) {
421 calculationResult[
"dPhiCmsC"] = std::abs(ROOT::Math::VectorUtil::DeltaPhi(negativeTrack.
p4CMS,
422 positiveTrack.
p4CMS)) * TMath::RadToDeg();
428 std::vector<SelectedECLCluster> selectedClusters;
431 const double time = cluster.getTime();
433 std::abs(time) > 200 or
438 const double dt99 = cluster.getDeltaTime99();
445 selectedCluster.
cluster = &cluster;
449 selectedCluster.
isTrack = cluster.isTrack();
451 selectedClusters.push_back(selectedCluster);
454 calculationResult[
"nElow"] += 1;
457 calculationResult[
"nEmedium"] += 1;
460 calculationResult[
"nEhigh"] += 1;
464 calculationResult[
"nVetoClust"] += 1;
469 const double thetaLab = selectedCluster.
p4Lab.Theta() * TMath::RadToDeg();
470 const double zmva = cluster.getZernikeMVA();
471 const bool photon = zmva > 0.5 and not selectedCluster.
isTrack;
472 const bool electron = zmva > 0.5 and selectedCluster.
isTrack;
476 calculationResult[
"chrgClust2GeV"] += 1;
479 const bool isInAcceptance = 17. < thetaLab and thetaLab < 150.;
480 if (isInAcceptance) {calculationResult[
"neutClust045GeVAcc"] += 1;}
481 const bool isInBarrel = 30. < thetaLab and thetaLab < 130.;
482 if (isInBarrel) {calculationResult[
"neutClust045GeVBarrel"] += 1;}
487 const bool notInHighBackgroundEndcapRegion = 18.5 < thetaLab and thetaLab < 139.3;
489 calculationResult[
"nE180Lab"] += 1;
493 calculationResult[
"nE300Lab"] += 1;
497 calculationResult[
"nE500Lab"] += 1;
500 if (selectedCluster.
energyCMS >
m_Ehigh and notInHighBackgroundEndcapRegion) {
501 calculationResult[
"nE2000CMS"] += 1;
506 calculationResult[
"nE250Lab"] += 1;
509 calculationResult[
"nE4000CMS"] += 1;
514 calculationResult[
"nEsingleClust"] += 1;
516 const bool barrelRegion = thetaLab > 45 and thetaLab < 115;
517 const bool extendedBarrelRegion = thetaLab > 30 and thetaLab < 130;
518 const bool endcapRegion = (thetaLab > 22 and thetaLab < 45) or (thetaLab > 115 and thetaLab < 145);
520 if (photon and barrelRegion) {
521 calculationResult[
"nEsinglePhotonBarrel"] += 1;
524 if (photon and extendedBarrelRegion) {
525 calculationResult[
"nEsinglePhotonExtendedBarrel"] += 1;
528 if (electron and barrelRegion) {
529 calculationResult[
"nEsingleElectronBarrel"] += 1;
532 if (electron and extendedBarrelRegion) {
533 calculationResult[
"nEsingleElectronExtendedBarrel"] += 1;
536 if (photon and endcapRegion) {
537 calculationResult[
"nEsinglePhotonEndcap"] += 1;
542 const bool reducedBarrelRegion = thetaLab > 44 and thetaLab < 98;
544 if (photon and reducedBarrelRegion) {
545 calculationResult[
"nReducedEsinglePhotonReducedBarrel"] += 1;
551 const bool barrelRegion = thetaLab > 32 and thetaLab < 130;
552 const bool endcapRegion = (thetaLab > 22 and thetaLab < 32) or (thetaLab > 130 and thetaLab < 145);
553 const bool lowAngleRegion = thetaLab < 22 or thetaLab > 145;
555 if (not selectedCluster.
isTrack and barrelRegion) {
556 calculationResult[
"n2GeVNeutBarrel"] += 1;
558 if (not selectedCluster.
isTrack and endcapRegion) {
559 calculationResult[
"n2GeVNeutEndcap"] += 1;
561 if (selectedCluster.
isTrack and not lowAngleRegion) {
562 calculationResult[
"n2GeVChrg"] += 1;
564 if (lowAngleRegion) {
565 calculationResult[
"nEhighLowAng"] += 1;
567 if (photon and barrelRegion) {
568 calculationResult[
"n2GeVPhotonBarrel"] += 1;
570 if (photon and endcapRegion) {
571 calculationResult[
"n2GeVPhotonEndcap"] += 1;
577 std::sort(selectedClusters.begin(), selectedClusters.end(), [](
const auto & lhs,
const auto & rhs) {
578 return lhs.energyCMS > rhs.energyCMS;
582 for (
short charge : { -1, 1}) {
583 auto& maximumPtTrack = maximumPtTracks.at(charge);
584 if (not maximumPtTrack) {
589 if (maximumPtTrack->clusterEnergySumCMS > 4.5) {
590 calculationResult[
"ee1leg"] = 1;
594 if (maximumPtTrack->pCMS > 3 and maximumPtTrack->clusterEnergySumLab > 0. and maximumPtTrack->clusterEnergySumLab < 1.) {
595 calculationResult[
"singleMuon"] = 1;
600 if (maximumPtTracks.at(-1) and maximumPtTracks.at(1)) {
604 double dphi = std::abs(negativeTrack.
p4CMS.Phi() - positiveTrack.
p4CMS.Phi()) * TMath::RadToDeg();
611 const double negativeP = negativeTrack.
pCMS;
614 const double positiveP = positiveTrack.
pCMS;
617 const double thetaSum = (negativeTrack.
p4CMS.Theta() + positiveTrack.
p4CMS.Theta()) * TMath::RadToDeg();
618 const double dthetaSum = std::abs(thetaSum - 180);
619 const auto back2back = dphi > 175 and dthetaSum < 15;
620 if (back2back and negativeClusterSum > 3 and positiveClusterSum > 3 and
621 (negativeClusterSum > 4.5 or positiveClusterSum > 4.5)) {
622 calculationResult[
"ee2leg"] = 1;
626 if (back2back and ((negativeClusterSum > 4.5 and positiveP > 3) or (positiveClusterSum > 4.5 and negativeP > 3))) {
627 calculationResult[
"ee1leg1trk"] = 1;
632 if ((negativeClusterSum > 4.5 and positiveClusterSum > 0.8 * positiveP) or
633 (positiveClusterSum > 4.5 and negativeClusterSum > 0.8 * negativeP)) {
634 calculationResult[
"ee1leg1e"] = 1;
638 const ROOT::Math::PxPyPzEVector p4Miss = p4ofCOM - negativeTrack.
p4CMS - positiveTrack.
p4CMS;
639 const double pmissTheta = p4Miss.Theta() * TMath::RadToDeg();
640 const double pmissp = p4Miss.P();
641 const double relMissAngle0 = ROOT::Math::VectorUtil::Angle(negativeTrack.
p4CMS, p4Miss) * TMath::RadToDeg();
642 const double relMissAngle1 = ROOT::Math::VectorUtil::Angle(positiveTrack.
p4CMS, p4Miss) * TMath::RadToDeg();
644 const bool electronEP = positiveClusterSum > 0.8 * positiveP or negativeClusterSum > 0.8 * negativeP;
645 const bool notMuonPair = negativeClusterSumLab > 1 or positiveClusterSumLab > 1;
646 const double highp = std::max(negativeP, positiveP);
647 const double lowp = std::min(negativeP, positiveP);
648 const bool lowEdep = negativeClusterSumLab < 0.5 and positiveClusterSumLab < 0.5;
651 if (calculationResult[
"maximumPCMS"] < 2 and dphi > 160 and (pmissTheta < 25. or pmissTheta > 155.)) {
652 calculationResult[
"eexxSelect"] = 1;
654 calculationResult[
"eeee"] = 1;
656 calculationResult[
"eemm"] = 1;
661 if ((pmissTheta < 20. or pmissTheta > 160.) and
662 ((calculationResult[
"maximumPCMS"] < 1.2 and dphi > 150.) or
663 (calculationResult[
"maximumPCMS"] < 2. and 175. < dphi))) {
664 calculationResult[
"eexx"] = 1;
668 if (negativeP > 1. and pmissTheta > 10. and pmissTheta < 170. and positiveP > 1. and dphi < 170. and pmissp > 1. and electronEP) {
669 if (calculationResult[
"nTrkLoose"] == 2 and calculationResult[
"nTrkTight"] >= 1) {
670 calculationResult[
"radBhabha"] = 1;
672 if (calculationResult[
"nTrkLooseB"] == 2 and calculationResult[
"nTrkTightB"] >= 1) {
673 calculationResult[
"radBhabhaB"] = 1;
678 if (negativeP > 2. and positiveP > 2. and 2 == calculationResult[
"nTrkLoose"] and
679 calculationResult[
"nTrkTight"] >= 1 and dphi > 175. and
680 (pmissTheta < 5. or pmissTheta > 175.) and electronEP) {
681 calculationResult[
"isrRadBhabha"] = 1;
685 if (highp > 4.5 and notMuonPair and pmissp > 1. and (relMissAngle0 < 5. or relMissAngle1 < 5.)) {
686 if (calculationResult[
"nTrkLoose"] == 2) { calculationResult[
"eeBrem"] = 1;}
687 if (calculationResult[
"nTrkLooseB"] >= 1) { calculationResult[
"eeBremB"] = 1;}
691 if (highp > 4.5 and lowEdep and thetaSum > 175. and thetaSum < 185. and dphi > 175.) {
692 if (calculationResult[
"nTrkLoose"] == 2) {calculationResult[
"muonPairV"] = 1;}
693 if (calculationResult[
"nTrkLooseB"] == 2) {calculationResult[
"muonPairVB"] = 1;}
697 if (highp > 3. and lowp > 2.5 and dphi > 165. and
698 ((negativeClusterSumLab > 0. and negativeClusterSumLab < 1.) or
699 (positiveClusterSumLab > 0. and positiveClusterSumLab < 1.))) {
700 calculationResult[
"selectmumu"] = 1;
705 if (selectedClusters.size() >= 2) {
709 double dphi = std::abs(firstCluster.
p4CMS.Phi() - secondCluster.
p4CMS.Phi()) * TMath::RadToDeg();
713 double thetaSum = (firstCluster.
p4CMS.Theta() + secondCluster.
p4CMS.Theta()) * TMath::RadToDeg();
714 double dthetaSum = std::abs(thetaSum - 180);
717 calculationResult[
"dphiCmsClust"] = dphi;
718 for (
int ic = 0; ic < 2; ic++) {
719 const double thetaLab = selectedClusters[ic].p4Lab.Theta() * TMath::RadToDeg();
720 const bool isInAcceptance = 17. < thetaLab and thetaLab < 150.;
721 const ECLCluster* cluster = selectedClusters[ic].cluster;
722 const double zmva = cluster->getZernikeMVA();
723 const bool photon = zmva > 0.5 and not selectedClusters[ic].isTrack;
724 if (isInAcceptance and photon) {calculationResult[
"nMaxEPhotonAcc"] += 1;}
727 const double firstEnergy = firstCluster.
p4CMS.E();
728 const double secondEnergy = secondCluster.
p4CMS.E();
730 const bool highEnergetic = firstEnergy > 3 and secondEnergy > 3 and (firstEnergy > 4.5 or secondEnergy > 4.5);
732 if (dphi > 160 and dphi < 178 and dthetaSum < 15 and highEnergetic) {
733 calculationResult[
"ee2clst"] = 1;
736 if (dphi > 178 and dthetaSum < 15 and highEnergetic) {
737 calculationResult[
"gg2clst"] = 1;
740 if ((calculationResult[
"ee2clst"] == 1 or calculationResult[
"gg2clst"] == 1) and
741 calculationResult[
"ee1leg"] == 1) {
742 calculationResult[
"ee1leg1clst"] = 1;
745 const double Elab0 = firstCluster.
p4Lab.E();
746 const double Elab1 = secondCluster.
p4Lab.E();
749 if (firstEnergy > 2 and secondEnergy > 2) {
750 const double thetaLab0 = firstCluster.
p4Lab.Theta() * TMath::RadToDeg();
751 const double thetaLab1 = secondCluster.
p4Lab.Theta() * TMath::RadToDeg();
753 const bool barrel0 = 32. < thetaLab0 and thetaLab0 < 130.;
754 const bool barrel1 = 32. < thetaLab1 and thetaLab1 < 130.;
755 const bool oneClustersAbove4 = firstEnergy > 4 or secondEnergy > 4;
756 const bool oneIsNeutral = not firstCluster.
isTrack or not secondCluster.
isTrack;
757 const bool bothAreNeutral = not firstCluster.
isTrack and not secondCluster.
isTrack;
758 const bool oneIsBarrel = barrel0 or barrel1;
759 const bool dphiCutExtraLoose = dphi > 175;
760 const bool dphiCutLoose = dphi > 177;
761 const bool dphiCutTight = dphi > 177.5;
763 if (dphiCutExtraLoose and oneIsNeutral and oneIsBarrel) {
764 calculationResult[
"ggBarrelVL"] = 1;
766 if (oneClustersAbove4 and dphiCutLoose and oneIsNeutral and oneIsBarrel) {
767 calculationResult[
"ggBarrelLoose"] = 1;
769 if (oneClustersAbove4 and dphiCutTight and bothAreNeutral and oneIsBarrel) {
770 calculationResult[
"ggBarrelTight"] = 1;
772 if (dphiCutExtraLoose and oneIsNeutral and not oneIsBarrel) {
773 calculationResult[
"ggEndcapVL"] = 1;
775 if (oneClustersAbove4 and dphiCutLoose and oneIsNeutral and not oneIsBarrel) {
776 calculationResult[
"ggEndcapLoose"] = 1;
778 if (oneClustersAbove4 and dphiCutTight and bothAreNeutral and not oneIsBarrel) {
779 calculationResult[
"ggEndcapTight"] = 1;
783 const double minEnergy = std::min(Elab0, Elab1);
784 const double maxEnergy = std::max(Elab0, Elab1);
785 if (dphi > 155 and thetaSum > 165 and thetaSum < 195 and minEnergy > 0.15 and minEnergy < 0.5 and
786 maxEnergy > 0.15 and maxEnergy < 0.5) {
787 calculationResult[
"muonPairECL"] = 1;
791 const double thetaLab0 = firstCluster.
p4Lab.Theta() * TMath::RadToDeg();
792 const double thetaLab1 = secondCluster.
p4Lab.Theta() * TMath::RadToDeg();
793 const bool inHieRegion0 = thetaLab0 > 26. and thetaLab0 < 130.;
794 const bool inHieRegion1 = thetaLab1 > 26. and thetaLab1 < 130.;
795 const bool firstIsNeutral = not firstCluster.
isTrack;
796 const bool secondIsNeutral = not secondCluster.
isTrack;
798 if (secondEnergy > 0.3 and inHieRegion0 and inHieRegion1 and firstIsNeutral and secondIsNeutral) {
799 const ROOT::Math::PxPyPzEVector ggP4CMS = firstCluster.
p4CMS + secondCluster.
p4CMS;
800 if (ggP4CMS.pt() > 1.) {calculationResult[
"ggHighPt"] = 1;}
807 double thetaFlatten = 0;
810 for (
short charge : { -1, 1}) {
811 const auto& maximumPtTrack = maximumPtTracks.at(charge);
812 if (not maximumPtTrack) {
816 if (maximumPtTrack->clusterEnergySumCMS > 1.5) {
818 double tempFlatten = 0.;
820 double tempInvMass = (maximumPtTrack->p4Lab + cluster.p4Lab).M();
821 if (tempInvMass > invMass) {
822 invMass = tempInvMass;
824 tempFlatten = cluster.p4Lab.Theta() * TMath::RadToDeg();
829 tempFlatten = maximumPtTrack->p4Lab.Theta() * TMath::RadToDeg();
831 if (invMass > 5.29) {
832 calculationResult[
"selectee1leg1clst"] = 1;
833 thetaFlatten = tempFlatten;
839 if (maximumPtTracks.at(-1) and maximumPtTracks.at(1)) {
842 const double invMass = (negativeTrack.
p4Lab + positiveTrack.
p4Lab).M();
844 calculationResult[
"selectee1leg1trk"] = 1;
847 thetaFlatten = negativeTrack.
p4Lab.Theta() * TMath::RadToDeg();
852 if ((invMass > 9.) and (missNegClust or missPosClust)) {
853 calculationResult[
"eeOneClust"] = 1;
857 if (calculationResult[
"selectee1leg1trk"] == 1 or calculationResult[
"selectee1leg1clst"] == 1) {
858 for (
int iflat = 0; iflat < 9; iflat++) {
859 const std::string& eeFlatName =
"eeFlat" + std::to_string(iflat);
860 calculationResult[eeFlatName] =
861 thetaFlatten >= flatBoundaries[iflat] and thetaFlatten < flatBoundaries[iflat + 1];
862 if (calculationResult[eeFlatName]) {
863 calculationResult[
"selectee"] = 1;
869 if (calculationResult[
"nTrkLoose"] == 1 and calculationResult[
"maximumPCMS"] > 0.8 and selectedClusters.size() >= 2) {
871 decltype(selectedClusters) selectedSingleTagClusters(selectedClusters.size());
872 auto lastItem = std::copy_if(selectedClusters.begin(), selectedClusters.end(), selectedSingleTagClusters.begin(),
874 const bool isNeutralCluster = not cluster.isTrack;
875 const bool hasEnoughEnergy = cluster.energyLab > 0.1;
876 const double clusterThetaLab = cluster.p4Lab.Theta() * TMath::RadToDeg();
877 const bool isInAcceptance = 17 < clusterThetaLab and clusterThetaLab < 150.;
878 return isNeutralCluster and hasEnoughEnergy and isInAcceptance;
880 selectedSingleTagClusters.resize(std::distance(selectedSingleTagClusters.begin(), lastItem));
882 if (selectedSingleTagClusters.size() >= 2) {
884 const auto& track = maximumPtTracks.at(-1) ? *maximumPtTracks.at(-1) : *maximumPtTracks.at(1);
886 const auto& firstCluster = selectedSingleTagClusters[0];
887 const auto& secondCluster = selectedSingleTagClusters[1];
889 const ROOT::Math::PxPyPzEVector trackP4CMS = track.p4CMS;
890 const ROOT::Math::PxPyPzEVector pi0P4CMS = firstCluster.p4CMS + secondCluster.p4CMS;
892 const bool passPi0ECMS = pi0P4CMS.E() > 1. and pi0P4CMS.E() < 0.525 * p4ofCOM.M();
893 const double thetaSumCMS = (pi0P4CMS.Theta() + trackP4CMS.Theta()) * TMath::RadToDeg();
894 const bool passThetaSum = thetaSumCMS < 170. or thetaSumCMS > 190.;
896 double dphiCMS = std::abs(trackP4CMS.Phi() - pi0P4CMS.Phi()) * TMath::RadToDeg();
898 dphiCMS = 360 - dphiCMS;
900 const bool passdPhi = dphiCMS > 160.;
902 if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() < 0.7) {
903 calculationResult[
"singleTagLowMass"] = 1;
904 }
else if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() > 0.7) {
905 calculationResult[
"singleTagHighMass"] = 1;
911 if (calculationResult[
"nTrkLooseB"] == 1 and calculationResult[
"maximumPCMSB"] > 0.8 and selectedClusters.size() >= 2) {
913 decltype(selectedClusters) selectedSingleTagClusters(selectedClusters.size());
914 auto lastItem = std::copy_if(selectedClusters.begin(), selectedClusters.end(), selectedSingleTagClusters.begin(),
916 const bool isNeutralCluster = not cluster.isTrack;
917 const bool hasEnoughEnergy = cluster.energyLab > 0.1;
918 const double clusterThetaLab = cluster.p4Lab.Theta() * TMath::RadToDeg();
919 const bool isInAcceptance = 17 < clusterThetaLab and clusterThetaLab < 150.;
920 return isNeutralCluster and hasEnoughEnergy and isInAcceptance;
922 selectedSingleTagClusters.resize(std::distance(selectedSingleTagClusters.begin(), lastItem));
924 if (selectedSingleTagClusters.size() >= 2) {
926 const auto& track = maximumPtTracks.at(-1) ? *maximumPtTracks.at(-1) : *maximumPtTracks.at(1);
928 const auto& firstCluster = selectedSingleTagClusters[0];
929 const auto& secondCluster = selectedSingleTagClusters[1];
931 const ROOT::Math::PxPyPzEVector trackP4CMS = track.p4CMS;
932 const ROOT::Math::PxPyPzEVector pi0P4CMS = firstCluster.p4CMS + secondCluster.p4CMS;
934 const bool passPi0ECMS = pi0P4CMS.E() > 1. and pi0P4CMS.E() < 0.525 * p4ofCOM.M();
935 const double thetaSumCMS = (pi0P4CMS.Theta() + trackP4CMS.Theta()) * TMath::RadToDeg();
936 const bool passThetaSum = thetaSumCMS < 170. or thetaSumCMS > 190.;
938 double dphiCMS = std::abs(trackP4CMS.Phi() - pi0P4CMS.Phi()) * TMath::RadToDeg();
940 dphiCMS = 360 - dphiCMS;
942 const bool passdPhi = dphiCMS > 160.;
944 if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() < 0.7) {
945 calculationResult[
"singleTagLowMassB"] = 1;
946 }
else if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() > 0.7) {
947 calculationResult[
"singleTagHighMassB"] = 1;
955 const auto negTrack = maximumPtTracksWithoutZCut.at(-1);
956 const auto posTrack = maximumPtTracksWithoutZCut.at(1);
958 if (negTrack and posTrack) {
960 const double maxNegpT = negTrack->pT;
961 const double maxPospT = posTrack->pT;
962 const double maxClusterENeg = negTrack->clusterEnergySumLab;
963 const double maxClusterEPos = posTrack->clusterEnergySumLab;
965 const ROOT::Math::PxPyPzEVector& momentumLabNeg(negTrack->p4Lab);
966 const ROOT::Math::PxPyPzEVector& momentumLabPos(posTrack->p4Lab);
968 const double& z0Neg = negTrack->track->getTrackFitResultWithClosestMass(
Const::pion)->getZ0();
969 const double& d0Neg = negTrack->track->getTrackFitResultWithClosestMass(
Const::pion)->getD0();
970 const double& z0Pos = posTrack->track->getTrackFitResultWithClosestMass(
Const::pion)->getZ0();
971 const double& d0Pos = posTrack->track->getTrackFitResultWithClosestMass(
Const::pion)->getD0();
978 double dphiLab = std::abs(momentumLabNeg.Phi() - momentumLabPos.Phi()) * TMath::RadToDeg();
980 dphiLab = 360 - dphiLab;
983 const double thetaSumLab = (momentumLabNeg.Theta() + momentumLabPos.Theta()) * TMath::RadToDeg();
985 constexpr double phiBackToBackTolerance = 2.;
986 constexpr double thetaBackToBackTolerance = 2.;
987 if ((180 - dphiLab) < phiBackToBackTolerance and std::abs(180 - thetaSumLab) < thetaBackToBackTolerance) {
988 calculationResult[
"cosmic"] = 1;
996 if (maximumPtTracksWithoutZCut.at(-1) and maximumPtTracksWithoutZCut.at(1)) {
999 const auto nTrack = maximumPtTracksWithoutZCut.at(-1)->track;
1002 const auto pTrack = maximumPtTracksWithoutZCut.at(1)->track;
1012 const double chisq = vertexFit->
getCHIsq();
1013 const int ndf = vertexFit->
getNDF();
1014 const double vertexProb = TMath::Prob(chisq, ndf);
1015 const auto vertexLocation = vertexFit->
getVertex();
1016 const double vertexXY = vertexLocation.perp();
1017 const double vertexTheta = vertexLocation.theta() * TMath::RadToDeg();
1023 const ROOT::Math::PxPyPzEVector& momentumLabNeg(maximumPtTracksWithoutZCut.at(-1)->p4Lab);
1024 const ROOT::Math::PxPyPzEVector& momentumLabPos(maximumPtTracksWithoutZCut.at(1)->p4Lab);
1025 const double thetaSumLab = (momentumLabNeg.Theta() + momentumLabPos.Theta()) * TMath::RadToDeg();
1026 double dPhiLab = std::abs(momentumLabNeg.Phi() - momentumLabPos.Phi()) * TMath::RadToDeg();
1027 if (dPhiLab > 180) {
1028 dPhiLab = 360 - dPhiLab;
1030 const double backToBackTolerance = 10.;
1031 const bool backToBackLab = std::abs(thetaSumLab - 180.) < backToBackTolerance and std::abs(dPhiLab - 180.) < backToBackTolerance;
1034 const double minProbChiVertex = 0.005;
1035 const double minXYVertex = 3.;
1036 const double maxXYVertex = 60.;
1037 const double minThetaVertex = 30.;
1038 const double maxThetaVertex = 120.;
1039 if (vertexProb > minProbChiVertex and vertexXY > minXYVertex and vertexXY < maxXYVertex and vertexTheta > minThetaVertex
1040 and vertexTheta < maxThetaVertex
1041 and not backToBackLab) {calculationResult[
"displacedVertex"] = 1;}