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)
bool MassCollinear (const xAOD::IParticle *p0, const xAOD::IParticle *p1, const xAOD::MissingET *met, const bool kMMCsynchronize, double &mass, double &xp1, double &xp2)
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 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 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
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
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 {}
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_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 {}
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_dRmax_tau {}
double m_MnuScanRange {}

Detailed Description

Definition at line 47 of file MissingMassCalculator.h.

Constructor & Destructor Documentation

◆ ~MissingMassCalculator()

MissingMassCalculator::~MissingMassCalculator ( )

Definition at line 188 of file MissingMassCalculator.cxx.

188{ delete Prob; }

◆ MissingMassCalculator() [1/2]

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

Definition at line 49 of file MissingMassCalculator.cxx.

51 : m_randomGen(), Prob(new MissingMassProb(aset, paramFilePath)) {
53 preparedInput.m_fUseVerbose = 0;
54 preparedInput.m_beamEnergy = 6500.0; // for now LHC default is sqrt(S)=7 TeV
55 m_niter_fit1 = 20;
56 m_niter_fit2 = 30;
57 m_niter_fit3 = 10;
58 m_NsucStop = -1;
59 m_NiterRandom = -1; // if the user does not set it to positive value,will be set
60 // in SpaceWalkerInit
61 m_niterRandomLocal = -1; // niterandom which is really used
62 // to be used with RMSSTOP NiterRandom=10000000; // number of random
63 // iterations for lh. Multiplied by 10 for ll, divided by 10 for hh (to be
64 // optimised)
65 // RMSStop=200;// Stop criteria depending of rms of histogram
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
82 Prob->SetUseTauProbability(true); // TauProbability is ON by default DRMERGE comment out for now
83 Prob->SetUseMnuProbability(false); // MnuProbability is OFF by default
84 Prob->SetUseDphiLL(false); // added by Tomas Davidek for lep-lep
85 preparedInput.m_METresSyst = 0; // no MET resolution systematics by default (+/-1: up/down 1 sigma)
86 preparedInput.m_dataType = 1; // set to "data" by default
87 preparedInput.m_fUseTailCleanup = 1; // cleanup by default for lep-had Moriond 2012 analysis
88 preparedInput.m_fUseDefaults = 0; // use pre-set defaults for various configurations; if set it to 0
89 // if need to study various options
90 m_fUseEfficiencyRecovery = 0; // no re-fit by default
95
96 preparedInput.m_METScanScheme = 1; // MET-scan scheme: 0- use JER; 1- use simple sumEt & missingHt
97 // for Njet=0 events in (lep-had winter 2012)
98 // MnuScanRange=ParticleConstants::tauMassInMeV / GEV; // range of M(nunu) scan
99 m_MnuScanRange = 1.5; // better value (sacha)
100 preparedInput.m_LFVmode = -1; // by default consider case of H->mu+tau(->ele)
101 preparedInput.ClearInput();
102
103 m_debugThisIteration = false;
104 m_lfvLeplepRefit = true;
105 m_SaveLlhHisto = false;
106
107 m_nsolmax = 4;
109
110 m_nuvecsol1.resize(m_nsolmax);
111 m_nuvecsol2.resize(m_nsolmax);
112 m_tauvecsol1.resize(m_nsolmax);
113 m_tauvecsol2.resize(m_nsolmax);
114 m_tauvecprob1.resize(m_nsolmax);
115 m_tauvecprob2.resize(m_nsolmax);
116
117 m_nsol = 0;
122
123 m_nsolOld = 0;
128
129 float hEmax = 3000.0; // maximum energy (GeV)
130 // number of bins
131 int hNbins = 1500; // original 2500 for mass, 10000 for P
132 // choice of hNbins also related to size of window for fitting (see
133 // maxFromHist)
134
135 //--- define histograms for histogram method
136 //--- upper limits need to be revisied in the future!!! It may be not enough
137 // for some analyses
138
139 m_fMfit_all = std::make_shared<TH1F>("MMC_h1", "M", hNbins, 0.0,
140 hEmax); // all solutions
141 m_fMfit_all->Sumw2(); // allow proper error bin calculation. Slightly slower but
142 // completely negligible
143
144 // histogram without weight. useful for debugging. negligibly slow until now
146 std::make_shared<TH1F>("MMC_h1NoW", "M no weight", hNbins, 0.0, hEmax); // all solutions
147
148 m_fPXfit1 = std::make_shared<TH1F>("MMC_h2", "Px1", 4 * hNbins, -hEmax,
149 hEmax); // Px for tau1
150 m_fPYfit1 = std::make_shared<TH1F>("MMC_h3", "Py1", 4 * hNbins, -hEmax,
151 hEmax); // Py for tau1
152 m_fPZfit1 = std::make_shared<TH1F>("MMC_h4", "Pz1", 4 * hNbins, -hEmax,
153 hEmax); // Pz for tau1
154 m_fPXfit2 = std::make_shared<TH1F>("MMC_h5", "Px2", 4 * hNbins, -hEmax,
155 hEmax); // Px for tau2
156 m_fPYfit2 = std::make_shared<TH1F>("MMC_h6", "Py2", 4 * hNbins, -hEmax,
157 hEmax); // Py for tau2
158 m_fPZfit2 = std::make_shared<TH1F>("MMC_h7", "Pz2", 4 * hNbins, -hEmax,
159 hEmax); // Pz for tau2
160
161 m_fMfit_all->SetDirectory(0);
162
163 m_fMfit_allNoWeight->SetDirectory(0);
164 m_fPXfit1->SetDirectory(0);
165 m_fPYfit1->SetDirectory(0);
166 m_fPZfit1->SetDirectory(0);
167 m_fPXfit2->SetDirectory(0);
168 m_fPYfit2->SetDirectory(0);
169 m_fPZfit2->SetDirectory(0);
170
171 // max hist fitting function
172 m_fFitting =
173 new TF1("MMC_maxFitting", this, &MissingMassCalculator::maxFitting, 0., hEmax, 3);
174 // Sets initial parameter names
175 m_fFitting->SetParNames("Max", "Mean", "InvWidth2");
176
177 if (preparedInput.m_fUseVerbose == 1) {
178 gDirectory->pwd();
179 gDirectory->ls();
180 }
181
182 if (preparedInput.m_fUseVerbose == 1) {
183 gDirectory->pwd();
184 gDirectory->ls();
185 }
186}
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 2657 of file MissingMassCalculator.cxx.

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

◆ checkMEtInRange()

bool MissingMassCalculator::checkMEtInRange ( )
inlineprotected

Definition at line 2698 of file MissingMassCalculator.cxx.

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

◆ CheckSolutions()

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

◆ ClearDitauStuff()

void MissingMassCalculator::ClearDitauStuff ( DitauStuff & fStuff)
private

Definition at line 285 of file MissingMassCalculator.cxx.

285 {
286 fStuff.Mditau_best = 0.0;
287 fStuff.Sign_best = 1.0E6;
288 fStuff.nutau1 = PtEtaPhiMVector(0., 0., 0., 0.);
289 fStuff.nutau2 = PtEtaPhiMVector(0., 0., 0., 0.);
290 fStuff.vistau1 = PtEtaPhiMVector(0., 0., 0., 0.);
291 fStuff.vistau2 = PtEtaPhiMVector(0., 0., 0., 0.);
292 fStuff.RMSoverMPV = 0.0;
293
294 return;
295}

◆ DitauMassCalculatorV9lfv()

int MissingMassCalculator::DitauMassCalculatorV9lfv ( bool refit)
inlineprotected

Definition at line 1006 of file MissingMassCalculator.cxx.

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

◆ DitauMassCalculatorV9walk()

int MissingMassCalculator::DitauMassCalculatorV9walk ( )
inlineprotected

Definition at line 786 of file MissingMassCalculator.cxx.

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

◆ DoOutputInfo()

void MissingMassCalculator::DoOutputInfo ( )
private

Definition at line 300 of file MissingMassCalculator.cxx.

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

◆ dTheta3DLimit()

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

Definition at line 2714 of file MissingMassCalculator.cxx.

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

◆ FinalizeSettings()

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

Definition at line 2819 of file MissingMassCalculator.cxx.

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

379{ return m_markovNAccept; }

◆ GetMarkovNFullscan()

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

Definition at line 380 of file MissingMassCalculator.h.

380{ return m_markovNFullScan;}

◆ GetMarkovNRejectMetropolis()

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

Definition at line 378 of file MissingMassCalculator.h.

◆ GetMarkovNRejectNoSol()

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

Definition at line 377 of file MissingMassCalculator.h.

377{ return m_markovNRejectNoSol;}

◆ GetMeanbinStop()

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

Definition at line 373 of file MissingMassCalculator.h.

373{ 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 366 of file MissingMassCalculator.h.

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

◆ GetNiterFit2()

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

Definition at line 367 of file MissingMassCalculator.h.

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

◆ GetNiterFit3()

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

Definition at line 368 of file MissingMassCalculator.h.

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

◆ GetNiterRandom()

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

Definition at line 369 of file MissingMassCalculator.h.

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

◆ GetNMetroReject()

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

Definition at line 403 of file MissingMassCalculator.h.

◆ GetNNoSol()

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

Definition at line 402 of file MissingMassCalculator.h.

402{return m_markovNRejectNoSol;}

◆ GetNSol()

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

Definition at line 404 of file MissingMassCalculator.h.

404{return m_markovNAccept;}

◆ GetNsucStop()

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

Definition at line 371 of file MissingMassCalculator.h.

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

◆ GetProposalTryMEt()

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

Definition at line 381 of file MissingMassCalculator.h.

381{return m_proposalTryMEt;}

◆ GetProposalTryMnu()

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

Definition at line 383 of file MissingMassCalculator.h.

383{return m_ProposalTryMnu;}

◆ GetProposalTryPhi()

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

Definition at line 382 of file MissingMassCalculator.h.

382{return m_ProposalTryPhi;}

◆ GetRMSStop()

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

Definition at line 372 of file MissingMassCalculator.h.

372{ return m_RMSStop; }

◆ GetRndmSeedAltering()

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

Definition at line 374 of file MissingMassCalculator.h.

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

◆ GetUseEfficiencyRecovery()

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

Definition at line 364 of file MissingMassCalculator.h.

364{ return m_fUseEfficiencyRecovery; }

◆ handleSolutions()

void MissingMassCalculator::handleSolutions ( )
inlineprotected

Definition at line 2039 of file MissingMassCalculator.cxx.

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

◆ MassCollinear()

bool MissingMassCalculator::MassCollinear ( const xAOD::IParticle * p0,
const xAOD::IParticle * p1,
const xAOD::MissingET * met,
const bool kMMCsynchronize,
double & mass,
double & xp1,
double & xp2 )

redefine tau vectors if necessary - MMC sychronization

Definition at line 3176 of file MissingMassCalculator.cxx.

3179 { // result
3180
3181 TLorentzVector k1 = p0->p4();
3182 TLorentzVector k2 = p1->p4();
3183
3185 if (kMMCsynchronize) {
3186 if (p0->type() == xAOD::Type::Tau) {
3187 const xAOD::TauJet *tau0 = static_cast<const xAOD::TauJet *>(p0);
3188 k1.SetPtEtaPhiM(k1.Pt(), k1.Eta(), k1.Phi(),
3189 tau0->nTracks() < 3 ? 800. : 1200.); // MeV
3190 }
3191
3192 if (p1->type() == xAOD::Type::Tau) {
3193 const xAOD::TauJet *tau1 = static_cast<const xAOD::TauJet *>(p1);
3194 k2.SetPtEtaPhiM(k2.Pt(), k2.Eta(), k2.Phi(),
3195 tau1->nTracks() < 3 ? 800. : 1200.); // MeV
3196 }
3197 }
3198
3199 TMatrixD K(2, 2);
3200 K(0, 0) = k1.Px();
3201 K(0, 1) = k2.Px();
3202 K(1, 0) = k1.Py();
3203 K(1, 1) = k2.Py();
3204
3205 if (K.Determinant() == 0)
3206 return false;
3207
3208 TMatrixD M(2, 1);
3209 M(0, 0) = met->mpx();
3210 M(1, 0) = met->mpy();
3211
3212 TMatrixD Kinv = K.Invert();
3213
3214 TMatrixD X(2, 1);
3215 X = Kinv * M;
3216
3217 double X1 = X(0, 0);
3218 double X2 = X(1, 0);
3219 double x1 = 1. / (1. + X1);
3220 double x2 = 1. / (1. + X2);
3221
3222 TLorentzVector par1 = k1 * (1 / x1);
3223 TLorentzVector par2 = k2 * (1 / x2);
3224
3225 double m = (par1 + par2).M();
3226
3227 // return to caller
3228 mass = m;
3229
3230 if (k1.Pt() > k2.Pt()) {
3231 xp1 = x1;
3232 xp2 = x2;
3233 } else {
3234 xp1 = x2;
3235 xp2 = x1;
3236 }
3237
3238 return true;
3239}
static Double_t tau0
size_t nTracks(TauJetParameters::TauTrackFlag flag=TauJetParameters::TauTrackFlag::classifiedCharged) const
@ Tau
The object is a tau (jet).
Definition ObjectType.h:49
TauJet_v3 TauJet
Definition of the current "tau version".
Definition TauJet.h:17

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

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

◆ 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 412 of file MissingMassCalculator.h.

412 {
413 return maxFromHist(theHist.get(), histInfo, maxHistStrategy, winHalfWidth, debug);
414 }
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 1511 of file MissingMassCalculator.cxx.

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

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

◆ NuPsolutionV3()

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

Definition at line 544 of file MissingMassCalculator.cxx.

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

2612 {
2613
2614 // copy tau 4 vect. If tau E scanning, these vectors will be modified
2615 m_tauVec1 = preparedInput.m_vistau1;
2616 m_tauVec2 = preparedInput.m_vistau2;
2617
2618 const XYVector &metVec = preparedInput.m_MetVec;
2619
2620 bool same = true;
2635
2637 same = updateDouble(std::pow(m_mTau, 2), m_mTau2) && same;
2641
2642 PtEtaPhiMVector Met4vec;
2643 Met4vec.SetPxPyPzE(preparedInput.m_MetVec.X(), preparedInput.m_MetVec.Y(), 0.0,
2644 preparedInput.m_MetVec.R());
2645 same = updateDouble((m_tauVec1 + m_tauVec2 + Met4vec).M(), m_Meff) && same;
2646
2647 same = updateDouble(preparedInput.m_HtOffset, preparedInput.m_htOffset) && same;
2648 // note that if useHT met_vec is actually -HT
2649 same = updateDouble(metVec.X(), preparedInput.m_inputMEtX) && same;
2650 same = updateDouble(metVec.Y(), preparedInput.m_inputMEtY) && same;
2651 same = updateDouble(metVec.R(), preparedInput.m_inputMEtT) && same;
2652
2653 return same;
2654}

◆ PrintOtherInput()

void MissingMassCalculator::PrintOtherInput ( )
private

Definition at line 397 of file MissingMassCalculator.cxx.

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

◆ PrintResults()

void MissingMassCalculator::PrintResults ( )
private

Definition at line 439 of file MissingMassCalculator.cxx.

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

◆ probCalculatorV9fast()

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

Definition at line 1787 of file MissingMassCalculator.cxx.

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

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

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

◆ SaveLlhHisto()

void MissingMassCalculator::SaveLlhHisto ( const bool val)

Definition at line 3088 of file MissingMassCalculator.cxx.

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

◆ SetBeamEnergy()

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

Definition at line 394 of file MissingMassCalculator.h.

394{ m_beamEnergy=val; }

◆ SetEventNumber()

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

Definition at line 355 of file MissingMassCalculator.h.

◆ SetFloatStoppingCheckFreq()

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

Definition at line 392 of file MissingMassCalculator.h.

◆ SetFloatStoppingComp()

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

Definition at line 393 of file MissingMassCalculator.h.

◆ SetFloatStoppingMinIter()

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

Definition at line 391 of file MissingMassCalculator.h.

◆ SetLFVLeplepRefit()

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

Definition at line 395 of file MissingMassCalculator.h.

◆ SetMeanbinStop()

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

Definition at line 353 of file MissingMassCalculator.h.

◆ SetMnuScanRange()

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

Definition at line 357 of file MissingMassCalculator.h.

◆ SetNiterFit1()

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

Definition at line 347 of file MissingMassCalculator.h.

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

◆ SetNiterFit2()

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

Definition at line 348 of file MissingMassCalculator.h.

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

◆ SetNiterFit3()

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

Definition at line 349 of file MissingMassCalculator.h.

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

◆ SetNiterRandom()

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

Definition at line 350 of file MissingMassCalculator.h.

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

◆ SetNsigmaMETscan()

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

Definition at line 388 of file MissingMassCalculator.h.

388{ 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 387 of file MissingMassCalculator.h.

387{ 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 386 of file MissingMassCalculator.h.

386{ 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 385 of file MissingMassCalculator.h.

385{ 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 351 of file MissingMassCalculator.h.

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

◆ SetProposalTryMEt()

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

Definition at line 359 of file MissingMassCalculator.h.

◆ SetProposalTryMnu()

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

Definition at line 361 of file MissingMassCalculator.h.

◆ SetProposalTryPhi()

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

Definition at line 360 of file MissingMassCalculator.h.

◆ SetRMSStop()

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

Definition at line 352 of file MissingMassCalculator.h.

352{ m_RMSStop=val;}

◆ SetRndmSeedAltering()

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

Definition at line 354 of file MissingMassCalculator.h.

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

◆ SetUseEfficiencyRecovery()

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

Definition at line 363 of file MissingMassCalculator.h.

◆ SetUseFloatStopping()

void MissingMassCalculator::SetUseFloatStopping ( const bool val)

Definition at line 3123 of file MissingMassCalculator.cxx.

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

◆ SpaceWalkerInit()

void MissingMassCalculator::SpaceWalkerInit ( )
inlineprotected

Definition at line 2308 of file MissingMassCalculator.cxx.

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

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

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

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

Member Data Documentation

◆ m_beamEnergy

double DiTauMassTools::MissingMassCalculator::m_beamEnergy {}
private

Definition at line 98 of file MissingMassCalculator.h.

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

◆ m_cosPhi1

double DiTauMassTools::MissingMassCalculator::m_cosPhi1 {}
private

Definition at line 172 of file MissingMassCalculator.h.

172{}, m_cosPhi2{}, m_sinPhi1{}, m_sinPhi2{};

◆ m_cosPhi2

double DiTauMassTools::MissingMassCalculator::m_cosPhi2 {}
private

Definition at line 172 of file MissingMassCalculator.h.

172{}, m_cosPhi2{}, m_sinPhi1{}, m_sinPhi2{};

◆ m_debugThisIteration

bool DiTauMassTools::MissingMassCalculator::m_debugThisIteration
private

Definition at line 94 of file MissingMassCalculator.h.

◆ m_dRmax_tau

double DiTauMassTools::MissingMassCalculator::m_dRmax_tau {}
private

Definition at line 255 of file MissingMassCalculator.h.

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

◆ m_E2v1

double DiTauMassTools::MissingMassCalculator::m_E2v1 {}
private

Definition at line 188 of file MissingMassCalculator.h.

188{};

◆ m_E2v2

double DiTauMassTools::MissingMassCalculator::m_E2v2 {}
private

Definition at line 189 of file MissingMassCalculator.h.

189{};

◆ m_ET2v1

double DiTauMassTools::MissingMassCalculator::m_ET2v1 {}
private

Definition at line 186 of file MissingMassCalculator.h.

186{};

◆ m_ET2v2

double DiTauMassTools::MissingMassCalculator::m_ET2v2 {}
private

Definition at line 187 of file MissingMassCalculator.h.

187{};

◆ m_eTau1

double DiTauMassTools::MissingMassCalculator::m_eTau1 {}
private

Definition at line 141 of file MissingMassCalculator.h.

141{}, m_eTau2{};

◆ m_eTau10

double DiTauMassTools::MissingMassCalculator::m_eTau10 {}
private

Definition at line 142 of file MissingMassCalculator.h.

142{}, m_eTau20{};

◆ m_eTau1Max

double DiTauMassTools::MissingMassCalculator::m_eTau1Max {}
private

◆ m_eTau1Min

double DiTauMassTools::MissingMassCalculator::m_eTau1Min {}
private

Definition at line 152 of file MissingMassCalculator.h.

◆ m_eTau1Proposal

double DiTauMassTools::MissingMassCalculator::m_eTau1Proposal {}
private

◆ m_eTau1Range

double DiTauMassTools::MissingMassCalculator::m_eTau1Range {}
private

Definition at line 152 of file MissingMassCalculator.h.

◆ m_eTau2

double DiTauMassTools::MissingMassCalculator::m_eTau2 {}
private

Definition at line 141 of file MissingMassCalculator.h.

141{}, m_eTau2{};

◆ m_eTau20

double DiTauMassTools::MissingMassCalculator::m_eTau20 {}
private

Definition at line 142 of file MissingMassCalculator.h.

142{}, m_eTau20{};

◆ m_eTau2Max

double DiTauMassTools::MissingMassCalculator::m_eTau2Max {}
private

◆ m_eTau2Min

double DiTauMassTools::MissingMassCalculator::m_eTau2Min {}
private

Definition at line 153 of file MissingMassCalculator.h.

◆ m_eTau2Proposal

double DiTauMassTools::MissingMassCalculator::m_eTau2Proposal {}
private

Definition at line 150 of file MissingMassCalculator.h.

150{},m_eTau2Proposal{};

◆ m_eTau2Range

double DiTauMassTools::MissingMassCalculator::m_eTau2Range {}
private

Definition at line 153 of file MissingMassCalculator.h.

◆ m_Ev1

double DiTauMassTools::MissingMassCalculator::m_Ev1 {}
private

Definition at line 191 of file MissingMassCalculator.h.

191{};

◆ m_Ev2

double DiTauMassTools::MissingMassCalculator::m_Ev2 {}
private

Definition at line 190 of file MissingMassCalculator.h.

190{};

◆ m_eventNumber

int DiTauMassTools::MissingMassCalculator::m_eventNumber {}
private

Definition at line 107 of file MissingMassCalculator.h.

107{};

◆ m_fDitauStuffFit

DitauStuff DiTauMassTools::MissingMassCalculator::m_fDitauStuffFit
private

Definition at line 244 of file MissingMassCalculator.h.

◆ m_fDitauStuffHisto

DitauStuff DiTauMassTools::MissingMassCalculator::m_fDitauStuffHisto
private

Definition at line 245 of file MissingMassCalculator.h.

◆ m_fFitting

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

Definition at line 229 of file MissingMassCalculator.h.

229{};

◆ m_fMEtL_all

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

Definition at line 198 of file MissingMassCalculator.h.

◆ m_fMEtL_split1

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

Definition at line 216 of file MissingMassCalculator.h.

◆ m_fMEtL_split2

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

Definition at line 223 of file MissingMassCalculator.h.

◆ m_fMEtP_all

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

Definition at line 197 of file MissingMassCalculator.h.

◆ m_fMEtP_split1

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

Definition at line 215 of file MissingMassCalculator.h.

◆ m_fMEtP_split2

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

Definition at line 222 of file MissingMassCalculator.h.

◆ m_fMetx

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

Definition at line 235 of file MissingMassCalculator.h.

235{};

◆ m_fMety

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

Definition at line 236 of file MissingMassCalculator.h.

236{};

◆ m_fMfit_all

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

Definition at line 196 of file MissingMassCalculator.h.

◆ m_fMfit_allGraph

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

Definition at line 203 of file MissingMassCalculator.h.

◆ m_fMfit_allNoWeight

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

Definition at line 204 of file MissingMassCalculator.h.

◆ m_fMmass_split1

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

Definition at line 214 of file MissingMassCalculator.h.

◆ m_fMmass_split2

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

Definition at line 221 of file MissingMassCalculator.h.

◆ m_fMnu1

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

Definition at line 233 of file MissingMassCalculator.h.

233{};

◆ m_fMnu1_all

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

Definition at line 199 of file MissingMassCalculator.h.

◆ m_fMnu1_split1

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

Definition at line 217 of file MissingMassCalculator.h.

◆ m_fMnu1_split2

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

Definition at line 224 of file MissingMassCalculator.h.

◆ m_fMnu2

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

Definition at line 234 of file MissingMassCalculator.h.

234{};

◆ m_fMnu2_all

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

Definition at line 200 of file MissingMassCalculator.h.

◆ m_fMnu2_split1

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

Definition at line 218 of file MissingMassCalculator.h.

◆ m_fMnu2_split2

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

Definition at line 225 of file MissingMassCalculator.h.

◆ m_fPhi1

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

Definition at line 231 of file MissingMassCalculator.h.

231{};

◆ m_fPhi1_all

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

Definition at line 201 of file MissingMassCalculator.h.

◆ m_fPhi1_split1

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

Definition at line 219 of file MissingMassCalculator.h.

◆ m_fPhi1_split2

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

Definition at line 226 of file MissingMassCalculator.h.

◆ m_fPhi2

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

Definition at line 232 of file MissingMassCalculator.h.

232{};

◆ m_fPhi2_all

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

Definition at line 202 of file MissingMassCalculator.h.

◆ m_fPhi2_split1

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

Definition at line 220 of file MissingMassCalculator.h.

◆ m_fPhi2_split2

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

Definition at line 227 of file MissingMassCalculator.h.

◆ m_fPXfit1

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

Definition at line 206 of file MissingMassCalculator.h.

◆ m_fPXfit2

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

Definition at line 209 of file MissingMassCalculator.h.

◆ m_fPYfit1

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

Definition at line 207 of file MissingMassCalculator.h.

◆ m_fPYfit2

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

Definition at line 210 of file MissingMassCalculator.h.

◆ m_fPZfit1

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

Definition at line 208 of file MissingMassCalculator.h.

◆ m_fPZfit2

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

Definition at line 211 of file MissingMassCalculator.h.

◆ m_fTauProb

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

Definition at line 238 of file MissingMassCalculator.h.

238{};

◆ m_fTheta3D

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

Definition at line 237 of file MissingMassCalculator.h.

237{};

◆ m_fullParamSpaceScan

bool DiTauMassTools::MissingMassCalculator::m_fullParamSpaceScan {}
private

Definition at line 154 of file MissingMassCalculator.h.

154{};

◆ m_fUseEfficiencyRecovery

bool DiTauMassTools::MissingMassCalculator::m_fUseEfficiencyRecovery {}
private

Definition at line 68 of file MissingMassCalculator.h.

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

◆ m_fUseFloatStopping

bool DiTauMassTools::MissingMassCalculator::m_fUseFloatStopping {}
private

Definition at line 69 of file MissingMassCalculator.h.

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

◆ m_fUseFloatStoppingCheckFreq

int DiTauMassTools::MissingMassCalculator::m_fUseFloatStoppingCheckFreq {}
private

Definition at line 71 of file MissingMassCalculator.h.

71{};

◆ m_fUseFloatStoppingComp

double DiTauMassTools::MissingMassCalculator::m_fUseFloatStoppingComp {}
private

Definition at line 72 of file MissingMassCalculator.h.

72{};

◆ m_fUseFloatStoppingMinIter

int DiTauMassTools::MissingMassCalculator::m_fUseFloatStoppingMinIter {}
private

Definition at line 70 of file MissingMassCalculator.h.

70{};

◆ 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 111 of file MissingMassCalculator.h.

111{};

◆ 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 119 of file MissingMassCalculator.h.

119{};

◆ m_iterNuPV3

int DiTauMassTools::MissingMassCalculator::m_iterNuPV3 {}
private

Definition at line 112 of file MissingMassCalculator.h.

112{};

◆ m_lfvLeplepRefit

bool DiTauMassTools::MissingMassCalculator::m_lfvLeplepRefit
private

Definition at line 94 of file MissingMassCalculator.h.

◆ m_m2Nu1

double DiTauMassTools::MissingMassCalculator::m_m2Nu1 {}
private

Definition at line 184 of file MissingMassCalculator.h.

184{};

◆ m_m2Nu2

double DiTauMassTools::MissingMassCalculator::m_m2Nu2 {}
private

Definition at line 185 of file MissingMassCalculator.h.

185{};

◆ m_markovCountDuplicate

int DiTauMassTools::MissingMassCalculator::m_markovCountDuplicate {}
private

Definition at line 125 of file MissingMassCalculator.h.

125{};

◆ m_markovNAccept

int DiTauMassTools::MissingMassCalculator::m_markovNAccept {}
private

Definition at line 129 of file MissingMassCalculator.h.

129{};

◆ m_markovNFullScan

int DiTauMassTools::MissingMassCalculator::m_markovNFullScan {}
private

Definition at line 126 of file MissingMassCalculator.h.

126{};

◆ m_markovNRejectMetropolis

int DiTauMassTools::MissingMassCalculator::m_markovNRejectMetropolis {}
private

Definition at line 128 of file MissingMassCalculator.h.

128{};

◆ m_markovNRejectNoSol

int DiTauMassTools::MissingMassCalculator::m_markovNRejectNoSol {}
private

Definition at line 127 of file MissingMassCalculator.h.

127{};

◆ m_meanbinStop

double DiTauMassTools::MissingMassCalculator::m_meanbinStop {}
private

Definition at line 78 of file MissingMassCalculator.h.

78{};

◆ m_meanbinToBeEvaluated

bool DiTauMassTools::MissingMassCalculator::m_meanbinToBeEvaluated {}
private

Definition at line 123 of file MissingMassCalculator.h.

123{};

◆ m_Meff

double DiTauMassTools::MissingMassCalculator::m_Meff {}
private

Definition at line 192 of file MissingMassCalculator.h.

192{},m_Meff{};

◆ m_metCovPhiCos

double DiTauMassTools::MissingMassCalculator::m_metCovPhiCos {}
private

Definition at line 149 of file MissingMassCalculator.h.

149{},m_metCovPhiSin{};

◆ m_metCovPhiSin

double DiTauMassTools::MissingMassCalculator::m_metCovPhiSin {}
private

Definition at line 149 of file MissingMassCalculator.h.

149{},m_metCovPhiSin{};

◆ m_MEtL

double DiTauMassTools::MissingMassCalculator::m_MEtL {}
private

Definition at line 140 of file MissingMassCalculator.h.

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

◆ m_MEtL0

double DiTauMassTools::MissingMassCalculator::m_MEtL0 {}
private

Definition at line 143 of file MissingMassCalculator.h.

◆ m_MEtLMax

double DiTauMassTools::MissingMassCalculator::m_MEtLMax {}
private

Definition at line 145 of file MissingMassCalculator.h.

◆ m_MEtLMin

double DiTauMassTools::MissingMassCalculator::m_MEtLMin {}
private

Definition at line 144 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 140 of file MissingMassCalculator.h.

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

◆ m_MEtP0

double DiTauMassTools::MissingMassCalculator::m_MEtP0 {}
private

Definition at line 143 of file MissingMassCalculator.h.

◆ m_MEtPMax

double DiTauMassTools::MissingMassCalculator::m_MEtPMax {}
private

Definition at line 145 of file MissingMassCalculator.h.

◆ m_MEtPMin

double DiTauMassTools::MissingMassCalculator::m_MEtPMin {}
private

Definition at line 144 of file MissingMassCalculator.h.

◆ m_MEtPRange

double DiTauMassTools::MissingMassCalculator::m_MEtPRange {}
private

◆ m_MEtProposal

double DiTauMassTools::MissingMassCalculator::m_MEtProposal {}
private

Definition at line 148 of file MissingMassCalculator.h.

◆ m_MEtPStep

double DiTauMassTools::MissingMassCalculator::m_MEtPStep {}
private

Definition at line 146 of file MissingMassCalculator.h.

◆ m_mmcCalibrationSet

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

Definition at line 66 of file MissingMassCalculator.h.

66{};

◆ m_Mnu1

double DiTauMassTools::MissingMassCalculator::m_Mnu1 {}
private

Definition at line 140 of file MissingMassCalculator.h.

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

◆ m_Mnu10

double DiTauMassTools::MissingMassCalculator::m_Mnu10 {}
private

Definition at line 143 of file MissingMassCalculator.h.

◆ m_Mnu1Exclude

bool DiTauMassTools::MissingMassCalculator::m_Mnu1Exclude {}
private

Definition at line 155 of file MissingMassCalculator.h.

155{};

◆ m_Mnu1ExcludeMax

double DiTauMassTools::MissingMassCalculator::m_Mnu1ExcludeMax {}
private

◆ m_Mnu1ExcludeMin

double DiTauMassTools::MissingMassCalculator::m_Mnu1ExcludeMin {}
private

Definition at line 169 of file MissingMassCalculator.h.

◆ m_Mnu1ExcludeRange

double DiTauMassTools::MissingMassCalculator::m_Mnu1ExcludeRange {}
private

Definition at line 169 of file MissingMassCalculator.h.

◆ m_Mnu1Max

double DiTauMassTools::MissingMassCalculator::m_Mnu1Max {}
private

Definition at line 145 of file MissingMassCalculator.h.

◆ m_Mnu1Min

double DiTauMassTools::MissingMassCalculator::m_Mnu1Min {}
private

Definition at line 144 of file MissingMassCalculator.h.

◆ m_Mnu1Range

double DiTauMassTools::MissingMassCalculator::m_Mnu1Range {}
private

◆ m_Mnu1Step

double DiTauMassTools::MissingMassCalculator::m_Mnu1Step {}
private

Definition at line 146 of file MissingMassCalculator.h.

◆ m_Mnu1XMax

double DiTauMassTools::MissingMassCalculator::m_Mnu1XMax {}
private

◆ m_Mnu1XMin

double DiTauMassTools::MissingMassCalculator::m_Mnu1XMin {}
private

Definition at line 170 of file MissingMassCalculator.h.

◆ m_Mnu1XRange

double DiTauMassTools::MissingMassCalculator::m_Mnu1XRange {}
private

Definition at line 170 of file MissingMassCalculator.h.

◆ m_Mnu2

double DiTauMassTools::MissingMassCalculator::m_Mnu2 {}
private

Definition at line 140 of file MissingMassCalculator.h.

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

◆ m_Mnu20

double DiTauMassTools::MissingMassCalculator::m_Mnu20 {}
private

Definition at line 143 of file MissingMassCalculator.h.

◆ m_Mnu2Max

double DiTauMassTools::MissingMassCalculator::m_Mnu2Max {}
private

Definition at line 145 of file MissingMassCalculator.h.

◆ m_Mnu2Min

double DiTauMassTools::MissingMassCalculator::m_Mnu2Min {}
private

Definition at line 144 of file MissingMassCalculator.h.

◆ m_Mnu2Range

double DiTauMassTools::MissingMassCalculator::m_Mnu2Range {}
private

◆ m_Mnu2Step

double DiTauMassTools::MissingMassCalculator::m_Mnu2Step {}
private

Definition at line 146 of file MissingMassCalculator.h.

◆ m_MnuProposal

double DiTauMassTools::MissingMassCalculator::m_MnuProposal {}
private

Definition at line 148 of file MissingMassCalculator.h.

◆ m_MnuScanRange

double DiTauMassTools::MissingMassCalculator::m_MnuScanRange {}
private

Definition at line 257 of file MissingMassCalculator.h.

257{}; // 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 139 of file MissingMassCalculator.h.

139{},m_mTau2{};

◆ m_mTau2

double DiTauMassTools::MissingMassCalculator::m_mTau2 {}
private

Definition at line 139 of file MissingMassCalculator.h.

139{},m_mTau2{};

◆ m_mtautauFinalSolOldVec

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

Definition at line 158 of file MissingMassCalculator.h.

◆ m_mtautauFinalSolVec

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

Definition at line 164 of file MissingMassCalculator.h.

◆ m_mtautauSum

double DiTauMassTools::MissingMassCalculator::m_mtautauSum {}
private

Definition at line 105 of file MissingMassCalculator.h.

105{};

◆ m_Mvis

double DiTauMassTools::MissingMassCalculator::m_Mvis {}
private

Definition at line 192 of file MissingMassCalculator.h.

192{},m_Meff{};

◆ m_niter_fit1

int DiTauMassTools::MissingMassCalculator::m_niter_fit1 {}
private

Definition at line 247 of file MissingMassCalculator.h.

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

◆ m_niter_fit2

int DiTauMassTools::MissingMassCalculator::m_niter_fit2 {}
private

Definition at line 248 of file MissingMassCalculator.h.

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

◆ m_niter_fit3

int DiTauMassTools::MissingMassCalculator::m_niter_fit3 {}
private

Definition at line 249 of file MissingMassCalculator.h.

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

◆ m_NiterRandom

int DiTauMassTools::MissingMassCalculator::m_NiterRandom {}
private

Definition at line 250 of file MissingMassCalculator.h.

250{}; // 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 75 of file MissingMassCalculator.h.

75{};

◆ m_nosol1

int DiTauMassTools::MissingMassCalculator::m_nosol1 {}
private

Definition at line 117 of file MissingMassCalculator.h.

117{};

◆ m_nosol2

int DiTauMassTools::MissingMassCalculator::m_nosol2 {}
private

Definition at line 118 of file MissingMassCalculator.h.

118{};

◆ m_nsigma_METscan

double DiTauMassTools::MissingMassCalculator::m_nsigma_METscan {}
private

Definition at line 96 of file MissingMassCalculator.h.

◆ m_nsigma_METscan2

double DiTauMassTools::MissingMassCalculator::m_nsigma_METscan2 {}
private

Definition at line 96 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 98 of file MissingMassCalculator.h.

98{}, 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 96 of file MissingMassCalculator.h.

◆ m_nsol

int DiTauMassTools::MissingMassCalculator::m_nsol {}
private

Definition at line 162 of file MissingMassCalculator.h.

162{};

◆ m_nsolfinalmax

int DiTauMassTools::MissingMassCalculator::m_nsolfinalmax {}
private

Definition at line 74 of file MissingMassCalculator.h.

74{};

◆ m_nsolmax

int DiTauMassTools::MissingMassCalculator::m_nsolmax
private

Definition at line 74 of file MissingMassCalculator.h.

◆ m_nsolOld

int DiTauMassTools::MissingMassCalculator::m_nsolOld {}
private

Definition at line 156 of file MissingMassCalculator.h.

156{};

◆ m_NsucStop

int DiTauMassTools::MissingMassCalculator::m_NsucStop {}
private

Definition at line 251 of file MissingMassCalculator.h.

251{};

◆ m_nsucStop

int DiTauMassTools::MissingMassCalculator::m_nsucStop {}
private

Definition at line 76 of file MissingMassCalculator.h.

76{};

◆ m_nu1FinalSolOldVec

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

Definition at line 159 of file MissingMassCalculator.h.

◆ m_nu1FinalSolVec

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

Definition at line 165 of file MissingMassCalculator.h.

◆ m_nu2FinalSolOldVec

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

Definition at line 160 of file MissingMassCalculator.h.

◆ m_nu2FinalSolVec

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

Definition at line 166 of file MissingMassCalculator.h.

◆ m_nuvec1_tmp

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

Definition at line 89 of file MissingMassCalculator.h.

◆ m_nuvec2_tmp

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

Definition at line 90 of file MissingMassCalculator.h.

◆ m_nuvecsol1

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

Definition at line 81 of file MissingMassCalculator.h.

◆ m_nuvecsol2

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

Definition at line 82 of file MissingMassCalculator.h.

◆ m_Phi1

double DiTauMassTools::MissingMassCalculator::m_Phi1 {}
private

Definition at line 140 of file MissingMassCalculator.h.

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

◆ m_Phi10

double DiTauMassTools::MissingMassCalculator::m_Phi10 {}
private

Definition at line 143 of file MissingMassCalculator.h.

◆ m_Phi1Max

double DiTauMassTools::MissingMassCalculator::m_Phi1Max {}
private

Definition at line 145 of file MissingMassCalculator.h.

◆ m_Phi1Min

double DiTauMassTools::MissingMassCalculator::m_Phi1Min {}
private

Definition at line 144 of file MissingMassCalculator.h.

◆ m_Phi1Range

double DiTauMassTools::MissingMassCalculator::m_Phi1Range {}
private

◆ m_Phi1Step

double DiTauMassTools::MissingMassCalculator::m_Phi1Step {}
private

Definition at line 146 of file MissingMassCalculator.h.

◆ m_Phi2

double DiTauMassTools::MissingMassCalculator::m_Phi2 {}
private

Definition at line 140 of file MissingMassCalculator.h.

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

◆ m_Phi20

double DiTauMassTools::MissingMassCalculator::m_Phi20 {}
private

Definition at line 143 of file MissingMassCalculator.h.

◆ m_Phi2Max

double DiTauMassTools::MissingMassCalculator::m_Phi2Max {}
private

Definition at line 145 of file MissingMassCalculator.h.

◆ m_Phi2Min

double DiTauMassTools::MissingMassCalculator::m_Phi2Min {}
private

Definition at line 144 of file MissingMassCalculator.h.

◆ m_Phi2Range

double DiTauMassTools::MissingMassCalculator::m_Phi2Range {}
private

◆ m_Phi2Step

double DiTauMassTools::MissingMassCalculator::m_Phi2Step {}
private

Definition at line 146 of file MissingMassCalculator.h.

◆ m_PhiProposal

double DiTauMassTools::MissingMassCalculator::m_PhiProposal {}
private

Definition at line 148 of file MissingMassCalculator.h.

◆ m_PrintmInvWidth2Error

double DiTauMassTools::MissingMassCalculator::m_PrintmInvWidth2Error {}
private

Definition at line 133 of file MissingMassCalculator.h.

133{};

◆ m_PrintmMaxError

double DiTauMassTools::MissingMassCalculator::m_PrintmMaxError {}
private

Definition at line 131 of file MissingMassCalculator.h.

131{};

◆ m_PrintmMeanError

double DiTauMassTools::MissingMassCalculator::m_PrintmMeanError {}
private

Definition at line 132 of file MissingMassCalculator.h.

132{};

◆ m_prob_tmp

double DiTauMassTools::MissingMassCalculator::m_prob_tmp {}
private

Definition at line 102 of file MissingMassCalculator.h.

102{};

◆ m_probFinalSolOldVec

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

Definition at line 157 of file MissingMassCalculator.h.

◆ m_probFinalSolVec

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

Definition at line 163 of file MissingMassCalculator.h.

◆ m_proposalTryMEt

double DiTauMassTools::MissingMassCalculator::m_proposalTryMEt {}
private

Definition at line 135 of file MissingMassCalculator.h.

135{};

◆ m_ProposalTryMnu

double DiTauMassTools::MissingMassCalculator::m_ProposalTryMnu {}
private

Definition at line 137 of file MissingMassCalculator.h.

137{};

◆ m_ProposalTryPhi

double DiTauMassTools::MissingMassCalculator::m_ProposalTryPhi {}
private

Definition at line 136 of file MissingMassCalculator.h.

136{};

◆ m_randomGen

TRandom2 DiTauMassTools::MissingMassCalculator::m_randomGen
private

Definition at line 64 of file MissingMassCalculator.h.

◆ m_RMSStop

int DiTauMassTools::MissingMassCalculator::m_RMSStop {}
private

Definition at line 252 of file MissingMassCalculator.h.

252{};

◆ m_rmsStop

int DiTauMassTools::MissingMassCalculator::m_rmsStop {}
private

Definition at line 77 of file MissingMassCalculator.h.

77{};

◆ m_RndmSeedAltering

int DiTauMassTools::MissingMassCalculator::m_RndmSeedAltering {}
private

Definition at line 253 of file MissingMassCalculator.h.

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

◆ m_SaveLlhHisto

bool DiTauMassTools::MissingMassCalculator::m_SaveLlhHisto
private

Definition at line 94 of file MissingMassCalculator.h.

◆ m_scanMnu1

bool DiTauMassTools::MissingMassCalculator::m_scanMnu1 {}
private

Definition at line 174 of file MissingMassCalculator.h.

174{},m_scanMnu2{};

◆ m_scanMnu2

bool DiTauMassTools::MissingMassCalculator::m_scanMnu2 {}
private

Definition at line 174 of file MissingMassCalculator.h.

174{},m_scanMnu2{};

◆ m_seed

int DiTauMassTools::MissingMassCalculator::m_seed {}
private

Definition at line 108 of file MissingMassCalculator.h.

108{};

◆ m_sinPhi1

double DiTauMassTools::MissingMassCalculator::m_sinPhi1 {}
private

Definition at line 172 of file MissingMassCalculator.h.

172{}, m_cosPhi2{}, m_sinPhi1{}, m_sinPhi2{};

◆ m_sinPhi2

double DiTauMassTools::MissingMassCalculator::m_sinPhi2 {}
private

Definition at line 172 of file MissingMassCalculator.h.

172{}, m_cosPhi2{}, m_sinPhi1{}, m_sinPhi2{};

◆ m_switch1

bool DiTauMassTools::MissingMassCalculator::m_switch1 {}
private

Definition at line 120 of file MissingMassCalculator.h.

120{};

◆ m_switch2

bool DiTauMassTools::MissingMassCalculator::m_switch2 {}
private

Definition at line 121 of file MissingMassCalculator.h.

121{};

◆ m_tautau_tmp

PtEtaPhiMVector DiTauMassTools::MissingMassCalculator::m_tautau_tmp
private

Definition at line 92 of file MissingMassCalculator.h.

◆ m_tauVec1

PtEtaPhiMVector DiTauMassTools::MissingMassCalculator::m_tauVec1
private

Definition at line 176 of file MissingMassCalculator.h.

◆ m_tauVec1E

double DiTauMassTools::MissingMassCalculator::m_tauVec1E {}
private

Definition at line 182 of file MissingMassCalculator.h.

182{};

◆ m_tauVec1M

double DiTauMassTools::MissingMassCalculator::m_tauVec1M {}
private

Definition at line 178 of file MissingMassCalculator.h.

178{}, m_tauVec2M{};

◆ m_tauVec1P

double DiTauMassTools::MissingMassCalculator::m_tauVec1P {}
private

Definition at line 181 of file MissingMassCalculator.h.

181{}, m_tauVec2P{};

◆ m_tauVec1Phi

double DiTauMassTools::MissingMassCalculator::m_tauVec1Phi {}
private

Definition at line 177 of file MissingMassCalculator.h.

177{}, m_tauVec2Phi{};

◆ m_tauVec1Px

double DiTauMassTools::MissingMassCalculator::m_tauVec1Px {}
private

Definition at line 179 of file MissingMassCalculator.h.

179{}, m_tauVec1Py{}, m_tauVec1Pz{};

◆ m_tauVec1Py

double DiTauMassTools::MissingMassCalculator::m_tauVec1Py {}
private

Definition at line 179 of file MissingMassCalculator.h.

179{}, m_tauVec1Py{}, m_tauVec1Pz{};

◆ m_tauVec1Pz

double DiTauMassTools::MissingMassCalculator::m_tauVec1Pz {}
private

Definition at line 179 of file MissingMassCalculator.h.

179{}, m_tauVec1Py{}, m_tauVec1Pz{};

◆ m_tauVec2

PtEtaPhiMVector DiTauMassTools::MissingMassCalculator::m_tauVec2
private

Definition at line 176 of file MissingMassCalculator.h.

◆ m_tauVec2E

double DiTauMassTools::MissingMassCalculator::m_tauVec2E {}
private

Definition at line 183 of file MissingMassCalculator.h.

183{};

◆ m_tauVec2M

double DiTauMassTools::MissingMassCalculator::m_tauVec2M {}
private

Definition at line 178 of file MissingMassCalculator.h.

178{}, m_tauVec2M{};

◆ m_tauVec2P

double DiTauMassTools::MissingMassCalculator::m_tauVec2P {}
private

Definition at line 181 of file MissingMassCalculator.h.

181{}, m_tauVec2P{};

◆ m_tauVec2Phi

double DiTauMassTools::MissingMassCalculator::m_tauVec2Phi {}
private

Definition at line 177 of file MissingMassCalculator.h.

177{}, m_tauVec2Phi{};

◆ m_tauVec2Px

double DiTauMassTools::MissingMassCalculator::m_tauVec2Px {}
private

Definition at line 180 of file MissingMassCalculator.h.

180{}, m_tauVec2Py{}, m_tauVec2Pz{};

◆ m_tauVec2Py

double DiTauMassTools::MissingMassCalculator::m_tauVec2Py {}
private

Definition at line 180 of file MissingMassCalculator.h.

180{}, m_tauVec2Py{}, m_tauVec2Pz{};

◆ m_tauVec2Pz

double DiTauMassTools::MissingMassCalculator::m_tauVec2Pz {}
private

Definition at line 180 of file MissingMassCalculator.h.

180{}, m_tauVec2Py{}, m_tauVec2Pz{};

◆ m_tauvecprob1

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

Definition at line 86 of file MissingMassCalculator.h.

◆ m_tauvecprob2

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

Definition at line 87 of file MissingMassCalculator.h.

◆ m_tauvecsol1

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

Definition at line 84 of file MissingMassCalculator.h.

◆ m_tauvecsol2

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

Definition at line 85 of file MissingMassCalculator.h.

◆ m_testdiscri1

int DiTauMassTools::MissingMassCalculator::m_testdiscri1 {}
private

Definition at line 115 of file MissingMassCalculator.h.

115{};

◆ m_testdiscri2

int DiTauMassTools::MissingMassCalculator::m_testdiscri2 {}
private

Definition at line 116 of file MissingMassCalculator.h.

116{};

◆ m_testptn1

int DiTauMassTools::MissingMassCalculator::m_testptn1 {}
private

Definition at line 113 of file MissingMassCalculator.h.

113{};

◆ m_testptn2

int DiTauMassTools::MissingMassCalculator::m_testptn2 {}
private

Definition at line 114 of file MissingMassCalculator.h.

114{};

◆ m_TLVdummy

PtEtaPhiMVector DiTauMassTools::MissingMassCalculator::m_TLVdummy
private

Definition at line 241 of file MissingMassCalculator.h.

◆ m_totalProbSum

double DiTauMassTools::MissingMassCalculator::m_totalProbSum {}
private

Definition at line 104 of file MissingMassCalculator.h.

104{};

◆ m_walkWeight

double DiTauMassTools::MissingMassCalculator::m_walkWeight {}
private

Definition at line 171 of file MissingMassCalculator.h.

171{};

◆ metvec_tmp

XYVector DiTauMassTools::MissingMassCalculator::metvec_tmp

Definition at line 416 of file MissingMassCalculator.h.

◆ OutputInfo

MissingMassOutput DiTauMassTools::MissingMassCalculator::OutputInfo

Definition at line 334 of file MissingMassCalculator.h.

◆ preparedInput

MissingMassInput DiTauMassTools::MissingMassCalculator::preparedInput

Definition at line 333 of file MissingMassCalculator.h.

◆ Prob

MissingMassProb* DiTauMassTools::MissingMassCalculator::Prob

Definition at line 335 of file MissingMassCalculator.h.


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