Belle II Software  release-08-00-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  calculationResult["bha3d"] = bha3d;
204  calculationResult["bhapur"] = bhapurPsnm;
205  calculationResult["bhapur_lml1"] = lml1 and bhapurFtdl;
206  } else {
207  calculationResult["l1_trigger_random"] = 1; // save every event if no L1 trigger info
208  calculationResult["l1_trigger_delayed_bhabha"] = 0;
209  calculationResult["l1_trigger_poisson"] = 0;
210  calculationResult["bha3d"] = 0;
211  calculationResult["bhapur"] = 0;
212  calculationResult["bhapur_lml1"] = 0;
213  }
214 
215  // Every 256th event has CDC NN trigger information
216  calculationResult["l1_trg_NN_info"] = 0;
217  if (m_bitsNN.isValid() and m_bitsNN.getEntries() > 0) {calculationResult["l1_trg_NN_info"] = 1;}
218 
219  calculationResult["true"] = 1;
220  calculationResult["false"] = 0;
221 
222  // Some utilities
223  // TODO: into constructor
224  ClusterUtils cUtil;
225  const ROOT::Math::XYZVector clustervertex = cUtil.GetIPPosition();
226  PCmsLabTransform boostrotate;
227  ROOT::Math::PxPyPzEVector p4ofCOM;
228  p4ofCOM.SetPxPyPzE(0, 0, 0, boostrotate.getCMSEnergy());
229 
230  // Pointers to the two tracks with the maximum pt
231  std::map<short, std::optional<MaximumPtTrack>> maximumPtTracks = {
232  { -1, {}}, {1, {}}
233  };
234 
235  // Pointer to the two tracks with maximum pt without a cut applied on z0 (used for cosmic trigger)
236  std::map<short, std::optional<MaximumPtTrack>> maximumPtTracksWithoutZCut = {
237  { -1, {}}, {1, {}}
238  };
239 
240  // --- Track variables -- //
241  for (const Track& track : m_tracks) {
242  const TrackFitResult* trackFitResult = track.getTrackFitResultWithClosestMass(Const::pion);
243  if (not trackFitResult) {
244  // Hu? No track fit result for pion? Ok, lets skip this track
245  continue;
246  }
247 
248  // TODO: Temporary cut on number of CDC hits
249  if (trackFitResult->getHitPatternCDC().getNHits() == 0) {
250  continue;
251  }
252 
253  // Count loose and tight tracks
254  const double z0 = trackFitResult->getZ0();
255  if (std::abs(z0) < m_tightTrkZ0) {
256  calculationResult["nTrkTight"] += 1;
257  }
258 
259  // From here on use only tracks with defined charge
260  const short charge = trackFitResult->getChargeSign();
261  if (charge == 0) {
262  continue;
263  }
264 
265  const ROOT::Math::PxPyPzEVector& momentumLab = trackFitResult->get4Momentum();
266  const ROOT::Math::PxPyPzEVector momentumCMS = boostrotate.rotateLabToCms() * momentumLab;
267  double pCMS = momentumCMS.P();
268  double pLab = momentumLab.P();
269 
270  // Find the maximum pt negative [0] and positive [1] tracks without z0 cut
271  const double pT = trackFitResult->getTransverseMomentum();
272  const auto& currentMaximum = maximumPtTracksWithoutZCut.at(charge);
273  if (not currentMaximum or pT > currentMaximum->pT) {
274  MaximumPtTrack newMaximum;
275  newMaximum.pT = pT;
276  newMaximum.track = &track;
277  newMaximum.pCMS = pCMS;
278  newMaximum.pLab = pLab;
279  newMaximum.p4CMS = momentumCMS;
280  newMaximum.p4Lab = momentumLab;
281  maximumPtTracksWithoutZCut[charge] = newMaximum;
282  }
283 
284  // Loose tracks
285  if (std::abs(z0) < m_looseTrkZ0) {
286  calculationResult["nTrkLoose"] += 1;
287  calculationResult["netChargeLoose"] += charge;
288 
289  if (std::isnan(calculationResult["maximumPCMS"]) or pCMS > calculationResult["maximumPCMS"]) {
290  calculationResult["maximumPCMS"] = pCMS;
291  }
292 
293  if (std::isnan(calculationResult["maximumPLab"]) or pLab > calculationResult["maximumPLab"]) {
294  calculationResult["maximumPLab"] = pLab;
295  }
296 
297  // Find the maximum pt negative [0] and positive [1] tracks
298  const double pTLoose = trackFitResult->getTransverseMomentum();
299  const auto& currentMaximumLoose = maximumPtTracks.at(charge);
300  if (not currentMaximumLoose or pTLoose > currentMaximumLoose->pT) {
301  MaximumPtTrack newMaximum;
302  newMaximum.pT = pTLoose;
303  newMaximum.track = &track;
304  newMaximum.pCMS = pCMS;
305  newMaximum.pLab = momentumLab.P();
306  newMaximum.p4CMS = momentumCMS;
307  newMaximum.p4Lab = momentumLab;
308  maximumPtTracks[charge] = newMaximum;
309  }
310  }
311 
312  //..New tighter track definitions.
313  // We can use pCMS and charge from above, because every looseB track is also a loose track
314  if (trackFitResult->getHitPatternCDC().getNHits() >= 5) {
315  if (std::abs(z0) < 1.) {calculationResult["nTrkTightB"] += 1;}
316  if (std::abs(z0) < 5.) {
317  calculationResult["nTrkLooseB"] += 1;
318  calculationResult["netChargeLooseB"] += charge;
319  if (std::isnan(calculationResult["maximumPCMSB"]) or pCMS > calculationResult["maximumPCMSB"]) {
320  calculationResult["maximumPCMSB"] = pCMS;
321 
322  }
323  }
324  }
325 
326  } // End loop over tracks
327 
328 
329  // -- Cluster variables -- //
330  std::vector<SelectedECLCluster> selectedClusters;
331  for (const ECLCluster& cluster : m_eclClusters) {
332  // Use the photon hypothesis, and only look at clusters with times < 200 ns
333  const double time = cluster.getTime();
334  if (not cluster.hasHypothesis(Belle2::ECLCluster::EHypothesisBit::c_nPhotons) or
335  std::abs(time) > 200 or
336  cluster.getEnergy(Belle2::ECLCluster::EHypothesisBit::c_nPhotons) < 0.1) {
337  continue;
338  }
339 
340  const double dt99 = cluster.getDeltaTime99();
341 
342  // Store infos if above the single photon threshold
343  SelectedECLCluster selectedCluster;
344  selectedCluster.p4Lab = cUtil.Get4MomentumFromCluster(&cluster, clustervertex,
345  Belle2::ECLCluster::EHypothesisBit::c_nPhotons);
346  selectedCluster.p4CMS = boostrotate.rotateLabToCms() * selectedCluster.p4Lab; // was clustp4COM
347  selectedCluster.cluster = &cluster;
348  selectedCluster.clusterTime = time / dt99; // was clustT
349  selectedCluster.energyLab = selectedCluster.p4Lab.E();
350  selectedCluster.energyCMS = selectedCluster.p4CMS.E(); // was first of ECOMPair
351  selectedCluster.isTrack = cluster.isTrack(); // was tempCharge
352 
353  selectedClusters.push_back(selectedCluster);
354 
355  if (selectedCluster.energyCMS > m_E2min) {
356  calculationResult["nElow"] += 1;
357  }
358  if (selectedCluster.energyCMS > m_E0min) {
359  calculationResult["nEmedium"] += 1;
360  }
361  if (selectedCluster.energyCMS > m_Ehigh) {
362  calculationResult["nEhigh"] += 1;
363  }
364 
365  if (selectedCluster.energyCMS > m_E0min and std::abs(selectedCluster.clusterTime) < 10) {
366  calculationResult["nVetoClust"] += 1;
367  }
368 
369 
370  //..Single cluster trigger objects use charge, Zernike moment, and thetaLab
371  const double thetaLab = selectedCluster.p4Lab.Theta() * TMath::RadToDeg();
372  const double zmva = cluster.getZernikeMVA();
373  const bool photon = zmva > 0.5 and not selectedCluster.isTrack;
374  const bool electron = zmva > 0.5 and selectedCluster.isTrack;
375 
376  //..For 1 track radiative Bhabha control sample
377  if (selectedCluster.energyCMS > 2. and selectedCluster.isTrack) {
378  calculationResult["chrgClust2GeV"] += 1;
379  }
380  if (selectedCluster.energyCMS > 0.45 and not selectedCluster.isTrack) {
381  const bool isInAcceptance = 17. < thetaLab and thetaLab < 150.;
382  if (isInAcceptance) {calculationResult["neutClust045GeVAcc"] += 1;}
383  const bool isInBarrel = 30. < thetaLab and thetaLab < 130.;
384  if (isInBarrel) {calculationResult["neutClust045GeVBarrel"] += 1;}
385  }
386 
387 
388  // improved 3 cluster (4 cluster) trigger
389  const bool notInHighBackgroundEndcapRegion = 18.5 < thetaLab and thetaLab < 139.3;
390  if (selectedCluster.energyLab > m_EminLab and notInHighBackgroundEndcapRegion) {
391  calculationResult["nE180Lab"] += 1;
392  }
393 
394  if (selectedCluster.energyLab > m_EminLab4Cluster and notInHighBackgroundEndcapRegion) {
395  calculationResult["nE300Lab"] += 1;
396  }
397 
398  if (selectedCluster.energyLab > m_EminLab3Cluster and notInHighBackgroundEndcapRegion) {
399  calculationResult["nE500Lab"] += 1;
400  }
401 
402  if (selectedCluster.energyCMS > m_Ehigh and notInHighBackgroundEndcapRegion) {
403  calculationResult["nE2000CMS"] += 1;
404  }
405 
406  //..For two-photon fusion ALP trigger
407  if (selectedCluster.energyLab > 0.25) {
408  calculationResult["nE250Lab"] += 1;
409  }
410  if (selectedCluster.energyCMS > 4.) {
411  calculationResult["nE4000CMS"] += 1;
412  }
413 
414  //..Single cluster triggers
415  if (selectedCluster.energyCMS > m_EsinglePhoton) {
416  calculationResult["nEsingleClust"] += 1;
417 
418  const bool barrelRegion = thetaLab > 45 and thetaLab < 115;
419  const bool extendedBarrelRegion = thetaLab > 30 and thetaLab < 130;
420  const bool endcapRegion = (thetaLab > 22 and thetaLab < 45) or (thetaLab > 115 and thetaLab < 145);
421 
422  if (photon and barrelRegion) {
423  calculationResult["nEsinglePhotonBarrel"] += 1;
424  }
425 
426  if (photon and extendedBarrelRegion) {
427  calculationResult["nEsinglePhotonExtendedBarrel"] += 1;
428  }
429 
430  if (electron and barrelRegion) {
431  calculationResult["nEsingleElectronBarrel"] += 1;
432  }
433 
434  if (electron and extendedBarrelRegion) {
435  calculationResult["nEsingleElectronExtendedBarrel"] += 1;
436  }
437 
438  if (photon and endcapRegion) {
439  calculationResult["nEsinglePhotonEndcap"] += 1;
440  }
441  }
442 
443  if (selectedCluster.energyCMS > m_reducedEsinglePhoton) {
444  const bool reducedBarrelRegion = thetaLab > 44 and thetaLab < 98;
445 
446  if (photon and reducedBarrelRegion) {
447  calculationResult["nReducedEsinglePhotonReducedBarrel"] += 1;
448  }
449  }
450 
451  if (selectedCluster.energyCMS > m_Ehigh) {
452  // TODO: different definition!
453  const bool barrelRegion = thetaLab > 32 and thetaLab < 130;
454  const bool endcapRegion = (thetaLab > 22 and thetaLab < 32) or (thetaLab > 130 and thetaLab < 145);
455  const bool lowAngleRegion = thetaLab < 22 or thetaLab > 145;
456 
457  if (not selectedCluster.isTrack and barrelRegion) {
458  calculationResult["n2GeVNeutBarrel"] += 1;
459  }
460  if (not selectedCluster.isTrack and endcapRegion) {
461  calculationResult["n2GeVNeutEndcap"] += 1;
462  }
463  if (selectedCluster.isTrack and not lowAngleRegion) {
464  calculationResult["n2GeVChrg"] += 1;
465  }
466  if (lowAngleRegion) {
467  calculationResult["nEhighLowAng"] += 1;
468  }
469  if (photon and barrelRegion) {
470  calculationResult["n2GeVPhotonBarrel"] += 1;
471  }
472  if (photon and endcapRegion) {
473  calculationResult["n2GeVPhotonEndcap"] += 1;
474  }
475  }
476  }
477  std::sort(selectedClusters.begin(), selectedClusters.end(), [](const auto & lhs, const auto & rhs) {
478  return lhs.energyCMS > rhs.energyCMS;
479  });
480 
481  // -- Bhabha and two-photon lepton preparation -- //
482  for (short charge : { -1, 1}) {
483  auto& maximumPtTrack = maximumPtTracks.at(charge);
484  if (not maximumPtTrack) {
485  continue;
486  }
487 
488  // Prep max two tracks for Bhabha and two-photon lepton selections
489  // Track / cluster matching
490  for (auto& cluster : maximumPtTrack->track->getRelationsTo<ECLCluster>()) {
491  if (cluster.hasHypothesis(Belle2::ECLCluster::EHypothesisBit::c_nPhotons)) {
492  maximumPtTrack->clusterEnergySumLab += cluster.getEnergy(Belle2::ECLCluster::EHypothesisBit::c_nPhotons);
493  }
494  }
495  maximumPtTrack->clusterEnergySumCMS =
496  maximumPtTrack->clusterEnergySumLab * maximumPtTrack->pCMS / maximumPtTrack->pLab;
497 
498  // Single leg Bhabha selections
499  if (maximumPtTrack->clusterEnergySumCMS > 4.5) {
500  calculationResult["ee1leg"] = 1;
501  }
502 
503  // single muon
504  if (maximumPtTrack->pCMS > 3 and maximumPtTrack->clusterEnergySumLab > 0. and maximumPtTrack->clusterEnergySumLab < 1.) {
505  calculationResult["singleMuon"] = 1;
506  }
507  }
508 
509  // Bhabha selections using two tracks
510  if (maximumPtTracks.at(-1) and maximumPtTracks.at(1)) {
511  const MaximumPtTrack& negativeTrack = *maximumPtTracks.at(-1);
512  const MaximumPtTrack& positiveTrack = *maximumPtTracks.at(1);
513 
514  double dphi = std::abs(negativeTrack.p4CMS.Phi() - positiveTrack.p4CMS.Phi()) * TMath::RadToDeg();
515  if (dphi > 180) {
516  dphi = 360 - dphi;
517  }
518 
519  const double negativeClusterSum = negativeTrack.clusterEnergySumCMS;
520  const double negativeClusterSumLab = negativeTrack.clusterEnergySumLab;
521  const double negativeP = negativeTrack.pCMS;
522  const double positiveClusterSum = positiveTrack.clusterEnergySumCMS;
523  const double positiveClusterSumLab = positiveTrack.clusterEnergySumLab;
524  const double positiveP = positiveTrack.pCMS;
525 
526  // two leg Bhabha
527  const double thetaSum = (negativeTrack.p4CMS.Theta() + positiveTrack.p4CMS.Theta()) * TMath::RadToDeg();
528  const double dthetaSum = std::abs(thetaSum - 180);
529  const auto back2back = dphi > 175 and dthetaSum < 15;
530  if (back2back and negativeClusterSum > 3 and positiveClusterSum > 3 and
531  (negativeClusterSum > 4.5 or positiveClusterSum > 4.5)) {
532  calculationResult["ee2leg"] = 1;
533  }
534 
535  // one leg one track Bhabha
536  if (back2back and ((negativeClusterSum > 4.5 and positiveP > 3) or (positiveClusterSum > 4.5 and negativeP > 3))) {
537  calculationResult["ee1leg1trk"] = 1;
538  }
539 
540 
541  // one leg plus one electron
542  if ((negativeClusterSum > 4.5 and positiveClusterSum > 0.8 * positiveP) or
543  (positiveClusterSum > 4.5 and negativeClusterSum > 0.8 * negativeP)) {
544  calculationResult["ee1leg1e"] = 1;
545  }
546 
547  // eeee, eemm, mu mu, and radiative Bhabha selection
548  const ROOT::Math::PxPyPzEVector p4Miss = p4ofCOM - negativeTrack.p4CMS - positiveTrack.p4CMS;
549  const double pmissTheta = p4Miss.Theta() * TMath::RadToDeg();
550  const double pmissp = p4Miss.P();
551  const double relMissAngle0 = B2Vector3D(negativeTrack.p4CMS.Vect()).Angle(B2Vector3D(p4Miss.Vect())) * TMath::RadToDeg();
552  const double relMissAngle1 = B2Vector3D(positiveTrack.p4CMS.Vect()).Angle(B2Vector3D(p4Miss.Vect())) * TMath::RadToDeg();
553 
554  const bool electronEP = positiveClusterSum > 0.8 * positiveP or negativeClusterSum > 0.8 * negativeP;
555  const bool notMuonPair = negativeClusterSumLab > 1 or positiveClusterSumLab > 1;
556  const double highp = std::max(negativeP, positiveP);
557  const double lowp = std::min(negativeP, positiveP);
558  const bool lowEdep = negativeClusterSumLab < 0.5 and positiveClusterSumLab < 0.5;
559 
560  //..Select two photon fusion lepton pairs
561  if (calculationResult["maximumPCMS"] < 2 and dphi > 160 and (pmissTheta < 25. or pmissTheta > 155.)) {
562  calculationResult["eexxSelect"] = 1;
563  if (electronEP) {
564  calculationResult["eeee"] = 1;
565  } else {
566  calculationResult["eemm"] = 1;
567  }
568  }
569 
570  //..Veto two photon fusion lepton pairs
571  if ((pmissTheta < 20. or pmissTheta > 160.) and
572  ((calculationResult["maximumPCMS"] < 1.2 and dphi > 150.) or
573  (calculationResult["maximumPCMS"] < 2. and 175. < dphi))) {
574  calculationResult["eexx"] = 1;
575  }
576 
577  //..Radiative Bhabhas with missing momentum in acceptance
578  if (negativeP > 1. and pmissTheta > 10. and pmissTheta < 170. and positiveP > 1. and dphi < 170. and pmissp > 1. and electronEP) {
579  if (calculationResult["nTrkLoose"] == 2 and calculationResult["nTrkTight"] >= 1) {
580  calculationResult["radBhabha"] = 1;
581  }
582  if (calculationResult["nTrkLooseB"] == 2 and calculationResult["nTrkTightB"] >= 1) {
583  calculationResult["radBhabhaB"] = 1;
584  }
585  }
586 
587  //..ISR radiative Bhabha with missing momentum down beam pipe
588  if (negativeP > 2. and positiveP > 2. and 2 == calculationResult["nTrkLoose"] and
589  calculationResult["nTrkTight"] >= 1 and dphi > 175. and
590  (pmissTheta < 5. or pmissTheta > 175.) and electronEP) {
591  calculationResult["isrRadBhabha"] = 1;
592  }
593 
594  //..Veto Bhabha with hard brem
595  if (highp > 4.5 and notMuonPair and pmissp > 1. and (relMissAngle0 < 5. or relMissAngle1 < 5.)) {
596  if (calculationResult["nTrkLoose"] == 2) { calculationResult["eeBrem"] = 1;}
597  if (calculationResult["nTrkLooseB"] >= 1) { calculationResult["eeBremB"] = 1;}
598  }
599 
600  //..Veto e+e- -> mu+mu-
601  if (highp > 4.5 and lowEdep and thetaSum > 175. and thetaSum < 185. and dphi > 175.) {
602  if (calculationResult["nTrkLoose"] == 2) {calculationResult["muonPairV"] = 1;}
603  if (calculationResult["nTrkLooseB"] == 2) {calculationResult["muonPairVB"] = 1;}
604  }
605 
606  //..Select e+e- -> mu+mu-
607  if (highp > 3. and lowp > 2.5 and dphi > 165. and
608  ((negativeClusterSumLab > 0. and negativeClusterSumLab < 1.) or
609  (positiveClusterSumLab > 0. and positiveClusterSumLab < 1.))) {
610  calculationResult["selectmumu"] = 1;
611  }
612  }
613 
614  //..Bhabha and gamma gamma vetoes using two maximum-energy clusters
615  if (selectedClusters.size() >= 2) {
616  const SelectedECLCluster& firstCluster = selectedClusters[0];
617  const SelectedECLCluster& secondCluster = selectedClusters[1];
618 
619  double dphi = std::abs(firstCluster.p4CMS.Phi() - secondCluster.p4CMS.Phi()) * TMath::RadToDeg();
620  if (dphi > 180) {
621  dphi = 360 - dphi;
622  }
623  double thetaSum = (firstCluster.p4CMS.Theta() + secondCluster.p4CMS.Theta()) * TMath::RadToDeg();
624  double dthetaSum = std::abs(thetaSum - 180);
625 
626  //..Quantities for two-photon fusion triggers
627  calculationResult["dphiCmsClust"] = dphi;
628  for (int ic = 0; ic < 2; ic++) {
629  const double thetaLab = selectedClusters[ic].p4Lab.Theta() * TMath::RadToDeg();
630  const bool isInAcceptance = 17. < thetaLab and thetaLab < 150.;
631  const ECLCluster* cluster = selectedClusters[ic].cluster;
632  const double zmva = cluster->getZernikeMVA();
633  const bool photon = zmva > 0.5 and not selectedClusters[ic].isTrack;
634  if (isInAcceptance and photon) {calculationResult["nMaxEPhotonAcc"] += 1;}
635  }
636 
637  const double firstEnergy = firstCluster.p4CMS.E();
638  const double secondEnergy = secondCluster.p4CMS.E();
639 
640  const bool highEnergetic = firstEnergy > 3 and secondEnergy > 3 and (firstEnergy > 4.5 or secondEnergy > 4.5);
641 
642  if (dphi > 160 and dphi < 178 and dthetaSum < 15 and highEnergetic) {
643  calculationResult["ee2clst"] = 1;
644  }
645 
646  if (dphi > 178 and dthetaSum < 15 and highEnergetic) {
647  calculationResult["gg2clst"] = 1;
648  }
649 
650  if ((calculationResult["ee2clst"] == 1 or calculationResult["gg2clst"] == 1) and
651  calculationResult["ee1leg"] == 1) {
652  calculationResult["ee1leg1clst"] = 1;
653  }
654 
655  const double Elab0 = firstCluster.p4Lab.E();
656  const double Elab1 = secondCluster.p4Lab.E();
657 
658  // gg and mumu accept triggers using ECL
659  if (firstEnergy > 2 and secondEnergy > 2) {
660  const double thetaLab0 = firstCluster.p4Lab.Theta() * TMath::RadToDeg();
661  const double thetaLab1 = secondCluster.p4Lab.Theta() * TMath::RadToDeg();
662 
663  const bool barrel0 = 32. < thetaLab0 and thetaLab0 < 130.;
664  const bool barrel1 = 32. < thetaLab1 and thetaLab1 < 130.;
665  const bool oneClustersAbove4 = firstEnergy > 4 or secondEnergy > 4;
666  const bool oneIsNeutral = not firstCluster.isTrack or not secondCluster.isTrack;
667  const bool bothAreNeutral = not firstCluster.isTrack and not secondCluster.isTrack;
668  const bool oneIsBarrel = barrel0 or barrel1;
669  const bool dphiCutExtraLoose = dphi > 175;
670  const bool dphiCutLoose = dphi > 177;
671  const bool dphiCutTight = dphi > 177.5;
672 
673  if (dphiCutExtraLoose and oneIsNeutral and oneIsBarrel) {
674  calculationResult["ggBarrelVL"] = 1;
675  }
676  if (oneClustersAbove4 and dphiCutLoose and oneIsNeutral and oneIsBarrel) {
677  calculationResult["ggBarrelLoose"] = 1;
678  }
679  if (oneClustersAbove4 and dphiCutTight and bothAreNeutral and oneIsBarrel) {
680  calculationResult["ggBarrelTight"] = 1;
681  }
682  if (dphiCutExtraLoose and oneIsNeutral and not oneIsBarrel) {
683  calculationResult["ggEndcapVL"] = 1;
684  }
685  if (oneClustersAbove4 and dphiCutLoose and oneIsNeutral and not oneIsBarrel) {
686  calculationResult["ggEndcapLoose"] = 1;
687  }
688  if (oneClustersAbove4 and dphiCutTight and bothAreNeutral and not oneIsBarrel) {
689  calculationResult["ggEndcapTight"] = 1;
690  }
691  }
692 
693  const double minEnergy = std::min(Elab0, Elab1);
694  const double maxEnergy = std::max(Elab0, Elab1);
695  if (dphi > 155 and thetaSum > 165 and thetaSum < 195 and minEnergy > 0.15 and minEnergy < 0.5 and
696  maxEnergy > 0.15 and maxEnergy < 0.5) {
697  calculationResult["muonPairECL"] = 1;
698  }
699 
700  }
701 
702  // Bhabha accept triggers.
703  // Use theta_lab of negative track if available; otherwise cluster.
704  double thetaFlatten = 0;
705 
706  // One leg, one cluster.
707  for (short charge : { -1, 1}) {
708  const auto& maximumPtTrack = maximumPtTracks.at(charge);
709  if (not maximumPtTrack) {
710  continue;
711  }
712 
713  if (maximumPtTrack->clusterEnergySumCMS > 1.5) {
714  double invMass = 0.;
715  double tempFlatten = 0.;
716  for (const SelectedECLCluster& cluster : selectedClusters) {
717  double tempInvMass = (maximumPtTrack->p4Lab + cluster.p4Lab).M();
718  if (tempInvMass > invMass) {
719  invMass = tempInvMass;
720  if (charge == 1) {
721  tempFlatten = cluster.p4Lab.Theta() * TMath::RadToDeg();
722  }
723  }
724  }
725  if (charge == -1) {
726  tempFlatten = maximumPtTrack->p4Lab.Theta() * TMath::RadToDeg();
727  }
728  if (invMass > 5.29) {
729  calculationResult["selectee1leg1clst"] = 1;
730  thetaFlatten = tempFlatten;
731  }
732  }
733  }
734 
735  // Two legs
736  if (maximumPtTracks.at(-1) and maximumPtTracks.at(1)) {
737  const MaximumPtTrack& negativeTrack = *maximumPtTracks.at(-1);
738  const MaximumPtTrack& positiveTrack = *maximumPtTracks.at(1);
739  const double invMass = (negativeTrack.p4Lab + positiveTrack.p4Lab).M();
740  if (invMass > 5.29 and (negativeTrack.clusterEnergySumCMS > 1.5 or positiveTrack.clusterEnergySumCMS > 1.5)) {
741  calculationResult["selectee1leg1trk"] = 1;
742  }
743 
744  thetaFlatten = negativeTrack.p4Lab.Theta() * TMath::RadToDeg();
745 
746  // Two tracks but (exactly) 1 high energy cluster
747  const bool missNegClust = positiveTrack.clusterEnergySumCMS > 4.5 and negativeTrack.clusterEnergySumCMS < 1.;
748  const bool missPosClust = negativeTrack.clusterEnergySumCMS > 4.5 and positiveTrack.clusterEnergySumCMS < 1.;
749  if ((invMass > 9.) and (missNegClust or missPosClust)) {
750  calculationResult["eeOneClust"] = 1;
751  }
752  }
753 
754  if (calculationResult["selectee1leg1trk"] == 1 or calculationResult["selectee1leg1clst"] == 1) {
755  for (int iflat = 0; iflat < 9; iflat++) {
756  const std::string& eeFlatName = "eeFlat" + std::to_string(iflat);
757  calculationResult[eeFlatName] =
758  thetaFlatten > flatBoundaries[iflat] and thetaFlatten < flatBoundaries[iflat + 1];
759  if (calculationResult[eeFlatName]) {
760  calculationResult["selectee"] = 1;
761  }
762  }
763  }
764 
765  // Single tag pi0 / eta dedicated lines
766  if (calculationResult["nTrkLoose"] == 1 and calculationResult["maximumPCMS"] > 0.8 and selectedClusters.size() >= 2) {
767 
768  decltype(selectedClusters) selectedSingleTagClusters(selectedClusters.size());
769  auto lastItem = std::copy_if(selectedClusters.begin(), selectedClusters.end(), selectedSingleTagClusters.begin(),
770  [](auto & cluster) {
771  const bool isNeutralCluster = not cluster.isTrack;
772  const bool hasEnoughEnergy = cluster.energyLab > 0.1;
773  const double clusterThetaLab = cluster.p4Lab.Theta() * TMath::RadToDeg();
774  const bool isInAcceptance = 17 < clusterThetaLab and clusterThetaLab < 150.;
775  return isNeutralCluster and hasEnoughEnergy and isInAcceptance;
776  });
777  selectedSingleTagClusters.resize(std::distance(selectedSingleTagClusters.begin(), lastItem));
778 
779  if (selectedSingleTagClusters.size() >= 2) { // One track and at least two clusters are found
780 
781  const auto& track = maximumPtTracks.at(-1) ? *maximumPtTracks.at(-1) : *maximumPtTracks.at(1);
782  // in real signal, the pi0 decay daughters are always the two most-energetic neutral clusters.
783  const auto& firstCluster = selectedSingleTagClusters[0];
784  const auto& secondCluster = selectedSingleTagClusters[1];
785 
786  const ROOT::Math::PxPyPzEVector trackP4CMS = track.p4CMS;
787  const ROOT::Math::PxPyPzEVector pi0P4CMS = firstCluster.p4CMS + secondCluster.p4CMS;
788 
789  const bool passPi0ECMS = pi0P4CMS.E() > 1. and pi0P4CMS.E() < 0.525 * p4ofCOM.M();
790  const double thetaSumCMS = (pi0P4CMS.Theta() + trackP4CMS.Theta()) * TMath::RadToDeg();
791  const bool passThetaSum = thetaSumCMS < 170. or thetaSumCMS > 190.;
792 
793  double dphiCMS = std::abs(trackP4CMS.Phi() - pi0P4CMS.Phi()) * TMath::RadToDeg();
794  if (dphiCMS > 180) {
795  dphiCMS = 360 - dphiCMS;
796  }
797  const bool passdPhi = dphiCMS > 160.;
798 
799  if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() < 0.7) {
800  calculationResult["singleTagLowMass"] = 1;
801  } else if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() > 0.7) {
802  calculationResult["singleTagHighMass"] = 1;
803  }
804  }
805  }
806 
807  //..Updated to use new track definitions
808  if (calculationResult["nTrkLooseB"] == 1 and calculationResult["maximumPCMSB"] > 0.8 and selectedClusters.size() >= 2) {
809 
810  decltype(selectedClusters) selectedSingleTagClusters(selectedClusters.size());
811  auto lastItem = std::copy_if(selectedClusters.begin(), selectedClusters.end(), selectedSingleTagClusters.begin(),
812  [](auto & cluster) {
813  const bool isNeutralCluster = not cluster.isTrack;
814  const bool hasEnoughEnergy = cluster.energyLab > 0.1;
815  const double clusterThetaLab = cluster.p4Lab.Theta() * TMath::RadToDeg();
816  const bool isInAcceptance = 17 < clusterThetaLab and clusterThetaLab < 150.;
817  return isNeutralCluster and hasEnoughEnergy and isInAcceptance;
818  });
819  selectedSingleTagClusters.resize(std::distance(selectedSingleTagClusters.begin(), lastItem));
820 
821  if (selectedSingleTagClusters.size() >= 2) { // One track and at least two clusters are found
822 
823  const auto& track = maximumPtTracks.at(-1) ? *maximumPtTracks.at(-1) : *maximumPtTracks.at(1);
824  // in real signal, the pi0 decay daughters are always the two most-energetic neutral clusters.
825  const auto& firstCluster = selectedSingleTagClusters[0];
826  const auto& secondCluster = selectedSingleTagClusters[1];
827 
828  const ROOT::Math::PxPyPzEVector trackP4CMS = track.p4CMS;
829  const ROOT::Math::PxPyPzEVector pi0P4CMS = firstCluster.p4CMS + secondCluster.p4CMS;
830 
831  const bool passPi0ECMS = pi0P4CMS.E() > 1. and pi0P4CMS.E() < 0.525 * p4ofCOM.M();
832  const double thetaSumCMS = (pi0P4CMS.Theta() + trackP4CMS.Theta()) * TMath::RadToDeg();
833  const bool passThetaSum = thetaSumCMS < 170. or thetaSumCMS > 190.;
834 
835  double dphiCMS = std::abs(trackP4CMS.Phi() - pi0P4CMS.Phi()) * TMath::RadToDeg();
836  if (dphiCMS > 180) {
837  dphiCMS = 360 - dphiCMS;
838  }
839  const bool passdPhi = dphiCMS > 160.;
840 
841  if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() < 0.7) {
842  calculationResult["singleTagLowMassB"] = 1;
843  } else if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() > 0.7) {
844  calculationResult["singleTagHighMassB"] = 1;
845  }
846  }
847  }
848 
849  //..Cosmic selection. Need exactly two tracks.
850  if (m_tracks.getEntries() == 2) {
851 
852  const auto negTrack = maximumPtTracksWithoutZCut.at(-1);
853  const auto posTrack = maximumPtTracksWithoutZCut.at(1);
854 
855  if (negTrack and posTrack) {
856 
857  const double maxNegpT = negTrack->pT;
858  const double maxPospT = posTrack->pT;
859 
860  auto accumulatePhotonEnergy = [](double result, const auto & cluster) {
861  return result + (cluster.hasHypothesis(Belle2::ECLCluster::EHypothesisBit::c_nPhotons) ? cluster.getEnergy(
862  Belle2::ECLCluster::EHypothesisBit::c_nPhotons) : 0);
863  };
864 
865  const auto& clustersOfNegTrack = negTrack->track->getRelationsTo<ECLCluster>();
866  const auto& clustersOfPosTrack = posTrack->track->getRelationsTo<ECLCluster>();
867 
868  const double maxClusterENeg = std::accumulate(clustersOfNegTrack.begin(), clustersOfNegTrack.end(), 0.0, accumulatePhotonEnergy);
869  const double maxClusterEPos = std::accumulate(clustersOfPosTrack.begin(), clustersOfPosTrack.end(), 0.0, accumulatePhotonEnergy);
870 
871  const ROOT::Math::PxPyPzEVector& momentumLabNeg(negTrack->p4Lab);
872  const ROOT::Math::PxPyPzEVector& momentumLabPos(posTrack->p4Lab);
873 
874  const double& z0Neg = negTrack->track->getTrackFitResultWithClosestMass(Const::pion)->getZ0();
875  const double& d0Neg = negTrack->track->getTrackFitResultWithClosestMass(Const::pion)->getD0();
876  const double& z0Pos = posTrack->track->getTrackFitResultWithClosestMass(Const::pion)->getZ0();
877  const double& d0Pos = posTrack->track->getTrackFitResultWithClosestMass(Const::pion)->getD0();
878 
879  // Select cosmic using these tracks
880  const bool goodMagneticRegion = (z0Neg<m_goodMagneticRegionZ0 or abs(d0Neg)>m_goodMagneticRegionD0)
881  and (z0Pos<m_goodMagneticRegionZ0 or abs(d0Pos)>m_goodMagneticRegionD0);
882  if (maxNegpT > m_cosmicMinPt and maxPospT > m_cosmicMinPt and maxClusterENeg < m_cosmicMaxClusterEnergy
883  and maxClusterEPos < m_cosmicMaxClusterEnergy and goodMagneticRegion) {
884  double dphiLab = std::abs(momentumLabNeg.Phi() - momentumLabPos.Phi()) * TMath::RadToDeg();
885  if (dphiLab > 180) {
886  dphiLab = 360 - dphiLab;
887  }
888 
889  const double thetaSumLab = (momentumLabNeg.Theta() + momentumLabPos.Theta()) * TMath::RadToDeg();
890 
891  constexpr double phiBackToBackTolerance = 2.;
892  constexpr double thetaBackToBackTolerance = 2.;
893  if ((180 - dphiLab) < phiBackToBackTolerance and std::abs(180 - thetaSumLab) < thetaBackToBackTolerance) {
894  calculationResult["cosmic"] = 1;
895  }
896  }
897  }
898  }
899 
900  //..Displaced vertex
901  // See https://indico.belle2.org/event/8973/ and references therein
902  if (maximumPtTracksWithoutZCut.at(-1) and maximumPtTracksWithoutZCut.at(1)) {
903 
904  //..Make particles of the two highest pt tracks (without IP cut)
905  const auto nTrack = maximumPtTracksWithoutZCut.at(-1)->track;
906  Particle* nParticle = new Particle(nTrack, Const::pion);
907 
908  const auto pTrack = maximumPtTracksWithoutZCut.at(1)->track;
909  Particle* pParticle = new Particle(pTrack, Const::pion);
910 
911  //..Make a vertex of these two
912  Belle2::analysis::VertexFitKFit* vertexFit = new Belle2::analysis::VertexFitKFit();
913  vertexFit->addParticle(nParticle);
914  vertexFit->addParticle(pParticle);
915  vertexFit->doFit();
916 
917  //..Vertex properties
918  const double chisq = vertexFit->getCHIsq();
919  const int ndf = vertexFit->getNDF();
920  const double vertexProb = TMath::Prob(chisq, ndf);
921  const auto vertexLocation = vertexFit->getVertex();
922  const double vertexXY = vertexLocation.perp();
923  const double vertexTheta = vertexLocation.theta() * TMath::RadToDeg();
924 
925  //..Angular differance of two tracks to reject cosmics
926  // Tolerance could be reduced from 10 deg to 2 deg if needed for physics reasons,
927  // for a 5% increase in the rate of selected displaced vertex triggers.
928  // See https://gitlab.desy.de/belle2/software/basf2/-/merge_requests/1867
929  const ROOT::Math::PxPyPzEVector& momentumLabNeg(maximumPtTracksWithoutZCut.at(-1)->p4Lab);
930  const ROOT::Math::PxPyPzEVector& momentumLabPos(maximumPtTracksWithoutZCut.at(1)->p4Lab);
931  const double thetaSumLab = (momentumLabNeg.Theta() + momentumLabPos.Theta()) * TMath::RadToDeg();
932  double dPhiLab = std::abs(momentumLabNeg.Phi() - momentumLabPos.Phi()) * TMath::RadToDeg();
933  if (dPhiLab > 180) {
934  dPhiLab = 360 - dPhiLab;
935  }
936  const double backToBackTolerance = 10.; // degrees
937  const bool backToBackLab = std::abs(thetaSumLab - 180.) < backToBackTolerance and std::abs(dPhiLab - 180.) < backToBackTolerance;
938 
939  //..Select a displaced vertex
940  const double minProbChiVertex = 0.005; // minimum probability chi square. Many backgrounds have bad chi sq.
941  const double minXYVertex = 3.; // minimum xy of vertex (cm). Should pass track filters below this.
942  const double maxXYVertex = 60.; // maximum xy. Insufficient CDC for good reconstruction beyond this.
943  const double minThetaVertex = 30.; // Large Bhabha background at low angles.
944  const double maxThetaVertex = 120.; // Large Bhabha background at low angles.
945  if (vertexProb > minProbChiVertex and vertexXY > minXYVertex and vertexXY < maxXYVertex and vertexTheta > minThetaVertex
946  and vertexTheta < maxThetaVertex
947  and not backToBackLab) {calculationResult["displacedVertex"] = 1;}
948 
949  //..Clean up
950  delete nParticle;
951  delete pParticle;
952  delete vertexFit;
953 
954  }
955 
956 }
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: