Belle II Software  release-05-01-25
FilterCalculator.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2019 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Chris Hearty, Nils Braun *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 #include <hlt/softwaretrigger/calculations/FilterCalculator.h>
11 #include <hlt/softwaretrigger/calculations/utilities.h>
12 
13 #include <TLorentzVector.h>
14 
15 #include <mdst/dataobjects/TrackFitResult.h>
16 #include <mdst/dataobjects/HitPatternCDC.h>
17 
18 #include <numeric>
19 
20 using namespace Belle2;
21 using namespace SoftwareTrigger;
22 
24 struct MaximumPtTrack {
26  double pT = NAN;
28  const Track* track = nullptr;
30  double clusterEnergySumCMS = 0;
32  double clusterEnergySumLab = 0;
34  double pCMS = NAN;
36  double pLab = NAN;
38  TLorentzVector p4CMS;
40  TLorentzVector p4Lab;
41 };
42 
44 struct SelectedECLCluster {
46  const ECLCluster* cluster = nullptr;
48  double energyCMS = NAN;
50  double energyLab = NAN;
52  TLorentzVector p4CMS;
54  TLorentzVector p4Lab;
56  bool isTrack = false;
58  double clusterTime = NAN;
59 };
60 
62 const double flatBoundaries[10] = {0., 19., 22., 25., 30., 35., 45., 60., 90., 180.};
63 
64 
65 FilterCalculator::FilterCalculator() : m_bitsNN("CDCTriggerNNBits")
66 {
67 
68 }
69 
70 void FilterCalculator::requireStoreArrays()
71 {
72  m_tracks.isRequired();
73  m_eclClusters.isRequired();
74  m_l1Trigger.isOptional();
75  m_bitsNN.isOptional();
76 }
77 
78 void FilterCalculator::doCalculation(SoftwareTriggerObject& calculationResult)
79 {
80  calculationResult["nTrkLoose"] = 0;
81  calculationResult["nTrkTight"] = 0;
82  calculationResult["ee2leg"] = 0;
83  calculationResult["nEmedium"] = 0;
84  calculationResult["nElow"] = 0;
85  calculationResult["nEhigh"] = 0;
86  calculationResult["nE180Lab"] = 0;
87  calculationResult["nE300Lab"] = 0;
88  calculationResult["nE500Lab"] = 0;
89  calculationResult["nE2000CMS"] = 0;
90  calculationResult["nE4000CMS"] = 0;
91  calculationResult["nE250Lab"] = 0;
92  calculationResult["nMaxEPhotonAcc"] = 0;
93  calculationResult["dphiCmsClust"] = NAN;
95  calculationResult["netChargeLoose"] = 0;
96  calculationResult["maximumPCMS"] = NAN;
97  calculationResult["eexx"] = 0;
98  calculationResult["ee1leg1trk"] = 0;
99  calculationResult["nEhighLowAng"] = 0;
100  calculationResult["nEsingleClust"] = 0;
101  calculationResult["nEsinglePhotonBarrel"] = 0;
102  calculationResult["nEsinglePhotonExtendedBarrel"] = 0;
103  calculationResult["nEsinglePhotonEndcap"] = 0;
104  calculationResult["nEsingleElectronBarrel"] = 0;
105  calculationResult["nEsingleElectronExtendedBarrel"] = 0;
106  calculationResult["nReducedEsinglePhotonReducedBarrel"] = 0;
107  calculationResult["nVetoClust"] = 0;
108  calculationResult["chrgClust2GeV"] = 0;
109  calculationResult["neutClust045GeVAcc"] = 0;
110  calculationResult["neutClust045GeVBarrel"] = 0;
111  calculationResult["singleTagLowMass"] = 0;
112  calculationResult["singleTagHighMass"] = 0;
113  calculationResult["n2GeVNeutBarrel"] = 0;
114  calculationResult["n2GeVNeutEndcap"] = 0;
115  calculationResult["n2GeVChrg"] = 0;
116  calculationResult["n2GeVPhotonBarrel"] = 0;
117  calculationResult["n2GeVPhotonEndcap"] = 0;
118  calculationResult["ee1leg"] = 0;
119  calculationResult["ee1leg1clst"] = 0;
120  calculationResult["ee1leg1e"] = 0;
121  calculationResult["ee2clst"] = 0;
122  calculationResult["gg2clst"] = 0;
123  calculationResult["eeee"] = 0;
124  calculationResult["eemm"] = 0;
125  calculationResult["eexxSelect"] = 0;
126  calculationResult["radBhabha"] = 0;
127  calculationResult["eeBrem"] = 0;
128  calculationResult["isrRadBhabha"] = 0;
129  calculationResult["muonPairECL"] = 0;
130  calculationResult["selectee1leg1trk"] = 0;
131  calculationResult["selectee1leg1clst"] = 0;
132  calculationResult["selectee"] = 0;
133  calculationResult["ggBarrelVL"] = 0;
134  calculationResult["ggBarrelLoose"] = 0;
135  calculationResult["ggBarrelTight"] = 0;
136  calculationResult["ggEndcapVL"] = 0;
137  calculationResult["ggEndcapLoose"] = 0;
138  calculationResult["ggEndcapTight"] = 0;
139  calculationResult["muonPairV"] = 0;
140  calculationResult["selectmumu"] = 0;
141  calculationResult["singleMuon"] = 0;
142  calculationResult["cosmic"] = 0;
143  calculationResult["eeFlat0"] = 0;
144  calculationResult["eeFlat1"] = 0;
145  calculationResult["eeFlat2"] = 0;
146  calculationResult["eeFlat3"] = 0;
147  calculationResult["eeFlat4"] = 0;
148  calculationResult["eeFlat5"] = 0;
149  calculationResult["eeFlat6"] = 0;
150  calculationResult["eeFlat7"] = 0;
151  calculationResult["eeFlat8"] = 0;
152  calculationResult["eeOneClust"] = 0;
153  // Passed on L1 information
154  if (m_l1Trigger.isValid()) {
155  calculationResult["l1_trigger_random"] = m_l1Trigger->getTimType() == TRGSummary::ETimingType::TTYP_RAND;
156  calculationResult["l1_trigger_delayed_bhabha"] = m_l1Trigger->getTimType() == TRGSummary::ETimingType::TTYP_DPHY;
157  calculationResult["l1_trigger_poisson"] = m_l1Trigger->getTimType() == TRGSummary::ETimingType::TTYP_POIS;
158 
159  calculationResult["bha3d"] = m_l1Trigger->testPsnm("bha3d");
160  calculationResult["bhapur"] = m_l1Trigger->testPsnm("bhapur");
161  calculationResult["bhapur_lml1"] = m_l1Trigger->testPsnm("lml1") and m_l1Trigger->testFtdl("bhapur");
162  } else {
163  calculationResult["l1_trigger_random"] = 1; // save every event if no L1 trigger info
164  calculationResult["l1_trigger_delayed_bhabha"] = 0;
165  calculationResult["l1_trigger_poisson"] = 0;
166  calculationResult["bha3d"] = 0;
167  calculationResult["bhapur"] = 0;
168  calculationResult["bhapur_lml1"] = 0;
169  }
170 
171  // Every 256th event has CDC NN trigger information
172  calculationResult["l1_trg_NN_info"] = 0;
173  if (m_bitsNN.isValid() and m_bitsNN.getEntries() > 0) {calculationResult["l1_trg_NN_info"] = 1;}
174 
175  calculationResult["true"] = 1;
176  calculationResult["false"] = 0;
177 
178  // Some utilities
179  // TODO: into constructor
180  ClusterUtils cUtil;
181  const TVector3 clustervertex = cUtil.GetIPPosition();
182  PCmsLabTransform boostrotate;
183  TLorentzVector p4ofCOM;
184  p4ofCOM.SetPxPyPzE(0, 0, 0, boostrotate.getCMSEnergy());
185 
186  // Pointers to the two tracks with the maximum pt
187  std::map<short, std::optional<MaximumPtTrack>> maximumPtTracks = {
188  { -1, {}}, {1, {}}
189  };
190 
191  // Pointer to the two tracks with maximum pt without a cut applied on z0 (used for cosmic trigger)
192  std::map<short, std::optional<MaximumPtTrack>> maximumPtTracksWithoutZCut = {
193  { -1, {}}, {1, {}}
194  };
195 
196  // --- Track variables -- //
197  for (const Track& track : m_tracks) {
198  const TrackFitResult* trackFitResult = track.getTrackFitResultWithClosestMass(Const::pion);
199  if (not trackFitResult) {
200  // Hu? No track fit result for pion? Ok, lets skip this track
201  continue;
202  }
203 
204  // TODO: Temporary cut on number of CDC hits
205  if (trackFitResult->getHitPatternCDC().getNHits() == 0) {
206  continue;
207  }
208 
209  // Count loose and tight tracks
210  const double z0 = trackFitResult->getZ0();
211  if (std::abs(z0) < m_tightTrkZ0) {
212  calculationResult["nTrkTight"] += 1;
213  }
214 
215  // From here on use only tracks with defined charge
216  const short charge = trackFitResult->getChargeSign();
217  if (charge == 0) {
218  continue;
219  }
220 
221  const TLorentzVector& momentumLab = trackFitResult->get4Momentum();
222  const TLorentzVector momentumCMS = boostrotate.rotateLabToCms() * momentumLab;
223  double pCMS = momentumCMS.Rho();
224 
225  // Find the maximum pt negative [0] and positive [1] tracks without z0 cut
226  const double pT = trackFitResult->getTransverseMomentum();
227  const auto& currentMaximum = maximumPtTracksWithoutZCut.at(charge);
228  if (not currentMaximum or pT > currentMaximum->pT) {
229  MaximumPtTrack newMaximum;
230  newMaximum.pT = pT;
231  newMaximum.track = &track;
232  newMaximum.pCMS = pCMS;
233  newMaximum.pLab = momentumLab.Rho();
234  newMaximum.p4CMS = momentumCMS;
235  newMaximum.p4Lab = momentumLab;
236  maximumPtTracksWithoutZCut[charge] = newMaximum;
237  }
238 
239  // Loose tracks
240  if (std::abs(z0) < m_looseTrkZ0) {
241  calculationResult["nTrkLoose"] += 1;
242  calculationResult["netChargeLoose"] += charge;
243 
244  if (std::isnan(calculationResult["maximumPCMS"]) or pCMS > calculationResult["maximumPCMS"]) {
245  calculationResult["maximumPCMS"] = pCMS;
246  }
247 
248  // Find the maximum pt negative [0] and positive [1] tracks
249  const double pTLoose = trackFitResult->getTransverseMomentum();
250  const auto& currentMaximumLoose = maximumPtTracks.at(charge);
251  if (not currentMaximumLoose or pTLoose > currentMaximumLoose->pT) {
252  MaximumPtTrack newMaximum;
253  newMaximum.pT = pTLoose;
254  newMaximum.track = &track;
255  newMaximum.pCMS = pCMS;
256  newMaximum.pLab = momentumLab.Rho();
257  newMaximum.p4CMS = momentumCMS;
258  newMaximum.p4Lab = momentumLab;
259  maximumPtTracks[charge] = newMaximum;
260  }
261  }
262  }
263 
264  // -- Cluster variables -- //
265  std::vector<SelectedECLCluster> selectedClusters;
266  for (const ECLCluster& cluster : m_eclClusters) {
267  // Use the photon hypothesis, and only look at clusters with times < 200 ns
268  const double time = cluster.getTime();
269  if (not cluster.hasHypothesis(Belle2::ECLCluster::EHypothesisBit::c_nPhotons) or
270  std::abs(time) > 200 or
271  cluster.getEnergy(Belle2::ECLCluster::EHypothesisBit::c_nPhotons) < 0.1) {
272  continue;
273  }
274 
275  const double dt99 = cluster.getDeltaTime99();
276 
277  // Store infos if above the single photon threshold
278  SelectedECLCluster selectedCluster;
279  selectedCluster.p4Lab = cUtil.Get4MomentumFromCluster(&cluster, clustervertex,
280  Belle2::ECLCluster::EHypothesisBit::c_nPhotons);
281  selectedCluster.p4CMS = boostrotate.rotateLabToCms() * selectedCluster.p4Lab; // was clustp4COM
282  selectedCluster.cluster = &cluster;
283  selectedCluster.clusterTime = time / dt99; // was clustT
284  selectedCluster.energyLab = selectedCluster.p4Lab.E();
285  selectedCluster.energyCMS = selectedCluster.p4CMS.E(); // was first of ECOMPair
286  selectedCluster.isTrack = cluster.isTrack(); // was tempCharge
287 
288  selectedClusters.push_back(selectedCluster);
289 
290  if (selectedCluster.energyCMS > m_E2min) {
291  calculationResult["nElow"] += 1;
292  }
293  if (selectedCluster.energyCMS > m_E0min) {
294  calculationResult["nEmedium"] += 1;
295  }
296  if (selectedCluster.energyCMS > m_Ehigh) {
297  calculationResult["nEhigh"] += 1;
298  }
299 
300  if (selectedCluster.energyCMS > m_E0min and std::abs(selectedCluster.clusterTime) < 10) {
301  calculationResult["nVetoClust"] += 1;
302  }
303 
304 
305  //..Single cluster trigger objects use charge, Zernike moment, and thetaLab
306  const double thetaLab = selectedCluster.p4Lab.Theta() * TMath::RadToDeg();
307  const double zmva = cluster.getZernikeMVA();
308  const bool photon = zmva > 0.5 and not selectedCluster.isTrack;
309  const bool electron = zmva > 0.5 and selectedCluster.isTrack;
310 
311  //..For 1 track radiative Bhabha control sample
312  if (selectedCluster.energyCMS > 2. and selectedCluster.isTrack) {
313  calculationResult["chrgClust2GeV"] += 1;
314  }
315  if (selectedCluster.energyCMS > 0.45 and not selectedCluster.isTrack) {
316  const bool isInAcceptance = 17. < thetaLab and thetaLab < 150.;
317  if (isInAcceptance) {calculationResult["neutClust045GeVAcc"] += 1;}
318  const bool isInBarrel = 30. < thetaLab and thetaLab < 130.;
319  if (isInBarrel) {calculationResult["neutClust045GeVBarrel"] += 1;}
320  }
321 
322 
323  // improved 3 cluster (4 cluster) trigger
324  const bool notInHighBackgroundEndcapRegion = 18.5 < thetaLab and thetaLab < 139.3;
325  if (selectedCluster.energyLab > m_EminLab and notInHighBackgroundEndcapRegion) {
326  calculationResult["nE180Lab"] += 1;
327  }
328 
329  if (selectedCluster.energyLab > m_EminLab4Cluster and notInHighBackgroundEndcapRegion) {
330  calculationResult["nE300Lab"] += 1;
331  }
332 
333  if (selectedCluster.energyLab > m_EminLab3Cluster and notInHighBackgroundEndcapRegion) {
334  calculationResult["nE500Lab"] += 1;
335  }
336 
337  if (selectedCluster.energyCMS > m_Ehigh and notInHighBackgroundEndcapRegion) {
338  calculationResult["nE2000CMS"] += 1;
339  }
340 
341  //..For two-photon fusion ALP trigger
342  if (selectedCluster.energyLab > 0.25) {
343  calculationResult["nE250Lab"] += 1;
344  }
345  if (selectedCluster.energyCMS > 4.) {
346  calculationResult["nE4000CMS"] += 1;
347  }
348 
349  //..Single cluster triggers
350  if (selectedCluster.energyCMS > m_EsinglePhoton) {
351  calculationResult["nEsingleClust"] += 1;
352 
353  const bool barrelRegion = thetaLab > 45 and thetaLab < 115;
354  const bool extendedBarrelRegion = thetaLab > 30 and thetaLab < 130;
355  const bool endcapRegion = (thetaLab > 22 and thetaLab < 45) or (thetaLab > 115 and thetaLab < 145);
356 
357  if (photon and barrelRegion) {
358  calculationResult["nEsinglePhotonBarrel"] += 1;
359  }
360 
361  if (photon and extendedBarrelRegion) {
362  calculationResult["nEsinglePhotonExtendedBarrel"] += 1;
363  }
364 
365  if (electron and barrelRegion) {
366  calculationResult["nEsingleElectronBarrel"] += 1;
367  }
368 
369  if (electron and extendedBarrelRegion) {
370  calculationResult["nEsingleElectronExtendedBarrel"] += 1;
371  }
372 
373  if (photon and endcapRegion) {
374  calculationResult["nEsinglePhotonEndcap"] += 1;
375  }
376  }
377 
378  if (selectedCluster.energyCMS > m_reducedEsinglePhoton) {
379  const bool reducedBarrelRegion = thetaLab > 44 and thetaLab < 98;
380 
381  if (photon and reducedBarrelRegion) {
382  calculationResult["nReducedEsinglePhotonReducedBarrel"] += 1;
383  }
384  }
385 
386  if (selectedCluster.energyCMS > m_Ehigh) {
387  // TODO: different definition!
388  const bool barrelRegion = thetaLab > 32 and thetaLab < 130;
389  const bool endcapRegion = (thetaLab > 22 and thetaLab < 32) or (thetaLab > 130 and thetaLab < 145);
390  const bool lowAngleRegion = thetaLab < 22 or thetaLab > 145;
391 
392  if (not selectedCluster.isTrack and barrelRegion) {
393  calculationResult["n2GeVNeutBarrel"] += 1;
394  }
395  if (not selectedCluster.isTrack and endcapRegion) {
396  calculationResult["n2GeVNeutEndcap"] += 1;
397  }
398  if (selectedCluster.isTrack and not lowAngleRegion) {
399  calculationResult["n2GeVChrg"] += 1;
400  }
401  if (lowAngleRegion) {
402  calculationResult["nEhighLowAng"] += 1;
403  }
404  if (photon and barrelRegion) {
405  calculationResult["n2GeVPhotonBarrel"] += 1;
406  }
407  if (photon and endcapRegion) {
408  calculationResult["n2GeVPhotonEndcap"] += 1;
409  }
410  }
411  }
412  std::sort(selectedClusters.begin(), selectedClusters.end(), [](const auto & lhs, const auto & rhs) {
413  return lhs.energyCMS > rhs.energyCMS;
414  });
415 
416  // -- Bhabha and two-photon lepton preparation -- //
417  for (short charge : { -1, 1}) {
418  auto& maximumPtTrack = maximumPtTracks.at(charge);
419  if (not maximumPtTrack) {
420  continue;
421  }
422 
423  // Prep max two tracks for Bhabha and two-photon lepton selections
424  // Track / cluster matching
425  for (auto& cluster : maximumPtTrack->track->getRelationsTo<ECLCluster>()) {
426  if (cluster.hasHypothesis(Belle2::ECLCluster::EHypothesisBit::c_nPhotons)) {
427  maximumPtTrack->clusterEnergySumLab += cluster.getEnergy(Belle2::ECLCluster::EHypothesisBit::c_nPhotons);
428  }
429  }
430  maximumPtTrack->clusterEnergySumCMS =
431  maximumPtTrack->clusterEnergySumLab * maximumPtTrack->pCMS / maximumPtTrack->pLab;
432 
433  // Single leg Bhabha selections
434  if (maximumPtTrack->clusterEnergySumCMS > 4.5) {
435  calculationResult["ee1leg"] = 1;
436  }
437 
438  // single muon
439  if (maximumPtTrack->pCMS > 3 and maximumPtTrack->clusterEnergySumLab > 0. and maximumPtTrack->clusterEnergySumLab < 1.) {
440  calculationResult["singleMuon"] = 1;
441  }
442  }
443 
444  // Bhabha selections using two tracks
445  if (maximumPtTracks.at(-1) and maximumPtTracks.at(1)) {
446  const MaximumPtTrack& negativeTrack = *maximumPtTracks.at(-1);
447  const MaximumPtTrack& positiveTrack = *maximumPtTracks.at(1);
448 
449  double dphi = std::abs(negativeTrack.p4CMS.Phi() - positiveTrack.p4CMS.Phi()) * TMath::RadToDeg();
450  if (dphi > 180) {
451  dphi = 360 - dphi;
452  }
453 
454  const double negativeClusterSum = negativeTrack.clusterEnergySumCMS;
455  const double negativeClusterSumLab = negativeTrack.clusterEnergySumLab;
456  const double negativeP = negativeTrack.pCMS;
457  const double positiveClusterSum = positiveTrack.clusterEnergySumCMS;
458  const double positiveClusterSumLab = positiveTrack.clusterEnergySumLab;
459  const double positiveP = positiveTrack.pCMS;
460 
461  // two leg Bhabha
462  const double thetaSum = (negativeTrack.p4CMS.Theta() + positiveTrack.p4CMS.Theta()) * TMath::RadToDeg();
463  const double dthetaSum = std::abs(thetaSum - 180);
464  const auto back2back = dphi > 175 and dthetaSum < 15;
465  if (back2back and negativeClusterSum > 3 and positiveClusterSum > 3 and
466  (negativeClusterSum > 4.5 or positiveClusterSum > 4.5)) {
467  calculationResult["ee2leg"] = 1;
468  }
469 
470  // one leg one track Bhabha
471  if (back2back and ((negativeClusterSum > 4.5 and positiveP > 3) or (positiveClusterSum > 4.5 and negativeP > 3))) {
472  calculationResult["ee1leg1trk"] = 1;
473  }
474 
475 
476  // one leg plus one electron
477  if ((negativeClusterSum > 4.5 and positiveClusterSum > 0.8 * positiveP) or
478  (positiveClusterSum > 4.5 and negativeClusterSum > 0.8 * negativeP)) {
479  calculationResult["ee1leg1e"] = 1;
480  }
481 
482  // eeee, eemm, mu mu, and radiative Bhabha selection
483  const TLorentzVector p4Miss = p4ofCOM - negativeTrack.p4CMS - positiveTrack.p4CMS;
484  const double pmissTheta = p4Miss.Theta() * TMath::RadToDeg();
485  const double pmissp = p4Miss.Rho();
486  const double relMissAngle0 = negativeTrack.p4CMS.Angle(p4Miss.Vect()) * TMath::RadToDeg();
487  const double relMissAngle1 = positiveTrack.p4CMS.Angle(p4Miss.Vect()) * TMath::RadToDeg();
488 
489  const bool electronEP = positiveClusterSum > 0.8 * positiveP or negativeClusterSum > 0.8 * negativeP;
490  const bool notMuonPair = negativeClusterSumLab > 1 or positiveClusterSumLab > 1;
491  const double highp = std::max(negativeP, positiveP);
492  const double lowp = std::min(negativeP, positiveP);
493  const bool lowEdep = negativeClusterSumLab < 0.5 and positiveClusterSumLab < 0.5;
494 
495 
496  if (calculationResult["maximumPCMS"] < 2 and dphi > 160 and (pmissTheta < 25. or pmissTheta > 155.)) {
497  calculationResult["eexxSelect"] = 1;
498  if (electronEP) {
499  calculationResult["eeee"] = 1;
500  } else {
501  calculationResult["eemm"] = 1;
502  }
503  }
504  if (1. < negativeP and 1. < positiveP and 2 == calculationResult["nTrkLoose"] and
505  calculationResult["nTrkTight"] >= 1 and dphi < 170. and
506  10. < pmissTheta and pmissTheta < 170. and 1. < pmissp and electronEP) {
507  calculationResult["radBhabha"] = 1;
508  }
509 
510  if (2. < negativeP and 2. < positiveP and 2 == calculationResult["nTrkLoose"] and
511  calculationResult["nTrkTight"] >= 1 and 175. < dphi and
512  (pmissTheta < 5. or pmissTheta > 175.) and electronEP) {
513  calculationResult["isrRadBhabha"] = 1;
514  }
515  if ((pmissTheta < 20. or pmissTheta > 160.) and
516  ((calculationResult["maximumPCMS"] < 1.2 and dphi > 150.) or
517  (calculationResult["maximumPCMS"] < 2. and 175. < dphi))) {
518  calculationResult["eexx"] = 1;
519  }
520  if (calculationResult["nTrkLoose"] == 2 and highp > 4.5 and notMuonPair and pmissp > 1. and
521  (relMissAngle0 < 5. or relMissAngle1 < 5.)) {
522  calculationResult["eeBrem"] = 1;
523  }
524 
525  if (calculationResult["nTrkLoose"] == 2 and highp > 4.5 and lowEdep and dphi > 175. and 175. < thetaSum and
526  thetaSum < 185.) {
527  calculationResult["muonPairV"] = 1;
528  }
529  if (3. < highp and 2.5 < lowp and 165. < dphi and
530  ((negativeClusterSumLab > 0. and negativeClusterSumLab < 1.) or
531  (positiveClusterSumLab > 0. and positiveClusterSumLab < 1.))) {
532  calculationResult["selectmumu"] = 1;
533  }
534  }
535 
536  if (selectedClusters.size() >= 2) {
537  const SelectedECLCluster& firstCluster = selectedClusters[0];
538  const SelectedECLCluster& secondCluster = selectedClusters[1];
539 
540  // Bhabha / gg vetoes using max two clusters
541  double dphi = std::abs(firstCluster.p4CMS.Phi() - secondCluster.p4CMS.Phi()) * TMath::RadToDeg();
542  if (dphi > 180) {
543  dphi = 360 - dphi;
544  }
545  double thetaSum = (firstCluster.p4CMS.Theta() + secondCluster.p4CMS.Theta()) * TMath::RadToDeg();
546  double dthetaSum = std::abs(thetaSum - 180);
547 
548  //..Quantities for two-photon fusion triggers
549  calculationResult["dphiCmsClust"] = dphi;
550  for (int ic = 0; ic < 2; ic++) {
551  const double thetaLab = selectedClusters[ic].p4Lab.Theta() * TMath::RadToDeg();
552  const bool isInAcceptance = 17. < thetaLab and thetaLab < 150.;
553  const ECLCluster* cluster = selectedClusters[ic].cluster;
554  const double zmva = cluster->getZernikeMVA();
555  const bool photon = zmva > 0.5 and not selectedClusters[ic].isTrack;
556  if (isInAcceptance and photon) {calculationResult["nMaxEPhotonAcc"] += 1;}
557  }
558 
559  const double firstEnergy = firstCluster.p4CMS.E();
560  const double secondEnergy = secondCluster.p4CMS.E();
561 
562  const bool highEnergetic = firstEnergy > 3 and secondEnergy > 3 and (firstEnergy > 4.5 or secondEnergy > 4.5);
563 
564  if (dphi > 160 and dphi < 178 and dthetaSum < 15 and highEnergetic) {
565  calculationResult["ee2clst"] = 1;
566  }
567 
568  if (dphi > 178 and dthetaSum < 15 and highEnergetic) {
569  calculationResult["gg2clst"] = 1;
570  }
571 
572  if ((calculationResult["ee2clst"] == 1 or calculationResult["gg2clst"] == 1) and
573  calculationResult["ee1leg"] == 1) {
574  calculationResult["ee1leg1clst"] = 1;
575  }
576 
577  const double Elab0 = firstCluster.p4Lab.E();
578  const double Elab1 = secondCluster.p4Lab.E();
579 
580  // gg and mumu accept triggers using ECL
581  if (firstEnergy > 2 and secondEnergy > 2) {
582  const double thetaLab0 = firstCluster.p4Lab.Theta() * TMath::RadToDeg();
583  const double thetaLab1 = secondCluster.p4Lab.Theta() * TMath::RadToDeg();
584 
585  const bool barrel0 = 32. < thetaLab0 and thetaLab0 < 130.;
586  const bool barrel1 = 32. < thetaLab1 and thetaLab1 < 130.;
587  const bool oneClustersAbove4 = firstEnergy > 4 or secondEnergy > 4;
588  const bool oneIsNeutral = not firstCluster.isTrack or not secondCluster.isTrack;
589  const bool bothAreNeutral = not firstCluster.isTrack and not secondCluster.isTrack;
590  const bool oneIsBarrel = barrel0 or barrel1;
591  const bool dphiCutExtraLoose = dphi > 175;
592  const bool dphiCutLoose = dphi > 177;
593  const bool dphiCutTight = dphi > 177.5;
594 
595  if (dphiCutExtraLoose and oneIsNeutral and oneIsBarrel) {
596  calculationResult["ggBarrelVL"] = 1;
597  }
598  if (oneClustersAbove4 and dphiCutLoose and oneIsNeutral and oneIsBarrel) {
599  calculationResult["ggBarrelLoose"] = 1;
600  }
601  if (oneClustersAbove4 and dphiCutTight and bothAreNeutral and oneIsBarrel) {
602  calculationResult["ggBarrelTight"] = 1;
603  }
604  if (dphiCutExtraLoose and oneIsNeutral and not oneIsBarrel) {
605  calculationResult["ggEndcapVL"] = 1;
606  }
607  if (oneClustersAbove4 and dphiCutLoose and oneIsNeutral and not oneIsBarrel) {
608  calculationResult["ggEndcapLoose"] = 1;
609  }
610  if (oneClustersAbove4 and dphiCutTight and bothAreNeutral and not oneIsBarrel) {
611  calculationResult["ggEndcapTight"] = 1;
612  }
613  }
614 
615  const double minEnergy = std::min(Elab0, Elab1);
616  const double maxEnergy = std::max(Elab0, Elab1);
617  if (dphi > 155 and thetaSum > 165 and thetaSum < 195 and minEnergy > 0.15 and minEnergy < 0.5 and
618  maxEnergy > 0.15 and maxEnergy < 0.5) {
619  calculationResult["muonPairECL"] = 1;
620  }
621 
622  }
623 
624  // Bhabha accept triggers.
625  // Use theta_lab of negative track if available; otherwise cluster.
626  double thetaFlatten = 0;
627 
628  // One leg, one cluster.
629  for (short charge : { -1, 1}) {
630  const auto& maximumPtTrack = maximumPtTracks.at(charge);
631  if (not maximumPtTrack) {
632  continue;
633  }
634 
635  if (maximumPtTrack->clusterEnergySumCMS > 1.5) {
636  double invMass = 0.;
637  double tempFlatten = 0.;
638  for (const SelectedECLCluster& cluster : selectedClusters) {
639  double tempInvMass = (maximumPtTrack->p4Lab + cluster.p4Lab).M();
640  if (tempInvMass > invMass) {
641  invMass = tempInvMass;
642  if (charge == 1) {
643  tempFlatten = cluster.p4Lab.Theta() * TMath::RadToDeg();
644  }
645  }
646  }
647  if (charge == -1) {
648  tempFlatten = maximumPtTrack->p4Lab.Theta() * TMath::RadToDeg();
649  }
650  if (invMass > 5.29) {
651  calculationResult["selectee1leg1clst"] = 1;
652  thetaFlatten = tempFlatten;
653  }
654  }
655  }
656 
657  // Two legs
658  if (maximumPtTracks.at(-1) and maximumPtTracks.at(1)) {
659  const MaximumPtTrack& negativeTrack = *maximumPtTracks.at(-1);
660  const MaximumPtTrack& positiveTrack = *maximumPtTracks.at(1);
661  const double invMass = (negativeTrack.p4Lab + positiveTrack.p4Lab).M();
662  if (invMass > 5.29 and (negativeTrack.clusterEnergySumCMS > 1.5 or positiveTrack.clusterEnergySumCMS > 1.5)) {
663  calculationResult["selectee1leg1trk"] = 1;
664  }
665 
666  thetaFlatten = negativeTrack.p4Lab.Theta() * TMath::RadToDeg();
667 
668  // Two tracks but (exactly) 1 high energy cluster
669  const bool missNegClust = positiveTrack.clusterEnergySumCMS > 4.5 and negativeTrack.clusterEnergySumCMS < 1.;
670  const bool missPosClust = negativeTrack.clusterEnergySumCMS > 4.5 and positiveTrack.clusterEnergySumCMS < 1.;
671  if ((invMass > 9.) and (missNegClust or missPosClust)) {
672  calculationResult["eeOneClust"] = 1;
673  }
674  }
675 
676  if (calculationResult["selectee1leg1trk"] == 1 or calculationResult["selectee1leg1clst"] == 1) {
677  for (int iflat = 0; iflat < 9; iflat++) {
678  const std::string& eeFlatName = "eeFlat" + std::to_string(iflat);
679  calculationResult[eeFlatName] =
680  thetaFlatten > flatBoundaries[iflat] and thetaFlatten < flatBoundaries[iflat + 1];
681  if (calculationResult[eeFlatName]) {
682  calculationResult["selectee"] = 1;
683  }
684  }
685  }
686 
687  // Single tag pi0 / eta dedicated lines
688  if (calculationResult["nTrkLoose"] == 1 and calculationResult["maximumPCMS"] > 0.8 and selectedClusters.size() >= 2) {
689 
690  decltype(selectedClusters) selectedSingleTagClusters(selectedClusters.size());
691  auto lastItem = std::copy_if(selectedClusters.begin(), selectedClusters.end(), selectedSingleTagClusters.begin(),
692  [](auto & cluster) {
693  const bool isNeutralCluster = not cluster.isTrack;
694  const bool hasEnoughEnergy = cluster.energyLab > 0.1;
695  const double clusterThetaLab = cluster.p4Lab.Theta() * TMath::RadToDeg();
696  const bool isInAcceptance = 17 < clusterThetaLab and clusterThetaLab < 150.;
697  return isNeutralCluster and hasEnoughEnergy and isInAcceptance;
698  });
699  selectedSingleTagClusters.resize(std::distance(selectedSingleTagClusters.begin(), lastItem));
700 
701  if (selectedSingleTagClusters.size() >= 2) { // One track and at least two clusters are found
702 
703  const auto& track = maximumPtTracks.at(-1) ? *maximumPtTracks.at(-1) : *maximumPtTracks.at(1);
704  // in real signal, the pi0 decay daughters are always the two most-energetic neutral clusters.
705  const auto& firstCluster = selectedSingleTagClusters[0];
706  const auto& secondCluster = selectedSingleTagClusters[1];
707 
708  const TLorentzVector trackP4CMS = track.p4CMS;
709  const TLorentzVector pi0P4CMS = firstCluster.p4CMS + secondCluster.p4CMS;
710 
711  const bool passPi0ECMS = pi0P4CMS.E() > 1. and pi0P4CMS.E() < 0.525 * p4ofCOM.M();
712  const double thetaSumCMS = (pi0P4CMS.Theta() + trackP4CMS.Theta()) * TMath::RadToDeg();
713  const bool passThetaSum = thetaSumCMS < 170. or thetaSumCMS > 190.;
714 
715  double dphiCMS = std::abs(trackP4CMS.Phi() - pi0P4CMS.Phi()) * TMath::RadToDeg();
716  if (dphiCMS > 180) {
717  dphiCMS = 360 - dphiCMS;
718  }
719  const bool passdPhi = dphiCMS > 160.;
720 
721  if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() < 0.7) {
722  calculationResult["singleTagLowMass"] = 1;
723  } else if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() > 0.7) {
724  calculationResult["singleTagHighMass"] = 1;
725  }
726  }
727  }
728 
729  //..Cosmic selection. Need exactly two tracks.
730  if (m_tracks.getEntries() == 2) {
731 
732  const auto negTrack = maximumPtTracksWithoutZCut.at(-1);
733  const auto posTrack = maximumPtTracksWithoutZCut.at(1);
734 
735  if (negTrack and posTrack) {
736 
737  const double maxNegpT = negTrack->pT;
738  const double maxPospT = posTrack->pT;
739 
740  auto accumulatePhotonEnergy = [](double result, const auto & cluster) {
741  return result + (cluster.hasHypothesis(Belle2::ECLCluster::EHypothesisBit::c_nPhotons) ? cluster.getEnergy(
742  Belle2::ECLCluster::EHypothesisBit::c_nPhotons) : 0);
743  };
744 
745  const auto& clustersOfNegTrack = negTrack->track->getRelationsTo<ECLCluster>();
746  const auto& clustersOfPosTrack = posTrack->track->getRelationsTo<ECLCluster>();
747 
748  const double maxClusterENeg = std::accumulate(clustersOfNegTrack.begin(), clustersOfNegTrack.end(), 0.0, accumulatePhotonEnergy);
749  const double maxClusterEPos = std::accumulate(clustersOfPosTrack.begin(), clustersOfPosTrack.end(), 0.0, accumulatePhotonEnergy);
750 
751  const TLorentzVector& momentumLabNeg(negTrack->p4Lab);
752  const TLorentzVector& momentumLabPos(posTrack->p4Lab);
753 
754  const double& z0Neg = negTrack->track->getTrackFitResultWithClosestMass(Const::pion)->getZ0();
755  const double& d0Neg = negTrack->track->getTrackFitResultWithClosestMass(Const::pion)->getD0();
756  const double& z0Pos = posTrack->track->getTrackFitResultWithClosestMass(Const::pion)->getZ0();
757  const double& d0Pos = posTrack->track->getTrackFitResultWithClosestMass(Const::pion)->getD0();
758 
759  // Select cosmic using these tracks
760  const bool goodMagneticRegion = (z0Neg<m_goodMagneticRegionZ0 or abs(d0Neg)>m_goodMagneticRegionD0)
761  and (z0Pos<m_goodMagneticRegionZ0 or abs(d0Pos)>m_goodMagneticRegionD0);
762  if (maxNegpT > m_cosmicMinPt and maxPospT > m_cosmicMinPt and maxClusterENeg < m_cosmicMaxClusterEnergy
763  and maxClusterEPos < m_cosmicMaxClusterEnergy and goodMagneticRegion) {
764  double dphiLab = std::abs(momentumLabNeg.Phi() - momentumLabPos.Phi()) * TMath::RadToDeg();
765  if (dphiLab > 180) {
766  dphiLab = 360 - dphiLab;
767  }
768 
769  const double thetaSumLab = (momentumLabNeg.Theta() + momentumLabPos.Theta()) * TMath::RadToDeg();
770 
771  constexpr double phiBackToBackTolerance = 2.;
772  constexpr double thetaBackToBackTolerance = 2.;
773  if ((180 - dphiLab) < phiBackToBackTolerance and std::abs(180 - thetaSumLab) < thetaBackToBackTolerance) {
774  calculationResult["cosmic"] = 1;
775  }
776  }
777  }
778  }
779 
780 }
MaximumPtTrack::track
const Track * track
the track
Definition: FilterCalculator.cc:28
MaximumPtTrack::pCMS
double pCMS
the momentum magnitude in CMS system
Definition: FilterCalculator.cc:34
Belle2::SoftwareTrigger::FilterCalculator::m_goodMagneticRegionD0
double m_goodMagneticRegionD0
minimum d0 for well understood magnetic field, if z0 is large (cm)
Definition: FilterCalculator.h:96
MaximumPtTrack::p4Lab
TLorentzVector p4Lab
the 4 momentum in lab system
Definition: FilterCalculator.cc:40
MaximumPtTrack::pT
double pT
the pT of the track
Definition: FilterCalculator.cc:26
Belle2::ECLCluster
ECL cluster data.
Definition: ECLCluster.h:39
Belle2::SoftwareTrigger::FilterCalculator::requireStoreArrays
void requireStoreArrays() override
Require the particle list. We do not need more here.
Definition: FilterCalculator.cc:70
Belle2::SoftwareTrigger::FilterCalculator::m_tightTrkZ0
double m_tightTrkZ0
which Z0 defines a tight track
Definition: FilterCalculator.h:72
MaximumPtTrack::p4CMS
TLorentzVector p4CMS
the 4 momentum in CMS system
Definition: FilterCalculator.cc:38
SelectedECLCluster::cluster
const ECLCluster * cluster
The ECL cluster.
Definition: FilterCalculator.cc:46
Belle2::PCmsLabTransform::getCMSEnergy
double getCMSEnergy() const
Returns CMS energy of e+e- (aka.
Definition: PCmsLabTransform.h:57
Belle2::ClusterUtils::Get4MomentumFromCluster
const TLorentzVector Get4MomentumFromCluster(const ECLCluster *cluster, ECLCluster::EHypothesisBit hypo)
Returns four momentum vector.
Definition: ClusterUtils.cc:26
MaximumPtTrack::clusterEnergySumLab
double clusterEnergySumLab
the sum of related cluster energies in lab system
Definition: FilterCalculator.cc:32
SelectedECLCluster
Temporary data structure holding the ECL clusters used for this analysis.
Definition: FilterCalculator.cc:44
Belle2::ECLCluster::EHypothesisBit::c_nPhotons
@ c_nPhotons
CR is split into n photons (N1)
Belle2::SoftwareTrigger::FilterCalculator::FilterCalculator
FilterCalculator()
Set the default names for the store object particle lists.
Definition: FilterCalculator.cc:65
MaximumPtTrack
Temporary data structure holding the track(s) with the maximum pT.
Definition: FilterCalculator.cc:24
Belle2::TrackFitResult
Values of the result of a track fit with a given particle hypothesis.
Definition: TrackFitResult.h:59
Belle2::SoftwareTrigger::FilterCalculator::m_EminLab4Cluster
double m_EminLab4Cluster
which lab energy defines nE300Lab
Definition: FilterCalculator.h:82
Belle2::RelationsInterface::getRelationsTo
RelationVector< TO > getRelationsTo(const std::string &name="", const std::string &namedRelation="") const
Get the relations that point from this object to another store array.
Definition: RelationsObject.h:199
Belle2::ClusterUtils::GetIPPosition
const TVector3 GetIPPosition()
Returns default IP position from beam parameters.
Definition: ClusterUtils.cc:182
Belle2::SoftwareTrigger::FilterCalculator::m_cosmicMaxClusterEnergy
double m_cosmicMaxClusterEnergy
which LAB cluster energy vetoes a cosmic candidate
Definition: FilterCalculator.h:92
SelectedECLCluster::clusterTime
double clusterTime
the time of the cluster
Definition: FilterCalculator.cc:58
Belle2::SoftwareTrigger::FilterCalculator::m_looseTrkZ0
double m_looseTrkZ0
which Z0 defines a loose track
Definition: FilterCalculator.h:70
Belle2::SoftwareTrigger::FilterCalculator::m_E0min
double m_E0min
which CMS energy defines nEmedium
Definition: FilterCalculator.h:76
SelectedECLCluster::p4Lab
TLorentzVector p4Lab
the 4 momentum in lab system
Definition: FilterCalculator.cc:54
Belle2::Const::pion
static const ChargedStable pion
charged pion particle
Definition: Const.h:535
Belle2::ClusterUtils
Class to provide momentum-related information from ECLClusters.
Definition: ClusterUtils.h:44
SelectedECLCluster::p4CMS
TLorentzVector p4CMS
the 4 momentum in CMS system
Definition: FilterCalculator.cc:52
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::SoftwareTrigger::FilterCalculator::m_EsinglePhoton
double m_EsinglePhoton
which CMS energy defines nEsingleClust
Definition: FilterCalculator.h:86
Belle2::SoftwareTrigger::FilterCalculator::m_l1Trigger
StoreObjPtr< TRGSummary > m_l1Trigger
Store Object with the trigger result.
Definition: FilterCalculator.h:65
Belle2::SoftwareTrigger::FilterCalculator::m_Ehigh
double m_Ehigh
which CMS energy defines nEhigh
Definition: FilterCalculator.h:78
SelectedECLCluster::isTrack
bool isTrack
is this ECL cluster likely from a track (or a photon) = is it charged?
Definition: FilterCalculator.cc:56
MaximumPtTrack::clusterEnergySumCMS
double clusterEnergySumCMS
the sum of related cluster energies in CMS system
Definition: FilterCalculator.cc:30
SelectedECLCluster::energyLab
double energyLab
the energy in Lab system
Definition: FilterCalculator.cc:50
Belle2::SoftwareTrigger::FilterCalculator::doCalculation
void doCalculation(SoftwareTriggerObject &calculationResult) override
Actually write out the variables into the map.
Definition: FilterCalculator.cc:78
Belle2::PCmsLabTransform
Class to hold Lorentz transformations from/to CMS and boost vector.
Definition: PCmsLabTransform.h:37
Belle2::PCmsLabTransform::rotateLabToCms
const TLorentzRotation rotateLabToCms() const
Returns Lorentz transformation from Lab to CMS.
Definition: PCmsLabTransform.h:74
Belle2::SoftwareTrigger::FilterCalculator::m_EminLab3Cluster
double m_EminLab3Cluster
which lab energy defines nE500Lab
Definition: FilterCalculator.h:84
Belle2::SoftwareTrigger::FilterCalculator::m_bitsNN
StoreArray< CDCTriggerUnpacker::NNBitStream > m_bitsNN
Store Object with the trigger NN bits.
Definition: FilterCalculator.h:67
Belle2::Track
Class that bundles various TrackFitResults.
Definition: Track.h:35
Belle2::SoftwareTrigger::FilterCalculator::m_reducedEsinglePhoton
double m_reducedEsinglePhoton
which CMS energy defines nReducedEsingle clusters
Definition: FilterCalculator.h:88
SelectedECLCluster::energyCMS
double energyCMS
the energy in CMS system
Definition: FilterCalculator.cc:48
MaximumPtTrack::pLab
double pLab
the momentum magnitude in lab system
Definition: FilterCalculator.cc:36
Belle2::SoftwareTrigger::FilterCalculator::m_cosmicMinPt
double m_cosmicMinPt
which LAB pt defines a cosmic
Definition: FilterCalculator.h:90
Belle2::SoftwareTrigger::FilterCalculator::m_eclClusters
StoreArray< ECLCluster > m_eclClusters
Store Array of the ecl clusters to be used.
Definition: FilterCalculator.h:63
Belle2::SoftwareTrigger::FilterCalculator::m_E2min
double m_E2min
which CMS energy defines nElow
Definition: FilterCalculator.h:74
Belle2::SoftwareTrigger::FilterCalculator::m_EminLab
double m_EminLab
which lab energy defines nE180Lab
Definition: FilterCalculator.h:80
Belle2::SoftwareTrigger::FilterCalculator::m_tracks
StoreArray< Track > m_tracks
Store Array of the tracks to be used.
Definition: FilterCalculator.h:61
Belle2::SoftwareTrigger::FilterCalculator::m_goodMagneticRegionZ0
double m_goodMagneticRegionZ0
maximum z0 for well understood magnetic field (cm)
Definition: FilterCalculator.h:94