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 2665 of file MissingMassCalculator.cxx.

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

◆ checkMEtInRange()

bool MissingMassCalculator::checkMEtInRange ( )
inlineprotected

Definition at line 2706 of file MissingMassCalculator.cxx.

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

◆ 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 1014 of file MissingMassCalculator.cxx.

1014 {
1015
1016 // debugThisIteration=false;
1017 m_debugThisIteration = true;
1018
1019 int fit_code = 0; // 0==bad, 1==good
1022 OutputInfo.m_NTrials = 0;
1023 OutputInfo.m_NSuccesses = 0;
1024 OutputInfo.m_AveSolRMS = 0.;
1025
1026 //------- Settings -------------------------------
1027 int NiterMET = m_niter_fit2; // number of iterations for each MET scan loop
1028 int NiterMnu = m_niter_fit3; // number of iterations for Mnu loop
1029 const double Mtau = ParticleConstants::tauMassInMeV / GEV;
1030 double Mnu_binSize = m_MnuScanRange / NiterMnu;
1031
1032 double METresX = preparedInput.m_METsigmaL; // MET resolution in direction parallel to
1033 // leading jet, for MET scan
1034 double METresY = preparedInput.m_METsigmaP; // MET resolution in direction perpendicular to
1035 // leading jet, for MET scan
1036
1037 //-------- end of Settings
1038
1039 // if m_nsigma_METscan was not set by user, set to default values
1040 if(m_nsigma_METscan == -1){
1041 if (preparedInput.m_tauTypes == TauTypes::ll) { // both tau's are leptonic
1043 } else if (preparedInput.m_tauTypes == TauTypes::lh) { // lep had
1045 }
1046 }
1047
1048 double N_METsigma = m_nsigma_METscan; // number of sigmas for MET scan
1049 double METresX_binSize = 2 * N_METsigma * METresX / NiterMET;
1050 double METresY_binSize = 2 * N_METsigma * METresY / NiterMET;
1051
1052 int solution = 0;
1053
1054 std::vector<PtEtaPhiMVector> nu_vec;
1055
1056 m_totalProbSum = 0;
1057 m_mtautauSum = 0;
1058
1059 double metprob = 1.0;
1060 double sign_tmp = 0.0;
1061 double tauprob = 1.0;
1062 double totalProb = 0.0;
1063
1064 m_prob_tmp = 0.0;
1065
1066 double met_smear_x = 0.0;
1067 double met_smear_y = 0.0;
1068 double met_smearL = 0.0;
1069 double met_smearP = 0.0;
1070
1071 double angle1 = 0.0;
1072
1073 if (m_fMfit_all) {
1074 m_fMfit_all->Reset();
1075 }
1076 if (m_fMfit_allNoWeight) {
1077 m_fMfit_allNoWeight->Reset();
1078 }
1079 if (m_fPXfit1) {
1080 m_fPXfit1->Reset();
1081 }
1082 if (m_fPYfit1) {
1083 m_fPYfit1->Reset();
1084 }
1085 if (m_fPZfit1) {
1086 m_fPZfit1->Reset();
1087 }
1088
1089 int iter0 = 0;
1090 m_iter1 = 0;
1091 m_iter2 = 0;
1092 m_iter3 = 0;
1093 m_iter4 = 0;
1094
1095 const double met_coscovphi = cos(preparedInput.m_METcovphi);
1096 const double met_sincovphi = sin(preparedInput.m_METcovphi);
1097
1098 m_iang1low = 0;
1099 m_iang1high = 0;
1100
1101 // double Mvis=(tau_vec1+tau_vec2).M();
1102 // PtEtaPhiMVector met4vec(0.0,0.0,0.0,0.0);
1103 // met4vec.SetPxPyPzE(met_vec.X(),met_vec.Y(),0.0,met_vec.R());
1104 // double Meff=(tau_vec1+tau_vec2+met4vec).M();
1105 // double met_det=met_vec.R();
1106
1107 //---------------------------------------------
1108 if (preparedInput.m_tauTypes == TauTypes::ll) // dilepton case
1109 {
1110 if (preparedInput.m_fUseVerbose == 1) {
1111 Info("DiTauMassTools", "Running in dilepton mode");
1112 }
1113 double input_metX = preparedInput.m_MetVec.X();
1114 double input_metY = preparedInput.m_MetVec.Y();
1115
1116 PtEtaPhiMVector tau_tmp(0.0, 0.0, 0.0, 0.0);
1117 PtEtaPhiMVector lep_tmp(0.0, 0.0, 0.0, 0.0);
1118 int tau_type_tmp;
1119 int tau_ind = 0;
1120
1121 if (preparedInput.m_LFVmode == 1) // muon case: H->mu+tau(->ele) decays
1122 {
1123 if ((preparedInput.m_vistau1.M() > 0.05 &&
1124 preparedInput.m_vistau2.M() < 0.05) != refit) // choosing lepton from Higgs decay
1125 //When the mass calculator is rerun with refit==true the alternative lepton ordering is used
1126 {
1127 tau_tmp = preparedInput.m_vistau2;
1128 lep_tmp = preparedInput.m_vistau1;
1129 tau_type_tmp = preparedInput.m_type_visTau2;
1130 tau_ind = 2;
1131 } else {
1132 tau_tmp = preparedInput.m_vistau1;
1133 lep_tmp = preparedInput.m_vistau2;
1134 tau_type_tmp = preparedInput.m_type_visTau1;
1135 tau_ind = 1;
1136 }
1137 }
1138 if (preparedInput.m_LFVmode == 0) // electron case: H->ele+tau(->mu) decays
1139 {
1140 if ((preparedInput.m_vistau1.M() < 0.05 &&
1141 preparedInput.m_vistau2.M() > 0.05) != refit) // choosing lepton from Higgs decay
1142 //When the mass calculator is rerun with refit=true the alternative lepton ordering is used
1143 {
1144 tau_tmp = preparedInput.m_vistau2;
1145 lep_tmp = preparedInput.m_vistau1;
1146 tau_type_tmp = preparedInput.m_type_visTau2;
1147 tau_ind = 2;
1148 } else {
1149 tau_tmp = preparedInput.m_vistau1;
1150 lep_tmp = preparedInput.m_vistau2;
1151 tau_type_tmp = preparedInput.m_type_visTau1;
1152 tau_ind = 1;
1153 }
1154 }
1155
1156 //------- Settings -------------------------------
1157 double Mlep = tau_tmp.M();
1158 // double dMnu_max=m_MnuScanRange-Mlep;
1159 // double Mnu_binSize=dMnu_max/NiterMnu;
1160 //-------- end of Settings
1161
1162 // double M=Mtau;
1163 double M_nu = 0.0;
1164 double MnuProb = 1.0;
1165 //---------------------------------------------
1166 for (int i3 = 0; i3 < NiterMnu; i3++) //---- loop-3: virtual neutrino mass
1167 {
1168 M_nu = Mnu_binSize * i3;
1169 if (M_nu >= (Mtau - Mlep))
1170 continue;
1171 // M=sqrt(Mtau*Mtau-M_nu*M_nu);
1172 MnuProb = Prob->MnuProbability(preparedInput, M_nu,
1173 Mnu_binSize); // Mnu probability
1174 //---------------------------------------------
1175 for (int i4 = 0; i4 < NiterMET + 1; i4++) // MET_X scan
1176 {
1177 met_smearL = METresX_binSize * i4 - N_METsigma * METresX;
1178 for (int i5 = 0; i5 < NiterMET + 1; i5++) // MET_Y scan
1179 {
1180 met_smearP = METresY_binSize * i5 - N_METsigma * METresY;
1181 if (pow(met_smearL / METresX, 2) + pow(met_smearP / METresY, 2) > pow(N_METsigma, 2))
1182 continue; // use ellipse instead of square
1183 met_smear_x = met_smearL * met_coscovphi - met_smearP * met_sincovphi;
1184 met_smear_y = met_smearL * met_sincovphi + met_smearP * met_coscovphi;
1185 metvec_tmp.SetXY(input_metX + met_smear_x, input_metY + met_smear_y);
1186
1187 solution = NuPsolutionLFV(metvec_tmp, tau_tmp, M_nu, nu_vec);
1188
1189 ++iter0;
1190
1191 if (solution < 1)
1192 continue;
1193 ++m_iter1;
1194
1195 // if fast sin cos, result to not match exactly nupsolutionv2, so skip
1196 // test
1197 // SpeedUp no nested loop to compute individual probability
1198 int ngoodsol1 = 0;
1199
1200 metprob = Prob->MetProbability(preparedInput, met_smearL, met_smearP, METresX, METresY);
1201 if (metprob <= 0)
1202 continue;
1203 for (unsigned int j1 = 0; j1 < nu_vec.size(); j1++) {
1204 if (tau_tmp.E() + nu_vec[j1].E() >= preparedInput.m_beamEnergy)
1205 continue;
1206 const double tau1_tmpp = (tau_tmp + nu_vec[j1]).P();
1207 angle1 = Angle(nu_vec[j1], tau_tmp);
1208
1209 if (angle1 < dTheta3DLimit(tau_type_tmp, 0, tau1_tmpp)) {
1210 ++m_iang1low;
1211 continue;
1212 } // lower 99% bound
1213 if (angle1 > dTheta3DLimit(tau_type_tmp, 1, tau1_tmpp)) {
1214 ++m_iang1high;
1215 continue;
1216 } // upper 99% bound
1217 double tauvecprob1j =
1218 Prob->dTheta3d_probabilityFast(preparedInput, tau_type_tmp, angle1, tau1_tmpp);
1219 if (tauvecprob1j == 0.)
1220 continue;
1221 tauprob = Prob->TauProbabilityLFV(preparedInput, tau_type_tmp, tau_tmp, nu_vec[j1]);
1222 totalProb = tauvecprob1j * metprob * MnuProb * tauprob;
1223
1224 m_tautau_tmp.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);
1225 m_tautau_tmp += tau_tmp;
1226 m_tautau_tmp += lep_tmp;
1227 m_tautau_tmp += nu_vec[j1];
1228
1229 const double mtautau = m_tautau_tmp.M();
1230
1231 m_totalProbSum += totalProb;
1232 m_mtautauSum += mtautau;
1233
1234 fit_code = 1; // at least one solution is found
1235
1236 m_fMfit_all->Fill(mtautau, totalProb);
1237 m_fMfit_allNoWeight->Fill(mtautau, 1.);
1238 //----------------- using P*fit to fill Px,y,z_tau
1239 m_fPXfit1->Fill((tau_tmp + nu_vec[j1]).Px(), totalProb);
1240 m_fPYfit1->Fill((tau_tmp + nu_vec[j1]).Py(), totalProb);
1241 m_fPZfit1->Fill((tau_tmp + nu_vec[j1]).Pz(), totalProb);
1242
1243 if (totalProb > m_prob_tmp) // fill solution with highest probability
1244 {
1245 sign_tmp = -log10(totalProb);
1246 m_prob_tmp = totalProb;
1247 m_fDitauStuffFit.Mditau_best = mtautau;
1248 m_fDitauStuffFit.Sign_best = sign_tmp;
1249 if (tau_ind == 1)
1250 m_fDitauStuffFit.nutau1 = nu_vec[j1];
1251 if (tau_ind == 2)
1252 m_fDitauStuffFit.nutau2 = nu_vec[j1];
1253 }
1254
1255 ++ngoodsol1;
1256 }
1257
1258 if (ngoodsol1 == 0)
1259 continue;
1260 m_iter2 += 1;
1261
1262 m_iter3 += 1;
1263 }
1264 }
1265 }
1266 } else if (preparedInput.m_tauTypes == TauTypes::lh) // lepton+tau case
1267 {
1268 if (preparedInput.m_fUseVerbose == 1) {
1269 Info("DiTauMassTools", "Running in lepton+tau mode");
1270 }
1271 //------- Settings -------------------------------
1272
1273 //----- Stuff below are for Winter 2012 lep-had analysis only; it has to be
1274 // replaced by a more common scheme once other channels are optimized
1275 // XYVector
1276 // mht_vec((tau_vec1+tau_vec2).Px(),(tau_vec1+tau_vec2).Py()); //
1277 // missing Ht vector for Njet25=0 events const double
1278 // mht=mht_vec.R();
1279 double input_metX = preparedInput.m_MetVec.X();
1280 double input_metY = preparedInput.m_MetVec.Y();
1281
1282 // double mht_offset=0.0;
1283 // if(InputInfo.UseHT) // use missing Ht (for 0-jet events only for
1284 // now)
1285 // {
1286 // input_metX=-mht_vec.X();
1287 // input_metY=-mht_vec.Y();
1288 // }
1289 // else // use MET (for 0-jet and 1-jet events)
1290 // {
1291 // input_metX=met_vec.X();
1292 // input_metY=met_vec.Y();
1293 // }
1294
1295 PtEtaPhiMVector tau_tmp(0.0, 0.0, 0.0, 0.0);
1296 PtEtaPhiMVector lep_tmp(0.0, 0.0, 0.0, 0.0);
1297 int tau_type_tmp;
1298 if (preparedInput.m_type_visTau1 == 8) {
1299 tau_tmp = preparedInput.m_vistau2;
1300 lep_tmp = preparedInput.m_vistau1;
1301 tau_type_tmp = preparedInput.m_type_visTau2;
1302 }
1303 if (preparedInput.m_type_visTau2 == 8) {
1304 tau_tmp = preparedInput.m_vistau1;
1305 lep_tmp = preparedInput.m_vistau2;
1306 tau_type_tmp = preparedInput.m_type_visTau1;
1307 }
1308
1309 //---------------------------------------------
1310 for (int i4 = 0; i4 < NiterMET + 1; i4++) // MET_X scan
1311 {
1312 met_smearL = METresX_binSize * i4 - N_METsigma * METresX;
1313 for (int i5 = 0; i5 < NiterMET + 1; i5++) // MET_Y scan
1314 {
1315 met_smearP = METresY_binSize * i5 - N_METsigma * METresY;
1316 if (pow(met_smearL / METresX, 2) + pow(met_smearP / METresY, 2) > pow(N_METsigma, 2))
1317 continue; // use ellipse instead of square
1318 met_smear_x = met_smearL * m_metCovPhiCos - met_smearP * m_metCovPhiSin;
1319 met_smear_y = met_smearL * m_metCovPhiSin + met_smearP * m_metCovPhiCos;
1320 metvec_tmp.SetXY(input_metX + met_smear_x, input_metY + met_smear_y);
1321
1322 solution = NuPsolutionLFV(metvec_tmp, tau_tmp, 0.0, nu_vec);
1323
1324 ++iter0;
1325
1326 if (solution < 1)
1327 continue;
1328 ++m_iter1;
1329
1330 // if fast sin cos, result to not match exactly nupsolutionv2, so skip
1331 // test
1332 // SpeedUp no nested loop to compute individual probability
1333 int ngoodsol1 = 0;
1334
1335 metprob = Prob->MetProbability(preparedInput, met_smearL, met_smearP, METresX, METresY);
1336 if (metprob <= 0)
1337 continue;
1338 for (unsigned int j1 = 0; j1 < nu_vec.size(); j1++) {
1339 if (tau_tmp.E() + nu_vec[j1].E() >= preparedInput.m_beamEnergy)
1340 continue;
1341 const double tau1_tmpp = (tau_tmp + nu_vec[j1]).P();
1342 angle1 = Angle(nu_vec[j1], tau_tmp);
1343
1344 if (angle1 < dTheta3DLimit(tau_type_tmp, 0, tau1_tmpp)) {
1345 ++m_iang1low;
1346 continue;
1347 } // lower 99% bound
1348 if (angle1 > dTheta3DLimit(tau_type_tmp, 1, tau1_tmpp)) {
1349 ++m_iang1high;
1350 continue;
1351 } // upper 99% bound
1352 double tauvecprob1j =
1353 Prob->dTheta3d_probabilityFast(preparedInput, tau_type_tmp, angle1, tau1_tmpp);
1354 if (tauvecprob1j == 0.)
1355 continue;
1356 tauprob = Prob->TauProbabilityLFV(preparedInput, tau_type_tmp, tau_tmp, nu_vec[j1]);
1357 totalProb = tauvecprob1j * metprob * tauprob;
1358
1359 m_tautau_tmp.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);
1360 m_tautau_tmp += tau_tmp;
1361 m_tautau_tmp += lep_tmp;
1362 m_tautau_tmp += nu_vec[j1];
1363
1364 const double mtautau = m_tautau_tmp.M();
1365
1366 m_totalProbSum += totalProb;
1367 m_mtautauSum += mtautau;
1368
1369 fit_code = 1; // at least one solution is found
1370
1371 m_fMfit_all->Fill(mtautau, totalProb);
1372 m_fMfit_allNoWeight->Fill(mtautau, 1.);
1374 // m_fPXfit1->Fill((tau_tmp+nu_vec[j1]).Px(),totalProb);
1375 // m_fPYfit1->Fill((tau_tmp+nu_vec[j1]).Py(),totalProb);
1376 // m_fPZfit1->Fill((tau_tmp+nu_vec[j1]).Pz(),totalProb);
1377
1378 if (totalProb > m_prob_tmp) // fill solution with highest probability
1379 {
1380 sign_tmp = -log10(totalProb);
1381 m_prob_tmp = totalProb;
1382 m_fDitauStuffFit.Mditau_best = mtautau;
1383 m_fDitauStuffFit.Sign_best = sign_tmp;
1384 if (preparedInput.m_type_visTau1 == 8) {
1385 m_fDitauStuffFit.vistau1 = lep_tmp;
1386 m_fDitauStuffFit.vistau2 = tau_tmp;
1387 m_fDitauStuffFit.nutau2 = nu_vec[j1];
1388 } else if (preparedInput.m_type_visTau2 == 8) {
1389 m_fDitauStuffFit.vistau2 = lep_tmp;
1390 m_fDitauStuffFit.vistau1 = tau_tmp;
1391 m_fDitauStuffFit.nutau1 = nu_vec[j1];
1392 }
1393 }
1394
1395 ++ngoodsol1;
1396 }
1397
1398 if (ngoodsol1 == 0)
1399 continue;
1400 m_iter2 += 1;
1401
1402 m_iter3 += 1;
1403 }
1404 }
1405 } else {
1406 Info("DiTauMassTools", "Running in an unknown mode?!?!");
1407 }
1408
1409 OutputInfo.m_NTrials = iter0;
1410 OutputInfo.m_NSuccesses = m_iter3;
1411
1412 if (preparedInput.m_fUseVerbose == 1) {
1413 Info("DiTauMassTools", "%s",
1414 ("SpeedUp niters=" + std::to_string(iter0) + " " + std::to_string(m_iter1) + " " +
1415 std::to_string(m_iter2) + " " + std::to_string(m_iter3) + "skip:" + std::to_string(m_iang1low) +
1416 " " + std::to_string(m_iang1high))
1417 .c_str());
1418 }
1419
1420 if (m_fMfit_all->GetEntries() > 0 && m_iter3 > 0) {
1421#ifdef SMOOTH
1422 m_fMfit_all->Smooth();
1423 m_fMfit_allNoWeight->Smooth();
1424 m_fPXfit1->Smooth();
1425 m_fPYfit1->Smooth();
1426 m_fPZfit1->Smooth();
1427#endif
1428
1429 // default max finding method defined in MissingMassCalculator.h
1430 // note that window defined in terms of number of bin, so depend on binning
1431 std::vector<double> histInfo(HistInfo::MAXHISTINFO);
1432 m_fDitauStuffHisto.Mditau_best = maxFromHist(m_fMfit_all, histInfo);
1433 double prob_hist = histInfo.at(HistInfo::PROB);
1434
1435 if (prob_hist != 0.0)
1436 m_fDitauStuffHisto.Sign_best = -log10(std::abs(prob_hist));
1437 else {
1438 // this mean the histogram is empty.
1439 // possible but very rare if all entries outside histogram range
1440 // fall back to maximum
1441 m_fDitauStuffHisto.Sign_best = -999.;
1442 m_fDitauStuffHisto.Mditau_best = m_fDitauStuffFit.Mditau_best;
1443 }
1444
1445 if (m_fDitauStuffHisto.Mditau_best > 0.0)
1446 m_fDitauStuffHisto.RMSoverMPV = m_fMfit_all->GetRMS() / m_fDitauStuffHisto.Mditau_best;
1447 std::vector<double> histInfoOther(HistInfo::MAXHISTINFO);
1448 //---- getting Nu1
1449 double Px1 = maxFromHist(m_fPXfit1, histInfoOther);
1450 double Py1 = maxFromHist(m_fPYfit1, histInfoOther);
1451 double Pz1 = maxFromHist(m_fPZfit1, histInfoOther);
1452 //---- setting 4-vecs
1453 PxPyPzMVector nu1_tmp(0.0, 0.0, 0.0, 0.0);
1454 PxPyPzMVector nu2_tmp(0.0, 0.0, 0.0, 0.0);
1455 if (preparedInput.m_type_visTau1 == 8) {
1456 nu1_tmp = preparedInput.m_vistau1;
1457 nu2_tmp.SetCoordinates(Px1, Py1, Pz1, ParticleConstants::tauMassInMeV / GEV);
1458 }
1459 if (preparedInput.m_type_visTau2 == 8) {
1460 nu2_tmp = preparedInput.m_vistau2;
1461 nu1_tmp.SetCoordinates(Px1, Py1, Pz1, ParticleConstants::tauMassInMeV / GEV);
1462 }
1463 m_fDitauStuffHisto.nutau1 = nu1_tmp - preparedInput.m_vistau1;
1464 m_fDitauStuffHisto.nutau2 = nu2_tmp - preparedInput.m_vistau2;
1465 }
1466 if (m_lfvLeplepRefit && fit_code==0 && !refit) {
1467 fit_code = DitauMassCalculatorV9lfv(true);
1468 return fit_code;
1469 }
1470
1471
1472
1473 if (preparedInput.m_fUseVerbose == 1) {
1474 if (fit_code == 0) {
1475 Info(
1476 "DiTauMassTools", "%s",
1477 ("!!!----> Warning-3 in MissingMassCalculator::DitauMassCalculatorV9lfv() : fit status=" +
1478 std::to_string(fit_code))
1479 .c_str());
1480 Info("DiTauMassTools", "....... No solution is found. Printing input info .......");
1481
1482 Info("DiTauMassTools", "%s", (" vis Tau-1: Pt="+std::to_string(preparedInput.m_vistau1.Pt())
1483 +" M="+std::to_string(preparedInput.m_vistau1.M())+" eta="+std::to_string(preparedInput.m_vistau1.Eta())
1484 +" phi="+std::to_string(preparedInput.m_vistau1.Phi())
1485 +" type="+std::to_string(preparedInput.m_type_visTau1)).c_str());
1486 Info("DiTauMassTools", "%s", (" vis Tau-2: Pt="+std::to_string(preparedInput.m_vistau2.Pt())
1487 +" M="+std::to_string(preparedInput.m_vistau2.M())+" eta="+std::to_string(preparedInput.m_vistau2.Eta())
1488 +" phi="+std::to_string(preparedInput.m_vistau2.Phi())
1489 +" type="+std::to_string(preparedInput.m_type_visTau2)).c_str());
1490 Info("DiTauMassTools", "%s", (" MET="+std::to_string(preparedInput.m_MetVec.R())+" Met_X="+std::to_string(preparedInput.m_MetVec.X())
1491 +" Met_Y="+std::to_string(preparedInput.m_MetVec.Y())).c_str());
1492 Info("DiTauMassTools", " ---------------------------------------------------------- ");
1493 }
1494 }
1495 return fit_code;
1496}
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 XYVector deltamet_vec;
846
847 // initialize a spacewalker, which walks the parameter space according to some
848 // algorithm
850
851 while (SpaceWalkerWalk()) {
852 bool paramInsideRange = false;
853 m_nsol = 0;
854
855 paramInsideRange = checkAllParamInRange();
856
857 // FIXME if no tau scanning, or symmetric matrices, rotatin is made twice
858 // which is inefficient
859 const double deltaMetx = m_MEtL * m_metCovPhiCos - m_MEtP * m_metCovPhiSin;
860 const double deltaMety = m_MEtL * m_metCovPhiSin + m_MEtP * m_metCovPhiCos;
861
862 // deltaMetVec.Set(met_smear_x,met_smear_y);
863 preparedInput.m_metVec.SetXY(preparedInput.m_inputMEtX + deltaMetx,
864 preparedInput.m_inputMEtY + deltaMety);
865
866 // save in global variable for speed sake
867 preparedInput.m_MEtX = preparedInput.m_metVec.X();
868 preparedInput.m_MEtY = preparedInput.m_metVec.Y();
869 preparedInput.m_MEtT = preparedInput.m_metVec.R();
870
871 if (paramInsideRange)
873
874 // DR for markov chain need to enter handleSolution also when zero solutions
876 // be careful that with markov, current solution is from now on stored in
877 // XYZOldSolVec
878
879 if (m_nsol <= 0)
880 continue;
881
882 // for markov, nsuccess more difficult to define. Decide this is the number
883 // of independent point accepted (hence without weight)
884 nsuccesses = m_markovNAccept;
886
887 m_iter1 += m_nsol;
888 fit_code = 1;
889
890 } // while loop
891
892 OutputInfo.m_NTrials = m_iter0;
893 OutputInfo.m_NSuccesses = nsuccesses;
894
895 if (nsuccesses > 0) {
896 OutputInfo.m_AveSolRMS /= nsuccesses;
897 } else {
898 OutputInfo.m_AveSolRMS = -1.;
899 }
900
901 double Px1, Py1, Pz1;
902 double Px2, Py2, Pz2;
903 if (nsuccesses > 0) {
904
905 // note that smoothing can slightly change the integral of the histogram
906
907#ifdef SMOOTH
908 m_fMfit_all->Smooth();
909 m_fMfit_allNoWeight->Smooth();
910 m_fPXfit1->Smooth();
911 m_fPYfit1->Smooth();
912 m_fPZfit1->Smooth();
913 m_fPXfit2->Smooth();
914 m_fPYfit2->Smooth();
915 m_fPZfit2->Smooth();
916#endif
917
918 // default max finding method defined in MissingMassCalculator.h
919 // note that window defined in terms of number of bin, so depend on binning
920 std::vector<double> histInfo(HistInfo::MAXHISTINFO);
921 m_fDitauStuffHisto.Mditau_best = maxFromHist(m_fMfit_all, histInfo);
922 double prob_hist = histInfo.at(HistInfo::PROB);
923
924 if (prob_hist != 0.0)
925 m_fDitauStuffHisto.Sign_best = -log10(std::abs(prob_hist));
926 else {
927 // this mean the histogram is empty.
928 // possible but very rare if all entries outside histogram range
929 // fall back to maximum
930 m_fDitauStuffHisto.Sign_best = -999.;
931 m_fDitauStuffHisto.Mditau_best = m_fDitauStuffFit.Mditau_best;
932 }
933
934 if (m_fDitauStuffHisto.Mditau_best > 0.0)
935 m_fDitauStuffHisto.RMSoverMPV = m_fMfit_all->GetRMS() / m_fDitauStuffHisto.Mditau_best;
936 std::vector<double> histInfoOther(HistInfo::MAXHISTINFO);
937 //---- getting full tau1 momentum
938 Px1 = maxFromHist(m_fPXfit1, histInfoOther);
939 Py1 = maxFromHist(m_fPYfit1, histInfoOther);
940 Pz1 = maxFromHist(m_fPZfit1, histInfoOther);
941
942 //---- getting full tau2 momentum
943 Px2 = maxFromHist(m_fPXfit2, histInfoOther);
944 Py2 = maxFromHist(m_fPYfit2, histInfoOther);
945 Pz2 = maxFromHist(m_fPZfit2, histInfoOther);
946
947 //---- setting 4-vecs
948 PxPyPzMVector fulltau1, fulltau2;
949 fulltau1.SetCoordinates(Px1, Py1, Pz1, ParticleConstants::tauMassInMeV / GEV);
950 fulltau2.SetCoordinates(Px2, Py2, Pz2, ParticleConstants::tauMassInMeV / GEV);
951 // PtEtaPhiMVector fulltau1(_fulltau1.Pt(), _fulltau1.Eta(), _fulltau1.Phi(), _fulltau1.M());
952 //PtEtaPhiMVector fulltau2(_fulltau2.Pt(), _fulltau2.Eta(), _fulltau2.Phi(), _fulltau2.M());
953
954 if (fulltau1.P() < preparedInput.m_vistau1.P())
955 fulltau1 = 1.01 * preparedInput.m_vistau1; // protection against cases when fitted tau
956 // momentum is smaller than visible tau momentum
957 if (fulltau2.P() < preparedInput.m_vistau2.P())
958 fulltau2 = 1.01 * preparedInput.m_vistau2; // protection against cases when fitted tau
959 // momentum is smaller than visible tau momentum
960 m_fDitauStuffHisto.vistau1 = preparedInput.m_vistau1; // FIXME should also be fitted if tau scan
961 m_fDitauStuffHisto.vistau2 = preparedInput.m_vistau2;
962 m_fDitauStuffHisto.nutau1 = fulltau1 - preparedInput.m_vistau1; // these are the original tau vis
963 m_fDitauStuffHisto.nutau2 =
964 fulltau2 - preparedInput.m_vistau2; // FIXME neutrino mass not necessarily zero
965 }
966
967 // Note that for v9walk, points outside the METx MEty disk are counted, while
968 // this was not the case for v9
969 if (preparedInput.m_fUseVerbose == 1) {
970 Info("DiTauMassTools", "Scanning ");
971 Info("DiTauMassTools", " Markov ");
972 Info("DiTauMassTools", "%s",
973 (" V9W niters=" + std::to_string(m_iter0) + " " + std::to_string(m_iter1)).c_str());
974 Info("DiTauMassTools", "%s", (" nFullScan " + std::to_string(m_markovNFullScan)).c_str());
975 Info("DiTauMassTools", "%s", (" nRejectNoSol " + std::to_string(m_markovNRejectNoSol)).c_str());
976 Info("DiTauMassTools", "%s", (" nRejectMetro " + std::to_string(m_markovNRejectMetropolis)).c_str());
977 Info("DiTauMassTools", "%s", (" nAccept " + std::to_string(m_markovNAccept)).c_str());
978 Info("DiTauMassTools", "%s",
979 (" probsum " + std::to_string(m_totalProbSum) + " msum " + std::to_string(m_mtautauSum))
980 .c_str());
981 }
982
983 if (preparedInput.m_fUseVerbose == 1) {
984 if (fit_code == 0) {
985 Info("DiTauMassTools", "%s", ("!!!----> Warning-3 in "
986 "MissingMassCalculator::DitauMassCalculatorV9Walk() : fit status=" +
987 std::to_string(fit_code))
988 .c_str());
989 Info("DiTauMassTools", "%s", "....... No solution is found. Printing input info .......");
990
991 Info("DiTauMassTools", "%s", (" vis Tau-1: Pt=" + std::to_string(preparedInput.m_vistau1.Pt()) +
992 " M=" + std::to_string(preparedInput.m_vistau1.M()) +
993 " eta=" + std::to_string(preparedInput.m_vistau1.Eta()) +
994 " phi=" + std::to_string(preparedInput.m_vistau1.Phi()) +
995 " type=" + std::to_string(preparedInput.m_type_visTau1))
996 .c_str());
997 Info("DiTauMassTools", "%s", (" vis Tau-2: Pt=" + std::to_string(preparedInput.m_vistau2.Pt()) +
998 " M=" + std::to_string(preparedInput.m_vistau2.M()) +
999 " eta=" + std::to_string(preparedInput.m_vistau2.Eta()) +
1000 " phi=" + std::to_string(preparedInput.m_vistau2.Phi()) +
1001 " type=" + std::to_string(preparedInput.m_type_visTau2))
1002 .c_str());
1003 Info("DiTauMassTools", "%s", (" MET=" + std::to_string(preparedInput.m_MetVec.R()) +
1004 " Met_X=" + std::to_string(preparedInput.m_MetVec.X()) +
1005 " Met_Y=" + std::to_string(preparedInput.m_MetVec.Y()))
1006 .c_str());
1007 Info("DiTauMassTools", " ---------------------------------------------------------- ");
1008 }
1009 }
1010
1011 return fit_code;
1012}
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 2722 of file MissingMassCalculator.cxx.

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

◆ FinalizeSettings()

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

Definition at line 2826 of file MissingMassCalculator.cxx.

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

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

◆ 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 1499 of file MissingMassCalculator.cxx.

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

◆ 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 1519 of file MissingMassCalculator.cxx.

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

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

◆ 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 1795 of file MissingMassCalculator.cxx.

1797 {
1798 // bool debug=true;
1799
1800 int nsol1;
1801 int nsol2;
1802
1803 const int solution = NuPsolutionV3(M_nu1, M_nu2, phi1, phi2, nsol1, nsol2);
1804
1805 if (solution != 1)
1806 return -4;
1807 // refineSolutions ( M_nu1,M_nu2,
1808 // met_smearL,met_smearP,metvec_tmp.R(),
1809 // nsol1, nsol2,m_Mvis,m_Meff);
1810 refineSolutions(M_nu1, M_nu2, nsol1, nsol2, m_Mvis, m_Meff);
1811
1812 if (m_nsol <= 0)
1813 return 0;
1814
1815 // success
1816
1817 return m_nsol; // for backward compatibility
1818}
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 1821 of file MissingMassCalculator.cxx.

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

◆ SaveLlhHisto()

void MissingMassCalculator::SaveLlhHisto ( const bool val)

Definition at line 3095 of file MissingMassCalculator.cxx.

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

◆ 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 3130 of file MissingMassCalculator.cxx.

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

◆ SpaceWalkerInit()

void MissingMassCalculator::SpaceWalkerInit ( )
inlineprotected

Definition at line 2316 of file MissingMassCalculator.cxx.

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

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

◆ 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 1992 of file MissingMassCalculator.cxx.

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

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: