Belle II Software development
FilterCalculator.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8#include <hlt/softwaretrigger/calculations/FilterCalculator.h>
9#include <hlt/softwaretrigger/calculations/utilities.h>
10
11#include <Math/Vector3D.h>
12#include <Math/Vector4D.h>
13#include <Math/VectorUtil.h>
14
15#include <mdst/dataobjects/TrackFitResult.h>
16#include <mdst/dataobjects/HitPatternCDC.h>
17
18#include <analysis/dataobjects/Particle.h>
19#include <analysis/VertexFitting/KFit/VertexFitKFit.h>
20
21#include <framework/logging/Logger.h>
22
23#include <cmath>
24#include <numeric>
25#include <stdexcept>
26#include <optional>
27
28using namespace Belle2;
29using namespace SoftwareTrigger;
30
34 double pT = NAN;
36 const Track* track = nullptr;
42 double pCMS = NAN;
44 double pLab = NAN;
46 ROOT::Math::PxPyPzEVector p4CMS;
48 ROOT::Math::PxPyPzEVector p4Lab;
49};
50
54 const ECLCluster* cluster = nullptr;
56 double energyCMS = NAN;
58 double energyLab = NAN;
60 ROOT::Math::PxPyPzEVector p4CMS;
62 ROOT::Math::PxPyPzEVector p4Lab;
64 bool isTrack = false;
66 double clusterTime = NAN;
67};
68
70const double flatBoundaries[10] = {0., 19., 22., 25., 30., 35., 45., 60., 90., 180.};
71
72
74{
75
76}
77
79{
80 m_tracks.isRequired();
81 m_eclClusters.isRequired();
82 m_l1Trigger.isOptional();
83 m_bitsNN.isOptional();
84}
85
86void FilterCalculator::doCalculation(SoftwareTriggerObject& calculationResult)
87{
88 calculationResult["nTrkLoose"] = 0;
89 calculationResult["nTrkTight"] = 0;
90 calculationResult["ee2leg"] = 0;
91 calculationResult["nEmedium"] = 0;
92 calculationResult["nElow"] = 0;
93 calculationResult["nEhigh"] = 0;
94 calculationResult["nE180Lab"] = 0;
95 calculationResult["nE300Lab"] = 0;
96 calculationResult["nE500Lab"] = 0;
97 calculationResult["nE2000CMS"] = 0;
98 calculationResult["nE4000CMS"] = 0;
99 calculationResult["nE250Lab"] = 0;
100 calculationResult["nMaxEPhotonAcc"] = 0;
101 calculationResult["dphiCmsClust"] = NAN;
102
103 calculationResult["netChargeLoose"] = 0;
104 calculationResult["maximumPCMS"] = NAN;
105 calculationResult["maximumPLab"] = NAN;
106 calculationResult["eexx"] = 0;
107 calculationResult["ee1leg1trk"] = 0;
108 calculationResult["nEhighLowAng"] = 0;
109 calculationResult["nEsingleClust"] = 0;
110 calculationResult["nEsinglePhotonBarrel"] = 0;
111 calculationResult["nEsinglePhotonExtendedBarrel"] = 0;
112 calculationResult["nEsinglePhotonEndcap"] = 0;
113 calculationResult["nEsingleElectronBarrel"] = 0;
114 calculationResult["nEsingleElectronExtendedBarrel"] = 0;
115 calculationResult["nReducedEsinglePhotonReducedBarrel"] = 0;
116 calculationResult["nVetoClust"] = 0;
117 calculationResult["chrgClust2GeV"] = 0;
118 calculationResult["neutClust045GeVAcc"] = 0;
119 calculationResult["neutClust045GeVBarrel"] = 0;
120 calculationResult["singleTagLowMass"] = 0;
121 calculationResult["singleTagHighMass"] = 0;
122 calculationResult["n2GeVNeutBarrel"] = 0;
123 calculationResult["n2GeVNeutEndcap"] = 0;
124 calculationResult["n2GeVChrg"] = 0;
125 calculationResult["n2GeVPhotonBarrel"] = 0;
126 calculationResult["n2GeVPhotonEndcap"] = 0;
127 calculationResult["ee1leg"] = 0;
128 calculationResult["ee1leg1clst"] = 0;
129 calculationResult["ee1leg1e"] = 0;
130 calculationResult["ee2clst"] = 0;
131 calculationResult["gg2clst"] = 0;
132 calculationResult["eeee"] = 0;
133 calculationResult["eemm"] = 0;
134 calculationResult["eexxSelect"] = 0;
135 calculationResult["radBhabha"] = 0;
136 calculationResult["eeBrem"] = 0;
137 calculationResult["isrRadBhabha"] = 0;
138 calculationResult["muonPairECL"] = 0;
139 calculationResult["ggHighPt"] = 0;
140 calculationResult["selectee1leg1trk"] = 0;
141 calculationResult["selectee1leg1clst"] = 0;
142 calculationResult["selectee"] = 0;
143 calculationResult["ggBarrelVL"] = 0;
144 calculationResult["ggBarrelLoose"] = 0;
145 calculationResult["ggBarrelTight"] = 0;
146 calculationResult["ggEndcapVL"] = 0;
147 calculationResult["ggEndcapLoose"] = 0;
148 calculationResult["ggEndcapTight"] = 0;
149 calculationResult["muonPairV"] = 0;
150 calculationResult["selectmumu"] = 0;
151 calculationResult["singleMuon"] = 0;
152 calculationResult["cosmic"] = 0;
153 calculationResult["displacedVertex"] = 0;
154 calculationResult["eeFlat0"] = 0;
155 calculationResult["eeFlat1"] = 0;
156 calculationResult["eeFlat2"] = 0;
157 calculationResult["eeFlat3"] = 0;
158 calculationResult["eeFlat4"] = 0;
159 calculationResult["eeFlat5"] = 0;
160 calculationResult["eeFlat6"] = 0;
161 calculationResult["eeFlat7"] = 0;
162 calculationResult["eeFlat8"] = 0;
163 calculationResult["eeOneClust"] = 0;
164
165 //..New track definitions
166 calculationResult["nTrkLooseB"] = 0;
167 calculationResult["nTrkTightB"] = 0;
168 calculationResult["maximumPCMSB"] = NAN;
169 calculationResult["netChargeLooseB"] = 0;
170 calculationResult["muonPairVB"] = 0;
171 calculationResult["eeBremB"] = 0;
172 calculationResult["singleTagLowMassB"] = 0;
173 calculationResult["singleTagHighMassB"] = 0;
174 calculationResult["radBhabhaB"] = 0;
175
176 //..Tracks without a requirement on nCDCHits
177 calculationResult["nTrackC"] = 0;
178 calculationResult["maximumPCMSC"] = NAN;
179 calculationResult["pCmsNegC"] = NAN;
180 calculationResult["clusterENegC"] = NAN;
181 calculationResult["pCmsPosC"] = NAN;
182 calculationResult["clusterEPosC"] = NAN;
183 calculationResult["dPhiCmsC"] = NAN;
184
185
186 // Passed on L1 information
187 if (m_l1Trigger.isValid()) {
188 calculationResult["l1_trigger_random"] = m_l1Trigger->getTimType() == TRGSummary::ETimingType::TTYP_RAND;
189 calculationResult["l1_trigger_delayed_bhabha"] = m_l1Trigger->getTimType() == TRGSummary::ETimingType::TTYP_DPHY;
190 calculationResult["l1_trigger_poisson"] = m_l1Trigger->getTimType() == TRGSummary::ETimingType::TTYP_POIS;
191 bool bha3d;
192 try {
193 bha3d = m_l1Trigger->testPsnm("bha3d");
194 } catch (const std::exception&) {
195 bha3d = false;
196 }
197 bool bhapurPsnm;
198 try {
199 bhapurPsnm = m_l1Trigger->testPsnm("bhapur");
200 } catch (const std::exception&) {
201 bhapurPsnm = false;
202 }
203 bool bhapurFtdl;
204 try {
205 bhapurFtdl = m_l1Trigger->testFtdl("bhapur");
206 } catch (const std::exception&) {
207 bhapurFtdl = false;
208 }
209 bool lml1;
210 try {
211 lml1 = m_l1Trigger->testPsnm("lml1");
212 } catch (const std::exception&) {
213 lml1 = false;
214 }
215 bool l1_bit_f;
216 try {
217 l1_bit_f = m_l1Trigger->testPsnm("fpre");
218 } catch (const std::exception&) {
219 l1_bit_f = false;
220 }
221 calculationResult["bha3d"] = bha3d;
222 calculationResult["bhapur"] = bhapurPsnm;
223 calculationResult["bhapur_lml1"] = lml1 and bhapurFtdl;
224 calculationResult["l1_bit_f"] = l1_bit_f;
225 } else {
226 calculationResult["l1_trigger_random"] = 1; // save every event if no L1 trigger info
227 calculationResult["l1_trigger_delayed_bhabha"] = 0;
228 calculationResult["l1_trigger_poisson"] = 0;
229 calculationResult["bha3d"] = 0;
230 calculationResult["bhapur"] = 0;
231 calculationResult["bhapur_lml1"] = 0;
232 calculationResult["l1_bit_f"] = 0;
233 }
234
235 // Every 256th event has CDC NN trigger information
236 calculationResult["l1_trg_NN_info"] = 0;
237 if (m_bitsNN.isValid() and m_bitsNN.getEntries() > 0) {calculationResult["l1_trg_NN_info"] = 1;}
238
239 calculationResult["true"] = 1;
240 calculationResult["false"] = 0;
241
242 // Some utilities
243 // TODO: into constructor
244 ClusterUtils cUtil;
245 const ROOT::Math::XYZVector clustervertex = cUtil.GetIPPosition();
246 PCmsLabTransform boostrotate;
247 ROOT::Math::PxPyPzEVector p4ofCOM;
248 p4ofCOM.SetPxPyPzE(0, 0, 0, boostrotate.getCMSEnergy());
249
250 // Pointers to the two tracks with the maximum pt
251 std::map<short, std::optional<MaximumPtTrack>> maximumPtTracks = {
252 { -1, {}}, {1, {}}
253 };
254
255 // Pointer to the two tracks with maximum pt without a cut applied on z0 (used for cosmic trigger)
256 std::map<short, std::optional<MaximumPtTrack>> maximumPtTracksWithoutZCut = {
257 { -1, {}}, {1, {}}
258 };
259
260 //..Pointers to the two "C" tracks with maximum pCms (no requirement on nCDC hits)
261 std::map<short, std::optional<MaximumPtTrack>> maximumPCmsTracksC = {
262 { -1, {}}, {1, {}}
263 };
264
265 //------------------------------------------------------------
266 // --- Track variables -- //
267 for (const Track& track : m_tracks) {
268 const TrackFitResult* trackFitResult = track.getTrackFitResultWithClosestMass(Const::pion);
269 if (not trackFitResult) {
270 // Huh? No track fit result for pion? Ok, lets skip this track
271 continue;
272 }
273
274 //..Track properties. Only use tracks with defined charge.
275 const short charge = trackFitResult->getChargeSign();
276 if (charge == 0) {continue;}
277 const double z0 = trackFitResult->getZ0();
278 const ROOT::Math::PxPyPzEVector& momentumLab = trackFitResult->get4Momentum();
279 const ROOT::Math::PxPyPzEVector momentumCMS = boostrotate.rotateLabToCms() * momentumLab;
280 const double pCMS = momentumCMS.P();
281 const double pLab = momentumLab.P();
282 const double trackTime = track.getTrackTime();
283 const double pT = trackFitResult->getTransverseMomentum();
284 const int nCDCHits = (int)(0.5 + trackFitResult->getHitPatternCDC().getNHits());
285 double clusterELab = 0.;
286 for (auto& cluster : track.getRelationsTo<ECLCluster>()) {
287 if (cluster.hasHypothesis(Belle2::ECLCluster::EHypothesisBit::c_nPhotons)) {
288 clusterELab += cluster.getEnergy(Belle2::ECLCluster::EHypothesisBit::c_nPhotons);
289 }
290 }
291 const double clusterECMS = clusterELab * pCMS / pLab;
292
293
294 //------------------------------------------------------------
295 //..New list without a requirement on number of CDC hits
296 bool goodTrackCTime = true;
297 if (std::abs(trackTime) > 15.) {goodTrackCTime = false;}
298 if (std::abs(z0) < 1. and goodTrackCTime and pCMS > 0.2) {
299 calculationResult["nTrackC"] += 1;
300 if (std::isnan(calculationResult["maximumPCMSC"]) or pCMS > calculationResult["maximumPCMSC"]) {
301 calculationResult["maximumPCMSC"] = pCMS;
302 }
303 }
304
305 //..Store maximum pCms positive and negative tracks
306 const auto& currentMaximumC = maximumPCmsTracksC.at(charge);
307 if (not currentMaximumC or pCMS > currentMaximumC->pCMS) {
308 MaximumPtTrack newMaximum;
309 newMaximum.pT = pT;
310 newMaximum.track = &track;
311 newMaximum.pCMS = pCMS;
312 newMaximum.pLab = pLab;
313 newMaximum.p4CMS = momentumCMS;
314 newMaximum.p4Lab = momentumLab;
315 newMaximum.clusterEnergySumLab = clusterELab;
316 newMaximum.clusterEnergySumCMS = clusterECMS;
317 maximumPCmsTracksC[charge] = newMaximum;
318 }
319
320 //------------------------------------------------------------
321 //..Now original loose and tight tracks, which require nCDC hits
322 // TODO: Temporary cut on number of CDC hits
323 if (nCDCHits == 0) {
324 continue;
325 }
326
327 // Count tight tracks
328 if (std::abs(z0) < m_tightTrkZ0) {
329 calculationResult["nTrkTight"] += 1;
330 }
331
332
333 //------------------------------------------------------------
334 // Find the maximum pt negative [0] and positive [1] tracks without z0 cut
335 const auto& currentMaximum = maximumPtTracksWithoutZCut.at(charge);
336 if (not currentMaximum or pT > currentMaximum->pT) {
337 MaximumPtTrack newMaximum;
338 newMaximum.pT = pT;
339 newMaximum.track = &track;
340 newMaximum.pCMS = pCMS;
341 newMaximum.pLab = pLab;
342 newMaximum.p4CMS = momentumCMS;
343 newMaximum.p4Lab = momentumLab;
344 newMaximum.clusterEnergySumLab = clusterELab;
345 newMaximum.clusterEnergySumCMS = clusterECMS;
346 maximumPtTracksWithoutZCut[charge] = newMaximum;
347 }
348
349
350 //------------------------------------------------------------
351 // Loose tracks
352 if (std::abs(z0) < m_looseTrkZ0) {
353 calculationResult["nTrkLoose"] += 1;
354 calculationResult["netChargeLoose"] += charge;
355
356 if (std::isnan(calculationResult["maximumPCMS"]) or pCMS > calculationResult["maximumPCMS"]) {
357 calculationResult["maximumPCMS"] = pCMS;
358 }
359
360 if (std::isnan(calculationResult["maximumPLab"]) or pLab > calculationResult["maximumPLab"]) {
361 calculationResult["maximumPLab"] = pLab;
362 }
363
364 // Find the maximum pt negative [0] and positive [1] tracks
365 const double pTLoose = trackFitResult->getTransverseMomentum();
366 const auto& currentMaximumLoose = maximumPtTracks.at(charge);
367 if (not currentMaximumLoose or pTLoose > currentMaximumLoose->pT) {
368 MaximumPtTrack newMaximum;
369 newMaximum.pT = pTLoose;
370 newMaximum.track = &track;
371 newMaximum.pCMS = pCMS;
372 newMaximum.pLab = momentumLab.P();
373 newMaximum.p4CMS = momentumCMS;
374 newMaximum.p4Lab = momentumLab;
375 newMaximum.clusterEnergySumLab = clusterELab;
376 newMaximum.clusterEnergySumCMS = clusterECMS;
377 maximumPtTracks[charge] = newMaximum;
378 }
379 }
380
381 //..New tighter track definitions.
382 // We can use pCMS and charge from above, because every looseB track is also a loose track
383 if (nCDCHits >= 5) {
384 if (std::abs(z0) < 1.) {calculationResult["nTrkTightB"] += 1;}
385 if (std::abs(z0) < 5.) {
386 calculationResult["nTrkLooseB"] += 1;
387 calculationResult["netChargeLooseB"] += charge;
388 if (std::isnan(calculationResult["maximumPCMSB"]) or pCMS > calculationResult["maximumPCMSB"]) {
389 calculationResult["maximumPCMSB"] = pCMS;
390
391 }
392 }
393 }
394
395 } // End loop over tracks
396
397
398 //------------------------------------------------------------
399 //..Store the "C" track calculations
400 for (short charge : { -1, 1}) {
401 auto& maximumPcmsTrackC = maximumPCmsTracksC.at(charge);
402 if (not maximumPcmsTrackC) {
403 continue;
404 }
405
406 //..Calculations for each track
407 if (charge == -1) {
408 calculationResult["pCmsNegC"] = maximumPcmsTrackC->pCMS;
409 calculationResult["clusterENegC"] = maximumPcmsTrackC->clusterEnergySumLab;
410 } else {
411 calculationResult["pCmsPosC"] = maximumPcmsTrackC->pCMS;
412 calculationResult["clusterEPosC"] = maximumPcmsTrackC->clusterEnergySumLab;
413 }
414 }
415
416 //..Calculations using both highest-pCms C tracks
417 if (maximumPCmsTracksC.at(-1) and maximumPCmsTracksC.at(1)) {
418 const MaximumPtTrack& negativeTrack = *maximumPCmsTracksC.at(-1);
419 const MaximumPtTrack& positiveTrack = *maximumPCmsTracksC.at(1);
420
421 calculationResult["dPhiCmsC"] = std::abs(ROOT::Math::VectorUtil::DeltaPhi(negativeTrack.p4CMS,
422 positiveTrack.p4CMS)) * TMath::RadToDeg();
423 }
424
425
426 //------------------------------------------------------------
427 // -- Cluster variables -- //
428 std::vector<SelectedECLCluster> selectedClusters;
429 for (const ECLCluster& cluster : m_eclClusters) {
430 // Use the photon hypothesis, and only look at clusters with times < 200 ns
431 const double time = cluster.getTime();
432 if (not cluster.hasHypothesis(Belle2::ECLCluster::EHypothesisBit::c_nPhotons) or
433 std::abs(time) > 200 or
434 cluster.getEnergy(Belle2::ECLCluster::EHypothesisBit::c_nPhotons) < 0.1) {
435 continue;
436 }
437
438 const double dt99 = cluster.getDeltaTime99();
439
440 // Store cluster information
441 SelectedECLCluster selectedCluster;
442 selectedCluster.p4Lab = cUtil.Get4MomentumFromCluster(&cluster, clustervertex,
444 selectedCluster.p4CMS = boostrotate.rotateLabToCms() * selectedCluster.p4Lab; // was clustp4COM
445 selectedCluster.cluster = &cluster;
446 selectedCluster.clusterTime = time / dt99; // was clustT
447 selectedCluster.energyLab = selectedCluster.p4Lab.E();
448 selectedCluster.energyCMS = selectedCluster.p4CMS.E(); // was first of ECOMPair
449 selectedCluster.isTrack = cluster.isTrack(); // was tempCharge
450
451 selectedClusters.push_back(selectedCluster);
452
453 if (selectedCluster.energyCMS > m_E2min) {
454 calculationResult["nElow"] += 1;
455 }
456 if (selectedCluster.energyCMS > m_E0min) {
457 calculationResult["nEmedium"] += 1;
458 }
459 if (selectedCluster.energyCMS > m_Ehigh) {
460 calculationResult["nEhigh"] += 1;
461 }
462
463 if (selectedCluster.energyCMS > m_E0min and std::abs(selectedCluster.clusterTime) < 10) {
464 calculationResult["nVetoClust"] += 1;
465 }
466
467
468 //..Single cluster trigger objects use charge, Zernike moment, and thetaLab
469 const double thetaLab = selectedCluster.p4Lab.Theta() * TMath::RadToDeg();
470 const double zmva = cluster.getZernikeMVA();
471 const bool photon = zmva > 0.5 and not selectedCluster.isTrack;
472 const bool electron = zmva > 0.5 and selectedCluster.isTrack;
473
474 //..For 1 track radiative Bhabha control sample
475 if (selectedCluster.energyCMS > 2. and selectedCluster.isTrack) {
476 calculationResult["chrgClust2GeV"] += 1;
477 }
478 if (selectedCluster.energyCMS > 0.45 and not selectedCluster.isTrack) {
479 const bool isInAcceptance = 17. < thetaLab and thetaLab < 150.;
480 if (isInAcceptance) {calculationResult["neutClust045GeVAcc"] += 1;}
481 const bool isInBarrel = 30. < thetaLab and thetaLab < 130.;
482 if (isInBarrel) {calculationResult["neutClust045GeVBarrel"] += 1;}
483 }
484
485
486 // improved 3 cluster (4 cluster) trigger
487 const bool notInHighBackgroundEndcapRegion = 18.5 < thetaLab and thetaLab < 139.3;
488 if (selectedCluster.energyLab > m_EminLab and notInHighBackgroundEndcapRegion) {
489 calculationResult["nE180Lab"] += 1;
490 }
491
492 if (selectedCluster.energyLab > m_EminLab4Cluster and notInHighBackgroundEndcapRegion) {
493 calculationResult["nE300Lab"] += 1;
494 }
495
496 if (selectedCluster.energyLab > m_EminLab3Cluster and notInHighBackgroundEndcapRegion) {
497 calculationResult["nE500Lab"] += 1;
498 }
499
500 if (selectedCluster.energyCMS > m_Ehigh and notInHighBackgroundEndcapRegion) {
501 calculationResult["nE2000CMS"] += 1;
502 }
503
504 //..For two-photon fusion ALP trigger
505 if (selectedCluster.energyLab > 0.25) {
506 calculationResult["nE250Lab"] += 1;
507 }
508 if (selectedCluster.energyCMS > 4.) {
509 calculationResult["nE4000CMS"] += 1;
510 }
511
512 //..Single cluster triggers
513 if (selectedCluster.energyCMS > m_EsinglePhoton) {
514 calculationResult["nEsingleClust"] += 1;
515
516 const bool barrelRegion = thetaLab > 45 and thetaLab < 115;
517 const bool extendedBarrelRegion = thetaLab > 30 and thetaLab < 130;
518 const bool endcapRegion = (thetaLab > 22 and thetaLab < 45) or (thetaLab > 115 and thetaLab < 145);
519
520 if (photon and barrelRegion) {
521 calculationResult["nEsinglePhotonBarrel"] += 1;
522 }
523
524 if (photon and extendedBarrelRegion) {
525 calculationResult["nEsinglePhotonExtendedBarrel"] += 1;
526 }
527
528 if (electron and barrelRegion) {
529 calculationResult["nEsingleElectronBarrel"] += 1;
530 }
531
532 if (electron and extendedBarrelRegion) {
533 calculationResult["nEsingleElectronExtendedBarrel"] += 1;
534 }
535
536 if (photon and endcapRegion) {
537 calculationResult["nEsinglePhotonEndcap"] += 1;
538 }
539 }
540
541 if (selectedCluster.energyCMS > m_reducedEsinglePhoton) {
542 const bool reducedBarrelRegion = thetaLab > 44 and thetaLab < 98;
543
544 if (photon and reducedBarrelRegion) {
545 calculationResult["nReducedEsinglePhotonReducedBarrel"] += 1;
546 }
547 }
548
549 if (selectedCluster.energyCMS > m_Ehigh) {
550 // TODO: different definition!
551 const bool barrelRegion = thetaLab > 32 and thetaLab < 130;
552 const bool endcapRegion = (thetaLab > 22 and thetaLab < 32) or (thetaLab > 130 and thetaLab < 145);
553 const bool lowAngleRegion = thetaLab < 22 or thetaLab > 145;
554
555 if (not selectedCluster.isTrack and barrelRegion) {
556 calculationResult["n2GeVNeutBarrel"] += 1;
557 }
558 if (not selectedCluster.isTrack and endcapRegion) {
559 calculationResult["n2GeVNeutEndcap"] += 1;
560 }
561 if (selectedCluster.isTrack and not lowAngleRegion) {
562 calculationResult["n2GeVChrg"] += 1;
563 }
564 if (lowAngleRegion) {
565 calculationResult["nEhighLowAng"] += 1;
566 }
567 if (photon and barrelRegion) {
568 calculationResult["n2GeVPhotonBarrel"] += 1;
569 }
570 if (photon and endcapRegion) {
571 calculationResult["n2GeVPhotonEndcap"] += 1;
572 }
573 }
574 } // end of loop over clusters
575
576 //..Order clusters by CMS energy
577 std::sort(selectedClusters.begin(), selectedClusters.end(), [](const auto & lhs, const auto & rhs) {
578 return lhs.energyCMS > rhs.energyCMS;
579 });
580
581 // -- Bhabha and two-photon lepton preparation -- //
582 for (short charge : { -1, 1}) {
583 auto& maximumPtTrack = maximumPtTracks.at(charge);
584 if (not maximumPtTrack) {
585 continue;
586 }
587
588 // Single leg Bhabha selections
589 if (maximumPtTrack->clusterEnergySumCMS > 4.5) {
590 calculationResult["ee1leg"] = 1;
591 }
592
593 // single muon
594 if (maximumPtTrack->pCMS > 3 and maximumPtTrack->clusterEnergySumLab > 0. and maximumPtTrack->clusterEnergySumLab < 1.) {
595 calculationResult["singleMuon"] = 1;
596 }
597 }
598
599 // Bhabha selections using two tracks
600 if (maximumPtTracks.at(-1) and maximumPtTracks.at(1)) {
601 const MaximumPtTrack& negativeTrack = *maximumPtTracks.at(-1);
602 const MaximumPtTrack& positiveTrack = *maximumPtTracks.at(1);
603
604 double dphi = std::abs(negativeTrack.p4CMS.Phi() - positiveTrack.p4CMS.Phi()) * TMath::RadToDeg();
605 if (dphi > 180) {
606 dphi = 360 - dphi;
607 }
608
609 const double negativeClusterSum = negativeTrack.clusterEnergySumCMS;
610 const double negativeClusterSumLab = negativeTrack.clusterEnergySumLab;
611 const double negativeP = negativeTrack.pCMS;
612 const double positiveClusterSum = positiveTrack.clusterEnergySumCMS;
613 const double positiveClusterSumLab = positiveTrack.clusterEnergySumLab;
614 const double positiveP = positiveTrack.pCMS;
615
616 // two leg Bhabha
617 const double thetaSum = (negativeTrack.p4CMS.Theta() + positiveTrack.p4CMS.Theta()) * TMath::RadToDeg();
618 const double dthetaSum = std::abs(thetaSum - 180);
619 const auto back2back = dphi > 175 and dthetaSum < 15;
620 if (back2back and negativeClusterSum > 3 and positiveClusterSum > 3 and
621 (negativeClusterSum > 4.5 or positiveClusterSum > 4.5)) {
622 calculationResult["ee2leg"] = 1;
623 }
624
625 // one leg one track Bhabha
626 if (back2back and ((negativeClusterSum > 4.5 and positiveP > 3) or (positiveClusterSum > 4.5 and negativeP > 3))) {
627 calculationResult["ee1leg1trk"] = 1;
628 }
629
630
631 // one leg plus one electron
632 if ((negativeClusterSum > 4.5 and positiveClusterSum > 0.8 * positiveP) or
633 (positiveClusterSum > 4.5 and negativeClusterSum > 0.8 * negativeP)) {
634 calculationResult["ee1leg1e"] = 1;
635 }
636
637 // eeee, eemm, mu mu, and radiative Bhabha selection
638 const ROOT::Math::PxPyPzEVector p4Miss = p4ofCOM - negativeTrack.p4CMS - positiveTrack.p4CMS;
639 const double pmissTheta = p4Miss.Theta() * TMath::RadToDeg();
640 const double pmissp = p4Miss.P();
641 const double relMissAngle0 = ROOT::Math::VectorUtil::Angle(negativeTrack.p4CMS, p4Miss) * TMath::RadToDeg();
642 const double relMissAngle1 = ROOT::Math::VectorUtil::Angle(positiveTrack.p4CMS, p4Miss) * TMath::RadToDeg();
643
644 const bool electronEP = positiveClusterSum > 0.8 * positiveP or negativeClusterSum > 0.8 * negativeP;
645 const bool notMuonPair = negativeClusterSumLab > 1 or positiveClusterSumLab > 1;
646 const double highp = std::max(negativeP, positiveP);
647 const double lowp = std::min(negativeP, positiveP);
648 const bool lowEdep = negativeClusterSumLab < 0.5 and positiveClusterSumLab < 0.5;
649
650 //..Select two photon fusion lepton pairs
651 if (calculationResult["maximumPCMS"] < 2 and dphi > 160 and (pmissTheta < 25. or pmissTheta > 155.)) {
652 calculationResult["eexxSelect"] = 1;
653 if (electronEP) {
654 calculationResult["eeee"] = 1;
655 } else {
656 calculationResult["eemm"] = 1;
657 }
658 }
659
660 //..Veto two photon fusion lepton pairs
661 if ((pmissTheta < 20. or pmissTheta > 160.) and
662 ((calculationResult["maximumPCMS"] < 1.2 and dphi > 150.) or
663 (calculationResult["maximumPCMS"] < 2. and 175. < dphi))) {
664 calculationResult["eexx"] = 1;
665 }
666
667 //..Radiative Bhabhas with missing momentum in acceptance
668 if (negativeP > 1. and pmissTheta > 10. and pmissTheta < 170. and positiveP > 1. and dphi < 170. and pmissp > 1. and electronEP) {
669 if (calculationResult["nTrkLoose"] == 2 and calculationResult["nTrkTight"] >= 1) {
670 calculationResult["radBhabha"] = 1;
671 }
672 if (calculationResult["nTrkLooseB"] == 2 and calculationResult["nTrkTightB"] >= 1) {
673 calculationResult["radBhabhaB"] = 1;
674 }
675 }
676
677 //..ISR radiative Bhabha with missing momentum down beam pipe
678 if (negativeP > 2. and positiveP > 2. and 2 == calculationResult["nTrkLoose"] and
679 calculationResult["nTrkTight"] >= 1 and dphi > 175. and
680 (pmissTheta < 5. or pmissTheta > 175.) and electronEP) {
681 calculationResult["isrRadBhabha"] = 1;
682 }
683
684 //..Veto Bhabha with hard brem
685 if (highp > 4.5 and notMuonPair and pmissp > 1. and (relMissAngle0 < 5. or relMissAngle1 < 5.)) {
686 if (calculationResult["nTrkLoose"] == 2) { calculationResult["eeBrem"] = 1;}
687 if (calculationResult["nTrkLooseB"] >= 1) { calculationResult["eeBremB"] = 1;}
688 }
689
690 //..Veto e+e- -> mu+mu-
691 if (highp > 4.5 and lowEdep and thetaSum > 175. and thetaSum < 185. and dphi > 175.) {
692 if (calculationResult["nTrkLoose"] == 2) {calculationResult["muonPairV"] = 1;}
693 if (calculationResult["nTrkLooseB"] == 2) {calculationResult["muonPairVB"] = 1;}
694 }
695
696 //..Select e+e- -> mu+mu-
697 if (highp > 3. and lowp > 2.5 and dphi > 165. and
698 ((negativeClusterSumLab > 0. and negativeClusterSumLab < 1.) or
699 (positiveClusterSumLab > 0. and positiveClusterSumLab < 1.))) {
700 calculationResult["selectmumu"] = 1;
701 }
702 }
703
704 //..Filters and vetoes using two maximum-energy clusters
705 if (selectedClusters.size() >= 2) {
706 const SelectedECLCluster& firstCluster = selectedClusters[0];
707 const SelectedECLCluster& secondCluster = selectedClusters[1];
708
709 double dphi = std::abs(firstCluster.p4CMS.Phi() - secondCluster.p4CMS.Phi()) * TMath::RadToDeg();
710 if (dphi > 180) {
711 dphi = 360 - dphi;
712 }
713 double thetaSum = (firstCluster.p4CMS.Theta() + secondCluster.p4CMS.Theta()) * TMath::RadToDeg();
714 double dthetaSum = std::abs(thetaSum - 180);
715
716 //..Quantities for two-photon fusion triggers
717 calculationResult["dphiCmsClust"] = dphi;
718 for (int ic = 0; ic < 2; ic++) {
719 const double thetaLab = selectedClusters[ic].p4Lab.Theta() * TMath::RadToDeg();
720 const bool isInAcceptance = 17. < thetaLab and thetaLab < 150.;
721 const ECLCluster* cluster = selectedClusters[ic].cluster;
722 const double zmva = cluster->getZernikeMVA();
723 const bool photon = zmva > 0.5 and not selectedClusters[ic].isTrack;
724 if (isInAcceptance and photon) {calculationResult["nMaxEPhotonAcc"] += 1;}
725 }
726
727 const double firstEnergy = firstCluster.p4CMS.E();
728 const double secondEnergy = secondCluster.p4CMS.E();
729
730 const bool highEnergetic = firstEnergy > 3 and secondEnergy > 3 and (firstEnergy > 4.5 or secondEnergy > 4.5);
731
732 if (dphi > 160 and dphi < 178 and dthetaSum < 15 and highEnergetic) {
733 calculationResult["ee2clst"] = 1;
734 }
735
736 if (dphi > 178 and dthetaSum < 15 and highEnergetic) {
737 calculationResult["gg2clst"] = 1;
738 }
739
740 if ((calculationResult["ee2clst"] == 1 or calculationResult["gg2clst"] == 1) and
741 calculationResult["ee1leg"] == 1) {
742 calculationResult["ee1leg1clst"] = 1;
743 }
744
745 const double Elab0 = firstCluster.p4Lab.E();
746 const double Elab1 = secondCluster.p4Lab.E();
747
748 // gg and mumu accept triggers using ECL
749 if (firstEnergy > 2 and secondEnergy > 2) {
750 const double thetaLab0 = firstCluster.p4Lab.Theta() * TMath::RadToDeg();
751 const double thetaLab1 = secondCluster.p4Lab.Theta() * TMath::RadToDeg();
752
753 const bool barrel0 = 32. < thetaLab0 and thetaLab0 < 130.;
754 const bool barrel1 = 32. < thetaLab1 and thetaLab1 < 130.;
755 const bool oneClustersAbove4 = firstEnergy > 4 or secondEnergy > 4;
756 const bool oneIsNeutral = not firstCluster.isTrack or not secondCluster.isTrack;
757 const bool bothAreNeutral = not firstCluster.isTrack and not secondCluster.isTrack;
758 const bool oneIsBarrel = barrel0 or barrel1;
759 const bool dphiCutExtraLoose = dphi > 175;
760 const bool dphiCutLoose = dphi > 177;
761 const bool dphiCutTight = dphi > 177.5;
762
763 if (dphiCutExtraLoose and oneIsNeutral and oneIsBarrel) {
764 calculationResult["ggBarrelVL"] = 1;
765 }
766 if (oneClustersAbove4 and dphiCutLoose and oneIsNeutral and oneIsBarrel) {
767 calculationResult["ggBarrelLoose"] = 1;
768 }
769 if (oneClustersAbove4 and dphiCutTight and bothAreNeutral and oneIsBarrel) {
770 calculationResult["ggBarrelTight"] = 1;
771 }
772 if (dphiCutExtraLoose and oneIsNeutral and not oneIsBarrel) {
773 calculationResult["ggEndcapVL"] = 1;
774 }
775 if (oneClustersAbove4 and dphiCutLoose and oneIsNeutral and not oneIsBarrel) {
776 calculationResult["ggEndcapLoose"] = 1;
777 }
778 if (oneClustersAbove4 and dphiCutTight and bothAreNeutral and not oneIsBarrel) {
779 calculationResult["ggEndcapTight"] = 1;
780 }
781 }
782
783 const double minEnergy = std::min(Elab0, Elab1);
784 const double maxEnergy = std::max(Elab0, Elab1);
785 if (dphi > 155 and thetaSum > 165 and thetaSum < 195 and minEnergy > 0.15 and minEnergy < 0.5 and
786 maxEnergy > 0.15 and maxEnergy < 0.5) {
787 calculationResult["muonPairECL"] = 1;
788 }
789
790 //..diPhoton line
791 const double thetaLab0 = firstCluster.p4Lab.Theta() * TMath::RadToDeg();
792 const double thetaLab1 = secondCluster.p4Lab.Theta() * TMath::RadToDeg();
793 const bool inHieRegion0 = thetaLab0 > 26. and thetaLab0 < 130.;
794 const bool inHieRegion1 = thetaLab1 > 26. and thetaLab1 < 130.;
795 const bool firstIsNeutral = not firstCluster.isTrack;
796 const bool secondIsNeutral = not secondCluster.isTrack;
797
798 if (secondEnergy > 0.3 and inHieRegion0 and inHieRegion1 and firstIsNeutral and secondIsNeutral) {
799 const ROOT::Math::PxPyPzEVector ggP4CMS = firstCluster.p4CMS + secondCluster.p4CMS;
800 if (ggP4CMS.pt() > 1.) {calculationResult["ggHighPt"] = 1;}
801 }
802
803 } // end of two-clusters lines
804
805 // Bhabha accept triggers.
806 // Use theta_lab of negative track if available; otherwise cluster.
807 double thetaFlatten = 0;
808
809 // One leg, one cluster.
810 for (short charge : { -1, 1}) {
811 const auto& maximumPtTrack = maximumPtTracks.at(charge);
812 if (not maximumPtTrack) {
813 continue;
814 }
815
816 if (maximumPtTrack->clusterEnergySumCMS > 1.5) {
817 double invMass = 0.;
818 double tempFlatten = 0.;
819 for (const SelectedECLCluster& cluster : selectedClusters) {
820 double tempInvMass = (maximumPtTrack->p4Lab + cluster.p4Lab).M();
821 if (tempInvMass > invMass) {
822 invMass = tempInvMass;
823 if (charge == 1) {
824 tempFlatten = cluster.p4Lab.Theta() * TMath::RadToDeg();
825 }
826 }
827 }
828 if (charge == -1) {
829 tempFlatten = maximumPtTrack->p4Lab.Theta() * TMath::RadToDeg();
830 }
831 if (invMass > 5.29) {
832 calculationResult["selectee1leg1clst"] = 1;
833 thetaFlatten = tempFlatten;
834 }
835 }
836 }
837
838 // Two legs
839 if (maximumPtTracks.at(-1) and maximumPtTracks.at(1)) {
840 const MaximumPtTrack& negativeTrack = *maximumPtTracks.at(-1);
841 const MaximumPtTrack& positiveTrack = *maximumPtTracks.at(1);
842 const double invMass = (negativeTrack.p4Lab + positiveTrack.p4Lab).M();
843 if (invMass > 5.29 and (negativeTrack.clusterEnergySumCMS > 1.5 or positiveTrack.clusterEnergySumCMS > 1.5)) {
844 calculationResult["selectee1leg1trk"] = 1;
845 }
846
847 thetaFlatten = negativeTrack.p4Lab.Theta() * TMath::RadToDeg();
848
849 // Two tracks but (exactly) 1 high energy cluster
850 const bool missNegClust = positiveTrack.clusterEnergySumCMS > 4.5 and negativeTrack.clusterEnergySumCMS < 1.;
851 const bool missPosClust = negativeTrack.clusterEnergySumCMS > 4.5 and positiveTrack.clusterEnergySumCMS < 1.;
852 if ((invMass > 9.) and (missNegClust or missPosClust)) {
853 calculationResult["eeOneClust"] = 1;
854 }
855 }
856
857 if (calculationResult["selectee1leg1trk"] == 1 or calculationResult["selectee1leg1clst"] == 1) {
858 for (int iflat = 0; iflat < 9; iflat++) {
859 const std::string& eeFlatName = "eeFlat" + std::to_string(iflat);
860 calculationResult[eeFlatName] =
861 thetaFlatten >= flatBoundaries[iflat] and thetaFlatten < flatBoundaries[iflat + 1];
862 if (calculationResult[eeFlatName]) {
863 calculationResult["selectee"] = 1;
864 }
865 }
866 }
867
868 // Single tag pi0 / eta dedicated lines
869 if (calculationResult["nTrkLoose"] == 1 and calculationResult["maximumPCMS"] > 0.8 and selectedClusters.size() >= 2) {
870
871 decltype(selectedClusters) selectedSingleTagClusters(selectedClusters.size());
872 auto lastItem = std::copy_if(selectedClusters.begin(), selectedClusters.end(), selectedSingleTagClusters.begin(),
873 [](auto & cluster) {
874 const bool isNeutralCluster = not cluster.isTrack;
875 const bool hasEnoughEnergy = cluster.energyLab > 0.1;
876 const double clusterThetaLab = cluster.p4Lab.Theta() * TMath::RadToDeg();
877 const bool isInAcceptance = 17 < clusterThetaLab and clusterThetaLab < 150.;
878 return isNeutralCluster and hasEnoughEnergy and isInAcceptance;
879 });
880 selectedSingleTagClusters.resize(std::distance(selectedSingleTagClusters.begin(), lastItem));
881
882 if (selectedSingleTagClusters.size() >= 2) { // One track and at least two clusters are found
883
884 const auto& track = maximumPtTracks.at(-1) ? *maximumPtTracks.at(-1) : *maximumPtTracks.at(1);
885 // in real signal, the pi0 decay daughters are always the two most-energetic neutral clusters.
886 const auto& firstCluster = selectedSingleTagClusters[0];
887 const auto& secondCluster = selectedSingleTagClusters[1];
888
889 const ROOT::Math::PxPyPzEVector trackP4CMS = track.p4CMS;
890 const ROOT::Math::PxPyPzEVector pi0P4CMS = firstCluster.p4CMS + secondCluster.p4CMS;
891
892 const bool passPi0ECMS = pi0P4CMS.E() > 1. and pi0P4CMS.E() < 0.525 * p4ofCOM.M();
893 const double thetaSumCMS = (pi0P4CMS.Theta() + trackP4CMS.Theta()) * TMath::RadToDeg();
894 const bool passThetaSum = thetaSumCMS < 170. or thetaSumCMS > 190.;
895
896 double dphiCMS = std::abs(trackP4CMS.Phi() - pi0P4CMS.Phi()) * TMath::RadToDeg();
897 if (dphiCMS > 180) {
898 dphiCMS = 360 - dphiCMS;
899 }
900 const bool passdPhi = dphiCMS > 160.;
901
902 if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() < 0.7) {
903 calculationResult["singleTagLowMass"] = 1;
904 } else if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() > 0.7) {
905 calculationResult["singleTagHighMass"] = 1;
906 }
907 }
908 }
909
910 //..Updated to use new track definitions
911 if (calculationResult["nTrkLooseB"] == 1 and calculationResult["maximumPCMSB"] > 0.8 and selectedClusters.size() >= 2) {
912
913 decltype(selectedClusters) selectedSingleTagClusters(selectedClusters.size());
914 auto lastItem = std::copy_if(selectedClusters.begin(), selectedClusters.end(), selectedSingleTagClusters.begin(),
915 [](auto & cluster) {
916 const bool isNeutralCluster = not cluster.isTrack;
917 const bool hasEnoughEnergy = cluster.energyLab > 0.1;
918 const double clusterThetaLab = cluster.p4Lab.Theta() * TMath::RadToDeg();
919 const bool isInAcceptance = 17 < clusterThetaLab and clusterThetaLab < 150.;
920 return isNeutralCluster and hasEnoughEnergy and isInAcceptance;
921 });
922 selectedSingleTagClusters.resize(std::distance(selectedSingleTagClusters.begin(), lastItem));
923
924 if (selectedSingleTagClusters.size() >= 2) { // One track and at least two clusters are found
925
926 const auto& track = maximumPtTracks.at(-1) ? *maximumPtTracks.at(-1) : *maximumPtTracks.at(1);
927 // in real signal, the pi0 decay daughters are always the two most-energetic neutral clusters.
928 const auto& firstCluster = selectedSingleTagClusters[0];
929 const auto& secondCluster = selectedSingleTagClusters[1];
930
931 const ROOT::Math::PxPyPzEVector trackP4CMS = track.p4CMS;
932 const ROOT::Math::PxPyPzEVector pi0P4CMS = firstCluster.p4CMS + secondCluster.p4CMS;
933
934 const bool passPi0ECMS = pi0P4CMS.E() > 1. and pi0P4CMS.E() < 0.525 * p4ofCOM.M();
935 const double thetaSumCMS = (pi0P4CMS.Theta() + trackP4CMS.Theta()) * TMath::RadToDeg();
936 const bool passThetaSum = thetaSumCMS < 170. or thetaSumCMS > 190.;
937
938 double dphiCMS = std::abs(trackP4CMS.Phi() - pi0P4CMS.Phi()) * TMath::RadToDeg();
939 if (dphiCMS > 180) {
940 dphiCMS = 360 - dphiCMS;
941 }
942 const bool passdPhi = dphiCMS > 160.;
943
944 if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() < 0.7) {
945 calculationResult["singleTagLowMassB"] = 1;
946 } else if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() > 0.7) {
947 calculationResult["singleTagHighMassB"] = 1;
948 }
949 }
950 }
951
952 //..Cosmic selection. Need exactly two tracks.
953 if (m_tracks.getEntries() == 2) {
954
955 const auto negTrack = maximumPtTracksWithoutZCut.at(-1);
956 const auto posTrack = maximumPtTracksWithoutZCut.at(1);
957
958 if (negTrack and posTrack) {
959
960 const double maxNegpT = negTrack->pT;
961 const double maxPospT = posTrack->pT;
962 const double maxClusterENeg = negTrack->clusterEnergySumLab;
963 const double maxClusterEPos = posTrack->clusterEnergySumLab;
964
965 const ROOT::Math::PxPyPzEVector& momentumLabNeg(negTrack->p4Lab);
966 const ROOT::Math::PxPyPzEVector& momentumLabPos(posTrack->p4Lab);
967
968 const double& z0Neg = negTrack->track->getTrackFitResultWithClosestMass(Const::pion)->getZ0();
969 const double& d0Neg = negTrack->track->getTrackFitResultWithClosestMass(Const::pion)->getD0();
970 const double& z0Pos = posTrack->track->getTrackFitResultWithClosestMass(Const::pion)->getZ0();
971 const double& d0Pos = posTrack->track->getTrackFitResultWithClosestMass(Const::pion)->getD0();
972
973 // Select cosmic using these tracks
974 const bool goodMagneticRegion = (z0Neg<m_goodMagneticRegionZ0 or std::abs(d0Neg)>m_goodMagneticRegionD0)
975 and (z0Pos<m_goodMagneticRegionZ0 or std::abs(d0Pos)>m_goodMagneticRegionD0);
976 if (maxNegpT > m_cosmicMinPt and maxPospT > m_cosmicMinPt and maxClusterENeg < m_cosmicMaxClusterEnergy
977 and maxClusterEPos < m_cosmicMaxClusterEnergy and goodMagneticRegion) {
978 double dphiLab = std::abs(momentumLabNeg.Phi() - momentumLabPos.Phi()) * TMath::RadToDeg();
979 if (dphiLab > 180) {
980 dphiLab = 360 - dphiLab;
981 }
982
983 const double thetaSumLab = (momentumLabNeg.Theta() + momentumLabPos.Theta()) * TMath::RadToDeg();
984
985 constexpr double phiBackToBackTolerance = 2.;
986 constexpr double thetaBackToBackTolerance = 2.;
987 if ((180 - dphiLab) < phiBackToBackTolerance and std::abs(180 - thetaSumLab) < thetaBackToBackTolerance) {
988 calculationResult["cosmic"] = 1;
989 }
990 }
991 }
992 }
993
994 //..Displaced vertex
995 // See https://indico.belle2.org/event/8973/ and references therein
996 if (maximumPtTracksWithoutZCut.at(-1) and maximumPtTracksWithoutZCut.at(1)) {
997
998 //..Make particles of the two highest pt tracks (without IP cut)
999 const auto nTrack = maximumPtTracksWithoutZCut.at(-1)->track;
1000 Particle* nParticle = new Particle(nTrack, Const::pion);
1001
1002 const auto pTrack = maximumPtTracksWithoutZCut.at(1)->track;
1003 Particle* pParticle = new Particle(pTrack, Const::pion);
1004
1005 //..Make a vertex of these two
1007 vertexFit->addParticle(nParticle);
1008 vertexFit->addParticle(pParticle);
1009 vertexFit->doFit();
1010
1011 //..Vertex properties
1012 const double chisq = vertexFit->getCHIsq();
1013 const int ndf = vertexFit->getNDF();
1014 const double vertexProb = TMath::Prob(chisq, ndf);
1015 const auto vertexLocation = vertexFit->getVertex();
1016 const double vertexXY = vertexLocation.perp();
1017 const double vertexTheta = vertexLocation.theta() * TMath::RadToDeg();
1018
1019 //..Angular differaence of two tracks to reject cosmics
1020 // Tolerance could be reduced from 10 deg to 2 deg if needed for physics reasons,
1021 // for a 5% increase in the rate of selected displaced vertex triggers.
1022 // See https://gitlab.desy.de/belle2/software/basf2/-/merge_requests/1867
1023 const ROOT::Math::PxPyPzEVector& momentumLabNeg(maximumPtTracksWithoutZCut.at(-1)->p4Lab);
1024 const ROOT::Math::PxPyPzEVector& momentumLabPos(maximumPtTracksWithoutZCut.at(1)->p4Lab);
1025 const double thetaSumLab = (momentumLabNeg.Theta() + momentumLabPos.Theta()) * TMath::RadToDeg();
1026 double dPhiLab = std::abs(momentumLabNeg.Phi() - momentumLabPos.Phi()) * TMath::RadToDeg();
1027 if (dPhiLab > 180) {
1028 dPhiLab = 360 - dPhiLab;
1029 }
1030 const double backToBackTolerance = 10.; // degrees
1031 const bool backToBackLab = std::abs(thetaSumLab - 180.) < backToBackTolerance and std::abs(dPhiLab - 180.) < backToBackTolerance;
1032
1033 //..Select a displaced vertex
1034 const double minProbChiVertex = 0.005; // minimum probability chi square. Many backgrounds have bad chi sq.
1035 const double minXYVertex = 3.; // minimum xy of vertex (cm). Should pass track filters below this.
1036 const double maxXYVertex = 60.; // maximum xy. Insufficient CDC for good reconstruction beyond this.
1037 const double minThetaVertex = 30.; // Large Bhabha background at low angles.
1038 const double maxThetaVertex = 120.; // Large Bhabha background at low angles.
1039 if (vertexProb > minProbChiVertex and vertexXY > minXYVertex and vertexXY < maxXYVertex and vertexTheta > minThetaVertex
1040 and vertexTheta < maxThetaVertex
1041 and not backToBackLab) {calculationResult["displacedVertex"] = 1;}
1042
1043 //..Clean up
1044 delete nParticle;
1045 delete pParticle;
1046 delete vertexFit;
1047
1048 }
1049
1050}
Class to provide momentum-related information from ECLClusters.
const ROOT::Math::PxPyPzEVector Get4MomentumFromCluster(const ECLCluster *cluster, ECLCluster::EHypothesisBit hypo)
Returns four momentum vector.
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)
Definition ECLCluster.h:41
unsigned short getNHits() const
Get the total Number of CDC hits in the fit.
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:76
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
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.
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
StoreArray< CDCTriggerUnpacker::NNBitStream > m_bitsNN
Store Object with the trigger NN bits.
FilterCalculator()
Set the default names for the store object particle lists.
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.
short getChargeSign() const
Return track charge (1 or -1).
ROOT::Math::PxPyPzEVector get4Momentum() const
Getter for the 4Momentum at the closest approach of the track in the r/phi projection.
double getTransverseMomentum() const
Getter for the absolute value of the transverse momentum at the perigee.
double getZ0() const
Getter for z0.
HitPatternCDC getHitPatternCDC() const
Getter for the hit pattern in the CDC;.
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.
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.
Abstract base class for different kinds of events.
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