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:

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

Detailed Description

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

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

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

Definition at line 30 of file FilterCalculator.h.

Constructor & Destructor Documentation

FilterCalculator()

FilterCalculator

(

)

Set the default names for the store object particle lists.

Definition at line 73 of file FilterCalculator.cc.

                                   : m_bitsNN("CDCTriggerNNBits")
{
 
}

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.

    {
      for (auto& identifierWithValue : m_calculationResult) {
        const std::string& identifier = identifierWithValue.first;
        const double value = identifierWithValue.second;
 
        storeObject->append(prefix + "_" + identifier, value);
      }
    }

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

< number of "C" tracks

< maximum pcms of "C" tracks GeV/c

< pCms of maximum p negatively charged track GeV/c

< clusterE lab of max p negatively charged track GeV

< pCms of maximum p positively charged track GeV/c

< clusterE lab of max p positively charged track GeV

< delta phi cms of max p positive and negative tracks deg

Implements SoftwareTriggerCalculation.

Definition at line 86 of file FilterCalculator.cc.

{
  calculationResult["nTrkLoose"] = 0; 
  calculationResult["nTrkTight"] = 0; 
  calculationResult["ee2leg"] = 0; 
  calculationResult["nEmedium"] = 0; 
  calculationResult["nElow"] = 0; 
  calculationResult["nEhigh"] = 0; 
  calculationResult["nE180Lab"] = 0; 
  calculationResult["nE300Lab"] = 0; 
  calculationResult["nE500Lab"] = 0; 
  calculationResult["nE2000CMS"] = 0; 
  calculationResult["nE4000CMS"] = 0; 
  calculationResult["nE250Lab"] = 0; 
  calculationResult["nMaxEPhotonAcc"] = 0; 
  calculationResult["dphiCmsClust"] = NAN; 
 
  calculationResult["netChargeLoose"] = 0; 
  calculationResult["maximumPCMS"] = NAN; 
  calculationResult["maximumPLab"] = NAN; 
  calculationResult["eexx"] = 0;
  calculationResult["ee1leg1trk"] = 0; 
  calculationResult["nEhighLowAng"] = 0; 
  calculationResult["nEsingleClust"] = 0; 
  calculationResult["nEsinglePhotonBarrel"] = 0; 
  calculationResult["nEsinglePhotonExtendedBarrel"] = 0; 
  calculationResult["nEsinglePhotonEndcap"] = 0; 
  calculationResult["nEsingleElectronBarrel"] = 0; 
  calculationResult["nEsingleElectronExtendedBarrel"] = 0; 
  calculationResult["nReducedEsinglePhotonReducedBarrel"] = 0; 
  calculationResult["nVetoClust"] = 0; 
  calculationResult["chrgClust2GeV"] = 0; 
  calculationResult["neutClust045GeVAcc"] = 0; 
  calculationResult["neutClust045GeVBarrel"] = 0; 
  calculationResult["singleTagLowMass"] = 0;
  calculationResult["singleTagHighMass"] = 0;
  calculationResult["n2GeVNeutBarrel"] = 0;
  calculationResult["n2GeVNeutEndcap"] = 0;
  calculationResult["n2GeVChrg"] = 0;
  calculationResult["n2GeVPhotonBarrel"] = 0;
  calculationResult["n2GeVPhotonEndcap"] = 0;
  calculationResult["ee1leg"] = 0; 
  calculationResult["ee1leg1clst"] = 0; 
  calculationResult["ee1leg1e"] = 0; 
  calculationResult["ee2clst"] = 0; 
  calculationResult["gg2clst"] = 0; 
  calculationResult["eeee"] = 0;
  calculationResult["eemm"] = 0;
  calculationResult["eexxSelect"] = 0;
  calculationResult["radBhabha"] = 0;
  calculationResult["eeBrem"] = 0;
  calculationResult["isrRadBhabha"] = 0;
  calculationResult["muonPairECL"] = 0;
  calculationResult["ggHighPt"] = 0;
  calculationResult["selectee1leg1trk"] = 0;
  calculationResult["selectee1leg1clst"] = 0;
  calculationResult["selectee"] = 0;
  calculationResult["ggBarrelVL"] = 0;
  calculationResult["ggBarrelLoose"] = 0;
  calculationResult["ggBarrelTight"] = 0;
  calculationResult["ggEndcapVL"] = 0;
  calculationResult["ggEndcapLoose"] = 0;
  calculationResult["ggEndcapTight"] = 0;
  calculationResult["muonPairV"] = 0;
  calculationResult["selectmumu"] = 0;
  calculationResult["singleMuon"] = 0;
  calculationResult["cosmic"] = 0;
  calculationResult["displacedVertex"] = 0;
  calculationResult["eeFlat0"] = 0;
  calculationResult["eeFlat1"] = 0;
  calculationResult["eeFlat2"] = 0;
  calculationResult["eeFlat3"] = 0;
  calculationResult["eeFlat4"] = 0;
  calculationResult["eeFlat5"] = 0;
  calculationResult["eeFlat6"] = 0;
  calculationResult["eeFlat7"] = 0;
  calculationResult["eeFlat8"] = 0;
  calculationResult["eeOneClust"] = 0;
 
  //..New track definitions
  calculationResult["nTrkLooseB"] = 0; 
  calculationResult["nTrkTightB"] = 0; 
  calculationResult["maximumPCMSB"] = NAN; 
  calculationResult["netChargeLooseB"] = 0; 
  calculationResult["muonPairVB"] = 0; 
  calculationResult["eeBremB"] = 0; 
  calculationResult["singleTagLowMassB"] = 0; 
  calculationResult["singleTagHighMassB"] = 0; 
  calculationResult["radBhabhaB"] = 0; 
 
  //..Tracks without a requirement on nCDCHits
  calculationResult["nTrackC"] = 0; 
  calculationResult["maximumPCMSC"] = NAN; 
  calculationResult["pCmsNegC"] = NAN; 
  calculationResult["clusterENegC"] = NAN; 
  calculationResult["pCmsPosC"] = NAN; 
  calculationResult["clusterEPosC"] = NAN; 
  calculationResult["dPhiCmsC"] = NAN; 
 
 
  // Passed on L1 information
  if (m_l1Trigger.isValid()) {
    calculationResult["l1_trigger_random"] = m_l1Trigger->getTimType() == TRGSummary::ETimingType::TTYP_RAND;
    calculationResult["l1_trigger_delayed_bhabha"] = m_l1Trigger->getTimType() == TRGSummary::ETimingType::TTYP_DPHY;
    calculationResult["l1_trigger_poisson"] = m_l1Trigger->getTimType() == TRGSummary::ETimingType::TTYP_POIS;
    bool bha3d;
    try {
      bha3d = m_l1Trigger->testPsnm("bha3d");
    } catch (const std::exception&) {
      bha3d = false;
    }
    bool bhapurPsnm;
    try {
      bhapurPsnm = m_l1Trigger->testPsnm("bhapur");
    } catch (const std::exception&) {
      bhapurPsnm = false;
    }
    bool bhapurFtdl;
    try {
      bhapurFtdl = m_l1Trigger->testFtdl("bhapur");
    } catch (const std::exception&) {
      bhapurFtdl = false;
    }
    bool lml1;
    try {
      lml1 = m_l1Trigger->testPsnm("lml1");
    } catch (const std::exception&) {
      lml1 = false;
    }
    bool l1_bit_f;
    try {
      l1_bit_f = m_l1Trigger->testPsnm("fpre");
    } catch (const std::exception&) {
      l1_bit_f = false;
    }
    calculationResult["bha3d"] = bha3d;
    calculationResult["bhapur"] = bhapurPsnm;
    calculationResult["bhapur_lml1"] = lml1 and bhapurFtdl;
    calculationResult["l1_bit_f"] = l1_bit_f;
  } else {
    calculationResult["l1_trigger_random"] = 1; // save every event if no L1 trigger info
    calculationResult["l1_trigger_delayed_bhabha"] = 0;
    calculationResult["l1_trigger_poisson"] = 0;
    calculationResult["bha3d"] = 0;
    calculationResult["bhapur"] = 0;
    calculationResult["bhapur_lml1"] = 0;
    calculationResult["l1_bit_f"] = 0;
  }
 
  // Every 256th event has CDC NN trigger information
  calculationResult["l1_trg_NN_info"] = 0;
  if (m_bitsNN.isValid() and m_bitsNN.getEntries() > 0) {calculationResult["l1_trg_NN_info"] = 1;}
 
  calculationResult["true"] = 1;
  calculationResult["false"] = 0;
 
  // Some utilities
  // TODO: into constructor
  ClusterUtils cUtil;
  const ROOT::Math::XYZVector clustervertex = cUtil.GetIPPosition();
  PCmsLabTransform boostrotate;
  ROOT::Math::PxPyPzEVector p4ofCOM;
  p4ofCOM.SetPxPyPzE(0, 0, 0, boostrotate.getCMSEnergy());
 
  // Pointers to the two tracks with the maximum pt
  std::map<short, std::optional<MaximumPtTrack>> maximumPtTracks = {
    { -1, {}}, {1,  {}}
  };
 
  // Pointer to the two tracks with maximum pt without a cut applied on z0 (used for cosmic trigger)
  std::map<short, std::optional<MaximumPtTrack>> maximumPtTracksWithoutZCut = {
    { -1, {}}, {1, {}}
  };
 
  //..Pointers to the two "C" tracks with maximum pCms (no requirement on nCDC hits)
  std::map<short, std::optional<MaximumPtTrack>> maximumPCmsTracksC = {
    { -1, {}}, {1, {}}
  };
 
  //------------------------------------------------------------
  // --- Track variables -- //
  for (const Track& track : m_tracks) {
    const TrackFitResult* trackFitResult = track.getTrackFitResultWithClosestMass(Const::pion);
    if (not trackFitResult) {
      // Huh? No track fit result for pion? Ok, lets skip this track
      continue;
    }
 
    //..Track properties. Only use tracks with defined charge.
    const short charge = trackFitResult->getChargeSign();
    if (charge == 0) {continue;}
    const double z0 = trackFitResult->getZ0();
    const ROOT::Math::PxPyPzEVector& momentumLab = trackFitResult->get4Momentum();
    const ROOT::Math::PxPyPzEVector momentumCMS = boostrotate.rotateLabToCms() * momentumLab;
    const double pCMS = momentumCMS.P();
    const double pLab = momentumLab.P();
    const double trackTime = track.getTrackTime();
    const double pT = trackFitResult->getTransverseMomentum();
    const int nCDCHits = (int)(0.5 + trackFitResult->getHitPatternCDC().getNHits());
    double clusterELab = 0.;
    for (auto& cluster : track.getRelationsTo<ECLCluster>()) {
      if (cluster.hasHypothesis(Belle2::ECLCluster::EHypothesisBit::c_nPhotons)) {
        clusterELab += cluster.getEnergy(Belle2::ECLCluster::EHypothesisBit::c_nPhotons);
      }
    }
    const double clusterECMS = clusterELab * pCMS / pLab;
 
 
    //------------------------------------------------------------
    //..New list without a requirement on number of CDC hits
    bool goodTrackCTime = true;
    if (std::abs(trackTime) > 15.) {goodTrackCTime = false;}
    if (std::abs(z0) < 1. and goodTrackCTime and pCMS > 0.2) {
      calculationResult["nTrackC"] += 1;
      if (std::isnan(calculationResult["maximumPCMSC"]) or pCMS > calculationResult["maximumPCMSC"]) {
        calculationResult["maximumPCMSC"] = pCMS;
      }
    }
 
    //..Store maximum pCms positive and negative tracks
    const auto& currentMaximumC = maximumPCmsTracksC.at(charge);
    if (not currentMaximumC or pCMS > currentMaximumC->pCMS) {
      MaximumPtTrack newMaximum;
      newMaximum.pT = pT;
      newMaximum.track = &track;
      newMaximum.pCMS = pCMS;
      newMaximum.pLab = pLab;
      newMaximum.p4CMS = momentumCMS;
      newMaximum.p4Lab = momentumLab;
      newMaximum.clusterEnergySumLab = clusterELab;
      newMaximum.clusterEnergySumCMS = clusterECMS;
      maximumPCmsTracksC[charge] = newMaximum;
    }
 
    //------------------------------------------------------------
    //..Now original loose and tight tracks, which require nCDC hits
    // TODO: Temporary cut on number of CDC hits
    if (nCDCHits == 0) {
      continue;
    }
 
    // Count tight tracks
    if (std::abs(z0) < m_tightTrkZ0) {
      calculationResult["nTrkTight"] += 1;
    }
 
 
    //------------------------------------------------------------
    // Find the maximum pt negative [0] and positive [1] tracks without z0 cut
    const auto& currentMaximum = maximumPtTracksWithoutZCut.at(charge);
    if (not currentMaximum or pT > currentMaximum->pT) {
      MaximumPtTrack newMaximum;
      newMaximum.pT = pT;
      newMaximum.track = &track;
      newMaximum.pCMS = pCMS;
      newMaximum.pLab = pLab;
      newMaximum.p4CMS = momentumCMS;
      newMaximum.p4Lab = momentumLab;
      newMaximum.clusterEnergySumLab = clusterELab;
      newMaximum.clusterEnergySumCMS = clusterECMS;
      maximumPtTracksWithoutZCut[charge] = newMaximum;
    }
 
 
    //------------------------------------------------------------
    // Loose tracks
    if (std::abs(z0) < m_looseTrkZ0) {
      calculationResult["nTrkLoose"] += 1;
      calculationResult["netChargeLoose"] += charge;
 
      if (std::isnan(calculationResult["maximumPCMS"]) or pCMS > calculationResult["maximumPCMS"]) {
        calculationResult["maximumPCMS"] = pCMS;
      }
 
      if (std::isnan(calculationResult["maximumPLab"]) or pLab > calculationResult["maximumPLab"]) {
        calculationResult["maximumPLab"] = pLab;
      }
 
      // Find the maximum pt negative [0] and positive [1] tracks
      const double pTLoose = trackFitResult->getTransverseMomentum();
      const auto& currentMaximumLoose = maximumPtTracks.at(charge);
      if (not currentMaximumLoose or pTLoose > currentMaximumLoose->pT) {
        MaximumPtTrack newMaximum;
        newMaximum.pT = pTLoose;
        newMaximum.track = &track;
        newMaximum.pCMS = pCMS;
        newMaximum.pLab = momentumLab.P();
        newMaximum.p4CMS = momentumCMS;
        newMaximum.p4Lab = momentumLab;
        newMaximum.clusterEnergySumLab = clusterELab;
        newMaximum.clusterEnergySumCMS = clusterECMS;
        maximumPtTracks[charge] = newMaximum;
      }
    }
 
    //..New tighter track definitions.
    //  We can use pCMS and charge from above, because every looseB track is also a loose track
    if (nCDCHits >= 5) {
      if (std::abs(z0) < 1.) {calculationResult["nTrkTightB"] += 1;}
      if (std::abs(z0) < 5.) {
        calculationResult["nTrkLooseB"] += 1;
        calculationResult["netChargeLooseB"] += charge;
        if (std::isnan(calculationResult["maximumPCMSB"]) or pCMS > calculationResult["maximumPCMSB"]) {
          calculationResult["maximumPCMSB"] = pCMS;
 
        }
      }
    }
 
  } // End loop over tracks
 
 
  //------------------------------------------------------------
  //..Store the "C" track calculations
  for (short charge : { -1, 1}) {
    auto& maximumPcmsTrackC = maximumPCmsTracksC.at(charge);
    if (not maximumPcmsTrackC) {
      continue;
    }
 
    //..Calculations for each track
    if (charge == -1) {
      calculationResult["pCmsNegC"] = maximumPcmsTrackC->pCMS;
      calculationResult["clusterENegC"] = maximumPcmsTrackC->clusterEnergySumLab;
    } else {
      calculationResult["pCmsPosC"] = maximumPcmsTrackC->pCMS;
      calculationResult["clusterEPosC"] = maximumPcmsTrackC->clusterEnergySumLab;
    }
  }
 
  //..Calculations using both highest-pCms C tracks
  if (maximumPCmsTracksC.at(-1) and maximumPCmsTracksC.at(1)) {
    const MaximumPtTrack& negativeTrack = *maximumPCmsTracksC.at(-1);
    const MaximumPtTrack& positiveTrack = *maximumPCmsTracksC.at(1);
 
    calculationResult["dPhiCmsC"] = std::abs(ROOT::Math::VectorUtil::DeltaPhi(negativeTrack.p4CMS,
                                             positiveTrack.p4CMS)) * TMath::RadToDeg();
  }
 
 
  //------------------------------------------------------------
  // -- Cluster variables -- //
  std::vector<SelectedECLCluster> selectedClusters;
  for (const ECLCluster& cluster : m_eclClusters) {
    // Use the photon hypothesis, and only look at clusters with times < 200 ns
    const double time = cluster.getTime();
    if (not cluster.hasHypothesis(Belle2::ECLCluster::EHypothesisBit::c_nPhotons) or
        std::abs(time) > 200 or
        cluster.getEnergy(Belle2::ECLCluster::EHypothesisBit::c_nPhotons) < 0.1) {
      continue;
    }
 
    const double dt99 = cluster.getDeltaTime99();
 
    // Store cluster information
    SelectedECLCluster selectedCluster;
    selectedCluster.p4Lab = cUtil.Get4MomentumFromCluster(&cluster, clustervertex,
                                                          Belle2::ECLCluster::EHypothesisBit::c_nPhotons);
    selectedCluster.p4CMS = boostrotate.rotateLabToCms() * selectedCluster.p4Lab; // was clustp4COM
    selectedCluster.cluster = &cluster;
    selectedCluster.clusterTime = time / dt99; // was clustT
    selectedCluster.energyLab = selectedCluster.p4Lab.E();
    selectedCluster.energyCMS = selectedCluster.p4CMS.E(); // was first of ECOMPair
    selectedCluster.isTrack = cluster.isTrack(); // was tempCharge
 
    selectedClusters.push_back(selectedCluster);
 
    if (selectedCluster.energyCMS > m_E2min) {
      calculationResult["nElow"] += 1;
    }
    if (selectedCluster.energyCMS > m_E0min) {
      calculationResult["nEmedium"] += 1;
    }
    if (selectedCluster.energyCMS > m_Ehigh) {
      calculationResult["nEhigh"] += 1;
    }
 
    if (selectedCluster.energyCMS > m_E0min and std::abs(selectedCluster.clusterTime) < 10) {
      calculationResult["nVetoClust"] += 1;
    }
 
 
    //..Single cluster trigger objects use charge, Zernike moment, and thetaLab
    const double thetaLab = selectedCluster.p4Lab.Theta() * TMath::RadToDeg();
    const double zmva = cluster.getZernikeMVA();
    const bool photon = zmva > 0.5 and not selectedCluster.isTrack;
    const bool electron = zmva > 0.5 and selectedCluster.isTrack;
 
    //..For 1 track radiative Bhabha control sample
    if (selectedCluster.energyCMS > 2. and selectedCluster.isTrack) {
      calculationResult["chrgClust2GeV"] += 1;
    }
    if (selectedCluster.energyCMS > 0.45 and not selectedCluster.isTrack) {
      const bool isInAcceptance = 17. < thetaLab and thetaLab < 150.;
      if (isInAcceptance) {calculationResult["neutClust045GeVAcc"] += 1;}
      const bool isInBarrel = 30. < thetaLab and thetaLab < 130.;
      if (isInBarrel) {calculationResult["neutClust045GeVBarrel"] += 1;}
    }
 
 
    // improved 3 cluster (4 cluster) trigger
    const bool notInHighBackgroundEndcapRegion = 18.5 < thetaLab and thetaLab < 139.3;
    if (selectedCluster.energyLab > m_EminLab and notInHighBackgroundEndcapRegion) {
      calculationResult["nE180Lab"] += 1;
    }
 
    if (selectedCluster.energyLab > m_EminLab4Cluster and notInHighBackgroundEndcapRegion) {
      calculationResult["nE300Lab"] += 1;
    }
 
    if (selectedCluster.energyLab > m_EminLab3Cluster and notInHighBackgroundEndcapRegion) {
      calculationResult["nE500Lab"] += 1;
    }
 
    if (selectedCluster.energyCMS > m_Ehigh and notInHighBackgroundEndcapRegion) {
      calculationResult["nE2000CMS"] += 1;
    }
 
    //..For two-photon fusion ALP trigger
    if (selectedCluster.energyLab > 0.25) {
      calculationResult["nE250Lab"] += 1;
    }
    if (selectedCluster.energyCMS > 4.) {
      calculationResult["nE4000CMS"] += 1;
    }
 
    //..Single cluster triggers
    if (selectedCluster.energyCMS > m_EsinglePhoton) {
      calculationResult["nEsingleClust"] += 1;
 
      const bool barrelRegion = thetaLab > 45 and thetaLab < 115;
      const bool extendedBarrelRegion = thetaLab > 30 and thetaLab < 130;
      const bool endcapRegion = (thetaLab > 22 and thetaLab < 45) or (thetaLab > 115 and thetaLab < 145);
 
      if (photon and barrelRegion) {
        calculationResult["nEsinglePhotonBarrel"] += 1;
      }
 
      if (photon and extendedBarrelRegion) {
        calculationResult["nEsinglePhotonExtendedBarrel"] += 1;
      }
 
      if (electron and barrelRegion) {
        calculationResult["nEsingleElectronBarrel"] += 1;
      }
 
      if (electron and extendedBarrelRegion) {
        calculationResult["nEsingleElectronExtendedBarrel"] += 1;
      }
 
      if (photon and endcapRegion) {
        calculationResult["nEsinglePhotonEndcap"] += 1;
      }
    }
 
    if (selectedCluster.energyCMS > m_reducedEsinglePhoton) {
      const bool reducedBarrelRegion = thetaLab > 44 and thetaLab < 98;
 
      if (photon and reducedBarrelRegion) {
        calculationResult["nReducedEsinglePhotonReducedBarrel"] += 1;
      }
    }
 
    if (selectedCluster.energyCMS > m_Ehigh) {
      // TODO: different definition!
      const bool barrelRegion = thetaLab > 32 and thetaLab < 130;
      const bool endcapRegion = (thetaLab > 22 and thetaLab < 32) or (thetaLab > 130 and thetaLab < 145);
      const bool lowAngleRegion = thetaLab < 22 or thetaLab > 145;
 
      if (not selectedCluster.isTrack and barrelRegion) {
        calculationResult["n2GeVNeutBarrel"] += 1;
      }
      if (not selectedCluster.isTrack and endcapRegion) {
        calculationResult["n2GeVNeutEndcap"] += 1;
      }
      if (selectedCluster.isTrack and not lowAngleRegion) {
        calculationResult["n2GeVChrg"] += 1;
      }
      if (lowAngleRegion) {
        calculationResult["nEhighLowAng"] += 1;
      }
      if (photon and barrelRegion) {
        calculationResult["n2GeVPhotonBarrel"] += 1;
      }
      if (photon and endcapRegion) {
        calculationResult["n2GeVPhotonEndcap"] += 1;
      }
    }
  } // end of loop over clusters
 
  //..Order clusters by CMS energy
  std::sort(selectedClusters.begin(), selectedClusters.end(), [](const auto & lhs, const auto & rhs) {
    return lhs.energyCMS > rhs.energyCMS;
  });
 
  // -- Bhabha and two-photon lepton preparation -- //
  for (short charge : { -1, 1}) {
    auto& maximumPtTrack = maximumPtTracks.at(charge);
    if (not maximumPtTrack) {
      continue;
    }
 
    // Single leg Bhabha selections
    if (maximumPtTrack->clusterEnergySumCMS > 4.5) {
      calculationResult["ee1leg"] = 1;
    }
 
    // single muon
    if (maximumPtTrack->pCMS > 3 and maximumPtTrack->clusterEnergySumLab > 0. and maximumPtTrack->clusterEnergySumLab < 1.) {
      calculationResult["singleMuon"] = 1;
    }
  }
 
  // Bhabha selections using two tracks
  if (maximumPtTracks.at(-1) and maximumPtTracks.at(1)) {
    const MaximumPtTrack& negativeTrack = *maximumPtTracks.at(-1);
    const MaximumPtTrack& positiveTrack = *maximumPtTracks.at(1);
 
    double dphi = std::abs(negativeTrack.p4CMS.Phi() - positiveTrack.p4CMS.Phi()) * TMath::RadToDeg();
    if (dphi > 180) {
      dphi = 360 - dphi;
    }
 
    const double negativeClusterSum = negativeTrack.clusterEnergySumCMS;
    const double negativeClusterSumLab = negativeTrack.clusterEnergySumLab;
    const double negativeP = negativeTrack.pCMS;
    const double positiveClusterSum = positiveTrack.clusterEnergySumCMS;
    const double positiveClusterSumLab = positiveTrack.clusterEnergySumLab;
    const double positiveP = positiveTrack.pCMS;
 
    // two leg Bhabha
    const double thetaSum = (negativeTrack.p4CMS.Theta() + positiveTrack.p4CMS.Theta()) * TMath::RadToDeg();
    const double dthetaSum = std::abs(thetaSum - 180);
    const auto back2back = dphi > 175 and dthetaSum < 15;
    if (back2back and negativeClusterSum > 3 and positiveClusterSum > 3 and
        (negativeClusterSum > 4.5 or positiveClusterSum > 4.5)) {
      calculationResult["ee2leg"] = 1;
    }
 
    // one leg one track Bhabha
    if (back2back and ((negativeClusterSum > 4.5 and positiveP > 3) or (positiveClusterSum > 4.5 and negativeP > 3))) {
      calculationResult["ee1leg1trk"] = 1;
    }
 
 
    // one leg plus one electron
    if ((negativeClusterSum > 4.5 and positiveClusterSum > 0.8 * positiveP) or
        (positiveClusterSum > 4.5 and negativeClusterSum > 0.8 * negativeP)) {
      calculationResult["ee1leg1e"] = 1;
    }
 
    // eeee, eemm, mu mu, and radiative Bhabha selection
    const ROOT::Math::PxPyPzEVector p4Miss = p4ofCOM - negativeTrack.p4CMS - positiveTrack.p4CMS;
    const double pmissTheta = p4Miss.Theta() * TMath::RadToDeg();
    const double pmissp = p4Miss.P();
    const double relMissAngle0 = ROOT::Math::VectorUtil::Angle(negativeTrack.p4CMS, p4Miss) * TMath::RadToDeg();
    const double relMissAngle1 = ROOT::Math::VectorUtil::Angle(positiveTrack.p4CMS, p4Miss) * TMath::RadToDeg();
 
    const bool electronEP = positiveClusterSum > 0.8 * positiveP or negativeClusterSum > 0.8 * negativeP;
    const bool notMuonPair = negativeClusterSumLab > 1 or positiveClusterSumLab > 1;
    const double highp = std::max(negativeP, positiveP);
    const double lowp = std::min(negativeP, positiveP);
    const bool lowEdep = negativeClusterSumLab < 0.5 and positiveClusterSumLab < 0.5;
 
    //..Select two photon fusion lepton pairs
    if (calculationResult["maximumPCMS"] < 2 and dphi > 160 and (pmissTheta < 25. or pmissTheta > 155.)) {
      calculationResult["eexxSelect"] = 1;
      if (electronEP) {
        calculationResult["eeee"] = 1;
      } else {
        calculationResult["eemm"] = 1;
      }
    }
 
    //..Veto two photon fusion lepton pairs
    if ((pmissTheta < 20. or pmissTheta > 160.) and
        ((calculationResult["maximumPCMS"] < 1.2 and dphi > 150.) or
         (calculationResult["maximumPCMS"] < 2. and 175. < dphi))) {
      calculationResult["eexx"] = 1;
    }
 
    //..Radiative Bhabhas with missing momentum in acceptance
    if (negativeP > 1. and pmissTheta > 10. and pmissTheta < 170. and positiveP > 1. and dphi < 170. and pmissp > 1. and electronEP) {
      if (calculationResult["nTrkLoose"] == 2 and calculationResult["nTrkTight"] >= 1) {
        calculationResult["radBhabha"] = 1;
      }
      if (calculationResult["nTrkLooseB"] == 2 and calculationResult["nTrkTightB"] >= 1) {
        calculationResult["radBhabhaB"] = 1;
      }
    }
 
    //..ISR radiative Bhabha with missing momentum down beam pipe
    if (negativeP > 2. and positiveP > 2. and 2 == calculationResult["nTrkLoose"] and
        calculationResult["nTrkTight"] >= 1 and dphi > 175. and
        (pmissTheta < 5. or pmissTheta > 175.) and electronEP) {
      calculationResult["isrRadBhabha"] = 1;
    }
 
    //..Veto Bhabha with hard brem
    if (highp > 4.5 and notMuonPair and pmissp > 1. and (relMissAngle0 < 5. or relMissAngle1 < 5.)) {
      if (calculationResult["nTrkLoose"] == 2) { calculationResult["eeBrem"] = 1;}
      if (calculationResult["nTrkLooseB"] >= 1) { calculationResult["eeBremB"] = 1;}
    }
 
    //..Veto e+e- -> mu+mu-
    if (highp > 4.5 and lowEdep and thetaSum > 175. and thetaSum < 185. and dphi > 175.) {
      if (calculationResult["nTrkLoose"] == 2) {calculationResult["muonPairV"] = 1;}
      if (calculationResult["nTrkLooseB"] == 2) {calculationResult["muonPairVB"] = 1;}
    }
 
    //..Select e+e- -> mu+mu-
    if (highp > 3. and lowp > 2.5 and dphi > 165. and
        ((negativeClusterSumLab > 0. and negativeClusterSumLab < 1.) or
         (positiveClusterSumLab > 0. and positiveClusterSumLab < 1.))) {
      calculationResult["selectmumu"] = 1;
    }
  }
 
  //..Filters and vetoes using two maximum-energy clusters
  if (selectedClusters.size() >= 2) {
    const SelectedECLCluster& firstCluster = selectedClusters[0];
    const SelectedECLCluster& secondCluster = selectedClusters[1];
 
    double dphi = std::abs(firstCluster.p4CMS.Phi() - secondCluster.p4CMS.Phi()) * TMath::RadToDeg();
    if (dphi > 180) {
      dphi = 360 - dphi;
    }
    double thetaSum = (firstCluster.p4CMS.Theta() + secondCluster.p4CMS.Theta()) * TMath::RadToDeg();
    double dthetaSum = std::abs(thetaSum - 180);
 
    //..Quantities for two-photon fusion triggers
    calculationResult["dphiCmsClust"] = dphi;
    for (int ic = 0; ic < 2; ic++) {
      const double thetaLab = selectedClusters[ic].p4Lab.Theta() * TMath::RadToDeg();
      const bool isInAcceptance = 17. < thetaLab and thetaLab < 150.;
      const ECLCluster* cluster = selectedClusters[ic].cluster;
      const double zmva = cluster->getZernikeMVA();
      const bool photon = zmva > 0.5 and not selectedClusters[ic].isTrack;
      if (isInAcceptance and photon) {calculationResult["nMaxEPhotonAcc"] += 1;}
    }
 
    const double firstEnergy = firstCluster.p4CMS.E();
    const double secondEnergy = secondCluster.p4CMS.E();
 
    const bool highEnergetic = firstEnergy > 3 and secondEnergy > 3 and (firstEnergy > 4.5 or secondEnergy > 4.5);
 
    if (dphi > 160 and dphi < 178 and dthetaSum < 15 and highEnergetic) {
      calculationResult["ee2clst"] = 1;
    }
 
    if (dphi > 178 and dthetaSum < 15 and highEnergetic) {
      calculationResult["gg2clst"] = 1;
    }
 
    if ((calculationResult["ee2clst"] == 1 or calculationResult["gg2clst"] == 1) and
        calculationResult["ee1leg"] == 1) {
      calculationResult["ee1leg1clst"] = 1;
    }
 
    const double Elab0 = firstCluster.p4Lab.E();
    const double Elab1 = secondCluster.p4Lab.E();
 
    // gg and mumu accept triggers using ECL
    if (firstEnergy > 2 and secondEnergy > 2) {
      const double thetaLab0 = firstCluster.p4Lab.Theta() * TMath::RadToDeg();
      const double thetaLab1 = secondCluster.p4Lab.Theta() * TMath::RadToDeg();
 
      const bool barrel0 = 32. < thetaLab0 and thetaLab0 < 130.;
      const bool barrel1 = 32. < thetaLab1 and thetaLab1 < 130.;
      const bool oneClustersAbove4 = firstEnergy > 4 or secondEnergy > 4;
      const bool oneIsNeutral = not firstCluster.isTrack or not secondCluster.isTrack;
      const bool bothAreNeutral = not firstCluster.isTrack and not secondCluster.isTrack;
      const bool oneIsBarrel = barrel0 or barrel1;
      const bool dphiCutExtraLoose = dphi > 175;
      const bool dphiCutLoose = dphi > 177;
      const bool dphiCutTight = dphi > 177.5;
 
      if (dphiCutExtraLoose and oneIsNeutral and oneIsBarrel) {
        calculationResult["ggBarrelVL"] = 1;
      }
      if (oneClustersAbove4 and dphiCutLoose and oneIsNeutral and oneIsBarrel) {
        calculationResult["ggBarrelLoose"] = 1;
      }
      if (oneClustersAbove4 and dphiCutTight and bothAreNeutral and oneIsBarrel) {
        calculationResult["ggBarrelTight"] = 1;
      }
      if (dphiCutExtraLoose and oneIsNeutral and not oneIsBarrel) {
        calculationResult["ggEndcapVL"] = 1;
      }
      if (oneClustersAbove4 and dphiCutLoose and oneIsNeutral and not oneIsBarrel) {
        calculationResult["ggEndcapLoose"] = 1;
      }
      if (oneClustersAbove4 and dphiCutTight and bothAreNeutral and not oneIsBarrel) {
        calculationResult["ggEndcapTight"] = 1;
      }
    }
 
    const double minEnergy = std::min(Elab0, Elab1);
    const double maxEnergy = std::max(Elab0, Elab1);
    if (dphi > 155 and thetaSum > 165 and thetaSum < 195 and minEnergy > 0.15 and minEnergy < 0.5 and
        maxEnergy > 0.15 and maxEnergy < 0.5) {
      calculationResult["muonPairECL"] = 1;
    }
 
    //..diPhoton line
    const double thetaLab0 = firstCluster.p4Lab.Theta() * TMath::RadToDeg();
    const double thetaLab1 = secondCluster.p4Lab.Theta() * TMath::RadToDeg();
    const bool inHieRegion0 = thetaLab0 > 26. and thetaLab0 < 130.;
    const bool inHieRegion1 = thetaLab1 > 26. and thetaLab1 < 130.;
    const bool firstIsNeutral = not firstCluster.isTrack;
    const bool secondIsNeutral = not secondCluster.isTrack;
 
    if (secondEnergy > 0.3 and inHieRegion0 and inHieRegion1 and firstIsNeutral and secondIsNeutral) {
      const ROOT::Math::PxPyPzEVector ggP4CMS = firstCluster.p4CMS + secondCluster.p4CMS;
      if (ggP4CMS.pt() > 1.) {calculationResult["ggHighPt"] = 1;}
    }
 
  } // end of two-clusters lines
 
  // Bhabha accept triggers.
  // Use theta_lab of negative track if available; otherwise cluster.
  double thetaFlatten = 0;
 
  // One leg, one cluster.
  for (short charge : { -1, 1}) {
    const auto& maximumPtTrack = maximumPtTracks.at(charge);
    if (not maximumPtTrack) {
      continue;
    }
 
    if (maximumPtTrack->clusterEnergySumCMS > 1.5) {
      double invMass = 0.;
      double tempFlatten = 0.;
      for (const SelectedECLCluster& cluster : selectedClusters) {
        double tempInvMass = (maximumPtTrack->p4Lab + cluster.p4Lab).M();
        if (tempInvMass > invMass) {
          invMass = tempInvMass;
          if (charge == 1) {
            tempFlatten = cluster.p4Lab.Theta() * TMath::RadToDeg();
          }
        }
      }
      if (charge == -1) {
        tempFlatten = maximumPtTrack->p4Lab.Theta() * TMath::RadToDeg();
      }
      if (invMass > 5.29) {
        calculationResult["selectee1leg1clst"] = 1;
        thetaFlatten = tempFlatten;
      }
    }
  }
 
  // Two legs
  if (maximumPtTracks.at(-1) and maximumPtTracks.at(1)) {
    const MaximumPtTrack& negativeTrack = *maximumPtTracks.at(-1);
    const MaximumPtTrack& positiveTrack = *maximumPtTracks.at(1);
    const double invMass = (negativeTrack.p4Lab + positiveTrack.p4Lab).M();
    if (invMass > 5.29 and (negativeTrack.clusterEnergySumCMS > 1.5 or positiveTrack.clusterEnergySumCMS > 1.5)) {
      calculationResult["selectee1leg1trk"] = 1;
    }
 
    thetaFlatten = negativeTrack.p4Lab.Theta() * TMath::RadToDeg();
 
    // Two tracks but (exactly) 1 high energy cluster
    const bool missNegClust = positiveTrack.clusterEnergySumCMS > 4.5 and negativeTrack.clusterEnergySumCMS < 1.;
    const bool missPosClust = negativeTrack.clusterEnergySumCMS > 4.5 and positiveTrack.clusterEnergySumCMS < 1.;
    if ((invMass > 9.) and (missNegClust or missPosClust)) {
      calculationResult["eeOneClust"] = 1;
    }
  }
 
  if (calculationResult["selectee1leg1trk"] == 1 or calculationResult["selectee1leg1clst"] == 1) {
    for (int iflat = 0; iflat < 9; iflat++) {
      const std::string& eeFlatName = "eeFlat" + std::to_string(iflat);
      calculationResult[eeFlatName] =
        thetaFlatten >= flatBoundaries[iflat] and thetaFlatten < flatBoundaries[iflat + 1];
      if (calculationResult[eeFlatName]) {
        calculationResult["selectee"] = 1;
      }
    }
  }
 
  // Single tag pi0 / eta dedicated lines
  if (calculationResult["nTrkLoose"] == 1 and calculationResult["maximumPCMS"] > 0.8 and selectedClusters.size() >= 2) {
 
    decltype(selectedClusters) selectedSingleTagClusters(selectedClusters.size());
    auto lastItem = std::copy_if(selectedClusters.begin(), selectedClusters.end(), selectedSingleTagClusters.begin(),
    [](auto & cluster) {
      const bool isNeutralCluster = not cluster.isTrack;
      const bool hasEnoughEnergy = cluster.energyLab > 0.1;
      const double clusterThetaLab = cluster.p4Lab.Theta() * TMath::RadToDeg();
      const bool isInAcceptance = 17 < clusterThetaLab and clusterThetaLab < 150.;
      return isNeutralCluster and hasEnoughEnergy and isInAcceptance;
    });
    selectedSingleTagClusters.resize(std::distance(selectedSingleTagClusters.begin(), lastItem));
 
    if (selectedSingleTagClusters.size() >= 2) {  // One track and at least two clusters are found
 
      const auto& track = maximumPtTracks.at(-1) ? *maximumPtTracks.at(-1) : *maximumPtTracks.at(1);
      // in real signal, the pi0 decay daughters are always the two most-energetic neutral clusters.
      const auto& firstCluster = selectedSingleTagClusters[0];
      const auto& secondCluster = selectedSingleTagClusters[1];
 
      const ROOT::Math::PxPyPzEVector trackP4CMS = track.p4CMS;
      const ROOT::Math::PxPyPzEVector pi0P4CMS = firstCluster.p4CMS + secondCluster.p4CMS;
 
      const bool passPi0ECMS = pi0P4CMS.E() > 1. and pi0P4CMS.E() < 0.525 * p4ofCOM.M();
      const double thetaSumCMS = (pi0P4CMS.Theta() + trackP4CMS.Theta()) * TMath::RadToDeg();
      const bool passThetaSum = thetaSumCMS < 170. or thetaSumCMS > 190.;
 
      double dphiCMS = std::abs(trackP4CMS.Phi() - pi0P4CMS.Phi()) * TMath::RadToDeg();
      if (dphiCMS > 180) {
        dphiCMS = 360 - dphiCMS;
      }
      const bool passdPhi = dphiCMS > 160.;
 
      if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() < 0.7) {
        calculationResult["singleTagLowMass"] = 1;
      } else if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() > 0.7) {
        calculationResult["singleTagHighMass"] = 1;
      }
    }
  }
 
  //..Updated to use new track definitions
  if (calculationResult["nTrkLooseB"] == 1 and calculationResult["maximumPCMSB"] > 0.8 and selectedClusters.size() >= 2) {
 
    decltype(selectedClusters) selectedSingleTagClusters(selectedClusters.size());
    auto lastItem = std::copy_if(selectedClusters.begin(), selectedClusters.end(), selectedSingleTagClusters.begin(),
    [](auto & cluster) {
      const bool isNeutralCluster = not cluster.isTrack;
      const bool hasEnoughEnergy = cluster.energyLab > 0.1;
      const double clusterThetaLab = cluster.p4Lab.Theta() * TMath::RadToDeg();
      const bool isInAcceptance = 17 < clusterThetaLab and clusterThetaLab < 150.;
      return isNeutralCluster and hasEnoughEnergy and isInAcceptance;
    });
    selectedSingleTagClusters.resize(std::distance(selectedSingleTagClusters.begin(), lastItem));
 
    if (selectedSingleTagClusters.size() >= 2) {  // One track and at least two clusters are found
 
      const auto& track = maximumPtTracks.at(-1) ? *maximumPtTracks.at(-1) : *maximumPtTracks.at(1);
      // in real signal, the pi0 decay daughters are always the two most-energetic neutral clusters.
      const auto& firstCluster = selectedSingleTagClusters[0];
      const auto& secondCluster = selectedSingleTagClusters[1];
 
      const ROOT::Math::PxPyPzEVector trackP4CMS = track.p4CMS;
      const ROOT::Math::PxPyPzEVector pi0P4CMS = firstCluster.p4CMS + secondCluster.p4CMS;
 
      const bool passPi0ECMS = pi0P4CMS.E() > 1. and pi0P4CMS.E() < 0.525 * p4ofCOM.M();
      const double thetaSumCMS = (pi0P4CMS.Theta() + trackP4CMS.Theta()) * TMath::RadToDeg();
      const bool passThetaSum = thetaSumCMS < 170. or thetaSumCMS > 190.;
 
      double dphiCMS = std::abs(trackP4CMS.Phi() - pi0P4CMS.Phi()) * TMath::RadToDeg();
      if (dphiCMS > 180) {
        dphiCMS = 360 - dphiCMS;
      }
      const bool passdPhi = dphiCMS > 160.;
 
      if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() < 0.7) {
        calculationResult["singleTagLowMassB"] = 1;
      } else if (passPi0ECMS and passThetaSum and passdPhi and pi0P4CMS.M() > 0.7) {
        calculationResult["singleTagHighMassB"] = 1;
      }
    }
  }
 
  //..Cosmic selection. Need exactly two tracks.
  if (m_tracks.getEntries() == 2) {
 
    const auto negTrack = maximumPtTracksWithoutZCut.at(-1);
    const auto posTrack = maximumPtTracksWithoutZCut.at(1);
 
    if (negTrack and posTrack) {
 
      const double maxNegpT = negTrack->pT;
      const double maxPospT = posTrack->pT;
      const double maxClusterENeg = negTrack->clusterEnergySumLab;
      const double maxClusterEPos = posTrack->clusterEnergySumLab;
 
      const ROOT::Math::PxPyPzEVector& momentumLabNeg(negTrack->p4Lab);
      const ROOT::Math::PxPyPzEVector& momentumLabPos(posTrack->p4Lab);
 
      const double& z0Neg = negTrack->track->getTrackFitResultWithClosestMass(Const::pion)->getZ0();
      const double& d0Neg = negTrack->track->getTrackFitResultWithClosestMass(Const::pion)->getD0();
      const double& z0Pos = posTrack->track->getTrackFitResultWithClosestMass(Const::pion)->getZ0();
      const double& d0Pos = posTrack->track->getTrackFitResultWithClosestMass(Const::pion)->getD0();
 
      // Select cosmic using these tracks
      const bool goodMagneticRegion = (z0Neg<m_goodMagneticRegionZ0 or std::abs(d0Neg)>m_goodMagneticRegionD0)
                                      and (z0Pos<m_goodMagneticRegionZ0 or std::abs(d0Pos)>m_goodMagneticRegionD0);
      if (maxNegpT > m_cosmicMinPt and maxPospT > m_cosmicMinPt and maxClusterENeg < m_cosmicMaxClusterEnergy
          and maxClusterEPos < m_cosmicMaxClusterEnergy and goodMagneticRegion) {
        double dphiLab = std::abs(momentumLabNeg.Phi() - momentumLabPos.Phi()) * TMath::RadToDeg();
        if (dphiLab > 180) {
          dphiLab = 360 - dphiLab;
        }
 
        const double thetaSumLab = (momentumLabNeg.Theta() + momentumLabPos.Theta()) * TMath::RadToDeg();
 
        constexpr double phiBackToBackTolerance = 2.;
        constexpr double thetaBackToBackTolerance = 2.;
        if ((180 - dphiLab) < phiBackToBackTolerance and std::abs(180 - thetaSumLab) < thetaBackToBackTolerance) {
          calculationResult["cosmic"] = 1;
        }
      }
    }
  }
 
  //..Displaced vertex
  // See https://indico.belle2.org/event/8973/ and references therein
  if (maximumPtTracksWithoutZCut.at(-1) and maximumPtTracksWithoutZCut.at(1)) {
 
    //..Make particles of the two highest pt tracks (without IP cut)
    const auto nTrack = maximumPtTracksWithoutZCut.at(-1)->track;
    Particle* nParticle = new Particle(nTrack, Const::pion);
 
    const auto pTrack = maximumPtTracksWithoutZCut.at(1)->track;
    Particle* pParticle = new Particle(pTrack, Const::pion);
 
    //..Make a vertex of these two
    Belle2::analysis::VertexFitKFit* vertexFit = new Belle2::analysis::VertexFitKFit();
    vertexFit->addParticle(nParticle);
    vertexFit->addParticle(pParticle);
    vertexFit->doFit();
 
    //..Vertex properties
    const double chisq = vertexFit->getCHIsq();
    const int ndf = vertexFit->getNDF();
    const double vertexProb = TMath::Prob(chisq, ndf);
    const auto vertexLocation = vertexFit->getVertex();
    const double vertexXY = vertexLocation.perp();
    const double vertexTheta = vertexLocation.theta() * TMath::RadToDeg();
 
    //..Angular differaence of two tracks to reject cosmics
    //  Tolerance could be reduced from 10 deg to 2 deg if needed for physics reasons,
    //  for a 5% increase in the rate of selected displaced vertex triggers.
    //  See https://gitlab.desy.de/belle2/software/basf2/-/merge_requests/1867
    const ROOT::Math::PxPyPzEVector& momentumLabNeg(maximumPtTracksWithoutZCut.at(-1)->p4Lab);
    const ROOT::Math::PxPyPzEVector& momentumLabPos(maximumPtTracksWithoutZCut.at(1)->p4Lab);
    const double thetaSumLab = (momentumLabNeg.Theta() + momentumLabPos.Theta()) * TMath::RadToDeg();
    double dPhiLab = std::abs(momentumLabNeg.Phi() - momentumLabPos.Phi()) * TMath::RadToDeg();
    if (dPhiLab > 180) {
      dPhiLab = 360 - dPhiLab;
    }
    const double backToBackTolerance = 10.; // degrees
    const bool backToBackLab = std::abs(thetaSumLab - 180.) < backToBackTolerance and std::abs(dPhiLab - 180.) < backToBackTolerance;
 
    //..Select a displaced vertex
    const double minProbChiVertex = 0.005; // minimum probability chi square. Many backgrounds have bad chi sq.
    const double minXYVertex = 3.; // minimum xy of vertex (cm). Should pass track filters below this.
    const double maxXYVertex = 60.; // maximum xy. Insufficient CDC for good reconstruction beyond this.
    const double minThetaVertex = 30.; // Large Bhabha background at low angles.
    const double maxThetaVertex = 120.; // Large Bhabha background at low angles.
    if (vertexProb > minProbChiVertex and vertexXY > minXYVertex and vertexXY < maxXYVertex and vertexTheta > minThetaVertex
        and vertexTheta < maxThetaVertex
        and not backToBackLab) {calculationResult["displacedVertex"] = 1;}
 
    //..Clean up
    delete nParticle;
    delete pParticle;
    delete vertexFit;
 
  }
 
}

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.

    {
      const unsigned int sizeBeforeCheck = m_calculationResult.size();
      doCalculation(m_calculationResult);
 
      if (m_calculationResult.size() != sizeBeforeCheck and sizeBeforeCheck > 0) {
        B2WARNING("Calculator added more variables (" << m_calculationResult.size() <<
                  ") than there were before (" << sizeBeforeCheck << "). Probably something strange is going on!");
      }
 
      return m_calculationResult;
    }

requireStoreArrays()

void requireStoreArrays ( )

overridevirtual

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

Implements SoftwareTriggerCalculation.

Definition at line 78 of file FilterCalculator.cc.

{
  m_tracks.isRequired();
  m_eclClusters.isRequired();
  m_l1Trigger.isOptional();
  m_bitsNN.isOptional();
}

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.

    {
      if (not m_debugPrepared) {
        for (auto& identifierWithValue : m_calculationResult) {
          const std::string& identifier = identifierWithValue.first;
          double& value = identifierWithValue.second;
 
          debugOutputTTree->Branch(identifier.c_str(), &value);
        }
        m_debugPrepared = true;
      }
 
      debugOutputTTree->Fill();
    }

Member Data Documentation

m_bitsNN

StoreArray<CDCTriggerUnpacker::NNBitStream> m_bitsNN

private

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:

• hlt/softwaretrigger/calculations/include/FilterCalculator.h
• hlt/softwaretrigger/calculations/src/FilterCalculator.cc