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

#include <MissingMassCalculator.h>

Collaboration diagram for DiTauMassTools::MissingMassCalculator:

Classes

struct  DitauStuff

Public Member Functions

 ~MissingMassCalculator ()
 MissingMassCalculator (MMCCalibrationSet::e aset, std::string paramFilePath)
 MissingMassCalculator (const MissingMassCalculator &)=delete
MissingMassCalculatoroperator= (const MissingMassCalculator &)=delete
int RunMissingMassCalculator (const xAOD::IParticle *part1, const xAOD::IParticle *part2, const xAOD::MissingET *met, const int &njets)
void FinalizeSettings (const xAOD::IParticle *part1, const xAOD::IParticle *part2, const xAOD::MissingET *met, const int &njets)
void SetNiterFit1 (const int val)
void SetNiterFit2 (const int val)
void SetNiterFit3 (const int val)
void SetNiterRandom (const int val)
void SetNsucStop (const int val)
void SetRMSStop (const int val)
void SetMeanbinStop (const double val)
void SetRndmSeedAltering (const int val)
void 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 46 of file MissingMassCalculator.h.

Constructor & Destructor Documentation

◆ ~MissingMassCalculator()

MissingMassCalculator::~MissingMassCalculator ( )

Definition at line 187 of file MissingMassCalculator.cxx.

187{ delete Prob; }

◆ MissingMassCalculator() [1/2]

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

Definition at line 48 of file MissingMassCalculator.cxx.

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

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

◆ checkMEtInRange()

bool MissingMassCalculator::checkMEtInRange ( )
inlineprotected

Definition at line 2697 of file MissingMassCalculator.cxx.

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

◆ CheckSolutions()

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

◆ ClearDitauStuff()

void MissingMassCalculator::ClearDitauStuff ( DitauStuff & fStuff)
private

Definition at line 284 of file MissingMassCalculator.cxx.

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

◆ DitauMassCalculatorV9lfv()

int MissingMassCalculator::DitauMassCalculatorV9lfv ( bool refit)
inlineprotected

Definition at line 1005 of file MissingMassCalculator.cxx.

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

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

◆ DoOutputInfo()

void MissingMassCalculator::DoOutputInfo ( )
private

Definition at line 299 of file MissingMassCalculator.cxx.

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

◆ dTheta3DLimit()

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

Definition at line 2713 of file MissingMassCalculator.cxx.

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

◆ FinalizeSettings()

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

Definition at line 2817 of file MissingMassCalculator.cxx.

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

372{ return m_markovNAccept; }

◆ GetMarkovNFullscan()

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

Definition at line 373 of file MissingMassCalculator.h.

373{ return m_markovNFullScan;}

◆ GetMarkovNRejectMetropolis()

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

Definition at line 371 of file MissingMassCalculator.h.

◆ GetMarkovNRejectNoSol()

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

Definition at line 370 of file MissingMassCalculator.h.

370{ return m_markovNRejectNoSol;}

◆ GetMeanbinStop()

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

Definition at line 366 of file MissingMassCalculator.h.

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

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

◆ GetNiterFit2()

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

Definition at line 360 of file MissingMassCalculator.h.

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

◆ GetNiterFit3()

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

Definition at line 361 of file MissingMassCalculator.h.

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

◆ GetNiterRandom()

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

Definition at line 362 of file MissingMassCalculator.h.

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

◆ GetNMetroReject()

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

Definition at line 396 of file MissingMassCalculator.h.

◆ GetNNoSol()

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

Definition at line 395 of file MissingMassCalculator.h.

395{return m_markovNRejectNoSol;}

◆ GetNSol()

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

Definition at line 397 of file MissingMassCalculator.h.

397{return m_markovNAccept;}

◆ GetNsucStop()

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

Definition at line 364 of file MissingMassCalculator.h.

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

◆ GetProposalTryMEt()

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

Definition at line 374 of file MissingMassCalculator.h.

374{return m_proposalTryMEt;}

◆ GetProposalTryMnu()

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

Definition at line 376 of file MissingMassCalculator.h.

376{return m_ProposalTryMnu;}

◆ GetProposalTryPhi()

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

Definition at line 375 of file MissingMassCalculator.h.

375{return m_ProposalTryPhi;}

◆ GetRMSStop()

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

Definition at line 365 of file MissingMassCalculator.h.

365{ return m_RMSStop; }

◆ GetRndmSeedAltering()

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

Definition at line 367 of file MissingMassCalculator.h.

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

◆ GetUseEfficiencyRecovery()

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

Definition at line 357 of file MissingMassCalculator.h.

357{ return m_fUseEfficiencyRecovery; }

◆ handleSolutions()

void MissingMassCalculator::handleSolutions ( )
inlineprotected

Definition at line 2038 of file MissingMassCalculator.cxx.

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

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

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

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

405 {
406 return maxFromHist(theHist.get(), histInfo, maxHistStrategy, winHalfWidth, debug);
407 }
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 1510 of file MissingMassCalculator.cxx.

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

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

◆ NuPsolutionV3()

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

Definition at line 543 of file MissingMassCalculator.cxx.

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

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

◆ PrintOtherInput()

void MissingMassCalculator::PrintOtherInput ( )
private

Definition at line 396 of file MissingMassCalculator.cxx.

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

◆ PrintResults()

void MissingMassCalculator::PrintResults ( )
private

Definition at line 438 of file MissingMassCalculator.cxx.

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

◆ probCalculatorV9fast()

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

Definition at line 1786 of file MissingMassCalculator.cxx.

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

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

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

◆ SaveLlhHisto()

void MissingMassCalculator::SaveLlhHisto ( const bool val)

Definition at line 3086 of file MissingMassCalculator.cxx.

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

◆ SetBeamEnergy()

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

Definition at line 387 of file MissingMassCalculator.h.

387{ m_beamEnergy=val; }

◆ SetEventNumber()

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

Definition at line 348 of file MissingMassCalculator.h.

◆ SetFloatStoppingCheckFreq()

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

Definition at line 385 of file MissingMassCalculator.h.

◆ SetFloatStoppingComp()

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

Definition at line 386 of file MissingMassCalculator.h.

◆ SetFloatStoppingMinIter()

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

Definition at line 384 of file MissingMassCalculator.h.

◆ SetLFVLeplepRefit()

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

Definition at line 388 of file MissingMassCalculator.h.

◆ SetMeanbinStop()

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

Definition at line 346 of file MissingMassCalculator.h.

◆ SetMnuScanRange()

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

Definition at line 350 of file MissingMassCalculator.h.

◆ SetNiterFit1()

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

Definition at line 340 of file MissingMassCalculator.h.

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

◆ SetNiterFit2()

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

Definition at line 341 of file MissingMassCalculator.h.

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

◆ SetNiterFit3()

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

Definition at line 342 of file MissingMassCalculator.h.

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

◆ SetNiterRandom()

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

Definition at line 343 of file MissingMassCalculator.h.

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

◆ SetNsigmaMETscan()

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

Definition at line 381 of file MissingMassCalculator.h.

381{ 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 380 of file MissingMassCalculator.h.

380{ 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 379 of file MissingMassCalculator.h.

379{ 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 378 of file MissingMassCalculator.h.

378{ 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 344 of file MissingMassCalculator.h.

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

◆ SetProposalTryMEt()

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

Definition at line 352 of file MissingMassCalculator.h.

◆ SetProposalTryMnu()

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

Definition at line 354 of file MissingMassCalculator.h.

◆ SetProposalTryPhi()

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

Definition at line 353 of file MissingMassCalculator.h.

◆ SetRMSStop()

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

Definition at line 345 of file MissingMassCalculator.h.

345{ m_RMSStop=val;}

◆ SetRndmSeedAltering()

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

Definition at line 347 of file MissingMassCalculator.h.

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

◆ SetUseEfficiencyRecovery()

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

Definition at line 356 of file MissingMassCalculator.h.

◆ SetUseFloatStopping()

void MissingMassCalculator::SetUseFloatStopping ( const bool val)

Definition at line 3121 of file MissingMassCalculator.cxx.

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

◆ SpaceWalkerInit()

void MissingMassCalculator::SpaceWalkerInit ( )
inlineprotected

Definition at line 2307 of file MissingMassCalculator.cxx.

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

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

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

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

Member Data Documentation

◆ m_beamEnergy

double DiTauMassTools::MissingMassCalculator::m_beamEnergy {}
private

Definition at line 97 of file MissingMassCalculator.h.

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

◆ m_cosPhi1

double DiTauMassTools::MissingMassCalculator::m_cosPhi1 {}
private

Definition at line 171 of file MissingMassCalculator.h.

171{}, m_cosPhi2{}, m_sinPhi1{}, m_sinPhi2{};

◆ m_cosPhi2

double DiTauMassTools::MissingMassCalculator::m_cosPhi2 {}
private

Definition at line 171 of file MissingMassCalculator.h.

171{}, m_cosPhi2{}, m_sinPhi1{}, m_sinPhi2{};

◆ m_debugThisIteration

bool DiTauMassTools::MissingMassCalculator::m_debugThisIteration
private

Definition at line 93 of file MissingMassCalculator.h.

◆ m_dRmax_tau

double DiTauMassTools::MissingMassCalculator::m_dRmax_tau {}
private

Definition at line 254 of file MissingMassCalculator.h.

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

◆ m_E2v1

double DiTauMassTools::MissingMassCalculator::m_E2v1 {}
private

Definition at line 187 of file MissingMassCalculator.h.

187{};

◆ m_E2v2

double DiTauMassTools::MissingMassCalculator::m_E2v2 {}
private

Definition at line 188 of file MissingMassCalculator.h.

188{};

◆ m_ET2v1

double DiTauMassTools::MissingMassCalculator::m_ET2v1 {}
private

Definition at line 185 of file MissingMassCalculator.h.

185{};

◆ m_ET2v2

double DiTauMassTools::MissingMassCalculator::m_ET2v2 {}
private

Definition at line 186 of file MissingMassCalculator.h.

186{};

◆ m_eTau1

double DiTauMassTools::MissingMassCalculator::m_eTau1 {}
private

Definition at line 140 of file MissingMassCalculator.h.

140{}, m_eTau2{};

◆ m_eTau10

double DiTauMassTools::MissingMassCalculator::m_eTau10 {}
private

Definition at line 141 of file MissingMassCalculator.h.

141{}, m_eTau20{};

◆ m_eTau1Max

double DiTauMassTools::MissingMassCalculator::m_eTau1Max {}
private

◆ m_eTau1Min

double DiTauMassTools::MissingMassCalculator::m_eTau1Min {}
private

Definition at line 151 of file MissingMassCalculator.h.

◆ m_eTau1Proposal

double DiTauMassTools::MissingMassCalculator::m_eTau1Proposal {}
private

◆ m_eTau1Range

double DiTauMassTools::MissingMassCalculator::m_eTau1Range {}
private

Definition at line 151 of file MissingMassCalculator.h.

◆ m_eTau2

double DiTauMassTools::MissingMassCalculator::m_eTau2 {}
private

Definition at line 140 of file MissingMassCalculator.h.

140{}, m_eTau2{};

◆ m_eTau20

double DiTauMassTools::MissingMassCalculator::m_eTau20 {}
private

Definition at line 141 of file MissingMassCalculator.h.

141{}, m_eTau20{};

◆ m_eTau2Max

double DiTauMassTools::MissingMassCalculator::m_eTau2Max {}
private

◆ m_eTau2Min

double DiTauMassTools::MissingMassCalculator::m_eTau2Min {}
private

Definition at line 152 of file MissingMassCalculator.h.

◆ m_eTau2Proposal

double DiTauMassTools::MissingMassCalculator::m_eTau2Proposal {}
private

Definition at line 149 of file MissingMassCalculator.h.

149{},m_eTau2Proposal{};

◆ m_eTau2Range

double DiTauMassTools::MissingMassCalculator::m_eTau2Range {}
private

Definition at line 152 of file MissingMassCalculator.h.

◆ m_Ev1

double DiTauMassTools::MissingMassCalculator::m_Ev1 {}
private

Definition at line 190 of file MissingMassCalculator.h.

190{};

◆ m_Ev2

double DiTauMassTools::MissingMassCalculator::m_Ev2 {}
private

Definition at line 189 of file MissingMassCalculator.h.

189{};

◆ m_eventNumber

int DiTauMassTools::MissingMassCalculator::m_eventNumber {}
private

Definition at line 106 of file MissingMassCalculator.h.

106{};

◆ m_fDitauStuffFit

DitauStuff DiTauMassTools::MissingMassCalculator::m_fDitauStuffFit
private

Definition at line 243 of file MissingMassCalculator.h.

◆ m_fDitauStuffHisto

DitauStuff DiTauMassTools::MissingMassCalculator::m_fDitauStuffHisto
private

Definition at line 244 of file MissingMassCalculator.h.

◆ m_fFitting

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

Definition at line 228 of file MissingMassCalculator.h.

228{};

◆ m_fMEtL_all

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

Definition at line 197 of file MissingMassCalculator.h.

◆ m_fMEtL_split1

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

Definition at line 215 of file MissingMassCalculator.h.

◆ m_fMEtL_split2

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

Definition at line 222 of file MissingMassCalculator.h.

◆ m_fMEtP_all

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

Definition at line 196 of file MissingMassCalculator.h.

◆ m_fMEtP_split1

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

Definition at line 214 of file MissingMassCalculator.h.

◆ m_fMEtP_split2

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

Definition at line 221 of file MissingMassCalculator.h.

◆ m_fMetx

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

Definition at line 234 of file MissingMassCalculator.h.

234{};

◆ m_fMety

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

Definition at line 235 of file MissingMassCalculator.h.

235{};

◆ m_fMfit_all

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

Definition at line 195 of file MissingMassCalculator.h.

◆ m_fMfit_allGraph

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

Definition at line 202 of file MissingMassCalculator.h.

◆ m_fMfit_allNoWeight

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

Definition at line 203 of file MissingMassCalculator.h.

◆ m_fMmass_split1

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

Definition at line 213 of file MissingMassCalculator.h.

◆ m_fMmass_split2

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

Definition at line 220 of file MissingMassCalculator.h.

◆ m_fMnu1

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

Definition at line 232 of file MissingMassCalculator.h.

232{};

◆ m_fMnu1_all

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

Definition at line 198 of file MissingMassCalculator.h.

◆ m_fMnu1_split1

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

Definition at line 216 of file MissingMassCalculator.h.

◆ m_fMnu1_split2

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

Definition at line 223 of file MissingMassCalculator.h.

◆ m_fMnu2

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

Definition at line 233 of file MissingMassCalculator.h.

233{};

◆ m_fMnu2_all

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

Definition at line 199 of file MissingMassCalculator.h.

◆ m_fMnu2_split1

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

Definition at line 217 of file MissingMassCalculator.h.

◆ m_fMnu2_split2

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

Definition at line 224 of file MissingMassCalculator.h.

◆ m_fPhi1

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

Definition at line 230 of file MissingMassCalculator.h.

230{};

◆ m_fPhi1_all

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

Definition at line 200 of file MissingMassCalculator.h.

◆ m_fPhi1_split1

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

Definition at line 218 of file MissingMassCalculator.h.

◆ m_fPhi1_split2

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

Definition at line 225 of file MissingMassCalculator.h.

◆ m_fPhi2

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

Definition at line 231 of file MissingMassCalculator.h.

231{};

◆ m_fPhi2_all

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

Definition at line 201 of file MissingMassCalculator.h.

◆ m_fPhi2_split1

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

Definition at line 219 of file MissingMassCalculator.h.

◆ m_fPhi2_split2

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

Definition at line 226 of file MissingMassCalculator.h.

◆ m_fPXfit1

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

Definition at line 205 of file MissingMassCalculator.h.

◆ m_fPXfit2

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

Definition at line 208 of file MissingMassCalculator.h.

◆ m_fPYfit1

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

Definition at line 206 of file MissingMassCalculator.h.

◆ m_fPYfit2

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

Definition at line 209 of file MissingMassCalculator.h.

◆ m_fPZfit1

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

Definition at line 207 of file MissingMassCalculator.h.

◆ m_fPZfit2

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

Definition at line 210 of file MissingMassCalculator.h.

◆ m_fTauProb

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

Definition at line 237 of file MissingMassCalculator.h.

237{};

◆ m_fTheta3D

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

Definition at line 236 of file MissingMassCalculator.h.

236{};

◆ m_fullParamSpaceScan

bool DiTauMassTools::MissingMassCalculator::m_fullParamSpaceScan {}
private

Definition at line 153 of file MissingMassCalculator.h.

153{};

◆ m_fUseEfficiencyRecovery

bool DiTauMassTools::MissingMassCalculator::m_fUseEfficiencyRecovery {}
private

Definition at line 67 of file MissingMassCalculator.h.

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

◆ m_fUseFloatStopping

bool DiTauMassTools::MissingMassCalculator::m_fUseFloatStopping {}
private

Definition at line 68 of file MissingMassCalculator.h.

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

◆ m_fUseFloatStoppingCheckFreq

int DiTauMassTools::MissingMassCalculator::m_fUseFloatStoppingCheckFreq {}
private

Definition at line 70 of file MissingMassCalculator.h.

70{};

◆ m_fUseFloatStoppingComp

double DiTauMassTools::MissingMassCalculator::m_fUseFloatStoppingComp {}
private

Definition at line 71 of file MissingMassCalculator.h.

71{};

◆ m_fUseFloatStoppingMinIter

int DiTauMassTools::MissingMassCalculator::m_fUseFloatStoppingMinIter {}
private

Definition at line 69 of file MissingMassCalculator.h.

69{};

◆ m_iang1high

int DiTauMassTools::MissingMassCalculator::m_iang1high {}
private

◆ m_iang1low

int DiTauMassTools::MissingMassCalculator::m_iang1low {}
private

◆ m_iang2high

int DiTauMassTools::MissingMassCalculator::m_iang2high {}
private

◆ m_iang2low

int DiTauMassTools::MissingMassCalculator::m_iang2low {}
private

◆ m_iter0

int DiTauMassTools::MissingMassCalculator::m_iter0 {}
private

Definition at line 110 of file MissingMassCalculator.h.

110{};

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

118{};

◆ m_iterNuPV3

int DiTauMassTools::MissingMassCalculator::m_iterNuPV3 {}
private

Definition at line 111 of file MissingMassCalculator.h.

111{};

◆ m_lfvLeplepRefit

bool DiTauMassTools::MissingMassCalculator::m_lfvLeplepRefit
private

Definition at line 93 of file MissingMassCalculator.h.

◆ m_m2Nu1

double DiTauMassTools::MissingMassCalculator::m_m2Nu1 {}
private

Definition at line 183 of file MissingMassCalculator.h.

183{};

◆ m_m2Nu2

double DiTauMassTools::MissingMassCalculator::m_m2Nu2 {}
private

Definition at line 184 of file MissingMassCalculator.h.

184{};

◆ m_markovCountDuplicate

int DiTauMassTools::MissingMassCalculator::m_markovCountDuplicate {}
private

Definition at line 124 of file MissingMassCalculator.h.

124{};

◆ m_markovNAccept

int DiTauMassTools::MissingMassCalculator::m_markovNAccept {}
private

Definition at line 128 of file MissingMassCalculator.h.

128{};

◆ m_markovNFullScan

int DiTauMassTools::MissingMassCalculator::m_markovNFullScan {}
private

Definition at line 125 of file MissingMassCalculator.h.

125{};

◆ m_markovNRejectMetropolis

int DiTauMassTools::MissingMassCalculator::m_markovNRejectMetropolis {}
private

Definition at line 127 of file MissingMassCalculator.h.

127{};

◆ m_markovNRejectNoSol

int DiTauMassTools::MissingMassCalculator::m_markovNRejectNoSol {}
private

Definition at line 126 of file MissingMassCalculator.h.

126{};

◆ m_meanbinStop

double DiTauMassTools::MissingMassCalculator::m_meanbinStop {}
private

Definition at line 77 of file MissingMassCalculator.h.

77{};

◆ m_meanbinToBeEvaluated

bool DiTauMassTools::MissingMassCalculator::m_meanbinToBeEvaluated {}
private

Definition at line 122 of file MissingMassCalculator.h.

122{};

◆ m_Meff

double DiTauMassTools::MissingMassCalculator::m_Meff {}
private

Definition at line 191 of file MissingMassCalculator.h.

191{},m_Meff{};

◆ m_metCovPhiCos

double DiTauMassTools::MissingMassCalculator::m_metCovPhiCos {}
private

Definition at line 148 of file MissingMassCalculator.h.

148{},m_metCovPhiSin{};

◆ m_metCovPhiSin

double DiTauMassTools::MissingMassCalculator::m_metCovPhiSin {}
private

Definition at line 148 of file MissingMassCalculator.h.

148{},m_metCovPhiSin{};

◆ m_MEtL

double DiTauMassTools::MissingMassCalculator::m_MEtL {}
private

Definition at line 139 of file MissingMassCalculator.h.

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

◆ m_MEtL0

double DiTauMassTools::MissingMassCalculator::m_MEtL0 {}
private

Definition at line 142 of file MissingMassCalculator.h.

◆ m_MEtLMax

double DiTauMassTools::MissingMassCalculator::m_MEtLMax {}
private

Definition at line 144 of file MissingMassCalculator.h.

◆ m_MEtLMin

double DiTauMassTools::MissingMassCalculator::m_MEtLMin {}
private

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

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

◆ m_MEtP0

double DiTauMassTools::MissingMassCalculator::m_MEtP0 {}
private

Definition at line 142 of file MissingMassCalculator.h.

◆ m_MEtPMax

double DiTauMassTools::MissingMassCalculator::m_MEtPMax {}
private

Definition at line 144 of file MissingMassCalculator.h.

◆ m_MEtPMin

double DiTauMassTools::MissingMassCalculator::m_MEtPMin {}
private

Definition at line 143 of file MissingMassCalculator.h.

◆ m_MEtPRange

double DiTauMassTools::MissingMassCalculator::m_MEtPRange {}
private

◆ m_MEtProposal

double DiTauMassTools::MissingMassCalculator::m_MEtProposal {}
private

Definition at line 147 of file MissingMassCalculator.h.

◆ m_MEtPStep

double DiTauMassTools::MissingMassCalculator::m_MEtPStep {}
private

Definition at line 145 of file MissingMassCalculator.h.

◆ m_mmcCalibrationSet

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

Definition at line 65 of file MissingMassCalculator.h.

65{};

◆ m_Mnu1

double DiTauMassTools::MissingMassCalculator::m_Mnu1 {}
private

Definition at line 139 of file MissingMassCalculator.h.

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

◆ m_Mnu10

double DiTauMassTools::MissingMassCalculator::m_Mnu10 {}
private

Definition at line 142 of file MissingMassCalculator.h.

◆ m_Mnu1Exclude

bool DiTauMassTools::MissingMassCalculator::m_Mnu1Exclude {}
private

Definition at line 154 of file MissingMassCalculator.h.

154{};

◆ m_Mnu1ExcludeMax

double DiTauMassTools::MissingMassCalculator::m_Mnu1ExcludeMax {}
private

◆ m_Mnu1ExcludeMin

double DiTauMassTools::MissingMassCalculator::m_Mnu1ExcludeMin {}
private

Definition at line 168 of file MissingMassCalculator.h.

◆ m_Mnu1ExcludeRange

double DiTauMassTools::MissingMassCalculator::m_Mnu1ExcludeRange {}
private

Definition at line 168 of file MissingMassCalculator.h.

◆ m_Mnu1Max

double DiTauMassTools::MissingMassCalculator::m_Mnu1Max {}
private

Definition at line 144 of file MissingMassCalculator.h.

◆ m_Mnu1Min

double DiTauMassTools::MissingMassCalculator::m_Mnu1Min {}
private

Definition at line 143 of file MissingMassCalculator.h.

◆ m_Mnu1Range

double DiTauMassTools::MissingMassCalculator::m_Mnu1Range {}
private

◆ m_Mnu1Step

double DiTauMassTools::MissingMassCalculator::m_Mnu1Step {}
private

Definition at line 145 of file MissingMassCalculator.h.

◆ m_Mnu1XMax

double DiTauMassTools::MissingMassCalculator::m_Mnu1XMax {}
private

◆ m_Mnu1XMin

double DiTauMassTools::MissingMassCalculator::m_Mnu1XMin {}
private

Definition at line 169 of file MissingMassCalculator.h.

◆ m_Mnu1XRange

double DiTauMassTools::MissingMassCalculator::m_Mnu1XRange {}
private

Definition at line 169 of file MissingMassCalculator.h.

◆ m_Mnu2

double DiTauMassTools::MissingMassCalculator::m_Mnu2 {}
private

Definition at line 139 of file MissingMassCalculator.h.

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

◆ m_Mnu20

double DiTauMassTools::MissingMassCalculator::m_Mnu20 {}
private

Definition at line 142 of file MissingMassCalculator.h.

◆ m_Mnu2Max

double DiTauMassTools::MissingMassCalculator::m_Mnu2Max {}
private

Definition at line 144 of file MissingMassCalculator.h.

◆ m_Mnu2Min

double DiTauMassTools::MissingMassCalculator::m_Mnu2Min {}
private

Definition at line 143 of file MissingMassCalculator.h.

◆ m_Mnu2Range

double DiTauMassTools::MissingMassCalculator::m_Mnu2Range {}
private

◆ m_Mnu2Step

double DiTauMassTools::MissingMassCalculator::m_Mnu2Step {}
private

Definition at line 145 of file MissingMassCalculator.h.

◆ m_MnuProposal

double DiTauMassTools::MissingMassCalculator::m_MnuProposal {}
private

Definition at line 147 of file MissingMassCalculator.h.

◆ m_MnuScanRange

double DiTauMassTools::MissingMassCalculator::m_MnuScanRange {}
private

Definition at line 256 of file MissingMassCalculator.h.

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

138{},m_mTau2{};

◆ m_mTau2

double DiTauMassTools::MissingMassCalculator::m_mTau2 {}
private

Definition at line 138 of file MissingMassCalculator.h.

138{},m_mTau2{};

◆ m_mtautauFinalSolOldVec

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

Definition at line 157 of file MissingMassCalculator.h.

◆ m_mtautauFinalSolVec

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

Definition at line 163 of file MissingMassCalculator.h.

◆ m_mtautauSum

double DiTauMassTools::MissingMassCalculator::m_mtautauSum {}
private

Definition at line 104 of file MissingMassCalculator.h.

104{};

◆ m_Mvis

double DiTauMassTools::MissingMassCalculator::m_Mvis {}
private

Definition at line 191 of file MissingMassCalculator.h.

191{},m_Meff{};

◆ m_niter_fit1

int DiTauMassTools::MissingMassCalculator::m_niter_fit1 {}
private

Definition at line 246 of file MissingMassCalculator.h.

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

◆ m_niter_fit2

int DiTauMassTools::MissingMassCalculator::m_niter_fit2 {}
private

Definition at line 247 of file MissingMassCalculator.h.

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

◆ m_niter_fit3

int DiTauMassTools::MissingMassCalculator::m_niter_fit3 {}
private

Definition at line 248 of file MissingMassCalculator.h.

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

◆ m_NiterRandom

int DiTauMassTools::MissingMassCalculator::m_NiterRandom {}
private

Definition at line 249 of file MissingMassCalculator.h.

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

◆ m_niterRandomLocal

int DiTauMassTools::MissingMassCalculator::m_niterRandomLocal {}
private

Definition at line 74 of file MissingMassCalculator.h.

74{};

◆ m_nosol1

int DiTauMassTools::MissingMassCalculator::m_nosol1 {}
private

Definition at line 116 of file MissingMassCalculator.h.

116{};

◆ m_nosol2

int DiTauMassTools::MissingMassCalculator::m_nosol2 {}
private

Definition at line 117 of file MissingMassCalculator.h.

117{};

◆ m_nsigma_METscan

double DiTauMassTools::MissingMassCalculator::m_nsigma_METscan {}
private

Definition at line 95 of file MissingMassCalculator.h.

◆ m_nsigma_METscan2

double DiTauMassTools::MissingMassCalculator::m_nsigma_METscan2 {}
private

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

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

◆ m_nsol

int DiTauMassTools::MissingMassCalculator::m_nsol {}
private

Definition at line 161 of file MissingMassCalculator.h.

161{};

◆ m_nsolfinalmax

int DiTauMassTools::MissingMassCalculator::m_nsolfinalmax {}
private

Definition at line 73 of file MissingMassCalculator.h.

73{};

◆ m_nsolmax

int DiTauMassTools::MissingMassCalculator::m_nsolmax
private

Definition at line 73 of file MissingMassCalculator.h.

◆ m_nsolOld

int DiTauMassTools::MissingMassCalculator::m_nsolOld {}
private

Definition at line 155 of file MissingMassCalculator.h.

155{};

◆ m_NsucStop

int DiTauMassTools::MissingMassCalculator::m_NsucStop {}
private

Definition at line 250 of file MissingMassCalculator.h.

250{};

◆ m_nsucStop

int DiTauMassTools::MissingMassCalculator::m_nsucStop {}
private

Definition at line 75 of file MissingMassCalculator.h.

75{};

◆ m_nu1FinalSolOldVec

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

Definition at line 158 of file MissingMassCalculator.h.

◆ m_nu1FinalSolVec

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

Definition at line 164 of file MissingMassCalculator.h.

◆ m_nu2FinalSolOldVec

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

Definition at line 159 of file MissingMassCalculator.h.

◆ m_nu2FinalSolVec

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

Definition at line 165 of file MissingMassCalculator.h.

◆ m_nuvec1_tmp

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

Definition at line 88 of file MissingMassCalculator.h.

◆ m_nuvec2_tmp

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

Definition at line 89 of file MissingMassCalculator.h.

◆ m_nuvecsol1

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

Definition at line 80 of file MissingMassCalculator.h.

◆ m_nuvecsol2

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

Definition at line 81 of file MissingMassCalculator.h.

◆ m_Phi1

double DiTauMassTools::MissingMassCalculator::m_Phi1 {}
private

Definition at line 139 of file MissingMassCalculator.h.

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

◆ m_Phi10

double DiTauMassTools::MissingMassCalculator::m_Phi10 {}
private

Definition at line 142 of file MissingMassCalculator.h.

◆ m_Phi1Max

double DiTauMassTools::MissingMassCalculator::m_Phi1Max {}
private

Definition at line 144 of file MissingMassCalculator.h.

◆ m_Phi1Min

double DiTauMassTools::MissingMassCalculator::m_Phi1Min {}
private

Definition at line 143 of file MissingMassCalculator.h.

◆ m_Phi1Range

double DiTauMassTools::MissingMassCalculator::m_Phi1Range {}
private

◆ m_Phi1Step

double DiTauMassTools::MissingMassCalculator::m_Phi1Step {}
private

Definition at line 145 of file MissingMassCalculator.h.

◆ m_Phi2

double DiTauMassTools::MissingMassCalculator::m_Phi2 {}
private

Definition at line 139 of file MissingMassCalculator.h.

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

◆ m_Phi20

double DiTauMassTools::MissingMassCalculator::m_Phi20 {}
private

Definition at line 142 of file MissingMassCalculator.h.

◆ m_Phi2Max

double DiTauMassTools::MissingMassCalculator::m_Phi2Max {}
private

Definition at line 144 of file MissingMassCalculator.h.

◆ m_Phi2Min

double DiTauMassTools::MissingMassCalculator::m_Phi2Min {}
private

Definition at line 143 of file MissingMassCalculator.h.

◆ m_Phi2Range

double DiTauMassTools::MissingMassCalculator::m_Phi2Range {}
private

◆ m_Phi2Step

double DiTauMassTools::MissingMassCalculator::m_Phi2Step {}
private

Definition at line 145 of file MissingMassCalculator.h.

◆ m_PhiProposal

double DiTauMassTools::MissingMassCalculator::m_PhiProposal {}
private

Definition at line 147 of file MissingMassCalculator.h.

◆ m_PrintmInvWidth2Error

double DiTauMassTools::MissingMassCalculator::m_PrintmInvWidth2Error {}
private

Definition at line 132 of file MissingMassCalculator.h.

132{};

◆ m_PrintmMaxError

double DiTauMassTools::MissingMassCalculator::m_PrintmMaxError {}
private

Definition at line 130 of file MissingMassCalculator.h.

130{};

◆ m_PrintmMeanError

double DiTauMassTools::MissingMassCalculator::m_PrintmMeanError {}
private

Definition at line 131 of file MissingMassCalculator.h.

131{};

◆ m_prob_tmp

double DiTauMassTools::MissingMassCalculator::m_prob_tmp {}
private

Definition at line 101 of file MissingMassCalculator.h.

101{};

◆ m_probFinalSolOldVec

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

Definition at line 156 of file MissingMassCalculator.h.

◆ m_probFinalSolVec

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

Definition at line 162 of file MissingMassCalculator.h.

◆ m_proposalTryMEt

double DiTauMassTools::MissingMassCalculator::m_proposalTryMEt {}
private

Definition at line 134 of file MissingMassCalculator.h.

134{};

◆ m_ProposalTryMnu

double DiTauMassTools::MissingMassCalculator::m_ProposalTryMnu {}
private

Definition at line 136 of file MissingMassCalculator.h.

136{};

◆ m_ProposalTryPhi

double DiTauMassTools::MissingMassCalculator::m_ProposalTryPhi {}
private

Definition at line 135 of file MissingMassCalculator.h.

135{};

◆ m_randomGen

TRandom2 DiTauMassTools::MissingMassCalculator::m_randomGen
private

Definition at line 63 of file MissingMassCalculator.h.

◆ m_RMSStop

int DiTauMassTools::MissingMassCalculator::m_RMSStop {}
private

Definition at line 251 of file MissingMassCalculator.h.

251{};

◆ m_rmsStop

int DiTauMassTools::MissingMassCalculator::m_rmsStop {}
private

Definition at line 76 of file MissingMassCalculator.h.

76{};

◆ m_RndmSeedAltering

int DiTauMassTools::MissingMassCalculator::m_RndmSeedAltering {}
private

Definition at line 252 of file MissingMassCalculator.h.

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

◆ m_SaveLlhHisto

bool DiTauMassTools::MissingMassCalculator::m_SaveLlhHisto
private

Definition at line 93 of file MissingMassCalculator.h.

◆ m_scanMnu1

bool DiTauMassTools::MissingMassCalculator::m_scanMnu1 {}
private

Definition at line 173 of file MissingMassCalculator.h.

173{},m_scanMnu2{};

◆ m_scanMnu2

bool DiTauMassTools::MissingMassCalculator::m_scanMnu2 {}
private

Definition at line 173 of file MissingMassCalculator.h.

173{},m_scanMnu2{};

◆ m_seed

int DiTauMassTools::MissingMassCalculator::m_seed {}
private

Definition at line 107 of file MissingMassCalculator.h.

107{};

◆ m_sinPhi1

double DiTauMassTools::MissingMassCalculator::m_sinPhi1 {}
private

Definition at line 171 of file MissingMassCalculator.h.

171{}, m_cosPhi2{}, m_sinPhi1{}, m_sinPhi2{};

◆ m_sinPhi2

double DiTauMassTools::MissingMassCalculator::m_sinPhi2 {}
private

Definition at line 171 of file MissingMassCalculator.h.

171{}, m_cosPhi2{}, m_sinPhi1{}, m_sinPhi2{};

◆ m_switch1

bool DiTauMassTools::MissingMassCalculator::m_switch1 {}
private

Definition at line 119 of file MissingMassCalculator.h.

119{};

◆ m_switch2

bool DiTauMassTools::MissingMassCalculator::m_switch2 {}
private

Definition at line 120 of file MissingMassCalculator.h.

120{};

◆ m_tautau_tmp

PtEtaPhiMVector DiTauMassTools::MissingMassCalculator::m_tautau_tmp
private

Definition at line 91 of file MissingMassCalculator.h.

◆ m_tauVec1

PtEtaPhiMVector DiTauMassTools::MissingMassCalculator::m_tauVec1
private

Definition at line 175 of file MissingMassCalculator.h.

◆ m_tauVec1E

double DiTauMassTools::MissingMassCalculator::m_tauVec1E {}
private

Definition at line 181 of file MissingMassCalculator.h.

181{};

◆ m_tauVec1M

double DiTauMassTools::MissingMassCalculator::m_tauVec1M {}
private

Definition at line 177 of file MissingMassCalculator.h.

177{}, m_tauVec2M{};

◆ m_tauVec1P

double DiTauMassTools::MissingMassCalculator::m_tauVec1P {}
private

Definition at line 180 of file MissingMassCalculator.h.

180{}, m_tauVec2P{};

◆ m_tauVec1Phi

double DiTauMassTools::MissingMassCalculator::m_tauVec1Phi {}
private

Definition at line 176 of file MissingMassCalculator.h.

176{}, m_tauVec2Phi{};

◆ m_tauVec1Px

double DiTauMassTools::MissingMassCalculator::m_tauVec1Px {}
private

Definition at line 178 of file MissingMassCalculator.h.

178{}, m_tauVec1Py{}, m_tauVec1Pz{};

◆ m_tauVec1Py

double DiTauMassTools::MissingMassCalculator::m_tauVec1Py {}
private

Definition at line 178 of file MissingMassCalculator.h.

178{}, m_tauVec1Py{}, m_tauVec1Pz{};

◆ m_tauVec1Pz

double DiTauMassTools::MissingMassCalculator::m_tauVec1Pz {}
private

Definition at line 178 of file MissingMassCalculator.h.

178{}, m_tauVec1Py{}, m_tauVec1Pz{};

◆ m_tauVec2

PtEtaPhiMVector DiTauMassTools::MissingMassCalculator::m_tauVec2
private

Definition at line 175 of file MissingMassCalculator.h.

◆ m_tauVec2E

double DiTauMassTools::MissingMassCalculator::m_tauVec2E {}
private

Definition at line 182 of file MissingMassCalculator.h.

182{};

◆ m_tauVec2M

double DiTauMassTools::MissingMassCalculator::m_tauVec2M {}
private

Definition at line 177 of file MissingMassCalculator.h.

177{}, m_tauVec2M{};

◆ m_tauVec2P

double DiTauMassTools::MissingMassCalculator::m_tauVec2P {}
private

Definition at line 180 of file MissingMassCalculator.h.

180{}, m_tauVec2P{};

◆ m_tauVec2Phi

double DiTauMassTools::MissingMassCalculator::m_tauVec2Phi {}
private

Definition at line 176 of file MissingMassCalculator.h.

176{}, m_tauVec2Phi{};

◆ m_tauVec2Px

double DiTauMassTools::MissingMassCalculator::m_tauVec2Px {}
private

Definition at line 179 of file MissingMassCalculator.h.

179{}, m_tauVec2Py{}, m_tauVec2Pz{};

◆ m_tauVec2Py

double DiTauMassTools::MissingMassCalculator::m_tauVec2Py {}
private

Definition at line 179 of file MissingMassCalculator.h.

179{}, m_tauVec2Py{}, m_tauVec2Pz{};

◆ m_tauVec2Pz

double DiTauMassTools::MissingMassCalculator::m_tauVec2Pz {}
private

Definition at line 179 of file MissingMassCalculator.h.

179{}, m_tauVec2Py{}, m_tauVec2Pz{};

◆ m_tauvecprob1

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

Definition at line 85 of file MissingMassCalculator.h.

◆ m_tauvecprob2

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

Definition at line 86 of file MissingMassCalculator.h.

◆ m_tauvecsol1

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

Definition at line 83 of file MissingMassCalculator.h.

◆ m_tauvecsol2

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

Definition at line 84 of file MissingMassCalculator.h.

◆ m_testdiscri1

int DiTauMassTools::MissingMassCalculator::m_testdiscri1 {}
private

Definition at line 114 of file MissingMassCalculator.h.

114{};

◆ m_testdiscri2

int DiTauMassTools::MissingMassCalculator::m_testdiscri2 {}
private

Definition at line 115 of file MissingMassCalculator.h.

115{};

◆ m_testptn1

int DiTauMassTools::MissingMassCalculator::m_testptn1 {}
private

Definition at line 112 of file MissingMassCalculator.h.

112{};

◆ m_testptn2

int DiTauMassTools::MissingMassCalculator::m_testptn2 {}
private

Definition at line 113 of file MissingMassCalculator.h.

113{};

◆ m_TLVdummy

PtEtaPhiMVector DiTauMassTools::MissingMassCalculator::m_TLVdummy
private

Definition at line 240 of file MissingMassCalculator.h.

◆ m_totalProbSum

double DiTauMassTools::MissingMassCalculator::m_totalProbSum {}
private

Definition at line 103 of file MissingMassCalculator.h.

103{};

◆ m_walkWeight

double DiTauMassTools::MissingMassCalculator::m_walkWeight {}
private

Definition at line 170 of file MissingMassCalculator.h.

170{};

◆ metvec_tmp

XYVector DiTauMassTools::MissingMassCalculator::metvec_tmp

Definition at line 409 of file MissingMassCalculator.h.

◆ OutputInfo

MissingMassOutput DiTauMassTools::MissingMassCalculator::OutputInfo

Definition at line 333 of file MissingMassCalculator.h.

◆ preparedInput

MissingMassInput DiTauMassTools::MissingMassCalculator::preparedInput

Definition at line 332 of file MissingMassCalculator.h.

◆ Prob

MissingMassProb* DiTauMassTools::MissingMassCalculator::Prob

Definition at line 334 of file MissingMassCalculator.h.


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