Belle II Software development
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:
SoftwareTriggerCalculation

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.
 
void writeDebugOutput (const std::unique_ptr< TTree > &debugOutputTTree)
 Function to write out debug output into the given TTree.
 
void addDebugOutput (const StoreObjPtr< SoftwareTriggerVariables > &storeObject, const std::string &prefix)
 Function to write out debug output into the given StoreObject.
 
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).
 

Private Attributes

StoreArray< Trackm_tracks
 Store Array of the tracks to be used.
 
StoreArray< ECLClusterm_eclClusters
 Store Array of the ecl clusters to be used.
 
StoreObjPtr< TRGSummarym_l1Trigger
 Store Object with the trigger result.
 
StoreArray< CDCTriggerUnpacker::NNBitStreamm_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.

Constructor & Destructor Documentation

◆ FilterCalculator()

Set the default names for the store object particle lists.

Definition at line 71 of file FilterCalculator.cc.

71 : m_bitsNN("CDCTriggerNNBits")
72{
73
74}
StoreArray< CDCTriggerUnpacker::NNBitStream > m_bitsNN
Store Object with the trigger NN bits.

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 }
SoftwareTriggerObject m_calculationResult
Internal storage of the result of the calculation.

◆ 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,
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(
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
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}
DataType Angle(const B2Vector3< DataType > &q) const
The angle w.r.t.
Definition: B2Vector3.h:302
Class to provide momentum-related information from ECLClusters.
Definition: ClusterUtils.h:36
const ROOT::Math::PxPyPzEVector Get4MomentumFromCluster(const ECLCluster *cluster, ECLCluster::EHypothesisBit hypo)
Returns four momentum vector.
Definition: ClusterUtils.cc:25
const ROOT::Math::XYZVector GetIPPosition()
Returns default IP position from beam parameters.
static const ChargedStable pion
charged pion particle
Definition: Const.h:661
ECL cluster data.
Definition: ECLCluster.h:27
@ c_nPhotons
CR is split into n photons (N1)
Class to hold Lorentz transformations from/to CMS and boost vector.
double getCMSEnergy() const
Returns CMS energy of e+e- (aka.
const ROOT::Math::LorentzRotation rotateLabToCms() const
Returns Lorentz transformation from Lab to CMS.
Class to store reconstructed particles.
Definition: Particle.h:75
RelationVector< TO > getRelationsTo(const std::string &name="", const std::string &namedRelation="") const
Get the relations that point from this object to another store array.
double m_EminLab
which lab energy defines nE180Lab
double m_reducedEsinglePhoton
which CMS energy defines nReducedEsingle clusters
StoreObjPtr< TRGSummary > m_l1Trigger
Store Object with the trigger result.
double m_E0min
which CMS energy defines nEmedium
double m_goodMagneticRegionZ0
maximum z0 for well understood magnetic field (cm)
double m_cosmicMinPt
which LAB pt defines a cosmic
double m_EsinglePhoton
which CMS energy defines nEsingleClust
double m_E2min
which CMS energy defines nElow
double m_tightTrkZ0
which Z0 defines a tight track
StoreArray< Track > m_tracks
Store Array of the tracks to be used.
double m_EminLab3Cluster
which lab energy defines nE500Lab
double m_Ehigh
which CMS energy defines nEhigh
StoreArray< ECLCluster > m_eclClusters
Store Array of the ecl clusters to be used.
double m_EminLab4Cluster
which lab energy defines nE300Lab
double m_cosmicMaxClusterEnergy
which LAB cluster energy vetoes a cosmic candidate
double m_goodMagneticRegionD0
minimum d0 for well understood magnetic field, if z0 is large (cm)
double m_looseTrkZ0
which Z0 defines a loose track
@ TTYP_DPHY
delayed physics events for background
Definition: TRGSummary.h:65
@ TTYP_POIS
poisson random trigger
Definition: TRGSummary.h:73
@ TTYP_RAND
random trigger events
Definition: TRGSummary.h:67
Values of the result of a track fit with a given particle hypothesis.
Class that bundles various TrackFitResults.
Definition: Track.h:25
virtual int getNDF(void) const
Get an NDF of the fit.
Definition: KFitBase.cc:114
enum KFitError::ECode addParticle(const Particle *particle)
Add a particle to the fitter.
Definition: KFitBase.cc:59
VertexFitKFit is a derived class from KFitBase to perform vertex-constraint kinematical fit.
Definition: VertexFitKFit.h:34
double getCHIsq(void) const override
Get a chi-square of the fit.
enum KFitError::ECode doFit(void)
Perform a vertex-constraint fit.
const HepPoint3D getVertex(const int flag=KFitConst::kAfterFit) const
Get a vertex position.
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:516
double charge(int pdgCode)
Returns electric charge of a particle with given pdg code.
Definition: EvtPDLUtil.cc:44
Temporary data structure holding the track(s) with the maximum pT.
double pLab
the momentum magnitude in lab system
double pT
the pT of the track
ROOT::Math::PxPyPzEVector p4CMS
the 4 momentum in CMS system
double clusterEnergySumLab
the sum of related cluster energies in lab system
ROOT::Math::PxPyPzEVector p4Lab
the 4 momentum in lab system
double clusterEnergySumCMS
the sum of related cluster energies in CMS system
double pCMS
the momentum magnitude in CMS system
const Track * track
the track
Temporary data structure holding the ECL clusters used for this analysis.
double energyLab
the energy in Lab system
double clusterTime
the time of the cluster
ROOT::Math::PxPyPzEVector p4CMS
the 4 momentum in CMS system
ROOT::Math::PxPyPzEVector p4Lab
the 4 momentum in lab system
const ECLCluster * cluster
The ECL cluster.
bool isTrack
is this ECL cluster likely from a track (or a photon) = is it charged?
double energyCMS
the energy in CMS system

◆ 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.

45 {
46 const unsigned int sizeBeforeCheck = m_calculationResult.size();
48
49 if (m_calculationResult.size() != sizeBeforeCheck and sizeBeforeCheck > 0) {
50 B2WARNING("Calculator added more variables (" << m_calculationResult.size() <<
51 ") than there were before (" << sizeBeforeCheck << "). Probably something strange is going on!");
52 }
53
55 }
virtual void doCalculation(SoftwareTriggerObject &m_calculationResult)=0
Override this function to implement your calculation.

◆ requireStoreArrays()

void requireStoreArrays ( )
overridevirtual

Require the particle list. We do not need more here.

Implements SoftwareTriggerCalculation.

Definition at line 76 of file FilterCalculator.cc.

77{
78 m_tracks.isRequired();
79 m_eclClusters.isRequired();
80 m_l1Trigger.isOptional();
81 m_bitsNN.isOptional();
82}

◆ 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.

20 {
21 if (not m_debugPrepared) {
22 for (auto& identifierWithValue : m_calculationResult) {
23 const std::string& identifier = identifierWithValue.first;
24 double& value = identifierWithValue.second;
25
26 debugOutputTTree->Branch(identifier.c_str(), &value);
27 }
28 m_debugPrepared = true;
29 }
30
31 debugOutputTTree->Fill();
32 }
bool m_debugPrepared
Flag to not add the branches twice to the TTree.

Member Data Documentation

◆ m_bitsNN

Store Object with the trigger NN bits.

Definition at line 49 of file FilterCalculator.h.

◆ m_calculationResult

SoftwareTriggerObject m_calculationResult
privateinherited

Internal storage of the result of the calculation.

Definition at line 74 of file SoftwareTriggerCalculation.h.

◆ m_cosmicMaxClusterEnergy

double m_cosmicMaxClusterEnergy = 1.0 * Unit::GeV
private

which LAB cluster energy vetoes a cosmic candidate

Definition at line 74 of file FilterCalculator.h.

◆ m_cosmicMinPt

double m_cosmicMinPt = 0.5 * Unit::GeV
private

which LAB pt defines a cosmic

Definition at line 72 of file FilterCalculator.h.

◆ m_debugPrepared

bool m_debugPrepared = false
privateinherited

Flag to not add the branches twice to the TTree.

Definition at line 76 of file SoftwareTriggerCalculation.h.

◆ m_E0min

double m_E0min = 0.3
private

which CMS energy defines nEmedium

Definition at line 58 of file FilterCalculator.h.

◆ m_E2min

double m_E2min = 0.2
private

which CMS energy defines nElow

Definition at line 56 of file FilterCalculator.h.

◆ m_eclClusters

StoreArray<ECLCluster> m_eclClusters
private

Store Array of the ecl clusters to be used.

Definition at line 45 of file FilterCalculator.h.

◆ m_Ehigh

double m_Ehigh = 2
private

which CMS energy defines nEhigh

Definition at line 60 of file FilterCalculator.h.

◆ m_EminLab

double m_EminLab = 0.18
private

which lab energy defines nE180Lab

Definition at line 62 of file FilterCalculator.h.

◆ m_EminLab3Cluster

double m_EminLab3Cluster = 0.5
private

which lab energy defines nE500Lab

Definition at line 66 of file FilterCalculator.h.

◆ m_EminLab4Cluster

double m_EminLab4Cluster = 0.3
private

which lab energy defines nE300Lab

Definition at line 64 of file FilterCalculator.h.

◆ m_EsinglePhoton

double m_EsinglePhoton = 1
private

which CMS energy defines nEsingleClust

Definition at line 68 of file FilterCalculator.h.

◆ m_goodMagneticRegionD0

double m_goodMagneticRegionD0 = 26.5
private

minimum d0 for well understood magnetic field, if z0 is large (cm)

Definition at line 78 of file FilterCalculator.h.

◆ m_goodMagneticRegionZ0

double m_goodMagneticRegionZ0 = 57.
private

maximum z0 for well understood magnetic field (cm)

Definition at line 76 of file FilterCalculator.h.

◆ m_l1Trigger

StoreObjPtr<TRGSummary> m_l1Trigger
private

Store Object with the trigger result.

Definition at line 47 of file FilterCalculator.h.

◆ m_looseTrkZ0

double m_looseTrkZ0 = 10 * Unit::cm
private

which Z0 defines a loose track

Definition at line 52 of file FilterCalculator.h.

◆ m_reducedEsinglePhoton

double m_reducedEsinglePhoton = 0.5
private

which CMS energy defines nReducedEsingle clusters

Definition at line 70 of file FilterCalculator.h.

◆ m_tightTrkZ0

double m_tightTrkZ0 = 2 * Unit::cm
private

which Z0 defines a tight track

Definition at line 54 of file FilterCalculator.h.

◆ m_tracks

StoreArray<Track> m_tracks
private

Store Array of the tracks to be used.

Definition at line 43 of file FilterCalculator.h.


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