Belle II Software  release-08-01-10
Public Member Functions | Private Attributes | List of all members
FilterCalculator Class Reference

Implementation of a calculator used in the SoftwareTriggerModule to fill a SoftwareTriggerObject for doing HLT cuts. More...

#include <FilterCalculator.h>

Inheritance diagram for FilterCalculator:
Inheritance graph
[legend]
Collaboration diagram for FilterCalculator:
Collaboration graph
[legend]

Public Member Functions

 FilterCalculator ()
 Set the default names for the store object particle lists.
 
void requireStoreArrays () override
 Require the particle list. We do not need more here.
 
void doCalculation (SoftwareTriggerObject &calculationResult) override
 Actually write out the variables into the map. More...
 
void writeDebugOutput (const std::unique_ptr< TTree > &debugOutputTTree)
 Function to write out debug output into the given TTree. More...
 
void addDebugOutput (const StoreObjPtr< SoftwareTriggerVariables > &storeObject, const std::string &prefix)
 Function to write out debug output into the given StoreObject. More...
 
const SoftwareTriggerObject & fillInCalculations ()
 Main function of this class: calculate the needed variables using the overwritten doCalculation function and write out the values into the results object (with their names). More...
 

Private Attributes

StoreArray< Track > m_tracks
 Store Array of the tracks to be used.
 
StoreArray< ECLCluster > m_eclClusters
 Store Array of the ecl clusters to be used.
 
StoreObjPtr< TRGSummary > m_l1Trigger
 Store Object with the trigger result.
 
StoreArray< CDCTriggerUnpacker::NNBitStream > m_bitsNN
 Store Object with the trigger NN bits.
 
double m_looseTrkZ0 = 10 * Unit::cm
 which Z0 defines a loose track
 
double m_tightTrkZ0 = 2 * Unit::cm
 which Z0 defines a tight track
 
double m_E2min = 0.2
 which CMS energy defines nElow
 
double m_E0min = 0.3
 which CMS energy defines nEmedium
 
double m_Ehigh = 2
 which CMS energy defines nEhigh
 
double m_EminLab = 0.18
 which lab energy defines nE180Lab
 
double m_EminLab4Cluster = 0.3
 which lab energy defines nE300Lab
 
double m_EminLab3Cluster = 0.5
 which lab energy defines nE500Lab
 
double m_EsinglePhoton = 1
 which CMS energy defines nEsingleClust
 
double m_reducedEsinglePhoton = 0.5
 which CMS energy defines nReducedEsingle clusters
 
double m_cosmicMinPt = 0.5 * Unit::GeV
 which LAB pt defines a cosmic
 
double m_cosmicMaxClusterEnergy = 1.0 * Unit::GeV
 which LAB cluster energy vetoes a cosmic candidate
 
double m_goodMagneticRegionZ0 = 57.
 maximum z0 for well understood magnetic field (cm)
 
double m_goodMagneticRegionD0 = 26.5
 minimum d0 for well understood magnetic field, if z0 is large (cm)
 
SoftwareTriggerObject m_calculationResult
 Internal storage of the result of the calculation.
 
bool m_debugPrepared = false
 Flag to not add the branches twice to the TTree.
 

Detailed Description

Implementation of a calculator used in the SoftwareTriggerModule to fill a SoftwareTriggerObject for doing HLT cuts.

This calculator exports variables needed for the trigger HLT part of the path ( = filtering out events)

This class implements the two main functions requireStoreArrays and doCalculation of the SoftwareTriggerCalculation class.

Definition at line 30 of file FilterCalculator.h.

Member Function Documentation

◆ addDebugOutput()

void addDebugOutput ( const StoreObjPtr< SoftwareTriggerVariables > &  storeObject,
const std::string &  prefix 
)
inherited

Function to write out debug output into the given StoreObject.

Needs an already prefilled calculationResult for this (probably using the fillInCalculations function). All added variables are prefixed with the given prefix string.

Definition at line 34 of file SoftwareTriggerCalculation.cc.

35  {
36  for (auto& identifierWithValue : m_calculationResult) {
37  const std::string& identifier = identifierWithValue.first;
38  const double value = identifierWithValue.second;
39 
40  storeObject->append(prefix + "_" + identifier, value);
41  }
42  }
Belle2::SoftwareTrigger::SoftwareTriggerCalculation::m_calculationResult
SoftwareTriggerObject m_calculationResult
Internal storage of the result of the calculation.
Definition: SoftwareTriggerCalculation.h:74

◆ doCalculation()

void doCalculation ( SoftwareTriggerObject &  calculationResult)
overridevirtual

Actually write out the variables into the map.

< number of loose tracks

< number of tight tracks

< Bhabha, 2 tracks with accompanying ECL info

< number of clusters with E*>m_Emedium (higher threshold)

< number of clusters with E*>m_Elow (lower threshold)

< number of clusters with E*>m_Ehigh (2 GeV)

< 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

< number of clusters with E*>4 GeV

< neutral clusters with Elab>250 MeV (anywhere)

< Neutral, zmva>0.5, in [17,150], max 2 energy clusters

< dphi* between 2 max E clusters

< net charge of loose tracks

< maximum p* of loose tracks (GeV/c)

< maximum pLab of loose tracks (GeV/c)

< Bhabha, 2 tracks, at least 1 of which has ECL info

< number of clusters with E*>m_Ehigh, low angles

< clusters with E*> m_EsinglePhoton (1 GeV)

< neutral clusters with E*> 1 GeV in [45,115]

< neutral clusters with E*> 1 GeV in [32,130]

< neutral clusters with E*> 1 GeV, not barrel or low

< charged clusters with E*> 1 GeV in [45,115]

< charged clusters with E*> 1 GeV in [32,130]

< charged clusters with E*> 0.5 GeV in [44,98]

< clusters with E>m_Emedium and |t|/dt99 < 10

< charged clusters with E*>2 GeV

< neutral clusters with E*>0.45 GeV in [17,150]

< neutral clusters with E*>0.45 GeV in [30,130]

< track + ECL info consistent with Bhabha

< Bhabha, 2 ECL clusters, one of which is charged

< Bhabha, 1 track has high energy cluster, other is electron

< Bhabha, 2 ECL clusters

< gg veto, 2 ECL clusters

< number of looseB tracks

< number of tightB tracks

< maximum pcms of looseB tracks (GeV/c)

< net charge of looseB tracks

< muon pair veto counting looseB tracks

< Bhabha veto hard Brem one leg

< single tag pi0 or eta

< single tag higher mass

< radiative Bhabha, recoil p in detector

Implements SoftwareTriggerCalculation.

Definition at line 84 of file FilterCalculator.cc.

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  //..New track definitions
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  // Passed on L1 information
175  if (m_l1Trigger.isValid()) {
176  calculationResult["l1_trigger_random"] = m_l1Trigger->getTimType() == TRGSummary::ETimingType::TTYP_RAND;
177  calculationResult["l1_trigger_delayed_bhabha"] = m_l1Trigger->getTimType() == TRGSummary::ETimingType::TTYP_DPHY;
178  calculationResult["l1_trigger_poisson"] = m_l1Trigger->getTimType() == TRGSummary::ETimingType::TTYP_POIS;
179  bool bha3d;
180  try {
181  bha3d = m_l1Trigger->testPsnm("bha3d");
182  } catch (const std::exception&) {
183  bha3d = false;
184  }
185  bool bhapurPsnm;
186  try {
187  bhapurPsnm = m_l1Trigger->testPsnm("bhapur");
188  } catch (const std::exception&) {
189  bhapurPsnm = false;
190  }
191  bool bhapurFtdl;
192  try {
193  bhapurFtdl = m_l1Trigger->testFtdl("bhapur");
194  } catch (const std::exception&) {
195  bhapurFtdl = false;
196  }
197  bool lml1;
198  try {
199  lml1 = m_l1Trigger->testPsnm("lml1");
200  } catch (const std::exception&) {
201  lml1 = false;
202  }
203  bool l1_bit_f;
204  try {
205  l1_bit_f = m_l1Trigger->testPsnm("fpre");
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; // save every event if no L1 trigger info
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  // Every 256th event has CDC NN trigger information
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  // Some utilities
231  // TODO: into constructor
232  ClusterUtils cUtil;
233  const ROOT::Math::XYZVector clustervertex = cUtil.GetIPPosition();
234  PCmsLabTransform boostrotate;
235  ROOT::Math::PxPyPzEVector p4ofCOM;
236  p4ofCOM.SetPxPyPzE(0, 0, 0, boostrotate.getCMSEnergy());
237 
238  // Pointers to the two tracks with the maximum pt
239  std::map<short, std::optional<MaximumPtTrack>> maximumPtTracks = {
240  { -1, {}}, {1, {}}
241  };
242 
243  // Pointer to the two tracks with maximum pt without a cut applied on z0 (used for cosmic trigger)
244  std::map<short, std::optional<MaximumPtTrack>> maximumPtTracksWithoutZCut = {
245  { -1, {}}, {1, {}}
246  };
247 
248  // --- Track variables -- //
249  for (const Track& track : m_tracks) {
250  const TrackFitResult* trackFitResult = track.getTrackFitResultWithClosestMass(Const::pion);
251  if (not trackFitResult) {
252  // Hu? No track fit result for pion? Ok, lets skip this track
253  continue;
254  }
255 
256  // TODO: Temporary cut on number of CDC hits
257  if (trackFitResult->getHitPatternCDC().getNHits() == 0) {
258  continue;
259  }
260 
261  // Count loose and tight tracks
262  const double z0 = trackFitResult->getZ0();
263  if (std::abs(z0) < m_tightTrkZ0) {
264  calculationResult["nTrkTight"] += 1;
265  }
266 
267  // From here on use only tracks with defined charge
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  // Find the maximum pt negative [0] and positive [1] tracks without z0 cut
279  const double pT = trackFitResult->getTransverseMomentum();
280  const auto& currentMaximum = maximumPtTracksWithoutZCut.at(charge);
281  if (not currentMaximum or pT > currentMaximum->pT) {
282  MaximumPtTrack newMaximum;
283  newMaximum.pT = 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  // Loose tracks
293  if (std::abs(z0) < m_looseTrkZ0) {
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  // Find the maximum pt negative [0] and positive [1] tracks
306  const double pTLoose = trackFitResult->getTransverseMomentum();
307  const auto& currentMaximumLoose = maximumPtTracks.at(charge);
308  if (not currentMaximumLoose or pTLoose > currentMaximumLoose->pT) {
309  MaximumPtTrack newMaximum;
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  //..New tighter track definitions.
321  // We can use pCMS and charge from above, because every looseB track is also a loose track
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  } // End loop over tracks
335 
336 
337  // -- Cluster variables -- //
338  std::vector<SelectedECLCluster> selectedClusters;
339  for (const ECLCluster& cluster : m_eclClusters) {
340  // Use the photon hypothesis, and only look at clusters with times < 200 ns
341  const double time = cluster.getTime();
342  if (not cluster.hasHypothesis(Belle2::ECLCluster::EHypothesisBit::c_nPhotons) or
343  std::abs(time) > 200 or
344  cluster.getEnergy(Belle2::ECLCluster::EHypothesisBit::c_nPhotons) < 0.1) {
345  continue;
346  }
347 
348  const double dt99 = cluster.getDeltaTime99();
349 
350  // Store infos if above the single photon threshold
351  SelectedECLCluster selectedCluster;
352  selectedCluster.p4Lab = cUtil.Get4MomentumFromCluster(&cluster, clustervertex,
353  Belle2::ECLCluster::EHypothesisBit::c_nPhotons);
354  selectedCluster.p4CMS = boostrotate.rotateLabToCms() * selectedCluster.p4Lab; // was clustp4COM
355  selectedCluster.cluster = &cluster;
356  selectedCluster.clusterTime = time / dt99; // was clustT
357  selectedCluster.energyLab = selectedCluster.p4Lab.E();
358  selectedCluster.energyCMS = selectedCluster.p4CMS.E(); // was first of ECOMPair
359  selectedCluster.isTrack = cluster.isTrack(); // was tempCharge
360 
361  selectedClusters.push_back(selectedCluster);
362 
363  if (selectedCluster.energyCMS > m_E2min) {
364  calculationResult["nElow"] += 1;
365  }
366  if (selectedCluster.energyCMS > m_E0min) {
367  calculationResult["nEmedium"] += 1;
368  }
369  if (selectedCluster.energyCMS > m_Ehigh) {
370  calculationResult["nEhigh"] += 1;
371  }
372 
373  if (selectedCluster.energyCMS > m_E0min and std::abs(selectedCluster.clusterTime) < 10) {
374  calculationResult["nVetoClust"] += 1;
375  }
376 
377 
378  //..Single cluster trigger objects use charge, Zernike moment, and thetaLab
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  //..For 1 track radiative Bhabha control sample
385  if (selectedCluster.energyCMS > 2. and selectedCluster.isTrack) {
386  calculationResult["chrgClust2GeV"] += 1;
387  }
388  if (selectedCluster.energyCMS > 0.45 and not selectedCluster.isTrack) {
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  // improved 3 cluster (4 cluster) trigger
397  const bool notInHighBackgroundEndcapRegion = 18.5 < thetaLab and thetaLab < 139.3;
398  if (selectedCluster.energyLab > m_EminLab and notInHighBackgroundEndcapRegion) {
399  calculationResult["nE180Lab"] += 1;
400  }
401 
402  if (selectedCluster.energyLab > m_EminLab4Cluster and notInHighBackgroundEndcapRegion) {
403  calculationResult["nE300Lab"] += 1;
404  }
405 
406  if (selectedCluster.energyLab > m_EminLab3Cluster and notInHighBackgroundEndcapRegion) {
407  calculationResult["nE500Lab"] += 1;
408  }
409 
410  if (selectedCluster.energyCMS > m_Ehigh and notInHighBackgroundEndcapRegion) {
411  calculationResult["nE2000CMS"] += 1;
412  }
413 
414  //..For two-photon fusion ALP trigger
415  if (selectedCluster.energyLab > 0.25) {
416  calculationResult["nE250Lab"] += 1;
417  }
418  if (selectedCluster.energyCMS > 4.) {
419  calculationResult["nE4000CMS"] += 1;
420  }
421 
422  //..Single cluster triggers
423  if (selectedCluster.energyCMS > m_EsinglePhoton) {
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 
451  if (selectedCluster.energyCMS > m_reducedEsinglePhoton) {
452  const bool reducedBarrelRegion = thetaLab > 44 and thetaLab < 98;
453 
454  if (photon and reducedBarrelRegion) {
455  calculationResult["nReducedEsinglePhotonReducedBarrel"] += 1;
456  }
457  }
458 
459  if (selectedCluster.energyCMS > m_Ehigh) {
460  // TODO: different definition!
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  // -- Bhabha and two-photon lepton preparation -- //
490  for (short charge : { -1, 1}) {
491  auto& maximumPtTrack = maximumPtTracks.at(charge);
492  if (not maximumPtTrack) {
493  continue;
494  }
495 
496  // Prep max two tracks for Bhabha and two-photon lepton selections
497  // Track / cluster matching
498  for (auto& cluster : maximumPtTrack->track->getRelationsTo<ECLCluster>()) {
499  if (cluster.hasHypothesis(Belle2::ECLCluster::EHypothesisBit::c_nPhotons)) {
500  maximumPtTrack->clusterEnergySumLab += cluster.getEnergy(Belle2::ECLCluster::EHypothesisBit::c_nPhotons);
501  }
502  }
503  maximumPtTrack->clusterEnergySumCMS =
504  maximumPtTrack->clusterEnergySumLab * maximumPtTrack->pCMS / maximumPtTrack->pLab;
505 
506  // Single leg Bhabha selections
507  if (maximumPtTrack->clusterEnergySumCMS > 4.5) {
508  calculationResult["ee1leg"] = 1;
509  }
510 
511  // single muon
512  if (maximumPtTrack->pCMS > 3 and maximumPtTrack->clusterEnergySumLab > 0. and maximumPtTrack->clusterEnergySumLab < 1.) {
513  calculationResult["singleMuon"] = 1;
514  }
515  }
516 
517  // Bhabha selections using two tracks
518  if (maximumPtTracks.at(-1) and maximumPtTracks.at(1)) {
519  const MaximumPtTrack& negativeTrack = *maximumPtTracks.at(-1);
520  const MaximumPtTrack& positiveTrack = *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 
527  const double negativeClusterSum = negativeTrack.clusterEnergySumCMS;
528  const double negativeClusterSumLab = negativeTrack.clusterEnergySumLab;
529  const double negativeP = negativeTrack.pCMS;
530  const double positiveClusterSum = positiveTrack.clusterEnergySumCMS;
531  const double positiveClusterSumLab = positiveTrack.clusterEnergySumLab;
532  const double positiveP = positiveTrack.pCMS;
533 
534  // two leg Bhabha
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  // one leg one track Bhabha
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  // one leg plus one electron
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  // eeee, eemm, mu mu, and radiative Bhabha selection
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();
559  const double relMissAngle0 = B2Vector3D(negativeTrack.p4CMS.Vect()).Angle(B2Vector3D(p4Miss.Vect())) * TMath::RadToDeg();
560  const double relMissAngle1 = B2Vector3D(positiveTrack.p4CMS.Vect()).Angle(B2Vector3D(p4Miss.Vect())) * TMath::RadToDeg();
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  //..Select two photon fusion lepton pairs
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  //..Veto two photon fusion lepton pairs
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  //..Radiative Bhabhas with missing momentum in acceptance
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  //..ISR radiative Bhabha with missing momentum down beam pipe
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  //..Veto Bhabha with hard brem
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  //..Veto e+e- -> mu+mu-
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  //..Select e+e- -> mu+mu-
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  //..Bhabha and gamma gamma vetoes using two maximum-energy clusters
623  if (selectedClusters.size() >= 2) {
624  const SelectedECLCluster& firstCluster = selectedClusters[0];
625  const SelectedECLCluster& secondCluster = selectedClusters[1];
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  //..Quantities for two-photon fusion triggers
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  // gg and mumu accept triggers using ECL
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  // Bhabha accept triggers.
711  // Use theta_lab of negative track if available; otherwise cluster.
712  double thetaFlatten = 0;
713 
714  // One leg, one cluster.
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.;
724  for (const SelectedECLCluster& cluster : selectedClusters) {
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  // Two legs
744  if (maximumPtTracks.at(-1) and maximumPtTracks.at(1)) {
745  const MaximumPtTrack& negativeTrack = *maximumPtTracks.at(-1);
746  const MaximumPtTrack& positiveTrack = *maximumPtTracks.at(1);
747  const double invMass = (negativeTrack.p4Lab + positiveTrack.p4Lab).M();
748  if (invMass > 5.29 and (negativeTrack.clusterEnergySumCMS > 1.5 or positiveTrack.clusterEnergySumCMS > 1.5)) {
749  calculationResult["selectee1leg1trk"] = 1;
750  }
751 
752  thetaFlatten = negativeTrack.p4Lab.Theta() * TMath::RadToDeg();
753 
754  // Two tracks but (exactly) 1 high energy cluster
755  const bool missNegClust = positiveTrack.clusterEnergySumCMS > 4.5 and negativeTrack.clusterEnergySumCMS < 1.;
756  const bool missPosClust = negativeTrack.clusterEnergySumCMS > 4.5 and positiveTrack.clusterEnergySumCMS < 1.;
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  // Single tag pi0 / eta dedicated lines
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) { // One track and at least two clusters are found
788 
789  const auto& track = maximumPtTracks.at(-1) ? *maximumPtTracks.at(-1) : *maximumPtTracks.at(1);
790  // in real signal, the pi0 decay daughters are always the two most-energetic neutral clusters.
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  //..Updated to use new track definitions
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) { // One track and at least two clusters are found
830 
831  const auto& track = maximumPtTracks.at(-1) ? *maximumPtTracks.at(-1) : *maximumPtTracks.at(1);
832  // in real signal, the pi0 decay daughters are always the two most-energetic neutral clusters.
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  //..Cosmic selection. Need exactly two tracks.
858  if (m_tracks.getEntries() == 2) {
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) {
869  return result + (cluster.hasHypothesis(Belle2::ECLCluster::EHypothesisBit::c_nPhotons) ? cluster.getEnergy(
870  Belle2::ECLCluster::EHypothesisBit::c_nPhotons) : 0);
871  };
872 
873  const auto& clustersOfNegTrack = negTrack->track->getRelationsTo<ECLCluster>();
874  const auto& clustersOfPosTrack = posTrack->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  // Select cosmic using these tracks
888  const bool goodMagneticRegion = (z0Neg<m_goodMagneticRegionZ0 or abs(d0Neg)>m_goodMagneticRegionD0)
889  and (z0Pos<m_goodMagneticRegionZ0 or abs(d0Pos)>m_goodMagneticRegionD0);
890  if (maxNegpT > m_cosmicMinPt and maxPospT > m_cosmicMinPt and maxClusterENeg < m_cosmicMaxClusterEnergy
891  and maxClusterEPos < m_cosmicMaxClusterEnergy and goodMagneticRegion) {
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  //..Displaced vertex
909  // See https://indico.belle2.org/event/8973/ and references therein
910  if (maximumPtTracksWithoutZCut.at(-1) and maximumPtTracksWithoutZCut.at(1)) {
911 
912  //..Make particles of the two highest pt tracks (without IP cut)
913  const auto nTrack = maximumPtTracksWithoutZCut.at(-1)->track;
914  Particle* nParticle = new Particle(nTrack, Const::pion);
915 
916  const auto pTrack = maximumPtTracksWithoutZCut.at(1)->track;
917  Particle* pParticle = new Particle(pTrack, Const::pion);
918 
919  //..Make a vertex of these two
920  Belle2::analysis::VertexFitKFit* vertexFit = new Belle2::analysis::VertexFitKFit();
921  vertexFit->addParticle(nParticle);
922  vertexFit->addParticle(pParticle);
923  vertexFit->doFit();
924 
925  //..Vertex properties
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  //..Angular differance of two tracks to reject cosmics
934  // Tolerance could be reduced from 10 deg to 2 deg if needed for physics reasons,
935  // for a 5% increase in the rate of selected displaced vertex triggers.
936  // See https://gitlab.desy.de/belle2/software/basf2/-/merge_requests/1867
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.; // degrees
945  const bool backToBackLab = std::abs(thetaSumLab - 180.) < backToBackTolerance and std::abs(dPhiLab - 180.) < backToBackTolerance;
946 
947  //..Select a displaced vertex
948  const double minProbChiVertex = 0.005; // minimum probability chi square. Many backgrounds have bad chi sq.
949  const double minXYVertex = 3.; // minimum xy of vertex (cm). Should pass track filters below this.
950  const double maxXYVertex = 60.; // maximum xy. Insufficient CDC for good reconstruction beyond this.
951  const double minThetaVertex = 30.; // Large Bhabha background at low angles.
952  const double maxThetaVertex = 120.; // Large Bhabha background at low angles.
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  //..Clean up
958  delete nParticle;
959  delete pParticle;
960  delete vertexFit;
961 
962  }
963 
964 }
Belle2::B2Vector3::Angle
DataType Angle(const B2Vector3< DataType > &q) const
The angle w.r.t.
Definition: B2Vector3.h:302
Belle2::ClusterUtils
Class to provide momentum-related information from ECLClusters.
Definition: ClusterUtils.h:35
Belle2::ClusterUtils::Get4MomentumFromCluster
const ROOT::Math::PxPyPzEVector Get4MomentumFromCluster(const ECLCluster *cluster, ECLCluster::EHypothesisBit hypo)
Returns four momentum vector.
Definition: ClusterUtils.cc:25
Belle2::ClusterUtils::GetIPPosition
const ROOT::Math::XYZVector GetIPPosition()
Returns default IP position from beam parameters.
Definition: ClusterUtils.cc:196
Belle2::Const::pion
static const ChargedStable pion
charged pion particle
Definition: Const.h:652
Belle2::ECLCluster
ECL cluster data.
Definition: ECLCluster.h:27
Belle2::ECLCluster::EHypothesisBit::c_nPhotons
@ c_nPhotons
CR is split into n photons (N1)
Belle2::PCmsLabTransform
Class to hold Lorentz transformations from/to CMS and boost vector.
Definition: PCmsLabTransform.h:31
Belle2::PCmsLabTransform::getCMSEnergy
double getCMSEnergy() const
Returns CMS energy of e+e- (aka.
Definition: PCmsLabTransform.h:51
Belle2::PCmsLabTransform::rotateLabToCms
const ROOT::Math::LorentzRotation rotateLabToCms() const
Returns Lorentz transformation from Lab to CMS.
Definition: PCmsLabTransform.h:69
Belle2::Particle
Class to store reconstructed particles.
Definition: Particle.h:75
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:197
Belle2::SoftwareTrigger::FilterCalculator::m_EminLab
double m_EminLab
which lab energy defines nE180Lab
Definition: FilterCalculator.h:62
Belle2::SoftwareTrigger::FilterCalculator::m_reducedEsinglePhoton
double m_reducedEsinglePhoton
which CMS energy defines nReducedEsingle clusters
Definition: FilterCalculator.h:70
Belle2::SoftwareTrigger::FilterCalculator::m_l1Trigger
StoreObjPtr< TRGSummary > m_l1Trigger
Store Object with the trigger result.
Definition: FilterCalculator.h:47
Belle2::SoftwareTrigger::FilterCalculator::m_E0min
double m_E0min
which CMS energy defines nEmedium
Definition: FilterCalculator.h:58
Belle2::SoftwareTrigger::FilterCalculator::m_goodMagneticRegionZ0
double m_goodMagneticRegionZ0
maximum z0 for well understood magnetic field (cm)
Definition: FilterCalculator.h:76
Belle2::SoftwareTrigger::FilterCalculator::m_cosmicMinPt
double m_cosmicMinPt
which LAB pt defines a cosmic
Definition: FilterCalculator.h:72
Belle2::SoftwareTrigger::FilterCalculator::m_EsinglePhoton
double m_EsinglePhoton
which CMS energy defines nEsingleClust
Definition: FilterCalculator.h:68
Belle2::SoftwareTrigger::FilterCalculator::m_E2min
double m_E2min
which CMS energy defines nElow
Definition: FilterCalculator.h:56
Belle2::SoftwareTrigger::FilterCalculator::m_tightTrkZ0
double m_tightTrkZ0
which Z0 defines a tight track
Definition: FilterCalculator.h:54
Belle2::SoftwareTrigger::FilterCalculator::m_tracks
StoreArray< Track > m_tracks
Store Array of the tracks to be used.
Definition: FilterCalculator.h:43
Belle2::SoftwareTrigger::FilterCalculator::m_EminLab3Cluster
double m_EminLab3Cluster
which lab energy defines nE500Lab
Definition: FilterCalculator.h:66
Belle2::SoftwareTrigger::FilterCalculator::m_Ehigh
double m_Ehigh
which CMS energy defines nEhigh
Definition: FilterCalculator.h:60
Belle2::SoftwareTrigger::FilterCalculator::m_eclClusters
StoreArray< ECLCluster > m_eclClusters
Store Array of the ecl clusters to be used.
Definition: FilterCalculator.h:45
Belle2::SoftwareTrigger::FilterCalculator::m_EminLab4Cluster
double m_EminLab4Cluster
which lab energy defines nE300Lab
Definition: FilterCalculator.h:64
Belle2::SoftwareTrigger::FilterCalculator::m_bitsNN
StoreArray< CDCTriggerUnpacker::NNBitStream > m_bitsNN
Store Object with the trigger NN bits.
Definition: FilterCalculator.h:49
Belle2::SoftwareTrigger::FilterCalculator::m_cosmicMaxClusterEnergy
double m_cosmicMaxClusterEnergy
which LAB cluster energy vetoes a cosmic candidate
Definition: FilterCalculator.h:74
Belle2::SoftwareTrigger::FilterCalculator::m_goodMagneticRegionD0
double m_goodMagneticRegionD0
minimum d0 for well understood magnetic field, if z0 is large (cm)
Definition: FilterCalculator.h:78
Belle2::SoftwareTrigger::FilterCalculator::m_looseTrkZ0
double m_looseTrkZ0
which Z0 defines a loose track
Definition: FilterCalculator.h:52
Belle2::TrackFitResult
Values of the result of a track fit with a given particle hypothesis.
Definition: TrackFitResult.h:49
Belle2::Track
Class that bundles various TrackFitResults.
Definition: Track.h:25
Belle2::analysis::KFitBase::getNDF
virtual int getNDF(void) const
Get an NDF of the fit.
Definition: KFitBase.cc:114
Belle2::analysis::KFitBase::addParticle
enum KFitError::ECode addParticle(const Particle *particle)
Add a particle to the fitter.
Definition: KFitBase.cc:59
Belle2::analysis::VertexFitKFit
VertexFitKFit is a derived class from KFitBase to perform vertex-constraint kinematical fit.
Definition: VertexFitKFit.h:34
Belle2::analysis::VertexFitKFit::getCHIsq
double getCHIsq(void) const override
Get a chi-square of the fit.
Definition: VertexFitKFit.cc:139
Belle2::analysis::VertexFitKFit::doFit
enum KFitError::ECode doFit(void)
Perform a vertex-constraint fit.
Definition: VertexFitKFit.cc:215
Belle2::analysis::VertexFitKFit::getVertex
const HepPoint3D getVertex(const int flag=KFitConst::kAfterFit) const
Get a vertex position.
Definition: VertexFitKFit.cc:113
Belle2::B2Vector3D
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:516
Belle2::EvtPDLUtil::charge
double charge(int pdgCode)
Returns electric charge of a particle with given pdg code.
Definition: EvtPDLUtil.cc:44
MaximumPtTrack
Temporary data structure holding the track(s) with the maximum pT.
Definition: FilterCalculator.cc:30
MaximumPtTrack::pLab
double pLab
the momentum magnitude in lab system
Definition: FilterCalculator.cc:42
MaximumPtTrack::pT
double pT
the pT of the track
Definition: FilterCalculator.cc:32
MaximumPtTrack::p4CMS
ROOT::Math::PxPyPzEVector p4CMS
the 4 momentum in CMS system
Definition: FilterCalculator.cc:44
MaximumPtTrack::clusterEnergySumLab
double clusterEnergySumLab
the sum of related cluster energies in lab system
Definition: FilterCalculator.cc:38
MaximumPtTrack::p4Lab
ROOT::Math::PxPyPzEVector p4Lab
the 4 momentum in lab system
Definition: FilterCalculator.cc:46
MaximumPtTrack::clusterEnergySumCMS
double clusterEnergySumCMS
the sum of related cluster energies in CMS system
Definition: FilterCalculator.cc:36
MaximumPtTrack::pCMS
double pCMS
the momentum magnitude in CMS system
Definition: FilterCalculator.cc:40
MaximumPtTrack::track
const Track * track
the track
Definition: FilterCalculator.cc:34
SelectedECLCluster
Temporary data structure holding the ECL clusters used for this analysis.
Definition: FilterCalculator.cc:50
SelectedECLCluster::energyLab
double energyLab
the energy in Lab system
Definition: FilterCalculator.cc:56
SelectedECLCluster::clusterTime
double clusterTime
the time of the cluster
Definition: FilterCalculator.cc:64
SelectedECLCluster::p4CMS
ROOT::Math::PxPyPzEVector p4CMS
the 4 momentum in CMS system
Definition: FilterCalculator.cc:58
SelectedECLCluster::p4Lab
ROOT::Math::PxPyPzEVector p4Lab
the 4 momentum in lab system
Definition: FilterCalculator.cc:60
SelectedECLCluster::cluster
const ECLCluster * cluster
The ECL cluster.
Definition: FilterCalculator.cc:52
SelectedECLCluster::isTrack
bool isTrack
is this ECL cluster likely from a track (or a photon) = is it charged?
Definition: FilterCalculator.cc:62
SelectedECLCluster::energyCMS
double energyCMS
the energy in CMS system
Definition: FilterCalculator.cc:54

◆ fillInCalculations()

const SoftwareTriggerObject & fillInCalculations ( )
inherited

Main function of this class: calculate the needed variables using the overwritten doCalculation function and write out the values into the results object (with their names).

Please make sure to override (or clear) the variables! Otherwise it can happen that their old values are still in the object.

What variables exactly are added to the result depends on the implementation details of the class.

Definition at line 44 of file SoftwareTriggerCalculation.cc.

◆ writeDebugOutput()

void writeDebugOutput ( const std::unique_ptr< TTree > &  debugOutputTTree)
inherited

Function to write out debug output into the given TTree.

Needs an already prefilled calculationResult for this (probably using the fillInCalculations function).

Definition at line 19 of file SoftwareTriggerCalculation.cc.

The documentation for this class was generated from the following files: