ATLAS Offline Software
Loading...
Searching...
No Matches
DiTauMassTools::MissingMassCalculator Class Reference

#include <MissingMassCalculator.h>

Collaboration diagram for DiTauMassTools::MissingMassCalculator:

Classes

struct  DitauStuff

Public Member Functions

 ~MissingMassCalculator ()
 MissingMassCalculator (MMCCalibrationSet::e aset, std::string paramFilePath)
 MissingMassCalculator (const MissingMassCalculator &)=delete
MissingMassCalculatoroperator= (const MissingMassCalculator &)=delete
int RunMissingMassCalculator (const xAOD::IParticle *part1, const xAOD::IParticle *part2, const xAOD::MissingET *met, const int &njets)
void FinalizeSettings (const xAOD::IParticle *part1, const xAOD::IParticle *part2, const xAOD::MissingET *met, const int &njets)
void SetNiterFit1 (const int val)
void SetNiterFit2 (const int val)
void SetNiterFit3 (const int val)
void SetNiterRandom (const int val)
void SetNsucStop (const int val)
void SetRMSStop (const int val)
void SetMeanbinStop (const double val)
void SetRndmSeedAltering (const int val)
void SetdTheta3d_binMax (const double val)
void SetdTheta3d_binMin (const double val)
void SetEventNumber (const int eventNumber)
void SetMnuScanRange (const double val)
void SetProposalTryMEt (const double val)
void SetProposalTryPhi (const double val)
void SetProposalTryMnu (const double val)
void SetProposalTryEtau (const double val)
void SetUseEfficiencyRecovery (const bool val)
bool GetUseEfficiencyRecovery () const
int GetNiterFit1 () const
int GetNiterFit2 () const
int GetNiterFit3 () const
int GetNiterRandom () const
int GetNsucStop () const
int GetRMSStop () const
double GetMeanbinStop () const
int GetRndmSeedAltering () const
int GetMarkovCountDuplicate () const
int GetMarkovNRejectNoSol () const
int GetMarkovNRejectMetropolis () const
int GetMarkovNAccept () const
int GetMarkovNFullscan () const
double GetProposalTryMEt () const
double GetProposalTryPhi () const
double GetProposalTryMnu () const
double GetProposalTryEtau () const
void SetNsigmaMETscan_ll (const double val)
void SetNsigmaMETscan_lh (const double val)
void SetNsigmaMETscan_hh (const double val)
void SetNsigmaMETscan (const double val)
void SetUseFloatStopping (const bool val)
void SetFloatStoppingMinIter (const int val)
void SetFloatStoppingCheckFreq (const int val)
void SetFloatStoppingComp (const double val)
void SetBeamEnergy (const double val)
void SetLFVLeplepRefit (const bool val)
void SaveLlhHisto (const bool val)
double GetmMaxError () const
double GetmMeanError () const
double GetmInvWidth2Error () const
int GetNNoSol () const
int GetNMetroReject () const
int GetNSol () const
Double_t maxFitting (Double_t *x, Double_t *par)
double maxFromHist (TH1F *theHist, std::vector< double > &histInfo, const MaxHistStrategy::e maxHistStrategy=MaxHistStrategy::FIT, const int winHalfWidth=2, bool debug=false)
double maxFromHist (const std::shared_ptr< TH1F > &theHist, std::vector< double > &histInfo, const MaxHistStrategy::e maxHistStrategy=MaxHistStrategy::FIT, const int winHalfWidth=2, bool debug=false)
double dTheta3DLimit (const int &tau_type, const int &limit_code, const double &P_tau)

Public Attributes

MissingMassInput preparedInput
MissingMassOutput OutputInfo
MissingMassProbProb
XYVector metvec_tmp

Protected Member Functions

int CheckSolutions (PtEtaPhiMVector nu_vec, PtEtaPhiMVector vis_vec, int decayType)
int TailCleanUp (const PtEtaPhiMVector &vis1, const PtEtaPhiMVector &nu1, const PtEtaPhiMVector &vis2, const PtEtaPhiMVector &nu2, const double &mmc_mass, const double &vis_mass, const double &eff_mass, const double &dphiTT)
int refineSolutions (const double &M_nu1, const double &M_nu2, const int nsol1, const int nsol2, const double &Mvis, const double &Meff)
void handleSolutions ()
double MassScale (int method, double mass, const int &tau_type1, const int &tau_type2)
int DitauMassCalculatorV9walk ()
int DitauMassCalculatorV9lfv (bool refit)
int probCalculatorV9fast (const double &phi1, const double &phi2, const double &M_nu1, const double &M_nu2)
void SpaceWalkerInit ()
bool SpaceWalkerWalk ()
bool precomputeCache ()
bool checkMEtInRange ()
bool checkAllParamInRange ()

Private Member Functions

void ClearDitauStuff (DitauStuff &fStuff)
void DoOutputInfo ()
void PrintOtherInput ()
void PrintResults ()
int NuPsolutionV3 (const double &mNu1, const double &mNu2, const double &phi1, const double &phi2, int &nsol1, int &nsol2)
int NuPsolutionLFV (const XYVector &met_vec, const PtEtaPhiMVector &tau, const double &m_nu, std::vector< PtEtaPhiMVector > &nu_vec)

Private Attributes

TRandom2 m_randomGen
MMCCalibrationSet::e m_mmcCalibrationSet {}
bool m_fUseEfficiencyRecovery {}
bool m_fUseFloatStopping {}
int m_fUseFloatStoppingMinIter {}
int m_fUseFloatStoppingCheckFreq {}
double m_fUseFloatStoppingComp {}
int m_nsolmax
int m_nsolfinalmax {}
int m_niterRandomLocal {}
int m_nsucStop {}
int m_rmsStop {}
double m_meanbinStop {}
std::vector< PtEtaPhiMVector > m_nuvecsol1
std::vector< PtEtaPhiMVector > m_nuvecsol2
std::vector< PtEtaPhiMVector > m_tauvecsol1
std::vector< PtEtaPhiMVector > m_tauvecsol2
std::vector< double > m_tauvecprob1
std::vector< double > m_tauvecprob2
std::vector< PtEtaPhiMVector > m_nuvec1_tmp
std::vector< PtEtaPhiMVector > m_nuvec2_tmp
PtEtaPhiMVector m_tautau_tmp
bool m_debugThisIteration
bool m_lfvLeplepRefit
bool m_SaveLlhHisto
int m_nCallprobCalculatorV9fast
double m_nsigma_METscan {}
double m_nsigma_METscan2 {}
double m_nsigma_METscan_ll {}
double m_nsigma_METscan_lh {}
double m_nsigma_METscan_hh {}
double m_nsigma_METscan_lfv_ll {}
double m_nsigma_METscan_lfv_lh {}
double m_beamEnergy {}
int m_iter1 {}
int m_iter2 {}
int m_iter3 {}
int m_iter4 {}
int m_iter5 {}
int m_iang1low {}
int m_iang1high {}
int m_iang2low {}
int m_iang2high {}
int m_iterTheta3d {}
double m_prob_tmp {}
double m_totalProbSum {}
double m_mtautauSum {}
int m_eventNumber {}
int m_seed {}
int m_iter0 {}
int m_iterNuPV3 {}
int m_testptn1 {}
int m_testptn2 {}
int m_testdiscri1 {}
int m_testdiscri2 {}
int m_nosol1 {}
int m_nosol2 {}
int m_iterNsuc {}
bool m_switch1 {}
bool m_switch2 {}
bool m_meanbinToBeEvaluated {}
int m_markovCountDuplicate {}
int m_markovNFullScan {}
int m_markovNRejectNoSol {}
int m_markovNRejectMetropolis {}
int m_markovNAccept {}
double m_PrintmMaxError {}
double m_PrintmMeanError {}
double m_PrintmInvWidth2Error {}
double m_proposalTryMEt {}
double m_ProposalTryPhi {}
double m_ProposalTryMnu {}
double m_ProposalTryEtau {}
double m_mTau {}
double m_mTau2 {}
double m_MEtL {}
double m_MEtP {}
double m_Phi1 {}
double m_Phi2 {}
double m_Mnu1 {}
double m_Mnu2 {}
double m_eTau1 {}
double m_eTau2 {}
double m_eTau10 {}
double m_eTau20 {}
double m_MEtL0 {}
double m_MEtP0 {}
double m_Phi10 {}
double m_Phi20 {}
double m_Mnu10 {}
double m_Mnu20 {}
double m_MEtLMin {}
double m_MEtPMin {}
double m_Phi1Min {}
double m_Phi2Min {}
double m_Mnu1Min {}
double m_Mnu2Min {}
double m_MEtLMax {}
double m_MEtPMax {}
double m_Phi1Max {}
double m_Phi2Max {}
double m_Mnu1Max {}
double m_Mnu2Max {}
double m_MEtLStep {}
double m_MEtPStep {}
double m_Phi1Step {}
double m_Phi2Step {}
double m_Mnu1Step {}
double m_Mnu2Step {}
double m_MEtLRange {}
double m_MEtPRange {}
double m_Phi1Range {}
double m_Phi2Range {}
double m_Mnu1Range {}
double m_Mnu2Range {}
double m_MEtProposal {}
double m_PhiProposal {}
double m_MnuProposal {}
double m_metCovPhiCos {}
double m_metCovPhiSin {}
double m_eTau1Proposal {}
double m_eTau2Proposal {}
double m_eTau1Min {}
double m_eTau1Max {}
double m_eTau1Range {}
double m_eTau2Min {}
double m_eTau2Max {}
double m_eTau2Range {}
bool m_fullParamSpaceScan {}
bool m_Mnu1Exclude {}
int m_nsolOld {}
std::vector< double > m_probFinalSolOldVec
std::vector< double > m_mtautauFinalSolOldVec
std::vector< PtEtaPhiMVector > m_nu1FinalSolOldVec
std::vector< PtEtaPhiMVector > m_nu2FinalSolOldVec
int m_nsol {}
std::vector< double > m_probFinalSolVec
std::vector< double > m_mtautauFinalSolVec
std::vector< PtEtaPhiMVector > m_nu1FinalSolVec
std::vector< PtEtaPhiMVector > m_nu2FinalSolVec
double m_Mnu1ExcludeMin {}
double m_Mnu1ExcludeMax {}
double m_Mnu1ExcludeRange {}
double m_Mnu1XMin {}
double m_Mnu1XMax {}
double m_Mnu1XRange {}
double m_walkWeight {}
double m_cosPhi1 {}
double m_cosPhi2 {}
double m_sinPhi1 {}
double m_sinPhi2 {}
bool m_scanMnu1 {}
bool m_scanMnu2 {}
PtEtaPhiMVector m_tauVec1
PtEtaPhiMVector m_tauVec2
double m_tauVec1Phi {}
double m_tauVec2Phi {}
double m_tauVec1M {}
double m_tauVec2M {}
double m_tauVec1Px {}
double m_tauVec1Py {}
double m_tauVec1Pz {}
double m_tauVec2Px {}
double m_tauVec2Py {}
double m_tauVec2Pz {}
double m_tauVec1P {}
double m_tauVec2P {}
double m_tauVec1E {}
double m_tauVec2E {}
double m_m2Nu1 {}
double m_m2Nu2 {}
double m_ET2v1 {}
double m_ET2v2 {}
double m_E2v1 {}
double m_E2v2 {}
double m_Ev2 {}
double m_Ev1 {}
double m_Mvis {}
double m_Meff {}
bool m_reRunWithBestMET {}
std::shared_ptr< TH1F > m_fMfit_all
std::shared_ptr< TH1F > m_fMEtP_all
std::shared_ptr< TH1F > m_fMEtL_all
std::shared_ptr< TH1F > m_fMnu1_all
std::shared_ptr< TH1F > m_fMnu2_all
std::shared_ptr< TH1F > m_fPhi1_all
std::shared_ptr< TH1F > m_fPhi2_all
std::shared_ptr< TGraph > m_fMfit_allGraph
std::shared_ptr< TH1F > m_fMfit_allNoWeight
std::shared_ptr< TH1F > m_fPXfit1
std::shared_ptr< TH1F > m_fPYfit1
std::shared_ptr< TH1F > m_fPZfit1
std::shared_ptr< TH1F > m_fPXfit2
std::shared_ptr< TH1F > m_fPYfit2
std::shared_ptr< TH1F > m_fPZfit2
std::shared_ptr< TH1F > m_fMmass_split1
std::shared_ptr< TH1F > m_fMEtP_split1
std::shared_ptr< TH1F > m_fMEtL_split1
std::shared_ptr< TH1F > m_fMnu1_split1
std::shared_ptr< TH1F > m_fMnu2_split1
std::shared_ptr< TH1F > m_fPhi1_split1
std::shared_ptr< TH1F > m_fPhi2_split1
std::shared_ptr< TH1F > m_fMmass_split2
std::shared_ptr< TH1F > m_fMEtP_split2
std::shared_ptr< TH1F > m_fMEtL_split2
std::shared_ptr< TH1F > m_fMnu1_split2
std::shared_ptr< TH1F > m_fMnu2_split2
std::shared_ptr< TH1F > m_fPhi1_split2
std::shared_ptr< TH1F > m_fPhi2_split2
TF1 * m_fFitting {}
TH1F * m_fPhi1 {}
TH1F * m_fPhi2 {}
TH1F * m_fMnu1 {}
TH1F * m_fMnu2 {}
TH1F * m_fMetx {}
TH1F * m_fMety {}
TH1F * m_fTheta3D {}
TH1F * m_fTauProb {}
PtEtaPhiMVector m_TLVdummy
DitauStuff m_fDitauStuffFit
DitauStuff m_fDitauStuffHisto
int m_niter_fit1 {}
int m_niter_fit2 {}
int m_niter_fit3 {}
int m_NiterRandom {}
int m_NsucStop {}
int m_RMSStop {}
int m_RndmSeedAltering {}
double m_dTheta3d_binMin {}
double m_dTheta3d_binMax {}
double m_dRmax_tau {}
double m_MnuScanRange {}

Detailed Description

Definition at line 46 of file MissingMassCalculator.h.

Constructor & Destructor Documentation

◆ ~MissingMassCalculator()

MissingMassCalculator::~MissingMassCalculator ( )

Definition at line 193 of file MissingMassCalculator.cxx.

193{ delete Prob; }

◆ MissingMassCalculator() [1/2]

MissingMassCalculator::MissingMassCalculator ( MMCCalibrationSet::e aset,
std::string paramFilePath )

Definition at line 48 of file MissingMassCalculator.cxx.

50 : m_randomGen(), Prob(new MissingMassProb(aset, paramFilePath)) {
52 preparedInput.m_fUseVerbose = 0;
53 preparedInput.m_beamEnergy = 6500.0; // for now LHC default is sqrt(S)=7 TeV
54 m_niter_fit1 = 20;
55 m_niter_fit2 = 30;
56 m_niter_fit3 = 10;
57 m_NsucStop = -1;
58 m_NiterRandom = -1; // if the user does not set it to positive value,will be set
59 // in SpaceWalkerInit
60 m_niterRandomLocal = -1; // niterandom which is really used
61 // to be used with RMSSTOP NiterRandom=10000000; // number of random
62 // iterations for lh. Multiplied by 10 for ll, divided by 10 for hh (to be
63 // optimised)
64 // RMSStop=200;// Stop criteria depending of rms of histogram
65 m_reRunWithBestMET = false;
66 m_RMSStop = -1; // disable
67
68 m_RndmSeedAltering = 0; // can be changed to re-compute with different random seed
69 m_dRmax_tau = 0.4; // changed from 0.2
70 m_nsigma_METscan = -1; // number of sigmas for MET scan
71 m_nsigma_METscan_ll = 3.0; // number of sigmas for MET scan
72 m_nsigma_METscan_lh = 3.0; // number of sigmas for MET scan
73 m_nsigma_METscan_hh = 4.0; // number of sigmas for MET scan (4 for hh 2013)
74 m_nsigma_METscan_lfv_ll = 5.0; // number of sigmas for MET scan (LFV leplep)
75 m_nsigma_METscan_lfv_lh = 5.0; // number of sigmas for MET scan (LFV lephad)
76
77 m_meanbinStop = -1; // meanbin stopping criterion (-1 if not used)
78 m_proposalTryMEt = -1; // loop on METproposal disable // FIXME should be cleaner
79 m_ProposalTryPhi = -1; // loop on Phiproposal disable
80 m_ProposalTryMnu = -1; // loop on MNuProposal disable
81 m_ProposalTryEtau = -1; // loop on ETauProposal disable
82
83 Prob->SetUseTauProbability(true); // TauProbability is ON by default DRMERGE comment out for now
84 Prob->SetUseMnuProbability(false); // MnuProbability is OFF by default
85 Prob->SetUseDphiLL(false); // added by Tomas Davidek for lep-lep
86 m_dTheta3d_binMin = 0.0025;
87 m_dTheta3d_binMax = 0.02;
88 preparedInput.m_METresSyst = 0; // no MET resolution systematics by default (+/-1: up/down 1 sigma)
89 preparedInput.m_dataType = 1; // set to "data" by default
90 preparedInput.m_fUseTailCleanup = 1; // cleanup by default for lep-had Moriond 2012 analysis
91 preparedInput.m_fUseDefaults = 0; // use pre-set defaults for various configurations; if set it to 0
92 // if need to study various options
93 m_fUseEfficiencyRecovery = 0; // no re-fit by default
98
99 preparedInput.m_METScanScheme = 1; // MET-scan scheme: 0- use JER; 1- use simple sumEt & missingHt
100 // for Njet=0 events in (lep-had winter 2012)
101 // MnuScanRange=ParticleConstants::tauMassInMeV / GEV; // range of M(nunu) scan
102 m_MnuScanRange = 1.5; // better value (sacha)
103 preparedInput.m_LFVmode = -1; // by default consider case of H->mu+tau(->ele)
104 preparedInput.ClearInput();
105
107 m_iterTheta3d = 0;
108 m_debugThisIteration = false;
109 m_lfvLeplepRefit = true;
110 m_SaveLlhHisto = false;
111
112 m_nsolmax = 4;
114
115 m_nuvecsol1.resize(m_nsolmax);
116 m_nuvecsol2.resize(m_nsolmax);
117 m_tauvecsol1.resize(m_nsolmax);
118 m_tauvecsol2.resize(m_nsolmax);
119 m_tauvecprob1.resize(m_nsolmax);
120 m_tauvecprob2.resize(m_nsolmax);
121
122 m_nsol = 0;
127
128 m_nsolOld = 0;
133
134 float hEmax = 3000.0; // maximum energy (GeV)
135 // number of bins
136 int hNbins = 1500; // original 2500 for mass, 10000 for P
137 // choice of hNbins also related to size of window for fitting (see
138 // maxFromHist)
139
140 //--- define histograms for histogram method
141 //--- upper limits need to be revisied in the future!!! It may be not enough
142 // for some analyses
143
144 m_fMfit_all = std::make_shared<TH1F>("MMC_h1", "M", hNbins, 0.0,
145 hEmax); // all solutions
146 m_fMfit_all->Sumw2(); // allow proper error bin calculation. Slightly slower but
147 // completely negligible
148
149 // histogram without weight. useful for debugging. negligibly slow until now
151 std::make_shared<TH1F>("MMC_h1NoW", "M no weight", hNbins, 0.0, hEmax); // all solutions
152
153 m_fPXfit1 = std::make_shared<TH1F>("MMC_h2", "Px1", 4 * hNbins, -hEmax,
154 hEmax); // Px for tau1
155 m_fPYfit1 = std::make_shared<TH1F>("MMC_h3", "Py1", 4 * hNbins, -hEmax,
156 hEmax); // Py for tau1
157 m_fPZfit1 = std::make_shared<TH1F>("MMC_h4", "Pz1", 4 * hNbins, -hEmax,
158 hEmax); // Pz for tau1
159 m_fPXfit2 = std::make_shared<TH1F>("MMC_h5", "Px2", 4 * hNbins, -hEmax,
160 hEmax); // Px for tau2
161 m_fPYfit2 = std::make_shared<TH1F>("MMC_h6", "Py2", 4 * hNbins, -hEmax,
162 hEmax); // Py for tau2
163 m_fPZfit2 = std::make_shared<TH1F>("MMC_h7", "Pz2", 4 * hNbins, -hEmax,
164 hEmax); // Pz for tau2
165
166 m_fMfit_all->SetDirectory(0);
167
168 m_fMfit_allNoWeight->SetDirectory(0);
169 m_fPXfit1->SetDirectory(0);
170 m_fPYfit1->SetDirectory(0);
171 m_fPZfit1->SetDirectory(0);
172 m_fPXfit2->SetDirectory(0);
173 m_fPYfit2->SetDirectory(0);
174 m_fPZfit2->SetDirectory(0);
175
176 // max hist fitting function
177 m_fFitting =
178 new TF1("MMC_maxFitting", this, &MissingMassCalculator::maxFitting, 0., hEmax, 3);
179 // Sets initial parameter names
180 m_fFitting->SetParNames("Max", "Mean", "InvWidth2");
181
182 if (preparedInput.m_fUseVerbose == 1) {
183 gDirectory->pwd();
184 gDirectory->ls();
185 }
186
187 if (preparedInput.m_fUseVerbose == 1) {
188 gDirectory->pwd();
189 gDirectory->ls();
190 }
191}
std::vector< PtEtaPhiMVector > m_nu2FinalSolOldVec
std::vector< PtEtaPhiMVector > m_nu1FinalSolOldVec
std::vector< PtEtaPhiMVector > m_tauvecsol1
std::vector< PtEtaPhiMVector > m_nuvecsol1
std::vector< PtEtaPhiMVector > m_nuvecsol2
std::vector< PtEtaPhiMVector > m_nu2FinalSolVec
Double_t maxFitting(Double_t *x, Double_t *par)
std::vector< PtEtaPhiMVector > m_tauvecsol2
std::vector< PtEtaPhiMVector > m_nu1FinalSolVec

◆ MissingMassCalculator() [2/2]

DiTauMassTools::MissingMassCalculator::MissingMassCalculator ( const MissingMassCalculator & )
delete

Member Function Documentation

◆ checkAllParamInRange()

bool MissingMassCalculator::checkAllParamInRange ( )
inlineprotected

Definition at line 2663 of file MissingMassCalculator.cxx.

2663 {
2664
2665 if (m_scanMnu1) {
2666 if (m_Mnu1 < m_Mnu1Min)
2667 return false;
2668 if (m_Mnu1 > m_Mnu1Max)
2669 return false;
2670 if (m_Mnu1 > m_mTau - m_tauVec1M)
2671 return false;
2672 }
2673
2674 if (m_scanMnu2) {
2675 if (m_Mnu2 < m_Mnu2Min)
2676 return false;
2677 if (m_Mnu2 > m_Mnu2Max)
2678 return false;
2679 if (m_Mnu2 > m_mTau - m_tauVec2M)
2680 return false;
2681 }
2682
2683 // FIXME note that since there is a coupling between Met and tau, should
2684 // rigorously test both together however since the 3 sigma range is just a
2685 // hack, it is probably OK
2686
2687 if (m_Phi1 < m_Phi1Min)
2688 return false;
2689 if (m_Phi1 > m_Phi1Max)
2690 return false;
2691
2692 if (m_Phi2 < m_Phi2Min)
2693 return false;
2694 if (m_Phi2 > m_Phi2Max)
2695 return false;
2696
2697 if (!checkMEtInRange())
2698 return false;
2699
2700 return true;
2701}

◆ checkMEtInRange()

bool MissingMassCalculator::checkMEtInRange ( )
inlineprotected

Definition at line 2704 of file MissingMassCalculator.cxx.

2704 {
2705 // check MEt is in allowed range
2706 // range is 3sigma disk ("cutting the corners")
2707 if (std::pow(m_MEtL / preparedInput.m_METsigmaL, 2) +
2708 std::pow(m_MEtP / preparedInput.m_METsigmaP, 2) >
2710 return false;
2711 } else {
2712 return true;
2713 }
2714}

◆ CheckSolutions()

int DiTauMassTools::MissingMassCalculator::CheckSolutions ( PtEtaPhiMVector nu_vec,
PtEtaPhiMVector vis_vec,
int decayType )
inlineprotected

◆ ClearDitauStuff()

void MissingMassCalculator::ClearDitauStuff ( DitauStuff & fStuff)
private

Definition at line 291 of file MissingMassCalculator.cxx.

291 {
292 fStuff.Mditau_best = 0.0;
293 fStuff.Sign_best = 1.0E6;
294 fStuff.nutau1 = PtEtaPhiMVector(0., 0., 0., 0.);
295 fStuff.nutau2 = PtEtaPhiMVector(0., 0., 0., 0.);
296 fStuff.vistau1 = PtEtaPhiMVector(0., 0., 0., 0.);
297 fStuff.vistau2 = PtEtaPhiMVector(0., 0., 0., 0.);
298 fStuff.RMSoverMPV = 0.0;
299
300 return;
301}

◆ DitauMassCalculatorV9lfv()

int MissingMassCalculator::DitauMassCalculatorV9lfv ( bool refit)
inlineprotected

Definition at line 1012 of file MissingMassCalculator.cxx.

1012 {
1013
1014 // debugThisIteration=false;
1015 m_debugThisIteration = true;
1016
1017 int fit_code = 0; // 0==bad, 1==good
1020 OutputInfo.m_NTrials = 0;
1021 OutputInfo.m_NSuccesses = 0;
1022 OutputInfo.m_AveSolRMS = 0.;
1023
1024 //------- Settings -------------------------------
1025 int NiterMET = m_niter_fit2; // number of iterations for each MET scan loop
1026 int NiterMnu = m_niter_fit3; // number of iterations for Mnu loop
1027 const double Mtau = ParticleConstants::tauMassInMeV / GEV;
1028 double Mnu_binSize = m_MnuScanRange / NiterMnu;
1029
1030 double METresX = preparedInput.m_METsigmaL; // MET resolution in direction parallel to
1031 // leading jet, for MET scan
1032 double METresY = preparedInput.m_METsigmaP; // MET resolution in direction perpendicular to
1033 // leading jet, for MET scan
1034
1035 //-------- end of Settings
1036
1037 // if m_nsigma_METscan was not set by user, set to default values
1038 if(m_nsigma_METscan == -1){
1039 if (preparedInput.m_tauTypes == TauTypes::ll) { // both tau's are leptonic
1041 } else if (preparedInput.m_tauTypes == TauTypes::lh) { // lep had
1043 }
1044 }
1045
1046 double N_METsigma = m_nsigma_METscan; // number of sigmas for MET scan
1047 double METresX_binSize = 2 * N_METsigma * METresX / NiterMET;
1048 double METresY_binSize = 2 * N_METsigma * METresY / NiterMET;
1049
1050 int solution = 0;
1051
1052 std::vector<PtEtaPhiMVector> nu_vec;
1053
1054 m_totalProbSum = 0;
1055 m_mtautauSum = 0;
1056
1057 double metprob = 1.0;
1058 double sign_tmp = 0.0;
1059 double tauprob = 1.0;
1060 double totalProb = 0.0;
1061
1062 m_prob_tmp = 0.0;
1063
1064 double met_smear_x = 0.0;
1065 double met_smear_y = 0.0;
1066 double met_smearL = 0.0;
1067 double met_smearP = 0.0;
1068
1069 double angle1 = 0.0;
1070
1071 if (m_fMfit_all) {
1072 m_fMfit_all->Reset();
1073 }
1074 if (m_fMfit_allNoWeight) {
1075 m_fMfit_allNoWeight->Reset();
1076 }
1077 if (m_fPXfit1) {
1078 m_fPXfit1->Reset();
1079 }
1080 if (m_fPYfit1) {
1081 m_fPYfit1->Reset();
1082 }
1083 if (m_fPZfit1) {
1084 m_fPZfit1->Reset();
1085 }
1086
1087 int iter0 = 0;
1088 m_iter1 = 0;
1089 m_iter2 = 0;
1090 m_iter3 = 0;
1091 m_iter4 = 0;
1092
1093 const double met_coscovphi = cos(preparedInput.m_METcovphi);
1094 const double met_sincovphi = sin(preparedInput.m_METcovphi);
1095
1096 m_iang1low = 0;
1097 m_iang1high = 0;
1098
1099 // double Mvis=(tau_vec1+tau_vec2).M();
1100 // PtEtaPhiMVector met4vec(0.0,0.0,0.0,0.0);
1101 // met4vec.SetPxPyPzE(met_vec.X(),met_vec.Y(),0.0,met_vec.R());
1102 // double Meff=(tau_vec1+tau_vec2+met4vec).M();
1103 // double met_det=met_vec.R();
1104
1105 //---------------------------------------------
1106 if (preparedInput.m_tauTypes == TauTypes::ll) // dilepton case
1107 {
1108 if (preparedInput.m_fUseVerbose == 1) {
1109 Info("DiTauMassTools", "Running in dilepton mode");
1110 }
1111 double input_metX = preparedInput.m_MetVec.X();
1112 double input_metY = preparedInput.m_MetVec.Y();
1113
1114 PtEtaPhiMVector tau_tmp(0.0, 0.0, 0.0, 0.0);
1115 PtEtaPhiMVector lep_tmp(0.0, 0.0, 0.0, 0.0);
1116 int tau_type_tmp;
1117 int tau_ind = 0;
1118
1119 if (preparedInput.m_LFVmode == 1) // muon case: H->mu+tau(->ele) decays
1120 {
1121 if ((preparedInput.m_vistau1.M() > 0.05 &&
1122 preparedInput.m_vistau2.M() < 0.05) != refit) // choosing lepton from Higgs decay
1123 //When the mass calculator is rerun with refit==true the alternative lepton ordering is used
1124 {
1125 tau_tmp = preparedInput.m_vistau2;
1126 lep_tmp = preparedInput.m_vistau1;
1127 tau_type_tmp = preparedInput.m_type_visTau2;
1128 tau_ind = 2;
1129 } else {
1130 tau_tmp = preparedInput.m_vistau1;
1131 lep_tmp = preparedInput.m_vistau2;
1132 tau_type_tmp = preparedInput.m_type_visTau1;
1133 tau_ind = 1;
1134 }
1135 }
1136 if (preparedInput.m_LFVmode == 0) // electron case: H->ele+tau(->mu) decays
1137 {
1138 if ((preparedInput.m_vistau1.M() < 0.05 &&
1139 preparedInput.m_vistau2.M() > 0.05) != refit) // choosing lepton from Higgs decay
1140 //When the mass calculator is rerun with refit=true the alternative lepton ordering is used
1141 {
1142 tau_tmp = preparedInput.m_vistau2;
1143 lep_tmp = preparedInput.m_vistau1;
1144 tau_type_tmp = preparedInput.m_type_visTau2;
1145 tau_ind = 2;
1146 } else {
1147 tau_tmp = preparedInput.m_vistau1;
1148 lep_tmp = preparedInput.m_vistau2;
1149 tau_type_tmp = preparedInput.m_type_visTau1;
1150 tau_ind = 1;
1151 }
1152 }
1153
1154 //------- Settings -------------------------------
1155 double Mlep = tau_tmp.M();
1156 // double dMnu_max=m_MnuScanRange-Mlep;
1157 // double Mnu_binSize=dMnu_max/NiterMnu;
1158 //-------- end of Settings
1159
1160 // double M=Mtau;
1161 double M_nu = 0.0;
1162 double MnuProb = 1.0;
1163 //---------------------------------------------
1164 for (int i3 = 0; i3 < NiterMnu; i3++) //---- loop-3: virtual neutrino mass
1165 {
1166 M_nu = Mnu_binSize * i3;
1167 if (M_nu >= (Mtau - Mlep))
1168 continue;
1169 // M=sqrt(Mtau*Mtau-M_nu*M_nu);
1170 MnuProb = Prob->MnuProbability(preparedInput, M_nu,
1171 Mnu_binSize); // Mnu probability
1172 //---------------------------------------------
1173 for (int i4 = 0; i4 < NiterMET + 1; i4++) // MET_X scan
1174 {
1175 met_smearL = METresX_binSize * i4 - N_METsigma * METresX;
1176 for (int i5 = 0; i5 < NiterMET + 1; i5++) // MET_Y scan
1177 {
1178 met_smearP = METresY_binSize * i5 - N_METsigma * METresY;
1179 if (pow(met_smearL / METresX, 2) + pow(met_smearP / METresY, 2) > pow(N_METsigma, 2))
1180 continue; // use ellipse instead of square
1181 met_smear_x = met_smearL * met_coscovphi - met_smearP * met_sincovphi;
1182 met_smear_y = met_smearL * met_sincovphi + met_smearP * met_coscovphi;
1183 metvec_tmp.SetXY(input_metX + met_smear_x, input_metY + met_smear_y);
1184
1185 solution = NuPsolutionLFV(metvec_tmp, tau_tmp, M_nu, nu_vec);
1186
1187 ++iter0;
1188
1189 if (solution < 1)
1190 continue;
1191 ++m_iter1;
1192
1193 // if fast sin cos, result to not match exactly nupsolutionv2, so skip
1194 // test
1195 // SpeedUp no nested loop to compute individual probability
1196 int ngoodsol1 = 0;
1197
1198 metprob = Prob->MetProbability(preparedInput, met_smearL, met_smearP, METresX, METresY);
1199 if (metprob <= 0)
1200 continue;
1201 for (unsigned int j1 = 0; j1 < nu_vec.size(); j1++) {
1202 if (tau_tmp.E() + nu_vec[j1].E() >= preparedInput.m_beamEnergy)
1203 continue;
1204 const double tau1_tmpp = (tau_tmp + nu_vec[j1]).P();
1205 angle1 = Angle(nu_vec[j1], tau_tmp);
1206
1207 if (angle1 < dTheta3DLimit(tau_type_tmp, 0, tau1_tmpp)) {
1208 ++m_iang1low;
1209 continue;
1210 } // lower 99% bound
1211 if (angle1 > dTheta3DLimit(tau_type_tmp, 1, tau1_tmpp)) {
1212 ++m_iang1high;
1213 continue;
1214 } // upper 99% bound
1215 double tauvecprob1j =
1216 Prob->dTheta3d_probabilityFast(preparedInput, tau_type_tmp, angle1, tau1_tmpp);
1217 if (tauvecprob1j == 0.)
1218 continue;
1219 tauprob = Prob->TauProbabilityLFV(preparedInput, tau_type_tmp, tau_tmp, nu_vec[j1]);
1220 totalProb = tauvecprob1j * metprob * MnuProb * tauprob;
1221
1222 m_tautau_tmp.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);
1223 m_tautau_tmp += tau_tmp;
1224 m_tautau_tmp += lep_tmp;
1225 m_tautau_tmp += nu_vec[j1];
1226
1227 const double mtautau = m_tautau_tmp.M();
1228
1229 m_totalProbSum += totalProb;
1230 m_mtautauSum += mtautau;
1231
1232 fit_code = 1; // at least one solution is found
1233
1234 m_fMfit_all->Fill(mtautau, totalProb);
1235 m_fMfit_allNoWeight->Fill(mtautau, 1.);
1236 //----------------- using P*fit to fill Px,y,z_tau
1237 m_fPXfit1->Fill((tau_tmp + nu_vec[j1]).Px(), totalProb);
1238 m_fPYfit1->Fill((tau_tmp + nu_vec[j1]).Py(), totalProb);
1239 m_fPZfit1->Fill((tau_tmp + nu_vec[j1]).Pz(), totalProb);
1240
1241 if (totalProb > m_prob_tmp) // fill solution with highest probability
1242 {
1243 sign_tmp = -log10(totalProb);
1244 m_prob_tmp = totalProb;
1245 m_fDitauStuffFit.Mditau_best = mtautau;
1246 m_fDitauStuffFit.Sign_best = sign_tmp;
1247 if (tau_ind == 1)
1248 m_fDitauStuffFit.nutau1 = nu_vec[j1];
1249 if (tau_ind == 2)
1250 m_fDitauStuffFit.nutau2 = nu_vec[j1];
1251 }
1252
1253 ++ngoodsol1;
1254 }
1255
1256 if (ngoodsol1 == 0)
1257 continue;
1258 m_iter2 += 1;
1259
1260 m_iter3 += 1;
1261 }
1262 }
1263 }
1264 } else if (preparedInput.m_tauTypes == TauTypes::lh) // lepton+tau case
1265 {
1266 if (preparedInput.m_fUseVerbose == 1) {
1267 Info("DiTauMassTools", "Running in lepton+tau mode");
1268 }
1269 //------- Settings -------------------------------
1270
1271 //----- Stuff below are for Winter 2012 lep-had analysis only; it has to be
1272 // replaced by a more common scheme once other channels are optimized
1273 // XYVector
1274 // mht_vec((tau_vec1+tau_vec2).Px(),(tau_vec1+tau_vec2).Py()); //
1275 // missing Ht vector for Njet25=0 events const double
1276 // mht=mht_vec.R();
1277 double input_metX = preparedInput.m_MetVec.X();
1278 double input_metY = preparedInput.m_MetVec.Y();
1279
1280 // double mht_offset=0.0;
1281 // if(InputInfo.UseHT) // use missing Ht (for 0-jet events only for
1282 // now)
1283 // {
1284 // input_metX=-mht_vec.X();
1285 // input_metY=-mht_vec.Y();
1286 // }
1287 // else // use MET (for 0-jet and 1-jet events)
1288 // {
1289 // input_metX=met_vec.X();
1290 // input_metY=met_vec.Y();
1291 // }
1292
1293 PtEtaPhiMVector tau_tmp(0.0, 0.0, 0.0, 0.0);
1294 PtEtaPhiMVector lep_tmp(0.0, 0.0, 0.0, 0.0);
1295 int tau_type_tmp;
1296 if (preparedInput.m_type_visTau1 == 8) {
1297 tau_tmp = preparedInput.m_vistau2;
1298 lep_tmp = preparedInput.m_vistau1;
1299 tau_type_tmp = preparedInput.m_type_visTau2;
1300 }
1301 if (preparedInput.m_type_visTau2 == 8) {
1302 tau_tmp = preparedInput.m_vistau1;
1303 lep_tmp = preparedInput.m_vistau2;
1304 tau_type_tmp = preparedInput.m_type_visTau1;
1305 }
1306
1307 //---------------------------------------------
1308 for (int i4 = 0; i4 < NiterMET + 1; i4++) // MET_X scan
1309 {
1310 met_smearL = METresX_binSize * i4 - N_METsigma * METresX;
1311 for (int i5 = 0; i5 < NiterMET + 1; i5++) // MET_Y scan
1312 {
1313 met_smearP = METresY_binSize * i5 - N_METsigma * METresY;
1314 if (pow(met_smearL / METresX, 2) + pow(met_smearP / METresY, 2) > pow(N_METsigma, 2))
1315 continue; // use ellipse instead of square
1316 met_smear_x = met_smearL * m_metCovPhiCos - met_smearP * m_metCovPhiSin;
1317 met_smear_y = met_smearL * m_metCovPhiSin + met_smearP * m_metCovPhiCos;
1318 metvec_tmp.SetXY(input_metX + met_smear_x, input_metY + met_smear_y);
1319
1320 solution = NuPsolutionLFV(metvec_tmp, tau_tmp, 0.0, nu_vec);
1321
1322 ++iter0;
1323
1324 if (solution < 1)
1325 continue;
1326 ++m_iter1;
1327
1328 // if fast sin cos, result to not match exactly nupsolutionv2, so skip
1329 // test
1330 // SpeedUp no nested loop to compute individual probability
1331 int ngoodsol1 = 0;
1332
1333 metprob = Prob->MetProbability(preparedInput, met_smearL, met_smearP, METresX, METresY);
1334 if (metprob <= 0)
1335 continue;
1336 for (unsigned int j1 = 0; j1 < nu_vec.size(); j1++) {
1337 if (tau_tmp.E() + nu_vec[j1].E() >= preparedInput.m_beamEnergy)
1338 continue;
1339 const double tau1_tmpp = (tau_tmp + nu_vec[j1]).P();
1340 angle1 = Angle(nu_vec[j1], tau_tmp);
1341
1342 if (angle1 < dTheta3DLimit(tau_type_tmp, 0, tau1_tmpp)) {
1343 ++m_iang1low;
1344 continue;
1345 } // lower 99% bound
1346 if (angle1 > dTheta3DLimit(tau_type_tmp, 1, tau1_tmpp)) {
1347 ++m_iang1high;
1348 continue;
1349 } // upper 99% bound
1350 double tauvecprob1j =
1351 Prob->dTheta3d_probabilityFast(preparedInput, tau_type_tmp, angle1, tau1_tmpp);
1352 if (tauvecprob1j == 0.)
1353 continue;
1354 tauprob = Prob->TauProbabilityLFV(preparedInput, tau_type_tmp, tau_tmp, nu_vec[j1]);
1355 totalProb = tauvecprob1j * metprob * tauprob;
1356
1357 m_tautau_tmp.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);
1358 m_tautau_tmp += tau_tmp;
1359 m_tautau_tmp += lep_tmp;
1360 m_tautau_tmp += nu_vec[j1];
1361
1362 const double mtautau = m_tautau_tmp.M();
1363
1364 m_totalProbSum += totalProb;
1365 m_mtautauSum += mtautau;
1366
1367 fit_code = 1; // at least one solution is found
1368
1369 m_fMfit_all->Fill(mtautau, totalProb);
1370 m_fMfit_allNoWeight->Fill(mtautau, 1.);
1372 // m_fPXfit1->Fill((tau_tmp+nu_vec[j1]).Px(),totalProb);
1373 // m_fPYfit1->Fill((tau_tmp+nu_vec[j1]).Py(),totalProb);
1374 // m_fPZfit1->Fill((tau_tmp+nu_vec[j1]).Pz(),totalProb);
1375
1376 if (totalProb > m_prob_tmp) // fill solution with highest probability
1377 {
1378 sign_tmp = -log10(totalProb);
1379 m_prob_tmp = totalProb;
1380 m_fDitauStuffFit.Mditau_best = mtautau;
1381 m_fDitauStuffFit.Sign_best = sign_tmp;
1382 if (preparedInput.m_type_visTau1 == 8) {
1383 m_fDitauStuffFit.vistau1 = lep_tmp;
1384 m_fDitauStuffFit.vistau2 = tau_tmp;
1385 m_fDitauStuffFit.nutau2 = nu_vec[j1];
1386 } else if (preparedInput.m_type_visTau2 == 8) {
1387 m_fDitauStuffFit.vistau2 = lep_tmp;
1388 m_fDitauStuffFit.vistau1 = tau_tmp;
1389 m_fDitauStuffFit.nutau1 = nu_vec[j1];
1390 }
1391 }
1392
1393 ++ngoodsol1;
1394 }
1395
1396 if (ngoodsol1 == 0)
1397 continue;
1398 m_iter2 += 1;
1399
1400 m_iter3 += 1;
1401 }
1402 }
1403 } else {
1404 Info("DiTauMassTools", "Running in an unknown mode?!?!");
1405 }
1406
1407 OutputInfo.m_NTrials = iter0;
1408 OutputInfo.m_NSuccesses = m_iter3;
1409
1410 if (preparedInput.m_fUseVerbose == 1) {
1411 Info("DiTauMassTools", "%s",
1412 ("SpeedUp niters=" + std::to_string(iter0) + " " + std::to_string(m_iter1) + " " +
1413 std::to_string(m_iter2) + " " + std::to_string(m_iter3) + "skip:" + std::to_string(m_iang1low) +
1414 " " + std::to_string(m_iang1high))
1415 .c_str());
1416 }
1417
1418 if (m_fMfit_all->GetEntries() > 0 && m_iter3 > 0) {
1419#ifdef SMOOTH
1420 m_fMfit_all->Smooth();
1421 m_fMfit_allNoWeight->Smooth();
1422 m_fPXfit1->Smooth();
1423 m_fPYfit1->Smooth();
1424 m_fPZfit1->Smooth();
1425#endif
1426
1427 // default max finding method defined in MissingMassCalculator.h
1428 // note that window defined in terms of number of bin, so depend on binning
1429 std::vector<double> histInfo(HistInfo::MAXHISTINFO);
1430 m_fDitauStuffHisto.Mditau_best = maxFromHist(m_fMfit_all, histInfo);
1431 double prob_hist = histInfo.at(HistInfo::PROB);
1432
1433 if (prob_hist != 0.0)
1434 m_fDitauStuffHisto.Sign_best = -log10(std::abs(prob_hist));
1435 else {
1436 // this mean the histogram is empty.
1437 // possible but very rare if all entries outside histogram range
1438 // fall back to maximum
1439 m_fDitauStuffHisto.Sign_best = -999.;
1440 m_fDitauStuffHisto.Mditau_best = m_fDitauStuffFit.Mditau_best;
1441 }
1442
1443 if (m_fDitauStuffHisto.Mditau_best > 0.0)
1444 m_fDitauStuffHisto.RMSoverMPV = m_fMfit_all->GetRMS() / m_fDitauStuffHisto.Mditau_best;
1445 std::vector<double> histInfoOther(HistInfo::MAXHISTINFO);
1446 //---- getting Nu1
1447 double Px1 = maxFromHist(m_fPXfit1, histInfoOther);
1448 double Py1 = maxFromHist(m_fPYfit1, histInfoOther);
1449 double Pz1 = maxFromHist(m_fPZfit1, histInfoOther);
1450 //---- setting 4-vecs
1451 PxPyPzMVector nu1_tmp(0.0, 0.0, 0.0, 0.0);
1452 PxPyPzMVector nu2_tmp(0.0, 0.0, 0.0, 0.0);
1453 if (preparedInput.m_type_visTau1 == 8) {
1454 nu1_tmp = preparedInput.m_vistau1;
1455 nu2_tmp.SetCoordinates(Px1, Py1, Pz1, ParticleConstants::tauMassInMeV / GEV);
1456 }
1457 if (preparedInput.m_type_visTau2 == 8) {
1458 nu2_tmp = preparedInput.m_vistau2;
1459 nu1_tmp.SetCoordinates(Px1, Py1, Pz1, ParticleConstants::tauMassInMeV / GEV);
1460 }
1461 m_fDitauStuffHisto.nutau1 = nu1_tmp - preparedInput.m_vistau1;
1462 m_fDitauStuffHisto.nutau2 = nu2_tmp - preparedInput.m_vistau2;
1463 }
1464 if (m_lfvLeplepRefit && fit_code==0 && !refit) {
1465 fit_code = DitauMassCalculatorV9lfv(true);
1466 return fit_code;
1467 }
1468
1469
1470
1471 if (preparedInput.m_fUseVerbose == 1) {
1472 if (fit_code == 0) {
1473 Info(
1474 "DiTauMassTools", "%s",
1475 ("!!!----> Warning-3 in MissingMassCalculator::DitauMassCalculatorV9lfv() : fit status=" +
1476 std::to_string(fit_code))
1477 .c_str());
1478 Info("DiTauMassTools", "....... No solution is found. Printing input info .......");
1479
1480 Info("DiTauMassTools", "%s", (" vis Tau-1: Pt="+std::to_string(preparedInput.m_vistau1.Pt())
1481 +" M="+std::to_string(preparedInput.m_vistau1.M())+" eta="+std::to_string(preparedInput.m_vistau1.Eta())
1482 +" phi="+std::to_string(preparedInput.m_vistau1.Phi())
1483 +" type="+std::to_string(preparedInput.m_type_visTau1)).c_str());
1484 Info("DiTauMassTools", "%s", (" vis Tau-2: Pt="+std::to_string(preparedInput.m_vistau2.Pt())
1485 +" M="+std::to_string(preparedInput.m_vistau2.M())+" eta="+std::to_string(preparedInput.m_vistau2.Eta())
1486 +" phi="+std::to_string(preparedInput.m_vistau2.Phi())
1487 +" type="+std::to_string(preparedInput.m_type_visTau2)).c_str());
1488 Info("DiTauMassTools", "%s", (" MET="+std::to_string(preparedInput.m_MetVec.R())+" Met_X="+std::to_string(preparedInput.m_MetVec.X())
1489 +" Met_Y="+std::to_string(preparedInput.m_MetVec.Y())).c_str());
1490 Info("DiTauMassTools", " ---------------------------------------------------------- ");
1491 }
1492 }
1493 return fit_code;
1494}
static Double_t P(Double_t *tt, Double_t *par)
#define GEV
constexpr int pow(int base, int exp) noexcept
double maxFromHist(TH1F *theHist, std::vector< double > &histInfo, const MaxHistStrategy::e maxHistStrategy=MaxHistStrategy::FIT, const int winHalfWidth=2, bool debug=false)
double dTheta3DLimit(const int &tau_type, const int &limit_code, const double &P_tau)
int NuPsolutionLFV(const XYVector &met_vec, const PtEtaPhiMVector &tau, const double &m_nu, std::vector< PtEtaPhiMVector > &nu_vec)
double Angle(const VectorType1 &vec1, const VectorType2 &vec2)
constexpr double tauMassInMeV
the mass of the tau (in MeV)
@ Info
Definition ZDCMsg.h:20

◆ DitauMassCalculatorV9walk()

int MissingMassCalculator::DitauMassCalculatorV9walk ( )
inlineprotected

Definition at line 792 of file MissingMassCalculator.cxx.

792 {
793
794 int nsuccesses = 0;
795
796 int fit_code = 0; // 0==bad, 1==good
799 OutputInfo.m_AveSolRMS = 0.;
800
801 m_fMfit_all->Reset();
802
803 if(m_SaveLlhHisto){
804 m_fMEtP_all->Reset();
805 m_fMEtL_all->Reset();
806 m_fMnu1_all->Reset();
807 m_fMnu2_all->Reset();
808 m_fPhi1_all->Reset();
809 m_fPhi2_all->Reset();
810 }
811
812 m_fMfit_allNoWeight->Reset();
813 m_fPXfit1->Reset();
814 m_fPYfit1->Reset();
815 m_fPZfit1->Reset();
816 m_fPXfit2->Reset();
817 m_fPYfit2->Reset();
818 m_fPZfit2->Reset();
819
820 // these histograms are used for the floating stopping criterion
822 m_fMmass_split1->Reset();
823 m_fMEtP_split1->Reset();
824 m_fMEtL_split1->Reset();
825 m_fMnu1_split1->Reset();
826 m_fMnu2_split1->Reset();
827 m_fPhi1_split1->Reset();
828 m_fPhi2_split1->Reset();
829 m_fMmass_split2->Reset();
830 m_fMEtP_split2->Reset();
831 m_fMEtL_split2->Reset();
832 m_fMnu1_split2->Reset();
833 m_fMnu2_split2->Reset();
834 m_fPhi1_split2->Reset();
835 m_fPhi2_split2->Reset();
836 }
837
838 m_prob_tmp = 0.0;
839
840 m_iter1 = 0;
841
842 m_totalProbSum = 0;
843 m_mtautauSum = 0;
844
845 // initialize a spacewalker, which walks the parameter space according to some
846 // algorithm
848
849 while (SpaceWalkerWalk()) {
850 bool paramInsideRange = false;
851 m_nsol = 0;
852
853 paramInsideRange = checkAllParamInRange();
854
855 // FIXME if no tau scanning, or symmetric matrices, rotatin is made twice
856 // which is inefficient
857 const double deltaMetx = m_MEtL * m_metCovPhiCos - m_MEtP * m_metCovPhiSin;
858 const double deltaMety = m_MEtL * m_metCovPhiSin + m_MEtP * m_metCovPhiCos;
859
860 // deltaMetVec.Set(met_smear_x,met_smear_y);
861 preparedInput.m_metVec.SetXY(preparedInput.m_inputMEtX + deltaMetx,
862 preparedInput.m_inputMEtY + deltaMety);
863
864 // save in global variable for speed sake
865 preparedInput.m_MEtX = preparedInput.m_metVec.X();
866 preparedInput.m_MEtY = preparedInput.m_metVec.Y();
867 preparedInput.m_MEtT = preparedInput.m_metVec.R();
868
869 if (paramInsideRange)
871
872 // DR for markov chain need to enter handleSolution also when zero solutions
874 // be careful that with markov, current solution is from now on stored in
875 // XYZOldSolVec
876
877 if (m_nsol <= 0)
878 continue;
879
880 // for markov, nsuccess more difficult to define. Decide this is the number
881 // of independent point accepted (hence without weight)
882 nsuccesses = m_markovNAccept;
884
885 m_iter1 += m_nsol;
886 fit_code = 1;
887
888 } // while loop
889
890 OutputInfo.m_NTrials = m_iter0;
891 OutputInfo.m_NSuccesses = nsuccesses;
892
893 if (nsuccesses > 0) {
894 OutputInfo.m_AveSolRMS /= nsuccesses;
895 } else {
896 OutputInfo.m_AveSolRMS = -1.;
897 }
898
899 double Px1, Py1, Pz1;
900 double Px2, Py2, Pz2;
901 if (nsuccesses > 0) {
902
903 // note that smoothing can slightly change the integral of the histogram
904
905#ifdef SMOOTH
906 m_fMfit_all->Smooth();
907 m_fMfit_allNoWeight->Smooth();
908 m_fPXfit1->Smooth();
909 m_fPYfit1->Smooth();
910 m_fPZfit1->Smooth();
911 m_fPXfit2->Smooth();
912 m_fPYfit2->Smooth();
913 m_fPZfit2->Smooth();
914#endif
915
916 // default max finding method defined in MissingMassCalculator.h
917 // note that window defined in terms of number of bin, so depend on binning
918 std::vector<double> histInfo(HistInfo::MAXHISTINFO);
919 m_fDitauStuffHisto.Mditau_best = maxFromHist(m_fMfit_all, histInfo);
920 double prob_hist = histInfo.at(HistInfo::PROB);
921
922 if (prob_hist != 0.0)
923 m_fDitauStuffHisto.Sign_best = -log10(std::abs(prob_hist));
924 else {
925 // this mean the histogram is empty.
926 // possible but very rare if all entries outside histogram range
927 // fall back to maximum
928 m_fDitauStuffHisto.Sign_best = -999.;
929 m_fDitauStuffHisto.Mditau_best = m_fDitauStuffFit.Mditau_best;
930 }
931
932 if (m_fDitauStuffHisto.Mditau_best > 0.0)
933 m_fDitauStuffHisto.RMSoverMPV = m_fMfit_all->GetRMS() / m_fDitauStuffHisto.Mditau_best;
934 std::vector<double> histInfoOther(HistInfo::MAXHISTINFO);
935 //---- getting full tau1 momentum
936 Px1 = maxFromHist(m_fPXfit1, histInfoOther);
937 Py1 = maxFromHist(m_fPYfit1, histInfoOther);
938 Pz1 = maxFromHist(m_fPZfit1, histInfoOther);
939
940 //---- getting full tau2 momentum
941 Px2 = maxFromHist(m_fPXfit2, histInfoOther);
942 Py2 = maxFromHist(m_fPYfit2, histInfoOther);
943 Pz2 = maxFromHist(m_fPZfit2, histInfoOther);
944
945 //---- setting 4-vecs
946 PxPyPzMVector fulltau1, fulltau2;
947 fulltau1.SetCoordinates(Px1, Py1, Pz1, ParticleConstants::tauMassInMeV / GEV);
948 fulltau2.SetCoordinates(Px2, Py2, Pz2, ParticleConstants::tauMassInMeV / GEV);
949 // PtEtaPhiMVector fulltau1(_fulltau1.Pt(), _fulltau1.Eta(), _fulltau1.Phi(), _fulltau1.M());
950 //PtEtaPhiMVector fulltau2(_fulltau2.Pt(), _fulltau2.Eta(), _fulltau2.Phi(), _fulltau2.M());
951
952 if (fulltau1.P() < preparedInput.m_vistau1.P())
953 fulltau1 = 1.01 * preparedInput.m_vistau1; // protection against cases when fitted tau
954 // momentum is smaller than visible tau momentum
955 if (fulltau2.P() < preparedInput.m_vistau2.P())
956 fulltau2 = 1.01 * preparedInput.m_vistau2; // protection against cases when fitted tau
957 // momentum is smaller than visible tau momentum
958 m_fDitauStuffHisto.vistau1 = preparedInput.m_vistau1; // FIXME should also be fitted if tau scan
959 m_fDitauStuffHisto.vistau2 = preparedInput.m_vistau2;
960 m_fDitauStuffHisto.nutau1 = fulltau1 - preparedInput.m_vistau1; // these are the original tau vis
961 m_fDitauStuffHisto.nutau2 =
962 fulltau2 - preparedInput.m_vistau2; // FIXME neutrino mass not necessarily zero
963 }
964
965 // Note that for v9walk, points outside the METx MEty disk are counted, while
966 // this was not the case for v9
967 if (preparedInput.m_fUseVerbose == 1) {
968 Info("DiTauMassTools", "Scanning ");
969 Info("DiTauMassTools", " Markov ");
970 Info("DiTauMassTools", "%s",
971 (" V9W niters=" + std::to_string(m_iter0) + " " + std::to_string(m_iter1)).c_str());
972 Info("DiTauMassTools", "%s", (" nFullScan " + std::to_string(m_markovNFullScan)).c_str());
973 Info("DiTauMassTools", "%s", (" nRejectNoSol " + std::to_string(m_markovNRejectNoSol)).c_str());
974 Info("DiTauMassTools", "%s", (" nRejectMetro " + std::to_string(m_markovNRejectMetropolis)).c_str());
975 Info("DiTauMassTools", "%s", (" nAccept " + std::to_string(m_markovNAccept)).c_str());
976 Info("DiTauMassTools", "%s",
977 (" probsum " + std::to_string(m_totalProbSum) + " msum " + std::to_string(m_mtautauSum))
978 .c_str());
979 }
980
981 if (preparedInput.m_fUseVerbose == 1) {
982 if (fit_code == 0) {
983 Info("DiTauMassTools", "%s", ("!!!----> Warning-3 in "
984 "MissingMassCalculator::DitauMassCalculatorV9Walk() : fit status=" +
985 std::to_string(fit_code))
986 .c_str());
987 Info("DiTauMassTools", "%s", "....... No solution is found. Printing input info .......");
988
989 Info("DiTauMassTools", "%s", (" vis Tau-1: Pt=" + std::to_string(preparedInput.m_vistau1.Pt()) +
990 " M=" + std::to_string(preparedInput.m_vistau1.M()) +
991 " eta=" + std::to_string(preparedInput.m_vistau1.Eta()) +
992 " phi=" + std::to_string(preparedInput.m_vistau1.Phi()) +
993 " type=" + std::to_string(preparedInput.m_type_visTau1))
994 .c_str());
995 Info("DiTauMassTools", "%s", (" vis Tau-2: Pt=" + std::to_string(preparedInput.m_vistau2.Pt()) +
996 " M=" + std::to_string(preparedInput.m_vistau2.M()) +
997 " eta=" + std::to_string(preparedInput.m_vistau2.Eta()) +
998 " phi=" + std::to_string(preparedInput.m_vistau2.Phi()) +
999 " type=" + std::to_string(preparedInput.m_type_visTau2))
1000 .c_str());
1001 Info("DiTauMassTools", "%s", (" MET=" + std::to_string(preparedInput.m_MetVec.R()) +
1002 " Met_X=" + std::to_string(preparedInput.m_MetVec.X()) +
1003 " Met_Y=" + std::to_string(preparedInput.m_MetVec.Y()))
1004 .c_str());
1005 Info("DiTauMassTools", " ---------------------------------------------------------- ");
1006 }
1007 }
1008
1009 return fit_code;
1010}
int probCalculatorV9fast(const double &phi1, const double &phi2, const double &M_nu1, const double &M_nu2)

◆ DoOutputInfo()

void MissingMassCalculator::DoOutputInfo ( )
private

Definition at line 306 of file MissingMassCalculator.cxx.

306 {
307 if (OutputInfo.m_FitStatus > 0) {
308 if (preparedInput.m_fUseVerbose == 1) {
309 Info("DiTauMassTools", "Retrieving output from fDitauStuffFit");
310 }
311 // MAXW method : get from fDittauStuffFit
312 OutputInfo.m_FitSignificance[MMCFitMethod::MAXW] = m_fDitauStuffFit.Sign_best;
313 OutputInfo.m_FittedMass[MMCFitMethod::MAXW] = m_fDitauStuffFit.Mditau_best;
314 double q1 = (1. - 0.68) / 2.;
315 double q2 = 1. - q1;
316 double xq[2], yq[2];
317 xq[0] = q1;
318 xq[1] = q2;
319 m_fMfit_all->GetQuantiles(2, yq, xq);
320 OutputInfo.m_FittedMassLowerError[MMCFitMethod::MAXW] = yq[0];
321 OutputInfo.m_FittedMassUpperError[MMCFitMethod::MAXW] = yq[1];
323 OutputInfo.m_objvec1[MMCFitMethod::MAXW] =
324 m_fDitauStuffFit.vistau1 + m_fDitauStuffFit.nutau1;
326 OutputInfo.m_objvec2[MMCFitMethod::MAXW] =
327 m_fDitauStuffFit.vistau2 + m_fDitauStuffFit.nutau2;
328 OutputInfo.m_totalvec[MMCFitMethod::MAXW] =
329 OutputInfo.m_objvec1[MMCFitMethod::MAXW] +
331 XYVector metmaxw(OutputInfo.m_nuvec1[MMCFitMethod::MAXW].Px() +
332 OutputInfo.m_nuvec2[MMCFitMethod::MAXW].Px(),
333 OutputInfo.m_nuvec1[MMCFitMethod::MAXW].Py() +
334 OutputInfo.m_nuvec2[MMCFitMethod::MAXW].Py());
335 OutputInfo.m_FittedMetVec[MMCFitMethod::MAXW] = metmaxw;
336
337 OutputInfo.m_FittedMass[MMCFitMethod::MLM] = m_fDitauStuffHisto.Mditau_best;
338 OutputInfo.m_FittedMassLowerError[MMCFitMethod::MLM] = yq[0];
339 OutputInfo.m_FittedMassUpperError[MMCFitMethod::MLM] = yq[1];
340
341 PtEtaPhiMVector tlvdummy(0., 0., 0., 0.);
342 XYVector metdummy(0., 0.);
343 OutputInfo.m_FitSignificance[MMCFitMethod::MLM] = -1.;
344 OutputInfo.m_nuvec1[MMCFitMethod::MLM] = tlvdummy;
345 OutputInfo.m_objvec1[MMCFitMethod::MLM] = tlvdummy;
346 OutputInfo.m_nuvec2[MMCFitMethod::MLM] = tlvdummy;
347 OutputInfo.m_objvec2[MMCFitMethod::MLM] = tlvdummy;
348 OutputInfo.m_totalvec[MMCFitMethod::MLM] = tlvdummy;
349 OutputInfo.m_FittedMetVec[MMCFitMethod::MLM] = metdummy;
350
351 // MLNU3P method : get from fDittauStuffHisto 4 momentum
354 m_fDitauStuffHisto.vistau1 + m_fDitauStuffHisto.nutau1;
357 m_fDitauStuffHisto.vistau2 + m_fDitauStuffHisto.nutau2;
358 OutputInfo.m_totalvec[MMCFitMethod::MLNU3P] =
361 OutputInfo.m_FittedMass[MMCFitMethod::MLNU3P] =
362 OutputInfo.m_totalvec[MMCFitMethod::MLNU3P].M();
363 OutputInfo.m_FittedMassUpperError[MMCFitMethod::MLNU3P] = 0.;
364 OutputInfo.m_FittedMassLowerError[MMCFitMethod::MLNU3P] = 0.;
365
366 XYVector metmlnu3p(OutputInfo.m_nuvec1[MMCFitMethod::MLNU3P].Px() +
367 OutputInfo.m_nuvec2[MMCFitMethod::MLNU3P].Px(),
368 OutputInfo.m_nuvec1[MMCFitMethod::MLNU3P].Py() +
369 OutputInfo.m_nuvec2[MMCFitMethod::MLNU3P].Py());
370 OutputInfo.m_FittedMetVec[MMCFitMethod::MLNU3P] = metmlnu3p;
371
372 OutputInfo.m_RMS2MPV = m_fDitauStuffHisto.RMSoverMPV;
373 }
374
375 OutputInfo.m_hMfit_all = m_fMfit_all;
376 OutputInfo.m_hMfit_allNoWeight = m_fMfit_allNoWeight;
377 OutputInfo.m_NSolutions = m_fMfit_all->GetEntries();
378 OutputInfo.m_SumW = m_fMfit_all->GetSumOfWeights();
379
380 //----------------- Check if input was re-ordered in FinalizeInputStuff() and
381 // restore the original order if needed
382 if (preparedInput.m_InputReorder == 1) {
383 PtEtaPhiMVector dummy_vec1(0.0, 0.0, 0.0, 0.0);
384 PtEtaPhiMVector dummy_vec2(0.0, 0.0, 0.0, 0.0);
385 for (int i = 0; i < 3; i++) {
386 // re-ordering neutrinos
387 dummy_vec1 = OutputInfo.m_nuvec1[i];
388 dummy_vec2 = OutputInfo.m_nuvec2[i];
389 OutputInfo.m_nuvec1[i] = dummy_vec2;
390 OutputInfo.m_nuvec2[i] = dummy_vec1;
391 // re-ordering tau's
392 dummy_vec1 = OutputInfo.m_objvec1[i];
393 dummy_vec2 = OutputInfo.m_objvec2[i];
394 OutputInfo.m_objvec1[i] = dummy_vec2;
395 OutputInfo.m_objvec2[i] = dummy_vec1;
396 }
397 }
398
399 return;
400}

◆ dTheta3DLimit()

double MissingMassCalculator::dTheta3DLimit ( const int & tau_type,
const int & limit_code,
const double & P_tau )
inline

Definition at line 2720 of file MissingMassCalculator.cxx.

2721 {
2722
2723#ifndef WITHDTHETA3DLIM
2724 // make the test ineffective if desired
2725 if (limit_code == 0)
2726 return 0.;
2727 if (limit_code == 1)
2728 return 10.;
2729 if (limit_code == 2)
2730 return 10.;
2731#endif
2732
2733 double limit = 1.0;
2734 if (limit_code == 0)
2735 limit = 0.0;
2736 double par[3] = {0.0, 0.0, 0.0};
2737 // ---- leptonic tau's
2738 if (tau_type == 8) {
2739 if (limit_code == 0) // lower 99% limit
2740 {
2741 par[0] = 0.3342;
2742 par[1] = -0.3376;
2743 par[2] = -0.001377;
2744 }
2745 if (limit_code == 1) // upper 99% limit
2746 {
2747 par[0] = 3.243;
2748 par[1] = -12.87;
2749 par[2] = 0.009656;
2750 }
2751 if (limit_code == 2) // upper 95% limit
2752 {
2753 par[0] = 2.927;
2754 par[1] = -7.911;
2755 par[2] = 0.007783;
2756 }
2757 }
2758 // ---- 1-prong tau's
2759 if (tau_type >= 0 && tau_type <= 2) {
2760 if (limit_code == 0) // lower 99% limit
2761 {
2762 par[0] = 0.2673;
2763 par[1] = -14.8;
2764 par[2] = -0.0004859;
2765 }
2766 if (limit_code == 1) // upper 99% limit
2767 {
2768 par[0] = 9.341;
2769 par[1] = -15.88;
2770 par[2] = 0.0333;
2771 }
2772 if (limit_code == 2) // upper 95% limit
2773 {
2774 par[0] = 6.535;
2775 par[1] = -8.649;
2776 par[2] = 0.00277;
2777 }
2778 }
2779 // ---- 3-prong tau's
2780 if (tau_type >= 3 && tau_type <= 5) {
2781 if (limit_code == 0) // lower 99% limit
2782 {
2783 par[0] = 0.2308;
2784 par[1] = -15.24;
2785 par[2] = -0.0009458;
2786 }
2787 if (limit_code == 1) // upper 99% limit
2788 {
2789 par[0] = 14.58;
2790 par[1] = -6.043;
2791 par[2] = -0.00928;
2792 }
2793 if (limit_code == 2) // upper 95% limit
2794 {
2795 par[0] = 8.233;
2796 par[1] = -0.3018;
2797 par[2] = -0.009399;
2798 }
2799 }
2800
2801 if (std::abs(P_tau + par[1]) > 0.0)
2802 limit = par[0] / (P_tau + par[1]) + par[2];
2803 if (limit_code == 0) {
2804 if (limit < 0.0) {
2805 limit = 0.0;
2806 } else if (limit > 0.03) {
2807 limit = 0.03;
2808 }
2809 } else {
2810 if (limit < 0.0 || limit > 0.5 * TMath::Pi()) {
2811 limit = 0.5 * TMath::Pi();
2812 } else if (limit < 0.05 && limit > 0.0) {
2813 limit = 0.05; // parameterization only runs up to P~220 GeV in this regime
2814 // will set an upper bound of 0.05
2815 }
2816 }
2817
2818 return limit;
2819}

◆ FinalizeSettings()

void MissingMassCalculator::FinalizeSettings ( const xAOD::IParticle * part1,
const xAOD::IParticle * part2,
const xAOD::MissingET * met,
const int & njets )

Definition at line 2824 of file MissingMassCalculator.cxx.

2827 {
2828 int mmcType1 = mmcType(part1);
2829 if (mmcType1 < 0)
2830 return; // return CP::CorrectionCode::Error;
2831
2832 int mmcType2 = mmcType(part2);
2833 if (mmcType2 < 0)
2834 return; // return CP::CorrectionCode::Error;
2835
2836 preparedInput.SetLFVmode(-2); // initialise LFV mode value for this event with being *not* LFV
2837 // if(getLFVMode(part1, part2, mmcType1, mmcType2) ==
2838 // CP::CorrectionCode::Error) {
2840 int LFVMode = getLFVMode(part1, part2, mmcType1, mmcType2);
2841 if (LFVMode == -1) {
2842 return; // return CP::CorrectionCode::Error;
2843 } else if (LFVMode != -2) {
2844 preparedInput.SetLFVmode(LFVMode);
2845 }
2846 }
2847
2848 // this will be in MeV but MMC allows MeV
2849 // assume the mass is correct as well
2850 PtEtaPhiMVector tlvTau1(part1->pt(), part1->eta(), part1->phi(), part1->m());
2851 PtEtaPhiMVector tlvTau2(part2->pt(), part2->eta(), part2->phi(), part2->m());
2852
2853 // Convert to GeV. In principle, MMC should cope with MeV but should check
2854 // thoroughly
2855 PtEtaPhiMVector fixedtau1;
2856 fixedtau1.SetCoordinates(tlvTau1.Pt() / GEV, tlvTau1.Eta(), tlvTau1.Phi(), tlvTau1.M() / GEV);
2857 PtEtaPhiMVector fixedtau2;
2858 fixedtau2.SetCoordinates(tlvTau2.Pt() / GEV, tlvTau2.Eta(), tlvTau2.Phi(), tlvTau2.M() / GEV);
2859
2860 preparedInput.SetVisTauType(0, mmcType1);
2861 preparedInput.SetVisTauType(1, mmcType2);
2862 preparedInput.SetVisTauVec(0, fixedtau1);
2863 preparedInput.SetVisTauVec(1, fixedtau2);
2864
2865 if (mmcType1 == 8 && mmcType2 == 8) {
2866 preparedInput.m_tauTypes = TauTypes::ll;
2867 } else if (mmcType1 >= 0 && mmcType1 <= 5 && mmcType2 >= 0 && mmcType2 <= 5) {
2868 preparedInput.m_tauTypes = TauTypes::hh;
2869 } else {
2870 preparedInput.m_tauTypes = TauTypes::lh;
2871 }
2872 if (preparedInput.m_fUseVerbose)
2873 Info("DiTauMassTools", "%s", ("running for tau types "+std::to_string(preparedInput.m_type_visTau1)+" "+std::to_string(preparedInput.m_type_visTau2)).c_str());
2874 XYVector met_vec(met->mpx() / GEV, met->mpy() / GEV);
2875 preparedInput.SetMetVec(met_vec);
2876 if (preparedInput.m_fUseVerbose)
2877 Info("DiTauMassTools", "%s", ("passing SumEt="+std::to_string(met->sumet() / GEV)).c_str());
2878 preparedInput.SetSumEt(met->sumet() / GEV);
2879 preparedInput.SetNjet25(njets);
2880
2881 // check that the calibration set has been chosen explicitly, otherwise abort
2883 Error("DiTauMassTools", "MMCCalibrationSet has not been set !. Please use "
2884 "fMMC.SetCalibrationSet(MMCCalibrationSet::MMC2019) or fMMC.SetCalibrationSet(MMCCalibrationSet::MMC2024)"
2885 ". Abort now. ");
2886 std::abort();
2887 }
2888 //----------- Re-ordering input info, to make sure there is no dependence of
2889 // results on input order
2890 // this might be needed because a random scan is used
2891 // highest pT tau is always first
2892 preparedInput.m_InputReorder = 0; // set flag to 0 by default, i.e. no re-ordering
2893 if ((preparedInput.m_type_visTau1 >= 0 && preparedInput.m_type_visTau1 <= 5) &&
2894 preparedInput.m_type_visTau2 == 8) // if hadron-lepton, reorder to have lepton first
2895 {
2896 preparedInput.m_InputReorder =
2897 1; // re-order to be done, this flag is to be checked in DoOutputInfo()
2898 } else if (!((preparedInput.m_type_visTau2 >= 0 && preparedInput.m_type_visTau2 <= 5) &&
2899 preparedInput.m_type_visTau1 == 8)) // if not lep-had nor had lep, reorder if tau1 is
2900 // after tau2 clockwise
2901 {
2902 if (fixPhiRange(preparedInput.m_vistau1.Phi() - preparedInput.m_vistau2.Phi()) > 0) {
2903 preparedInput.m_InputReorder = 1; // re-order to be done, this flag is to be
2904 // checked in DoOutputInfo()
2905 }
2906 }
2907
2908 if (preparedInput.m_InputReorder == 1) // copy and re-order
2909 {
2910 std::swap(preparedInput.m_vistau1, preparedInput.m_vistau2);
2911 std::swap(preparedInput.m_type_visTau1, preparedInput.m_type_visTau2);
2912 std::swap(preparedInput.m_Nprong_tau1, preparedInput.m_Nprong_tau2);
2913 }
2914 //--------- re-ordering is done ---------------------------------------
2915
2916 preparedInput.m_DelPhiTT =
2917 std::abs(Phi_mpi_pi(preparedInput.m_vistau1.Phi() - preparedInput.m_vistau2.Phi()));
2918
2919 for (unsigned int i = 0; i < preparedInput.m_jet4vecs.size(); i++) {
2920 // correcting sumEt, give priority to SetMetScanParamsUE()
2921 if (preparedInput.m_METScanScheme == 0) {
2922 if ((preparedInput.m_METsigmaP < 0.1 || preparedInput.m_METsigmaL < 0.1) &&
2923 preparedInput.m_SumEt > preparedInput.m_jet4vecs[i].Pt() &&
2924 preparedInput.m_jet4vecs[i].Pt() > 20.0) {
2925 if (preparedInput.m_fUseVerbose == 1) {
2926 Info("DiTauMassTools", "correcting sumET");
2927 }
2928 preparedInput.m_SumEt -= preparedInput.m_jet4vecs[i].Pt();
2929 }
2930 }
2931 }
2932
2933 // give priority to SetVisTauType, only do this if type_visTau1 and
2934 // type_visTau2 are not set
2935 /*if(type_visTau1<0 && type_visTau2<0 && Nprong_tau1>-1 && Nprong_tau2>-1)
2936 {
2937 if(Nprong_tau1==0) type_visTau1 = 8; // leptonic tau
2938 else if( Nprong_tau1==1) type_visTau1 = 0; // set to 1p0n for now, may use
2939different solution later like explicit integer for this case that pantau info is
2940not available? else if( Nprong_tau1==3) type_visTau1 = 3; // set to 3p0n for
2941now, see above if(Nprong_tau2==0) type_visTau2 = 8; // leptonic tau else if(
2942Nprong_tau2==1) type_visTau2 = 0; // set to 1p0n for now, see above else if(
2943Nprong_tau2==3) type_visTau2=3; // set to 3p0n for now, see above
2944 }
2945 */
2946 // checking input mass of hadronic tau-1
2947 // one prong
2948 // // checking input mass of hadronic tau-1
2949 // DRMERGE LFV addition
2951 if ((preparedInput.m_type_visTau1 >= 0 && preparedInput.m_type_visTau1 <= 2) &&
2952 preparedInput.m_vistau1.M() != 1.1) {
2953 preparedInput.m_vistau1.SetCoordinates(preparedInput.m_vistau1.Pt(), preparedInput.m_vistau1.Eta(),
2954 preparedInput.m_vistau1.Phi(), 1.1);
2955 }
2956 if ((preparedInput.m_type_visTau1 >= 3 && preparedInput.m_type_visTau1 <= 5) &&
2957 preparedInput.m_vistau1.M() != 1.35) {
2958 preparedInput.m_vistau1.SetCoordinates(preparedInput.m_vistau1.Pt(), preparedInput.m_vistau1.Eta(),
2959 preparedInput.m_vistau1.Phi(), 1.35);
2960 }
2961 // checking input mass of hadronic tau-2
2962 if ((preparedInput.m_type_visTau2 >= 0 && preparedInput.m_type_visTau2 <= 2) &&
2963 preparedInput.m_vistau2.M() != 1.1) {
2964 preparedInput.m_vistau2.SetCoordinates(preparedInput.m_vistau2.Pt(), preparedInput.m_vistau2.Eta(),
2965 preparedInput.m_vistau2.Phi(), 1.1);
2966 }
2967 if ((preparedInput.m_type_visTau2 >= 3 && preparedInput.m_type_visTau2 <= 5) &&
2968 preparedInput.m_vistau2.M() != 1.35) {
2969 preparedInput.m_vistau2.SetCoordinates(preparedInput.m_vistau2.Pt(), preparedInput.m_vistau2.Eta(),
2970 preparedInput.m_vistau2.Phi(), 1.35);
2971 }
2972 } else {
2973 // DRMERGE end LFV addition
2974 if ((preparedInput.m_type_visTau1 >= 0 && preparedInput.m_type_visTau1 <= 2) &&
2975 preparedInput.m_vistau1.M() != 0.8) {
2976 preparedInput.m_vistau1.SetCoordinates(preparedInput.m_vistau1.Pt(), preparedInput.m_vistau1.Eta(),
2977 preparedInput.m_vistau1.Phi(), 0.8);
2978 }
2979 // 3 prong
2980 if ((preparedInput.m_type_visTau1 >= 3 && preparedInput.m_type_visTau1 <= 5) &&
2981 preparedInput.m_vistau1.M() != 1.2) {
2982 preparedInput.m_vistau1.SetCoordinates(preparedInput.m_vistau1.Pt(), preparedInput.m_vistau1.Eta(),
2983 preparedInput.m_vistau1.Phi(), 1.2);
2984 }
2985 // checking input mass of hadronic tau-2
2986 // one prong
2987 if ((preparedInput.m_type_visTau2 >= 0 && preparedInput.m_type_visTau2 <= 2) &&
2988 preparedInput.m_vistau2.M() != 0.8) {
2989 preparedInput.m_vistau2.SetCoordinates(preparedInput.m_vistau2.Pt(), preparedInput.m_vistau2.Eta(),
2990 preparedInput.m_vistau2.Phi(), 0.8);
2991 }
2992 // 3 prong
2993 if ((preparedInput.m_type_visTau2 >= 3 && preparedInput.m_type_visTau2 <= 5) &&
2994 preparedInput.m_vistau2.M() != 1.2) {
2995 preparedInput.m_vistau2.SetCoordinates(preparedInput.m_vistau2.Pt(), preparedInput.m_vistau2.Eta(),
2996 preparedInput.m_vistau2.Phi(), 1.2);
2997 }
2998 } // DRDRMERGE LFV else closing
2999
3000 // correcting sumEt for electron pt, give priority to SetMetScanParamsUE()
3001 // DR20150615 in tag 00-00-11 and before. The following was done before the
3002 // mass of the hadronic tau was set which mean that sumEt was wrongly
3003 // corrected for the hadronic tau pt if the hadronic tau mass was set to zero
3004 // Sasha 08/12/15: don't do electron Pt subtraction for high mass studies; in
3005 // the future, need to check if lepton Pt needs to be subtracted for both ele
3006 // and muon
3007 if (preparedInput.m_METsigmaP < 0.1 || preparedInput.m_METsigmaL < 0.1) {
3008
3009 // T. Davidek: hack for lep-lep -- subtract lepton pT both for muon and
3010 // electron
3013 preparedInput.m_vistau1.M() < 0.12 && preparedInput.m_vistau2.M() < 0.12) { // lep-lep channel
3014 if (preparedInput.m_SumEt > preparedInput.m_vistau1.Pt())
3015 preparedInput.m_SumEt -= preparedInput.m_vistau1.Pt();
3016 if (preparedInput.m_SumEt > preparedInput.m_vistau2.Pt())
3017 preparedInput.m_SumEt -= preparedInput.m_vistau2.Pt();
3018 } else {
3019 // continue with the original code
3020 if (preparedInput.m_SumEt > preparedInput.m_vistau1.Pt() && preparedInput.m_vistau1.M() < 0.05 &&
3022 if (preparedInput.m_fUseVerbose == 1) {
3023 Info("DiTauMassTools", "Substracting pt1 from sumEt");
3024 }
3025 preparedInput.m_SumEt -= preparedInput.m_vistau1.Pt();
3026 }
3027 if (preparedInput.m_SumEt > preparedInput.m_vistau2.Pt() && preparedInput.m_vistau2.M() < 0.05 &&
3029 if (preparedInput.m_fUseVerbose == 1) {
3030 Info("DiTauMassTools", "Substracting pt2 from sumEt");
3031 }
3032 preparedInput.m_SumEt -= preparedInput.m_vistau2.Pt();
3033 }
3034 }
3035 }
3036
3037 // controling TauProbability settings for UPGRADE studies
3039 preparedInput.m_fUseDefaults == 1) {
3040 if ((preparedInput.m_vistau1.M() < 0.12 && preparedInput.m_vistau2.M() > 0.12) ||
3041 (preparedInput.m_vistau2.M() < 0.12 && preparedInput.m_vistau1.M() > 0.12)) {
3042 Prob->SetUseTauProbability(true); // lep-had case
3043 }
3044 if (preparedInput.m_vistau1.M() > 0.12 && preparedInput.m_vistau2.M() > 0.12) {
3045 Prob->SetUseTauProbability(false); // had-had case
3046 }
3047 }
3048
3049 // change Beam Energy for different running conditions
3050 preparedInput.m_beamEnergy = m_beamEnergy;
3051
3052 //--------------------- pre-set defaults for Run-2. To disable pre-set
3053 // defaults set fUseDefaults=0
3054 if (preparedInput.m_fUseDefaults == 1) {
3059 preparedInput.m_fUseTailCleanup = 0;
3060 if ((preparedInput.m_vistau1.M() < 0.12 && preparedInput.m_vistau2.M() > 0.12) ||
3061 (preparedInput.m_vistau2.M() < 0.12 && preparedInput.m_vistau1.M() > 0.12))
3062 Prob->SetUseTauProbability(false); // lep-had
3063 if (preparedInput.m_tauTypes == TauTypes::hh)
3064 Prob->SetUseTauProbability(true); // had-had
3065 Prob->SetUseMnuProbability(false);
3066 }
3067 }
3068
3069 // compute HTOffset if relevant
3070 if (Prob->GetUseHT()) // use missing Ht for Njet25=0 events
3071 {
3072 // dPhi(l-t) dependence of misHt-trueMET
3073 double HtOffset = 0.;
3074 // proper for hh
3075 if (preparedInput.m_tauTypes == TauTypes::hh) {
3076 // hh
3077 double x = preparedInput.m_DelPhiTT;
3078 HtOffset = 87.5 - 27.0 * x;
3079 }
3080
3081 preparedInput.m_HtOffset = HtOffset;
3082
3083 // if use HT, replace MET with HT
3084 preparedInput.m_METsigmaP =
3085 preparedInput.m_MHtSigma2; // sigma of 2nd Gaussian for missing Ht resolution
3086 preparedInput.m_METsigmaL = preparedInput.m_MHtSigma2;
3087
3088 PtEtaPhiMVector tauSum = preparedInput.m_vistau1 + preparedInput.m_vistau2;
3089 preparedInput.m_MetVec.SetXY(-tauSum.Px(), -tauSum.Py()); // WARNING this replace metvec by -mht
3090 }
3091}
__HOSTDEV__ double Phi_mpi_pi(double)
Definition GeoRegion.cxx:10
#define x
virtual double eta() const =0
The pseudorapidity ( ) of the particle.
virtual double pt() const =0
The transverse momentum ( ) of the particle.
virtual double m() const =0
The invariant mass of the particle.
virtual double phi() const =0
The azimuthal angle ( ) of the particle.
float sumet() const
Returns.
float mpx() const
Returns .
float mpy() const
Returns .
int getLFVMode(const xAOD::IParticle *p1, const xAOD::IParticle *p2, int mmcType1, int mmcType2)
Error
The different types of error that can be flagged in the L1TopoRDO.
Definition Error.h:16
void swap(ElementLinkVector< DOBJ > &lhs, ElementLinkVector< DOBJ > &rhs)

◆ GetMarkovCountDuplicate()

int DiTauMassTools::MissingMassCalculator::GetMarkovCountDuplicate ( ) const
inline

◆ GetMarkovNAccept()

int DiTauMassTools::MissingMassCalculator::GetMarkovNAccept ( ) const
inline

Definition at line 383 of file MissingMassCalculator.h.

383{ return m_markovNAccept; }

◆ GetMarkovNFullscan()

int DiTauMassTools::MissingMassCalculator::GetMarkovNFullscan ( ) const
inline

Definition at line 384 of file MissingMassCalculator.h.

384{ return m_markovNFullScan;}

◆ GetMarkovNRejectMetropolis()

int DiTauMassTools::MissingMassCalculator::GetMarkovNRejectMetropolis ( ) const
inline

Definition at line 382 of file MissingMassCalculator.h.

◆ GetMarkovNRejectNoSol()

int DiTauMassTools::MissingMassCalculator::GetMarkovNRejectNoSol ( ) const
inline

Definition at line 381 of file MissingMassCalculator.h.

381{ return m_markovNRejectNoSol;}

◆ GetMeanbinStop()

double DiTauMassTools::MissingMassCalculator::GetMeanbinStop ( ) const
inline

Definition at line 377 of file MissingMassCalculator.h.

377{ return m_meanbinStop;}

◆ GetmInvWidth2Error()

double DiTauMassTools::MissingMassCalculator::GetmInvWidth2Error ( ) const
inline

◆ GetmMaxError()

double DiTauMassTools::MissingMassCalculator::GetmMaxError ( ) const
inline

◆ GetmMeanError()

double DiTauMassTools::MissingMassCalculator::GetmMeanError ( ) const
inline

◆ GetNiterFit1()

int DiTauMassTools::MissingMassCalculator::GetNiterFit1 ( ) const
inline

Definition at line 370 of file MissingMassCalculator.h.

370{ return m_niter_fit1; } // number of iterations per loop in dPhi loop

◆ GetNiterFit2()

int DiTauMassTools::MissingMassCalculator::GetNiterFit2 ( ) const
inline

Definition at line 371 of file MissingMassCalculator.h.

371{ return m_niter_fit2; } // number of iterations per loop in MET loop

◆ GetNiterFit3()

int DiTauMassTools::MissingMassCalculator::GetNiterFit3 ( ) const
inline

Definition at line 372 of file MissingMassCalculator.h.

372{ return m_niter_fit3; } // number of iterations per loop in Mnu loop

◆ GetNiterRandom()

int DiTauMassTools::MissingMassCalculator::GetNiterRandom ( ) const
inline

Definition at line 373 of file MissingMassCalculator.h.

373{ return m_niterRandomLocal; } // number of random iterations

◆ GetNMetroReject()

int DiTauMassTools::MissingMassCalculator::GetNMetroReject ( ) const
inline

Definition at line 408 of file MissingMassCalculator.h.

◆ GetNNoSol()

int DiTauMassTools::MissingMassCalculator::GetNNoSol ( ) const
inline

Definition at line 407 of file MissingMassCalculator.h.

407{return m_markovNRejectNoSol;}

◆ GetNSol()

int DiTauMassTools::MissingMassCalculator::GetNSol ( ) const
inline

Definition at line 409 of file MissingMassCalculator.h.

409{return m_markovNAccept;}

◆ GetNsucStop()

int DiTauMassTools::MissingMassCalculator::GetNsucStop ( ) const
inline

Definition at line 375 of file MissingMassCalculator.h.

375{ return m_NsucStop; } // Arrest criteria for NSuc

◆ GetProposalTryEtau()

double DiTauMassTools::MissingMassCalculator::GetProposalTryEtau ( ) const
inline

Definition at line 388 of file MissingMassCalculator.h.

388{return m_ProposalTryEtau;}

◆ GetProposalTryMEt()

double DiTauMassTools::MissingMassCalculator::GetProposalTryMEt ( ) const
inline

Definition at line 385 of file MissingMassCalculator.h.

385{return m_proposalTryMEt;}

◆ GetProposalTryMnu()

double DiTauMassTools::MissingMassCalculator::GetProposalTryMnu ( ) const
inline

Definition at line 387 of file MissingMassCalculator.h.

387{return m_ProposalTryMnu;}

◆ GetProposalTryPhi()

double DiTauMassTools::MissingMassCalculator::GetProposalTryPhi ( ) const
inline

Definition at line 386 of file MissingMassCalculator.h.

386{return m_ProposalTryPhi;}

◆ GetRMSStop()

int DiTauMassTools::MissingMassCalculator::GetRMSStop ( ) const
inline

Definition at line 376 of file MissingMassCalculator.h.

376{ return m_RMSStop; }

◆ GetRndmSeedAltering()

int DiTauMassTools::MissingMassCalculator::GetRndmSeedAltering ( ) const
inline

Definition at line 378 of file MissingMassCalculator.h.

378{ return m_RndmSeedAltering; } // number of iterations per loop in Mnu loop

◆ GetUseEfficiencyRecovery()

bool DiTauMassTools::MissingMassCalculator::GetUseEfficiencyRecovery ( ) const
inline

Definition at line 368 of file MissingMassCalculator.h.

368{ return m_fUseEfficiencyRecovery; }

◆ handleSolutions()

void MissingMassCalculator::handleSolutions ( )
inlineprotected

Definition at line 2045 of file MissingMassCalculator.cxx.

2047{
2048
2049 bool reject = true;
2050 double totalProbSumSol = 0.;
2051 double totalProbSumSolOld = 0.;
2052 bool firstPointWithSol = false;
2053
2054 for (int isol = 0; isol < m_nsol; ++isol) {
2055 totalProbSumSol += m_probFinalSolVec[isol];
2056 }
2057
2058 double uMC = -1.;
2059 bool notSureToKeep = true;
2060 // note : if no solution, the point is treated as having a zero probability
2062 reject = false; // accept anyway in this mode
2063 notSureToKeep = false; // do not need to test on prob
2064 if (m_nsol <= 0) {
2065 // if initial full scaning and no sol : continue
2066 m_markovNFullScan += 1;
2067 } else {
2068 // if we were in in full scan mode and we have a solution, switch it off
2069 m_fullParamSpaceScan = false;
2070 firstPointWithSol = true; // as this is the first point without a solution
2071 // there is no old sol
2072 m_iter0 = 0; // reset the counter so that separately the full scan pphase
2073 // and the markov phase use m_niterRandomLocal points
2074 // hack for hh : allow 10 times less iteration for markov than for the
2075 // fullscan phase
2076 if (preparedInput.m_tauTypes == TauTypes::hh) {
2077 m_niterRandomLocal /= 10;
2078 }
2079 }
2080 }
2081
2082 if (notSureToKeep) {
2083 // apply Metropolis algorithm to decide to keep this point.
2084 // compute the probability of the previous point and the current one
2085 for (int isol = 0; isol < m_nsolOld; ++isol) {
2086 totalProbSumSolOld += m_probFinalSolOldVec[isol];
2087 }
2088
2089 // accept anyway if null old probability (should only happen for the very
2090 // first point with a solution)
2091 if (!firstPointWithSol && totalProbSumSolOld <= 0.) {
2092 Error("DiTauMassTools", "%s",
2093 (" ERROR null old probability !!! " + std::to_string(totalProbSumSolOld) + " nsolOld " +
2094 std::to_string(m_nsolOld))
2095 .c_str());
2096 reject = false;
2097 } else if (totalProbSumSol > totalProbSumSolOld) {
2098 // if going up, accept anyway
2099 reject = false;
2100 // else if (totalProbSumSol < 1E-16) { // if null target probability,
2101 // reject anyway
2102 } else if (totalProbSumSol < totalProbSumSolOld * 1E-6) { // if ratio of probability <1e6, point
2103 // will be accepted only every 1E6
2104 // iteration, so can reject anyway
2105 reject = true;
2106 } else if (m_nsol <= 0) { // new parametrisation give prob too small to
2107 // trigger above condition if no solution is found
2108 reject = true;
2109 } else {
2110 // if going down, reject with a probability
2111 // 1-totalProbSum/totalProbSumOld)
2112 uMC = m_randomGen.Rndm();
2113 reject = (uMC > totalProbSumSol / totalProbSumSolOld);
2114 }
2115 } // if reject
2116
2117 // proceed with the handling of the solutions wether the old or the new ones
2118
2119 // optionally fill the vectors with the complete list of points (for all
2120 // walkstrategy)
2121
2122 if (reject) {
2123 // current point reset to the previous one
2124 // Note : only place where m_MEtP etc... are modified outside spacewalkerXYZ
2125 m_MEtP = m_MEtP0;
2126 m_MEtL = m_MEtL0;
2127 m_Phi1 = m_Phi10;
2128 m_Phi2 = m_Phi20;
2129 m_eTau1 = m_eTau10;
2130 m_eTau2 = m_eTau20;
2131 if (m_scanMnu1)
2132 m_Mnu1 = m_Mnu10;
2133 if (m_scanMnu2)
2134 m_Mnu2 = m_Mnu20;
2135 }
2136
2137 // default case : fill the histogram with solution, using current point
2138 bool fillSolution = true;
2139 bool oldToBeUsed = false;
2140
2141 // now handle the reject or accept cases
2142 // the tricky thing is that for markov, we accept the old point as soon as a
2143 // new accepted point is found with a weight equal to one plus the number of
2144 // rejected point inbetween
2145
2146 if (reject) {
2147 fillSolution = false; // do not fill solution, just count number of replication
2149 if (m_nsol <= 0) {
2151 } else {
2153 }
2154
2155 } else {
2156 // if accept, will fill solution (except for very first point) but taking
2157 // the values from the previous point
2158 if (!m_fullParamSpaceScan) {
2159 m_markovNAccept += 1;
2160 }
2161 if (!firstPointWithSol) {
2162 fillSolution = true;
2163 oldToBeUsed = true;
2164 } else {
2165 fillSolution = false;
2166 }
2167 } // else reject
2168
2169 // if do not fill solution exit now
2170 // for the first point with solution we need to copy the new sol into the old
2171 // one before leaving
2172 if (!fillSolution) {
2173 if (firstPointWithSol) {
2174 // current point is the future previous one
2175 m_nsolOld = m_nsol;
2176 for (int isol = 0; isol < m_nsol; ++isol) {
2181 }
2182 }
2183 return;
2184 }
2185
2186 // compute RMS of the different solutions
2187 double solSum = 0.;
2188 double solSum2 = 0.;
2189
2190 for (int isol = 0; isol < m_nsol; ++isol) {
2191 ++m_iter5;
2192 double totalProb;
2193 double mtautau;
2194 const PtEtaPhiMVector *pnuvec1_tmpj;
2195 const PtEtaPhiMVector *pnuvec2_tmpj;
2196
2197 if (oldToBeUsed) {
2198 totalProb = m_probFinalSolOldVec[isol];
2199 mtautau = m_mtautauFinalSolOldVec[isol];
2200 pnuvec1_tmpj = &m_nu1FinalSolOldVec[isol];
2201 pnuvec2_tmpj = &m_nu2FinalSolOldVec[isol];
2202 } else {
2203 totalProb = m_probFinalSolVec[isol];
2204 mtautau = m_mtautauFinalSolVec[isol];
2205 pnuvec1_tmpj = &m_nu1FinalSolVec[isol];
2206 pnuvec2_tmpj = &m_nu2FinalSolVec[isol];
2207 }
2208 const PtEtaPhiMVector &nuvec1_tmpj = *pnuvec1_tmpj;
2209 const PtEtaPhiMVector &nuvec2_tmpj = *pnuvec2_tmpj;
2210
2211 solSum += mtautau;
2212 solSum2 += mtautau * mtautau;
2213
2214 double weight;
2215 // MarkovChain : accepted events already distributed according to
2216 // probability distribution, so weight is 1. acutally to have a proper
2217 // estimate of per bin error, instead of putting several time the same point
2218 // when metropolis alg reject one (or no solution), rather put it with the
2219 // multiplicity weight. Should only change the error bars might change if
2220 // weighted markov chain are used there is also an issue with the 4 very
2221 // close nearly identical solution
2223 1; // incremented only when a point is rejected, hence need to add 1
2224
2225 m_fMfit_all->Fill(mtautau, weight);
2226
2227 if(m_SaveLlhHisto){
2228 m_fMEtP_all->Fill(m_MEtP, weight);
2229 m_fMEtL_all->Fill(m_MEtL, weight);
2230 m_fMnu1_all->Fill(m_Mnu1, weight);
2231 m_fMnu2_all->Fill(m_Mnu2, weight);
2232 m_fPhi1_all->Fill(m_Phi1, weight);
2233 m_fPhi2_all->Fill(m_Phi2, weight);
2234 if (mtautau != 0. && weight != 0.)
2235 m_fMfit_allGraph->SetPoint(m_iter0, mtautau, -TMath::Log(weight));
2236 }
2237
2238 m_fMfit_allNoWeight->Fill(mtautau, 1.);
2239
2240 // m_fPXfit1->Fill(nuvec1_tmpj.Px(),weight);
2241 // m_fPYfit1->Fill(nuvec1_tmpj.Py(),weight);
2242 // m_fPZfit1->Fill(nuvec1_tmpj.Pz(),weight);
2243 // m_fPXfit2->Fill(nuvec2_tmpj.Px(),weight);
2244 // m_fPYfit2->Fill(nuvec2_tmpj.Py(),weight);
2245 // m_fPZfit2->Fill(nuvec2_tmpj.Pz(),weight);
2246
2247 //----------------- using P*fit to fill Px,y,z_tau
2248 // Note that the original vistau are used there deliberately,
2249 // since they will be subtracted after histogram fitting
2250 // DR, kudos Antony Lesage : do not create temporary TLV within each Fill,
2251 // saves 10% CPU
2252 m_fPXfit1->Fill(preparedInput.m_vistau1.Px() + nuvec1_tmpj.Px(), totalProb);
2253 m_fPYfit1->Fill(preparedInput.m_vistau1.Py() + nuvec1_tmpj.Py(), totalProb);
2254 m_fPZfit1->Fill(preparedInput.m_vistau1.Pz() + nuvec1_tmpj.Pz(), totalProb);
2255 m_fPXfit2->Fill(preparedInput.m_vistau2.Px() + nuvec2_tmpj.Px(), totalProb);
2256 m_fPYfit2->Fill(preparedInput.m_vistau2.Py() + nuvec2_tmpj.Py(), totalProb);
2257 m_fPZfit2->Fill(preparedInput.m_vistau2.Pz() + nuvec2_tmpj.Pz(), totalProb);
2258
2259 // fill histograms for floating stopping criterion, split randomly
2260 if (m_fUseFloatStopping) {
2261 if (m_randomGen.Rndm() <= 0.5) {
2262 m_fMmass_split1->Fill(mtautau, weight);
2263 m_fMEtP_split1->Fill(m_MEtP, weight);
2264 m_fMEtL_split1->Fill(m_MEtL, weight);
2265 m_fMnu1_split1->Fill(m_Mnu1, weight);
2266 m_fMnu2_split1->Fill(m_Mnu2, weight);
2267 m_fPhi1_split1->Fill(m_Phi1, weight);
2268 m_fPhi2_split1->Fill(m_Phi2, weight);
2269 } else {
2270 m_fMmass_split2->Fill(mtautau, weight);
2271 m_fMEtP_split2->Fill(m_MEtP, weight);
2272 m_fMEtL_split2->Fill(m_MEtL, weight);
2273 m_fMnu1_split2->Fill(m_Mnu1, weight);
2274 m_fMnu2_split2->Fill(m_Mnu2, weight);
2275 m_fPhi1_split2->Fill(m_Phi1, weight);
2276 m_fPhi2_split2->Fill(m_Phi2, weight);
2277 }
2278 }
2279
2280 if (totalProb > m_prob_tmp) // fill solution with highest probability
2281 {
2282 m_prob_tmp = totalProb;
2283 m_fDitauStuffFit.Mditau_best = mtautau;
2284 m_fDitauStuffFit.Sign_best = -log10(totalProb);
2285 ;
2286 m_fDitauStuffFit.nutau1 = nuvec1_tmpj;
2287 m_fDitauStuffFit.nutau2 = nuvec2_tmpj;
2288 m_fDitauStuffFit.vistau1 = m_tauVec1;
2289 m_fDitauStuffFit.vistau2 = m_tauVec2;
2290 }
2291 } // loop on solutions
2292
2293 m_markovCountDuplicate = 0; // now can reset the duplicate count
2294
2295 if (oldToBeUsed) {
2296 // current point is the future previous one
2297 // TLV copy not super efficient but not dramatic
2298 m_nsolOld = m_nsol;
2299 for (int isol = 0; isol < m_nsol; ++isol) {
2304 }
2305 }
2306
2307 // compute rms of solutions
2308 const double solRMS = sqrt(solSum2 / m_nsol - std::pow(solSum / m_nsol, 2));
2309 OutputInfo.m_AveSolRMS += solRMS;
2310
2311 return;
2312}

◆ MassScale()

double DiTauMassTools::MissingMassCalculator::MassScale ( int method,
double mass,
const int & tau_type1,
const int & tau_type2 )
inlineprotected

◆ maxFitting()

Double_t MissingMassCalculator::maxFitting ( Double_t * x,
Double_t * par )

Definition at line 1497 of file MissingMassCalculator.cxx.

1499{
1500 // parabola with parameters max, mean and invwidth
1501 const double mM = x[0];
1502 const double mMax = par[0];
1503 const double mMean = par[1];
1504 const double mInvWidth2 = par[2]; // if param positif distance between intersection of the
1505 // parabola with x axis: 1/Sqrt(mInvWidth2)
1506 const double fitval = mMax * (1 - 4 * mInvWidth2 * std::pow(mM - mMean, 2));
1507 return fitval;
1508}

◆ maxFromHist() [1/2]

double DiTauMassTools::MissingMassCalculator::maxFromHist ( const std::shared_ptr< TH1F > & theHist,
std::vector< double > & histInfo,
const MaxHistStrategy::e maxHistStrategy = MaxHistStrategy::FIT,
const int winHalfWidth = 2,
bool debug = false )
inline

Definition at line 417 of file MissingMassCalculator.h.

417 {
418 return maxFromHist(theHist.get(), histInfo, maxHistStrategy, winHalfWidth, debug);
419 }
const bool debug

◆ maxFromHist() [2/2]

double MissingMassCalculator::maxFromHist ( TH1F * theHist,
std::vector< double > & histInfo,
const MaxHistStrategy::e maxHistStrategy = MaxHistStrategy::FIT,
const int winHalfWidth = 2,
bool debug = false )

Definition at line 1517 of file MissingMassCalculator.cxx.

1519 {
1520 // namespace HistInfo
1521 // enum e {
1522 // PROB=0,INTEGRAL,CHI2,DISCRI,TANTHETA,TANTHETAW,FITLENGTH,RMS,RMSVSDISCRI,MAXHISTINFO
1523 // };
1524 double maxPos = 0.;
1525 double prob = 0.;
1526
1527 for (std::vector<double>::iterator itr = histInfo.begin(); itr != histInfo.end(); ++itr) {
1528 *itr = -1;
1529 }
1530
1531 histInfo[HistInfo::INTEGRAL] = theHist->Integral();
1532
1533 if (maxHistStrategy == MaxHistStrategy::MAXBIN ||
1534 ((maxHistStrategy == MaxHistStrategy::MAXBINWINDOW ||
1535 maxHistStrategy == MaxHistStrategy::SLIDINGWINDOW) &&
1536 winHalfWidth == 0)) {
1537
1538 // simple max search
1539 // original version, simple bin maximum
1540 int max_bin = theHist->GetMaximumBin();
1541 maxPos = theHist->GetBinCenter(max_bin);
1542
1543 // FIXME GetEntries is unweighted
1544 prob = theHist->GetBinContent(max_bin) / double(theHist->GetEntries());
1545 if (prob > 1.)
1546 prob = 1.;
1547 histInfo[HistInfo::PROB] = prob;
1548 return maxPos;
1549 }
1550
1551 int hNbins = theHist->GetNbinsX();
1552
1553 if (maxHistStrategy == MaxHistStrategy::MAXBINWINDOW) {
1554 // average around maximum bin (nearly useless in fact)
1555 // could be faster
1556 int max_bin = theHist->GetMaximumBin();
1557 int iBinMin = max_bin - winHalfWidth;
1558 if (iBinMin < 0)
1559 iBinMin = 0;
1560 int iBinMax = max_bin + winHalfWidth;
1561 if (iBinMax > hNbins)
1562 iBinMax = hNbins - 1;
1563 double sumw = 0;
1564 double sumx = 0;
1565 for (int iBin = iBinMin; iBin <= iBinMax; ++iBin) {
1566 const double weight = theHist->GetBinContent(iBin);
1567 sumw += weight;
1568 sumx += weight * theHist->GetBinCenter(iBin);
1569 }
1570 maxPos = sumx / sumw;
1571
1572 // FIXME GetEntries is unweighted
1573 prob = sumw / theHist->GetEntries();
1574 if (prob > 1.)
1575 prob = 1.;
1576
1577 return maxPos;
1578 }
1579
1580 // now compute sliding window anyway
1581 if (maxHistStrategy != MaxHistStrategy::SLIDINGWINDOW &&
1582 maxHistStrategy != MaxHistStrategy::FIT) {
1583 Error("DiTauMassTools", "%s",
1584 ("ERROR undefined maxHistStrategy:" + std::to_string(maxHistStrategy)).c_str());
1585 return -10.;
1586 }
1587
1588 // first iteration to find the first and last non zero bin, and the histogram
1589 // integral (not same as Entries because of weights)
1590 int lastNonZeroBin = -1;
1591 int firstNonZeroBin = -1;
1592 double totalSumw = 0.;
1593 bool firstNullPart = true;
1594 for (int iBin = 0; iBin < hNbins; ++iBin) {
1595 const double weight = theHist->GetBinContent(iBin);
1596 if (weight > 0) {
1597 totalSumw += weight;
1598 lastNonZeroBin = iBin;
1599 if (firstNullPart) {
1600 firstNullPart = false;
1601 firstNonZeroBin = iBin;
1602 }
1603 }
1604 }
1605
1606 // enlarge first and last non zero bin with window width to avoid side effect
1607 // (maximum close to the edge)
1608 firstNonZeroBin = std::max(0, firstNonZeroBin - winHalfWidth - 1);
1609 lastNonZeroBin = std::min(hNbins - 1, lastNonZeroBin + winHalfWidth + 1);
1610
1611 // if null histogram quit
1612 if (firstNullPart)
1613 return maxPos;
1614
1615 // determine the size of the sliding window in the fit case
1616
1617 // sliding window
1618 const int nwidth = 2 * winHalfWidth + 1;
1619 double winsum = 0.;
1620
1621 for (int ibin = 0; ibin < nwidth; ++ibin) {
1622 winsum += theHist->GetBinContent(ibin);
1623 }
1624 double winmax = winsum;
1625
1626 int max_bin = 0.;
1627 int iBinL = firstNonZeroBin;
1628 int iBinR = iBinL + 2 * winHalfWidth;
1629 bool goingUp = true;
1630
1631 do {
1632 ++iBinL;
1633 ++iBinR;
1634 const double deltawin = theHist->GetBinContent(iBinR) - theHist->GetBinContent(iBinL - 1);
1635
1636 if (deltawin < 0) {
1637 if (goingUp) {
1638 // if were climbing and now loose more on the left
1639 // than win on the right. This was a local maxima
1640 if (winsum > winmax) {
1641 // global maximum one so far
1642 winmax = winsum;
1643 max_bin = (iBinR + iBinL) / 2 - 1;
1644 }
1645 goingUp = false; // now going down
1646 }
1647 } else {
1648 // do not care about minima, simply indicate we are going down
1649 goingUp = true;
1650 }
1651
1652 winsum += deltawin;
1653
1654 } while (iBinR < lastNonZeroBin);
1655
1656 // now compute average
1657 int iBinMin = max_bin - winHalfWidth;
1658 if (iBinMin < 0)
1659 iBinMin = 0;
1660 int iBinMax = max_bin + winHalfWidth;
1661 if (iBinMax >= hNbins)
1662 iBinMax = hNbins - 1;
1663 double sumw = 0;
1664 double sumx = 0;
1665 for (int iBin = iBinMin; iBin <= iBinMax; ++iBin) {
1666 const double weight = theHist->GetBinContent(iBin);
1667 sumw += weight;
1668 sumx += weight * theHist->GetBinCenter(iBin);
1669 }
1670
1671 double maxPosWin = -1.;
1672
1673 if (sumw > 0.) {
1674 maxPosWin = sumx / sumw;
1675 }
1676 // prob if the fraction of events in the window
1677 prob = sumw / totalSumw;
1678
1679 // Definitions of some useful parameters
1680
1681 const double h_rms = theHist->GetRMS(1);
1682 histInfo[HistInfo::RMS] = h_rms;
1683
1684 double num = 0;
1685 double numerator = 0;
1686 double denominator = 0;
1687 bool nullBin = false;
1688
1689 for (int i = iBinMin; i < iBinMax; ++i) {
1690 double binError = theHist->GetBinError(i);
1691 if (binError < 1e-10) {
1692 nullBin = true;
1693 }
1694 double binErrorSquare = std::pow(binError, 2);
1695 num = theHist->GetBinContent(i) / (binErrorSquare);
1696 numerator = numerator + num;
1697 denominator = denominator + (1 / (binErrorSquare));
1698 }
1699 if (numerator < 1e-10 || denominator < 1e-10 || nullBin == true) {
1700 histInfo[HistInfo::MEANBIN] = -1;
1701 } else {
1702 histInfo[HistInfo::MEANBIN] = sqrt(1 / denominator) / (numerator / denominator);
1703 }
1704
1705 // stop here if only looking for sliding window
1706 if (maxHistStrategy == MaxHistStrategy::SLIDINGWINDOW) {
1707 return maxPosWin;
1708 }
1709
1710 maxPos = maxPosWin;
1711 // now FIT maxHistStrategy==MaxHistStrategy::FIT
1712
1713 // now mass fit in range defined by sliding window
1714 // window will be around maxPos
1715 const double binWidth = theHist->GetBinCenter(2) - theHist->GetBinCenter(1);
1716 double fitWidth = (winHalfWidth + 0.5) * binWidth;
1717 // fit range 2 larger than original window range, 3 if less than 20% of the
1718 // histogram in slinding window
1719
1720 if (prob > 0.2) {
1721 fitWidth *= 2.;
1722 } else {
1723 fitWidth *= 3.;
1724 }
1725 // fit option : Q == Quiet, no printout S result of the fit returned in
1726 // TFitResultPtr N do not draw the resulting function
1727
1728 // if debug plot the fitted function
1729 TString fitOption = debug ? "QS" : "QNS";
1730 // root fit
1731 // Sets initial values
1732 m_fFitting->SetParameters(sumw / winHalfWidth, maxPos, 0.0025);
1733 // TFitResultPtr
1734 // fitRes=theHist->Fit("pol2",fitOption,"",maxPos-fitWidth,maxPos+fitWidth);
1735 TFitResultPtr fitRes =
1736 theHist->Fit(m_fFitting, fitOption, "", maxPos - fitWidth, maxPos + fitWidth);
1737
1738 double maxPosFit = -1.;
1739
1740 if (int(fitRes) == 0) {
1741 // root fit
1742 histInfo[HistInfo::CHI2] = fitRes->Chi2();
1743 const double mMax = fitRes->Parameter(0);
1744 const double mMean = fitRes->Parameter(1);
1745 const double mInvWidth2 = fitRes->Parameter(2);
1746 double mMaxError = fitRes->ParError(0);
1747 m_PrintmMaxError = mMaxError;
1748 double mMeanError = fitRes->ParError(1);
1749 m_PrintmMeanError = mMeanError;
1750 double mInvWidth2Error = fitRes->ParError(2);
1751 m_PrintmInvWidth2Error = mInvWidth2Error;
1752 mMeanError = 0.; // avoid warning
1753 mInvWidth2Error = 0.; // avoid warning
1754 const double c = mMax * (1 - 4 * mMean * mMean * mInvWidth2);
1755 const double b = 8 * mMax * mMean * mInvWidth2;
1756 const double a = -4 * mMax * mInvWidth2;
1757 // when built in polynomial fit
1758 // const double c=fitRes->Parameter(0);
1759 // const double b=fitRes->Parameter(1);
1760 // const double a=fitRes->Parameter(2);
1761
1762 const double h_discri = b * b - 4 * a * c;
1763 histInfo[HistInfo::DISCRI] = h_discri;
1764 const double sqrth_discri = sqrt(h_discri);
1765 const double h_fitLength = sqrth_discri / a;
1766 histInfo[HistInfo::FITLENGTH] = h_fitLength;
1767 histInfo[HistInfo::TANTHETA] = 2 * a / sqrth_discri;
1768 histInfo[HistInfo::TANTHETAW] = 2 * a * sumw / sqrth_discri;
1769 histInfo[HistInfo::RMSVSDISCRI] = h_rms / h_fitLength;
1770 // compute maximum position (only if inverted parabola)
1771 if (a < 0)
1772 maxPosFit = -b / (2 * a);
1773 }
1774
1775 // keep fit result only if within 80% of fit window, and fit succeeded
1776 if (maxPosFit >= 0. and std::abs(maxPosFit - maxPosWin) < 0.8 * fitWidth) {
1777 histInfo[HistInfo::PROB] = prob;
1778 return maxPosFit;
1779 } else {
1780 // otherwise keep the weighted average
1781 // negate prob just to flag such event
1782 prob = -prob;
1783 histInfo[HistInfo::PROB] = prob;
1784 return maxPosWin;
1785 }
1786}
static Double_t a
void binWidth(TH1 *h)
Definition listroot.cxx:80

◆ NuPsolutionLFV()

int MissingMassCalculator::NuPsolutionLFV ( const XYVector & met_vec,
const PtEtaPhiMVector & tau,
const double & m_nu,
std::vector< PtEtaPhiMVector > & nu_vec )
inlineprivate

Definition at line 755 of file MissingMassCalculator.cxx.

757 {
758 int solution_code = 0; // 0 with no solution, 1 with solution
759
760 nu_vec.clear();
761 PxPyPzMVector nu(met_vec.X(), met_vec.Y(), 0.0, l_nu);
762 PxPyPzMVector nu2(met_vec.X(), met_vec.Y(), 0.0, l_nu);
763
764 const double Mtau = ParticleConstants::tauMassInMeV / GEV;
765 // double msq = (Mtau*Mtau-tau.M()*tau.M())/2;
766 double msq = (Mtau * Mtau - tau.M() * tau.M() - l_nu * l_nu) /
767 2; // to take into account the fact that 2-nu systema has mass
768 double gamma = nu.Px() * nu.Px() + nu.Py() * nu.Py();
769 double beta = tau.Px() * nu.Px() + tau.Py() * nu.Py() + msq;
770 double a = tau.E() * tau.E() - tau.Pz() * tau.Pz();
771 double b = -2 * tau.Pz() * beta;
772 double c = tau.E() * tau.E() * gamma - beta * beta;
773 if ((b * b - 4 * a * c) < 0)
774 return solution_code; // no solution found
775 else
776 solution_code = 2;
777 double pvz1 = (-b + sqrt(b * b - 4 * a * c)) / (2 * a);
778 double pvz2 = (-b - sqrt(b * b - 4 * a * c)) / (2 * a);
779
780 nu.SetCoordinates(met_vec.X(), met_vec.Y(), pvz1, l_nu);
781 nu2.SetCoordinates(met_vec.X(), met_vec.Y(), pvz2, l_nu);
782
783 PtEtaPhiMVector return_nu(nu.Pt(), nu.Eta(), nu.Phi(), nu.M());
784 PtEtaPhiMVector return_nu2(nu2.Pt(), nu2.Eta(), nu2.Phi(), nu2.M());
785 nu_vec.push_back(return_nu);
786 nu_vec.push_back(return_nu2);
787 return solution_code;
788}

◆ NuPsolutionV3()

int MissingMassCalculator::NuPsolutionV3 ( const double & mNu1,
const double & mNu2,
const double & phi1,
const double & phi2,
int & nsol1,
int & nsol2 )
inlineprivate

Definition at line 550 of file MissingMassCalculator.cxx.

552 {
553
554 // Pv1, Pv2 : visible tau decay product momentum
555 // Pn1 Pn2 : neutrino momentum
556 // phi1, phi2 : neutrino azymutal angles
557 // PTmiss2=PTmissy Cos[phi2] - PTmissx Sin[phi2]
558 // PTmiss2cscdphi=PTmiss2/Sin[phi1-phi2]
559 // Pv1proj=Pv1x Cos[phi1] + Pv1y Sin[phi1]
560 // M2noma1=Mtau^2-Mv1^2-Mn1^2
561 // ETv1^2=Ev1^2-Pv1z^2
562
563 // discriminant : 16 Ev1^2 (M2noma1^2 + 4 M2noma1 PTmiss2cscdphi Pv1proj - 4
564 // (ETv1^2 (Mn1^2 + PTmiss2cscdphi^2) - PTmiss2cscdphi^2 Pv1proj^2))
565 // two solutions for epsilon = +/-1
566 // Pn1z=(1/(2 ETv1^2))(epsilon Ev1 Sqrt[ M2noma1^2 + 4 M2noma1 PTmiss2cscdphi
567 // Pv1proj - 4 (ETv1^2 (Mn1^2 + qPTmiss2cscdphi^2) - PTmiss2cscdphi^2
568 // Pv1proj^2)] + M2noma1 Pv1z + 2 PTmiss2cscdphi Pv1proj Pv1z)
569 // with conditions: M2noma1 + 2 PTmiss2cscdphi Pv1proj + 2 Pn1z Pv1z > 0
570 // PTn1 -> PTmiss2 Csc[phi1 - phi2]
571
572 // if initialisation precompute some quantities
573 int solution_code = 0; // 0 with no solution, 1 with solution
574 nsol1 = 0;
575 nsol2 = 0;
576
577 // Variables used to test PTn1 and PTn2 > 0
578
579 const double &pTmissx = preparedInput.m_MEtX;
580 const double &pTmissy = preparedInput.m_MEtY;
581
583 double pTmiss2 = pTmissy * m_cosPhi2 - pTmissx * m_sinPhi2;
584
585 int dPhiSign = 0;
586 dPhiSign = fixPhiRange(phi1 - phi2) > 0 ? +1 : -1;
587
588 // Test if PTn1 and PTn2 > 0. Then MET vector is between the two neutrino
589 // vector
590
591 if (pTmiss2 * dPhiSign < 0) {
592 ++m_testptn1;
593 return solution_code;
594 }
595
597 double pTmiss1 = pTmissy * m_cosPhi1 - pTmissx * m_sinPhi1;
598
599 if (pTmiss1 * (-dPhiSign) < 0) {
600 ++m_testptn2;
601 return solution_code;
602 }
603
604 // Variables used to calculate discri1
605
606 double m2Vis1 = m_tauVec1M * m_tauVec1M;
607 m_ET2v1 = std::pow(m_tauVec1E, 2) - std::pow(m_tauVec1Pz, 2);
608 m_m2Nu1 = mNu1 * mNu1;
609 double m2noma1 = m_mTau2 - m_m2Nu1 - m2Vis1;
610 double m4noma1 = m2noma1 * m2noma1;
611 double pv1proj = m_tauVec1Px * m_cosPhi1 + m_tauVec1Py * m_sinPhi1;
612 double p2v1proj = std::pow(pv1proj, 2);
613 double sinDPhi2 = m_cosPhi2 * m_sinPhi1 - m_sinPhi2 * m_cosPhi1; // sin(Phi1-Phi2)
614 double pTmiss2CscDPhi = pTmiss2 / sinDPhi2;
615 double &pTn1 = pTmiss2CscDPhi;
616 double pT2miss2CscDPhi = pTmiss2CscDPhi * pTmiss2CscDPhi;
617
618 // Test on discri1
619 const double discri1 = m4noma1 + 4 * m2noma1 * pTmiss2CscDPhi * pv1proj -
620 4 * (m_ET2v1 * (m_m2Nu1 + pT2miss2CscDPhi) - (pT2miss2CscDPhi * p2v1proj));
621
622 if (discri1 < 0) // discriminant negative -> no solution
623 {
625 return solution_code;
626 }
627
628 // Variables used to calculate discri2
629 double m2Vis2 = m_tauVec2M * m_tauVec2M;
630 m_ET2v2 = std::pow(m_tauVec2E, 2) - std::pow(m_tauVec2Pz, 2);
631 m_m2Nu2 = mNu2 * mNu2;
632 double m2noma2 = m_mTau2 - m_m2Nu2 - m2Vis2;
633 double m4noma2 = m2noma2 * m2noma2;
634 double pv2proj = m_tauVec2Px * m_cosPhi2 + m_tauVec2Py * m_sinPhi2;
635 double p2v2proj = std::pow(pv2proj, 2);
636 double sinDPhi1 = -sinDPhi2;
637 double pTmiss1CscDPhi = pTmiss1 / sinDPhi1;
638 double &pTn2 = pTmiss1CscDPhi;
639 double pT2miss1CscDPhi = pTmiss1CscDPhi * pTmiss1CscDPhi;
640
641 const double discri2 = m4noma2 + 4 * m2noma2 * pTmiss1CscDPhi * pv2proj -
642 4 * (m_ET2v2 * (m_m2Nu2 + pT2miss1CscDPhi) - (pT2miss1CscDPhi * p2v2proj));
643
644 if (discri2 < 0) // discriminant negative -> no solution
645 {
647 return solution_code;
648 }
649
650 // this should be done only once we know there are solutions for nu2
652 m_Ev1 = sqrt(m_E2v1);
653 double sqdiscri1 = sqrt(discri1);
654 double first1 =
655 (m2noma1 * m_tauVec1Pz + 2 * pTmiss2CscDPhi * pv1proj * m_tauVec1Pz) / (2 * m_ET2v1);
656 double second1 = sqdiscri1 * m_Ev1 / (2 * m_ET2v1);
657
658 // first solution
659 double pn1Z = first1 + second1;
660
661 if (m2noma1 + 2 * pTmiss2CscDPhi * pv1proj + 2 * pn1Z * m_tauVec1Pz >
662 0) // Condition for solution to exist
663 {
664 m_nuvecsol1[nsol1].SetPxPyPzE(pTn1 * m_cosPhi1, pTn1 * m_sinPhi1, pn1Z,
665 sqrt(std::pow(pTn1, 2) + std::pow(pn1Z, 2) + m_m2Nu1));
666
667 ++nsol1;
668 }
669
670 pn1Z = first1 - second1;
671
672 if (m2noma1 + 2 * pTmiss2CscDPhi * pv1proj + 2 * pn1Z * m_tauVec1Pz >
673 0) // Condition for solution to exist
674 {
675
676 m_nuvecsol1[nsol1].SetPxPyPzE(pTn1 * m_cosPhi1, pTn1 * m_sinPhi1, pn1Z,
677 sqrt(std::pow(pTn1, 2) + std::pow(pn1Z, 2) + m_m2Nu1));
678
679 ++nsol1;
680 }
681
682 if (nsol1 == 0) {
683 ++m_nosol1;
684 return solution_code;
685 }
686
688 m_Ev2 = sqrt(m_E2v2);
689 double sqdiscri2 = sqrt(discri2);
690 double first2 =
691 (m2noma2 * m_tauVec2Pz + 2 * pTmiss1CscDPhi * pv2proj * m_tauVec2Pz) / (2 * m_ET2v2);
692 double second2 = sqdiscri2 * m_Ev2 / (2 * m_ET2v2);
693
694 // second solution
695 double pn2Z = first2 + second2;
696
697 if (m2noma2 + 2 * pTmiss1CscDPhi * pv2proj + 2 * pn2Z * m_tauVec2Pz >
698 0) // Condition for solution to exist
699 {
700 m_nuvecsol2[nsol2].SetPxPyPzE(pTn2 * m_cosPhi2, pTn2 * m_sinPhi2, pn2Z,
701 sqrt(std::pow(pTn2, 2) + std::pow(pn2Z, 2) + m_m2Nu2));
702
703 ++nsol2;
704 }
705
706 pn2Z = first2 - second2;
707 ;
708
709 if (m2noma2 + 2 * pTmiss1CscDPhi * pv2proj + 2 * pn2Z * m_tauVec2Pz >
710 0) // Condition for solution to exist
711 {
712 m_nuvecsol2[nsol2].SetPxPyPzE(pTn2 * m_cosPhi2, pTn2 * m_sinPhi2, pn2Z,
713 sqrt(std::pow(pTn2, 2) + std::pow(pn2Z, 2) + m_m2Nu2));
714
715 ++nsol2;
716 }
717
718 if (nsol2 == 0) {
719 ++m_nosol2;
720 return solution_code;
721 }
722
723 // Verification if solution exist
724
725 solution_code = 1;
726 ++m_iterNuPV3;
727
728 // double check solutions from time to time
729 if (m_iterNuPV3 % 1000 == 1) {
730 double pnux = m_nuvecsol1[0].Px() + m_nuvecsol2[0].Px();
731 double pnuy = m_nuvecsol1[0].Py() + m_nuvecsol2[0].Py();
732 double mtau1plus = (m_nuvecsol1[0] + m_tauVec1).M();
733 double mtau1moins = (m_nuvecsol1[1] + m_tauVec1).M();
734 double mtau2plus = (m_nuvecsol2[0] + m_tauVec2).M();
735 double mtau2moins = (m_nuvecsol2[1] + m_tauVec2).M();
736 if (std::abs(pnux - pTmissx) > 0.001 || std::abs(pnuy - pTmissy) > 0.001) {
737 Info("DiTauMassTools", "%s", ("NuPsolutionV3 ERROR Pnux-Met.X or Pnuy-Met.Y > 0.001 : " +
738 std::to_string(pnux - pTmissx) + " and " +
739 std::to_string(pnuy - pTmissx) + " " + "Invalid solutions")
740 .c_str());
741 }
742 if (std::abs(mtau1plus - m_mTau) > 0.001 || std::abs(mtau1moins - m_mTau) > 0.001 ||
743 std::abs(mtau2plus - m_mTau) > 0.001 || std::abs(mtau2moins - m_mTau) > 0.001) {
744 Info("DiTauMassTools", "%s", ("NuPsolutionV3 ERROR tau mass not recovered : " +
745 std::to_string(mtau1plus) + " " + std::to_string(mtau1moins) + " " +
746 std::to_string(mtau2plus) + " " + std::to_string(mtau2moins))
747 .c_str());
748 }
749 }
750
751 return solution_code;
752}
void fastSinCos(const double &phi, double &sinPhi, double &cosPhi)

◆ operator=()

MissingMassCalculator & DiTauMassTools::MissingMassCalculator::operator= ( const MissingMassCalculator & )
delete

◆ precomputeCache()

bool MissingMassCalculator::precomputeCache ( )
inlineprotected

Definition at line 2618 of file MissingMassCalculator.cxx.

2618 {
2619
2620 // copy tau 4 vect. If tau E scanning, these vectors will be modified
2621 m_tauVec1 = preparedInput.m_vistau1;
2622 m_tauVec2 = preparedInput.m_vistau2;
2623
2624 const XYVector &metVec = preparedInput.m_MetVec;
2625
2626 bool same = true;
2641
2643 same = updateDouble(std::pow(m_mTau, 2), m_mTau2) && same;
2647
2648 PtEtaPhiMVector Met4vec;
2649 Met4vec.SetPxPyPzE(preparedInput.m_MetVec.X(), preparedInput.m_MetVec.Y(), 0.0,
2650 preparedInput.m_MetVec.R());
2651 same = updateDouble((m_tauVec1 + m_tauVec2 + Met4vec).M(), m_Meff) && same;
2652
2653 same = updateDouble(preparedInput.m_HtOffset, preparedInput.m_htOffset) && same;
2654 // note that if useHT met_vec is actually -HT
2655 same = updateDouble(metVec.X(), preparedInput.m_inputMEtX) && same;
2656 same = updateDouble(metVec.Y(), preparedInput.m_inputMEtY) && same;
2657 same = updateDouble(metVec.R(), preparedInput.m_inputMEtT) && same;
2658
2659 return same;
2660}

◆ PrintOtherInput()

void MissingMassCalculator::PrintOtherInput ( )
private

Definition at line 403 of file MissingMassCalculator.cxx.

403 {
404 if (preparedInput.m_fUseVerbose != 1)
405 return;
406
407 Info("DiTauMassTools",
408 ".........................Other input.....................................");
409 Info("DiTauMassTools", "%s",
410 ("Beam energy =" + std::to_string(preparedInput.m_beamEnergy) +
411 " sqrt(S) for collisions =" + std::to_string(2.0 * preparedInput.m_beamEnergy))
412 .c_str());
413 Info("DiTauMassTools", "%s",
414 ("CalibrationSet " + MMCCalibrationSet::name[m_mmcCalibrationSet])
415 .c_str());
416 Info("DiTauMassTools", "%s",
417 ("LFV mode " + std::to_string(preparedInput.m_LFVmode) + " seed=" + std::to_string(m_seed))
418 .c_str());
419 Info("DiTauMassTools", "%s", ("usetauProbability=" + std::to_string(Prob->GetUseTauProbability()) +
420 " useTailCleanup=" + std::to_string(preparedInput.m_fUseTailCleanup))
421 .c_str());
422
423 if (preparedInput.m_InputReorder != 0) {
424 Info("DiTauMassTools",
425 "tau1 and tau2 were internally swapped (visible on prepared input printout)");
426 } else {
427 Info("DiTauMassTools", "tau1 and tau2 were NOT internally swapped");
428 }
429
430 Info("DiTauMassTools", "%s",
431 (" MEtLMin=" + std::to_string(m_MEtLMin) + " MEtLMax=" + std::to_string(m_MEtLMax)).c_str());
432 Info("DiTauMassTools", "%s",
433 (" MEtPMin=" + std::to_string(m_MEtPMin) + " MEtPMax=" + std::to_string(m_MEtPMax)).c_str());
434 Info("DiTauMassTools", "%s",
435 (" Phi1Min=" + std::to_string(m_Phi1Min) + " Phi1Max=" + std::to_string(m_Phi1Max)).c_str());
436 Info("DiTauMassTools", "%s",
437 (" Phi2Min=" + std::to_string(m_Phi2Min) + " Phi2Max=" + std::to_string(m_Phi2Max)).c_str());
438 Info("DiTauMassTools", "%s",
439 (" Mnu1Min=" + std::to_string(m_Mnu1Min) + " Mnu1Max=" + std::to_string(m_Mnu1Max)).c_str());
440 Info("DiTauMassTools", "%s",
441 (" Mnu2Min=" + std::to_string(m_Mnu2Min) + " Mnu2Max=" + std::to_string(m_Mnu2Max)).c_str());
442}

◆ PrintResults()

void MissingMassCalculator::PrintResults ( )
private

Definition at line 445 of file MissingMassCalculator.cxx.

445 {
446
447 if (preparedInput.m_fUseVerbose != 1)
448 return;
449
450 const PtEtaPhiMVector *origVisTau1 = 0;
451 const PtEtaPhiMVector *origVisTau2 = 0;
452
453 if (preparedInput.m_InputReorder == 0) {
454 origVisTau1 = &preparedInput.m_vistau1;
455 origVisTau2 = &preparedInput.m_vistau2;
456 } else // input order was flipped
457 {
458 origVisTau1 = &preparedInput.m_vistau2;
459 origVisTau2 = &preparedInput.m_vistau1;
460 }
461
463
464 Info("DiTauMassTools",
465 "------------- Printing Final Results for MissingMassCalculator --------------");
466 Info("DiTauMassTools",
467 ".............................................................................");
468 Info("DiTauMassTools", "%s", ("Fit status=" + std::to_string(OutputInfo.m_FitStatus)).c_str());
469
470 for (int imeth = 0; imeth < MMCFitMethod::MAX; ++imeth) {
471 Info("DiTauMassTools", "%s",
472 ("___ Results for " + MMCFitMethod::name[imeth] + "Method ___")
473 .c_str());
474 Info("DiTauMassTools", "%s",
475 (" signif=" + std::to_string(OutputInfo.m_FitSignificance[imeth])).c_str());
476 Info("DiTauMassTools", "%s", (" mass=" + std::to_string(OutputInfo.m_FittedMass[imeth])).c_str());
477 Info("DiTauMassTools", "%s", (" rms/mpv=" + std::to_string(OutputInfo.m_RMS2MPV)).c_str());
478
479 if (imeth == MMCFitMethod::MLM) {
480 Info("DiTauMassTools", " no 4-momentum or MET from this method ");
481 continue;
482 }
483
484 if (OutputInfo.m_FitStatus <= 0) {
485 Info("DiTauMassTools", " fit failed ");
486 }
487
488 const PtEtaPhiMVector &tlvnu1 = OutputInfo.m_nuvec1[imeth];
489 const PtEtaPhiMVector &tlvnu2 = OutputInfo.m_nuvec2[imeth];
490 const PtEtaPhiMVector &tlvo1 = OutputInfo.m_objvec1[imeth];
491 const PtEtaPhiMVector &tlvo2 = OutputInfo.m_objvec2[imeth];
492 const XYVector &tvmet = OutputInfo.m_FittedMetVec[imeth];
493
494 Info("DiTauMassTools", "%s",
495 (" Neutrino-1: P=" + std::to_string(tlvnu1.P()) + " Pt=" + std::to_string(tlvnu1.Pt()) +
496 " Eta=" + std::to_string(tlvnu1.Eta()) + " Phi=" + std::to_string(tlvnu1.Phi()) +
497 " M=" + std::to_string(tlvnu1.M()) + " Px=" + std::to_string(tlvnu1.Px()) +
498 " Py=" + std::to_string(tlvnu1.Py()) + " Pz=" + std::to_string(tlvnu1.Pz()))
499 .c_str());
500 Info("DiTauMassTools", "%s",
501 (" Neutrino-2: P=" + std::to_string(tlvnu2.P()) + " Pt=" + std::to_string(tlvnu2.Pt()) +
502 " Eta=" + std::to_string(tlvnu2.Eta()) + " Phi=" + std::to_string(tlvnu2.Phi()) +
503 " M=" + std::to_string(tlvnu2.M()) + " Px=" + std::to_string(tlvnu2.Px()) +
504 " Py=" + std::to_string(tlvnu2.Py()) + " Pz=" + std::to_string(tlvnu2.Pz()))
505 .c_str());
506 Info("DiTauMassTools", "%s",
507 (" Tau-1: P=" + std::to_string(tlvo1.P()) + " Pt=" + std::to_string(tlvo1.Pt()) +
508 " Eta=" + std::to_string(tlvo1.Eta()) + " Phi=" + std::to_string(tlvo1.Phi()) +
509 " M=" + std::to_string(tlvo1.M()) + " Px=" + std::to_string(tlvo1.Px()) +
510 " Py=" + std::to_string(tlvo1.Py()) + " Pz=" + std::to_string(tlvo1.Pz()))
511 .c_str());
512 Info("DiTauMassTools", "%s",
513 (" Tau-2: P=" + std::to_string(tlvo2.P()) + " Pt=" + std::to_string(tlvo2.Pt()) +
514 " Eta=" + std::to_string(tlvo2.Eta()) + " Phi=" + std::to_string(tlvo2.Phi()) +
515 " M=" + std::to_string(tlvo2.M()) + " Px=" + std::to_string(tlvo2.Px()) +
516 " Py=" + std::to_string(tlvo2.Py()) + " Pz=" + std::to_string(tlvo2.Pz()))
517 .c_str());
518
519 Info("DiTauMassTools", "%s",
520 (" dR(nu1-visTau1)=" + std::to_string(DeltaR(tlvnu1,*origVisTau1))).c_str());
521 Info("DiTauMassTools", "%s",
522 (" dR(nu2-visTau2)=" + std::to_string(DeltaR(tlvnu2,*origVisTau2))).c_str());
523
524 Info("DiTauMassTools", "%s",
525 (" Fitted MET =" + std::to_string(tvmet.R()) + " Phi=" + std::to_string(tlvnu1.Phi()) +
526 " Px=" + std::to_string(tvmet.X()) + " Py=" + std::to_string(tvmet.Y()))
527 .c_str());
528
529 Info("DiTauMassTools", "%s", (" Resonance: P=" + std::to_string(OutputInfo.m_totalvec[imeth].P()) +
530 " Pt=" + std::to_string(OutputInfo.m_totalvec[imeth].Pt()) +
531 " Eta=" + std::to_string(OutputInfo.m_totalvec[imeth].Eta()) +
532 " Phi=" + std::to_string(OutputInfo.m_totalvec[imeth].Phi()) +
533 " M=" + std::to_string(OutputInfo.m_totalvec[imeth].M()) +
534 " Px=" + std::to_string(OutputInfo.m_totalvec[imeth].Px()) +
535 " Py=" + std::to_string(OutputInfo.m_totalvec[imeth].Py()) +
536 " Pz=" + std::to_string(OutputInfo.m_totalvec[imeth].Pz()))
537 .c_str());
538 }
539
540 return;
541}

◆ probCalculatorV9fast()

int MissingMassCalculator::probCalculatorV9fast ( const double & phi1,
const double & phi2,
const double & M_nu1,
const double & M_nu2 )
inlineprotected

Definition at line 1793 of file MissingMassCalculator.cxx.

1795 {
1796 // bool debug=true;
1797
1798 int nsol1;
1799 int nsol2;
1800
1801 const int solution = NuPsolutionV3(M_nu1, M_nu2, phi1, phi2, nsol1, nsol2);
1802
1803 if (solution != 1)
1804 return -4;
1805 // refineSolutions ( M_nu1,M_nu2,
1806 // met_smearL,met_smearP,metvec_tmp.R(),
1807 // nsol1, nsol2,m_Mvis,m_Meff);
1808 refineSolutions(M_nu1, M_nu2, nsol1, nsol2, m_Mvis, m_Meff);
1809
1810 if (m_nsol <= 0)
1811 return 0;
1812
1813 // success
1814
1815 return m_nsol; // for backward compatibility
1816}
int refineSolutions(const double &M_nu1, const double &M_nu2, const int nsol1, const int nsol2, const double &Mvis, const double &Meff)
int NuPsolutionV3(const double &mNu1, const double &mNu2, const double &phi1, const double &phi2, int &nsol1, int &nsol2)

◆ refineSolutions()

int MissingMassCalculator::refineSolutions ( const double & M_nu1,
const double & M_nu2,
const int nsol1,
const int nsol2,
const double & Mvis,
const double & Meff )
inlineprotected

Definition at line 1819 of file MissingMassCalculator.cxx.

1823{
1824 m_nsol = 0;
1825
1826 if (int(m_probFinalSolVec.size()) < m_nsolfinalmax)
1827 Error("DiTauMassTools", "%s",
1828 ("refineSolutions ERROR probFinalSolVec.size() should be " + std::to_string(m_nsolfinalmax))
1829 .c_str());
1830 if (int(m_mtautauFinalSolVec.size()) < m_nsolfinalmax)
1831 Error("DiTauMassTools", "%s",
1832 ("refineSolutions ERROR mtautauSolVec.size() should be " + std::to_string(m_nsolfinalmax))
1833 .c_str());
1834 if (int(m_nu1FinalSolVec.size()) < m_nsolfinalmax)
1835 Error("DiTauMassTools", "%s",
1836 ("refineSolutions ERROR nu1FinalSolVec.size() should be " + std::to_string(m_nsolfinalmax))
1837 .c_str());
1838 if (int(m_nu2FinalSolVec.size()) < m_nsolfinalmax)
1839 Error("DiTauMassTools", "%s",
1840 ("refineSolutions ERROR nu2FinalSolVec.size() should be " + std::to_string(m_nsolfinalmax))
1841 .c_str());
1842 if (nsol1 > int(m_nsolmax))
1843 Error("DiTauMassTools", "%s", ("refineSolutions ERROR nsol1 " + std::to_string(nsol1) +
1844 " > nsolmax !" + std::to_string(m_nsolmax))
1845 .c_str());
1846 if (nsol2 > int(m_nsolmax))
1847 Error("DiTauMassTools", "%s", ("refineSolutions ERROR nsol1 " + std::to_string(nsol2) +
1848 " > nsolmax !" + std::to_string(m_nsolmax))
1849 .c_str());
1850
1851 int ngoodsol1 = 0;
1852 int ngoodsol2 = 0;
1853 double constProb =
1854 Prob->apply(preparedInput, -99, -99, PtEtaPhiMVector(0, 0, 0, 0), PtEtaPhiMVector(0, 0, 0, 0),
1855 PtEtaPhiMVector(0, 0, 0, 0), PtEtaPhiMVector(0, 0, 0, 0), true, false, false);
1856
1857 for (int j1 = 0; j1 < nsol1; ++j1) {
1858 PtEtaPhiMVector &nuvec1_tmpj = m_nuvecsol1[j1];
1859 PtEtaPhiMVector &tauvecsol1j = m_tauvecsol1[j1];
1860 double &tauvecprob1j = m_tauvecprob1[j1];
1861 tauvecprob1j = 0.;
1862 // take first or second solution
1863 // no time to call rndm, switch more or less randomely, according to an
1864 // oscillating switch perturbed by m_phi1
1865 if (nsol1 > 1) {
1866 if (j1 == 0) { // decide at the first solution which one we will take
1867 const int pickInt = std::abs(10000 * m_Phi1);
1868 const int pickDigit = pickInt - 10 * (pickInt / 10);
1869 if (pickDigit < 5)
1871 }
1873 }
1874
1875 if (!m_switch1) {
1876 nuvec1_tmpj.SetCoordinates(nuvec1_tmpj.Pt(), nuvec1_tmpj.Eta(), nuvec1_tmpj.Phi(), M_nu1);
1877 tauvecsol1j.SetPxPyPzE(0., 0., 0., 0.);
1878 tauvecsol1j += nuvec1_tmpj;
1879 tauvecsol1j += m_tauVec1;
1880 if (tauvecsol1j.E() >= preparedInput.m_beamEnergy)
1881 continue;
1882 tauvecprob1j = Prob->apply(preparedInput, preparedInput.m_type_visTau1, -99, m_tauVec1,
1883 PtEtaPhiMVector(0, 0, 0, 0), nuvec1_tmpj,
1884 PtEtaPhiMVector(0, 0, 0, 0), false, true, false);
1885 ++ngoodsol1;
1886 }
1887
1888 for (int j2 = 0; j2 < nsol2; ++j2) {
1889 PtEtaPhiMVector &nuvec2_tmpj = m_nuvecsol2[j2];
1890 PtEtaPhiMVector &tauvecsol2j = m_tauvecsol2[j2];
1891 double &tauvecprob2j = m_tauvecprob2[j2];
1892 if (j1 == 0) {
1893 tauvecprob2j = 0.;
1894 // take first or second solution
1895 // no time to call rndm, switch more or less randomely, according to an
1896 // oscillating switch perturbed by m_phi2
1897 if (nsol2 > 1) {
1898 if (j2 == 0) { // decide at the first solution which one we will take
1899 const int pickInt = std::abs(10000 * m_Phi2);
1900 const int pickDigit = pickInt - 10 * int(pickInt / 10);
1901 if (pickDigit < 5)
1903 }
1905 }
1906
1907 if (!m_switch2) {
1908 nuvec2_tmpj.SetCoordinates(nuvec2_tmpj.Pt(), nuvec2_tmpj.Eta(), nuvec2_tmpj.Phi(), M_nu2);
1909 tauvecsol2j.SetPxPyPzE(0., 0., 0., 0.);
1910 tauvecsol2j += nuvec2_tmpj;
1911 tauvecsol2j += m_tauVec2;
1912 if (tauvecsol2j.E() >= preparedInput.m_beamEnergy)
1913 continue;
1914 tauvecprob2j = Prob->apply(preparedInput, -99, preparedInput.m_type_visTau2,
1915 PtEtaPhiMVector(0, 0, 0, 0), m_tauVec2,
1916 PtEtaPhiMVector(0, 0, 0, 0), nuvec2_tmpj, false, true, false);
1917 ++ngoodsol2;
1918 }
1919 }
1920 if (tauvecprob1j == 0.)
1921 continue;
1922 if (tauvecprob2j == 0.)
1923 continue;
1924
1925 double totalProb = 1.;
1926
1927 m_tautau_tmp.SetPxPyPzE(0., 0., 0., 0.);
1928 m_tautau_tmp += tauvecsol1j;
1929 m_tautau_tmp += tauvecsol2j;
1930 const double mtautau = m_tautau_tmp.M();
1931
1932 if (TailCleanUp(m_tauVec1, nuvec1_tmpj, m_tauVec2, nuvec2_tmpj, mtautau, Mvis, Meff,
1933 preparedInput.m_DelPhiTT) == 0) {
1934 continue;
1935 }
1936
1937 totalProb *=
1938 (constProb * tauvecprob1j * tauvecprob2j *
1939 Prob->apply(preparedInput, preparedInput.m_type_visTau1, preparedInput.m_type_visTau2,
1940 m_tauVec1, m_tauVec2, nuvec1_tmpj, nuvec2_tmpj, false, false, true));
1941
1942 if (totalProb <= 0) {
1943 if (preparedInput.m_fUseVerbose)
1944 Warning("DiTauMassTools", "%s",
1945 ("null proba solution, rejected "+std::to_string(totalProb)).c_str());
1946 } else {
1947 // only count solution with non zero probability
1948 m_totalProbSum += totalProb;
1949 m_mtautauSum += mtautau;
1950
1951 if (m_nsol >= int(m_nsolfinalmax)) {
1952 Error("DiTauMassTools", "%s",
1953 ("refineSolutions ERROR nsol getting larger than nsolfinalmax!!! " +
1954 std::to_string(m_nsol))
1955 .c_str());
1956 Error("DiTauMassTools", "%s",
1957 (" j1 " + std::to_string(j1) + " j2 " + std::to_string(j2) + " nsol1 " +
1958 std::to_string(nsol1) + " nsol2 " + std::to_string(nsol2))
1959 .c_str());
1960 --m_nsol; // overwrite last solution. However this should really never
1961 // happen
1962 }
1963
1964 // good solution found, copy in vector
1965 m_mtautauFinalSolVec[m_nsol] = mtautau;
1966 m_probFinalSolVec[m_nsol] = totalProb;
1967
1968 PtEtaPhiMVector &nu1Final = m_nu1FinalSolVec[m_nsol];
1969 PtEtaPhiMVector &nu2Final = m_nu2FinalSolVec[m_nsol];
1970 // for (int iv=0;iv<4;++iv){
1971
1972 nu1Final.SetPxPyPzE(nuvec1_tmpj.Px(), nuvec1_tmpj.Py(), nuvec1_tmpj.Pz(), nuvec1_tmpj.E());
1973 nu2Final.SetPxPyPzE(nuvec2_tmpj.Px(), nuvec2_tmpj.Py(), nuvec2_tmpj.Pz(), nuvec2_tmpj.E());
1974 // }
1975
1976 ++m_nsol;
1977 } // else totalProb<=0
1978
1979 } // loop j2
1980 } // loop j1
1981 if (ngoodsol1 == 0) {
1982 return -1;
1983 }
1984 if (ngoodsol2 == 0) {
1985 return -2;
1986 }
1987 return m_nsol;
1988}
int TailCleanUp(const PtEtaPhiMVector &vis1, const PtEtaPhiMVector &nu1, const PtEtaPhiMVector &vis2, const PtEtaPhiMVector &nu2, const double &mmc_mass, const double &vis_mass, const double &eff_mass, const double &dphiTT)

◆ RunMissingMassCalculator()

int MissingMassCalculator::RunMissingMassCalculator ( const xAOD::IParticle * part1,
const xAOD::IParticle * part2,
const xAOD::MissingET * met,
const int & njets )

Definition at line 197 of file MissingMassCalculator.cxx.

200 {
201 m_reRunWithBestMET = false;
202
203 OutputInfo.ClearOutput(preparedInput.m_fUseVerbose);
204 if (preparedInput.m_fUseVerbose == 1) {
205 Info("DiTauMassTools", "------------- Raw Input for MissingMassCalculator --------------");
206 }
207 FinalizeSettings(part1, part2, met, njets); // rawInput, preparedInput );
208 Prob->MET(preparedInput);
209 if (preparedInput.m_fUseVerbose == 1) {
210 Info("DiTauMassTools", "------------- Prepared Input for MissingMassCalculator--------------");
211 preparedInput.PrintInputInfo();
212 }
213
214 if (preparedInput.m_LFVmode < 0) {
215 // remove argument DiTauMassCalculatorV9Walk work directly on preparedInput
217
218 // re-running MMC for on failed events
219 if (m_fUseEfficiencyRecovery == 1 && OutputInfo.m_FitStatus != 1) {
220 // most events where MMC failed happened to have dPhi>2.9. Run re-fit only
221 // on these events
222 if (preparedInput.m_DelPhiTT > 2.9) {
223 // preparedInput.MetVec.Set(-(preparedInput.vistau1+preparedInput.vistau2).Px(),-(preparedInput.vistau1+preparedInput.vistau2).Py());
224 // // replace MET by MPT
225
226 XYVector dummy_met(-(preparedInput.m_vistau1 + preparedInput.m_vistau2).Px(),
227 -(preparedInput.m_vistau1 + preparedInput.m_vistau2).Py());
228 preparedInput.m_METcovphi = dummy_met.Phi();
229 double dummy_METres =
230 sqrt(pow(preparedInput.m_METsigmaL, 2) + pow(preparedInput.m_METsigmaP, 2));
231 preparedInput.m_METsigmaL =
232 dummy_METres * std::abs(cos(dummy_met.Phi() - preparedInput.m_MetVec.Phi()));
233 preparedInput.m_METsigmaP =
234 dummy_METres * std::abs(sin(dummy_met.Phi() - preparedInput.m_MetVec.Phi()));
235 if (preparedInput.m_METsigmaP < 5.0)
236 preparedInput.m_METsigmaP = 5.0;
237 m_nsigma_METscan_lh = 6.0; // increase range of MET scan
238 m_nsigma_METscan_hh = 6.0; // increase range of MET scan
239
240 OutputInfo.ClearOutput(preparedInput.m_fUseVerbose); // clear output stuff before re-running
241 OutputInfo.m_FitStatus = DitauMassCalculatorV9walk(); // run MMC again
242 }
243 }
244
245 }
246
247 // running MMC in LFV mode for reconstructing mass of X->lep+tau
248 else {
249 if (preparedInput.m_fUseVerbose == 1) {
250 Info("DiTauMassTools", "Calling DitauMassCalculatorV9lfv");
251 }
252 OutputInfo.m_FitStatus = DitauMassCalculatorV9lfv(false);
253 }
254
255 if(m_SaveLlhHisto){
256 TFile *outFile = TFile::Open("MMC_likelihoods.root", "UPDATE");
257 outFile->cd();
258 auto path = std::to_string(m_eventNumber);
259 if (!outFile->GetDirectory(path.c_str()))
260 outFile->mkdir(path.c_str());
261 outFile->cd(path.c_str());
262 m_fMfit_all->Write(m_fMfit_all->GetName(), TObject::kOverwrite);
263 m_fMEtP_all->Write(m_fMEtP_all->GetName(), TObject::kOverwrite);
264 m_fMEtL_all->Write(m_fMEtL_all->GetName(), TObject::kOverwrite);
265 m_fMnu1_all->Write(m_fMnu1_all->GetName(), TObject::kOverwrite);
266 m_fMnu2_all->Write(m_fMnu2_all->GetName(), TObject::kOverwrite);
267 m_fPhi1_all->Write(m_fPhi1_all->GetName(), TObject::kOverwrite);
268 m_fPhi2_all->Write(m_fPhi2_all->GetName(), TObject::kOverwrite);
269 m_fMfit_allNoWeight->Write(m_fMfit_allNoWeight->GetName(), TObject::kOverwrite);
270 m_fMfit_allGraph->Write("Graph", TObject::kOverwrite);
271 TH1D *nosol = new TH1D("nosol", "nosol", 7, 0, 7);
272 nosol->SetBinContent(1, m_testptn1);
273 nosol->SetBinContent(2, m_testptn2);
274 nosol->SetBinContent(3, m_testdiscri1);
275 nosol->SetBinContent(4, m_testdiscri2);
276 nosol->SetBinContent(5, m_nosol1);
277 nosol->SetBinContent(6, m_nosol1);
278 nosol->SetBinContent(7, m_iterNuPV3);
279 nosol->Write(nosol->GetName(), TObject::kOverwrite);
280 outFile->Write();
281 outFile->Close();
282 }
283
284 DoOutputInfo();
285 PrintResults();
286 preparedInput.ClearInput();
287 return 1;
288}
void FinalizeSettings(const xAOD::IParticle *part1, const xAOD::IParticle *part2, const xAOD::MissingET *met, const int &njets)
outFile
Comment Out Those You do not wish to run.
path
python interpreter configuration --------------------------------------—
Definition athena.py:126

◆ SaveLlhHisto()

void MissingMassCalculator::SaveLlhHisto ( const bool val)

Definition at line 3093 of file MissingMassCalculator.cxx.

3093 {
3095 if(!m_SaveLlhHisto) return;
3096
3097 float hEmax = 3000.0; // maximum energy (GeV)
3098 int hNbins = 1500;
3099 m_fMEtP_all = std::make_shared<TH1F>("MEtP_h1", "M", hNbins, -100.0,
3100 100.); // all solutions
3101 m_fMEtL_all = std::make_shared<TH1F>("MEtL_h1", "M", hNbins, -100.0,
3102 100.); // all solutions
3103 m_fMnu1_all = std::make_shared<TH1F>("Mnu1_h1", "M", hNbins, 0.0,
3104 hEmax); // all solutions
3105 m_fMnu2_all = std::make_shared<TH1F>("Mnu2_h1", "M", hNbins, 0.0,
3106 hEmax); // all solutions
3107 m_fPhi1_all = std::make_shared<TH1F>("Phi1_h1", "M", hNbins, -10.0,
3108 10.); // all solutions
3109 m_fPhi2_all = std::make_shared<TH1F>("Phi2_h1", "M", hNbins, -10.0,
3110 10.); // all solutions
3111 m_fMfit_allGraph = std::make_shared<TGraph>(); // all solutions
3112
3113 m_fMEtP_all->Sumw2();
3114 m_fMEtL_all->Sumw2();
3115 m_fMnu1_all->Sumw2();
3116 m_fMnu2_all->Sumw2();
3117 m_fPhi1_all->Sumw2();
3118 m_fPhi2_all->Sumw2();
3119
3120 m_fMEtP_all->SetDirectory(0);
3121 m_fMEtL_all->SetDirectory(0);
3122 m_fMnu1_all->SetDirectory(0);
3123 m_fMnu2_all->SetDirectory(0);
3124 m_fPhi1_all->SetDirectory(0);
3125 m_fPhi2_all->SetDirectory(0);
3126}

◆ SetBeamEnergy()

void DiTauMassTools::MissingMassCalculator::SetBeamEnergy ( const double val)
inline

Definition at line 399 of file MissingMassCalculator.h.

399{ m_beamEnergy=val; }

◆ SetdTheta3d_binMax()

void DiTauMassTools::MissingMassCalculator::SetdTheta3d_binMax ( const double val)
inline

Definition at line 356 of file MissingMassCalculator.h.

356{ m_dTheta3d_binMax=val; } // maximum step size for dTheta3D

◆ SetdTheta3d_binMin()

void DiTauMassTools::MissingMassCalculator::SetdTheta3d_binMin ( const double val)
inline

Definition at line 357 of file MissingMassCalculator.h.

357{ m_dTheta3d_binMin=val; } // minimal step size for dTheta3D

◆ SetEventNumber()

void DiTauMassTools::MissingMassCalculator::SetEventNumber ( const int eventNumber)
inline

Definition at line 358 of file MissingMassCalculator.h.

◆ SetFloatStoppingCheckFreq()

void DiTauMassTools::MissingMassCalculator::SetFloatStoppingCheckFreq ( const int val)
inline

Definition at line 397 of file MissingMassCalculator.h.

◆ SetFloatStoppingComp()

void DiTauMassTools::MissingMassCalculator::SetFloatStoppingComp ( const double val)
inline

Definition at line 398 of file MissingMassCalculator.h.

◆ SetFloatStoppingMinIter()

void DiTauMassTools::MissingMassCalculator::SetFloatStoppingMinIter ( const int val)
inline

Definition at line 396 of file MissingMassCalculator.h.

◆ SetLFVLeplepRefit()

void DiTauMassTools::MissingMassCalculator::SetLFVLeplepRefit ( const bool val)
inline

Definition at line 400 of file MissingMassCalculator.h.

◆ SetMeanbinStop()

void DiTauMassTools::MissingMassCalculator::SetMeanbinStop ( const double val)
inline

Definition at line 354 of file MissingMassCalculator.h.

◆ SetMnuScanRange()

void DiTauMassTools::MissingMassCalculator::SetMnuScanRange ( const double val)
inline

Definition at line 360 of file MissingMassCalculator.h.

◆ SetNiterFit1()

void DiTauMassTools::MissingMassCalculator::SetNiterFit1 ( const int val)
inline

Definition at line 348 of file MissingMassCalculator.h.

348{ m_niter_fit1=val; } // number of iterations per loop in dPhi loop

◆ SetNiterFit2()

void DiTauMassTools::MissingMassCalculator::SetNiterFit2 ( const int val)
inline

Definition at line 349 of file MissingMassCalculator.h.

349{ m_niter_fit2=val; } // number of iterations per loop in MET loop

◆ SetNiterFit3()

void DiTauMassTools::MissingMassCalculator::SetNiterFit3 ( const int val)
inline

Definition at line 350 of file MissingMassCalculator.h.

350{ m_niter_fit3=val; } // number of iterations per loop in Mnu loop

◆ SetNiterRandom()

void DiTauMassTools::MissingMassCalculator::SetNiterRandom ( const int val)
inline

Definition at line 351 of file MissingMassCalculator.h.

351{ m_NiterRandom=val; } // number of random iterations

◆ SetNsigmaMETscan()

void DiTauMassTools::MissingMassCalculator::SetNsigmaMETscan ( const double val)
inline

Definition at line 393 of file MissingMassCalculator.h.

393{ m_nsigma_METscan=val; } // number of sigma's for MET-scan

◆ SetNsigmaMETscan_hh()

void DiTauMassTools::MissingMassCalculator::SetNsigmaMETscan_hh ( const double val)
inline

Definition at line 392 of file MissingMassCalculator.h.

392{ m_nsigma_METscan_hh=val; } // number of sigma's for MET-scan in hh events

◆ SetNsigmaMETscan_lh()

void DiTauMassTools::MissingMassCalculator::SetNsigmaMETscan_lh ( const double val)
inline

Definition at line 391 of file MissingMassCalculator.h.

391{ m_nsigma_METscan_lh=val; } // number of sigma's for MET-scan in lh events

◆ SetNsigmaMETscan_ll()

void DiTauMassTools::MissingMassCalculator::SetNsigmaMETscan_ll ( const double val)
inline

Definition at line 390 of file MissingMassCalculator.h.

390{ m_nsigma_METscan_ll=val; } // number of sigma's for MET-scan in ll events

◆ SetNsucStop()

void DiTauMassTools::MissingMassCalculator::SetNsucStop ( const int val)
inline

Definition at line 352 of file MissingMassCalculator.h.

352{ m_NsucStop=val; } // Arrest criteria for Nsuccesses

◆ SetProposalTryEtau()

void DiTauMassTools::MissingMassCalculator::SetProposalTryEtau ( const double val)
inline

Definition at line 365 of file MissingMassCalculator.h.

◆ SetProposalTryMEt()

void DiTauMassTools::MissingMassCalculator::SetProposalTryMEt ( const double val)
inline

Definition at line 362 of file MissingMassCalculator.h.

◆ SetProposalTryMnu()

void DiTauMassTools::MissingMassCalculator::SetProposalTryMnu ( const double val)
inline

Definition at line 364 of file MissingMassCalculator.h.

◆ SetProposalTryPhi()

void DiTauMassTools::MissingMassCalculator::SetProposalTryPhi ( const double val)
inline

Definition at line 363 of file MissingMassCalculator.h.

◆ SetRMSStop()

void DiTauMassTools::MissingMassCalculator::SetRMSStop ( const int val)
inline

Definition at line 353 of file MissingMassCalculator.h.

353{ m_RMSStop=val;}

◆ SetRndmSeedAltering()

void DiTauMassTools::MissingMassCalculator::SetRndmSeedAltering ( const int val)
inline

Definition at line 355 of file MissingMassCalculator.h.

355{ m_RndmSeedAltering=val; } // number of iterations per loop in Mnu loop

◆ SetUseEfficiencyRecovery()

void DiTauMassTools::MissingMassCalculator::SetUseEfficiencyRecovery ( const bool val)
inline

Definition at line 367 of file MissingMassCalculator.h.

◆ SetUseFloatStopping()

void MissingMassCalculator::SetUseFloatStopping ( const bool val)

Definition at line 3128 of file MissingMassCalculator.cxx.

3128 {
3130 if(!m_fUseFloatStopping) return;
3131
3132 float hEmax = 3000.0; // maximum energy (GeV)
3133 int hNbins = 1500;
3134 m_fMmass_split1 = std::make_shared<TH1F>("mass_h1_1", "M", hNbins, 0.0, hEmax);
3135 m_fMEtP_split1 = std::make_shared<TH1F>("MEtP_h1_1", "M", hNbins, -100.0, 100.0);
3136 m_fMEtL_split1 = std::make_shared<TH1F>("MEtL_h1_1", "M", hNbins, -100.0, 100.0);
3137 m_fMnu1_split1 = std::make_shared<TH1F>("Mnu1_h1_1", "M", hNbins, 0.0, hEmax);
3138 m_fMnu2_split1 = std::make_shared<TH1F>("Mnu2_h1_1", "M", hNbins, 0.0, hEmax);
3139 m_fPhi1_split1 = std::make_shared<TH1F>("Phi1_h1_1", "M", hNbins, -10.0, 10.0);
3140 m_fPhi2_split1 = std::make_shared<TH1F>("Phi2_h1_1", "M", hNbins, -10.0, 10.0);
3141 m_fMmass_split2 = std::make_shared<TH1F>("mass_h1_2", "M", hNbins, 0.0, hEmax);
3142 m_fMEtP_split2 = std::make_shared<TH1F>("MEtP_h1_2", "M", hNbins, -100.0, 100.0);
3143 m_fMEtL_split2 = std::make_shared<TH1F>("MEtL_h1_2", "M", hNbins, -100.0, 100.0);
3144 m_fMnu1_split2 = std::make_shared<TH1F>("Mnu1_h1_2", "M", hNbins, 0.0, hEmax);
3145 m_fMnu2_split2 = std::make_shared<TH1F>("Mnu2_h1_2", "M", hNbins, 0.0, hEmax);
3146 m_fPhi1_split2 = std::make_shared<TH1F>("Phi1_h1_2", "M", hNbins, -10.0, 10.0);
3147 m_fPhi2_split2 = std::make_shared<TH1F>("Phi2_h1_2", "M", hNbins, -10.0, 10.0);
3148
3149 m_fMmass_split1->Sumw2();
3150 m_fMEtP_split1->Sumw2();
3151 m_fMEtL_split1->Sumw2();
3152 m_fMnu1_split1->Sumw2();
3153 m_fMnu2_split1->Sumw2();
3154 m_fPhi1_split1->Sumw2();
3155 m_fPhi2_split1->Sumw2();
3156 m_fMmass_split2->Sumw2();
3157 m_fMEtP_split2->Sumw2();
3158 m_fMEtL_split2->Sumw2();
3159 m_fMnu1_split2->Sumw2();
3160 m_fMnu2_split2->Sumw2();
3161 m_fPhi1_split2->Sumw2();
3162 m_fPhi2_split2->Sumw2();
3163
3164 m_fMmass_split1->SetDirectory(0);
3165 m_fMEtP_split1->SetDirectory(0);
3166 m_fMEtL_split1->SetDirectory(0);
3167 m_fMnu1_split1->SetDirectory(0);
3168 m_fMnu2_split1->SetDirectory(0);
3169 m_fPhi1_split1->SetDirectory(0);
3170 m_fPhi2_split1->SetDirectory(0);
3171 m_fMmass_split2->SetDirectory(0);
3172 m_fMEtP_split2->SetDirectory(0);
3173 m_fMEtL_split2->SetDirectory(0);
3174 m_fMnu1_split2->SetDirectory(0);
3175 m_fMnu2_split2->SetDirectory(0);
3176 m_fPhi1_split2->SetDirectory(0);
3177 m_fPhi2_split2->SetDirectory(0);
3178}

◆ SpaceWalkerInit()

void MissingMassCalculator::SpaceWalkerInit ( )
inlineprotected

Definition at line 2314 of file MissingMassCalculator.cxx.

2314 {
2315 // FIXME could use function pointer to switch between functions
2316 m_nsolOld = 0;
2317
2318 double METresX = preparedInput.m_METsigmaL; // MET resolution in direction parallel to MET
2319 // resolution major axis, for MET scan
2320 double METresY = preparedInput.m_METsigmaP; // MET resolution in direction perpendicular to
2321 // to MET resolution major axis, for MET scan
2322
2323 // precompute some quantities and store in m_ data members
2326 if (Prob->GetUseMnuProbability() == true && (preparedInput.m_tauTypes == TauTypes::ll || preparedInput.m_tauTypes == TauTypes::lh) ) Prob->setParamNuMass();
2327 Prob->setParamAngle(m_tauVec1, 1, preparedInput.m_type_visTau1);
2328 Prob->setParamAngle(m_tauVec2, 2, preparedInput.m_type_visTau2);
2329 Prob->setParamRatio(1, preparedInput.m_type_visTau1);
2330 Prob->setParamRatio(2, preparedInput.m_type_visTau2);
2331 }
2332
2333 // if m_nsigma_METscan was not set by user, set to default values
2334 if(m_nsigma_METscan == -1){
2335 if (preparedInput.m_tauTypes == TauTypes::ll) // both tau's are leptonic
2336 {
2338 } else if (preparedInput.m_tauTypes == TauTypes::lh) // lep had
2339 {
2341 } else // hh
2342 {
2344 }
2345 }
2346
2347 m_nsigma_METscan2 = std::pow(m_nsigma_METscan, 2);
2348
2349 const double deltaPhi1 = MaxDelPhi(preparedInput.m_type_visTau1, m_tauVec1P, m_dRmax_tau);
2350 const double deltaPhi2 = MaxDelPhi(preparedInput.m_type_visTau2, m_tauVec2P, m_dRmax_tau);
2351
2352 m_walkWeight = 1.;
2353
2354 // dummy initial value to avoid printout with random values
2355 m_Phi10 = 0.;
2356 m_Phi20 = 0.;
2357 m_MEtL0 = 0.;
2358 m_MEtP0 = 0.;
2359 m_Mnu10 = 0.;
2360 m_Mnu20 = 0.;
2361
2363
2364 // seeds the random generator in a reproducible way from the phi of both tau;
2365 double aux = std::abs(m_tauVec1Phi + double(m_tauVec2Phi) / 100. / TMath::Pi()) * 100;
2366 m_seed = (aux - floor(aux)) * 1E6 * (1 + m_RndmSeedAltering) + 13;
2367
2368 m_randomGen.SetSeed(m_seed);
2369 // int Niter=Niter_fit1; // number of points for each dR loop
2370 // int NiterMET=Niter_fit2; // number of iterations for each MET scan loop
2371 // int NiterMnu=Niter_fit3; // number of iterations for Mnu loop
2372
2373 // approximately compute the number of points from the grid scanning
2374 // divide by abritry number to recover timing with still better results
2375 // m_NiterRandom=(NiterMET+1)*(NiterMET+1)*4*Niter*Niter/10;
2376
2380
2384
2385 m_Mnu1Min = 0.;
2386 m_scanMnu1 = false;
2387 m_Mnu1 = m_Mnu1Min;
2388
2389 // for markov chain use factor 2
2391
2392 // NiterRandom set by user (default is -1). If negative, defines the default
2393 // here. no more automatic scaling for ll hl hh
2394 if (m_NiterRandom <= 0) {
2395 m_niterRandomLocal = 100000; // number of iterations for Markov for lh
2396 if (preparedInput.m_tauTypes == TauTypes::ll)
2397 m_niterRandomLocal *= 2; // multiplied for ll , unchecked
2398 if (preparedInput.m_tauTypes == TauTypes::hh)
2399 m_niterRandomLocal *= 5; // divided for hh ,checked
2400 } else {
2402 }
2403
2404 if (preparedInput.m_type_visTau1 == 8) {
2405 // m_Mnu1Max=m_mTau-m_tauVec1M;
2408 m_scanMnu1 = true;
2409 }
2410
2411 m_Mnu2Min = 0.;
2412 m_scanMnu2 = false;
2413 m_Mnu2 = m_Mnu2Min;
2414 if (preparedInput.m_type_visTau2 == 8) {
2415 // m_Mnu2Max=m_mTau-m_tauVec2M;
2418 m_scanMnu2 = true;
2419 }
2420
2421 m_MEtLMin = -m_nsigma_METscan * METresX;
2422 m_MEtLMax = +m_nsigma_METscan * METresX;
2424
2425 m_MEtPMin = -m_nsigma_METscan * METresY;
2426 m_MEtPMax = +m_nsigma_METscan * METresY;
2428
2429 m_eTau1Min = -1;
2430 m_eTau1Max = -1;
2431 m_eTau2Min = -1;
2432 m_eTau2Max = -1;
2433
2434 m_switch1 = true;
2435 m_switch2 = true;
2436
2439
2440 m_iter0 = -1;
2441 m_iterNuPV3 = 0;
2442 m_testptn1 = 0;
2443 m_testptn2 = 0;
2444 m_testdiscri1 = 0;
2445 m_testdiscri2 = 0;
2446 m_nosol1 = 0;
2447 m_nosol2 = 0;
2448 m_iterNsuc = 0;
2449 if (m_meanbinStop > 0) {
2451 } else {
2452 m_meanbinToBeEvaluated = false;
2453 }
2454
2458 m_markovNAccept = 0;
2460 // set full parameter space scannning for the first steps, until a solution is
2461 // found
2462 m_fullParamSpaceScan = true;
2463 // size of step. Needs to be tune. Start with simple heuristic.
2464 if (m_proposalTryMEt < 0) {
2465 m_MEtProposal = m_MEtPRange / 30.;
2466 } else {
2468 }
2469 if (m_ProposalTryPhi < 0) {
2470 m_PhiProposal = 0.04;
2471 } else {
2473 }
2474 // FIXME if m_Mnu1Range !ne m_Mnu2Range same proposal will be done
2475 if (m_scanMnu1) {
2476 if (m_ProposalTryMnu < 0) {
2477 m_MnuProposal = m_Mnu1Range / 10.;
2478 } else {
2480 }
2481 }
2482 if (m_scanMnu2) {
2483 if (m_ProposalTryMnu < 0) {
2484 m_MnuProposal = m_Mnu2Range / 10.;
2485 } else {
2487 }
2488 }
2489}
double MaxDelPhi(int tau_type, double Pvis, double dRmax_tau)
@ deltaPhi1
difference between the cluster eta (1st sampling) and the eta of the track extrapolated to the 1st sa...
@ deltaPhi2
difference between the cluster phi (second sampling) and the phi of the track extrapolated to the sec...

◆ SpaceWalkerWalk()

bool MissingMassCalculator::SpaceWalkerWalk ( )
inlineprotected

Definition at line 2494 of file MissingMassCalculator.cxx.

2494 {
2495 preparedInput.m_MEtX = -999.;
2496 preparedInput.m_MEtY = -999.;
2497
2498 ++m_iter0;
2499
2500 if (m_meanbinToBeEvaluated && m_iterNsuc == 500) {
2501 Info("DiTauMassTools", " in m_meanbinToBeEvaluated && m_iterNsuc==500 ");
2502 // for markov chain m_iterNsuc is the number of *accepted* points, so there
2503 // can be several iterations without any increment of m_iterNsuc. Hence need
2504 // to make sure meanbin is evaluated only once
2505 m_meanbinToBeEvaluated = false;
2506
2507 // Meanbin stopping criterion
2508 std::vector<double> histInfo(HistInfo::MAXHISTINFO);
2509 // SLIDINGWINDOW strategy to avoid doing the parabola fit now given it will
2510 // not be use
2512 double meanbin = histInfo.at(HistInfo::MEANBIN);
2513 if (meanbin < 0) {
2514 m_nsucStop = -1; // no meaningful meanbin switch back to niter criterion
2515 } else {
2516 double stopdouble = 500 * std::pow((meanbin / m_meanbinStop), 2);
2517 int stopint = stopdouble;
2518 m_nsucStop = stopint;
2519 }
2520 if (m_nsucStop < 500)
2521 return false;
2522 }
2523 // should be outside m_meanbinStop test
2524 if (m_iterNsuc == m_nsucStop)
2525 return false; // Critere d'arret pour nombre de succes
2526
2528 return false; // for now simple stopping criterion on number of iteration
2529
2530 // floating stopping criterion, reduces run-time for lh, hh by a factor ~2 and ll by roughly
2531 // factor ~3 check if every scanned variable and resulting mass thermalised after N (default 10k) iterations
2532 // and then every M (default 1k) iterations do this by checking that the means of the split distributions is
2533 // comparable within X% (default 5%) of their sigma
2535 if (std::abs(m_fMEtP_split1->GetMean() - m_fMEtP_split2->GetMean()) <= m_fUseFloatStoppingComp * m_fMEtP_split1->GetRMS()) {
2536 if (std::abs(m_fMEtL_split1->GetMean() - m_fMEtL_split2->GetMean()) <=
2538 if (std::abs(m_fMnu1_split1->GetMean() - m_fMnu1_split2->GetMean()) <=
2540 if (std::abs(m_fMnu2_split1->GetMean() - m_fMnu2_split2->GetMean()) <=
2542 if (std::abs(m_fPhi1_split1->GetMean() - m_fPhi1_split2->GetMean()) <=
2544 if (std::abs(m_fPhi2_split1->GetMean() - m_fPhi2_split2->GetMean()) <=
2546 if (std::abs(m_fMmass_split1->GetMean() - m_fMmass_split2->GetMean()) <=
2548 return false;
2549 }
2550 }
2551 }
2552 }
2553 }
2554 }
2555 }
2556 }
2557
2559 // as long as no solution found need to randomise on the full parameter
2560 // space
2561
2562 // cut the corners in MissingET (not optimised at all)
2563 // not needed if distribution is already gaussian
2564 do {
2567 } while (!checkMEtInRange());
2568
2569 if (m_scanMnu1) {
2571 }
2572
2573 if (m_scanMnu2) {
2575 }
2576
2579
2580 return true;
2581 }
2582
2583 // here the real markov chain takes place : "propose" the new point
2584 // note that if one parameter goes outside range, this should not be fixed
2585 // here but later in handleSolution, otherwise would cause a bias
2586
2587 // m_MEtP0 etc... also store the position of the previous Markov Chain step,
2588 // which is needed by the algorithm
2589 m_MEtP0 = m_MEtP;
2590 m_MEtL0 = m_MEtL;
2591
2593
2595
2596 if (m_scanMnu1) {
2597 m_Mnu10 = m_Mnu1;
2599 }
2600
2601 if (m_scanMnu2) {
2602 m_Mnu20 = m_Mnu2;
2604 }
2605
2606 m_Phi10 = m_Phi1;
2608
2609 m_Phi20 = m_Phi2;
2610
2612
2613 return true;
2614}

◆ TailCleanUp()

int MissingMassCalculator::TailCleanUp ( const PtEtaPhiMVector & vis1,
const PtEtaPhiMVector & nu1,
const PtEtaPhiMVector & vis2,
const PtEtaPhiMVector & nu2,
const double & mmc_mass,
const double & vis_mass,
const double & eff_mass,
const double & dphiTT )
inlineprotected

Definition at line 1990 of file MissingMassCalculator.cxx.

1995 {
1996
1997 int pass_code = 1;
1998 if (preparedInput.m_fUseTailCleanup == 0)
1999 return pass_code;
2000
2001 // the Clean-up cuts are specifically for rel16 analyses.
2002 // the will change in rel17 analyses and after the MMC is updated
2003
2004 if (preparedInput.m_tauTypes == TauTypes::ll) // lepton-lepton channel
2005 {
2006 const double MrecoMvis = mmc_mass / vis_mass;
2007 if (MrecoMvis > 2.6)
2008 return 0;
2009 const double MrecoMeff = mmc_mass / eff_mass;
2010 if (MrecoMeff > 1.9)
2011 return 0;
2012 const double e1p1 = nu1.E() / vis1.P();
2013 const double e2p2 = nu2.E() / vis2.P();
2014 if ((e1p1 + e2p2) > 4.5)
2015 return 0;
2016 if (e2p2 > 4.0)
2017 return 0;
2018 if (e1p1 > 3.0)
2019 return 0;
2020 }
2021
2022 //-------- these are new cuts for lep-had analysis for Moriond
2023 if (preparedInput.m_tauTypes == TauTypes::lh) // lepton-hadron channel
2024 {
2025
2030 return pass_code; // don't use TailCleanup for 8 & 13 TeV data
2031
2032 //--------- leave code uncommented to avoid Compilation warnings
2033 if (Prob->GetUseHT()) {
2034 const double MrecoMvis = mmc_mass / vis_mass;
2035 const double MrecoMeff = mmc_mass / eff_mass;
2036 const double x = dphiTT > 1.5 ? dphiTT : 1.5;
2037 if ((MrecoMeff + MrecoMvis) > 5.908 - 1.881 * x + 0.2995 * x * x)
2038 return 0;
2039 }
2040 }
2041 return pass_code;
2042}

Member Data Documentation

◆ m_beamEnergy

double DiTauMassTools::MissingMassCalculator::m_beamEnergy {}
private

Definition at line 99 of file MissingMassCalculator.h.

99{}, m_beamEnergy{}; // number of sigmas for MET-scan

◆ m_cosPhi1

double DiTauMassTools::MissingMassCalculator::m_cosPhi1 {}
private

Definition at line 176 of file MissingMassCalculator.h.

176{}, m_cosPhi2{}, m_sinPhi1{}, m_sinPhi2{};

◆ m_cosPhi2

double DiTauMassTools::MissingMassCalculator::m_cosPhi2 {}
private

Definition at line 176 of file MissingMassCalculator.h.

176{}, m_cosPhi2{}, m_sinPhi1{}, m_sinPhi2{};

◆ m_debugThisIteration

bool DiTauMassTools::MissingMassCalculator::m_debugThisIteration
private

Definition at line 93 of file MissingMassCalculator.h.

◆ m_dRmax_tau

double DiTauMassTools::MissingMassCalculator::m_dRmax_tau {}
private

Definition at line 262 of file MissingMassCalculator.h.

262{}; // maximum dR(nu-visTau)

◆ m_dTheta3d_binMax

double DiTauMassTools::MissingMassCalculator::m_dTheta3d_binMax {}
private

Definition at line 261 of file MissingMassCalculator.h.

261{}; // maximum step size for dTheta3D

◆ m_dTheta3d_binMin

double DiTauMassTools::MissingMassCalculator::m_dTheta3d_binMin {}
private

Definition at line 260 of file MissingMassCalculator.h.

260{}; // minimal step size for dTheta3D

◆ m_E2v1

double DiTauMassTools::MissingMassCalculator::m_E2v1 {}
private

Definition at line 192 of file MissingMassCalculator.h.

192{};

◆ m_E2v2

double DiTauMassTools::MissingMassCalculator::m_E2v2 {}
private

Definition at line 193 of file MissingMassCalculator.h.

193{};

◆ m_ET2v1

double DiTauMassTools::MissingMassCalculator::m_ET2v1 {}
private

Definition at line 190 of file MissingMassCalculator.h.

190{};

◆ m_ET2v2

double DiTauMassTools::MissingMassCalculator::m_ET2v2 {}
private

Definition at line 191 of file MissingMassCalculator.h.

191{};

◆ m_eTau1

double DiTauMassTools::MissingMassCalculator::m_eTau1 {}
private

Definition at line 145 of file MissingMassCalculator.h.

145{}, m_eTau2{};

◆ m_eTau10

double DiTauMassTools::MissingMassCalculator::m_eTau10 {}
private

Definition at line 146 of file MissingMassCalculator.h.

146{}, m_eTau20{};

◆ m_eTau1Max

double DiTauMassTools::MissingMassCalculator::m_eTau1Max {}
private

◆ m_eTau1Min

double DiTauMassTools::MissingMassCalculator::m_eTau1Min {}
private

Definition at line 156 of file MissingMassCalculator.h.

◆ m_eTau1Proposal

double DiTauMassTools::MissingMassCalculator::m_eTau1Proposal {}
private

◆ m_eTau1Range

double DiTauMassTools::MissingMassCalculator::m_eTau1Range {}
private

Definition at line 156 of file MissingMassCalculator.h.

◆ m_eTau2

double DiTauMassTools::MissingMassCalculator::m_eTau2 {}
private

Definition at line 145 of file MissingMassCalculator.h.

145{}, m_eTau2{};

◆ m_eTau20

double DiTauMassTools::MissingMassCalculator::m_eTau20 {}
private

Definition at line 146 of file MissingMassCalculator.h.

146{}, m_eTau20{};

◆ m_eTau2Max

double DiTauMassTools::MissingMassCalculator::m_eTau2Max {}
private

◆ m_eTau2Min

double DiTauMassTools::MissingMassCalculator::m_eTau2Min {}
private

Definition at line 157 of file MissingMassCalculator.h.

◆ m_eTau2Proposal

double DiTauMassTools::MissingMassCalculator::m_eTau2Proposal {}
private

Definition at line 154 of file MissingMassCalculator.h.

154{},m_eTau2Proposal{};

◆ m_eTau2Range

double DiTauMassTools::MissingMassCalculator::m_eTau2Range {}
private

Definition at line 157 of file MissingMassCalculator.h.

◆ m_Ev1

double DiTauMassTools::MissingMassCalculator::m_Ev1 {}
private

Definition at line 195 of file MissingMassCalculator.h.

195{};

◆ m_Ev2

double DiTauMassTools::MissingMassCalculator::m_Ev2 {}
private

Definition at line 194 of file MissingMassCalculator.h.

194{};

◆ m_eventNumber

int DiTauMassTools::MissingMassCalculator::m_eventNumber {}
private

Definition at line 109 of file MissingMassCalculator.h.

109{};

◆ m_fDitauStuffFit

DitauStuff DiTauMassTools::MissingMassCalculator::m_fDitauStuffFit
private

Definition at line 249 of file MissingMassCalculator.h.

◆ m_fDitauStuffHisto

DitauStuff DiTauMassTools::MissingMassCalculator::m_fDitauStuffHisto
private

Definition at line 250 of file MissingMassCalculator.h.

◆ m_fFitting

TF1* DiTauMassTools::MissingMassCalculator::m_fFitting {}
private

Definition at line 234 of file MissingMassCalculator.h.

234{};

◆ m_fMEtL_all

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fMEtL_all
private

Definition at line 203 of file MissingMassCalculator.h.

◆ m_fMEtL_split1

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fMEtL_split1
private

Definition at line 221 of file MissingMassCalculator.h.

◆ m_fMEtL_split2

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fMEtL_split2
private

Definition at line 228 of file MissingMassCalculator.h.

◆ m_fMEtP_all

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fMEtP_all
private

Definition at line 202 of file MissingMassCalculator.h.

◆ m_fMEtP_split1

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fMEtP_split1
private

Definition at line 220 of file MissingMassCalculator.h.

◆ m_fMEtP_split2

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fMEtP_split2
private

Definition at line 227 of file MissingMassCalculator.h.

◆ m_fMetx

TH1F* DiTauMassTools::MissingMassCalculator::m_fMetx {}
private

Definition at line 240 of file MissingMassCalculator.h.

240{};

◆ m_fMety

TH1F* DiTauMassTools::MissingMassCalculator::m_fMety {}
private

Definition at line 241 of file MissingMassCalculator.h.

241{};

◆ m_fMfit_all

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fMfit_all
private

Definition at line 201 of file MissingMassCalculator.h.

◆ m_fMfit_allGraph

std::shared_ptr<TGraph> DiTauMassTools::MissingMassCalculator::m_fMfit_allGraph
private

Definition at line 208 of file MissingMassCalculator.h.

◆ m_fMfit_allNoWeight

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fMfit_allNoWeight
private

Definition at line 209 of file MissingMassCalculator.h.

◆ m_fMmass_split1

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fMmass_split1
private

Definition at line 219 of file MissingMassCalculator.h.

◆ m_fMmass_split2

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fMmass_split2
private

Definition at line 226 of file MissingMassCalculator.h.

◆ m_fMnu1

TH1F* DiTauMassTools::MissingMassCalculator::m_fMnu1 {}
private

Definition at line 238 of file MissingMassCalculator.h.

238{};

◆ m_fMnu1_all

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fMnu1_all
private

Definition at line 204 of file MissingMassCalculator.h.

◆ m_fMnu1_split1

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fMnu1_split1
private

Definition at line 222 of file MissingMassCalculator.h.

◆ m_fMnu1_split2

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fMnu1_split2
private

Definition at line 229 of file MissingMassCalculator.h.

◆ m_fMnu2

TH1F* DiTauMassTools::MissingMassCalculator::m_fMnu2 {}
private

Definition at line 239 of file MissingMassCalculator.h.

239{};

◆ m_fMnu2_all

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fMnu2_all
private

Definition at line 205 of file MissingMassCalculator.h.

◆ m_fMnu2_split1

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fMnu2_split1
private

Definition at line 223 of file MissingMassCalculator.h.

◆ m_fMnu2_split2

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fMnu2_split2
private

Definition at line 230 of file MissingMassCalculator.h.

◆ m_fPhi1

TH1F* DiTauMassTools::MissingMassCalculator::m_fPhi1 {}
private

Definition at line 236 of file MissingMassCalculator.h.

236{};

◆ m_fPhi1_all

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fPhi1_all
private

Definition at line 206 of file MissingMassCalculator.h.

◆ m_fPhi1_split1

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fPhi1_split1
private

Definition at line 224 of file MissingMassCalculator.h.

◆ m_fPhi1_split2

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fPhi1_split2
private

Definition at line 231 of file MissingMassCalculator.h.

◆ m_fPhi2

TH1F* DiTauMassTools::MissingMassCalculator::m_fPhi2 {}
private

Definition at line 237 of file MissingMassCalculator.h.

237{};

◆ m_fPhi2_all

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fPhi2_all
private

Definition at line 207 of file MissingMassCalculator.h.

◆ m_fPhi2_split1

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fPhi2_split1
private

Definition at line 225 of file MissingMassCalculator.h.

◆ m_fPhi2_split2

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fPhi2_split2
private

Definition at line 232 of file MissingMassCalculator.h.

◆ m_fPXfit1

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fPXfit1
private

Definition at line 211 of file MissingMassCalculator.h.

◆ m_fPXfit2

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fPXfit2
private

Definition at line 214 of file MissingMassCalculator.h.

◆ m_fPYfit1

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fPYfit1
private

Definition at line 212 of file MissingMassCalculator.h.

◆ m_fPYfit2

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fPYfit2
private

Definition at line 215 of file MissingMassCalculator.h.

◆ m_fPZfit1

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fPZfit1
private

Definition at line 213 of file MissingMassCalculator.h.

◆ m_fPZfit2

std::shared_ptr<TH1F> DiTauMassTools::MissingMassCalculator::m_fPZfit2
private

Definition at line 216 of file MissingMassCalculator.h.

◆ m_fTauProb

TH1F* DiTauMassTools::MissingMassCalculator::m_fTauProb {}
private

Definition at line 243 of file MissingMassCalculator.h.

243{};

◆ m_fTheta3D

TH1F* DiTauMassTools::MissingMassCalculator::m_fTheta3D {}
private

Definition at line 242 of file MissingMassCalculator.h.

242{};

◆ m_fullParamSpaceScan

bool DiTauMassTools::MissingMassCalculator::m_fullParamSpaceScan {}
private

Definition at line 158 of file MissingMassCalculator.h.

158{};

◆ m_fUseEfficiencyRecovery

bool DiTauMassTools::MissingMassCalculator::m_fUseEfficiencyRecovery {}
private

Definition at line 67 of file MissingMassCalculator.h.

67{}; // switch to turn ON/OFF re-fit in order to recover efficiency

◆ m_fUseFloatStopping

bool DiTauMassTools::MissingMassCalculator::m_fUseFloatStopping {}
private

Definition at line 68 of file MissingMassCalculator.h.

68{}; // switch to turn ON/OFF floating stopping criterion

◆ m_fUseFloatStoppingCheckFreq

int DiTauMassTools::MissingMassCalculator::m_fUseFloatStoppingCheckFreq {}
private

Definition at line 70 of file MissingMassCalculator.h.

70{};

◆ m_fUseFloatStoppingComp

double DiTauMassTools::MissingMassCalculator::m_fUseFloatStoppingComp {}
private

Definition at line 71 of file MissingMassCalculator.h.

71{};

◆ m_fUseFloatStoppingMinIter

int DiTauMassTools::MissingMassCalculator::m_fUseFloatStoppingMinIter {}
private

Definition at line 69 of file MissingMassCalculator.h.

69{};

◆ m_iang1high

int DiTauMassTools::MissingMassCalculator::m_iang1high {}
private

◆ m_iang1low

int DiTauMassTools::MissingMassCalculator::m_iang1low {}
private

◆ m_iang2high

int DiTauMassTools::MissingMassCalculator::m_iang2high {}
private

◆ m_iang2low

int DiTauMassTools::MissingMassCalculator::m_iang2low {}
private

◆ m_iter0

int DiTauMassTools::MissingMassCalculator::m_iter0 {}
private

Definition at line 113 of file MissingMassCalculator.h.

113{};

◆ m_iter1

int DiTauMassTools::MissingMassCalculator::m_iter1 {}
private

◆ m_iter2

int DiTauMassTools::MissingMassCalculator::m_iter2 {}
private

◆ m_iter3

int DiTauMassTools::MissingMassCalculator::m_iter3 {}
private

◆ m_iter4

int DiTauMassTools::MissingMassCalculator::m_iter4 {}
private

◆ m_iter5

int DiTauMassTools::MissingMassCalculator::m_iter5 {}
private

◆ m_iterNsuc

int DiTauMassTools::MissingMassCalculator::m_iterNsuc {}
private

Definition at line 121 of file MissingMassCalculator.h.

121{};

◆ m_iterNuPV3

int DiTauMassTools::MissingMassCalculator::m_iterNuPV3 {}
private

Definition at line 114 of file MissingMassCalculator.h.

114{};

◆ m_iterTheta3d

int DiTauMassTools::MissingMassCalculator::m_iterTheta3d {}
private

Definition at line 102 of file MissingMassCalculator.h.

102{};

◆ m_lfvLeplepRefit

bool DiTauMassTools::MissingMassCalculator::m_lfvLeplepRefit
private

Definition at line 93 of file MissingMassCalculator.h.

◆ m_m2Nu1

double DiTauMassTools::MissingMassCalculator::m_m2Nu1 {}
private

Definition at line 188 of file MissingMassCalculator.h.

188{};

◆ m_m2Nu2

double DiTauMassTools::MissingMassCalculator::m_m2Nu2 {}
private

Definition at line 189 of file MissingMassCalculator.h.

189{};

◆ m_markovCountDuplicate

int DiTauMassTools::MissingMassCalculator::m_markovCountDuplicate {}
private

Definition at line 127 of file MissingMassCalculator.h.

127{};

◆ m_markovNAccept

int DiTauMassTools::MissingMassCalculator::m_markovNAccept {}
private

Definition at line 131 of file MissingMassCalculator.h.

131{};

◆ m_markovNFullScan

int DiTauMassTools::MissingMassCalculator::m_markovNFullScan {}
private

Definition at line 128 of file MissingMassCalculator.h.

128{};

◆ m_markovNRejectMetropolis

int DiTauMassTools::MissingMassCalculator::m_markovNRejectMetropolis {}
private

Definition at line 130 of file MissingMassCalculator.h.

130{};

◆ m_markovNRejectNoSol

int DiTauMassTools::MissingMassCalculator::m_markovNRejectNoSol {}
private

Definition at line 129 of file MissingMassCalculator.h.

129{};

◆ m_meanbinStop

double DiTauMassTools::MissingMassCalculator::m_meanbinStop {}
private

Definition at line 77 of file MissingMassCalculator.h.

77{};

◆ m_meanbinToBeEvaluated

bool DiTauMassTools::MissingMassCalculator::m_meanbinToBeEvaluated {}
private

Definition at line 125 of file MissingMassCalculator.h.

125{};

◆ m_Meff

double DiTauMassTools::MissingMassCalculator::m_Meff {}
private

Definition at line 196 of file MissingMassCalculator.h.

196{},m_Meff{};

◆ m_metCovPhiCos

double DiTauMassTools::MissingMassCalculator::m_metCovPhiCos {}
private

Definition at line 153 of file MissingMassCalculator.h.

153{},m_metCovPhiSin{};

◆ m_metCovPhiSin

double DiTauMassTools::MissingMassCalculator::m_metCovPhiSin {}
private

Definition at line 153 of file MissingMassCalculator.h.

153{},m_metCovPhiSin{};

◆ m_MEtL

double DiTauMassTools::MissingMassCalculator::m_MEtL {}
private

Definition at line 144 of file MissingMassCalculator.h.

144{},m_MEtP{},m_Phi1{},m_Phi2{},m_Mnu1{},m_Mnu2{};

◆ m_MEtL0

double DiTauMassTools::MissingMassCalculator::m_MEtL0 {}
private

Definition at line 147 of file MissingMassCalculator.h.

◆ m_MEtLMax

double DiTauMassTools::MissingMassCalculator::m_MEtLMax {}
private

Definition at line 149 of file MissingMassCalculator.h.

◆ m_MEtLMin

double DiTauMassTools::MissingMassCalculator::m_MEtLMin {}
private

Definition at line 148 of file MissingMassCalculator.h.

◆ m_MEtLRange

double DiTauMassTools::MissingMassCalculator::m_MEtLRange {}
private

◆ m_MEtLStep

◆ m_MEtP

double DiTauMassTools::MissingMassCalculator::m_MEtP {}
private

Definition at line 144 of file MissingMassCalculator.h.

144{},m_MEtP{},m_Phi1{},m_Phi2{},m_Mnu1{},m_Mnu2{};

◆ m_MEtP0

double DiTauMassTools::MissingMassCalculator::m_MEtP0 {}
private

Definition at line 147 of file MissingMassCalculator.h.

◆ m_MEtPMax

double DiTauMassTools::MissingMassCalculator::m_MEtPMax {}
private

Definition at line 149 of file MissingMassCalculator.h.

◆ m_MEtPMin

double DiTauMassTools::MissingMassCalculator::m_MEtPMin {}
private

Definition at line 148 of file MissingMassCalculator.h.

◆ m_MEtPRange

double DiTauMassTools::MissingMassCalculator::m_MEtPRange {}
private

◆ m_MEtProposal

double DiTauMassTools::MissingMassCalculator::m_MEtProposal {}
private

Definition at line 152 of file MissingMassCalculator.h.

◆ m_MEtPStep

double DiTauMassTools::MissingMassCalculator::m_MEtPStep {}
private

Definition at line 150 of file MissingMassCalculator.h.

◆ m_mmcCalibrationSet

MMCCalibrationSet::e DiTauMassTools::MissingMassCalculator::m_mmcCalibrationSet {}
private

Definition at line 65 of file MissingMassCalculator.h.

65{};

◆ m_Mnu1

double DiTauMassTools::MissingMassCalculator::m_Mnu1 {}
private

Definition at line 144 of file MissingMassCalculator.h.

144{},m_MEtP{},m_Phi1{},m_Phi2{},m_Mnu1{},m_Mnu2{};

◆ m_Mnu10

double DiTauMassTools::MissingMassCalculator::m_Mnu10 {}
private

Definition at line 147 of file MissingMassCalculator.h.

◆ m_Mnu1Exclude

bool DiTauMassTools::MissingMassCalculator::m_Mnu1Exclude {}
private

Definition at line 159 of file MissingMassCalculator.h.

159{};

◆ m_Mnu1ExcludeMax

double DiTauMassTools::MissingMassCalculator::m_Mnu1ExcludeMax {}
private

◆ m_Mnu1ExcludeMin

double DiTauMassTools::MissingMassCalculator::m_Mnu1ExcludeMin {}
private

Definition at line 173 of file MissingMassCalculator.h.

◆ m_Mnu1ExcludeRange

double DiTauMassTools::MissingMassCalculator::m_Mnu1ExcludeRange {}
private

Definition at line 173 of file MissingMassCalculator.h.

◆ m_Mnu1Max

double DiTauMassTools::MissingMassCalculator::m_Mnu1Max {}
private

Definition at line 149 of file MissingMassCalculator.h.

◆ m_Mnu1Min

double DiTauMassTools::MissingMassCalculator::m_Mnu1Min {}
private

Definition at line 148 of file MissingMassCalculator.h.

◆ m_Mnu1Range

double DiTauMassTools::MissingMassCalculator::m_Mnu1Range {}
private

◆ m_Mnu1Step

double DiTauMassTools::MissingMassCalculator::m_Mnu1Step {}
private

Definition at line 150 of file MissingMassCalculator.h.

◆ m_Mnu1XMax

double DiTauMassTools::MissingMassCalculator::m_Mnu1XMax {}
private

◆ m_Mnu1XMin

double DiTauMassTools::MissingMassCalculator::m_Mnu1XMin {}
private

Definition at line 174 of file MissingMassCalculator.h.

◆ m_Mnu1XRange

double DiTauMassTools::MissingMassCalculator::m_Mnu1XRange {}
private

Definition at line 174 of file MissingMassCalculator.h.

◆ m_Mnu2

double DiTauMassTools::MissingMassCalculator::m_Mnu2 {}
private

Definition at line 144 of file MissingMassCalculator.h.

144{},m_MEtP{},m_Phi1{},m_Phi2{},m_Mnu1{},m_Mnu2{};

◆ m_Mnu20

double DiTauMassTools::MissingMassCalculator::m_Mnu20 {}
private

Definition at line 147 of file MissingMassCalculator.h.

◆ m_Mnu2Max

double DiTauMassTools::MissingMassCalculator::m_Mnu2Max {}
private

Definition at line 149 of file MissingMassCalculator.h.

◆ m_Mnu2Min

double DiTauMassTools::MissingMassCalculator::m_Mnu2Min {}
private

Definition at line 148 of file MissingMassCalculator.h.

◆ m_Mnu2Range

double DiTauMassTools::MissingMassCalculator::m_Mnu2Range {}
private

◆ m_Mnu2Step

double DiTauMassTools::MissingMassCalculator::m_Mnu2Step {}
private

Definition at line 150 of file MissingMassCalculator.h.

◆ m_MnuProposal

double DiTauMassTools::MissingMassCalculator::m_MnuProposal {}
private

Definition at line 152 of file MissingMassCalculator.h.

◆ m_MnuScanRange

double DiTauMassTools::MissingMassCalculator::m_MnuScanRange {}
private

Definition at line 264 of file MissingMassCalculator.h.

264{}; // range of M(nunu) scan; M(nunu) range can be affected by selection cuts

◆ m_mTau

double DiTauMassTools::MissingMassCalculator::m_mTau {}
private

Definition at line 143 of file MissingMassCalculator.h.

143{},m_mTau2{};

◆ m_mTau2

double DiTauMassTools::MissingMassCalculator::m_mTau2 {}
private

Definition at line 143 of file MissingMassCalculator.h.

143{},m_mTau2{};

◆ m_mtautauFinalSolOldVec

std::vector<double> DiTauMassTools::MissingMassCalculator::m_mtautauFinalSolOldVec
private

Definition at line 162 of file MissingMassCalculator.h.

◆ m_mtautauFinalSolVec

std::vector<double> DiTauMassTools::MissingMassCalculator::m_mtautauFinalSolVec
private

Definition at line 168 of file MissingMassCalculator.h.

◆ m_mtautauSum

double DiTauMassTools::MissingMassCalculator::m_mtautauSum {}
private

Definition at line 107 of file MissingMassCalculator.h.

107{};

◆ m_Mvis

double DiTauMassTools::MissingMassCalculator::m_Mvis {}
private

Definition at line 196 of file MissingMassCalculator.h.

196{},m_Meff{};

◆ m_nCallprobCalculatorV9fast

int DiTauMassTools::MissingMassCalculator::m_nCallprobCalculatorV9fast
private

Definition at line 95 of file MissingMassCalculator.h.

◆ m_niter_fit1

int DiTauMassTools::MissingMassCalculator::m_niter_fit1 {}
private

Definition at line 252 of file MissingMassCalculator.h.

252{}; // number of iterations for dR-dPhi scan

◆ m_niter_fit2

int DiTauMassTools::MissingMassCalculator::m_niter_fit2 {}
private

Definition at line 253 of file MissingMassCalculator.h.

253{}; // number of iterations for MET-scan

◆ m_niter_fit3

int DiTauMassTools::MissingMassCalculator::m_niter_fit3 {}
private

Definition at line 254 of file MissingMassCalculator.h.

254{}; // number of iterations for Mnu-scan

◆ m_NiterRandom

int DiTauMassTools::MissingMassCalculator::m_NiterRandom {}
private

Definition at line 255 of file MissingMassCalculator.h.

255{}; // number of random iterations (for lh, multiply or divide by 10 for ll and hh)

◆ m_niterRandomLocal

int DiTauMassTools::MissingMassCalculator::m_niterRandomLocal {}
private

Definition at line 74 of file MissingMassCalculator.h.

74{};

◆ m_nosol1

int DiTauMassTools::MissingMassCalculator::m_nosol1 {}
private

Definition at line 119 of file MissingMassCalculator.h.

119{};

◆ m_nosol2

int DiTauMassTools::MissingMassCalculator::m_nosol2 {}
private

Definition at line 120 of file MissingMassCalculator.h.

120{};

◆ m_nsigma_METscan

double DiTauMassTools::MissingMassCalculator::m_nsigma_METscan {}
private

Definition at line 97 of file MissingMassCalculator.h.

◆ m_nsigma_METscan2

double DiTauMassTools::MissingMassCalculator::m_nsigma_METscan2 {}
private

Definition at line 97 of file MissingMassCalculator.h.

◆ m_nsigma_METscan_hh

double DiTauMassTools::MissingMassCalculator::m_nsigma_METscan_hh {}
private

◆ m_nsigma_METscan_lfv_lh

double DiTauMassTools::MissingMassCalculator::m_nsigma_METscan_lfv_lh {}
private

Definition at line 99 of file MissingMassCalculator.h.

99{}, m_beamEnergy{}; // number of sigmas for MET-scan

◆ m_nsigma_METscan_lfv_ll

double DiTauMassTools::MissingMassCalculator::m_nsigma_METscan_lfv_ll {}
private

◆ m_nsigma_METscan_lh

double DiTauMassTools::MissingMassCalculator::m_nsigma_METscan_lh {}
private

◆ m_nsigma_METscan_ll

double DiTauMassTools::MissingMassCalculator::m_nsigma_METscan_ll {}
private

Definition at line 97 of file MissingMassCalculator.h.

◆ m_nsol

int DiTauMassTools::MissingMassCalculator::m_nsol {}
private

Definition at line 166 of file MissingMassCalculator.h.

166{};

◆ m_nsolfinalmax

int DiTauMassTools::MissingMassCalculator::m_nsolfinalmax {}
private

Definition at line 73 of file MissingMassCalculator.h.

73{};

◆ m_nsolmax

int DiTauMassTools::MissingMassCalculator::m_nsolmax
private

Definition at line 73 of file MissingMassCalculator.h.

◆ m_nsolOld

int DiTauMassTools::MissingMassCalculator::m_nsolOld {}
private

Definition at line 160 of file MissingMassCalculator.h.

160{};

◆ m_NsucStop

int DiTauMassTools::MissingMassCalculator::m_NsucStop {}
private

Definition at line 256 of file MissingMassCalculator.h.

256{};

◆ m_nsucStop

int DiTauMassTools::MissingMassCalculator::m_nsucStop {}
private

Definition at line 75 of file MissingMassCalculator.h.

75{};

◆ m_nu1FinalSolOldVec

std::vector<PtEtaPhiMVector> DiTauMassTools::MissingMassCalculator::m_nu1FinalSolOldVec
private

Definition at line 163 of file MissingMassCalculator.h.

◆ m_nu1FinalSolVec

std::vector<PtEtaPhiMVector> DiTauMassTools::MissingMassCalculator::m_nu1FinalSolVec
private

Definition at line 169 of file MissingMassCalculator.h.

◆ m_nu2FinalSolOldVec

std::vector<PtEtaPhiMVector> DiTauMassTools::MissingMassCalculator::m_nu2FinalSolOldVec
private

Definition at line 164 of file MissingMassCalculator.h.

◆ m_nu2FinalSolVec

std::vector<PtEtaPhiMVector> DiTauMassTools::MissingMassCalculator::m_nu2FinalSolVec
private

Definition at line 170 of file MissingMassCalculator.h.

◆ m_nuvec1_tmp

std::vector<PtEtaPhiMVector> DiTauMassTools::MissingMassCalculator::m_nuvec1_tmp
private

Definition at line 88 of file MissingMassCalculator.h.

◆ m_nuvec2_tmp

std::vector<PtEtaPhiMVector> DiTauMassTools::MissingMassCalculator::m_nuvec2_tmp
private

Definition at line 89 of file MissingMassCalculator.h.

◆ m_nuvecsol1

std::vector<PtEtaPhiMVector> DiTauMassTools::MissingMassCalculator::m_nuvecsol1
private

Definition at line 80 of file MissingMassCalculator.h.

◆ m_nuvecsol2

std::vector<PtEtaPhiMVector> DiTauMassTools::MissingMassCalculator::m_nuvecsol2
private

Definition at line 81 of file MissingMassCalculator.h.

◆ m_Phi1

double DiTauMassTools::MissingMassCalculator::m_Phi1 {}
private

Definition at line 144 of file MissingMassCalculator.h.

144{},m_MEtP{},m_Phi1{},m_Phi2{},m_Mnu1{},m_Mnu2{};

◆ m_Phi10

double DiTauMassTools::MissingMassCalculator::m_Phi10 {}
private

Definition at line 147 of file MissingMassCalculator.h.

◆ m_Phi1Max

double DiTauMassTools::MissingMassCalculator::m_Phi1Max {}
private

Definition at line 149 of file MissingMassCalculator.h.

◆ m_Phi1Min

double DiTauMassTools::MissingMassCalculator::m_Phi1Min {}
private

Definition at line 148 of file MissingMassCalculator.h.

◆ m_Phi1Range

double DiTauMassTools::MissingMassCalculator::m_Phi1Range {}
private

◆ m_Phi1Step

double DiTauMassTools::MissingMassCalculator::m_Phi1Step {}
private

Definition at line 150 of file MissingMassCalculator.h.

◆ m_Phi2

double DiTauMassTools::MissingMassCalculator::m_Phi2 {}
private

Definition at line 144 of file MissingMassCalculator.h.

144{},m_MEtP{},m_Phi1{},m_Phi2{},m_Mnu1{},m_Mnu2{};

◆ m_Phi20

double DiTauMassTools::MissingMassCalculator::m_Phi20 {}
private

Definition at line 147 of file MissingMassCalculator.h.

◆ m_Phi2Max

double DiTauMassTools::MissingMassCalculator::m_Phi2Max {}
private

Definition at line 149 of file MissingMassCalculator.h.

◆ m_Phi2Min

double DiTauMassTools::MissingMassCalculator::m_Phi2Min {}
private

Definition at line 148 of file MissingMassCalculator.h.

◆ m_Phi2Range

double DiTauMassTools::MissingMassCalculator::m_Phi2Range {}
private

◆ m_Phi2Step

double DiTauMassTools::MissingMassCalculator::m_Phi2Step {}
private

Definition at line 150 of file MissingMassCalculator.h.

◆ m_PhiProposal

double DiTauMassTools::MissingMassCalculator::m_PhiProposal {}
private

Definition at line 152 of file MissingMassCalculator.h.

◆ m_PrintmInvWidth2Error

double DiTauMassTools::MissingMassCalculator::m_PrintmInvWidth2Error {}
private

Definition at line 135 of file MissingMassCalculator.h.

135{};

◆ m_PrintmMaxError

double DiTauMassTools::MissingMassCalculator::m_PrintmMaxError {}
private

Definition at line 133 of file MissingMassCalculator.h.

133{};

◆ m_PrintmMeanError

double DiTauMassTools::MissingMassCalculator::m_PrintmMeanError {}
private

Definition at line 134 of file MissingMassCalculator.h.

134{};

◆ m_prob_tmp

double DiTauMassTools::MissingMassCalculator::m_prob_tmp {}
private

Definition at line 104 of file MissingMassCalculator.h.

104{};

◆ m_probFinalSolOldVec

std::vector<double> DiTauMassTools::MissingMassCalculator::m_probFinalSolOldVec
private

Definition at line 161 of file MissingMassCalculator.h.

◆ m_probFinalSolVec

std::vector<double> DiTauMassTools::MissingMassCalculator::m_probFinalSolVec
private

Definition at line 167 of file MissingMassCalculator.h.

◆ m_ProposalTryEtau

double DiTauMassTools::MissingMassCalculator::m_ProposalTryEtau {}
private

Definition at line 140 of file MissingMassCalculator.h.

140{};

◆ m_proposalTryMEt

double DiTauMassTools::MissingMassCalculator::m_proposalTryMEt {}
private

Definition at line 137 of file MissingMassCalculator.h.

137{};

◆ m_ProposalTryMnu

double DiTauMassTools::MissingMassCalculator::m_ProposalTryMnu {}
private

Definition at line 139 of file MissingMassCalculator.h.

139{};

◆ m_ProposalTryPhi

double DiTauMassTools::MissingMassCalculator::m_ProposalTryPhi {}
private

Definition at line 138 of file MissingMassCalculator.h.

138{};

◆ m_randomGen

TRandom2 DiTauMassTools::MissingMassCalculator::m_randomGen
private

Definition at line 63 of file MissingMassCalculator.h.

◆ m_reRunWithBestMET

bool DiTauMassTools::MissingMassCalculator::m_reRunWithBestMET {}
private

Definition at line 197 of file MissingMassCalculator.h.

197{};

◆ m_RMSStop

int DiTauMassTools::MissingMassCalculator::m_RMSStop {}
private

Definition at line 257 of file MissingMassCalculator.h.

257{};

◆ m_rmsStop

int DiTauMassTools::MissingMassCalculator::m_rmsStop {}
private

Definition at line 76 of file MissingMassCalculator.h.

76{};

◆ m_RndmSeedAltering

int DiTauMassTools::MissingMassCalculator::m_RndmSeedAltering {}
private

Definition at line 258 of file MissingMassCalculator.h.

258{}; // reset seed (not necessary by default)

◆ m_SaveLlhHisto

bool DiTauMassTools::MissingMassCalculator::m_SaveLlhHisto
private

Definition at line 93 of file MissingMassCalculator.h.

◆ m_scanMnu1

bool DiTauMassTools::MissingMassCalculator::m_scanMnu1 {}
private

Definition at line 178 of file MissingMassCalculator.h.

178{},m_scanMnu2{};

◆ m_scanMnu2

bool DiTauMassTools::MissingMassCalculator::m_scanMnu2 {}
private

Definition at line 178 of file MissingMassCalculator.h.

178{},m_scanMnu2{};

◆ m_seed

int DiTauMassTools::MissingMassCalculator::m_seed {}
private

Definition at line 110 of file MissingMassCalculator.h.

110{};

◆ m_sinPhi1

double DiTauMassTools::MissingMassCalculator::m_sinPhi1 {}
private

Definition at line 176 of file MissingMassCalculator.h.

176{}, m_cosPhi2{}, m_sinPhi1{}, m_sinPhi2{};

◆ m_sinPhi2

double DiTauMassTools::MissingMassCalculator::m_sinPhi2 {}
private

Definition at line 176 of file MissingMassCalculator.h.

176{}, m_cosPhi2{}, m_sinPhi1{}, m_sinPhi2{};

◆ m_switch1

bool DiTauMassTools::MissingMassCalculator::m_switch1 {}
private

Definition at line 122 of file MissingMassCalculator.h.

122{};

◆ m_switch2

bool DiTauMassTools::MissingMassCalculator::m_switch2 {}
private

Definition at line 123 of file MissingMassCalculator.h.

123{};

◆ m_tautau_tmp

PtEtaPhiMVector DiTauMassTools::MissingMassCalculator::m_tautau_tmp
private

Definition at line 91 of file MissingMassCalculator.h.

◆ m_tauVec1

PtEtaPhiMVector DiTauMassTools::MissingMassCalculator::m_tauVec1
private

Definition at line 180 of file MissingMassCalculator.h.

◆ m_tauVec1E

double DiTauMassTools::MissingMassCalculator::m_tauVec1E {}
private

Definition at line 186 of file MissingMassCalculator.h.

186{};

◆ m_tauVec1M

double DiTauMassTools::MissingMassCalculator::m_tauVec1M {}
private

Definition at line 182 of file MissingMassCalculator.h.

182{}, m_tauVec2M{};

◆ m_tauVec1P

double DiTauMassTools::MissingMassCalculator::m_tauVec1P {}
private

Definition at line 185 of file MissingMassCalculator.h.

185{}, m_tauVec2P{};

◆ m_tauVec1Phi

double DiTauMassTools::MissingMassCalculator::m_tauVec1Phi {}
private

Definition at line 181 of file MissingMassCalculator.h.

181{}, m_tauVec2Phi{};

◆ m_tauVec1Px

double DiTauMassTools::MissingMassCalculator::m_tauVec1Px {}
private

Definition at line 183 of file MissingMassCalculator.h.

183{}, m_tauVec1Py{}, m_tauVec1Pz{};

◆ m_tauVec1Py

double DiTauMassTools::MissingMassCalculator::m_tauVec1Py {}
private

Definition at line 183 of file MissingMassCalculator.h.

183{}, m_tauVec1Py{}, m_tauVec1Pz{};

◆ m_tauVec1Pz

double DiTauMassTools::MissingMassCalculator::m_tauVec1Pz {}
private

Definition at line 183 of file MissingMassCalculator.h.

183{}, m_tauVec1Py{}, m_tauVec1Pz{};

◆ m_tauVec2

PtEtaPhiMVector DiTauMassTools::MissingMassCalculator::m_tauVec2
private

Definition at line 180 of file MissingMassCalculator.h.

◆ m_tauVec2E

double DiTauMassTools::MissingMassCalculator::m_tauVec2E {}
private

Definition at line 187 of file MissingMassCalculator.h.

187{};

◆ m_tauVec2M

double DiTauMassTools::MissingMassCalculator::m_tauVec2M {}
private

Definition at line 182 of file MissingMassCalculator.h.

182{}, m_tauVec2M{};

◆ m_tauVec2P

double DiTauMassTools::MissingMassCalculator::m_tauVec2P {}
private

Definition at line 185 of file MissingMassCalculator.h.

185{}, m_tauVec2P{};

◆ m_tauVec2Phi

double DiTauMassTools::MissingMassCalculator::m_tauVec2Phi {}
private

Definition at line 181 of file MissingMassCalculator.h.

181{}, m_tauVec2Phi{};

◆ m_tauVec2Px

double DiTauMassTools::MissingMassCalculator::m_tauVec2Px {}
private

Definition at line 184 of file MissingMassCalculator.h.

184{}, m_tauVec2Py{}, m_tauVec2Pz{};

◆ m_tauVec2Py

double DiTauMassTools::MissingMassCalculator::m_tauVec2Py {}
private

Definition at line 184 of file MissingMassCalculator.h.

184{}, m_tauVec2Py{}, m_tauVec2Pz{};

◆ m_tauVec2Pz

double DiTauMassTools::MissingMassCalculator::m_tauVec2Pz {}
private

Definition at line 184 of file MissingMassCalculator.h.

184{}, m_tauVec2Py{}, m_tauVec2Pz{};

◆ m_tauvecprob1

std::vector<double> DiTauMassTools::MissingMassCalculator::m_tauvecprob1
private

Definition at line 85 of file MissingMassCalculator.h.

◆ m_tauvecprob2

std::vector<double> DiTauMassTools::MissingMassCalculator::m_tauvecprob2
private

Definition at line 86 of file MissingMassCalculator.h.

◆ m_tauvecsol1

std::vector<PtEtaPhiMVector> DiTauMassTools::MissingMassCalculator::m_tauvecsol1
private

Definition at line 83 of file MissingMassCalculator.h.

◆ m_tauvecsol2

std::vector<PtEtaPhiMVector> DiTauMassTools::MissingMassCalculator::m_tauvecsol2
private

Definition at line 84 of file MissingMassCalculator.h.

◆ m_testdiscri1

int DiTauMassTools::MissingMassCalculator::m_testdiscri1 {}
private

Definition at line 117 of file MissingMassCalculator.h.

117{};

◆ m_testdiscri2

int DiTauMassTools::MissingMassCalculator::m_testdiscri2 {}
private

Definition at line 118 of file MissingMassCalculator.h.

118{};

◆ m_testptn1

int DiTauMassTools::MissingMassCalculator::m_testptn1 {}
private

Definition at line 115 of file MissingMassCalculator.h.

115{};

◆ m_testptn2

int DiTauMassTools::MissingMassCalculator::m_testptn2 {}
private

Definition at line 116 of file MissingMassCalculator.h.

116{};

◆ m_TLVdummy

PtEtaPhiMVector DiTauMassTools::MissingMassCalculator::m_TLVdummy
private

Definition at line 246 of file MissingMassCalculator.h.

◆ m_totalProbSum

double DiTauMassTools::MissingMassCalculator::m_totalProbSum {}
private

Definition at line 106 of file MissingMassCalculator.h.

106{};

◆ m_walkWeight

double DiTauMassTools::MissingMassCalculator::m_walkWeight {}
private

Definition at line 175 of file MissingMassCalculator.h.

175{};

◆ metvec_tmp

XYVector DiTauMassTools::MissingMassCalculator::metvec_tmp

Definition at line 421 of file MissingMassCalculator.h.

◆ OutputInfo

MissingMassOutput DiTauMassTools::MissingMassCalculator::OutputInfo

Definition at line 341 of file MissingMassCalculator.h.

◆ preparedInput

MissingMassInput DiTauMassTools::MissingMassCalculator::preparedInput

Definition at line 340 of file MissingMassCalculator.h.

◆ Prob

MissingMassProb* DiTauMassTools::MissingMassCalculator::Prob

Definition at line 342 of file MissingMassCalculator.h.


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