ATLAS Offline Software
Loading...
Searching...
No Matches
MissingMassCalculator.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
3*/
4
5// vim: ts=8 sw=2
6/*
7 Missing Mass Calculator
8*/
9//
10// to be done : tau 4-vect and type should be data member of MMC.
11
12// if histogram smoothing
13//#define SMOOTH
14
15#include "DiTauMassTools/MissingMassCalculator.h" // this is for RootCore package
16#include "DiTauMassTools/MissingMassProb.h"
17#include "xAODMissingET/MissingET.h"
18#include "xAODTau/TauJet.h"
19
20#include <TObject.h>
21// SpeedUp committed from revision 163876
22#include <TF1.h>
23#include <TFile.h>
24#include <TFitResult.h>
25#include <TFitResultPtr.h>
26#include <TMatrixDSym.h>
27#include <TMatrixT.h>
28#include <TObject.h>
29#include <TVectorD.h>
30#include "Math/VectorUtil.h"
31
32#include "TruthUtils/ParticleConstants.h"
33
34#include <fstream>
35#include <iomanip>
36#include <iostream>
37#include <sstream>
38#include <stdexcept>
39namespace {
40 constexpr double GEV = 1000.0;
41}
42
43
44using namespace DiTauMassTools;
45using ROOT::Math::PtEtaPhiMVector;
46using ROOT::Math::PxPyPzMVector;
47using ROOT::Math::XYVector;
48using ROOT::Math::VectorUtil::DeltaR;
49using ROOT::Math::VectorUtil::Phi_mpi_pi;
50
51//______________________________constructor________________________________
52MissingMassCalculator::MissingMassCalculator(
53 MMCCalibrationSet::e aset, std::string paramFilePath)
54 : m_randomGen(), Prob(new MissingMassProb(aset, paramFilePath)) {
55 m_mmcCalibrationSet = aset;
56 preparedInput.m_fUseVerbose = 0;
57 preparedInput.m_beamEnergy = 6500.0; // for now LHC default is sqrt(S)=7 TeV
58 m_niter_fit1 = 20;
59 m_niter_fit2 = 30;
60 m_niter_fit3 = 10;
61 m_NsucStop = -1;
62 m_NiterRandom = -1; // if the user does not set it to positive value,will be set
63 // in SpaceWalkerInit
64 m_niterRandomLocal = -1; // niterandom which is really used
65 // to be used with RMSSTOP NiterRandom=10000000; // number of random
66 // iterations for lh. Multiplied by 10 for ll, divided by 10 for hh (to be
67 // optimised)
68 // RMSStop=200;// Stop criteria depending of rms of histogram
69 m_RMSStop = -1; // disable
70
71 m_RndmSeedAltering = 0; // can be changed to re-compute with different random seed
72 m_dRmax_tau = 0.4; // changed from 0.2
73 m_nsigma_METscan = -1; // number of sigmas for MET scan
74 m_nsigma_METscan_ll = 3.0; // number of sigmas for MET scan
75 m_nsigma_METscan_lh = 3.0; // number of sigmas for MET scan
76 m_nsigma_METscan_hh = 4.0; // number of sigmas for MET scan (4 for hh 2013)
77 m_nsigma_METscan_lfv_ll = 5.0; // number of sigmas for MET scan (LFV leplep)
78 m_nsigma_METscan_lfv_lh = 5.0; // number of sigmas for MET scan (LFV lephad)
79
80 m_meanbinStop = -1; // meanbin stopping criterion (-1 if not used)
81 m_proposalTryMEt = -1; // loop on METproposal disable // FIXME should be cleaner
82 m_ProposalTryPhi = -1; // loop on Phiproposal disable
83 m_ProposalTryMnu = -1; // loop on MNuProposal disable
84
85 Prob->SetUseTauProbability(true); // TauProbability is ON by default DRMERGE comment out for now
86 Prob->SetUseMnuProbability(false); // MnuProbability is OFF by default
87 Prob->SetUseDphiLL(false); // added by Tomas Davidek for lep-lep
88 preparedInput.m_METresSyst = 0; // no MET resolution systematics by default (+/-1: up/down 1 sigma)
89 preparedInput.m_dataType = 1; // set to "data" by default
90 preparedInput.m_fUseTailCleanup = 1; // cleanup by default for lep-had Moriond 2012 analysis
91 preparedInput.m_fUseDefaults = 0; // use pre-set defaults for various configurations; if set it to 0
92 // if need to study various options
93 m_fUseEfficiencyRecovery = 0; // no re-fit by default
94 m_fUseFloatStopping = 0;
95 m_fUseFloatStoppingMinIter = 10000;
96 m_fUseFloatStoppingCheckFreq = 1000;
97 m_fUseFloatStoppingComp = 0.05;
98
99 preparedInput.m_METScanScheme = 1; // MET-scan scheme: 0- use JER; 1- use simple sumEt & missingHt
100 // for Njet=0 events in (lep-had winter 2012)
101 // MnuScanRange=ParticleConstants::tauMassInMeV / GEV; // range of M(nunu) scan
102 m_MnuScanRange = 1.5; // better value (sacha)
103 preparedInput.m_LFVmode = -1; // by default consider case of H->mu+tau(->ele)
104 preparedInput.ClearInput();
105
106 m_debugThisIteration = false;
107 m_lfvLeplepRefit = true;
108 m_SaveLlhHisto = false;
109
110 m_nsolmax = 4;
111 m_nsolfinalmax = m_nsolmax * m_nsolmax;
112
113 m_nuvecsol1.resize(m_nsolmax);
114 m_nuvecsol2.resize(m_nsolmax);
115 m_tauvecsol1.resize(m_nsolmax);
116 m_tauvecsol2.resize(m_nsolmax);
117 m_tauvecprob1.resize(m_nsolmax);
118 m_tauvecprob2.resize(m_nsolmax);
119
120 m_nsol = 0;
121 m_probFinalSolVec.resize(m_nsolfinalmax);
122 m_mtautauFinalSolVec.resize(m_nsolfinalmax);
123 m_nu1FinalSolVec.resize(m_nsolfinalmax);
124 m_nu2FinalSolVec.resize(m_nsolfinalmax);
125
126 m_nsolOld = 0;
127 m_probFinalSolOldVec.resize(m_nsolfinalmax);
128 m_mtautauFinalSolOldVec.resize(m_nsolfinalmax);
129 m_nu1FinalSolOldVec.resize(m_nsolfinalmax);
130 m_nu2FinalSolOldVec.resize(m_nsolfinalmax);
131
132 float hEmax = 3000.0; // maximum energy (GeV)
133 // number of bins
134 int hNbins = 1500; // original 2500 for mass, 10000 for P
135 // choice of hNbins also related to size of window for fitting (see
136 // maxFromHist)
137
138 //--- define histograms for histogram method
139 //--- upper limits need to be revisied in the future!!! It may be not enough
140 // for some analyses
141
142 m_fMfit_all = std::make_shared<TH1F>("MMC_h1", "M", hNbins, 0.0,
143 hEmax); // all solutions
144 m_fMfit_all->Sumw2(); // allow proper error bin calculation. Slightly slower but
145 // completely negligible
146
147 // histogram without weight. useful for debugging. negligibly slow until now
148 m_fMfit_allNoWeight =
149 std::make_shared<TH1F>("MMC_h1NoW", "M no weight", hNbins, 0.0, hEmax); // all solutions
150
151 m_fPXfit1 = std::make_shared<TH1F>("MMC_h2", "Px1", 4 * hNbins, -hEmax,
152 hEmax); // Px for tau1
153 m_fPYfit1 = std::make_shared<TH1F>("MMC_h3", "Py1", 4 * hNbins, -hEmax,
154 hEmax); // Py for tau1
155 m_fPZfit1 = std::make_shared<TH1F>("MMC_h4", "Pz1", 4 * hNbins, -hEmax,
156 hEmax); // Pz for tau1
157 m_fPXfit2 = std::make_shared<TH1F>("MMC_h5", "Px2", 4 * hNbins, -hEmax,
158 hEmax); // Px for tau2
159 m_fPYfit2 = std::make_shared<TH1F>("MMC_h6", "Py2", 4 * hNbins, -hEmax,
160 hEmax); // Py for tau2
161 m_fPZfit2 = std::make_shared<TH1F>("MMC_h7", "Pz2", 4 * hNbins, -hEmax,
162 hEmax); // Pz for tau2
163
164 m_fMfit_all->SetDirectory(0);
165
166 m_fMfit_allNoWeight->SetDirectory(0);
167 m_fPXfit1->SetDirectory(0);
168 m_fPYfit1->SetDirectory(0);
169 m_fPZfit1->SetDirectory(0);
170 m_fPXfit2->SetDirectory(0);
171 m_fPYfit2->SetDirectory(0);
172 m_fPZfit2->SetDirectory(0);
173
174 // max hist fitting function
175 m_fFitting =
176 new TF1("MMC_maxFitting", this, &MissingMassCalculator::maxFitting, 0., hEmax, 3);
177 // Sets initial parameter names
178 m_fFitting->SetParNames("Max", "Mean", "InvWidth2");
179
180 if (preparedInput.m_fUseVerbose == 1) {
181 gDirectory->pwd();
182 gDirectory->ls();
183 }
184
185 if (preparedInput.m_fUseVerbose == 1) {
186 gDirectory->pwd();
187 gDirectory->ls();
188 }
189}
190
191MissingMassCalculator::~MissingMassCalculator() { delete Prob; }
192
193//_____________________________________________________________________________
194// Main Method to run MissingMassCalculator
195int MissingMassCalculator::RunMissingMassCalculator(const xAOD::IParticle *part1,
196 const xAOD::IParticle *part2,
197 const xAOD::MissingET *met,
198 const int &njets) {
199
200 OutputInfo.ClearOutput(preparedInput.m_fUseVerbose);
201 if (preparedInput.m_fUseVerbose == 1) {
202 Info("DiTauMassTools", "------------- Raw Input for MissingMassCalculator --------------");
203 }
204 FinalizeSettings(part1, part2, met, njets); // rawInput, preparedInput );
205 Prob->MET(preparedInput);
206 if (preparedInput.m_fUseVerbose == 1) {
207 Info("DiTauMassTools", "------------- Prepared Input for MissingMassCalculator--------------");
208 preparedInput.PrintInputInfo();
209 }
210
211 if (preparedInput.m_LFVmode < 0) {
212 // remove argument DiTauMassCalculatorV9Walk work directly on preparedInput
213 OutputInfo.m_FitStatus = DitauMassCalculatorV9walk();
214
215 // re-running MMC for on failed events
216 if (m_fUseEfficiencyRecovery == 1 && OutputInfo.m_FitStatus != 1) {
217 // most events where MMC failed happened to have dPhi>2.9. Run re-fit only
218 // on these events
219 if (preparedInput.m_DelPhiTT > 2.9) {
220 // preparedInput.MetVec.Set(-(preparedInput.vistau1+preparedInput.vistau2).Px(),-(preparedInput.vistau1+preparedInput.vistau2).Py());
221 // // replace MET by MPT
222
223 XYVector dummy_met(-(preparedInput.m_vistau1 + preparedInput.m_vistau2).Px(),
224 -(preparedInput.m_vistau1 + preparedInput.m_vistau2).Py());
225 preparedInput.m_METcovphi = dummy_met.Phi();
226 double dummy_METres =
227 sqrt(pow(preparedInput.m_METsigmaL, 2) + pow(preparedInput.m_METsigmaP, 2));
228 preparedInput.m_METsigmaL =
229 dummy_METres * std::abs(cos(dummy_met.Phi() - preparedInput.m_MetVec.Phi()));
230 preparedInput.m_METsigmaP =
231 dummy_METres * std::abs(sin(dummy_met.Phi() - preparedInput.m_MetVec.Phi()));
232 if (preparedInput.m_METsigmaP < 5.0)
233 preparedInput.m_METsigmaP = 5.0;
234 m_nsigma_METscan_lh = 6.0; // increase range of MET scan
235 m_nsigma_METscan_hh = 6.0; // increase range of MET scan
236
237 OutputInfo.ClearOutput(preparedInput.m_fUseVerbose); // clear output stuff before re-running
238 OutputInfo.m_FitStatus = DitauMassCalculatorV9walk(); // run MMC again
239 }
240 }
241
242 }
243
244 // running MMC in LFV mode for reconstructing mass of X->lep+tau
245 else {
246 if (preparedInput.m_fUseVerbose == 1) {
247 Info("DiTauMassTools", "Calling DitauMassCalculatorV9lfv");
248 }
249 OutputInfo.m_FitStatus = DitauMassCalculatorV9lfv(false);
250 }
251
252 if(m_SaveLlhHisto){
253 TFile *outFile = TFile::Open("MMC_likelihoods.root", "UPDATE");
254 outFile->cd();
255 auto path = std::to_string(m_eventNumber);
256 if (!outFile->GetDirectory(path.c_str()))
257 outFile->mkdir(path.c_str());
258 outFile->cd(path.c_str());
259 m_fMfit_all->Write(m_fMfit_all->GetName(), TObject::kOverwrite);
260 m_fMEtP_all->Write(m_fMEtP_all->GetName(), TObject::kOverwrite);
261 m_fMEtL_all->Write(m_fMEtL_all->GetName(), TObject::kOverwrite);
262 m_fMnu1_all->Write(m_fMnu1_all->GetName(), TObject::kOverwrite);
263 m_fMnu2_all->Write(m_fMnu2_all->GetName(), TObject::kOverwrite);
264 m_fPhi1_all->Write(m_fPhi1_all->GetName(), TObject::kOverwrite);
265 m_fPhi2_all->Write(m_fPhi2_all->GetName(), TObject::kOverwrite);
266 m_fMfit_allNoWeight->Write(m_fMfit_allNoWeight->GetName(), TObject::kOverwrite);
267 m_fMfit_allGraph->Write("Graph", TObject::kOverwrite);
268 TH1D *nosol = new TH1D("nosol", "nosol", 7, 0, 7);
269 nosol->SetBinContent(1, m_testptn1);
270 nosol->SetBinContent(2, m_testptn2);
271 nosol->SetBinContent(3, m_testdiscri1);
272 nosol->SetBinContent(4, m_testdiscri2);
273 nosol->SetBinContent(5, m_nosol1);
274 nosol->SetBinContent(6, m_nosol1);
275 nosol->SetBinContent(7, m_iterNuPV3);
276 nosol->Write(nosol->GetName(), TObject::kOverwrite);
277 outFile->Write();
278 outFile->Close();
279 }
280
281 DoOutputInfo();
282 PrintResults();
283 preparedInput.ClearInput();
284 return 1;
285}
286
287//-------- clearing ditau container
288void MissingMassCalculator::ClearDitauStuff(DitauStuff &fStuff) {
289 fStuff.Mditau_best = 0.0;
290 fStuff.Sign_best = 1.0E6;
291 fStuff.nutau1 = PtEtaPhiMVector(0., 0., 0., 0.);
292 fStuff.nutau2 = PtEtaPhiMVector(0., 0., 0., 0.);
293 fStuff.vistau1 = PtEtaPhiMVector(0., 0., 0., 0.);
294 fStuff.vistau2 = PtEtaPhiMVector(0., 0., 0., 0.);
295 fStuff.RMSoverMPV = 0.0;
296
297 return;
298}
299
300//---------------------------- Accessors to output parameters
301//------------------------
302// finalizes output information
303void MissingMassCalculator::DoOutputInfo() {
304 if (OutputInfo.m_FitStatus > 0) {
305 if (preparedInput.m_fUseVerbose == 1) {
306 Info("DiTauMassTools", "Retrieving output from fDitauStuffFit");
307 }
308 // MAXW method : get from fDittauStuffFit
309 OutputInfo.m_FitSignificance[MMCFitMethod::MAXW] = m_fDitauStuffFit.Sign_best;
310 OutputInfo.m_FittedMass[MMCFitMethod::MAXW] = m_fDitauStuffFit.Mditau_best;
311 double q1 = (1. - 0.68) / 2.;
312 double q2 = 1. - q1;
313 double xq[2], yq[2];
314 xq[0] = q1;
315 xq[1] = q2;
316 m_fMfit_all->GetQuantiles(2, yq, xq);
317 OutputInfo.m_FittedMassLowerError[MMCFitMethod::MAXW] = yq[0];
318 OutputInfo.m_FittedMassUpperError[MMCFitMethod::MAXW] = yq[1];
319 OutputInfo.m_nuvec1[MMCFitMethod::MAXW] = m_fDitauStuffFit.nutau1;
320 OutputInfo.m_objvec1[MMCFitMethod::MAXW] =
321 m_fDitauStuffFit.vistau1 + m_fDitauStuffFit.nutau1;
322 OutputInfo.m_nuvec2[MMCFitMethod::MAXW] = m_fDitauStuffFit.nutau2;
323 OutputInfo.m_objvec2[MMCFitMethod::MAXW] =
324 m_fDitauStuffFit.vistau2 + m_fDitauStuffFit.nutau2;
325 OutputInfo.m_totalvec[MMCFitMethod::MAXW] =
326 OutputInfo.m_objvec1[MMCFitMethod::MAXW] +
327 OutputInfo.m_objvec2[MMCFitMethod::MAXW];
328 XYVector metmaxw(OutputInfo.m_nuvec1[MMCFitMethod::MAXW].Px() +
329 OutputInfo.m_nuvec2[MMCFitMethod::MAXW].Px(),
330 OutputInfo.m_nuvec1[MMCFitMethod::MAXW].Py() +
331 OutputInfo.m_nuvec2[MMCFitMethod::MAXW].Py());
332 OutputInfo.m_FittedMetVec[MMCFitMethod::MAXW] = metmaxw;
333
334 OutputInfo.m_FittedMass[MMCFitMethod::MLM] = m_fDitauStuffHisto.Mditau_best;
335 OutputInfo.m_FittedMassLowerError[MMCFitMethod::MLM] = yq[0];
336 OutputInfo.m_FittedMassUpperError[MMCFitMethod::MLM] = yq[1];
337
338 PtEtaPhiMVector tlvdummy(0., 0., 0., 0.);
339 XYVector metdummy(0., 0.);
340 OutputInfo.m_FitSignificance[MMCFitMethod::MLM] = -1.;
341 OutputInfo.m_nuvec1[MMCFitMethod::MLM] = tlvdummy;
342 OutputInfo.m_objvec1[MMCFitMethod::MLM] = tlvdummy;
343 OutputInfo.m_nuvec2[MMCFitMethod::MLM] = tlvdummy;
344 OutputInfo.m_objvec2[MMCFitMethod::MLM] = tlvdummy;
345 OutputInfo.m_totalvec[MMCFitMethod::MLM] = tlvdummy;
346 OutputInfo.m_FittedMetVec[MMCFitMethod::MLM] = metdummy;
347
348 // MLNU3P method : get from fDittauStuffHisto 4 momentum
349 OutputInfo.m_nuvec1[MMCFitMethod::MLNU3P] = m_fDitauStuffHisto.nutau1;
350 OutputInfo.m_objvec1[MMCFitMethod::MLNU3P] =
351 m_fDitauStuffHisto.vistau1 + m_fDitauStuffHisto.nutau1;
352 OutputInfo.m_nuvec2[MMCFitMethod::MLNU3P] = m_fDitauStuffHisto.nutau2;
353 OutputInfo.m_objvec2[MMCFitMethod::MLNU3P] =
354 m_fDitauStuffHisto.vistau2 + m_fDitauStuffHisto.nutau2;
355 OutputInfo.m_totalvec[MMCFitMethod::MLNU3P] =
356 OutputInfo.m_objvec1[MMCFitMethod::MLNU3P] +
357 OutputInfo.m_objvec2[MMCFitMethod::MLNU3P];
358 OutputInfo.m_FittedMass[MMCFitMethod::MLNU3P] =
359 OutputInfo.m_totalvec[MMCFitMethod::MLNU3P].M();
360 OutputInfo.m_FittedMassUpperError[MMCFitMethod::MLNU3P] = 0.;
361 OutputInfo.m_FittedMassLowerError[MMCFitMethod::MLNU3P] = 0.;
362
363 XYVector metmlnu3p(OutputInfo.m_nuvec1[MMCFitMethod::MLNU3P].Px() +
364 OutputInfo.m_nuvec2[MMCFitMethod::MLNU3P].Px(),
365 OutputInfo.m_nuvec1[MMCFitMethod::MLNU3P].Py() +
366 OutputInfo.m_nuvec2[MMCFitMethod::MLNU3P].Py());
367 OutputInfo.m_FittedMetVec[MMCFitMethod::MLNU3P] = metmlnu3p;
368
369 OutputInfo.m_RMS2MPV = m_fDitauStuffHisto.RMSoverMPV;
370 }
371
372 OutputInfo.m_hMfit_all = m_fMfit_all;
373 OutputInfo.m_hMfit_allNoWeight = m_fMfit_allNoWeight;
374 OutputInfo.m_NSolutions = m_fMfit_all->GetEntries();
375 OutputInfo.m_SumW = m_fMfit_all->GetSumOfWeights();
376
377 //----------------- Check if input was re-ordered in FinalizeInputStuff() and
378 // restore the original order if needed
379 if (preparedInput.m_InputReorder == 1) {
380 PtEtaPhiMVector dummy_vec1(0.0, 0.0, 0.0, 0.0);
381 PtEtaPhiMVector dummy_vec2(0.0, 0.0, 0.0, 0.0);
382 for (int i = 0; i < 3; i++) {
383 // re-ordering neutrinos
384 dummy_vec1 = OutputInfo.m_nuvec1[i];
385 dummy_vec2 = OutputInfo.m_nuvec2[i];
386 OutputInfo.m_nuvec1[i] = dummy_vec2;
387 OutputInfo.m_nuvec2[i] = dummy_vec1;
388 // re-ordering tau's
389 dummy_vec1 = OutputInfo.m_objvec1[i];
390 dummy_vec2 = OutputInfo.m_objvec2[i];
391 OutputInfo.m_objvec1[i] = dummy_vec2;
392 OutputInfo.m_objvec2[i] = dummy_vec1;
393 }
394 }
395
396 return;
397}
398
399// Printout of final results
400void MissingMassCalculator::PrintOtherInput() {
401 if (preparedInput.m_fUseVerbose != 1)
402 return;
403
404 Info("DiTauMassTools",
405 ".........................Other input.....................................");
406 Info("DiTauMassTools", "%s",
407 ("Beam energy =" + std::to_string(preparedInput.m_beamEnergy) +
408 " sqrt(S) for collisions =" + std::to_string(2.0 * preparedInput.m_beamEnergy))
409 .c_str());
410 Info("DiTauMassTools", "%s",
411 ("CalibrationSet " + MMCCalibrationSet::name[m_mmcCalibrationSet])
412 .c_str());
413 Info("DiTauMassTools", "%s",
414 ("LFV mode " + std::to_string(preparedInput.m_LFVmode) + " seed=" + std::to_string(m_seed))
415 .c_str());
416 Info("DiTauMassTools", "%s", ("usetauProbability=" + std::to_string(Prob->GetUseTauProbability()) +
417 " useTailCleanup=" + std::to_string(preparedInput.m_fUseTailCleanup))
418 .c_str());
419
420 if (preparedInput.m_InputReorder != 0) {
421 Info("DiTauMassTools",
422 "tau1 and tau2 were internally swapped (visible on prepared input printout)");
423 } else {
424 Info("DiTauMassTools", "tau1 and tau2 were NOT internally swapped");
425 }
426
427 Info("DiTauMassTools", "%s",
428 (" MEtLMin=" + std::to_string(m_MEtLMin) + " MEtLMax=" + std::to_string(m_MEtLMax)).c_str());
429 Info("DiTauMassTools", "%s",
430 (" MEtPMin=" + std::to_string(m_MEtPMin) + " MEtPMax=" + std::to_string(m_MEtPMax)).c_str());
431 Info("DiTauMassTools", "%s",
432 (" Phi1Min=" + std::to_string(m_Phi1Min) + " Phi1Max=" + std::to_string(m_Phi1Max)).c_str());
433 Info("DiTauMassTools", "%s",
434 (" Phi2Min=" + std::to_string(m_Phi2Min) + " Phi2Max=" + std::to_string(m_Phi2Max)).c_str());
435 Info("DiTauMassTools", "%s",
436 (" Mnu1Min=" + std::to_string(m_Mnu1Min) + " Mnu1Max=" + std::to_string(m_Mnu1Max)).c_str());
437 Info("DiTauMassTools", "%s",
438 (" Mnu2Min=" + std::to_string(m_Mnu2Min) + " Mnu2Max=" + std::to_string(m_Mnu2Max)).c_str());
439}
440
441// Printout of final results
442void MissingMassCalculator::PrintResults() {
443
444 if (preparedInput.m_fUseVerbose != 1)
445 return;
446
447 const PtEtaPhiMVector *origVisTau1 = 0;
448 const PtEtaPhiMVector *origVisTau2 = 0;
449
450 if (preparedInput.m_InputReorder == 0) {
451 origVisTau1 = &preparedInput.m_vistau1;
452 origVisTau2 = &preparedInput.m_vistau2;
453 } else // input order was flipped
454 {
455 origVisTau1 = &preparedInput.m_vistau2;
456 origVisTau2 = &preparedInput.m_vistau1;
457 }
458
459 PrintOtherInput();
460
461 Info("DiTauMassTools",
462 "------------- Printing Final Results for MissingMassCalculator --------------");
463 Info("DiTauMassTools",
464 ".............................................................................");
465 Info("DiTauMassTools", "%s", ("Fit status=" + std::to_string(OutputInfo.m_FitStatus)).c_str());
466
467 for (int imeth = 0; imeth < MMCFitMethod::MAX; ++imeth) {
468 Info("DiTauMassTools", "%s",
469 ("___ Results for " + MMCFitMethod::name[imeth] + "Method ___")
470 .c_str());
471 Info("DiTauMassTools", "%s",
472 (" signif=" + std::to_string(OutputInfo.m_FitSignificance[imeth])).c_str());
473 Info("DiTauMassTools", "%s", (" mass=" + std::to_string(OutputInfo.m_FittedMass[imeth])).c_str());
474 Info("DiTauMassTools", "%s", (" rms/mpv=" + std::to_string(OutputInfo.m_RMS2MPV)).c_str());
475
476 if (imeth == MMCFitMethod::MLM) {
477 Info("DiTauMassTools", " no 4-momentum or MET from this method ");
478 continue;
479 }
480
481 if (OutputInfo.m_FitStatus <= 0) {
482 Info("DiTauMassTools", " fit failed ");
483 }
484
485 const PtEtaPhiMVector &tlvnu1 = OutputInfo.m_nuvec1[imeth];
486 const PtEtaPhiMVector &tlvnu2 = OutputInfo.m_nuvec2[imeth];
487 const PtEtaPhiMVector &tlvo1 = OutputInfo.m_objvec1[imeth];
488 const PtEtaPhiMVector &tlvo2 = OutputInfo.m_objvec2[imeth];
489 const XYVector &tvmet = OutputInfo.m_FittedMetVec[imeth];
490
491 Info("DiTauMassTools", "%s",
492 (" Neutrino-1: P=" + std::to_string(tlvnu1.P()) + " Pt=" + std::to_string(tlvnu1.Pt()) +
493 " Eta=" + std::to_string(tlvnu1.Eta()) + " Phi=" + std::to_string(tlvnu1.Phi()) +
494 " M=" + std::to_string(tlvnu1.M()) + " Px=" + std::to_string(tlvnu1.Px()) +
495 " Py=" + std::to_string(tlvnu1.Py()) + " Pz=" + std::to_string(tlvnu1.Pz()))
496 .c_str());
497 Info("DiTauMassTools", "%s",
498 (" Neutrino-2: P=" + std::to_string(tlvnu2.P()) + " Pt=" + std::to_string(tlvnu2.Pt()) +
499 " Eta=" + std::to_string(tlvnu2.Eta()) + " Phi=" + std::to_string(tlvnu2.Phi()) +
500 " M=" + std::to_string(tlvnu2.M()) + " Px=" + std::to_string(tlvnu2.Px()) +
501 " Py=" + std::to_string(tlvnu2.Py()) + " Pz=" + std::to_string(tlvnu2.Pz()))
502 .c_str());
503 Info("DiTauMassTools", "%s",
504 (" Tau-1: P=" + std::to_string(tlvo1.P()) + " Pt=" + std::to_string(tlvo1.Pt()) +
505 " Eta=" + std::to_string(tlvo1.Eta()) + " Phi=" + std::to_string(tlvo1.Phi()) +
506 " M=" + std::to_string(tlvo1.M()) + " Px=" + std::to_string(tlvo1.Px()) +
507 " Py=" + std::to_string(tlvo1.Py()) + " Pz=" + std::to_string(tlvo1.Pz()))
508 .c_str());
509 Info("DiTauMassTools", "%s",
510 (" Tau-2: P=" + std::to_string(tlvo2.P()) + " Pt=" + std::to_string(tlvo2.Pt()) +
511 " Eta=" + std::to_string(tlvo2.Eta()) + " Phi=" + std::to_string(tlvo2.Phi()) +
512 " M=" + std::to_string(tlvo2.M()) + " Px=" + std::to_string(tlvo2.Px()) +
513 " Py=" + std::to_string(tlvo2.Py()) + " Pz=" + std::to_string(tlvo2.Pz()))
514 .c_str());
515
516 Info("DiTauMassTools", "%s",
517 (" dR(nu1-visTau1)=" + std::to_string(DeltaR(tlvnu1,*origVisTau1))).c_str());
518 Info("DiTauMassTools", "%s",
519 (" dR(nu2-visTau2)=" + std::to_string(DeltaR(tlvnu2,*origVisTau2))).c_str());
520
521 Info("DiTauMassTools", "%s",
522 (" Fitted MET =" + std::to_string(tvmet.R()) + " Phi=" + std::to_string(tlvnu1.Phi()) +
523 " Px=" + std::to_string(tvmet.X()) + " Py=" + std::to_string(tvmet.Y()))
524 .c_str());
525
526 Info("DiTauMassTools", "%s", (" Resonance: P=" + std::to_string(OutputInfo.m_totalvec[imeth].P()) +
527 " Pt=" + std::to_string(OutputInfo.m_totalvec[imeth].Pt()) +
528 " Eta=" + std::to_string(OutputInfo.m_totalvec[imeth].Eta()) +
529 " Phi=" + std::to_string(OutputInfo.m_totalvec[imeth].Phi()) +
530 " M=" + std::to_string(OutputInfo.m_totalvec[imeth].M()) +
531 " Px=" + std::to_string(OutputInfo.m_totalvec[imeth].Px()) +
532 " Py=" + std::to_string(OutputInfo.m_totalvec[imeth].Py()) +
533 " Pz=" + std::to_string(OutputInfo.m_totalvec[imeth].Pz()))
534 .c_str());
535 }
536
537 return;
538}
539
540// returns P1, P2, and theta1 & theta2 solutions
541// This compute the nu1 nu2 solution in the most efficient way. Wrt to
542// NuPsolutionV2, the output nu1 nu2 4-vector have non zero mass (if relevant).
543// It is not optimised for grid running so much less caching is done (which
544// makes it more readable). Only quantities fixed within an event are cached.
545// relies on a number of these variables to be initialised before the loop.
546
547int MissingMassCalculator::NuPsolutionV3(const double &mNu1, const double &mNu2,
548 const double &phi1, const double &phi2,
549 int &nsol1, int &nsol2) {
550
551 // Pv1, Pv2 : visible tau decay product momentum
552 // Pn1 Pn2 : neutrino momentum
553 // phi1, phi2 : neutrino azymutal angles
554 // PTmiss2=PTmissy Cos[phi2] - PTmissx Sin[phi2]
555 // PTmiss2cscdphi=PTmiss2/Sin[phi1-phi2]
556 // Pv1proj=Pv1x Cos[phi1] + Pv1y Sin[phi1]
557 // M2noma1=Mtau^2-Mv1^2-Mn1^2
558 // ETv1^2=Ev1^2-Pv1z^2
559
560 // discriminant : 16 Ev1^2 (M2noma1^2 + 4 M2noma1 PTmiss2cscdphi Pv1proj - 4
561 // (ETv1^2 (Mn1^2 + PTmiss2cscdphi^2) - PTmiss2cscdphi^2 Pv1proj^2))
562 // two solutions for epsilon = +/-1
563 // Pn1z=(1/(2 ETv1^2))(epsilon Ev1 Sqrt[ M2noma1^2 + 4 M2noma1 PTmiss2cscdphi
564 // Pv1proj - 4 (ETv1^2 (Mn1^2 + qPTmiss2cscdphi^2) - PTmiss2cscdphi^2
565 // Pv1proj^2)] + M2noma1 Pv1z + 2 PTmiss2cscdphi Pv1proj Pv1z)
566 // with conditions: M2noma1 + 2 PTmiss2cscdphi Pv1proj + 2 Pn1z Pv1z > 0
567 // PTn1 -> PTmiss2 Csc[phi1 - phi2]
568
569 // if initialisation precompute some quantities
570 int solution_code = 0; // 0 with no solution, 1 with solution
571 nsol1 = 0;
572 nsol2 = 0;
573
574 // Variables used to test PTn1 and PTn2 > 0
575
576 const double &pTmissx = preparedInput.m_MEtX;
577 const double &pTmissy = preparedInput.m_MEtY;
578
579 fastSinCos(phi2, m_sinPhi2, m_cosPhi2);
580 double pTmiss2 = pTmissy * m_cosPhi2 - pTmissx * m_sinPhi2;
581
582 int dPhiSign = 0;
583 dPhiSign = fixPhiRange(phi1 - phi2) > 0 ? +1 : -1;
584
585 // Test if PTn1 and PTn2 > 0. Then MET vector is between the two neutrino
586 // vector
587
588 if (pTmiss2 * dPhiSign < 0) {
589 ++m_testptn1;
590 return solution_code;
591 }
592
593 fastSinCos(phi1, m_sinPhi1, m_cosPhi1);
594 double pTmiss1 = pTmissy * m_cosPhi1 - pTmissx * m_sinPhi1;
595
596 if (pTmiss1 * (-dPhiSign) < 0) {
597 ++m_testptn2;
598 return solution_code;
599 }
600
601 // Variables used to calculate discri1
602
603 double m2Vis1 = m_tauVec1M * m_tauVec1M;
604 m_ET2v1 = std::pow(m_tauVec1E, 2) - std::pow(m_tauVec1Pz, 2);
605 m_m2Nu1 = mNu1 * mNu1;
606 double m2noma1 = m_mTau2 - m_m2Nu1 - m2Vis1;
607 double m4noma1 = m2noma1 * m2noma1;
608 double pv1proj = m_tauVec1Px * m_cosPhi1 + m_tauVec1Py * m_sinPhi1;
609 double p2v1proj = std::pow(pv1proj, 2);
610 double sinDPhi2 = m_cosPhi2 * m_sinPhi1 - m_sinPhi2 * m_cosPhi1; // sin(Phi1-Phi2)
611 double pTmiss2CscDPhi = pTmiss2 / sinDPhi2;
612 double &pTn1 = pTmiss2CscDPhi;
613 double pT2miss2CscDPhi = pTmiss2CscDPhi * pTmiss2CscDPhi;
614
615 // Test on discri1
616 const double discri1 = m4noma1 + 4 * m2noma1 * pTmiss2CscDPhi * pv1proj -
617 4 * (m_ET2v1 * (m_m2Nu1 + pT2miss2CscDPhi) - (pT2miss2CscDPhi * p2v1proj));
618
619 if (discri1 < 0) // discriminant negative -> no solution
620 {
621 ++m_testdiscri1;
622 return solution_code;
623 }
624
625 // Variables used to calculate discri2
626 double m2Vis2 = m_tauVec2M * m_tauVec2M;
627 m_ET2v2 = std::pow(m_tauVec2E, 2) - std::pow(m_tauVec2Pz, 2);
628 m_m2Nu2 = mNu2 * mNu2;
629 double m2noma2 = m_mTau2 - m_m2Nu2 - m2Vis2;
630 double m4noma2 = m2noma2 * m2noma2;
631 double pv2proj = m_tauVec2Px * m_cosPhi2 + m_tauVec2Py * m_sinPhi2;
632 double p2v2proj = std::pow(pv2proj, 2);
633 double sinDPhi1 = -sinDPhi2;
634 double pTmiss1CscDPhi = pTmiss1 / sinDPhi1;
635 double &pTn2 = pTmiss1CscDPhi;
636 double pT2miss1CscDPhi = pTmiss1CscDPhi * pTmiss1CscDPhi;
637
638 const double discri2 = m4noma2 + 4 * m2noma2 * pTmiss1CscDPhi * pv2proj -
639 4 * (m_ET2v2 * (m_m2Nu2 + pT2miss1CscDPhi) - (pT2miss1CscDPhi * p2v2proj));
640
641 if (discri2 < 0) // discriminant negative -> no solution
642 {
643 ++m_testdiscri2;
644 return solution_code;
645 }
646
647 // this should be done only once we know there are solutions for nu2
648 m_E2v1 = m_tauVec1Pz * m_tauVec1Pz + m_ET2v1;
649 m_Ev1 = sqrt(m_E2v1);
650 double sqdiscri1 = sqrt(discri1);
651 double first1 =
652 (m2noma1 * m_tauVec1Pz + 2 * pTmiss2CscDPhi * pv1proj * m_tauVec1Pz) / (2 * m_ET2v1);
653 double second1 = sqdiscri1 * m_Ev1 / (2 * m_ET2v1);
654
655 // first solution
656 double pn1Z = first1 + second1;
657
658 if (m2noma1 + 2 * pTmiss2CscDPhi * pv1proj + 2 * pn1Z * m_tauVec1Pz >
659 0) // Condition for solution to exist
660 {
661 m_nuvecsol1[nsol1].SetPxPyPzE(pTn1 * m_cosPhi1, pTn1 * m_sinPhi1, pn1Z,
662 sqrt(std::pow(pTn1, 2) + std::pow(pn1Z, 2) + m_m2Nu1));
663
664 ++nsol1;
665 }
666
667 pn1Z = first1 - second1;
668
669 if (m2noma1 + 2 * pTmiss2CscDPhi * pv1proj + 2 * pn1Z * m_tauVec1Pz >
670 0) // Condition for solution to exist
671 {
672
673 m_nuvecsol1[nsol1].SetPxPyPzE(pTn1 * m_cosPhi1, pTn1 * m_sinPhi1, pn1Z,
674 sqrt(std::pow(pTn1, 2) + std::pow(pn1Z, 2) + m_m2Nu1));
675
676 ++nsol1;
677 }
678
679 if (nsol1 == 0) {
680 ++m_nosol1;
681 return solution_code;
682 }
683
684 m_E2v2 = m_tauVec2Pz * m_tauVec2Pz + m_ET2v2;
685 m_Ev2 = sqrt(m_E2v2);
686 double sqdiscri2 = sqrt(discri2);
687 double first2 =
688 (m2noma2 * m_tauVec2Pz + 2 * pTmiss1CscDPhi * pv2proj * m_tauVec2Pz) / (2 * m_ET2v2);
689 double second2 = sqdiscri2 * m_Ev2 / (2 * m_ET2v2);
690
691 // second solution
692 double pn2Z = first2 + second2;
693
694 if (m2noma2 + 2 * pTmiss1CscDPhi * pv2proj + 2 * pn2Z * m_tauVec2Pz >
695 0) // Condition for solution to exist
696 {
697 m_nuvecsol2[nsol2].SetPxPyPzE(pTn2 * m_cosPhi2, pTn2 * m_sinPhi2, pn2Z,
698 sqrt(std::pow(pTn2, 2) + std::pow(pn2Z, 2) + m_m2Nu2));
699
700 ++nsol2;
701 }
702
703 pn2Z = first2 - second2;
704 ;
705
706 if (m2noma2 + 2 * pTmiss1CscDPhi * pv2proj + 2 * pn2Z * m_tauVec2Pz >
707 0) // Condition for solution to exist
708 {
709 m_nuvecsol2[nsol2].SetPxPyPzE(pTn2 * m_cosPhi2, pTn2 * m_sinPhi2, pn2Z,
710 sqrt(std::pow(pTn2, 2) + std::pow(pn2Z, 2) + m_m2Nu2));
711
712 ++nsol2;
713 }
714
715 if (nsol2 == 0) {
716 ++m_nosol2;
717 return solution_code;
718 }
719
720 // Verification if solution exist
721
722 solution_code = 1;
723 ++m_iterNuPV3;
724
725 // double check solutions from time to time
726 if (m_iterNuPV3 % 1000 == 1) {
727 double pnux = m_nuvecsol1[0].Px() + m_nuvecsol2[0].Px();
728 double pnuy = m_nuvecsol1[0].Py() + m_nuvecsol2[0].Py();
729 double mtau1plus = (m_nuvecsol1[0] + m_tauVec1).M();
730 double mtau1moins = (m_nuvecsol1[1] + m_tauVec1).M();
731 double mtau2plus = (m_nuvecsol2[0] + m_tauVec2).M();
732 double mtau2moins = (m_nuvecsol2[1] + m_tauVec2).M();
733 if (std::abs(pnux - pTmissx) > 0.001 || std::abs(pnuy - pTmissy) > 0.001) {
734 Info("DiTauMassTools", "%s", ("NuPsolutionV3 ERROR Pnux-Met.X or Pnuy-Met.Y > 0.001 : " +
735 std::to_string(pnux - pTmissx) + " and " +
736 std::to_string(pnuy - pTmissx) + " " + "Invalid solutions")
737 .c_str());
738 }
739 if (std::abs(mtau1plus - m_mTau) > 0.001 || std::abs(mtau1moins - m_mTau) > 0.001 ||
740 std::abs(mtau2plus - m_mTau) > 0.001 || std::abs(mtau2moins - m_mTau) > 0.001) {
741 Info("DiTauMassTools", "%s", ("NuPsolutionV3 ERROR tau mass not recovered : " +
742 std::to_string(mtau1plus) + " " + std::to_string(mtau1moins) + " " +
743 std::to_string(mtau2plus) + " " + std::to_string(mtau2moins))
744 .c_str());
745 }
746 }
747
748 return solution_code;
749}
750
751// returns solution for Lepton Flavor Violating X->lep+tau studies
752int MissingMassCalculator::NuPsolutionLFV(const XYVector &met_vec,
753 const PtEtaPhiMVector &tau, const double &l_nu,
754 std::vector<PtEtaPhiMVector> &nu_vec) {
755 int solution_code = 0; // 0 with no solution, 1 with solution
756
757 nu_vec.clear();
758 PxPyPzMVector nu(met_vec.X(), met_vec.Y(), 0.0, l_nu);
759 PxPyPzMVector nu2(met_vec.X(), met_vec.Y(), 0.0, l_nu);
760
761 const double Mtau = ParticleConstants::tauMassInMeV / GEV;
762 // double msq = (Mtau*Mtau-tau.M()*tau.M())/2;
763 double msq = (Mtau * Mtau - tau.M() * tau.M() - l_nu * l_nu) /
764 2; // to take into account the fact that 2-nu systema has mass
765 double gamma = nu.Px() * nu.Px() + nu.Py() * nu.Py();
766 double beta = tau.Px() * nu.Px() + tau.Py() * nu.Py() + msq;
767 double a = tau.E() * tau.E() - tau.Pz() * tau.Pz();
768 double b = -2 * tau.Pz() * beta;
769 double c = tau.E() * tau.E() * gamma - beta * beta;
770 if ((b * b - 4 * a * c) < 0)
771 return solution_code; // no solution found
772 else
773 solution_code = 2;
774 double pvz1 = (-b + sqrt(b * b - 4 * a * c)) / (2 * a);
775 double pvz2 = (-b - sqrt(b * b - 4 * a * c)) / (2 * a);
776
777 nu.SetCoordinates(met_vec.X(), met_vec.Y(), pvz1, l_nu);
778 nu2.SetCoordinates(met_vec.X(), met_vec.Y(), pvz2, l_nu);
779
780 PtEtaPhiMVector return_nu(nu.Pt(), nu.Eta(), nu.Phi(), nu.M());
781 PtEtaPhiMVector return_nu2(nu2.Pt(), nu2.Eta(), nu2.Phi(), nu2.M());
782 nu_vec.push_back(return_nu);
783 nu_vec.push_back(return_nu2);
784 return solution_code;
785}
786
787// like v9fast, but the parameter space scanning is now factorised out, to allow
788// flexibility
789int MissingMassCalculator::DitauMassCalculatorV9walk() {
790
791 int nsuccesses = 0;
792
793 int fit_code = 0; // 0==bad, 1==good
794 ClearDitauStuff(m_fDitauStuffFit);
795 ClearDitauStuff(m_fDitauStuffHisto);
796 OutputInfo.m_AveSolRMS = 0.;
797
798 m_fMfit_all->Reset();
799
800 if(m_SaveLlhHisto){
801 m_fMEtP_all->Reset();
802 m_fMEtL_all->Reset();
803 m_fMnu1_all->Reset();
804 m_fMnu2_all->Reset();
805 m_fPhi1_all->Reset();
806 m_fPhi2_all->Reset();
807 }
808
809 m_fMfit_allNoWeight->Reset();
810 m_fPXfit1->Reset();
811 m_fPYfit1->Reset();
812 m_fPZfit1->Reset();
813 m_fPXfit2->Reset();
814 m_fPYfit2->Reset();
815 m_fPZfit2->Reset();
816
817 // these histograms are used for the floating stopping criterion
818 if (m_fUseFloatStopping) {
819 m_fMmass_split1->Reset();
820 m_fMEtP_split1->Reset();
821 m_fMEtL_split1->Reset();
822 m_fMnu1_split1->Reset();
823 m_fMnu2_split1->Reset();
824 m_fPhi1_split1->Reset();
825 m_fPhi2_split1->Reset();
826 m_fMmass_split2->Reset();
827 m_fMEtP_split2->Reset();
828 m_fMEtL_split2->Reset();
829 m_fMnu1_split2->Reset();
830 m_fMnu2_split2->Reset();
831 m_fPhi1_split2->Reset();
832 m_fPhi2_split2->Reset();
833 }
834
835 m_prob_tmp = 0.0;
836
837 m_iter1 = 0;
838
839 m_totalProbSum = 0;
840 m_mtautauSum = 0;
841
842 // initialize a spacewalker, which walks the parameter space according to some
843 // algorithm
844 SpaceWalkerInit();
845
846 while (SpaceWalkerWalk()) {
847 bool paramInsideRange = false;
848 m_nsol = 0;
849
850 paramInsideRange = checkAllParamInRange();
851
852 // FIXME if no tau scanning, or symmetric matrices, rotatin is made twice
853 // which is inefficient
854 const double deltaMetx = m_MEtL * m_metCovPhiCos - m_MEtP * m_metCovPhiSin;
855 const double deltaMety = m_MEtL * m_metCovPhiSin + m_MEtP * m_metCovPhiCos;
856
857 // deltaMetVec.Set(met_smear_x,met_smear_y);
858 preparedInput.m_metVec.SetXY(preparedInput.m_inputMEtX + deltaMetx,
859 preparedInput.m_inputMEtY + deltaMety);
860
861 // save in global variable for speed sake
862 preparedInput.m_MEtX = preparedInput.m_metVec.X();
863 preparedInput.m_MEtY = preparedInput.m_metVec.Y();
864 preparedInput.m_MEtT = preparedInput.m_metVec.R();
865
866 if (paramInsideRange)
867 probCalculatorV9fast(m_Phi1, m_Phi2, m_Mnu1, m_Mnu2);
868
869 // DR for markov chain need to enter handleSolution also when zero solutions
870 handleSolutions();
871 // be careful that with markov, current solution is from now on stored in
872 // XYZOldSolVec
873
874 if (m_nsol <= 0)
875 continue;
876
877 // for markov, nsuccess more difficult to define. Decide this is the number
878 // of independent point accepted (hence without weight)
879 nsuccesses = m_markovNAccept;
880 m_iterNsuc = m_markovNAccept;
881
882 m_iter1 += m_nsol;
883 fit_code = 1;
884
885 } // while loop
886
887 OutputInfo.m_NTrials = m_iter0;
888 OutputInfo.m_NSuccesses = nsuccesses;
889
890 if (nsuccesses > 0) {
891 OutputInfo.m_AveSolRMS /= nsuccesses;
892 } else {
893 OutputInfo.m_AveSolRMS = -1.;
894 }
895
896 double Px1, Py1, Pz1;
897 double Px2, Py2, Pz2;
898 if (nsuccesses > 0) {
899
900 // note that smoothing can slightly change the integral of the histogram
901
902#ifdef SMOOTH
903 m_fMfit_all->Smooth();
904 m_fMfit_allNoWeight->Smooth();
905 m_fPXfit1->Smooth();
906 m_fPYfit1->Smooth();
907 m_fPZfit1->Smooth();
908 m_fPXfit2->Smooth();
909 m_fPYfit2->Smooth();
910 m_fPZfit2->Smooth();
911#endif
912
913 // default max finding method defined in MissingMassCalculator.h
914 // note that window defined in terms of number of bin, so depend on binning
915 std::vector<double> histInfo(HistInfo::MAXHISTINFO);
916 m_fDitauStuffHisto.Mditau_best = maxFromHist(m_fMfit_all, histInfo);
917 double prob_hist = histInfo.at(HistInfo::PROB);
918
919 if (prob_hist != 0.0)
920 m_fDitauStuffHisto.Sign_best = -log10(std::abs(prob_hist));
921 else {
922 // this mean the histogram is empty.
923 // possible but very rare if all entries outside histogram range
924 // fall back to maximum
925 m_fDitauStuffHisto.Sign_best = -999.;
926 m_fDitauStuffHisto.Mditau_best = m_fDitauStuffFit.Mditau_best;
927 }
928
929 if (m_fDitauStuffHisto.Mditau_best > 0.0)
930 m_fDitauStuffHisto.RMSoverMPV = m_fMfit_all->GetRMS() / m_fDitauStuffHisto.Mditau_best;
931 std::vector<double> histInfoOther(HistInfo::MAXHISTINFO);
932 //---- getting full tau1 momentum
933 Px1 = maxFromHist(m_fPXfit1, histInfoOther);
934 Py1 = maxFromHist(m_fPYfit1, histInfoOther);
935 Pz1 = maxFromHist(m_fPZfit1, histInfoOther);
936
937 //---- getting full tau2 momentum
938 Px2 = maxFromHist(m_fPXfit2, histInfoOther);
939 Py2 = maxFromHist(m_fPYfit2, histInfoOther);
940 Pz2 = maxFromHist(m_fPZfit2, histInfoOther);
941
942 //---- setting 4-vecs
943 PxPyPzMVector fulltau1, fulltau2;
944 fulltau1.SetCoordinates(Px1, Py1, Pz1, ParticleConstants::tauMassInMeV / GEV);
945 fulltau2.SetCoordinates(Px2, Py2, Pz2, ParticleConstants::tauMassInMeV / GEV);
946 // PtEtaPhiMVector fulltau1(_fulltau1.Pt(), _fulltau1.Eta(), _fulltau1.Phi(), _fulltau1.M());
947 //PtEtaPhiMVector fulltau2(_fulltau2.Pt(), _fulltau2.Eta(), _fulltau2.Phi(), _fulltau2.M());
948
949 if (fulltau1.P() < preparedInput.m_vistau1.P())
950 fulltau1 = 1.01 * preparedInput.m_vistau1; // protection against cases when fitted tau
951 // momentum is smaller than visible tau momentum
952 if (fulltau2.P() < preparedInput.m_vistau2.P())
953 fulltau2 = 1.01 * preparedInput.m_vistau2; // protection against cases when fitted tau
954 // momentum is smaller than visible tau momentum
955 m_fDitauStuffHisto.vistau1 = preparedInput.m_vistau1; // FIXME should also be fitted if tau scan
956 m_fDitauStuffHisto.vistau2 = preparedInput.m_vistau2;
957 m_fDitauStuffHisto.nutau1 = fulltau1 - preparedInput.m_vistau1; // these are the original tau vis
958 m_fDitauStuffHisto.nutau2 =
959 fulltau2 - preparedInput.m_vistau2; // FIXME neutrino mass not necessarily zero
960 }
961
962 // Note that for v9walk, points outside the METx MEty disk are counted, while
963 // this was not the case for v9
964 if (preparedInput.m_fUseVerbose == 1) {
965 Info("DiTauMassTools", "Scanning ");
966 Info("DiTauMassTools", " Markov ");
967 Info("DiTauMassTools", "%s",
968 (" V9W niters=" + std::to_string(m_iter0) + " " + std::to_string(m_iter1)).c_str());
969 Info("DiTauMassTools", "%s", (" nFullScan " + std::to_string(m_markovNFullScan)).c_str());
970 Info("DiTauMassTools", "%s", (" nRejectNoSol " + std::to_string(m_markovNRejectNoSol)).c_str());
971 Info("DiTauMassTools", "%s", (" nRejectMetro " + std::to_string(m_markovNRejectMetropolis)).c_str());
972 Info("DiTauMassTools", "%s", (" nAccept " + std::to_string(m_markovNAccept)).c_str());
973 Info("DiTauMassTools", "%s",
974 (" probsum " + std::to_string(m_totalProbSum) + " msum " + std::to_string(m_mtautauSum))
975 .c_str());
976 }
977
978 if (preparedInput.m_fUseVerbose == 1) {
979 if (fit_code == 0) {
980 Info("DiTauMassTools", "%s", ("!!!----> Warning-3 in "
981 "MissingMassCalculator::DitauMassCalculatorV9Walk() : fit status=" +
982 std::to_string(fit_code))
983 .c_str());
984 Info("DiTauMassTools", "%s", "....... No solution is found. Printing input info .......");
985
986 Info("DiTauMassTools", "%s", (" vis Tau-1: Pt=" + std::to_string(preparedInput.m_vistau1.Pt()) +
987 " M=" + std::to_string(preparedInput.m_vistau1.M()) +
988 " eta=" + std::to_string(preparedInput.m_vistau1.Eta()) +
989 " phi=" + std::to_string(preparedInput.m_vistau1.Phi()) +
990 " type=" + std::to_string(preparedInput.m_type_visTau1))
991 .c_str());
992 Info("DiTauMassTools", "%s", (" vis Tau-2: Pt=" + std::to_string(preparedInput.m_vistau2.Pt()) +
993 " M=" + std::to_string(preparedInput.m_vistau2.M()) +
994 " eta=" + std::to_string(preparedInput.m_vistau2.Eta()) +
995 " phi=" + std::to_string(preparedInput.m_vistau2.Phi()) +
996 " type=" + std::to_string(preparedInput.m_type_visTau2))
997 .c_str());
998 Info("DiTauMassTools", "%s", (" MET=" + std::to_string(preparedInput.m_MetVec.R()) +
999 " Met_X=" + std::to_string(preparedInput.m_MetVec.X()) +
1000 " Met_Y=" + std::to_string(preparedInput.m_MetVec.Y()))
1001 .c_str());
1002 Info("DiTauMassTools", " ---------------------------------------------------------- ");
1003 }
1004 }
1005
1006 return fit_code;
1007}
1008
1009int MissingMassCalculator::DitauMassCalculatorV9lfv(bool refit) {
1010
1011 // debugThisIteration=false;
1012 m_debugThisIteration = true;
1013
1014 int fit_code = 0; // 0==bad, 1==good
1015 ClearDitauStuff(m_fDitauStuffFit);
1016 ClearDitauStuff(m_fDitauStuffHisto);
1017 OutputInfo.m_NTrials = 0;
1018 OutputInfo.m_NSuccesses = 0;
1019 OutputInfo.m_AveSolRMS = 0.;
1020
1021 //------- Settings -------------------------------
1022 int NiterMET = m_niter_fit2; // number of iterations for each MET scan loop
1023 int NiterMnu = m_niter_fit3; // number of iterations for Mnu loop
1024 const double Mtau = ParticleConstants::tauMassInMeV / GEV;
1025 double Mnu_binSize = m_MnuScanRange / NiterMnu;
1026
1027 double METresX = preparedInput.m_METsigmaL; // MET resolution in direction parallel to
1028 // leading jet, for MET scan
1029 double METresY = preparedInput.m_METsigmaP; // MET resolution in direction perpendicular to
1030 // leading jet, for MET scan
1031
1032 //-------- end of Settings
1033
1034 // if m_nsigma_METscan was not set by user, set to default values
1035 if(m_nsigma_METscan == -1){
1036 if (preparedInput.m_tauTypes == TauTypes::ll) { // both tau's are leptonic
1037 m_nsigma_METscan = m_nsigma_METscan_lfv_ll;
1038 } else if (preparedInput.m_tauTypes == TauTypes::lh) { // lep had
1039 m_nsigma_METscan = m_nsigma_METscan_lfv_lh;
1040 }
1041 }
1042
1043 double N_METsigma = m_nsigma_METscan; // number of sigmas for MET scan
1044 double METresX_binSize = 2 * N_METsigma * METresX / NiterMET;
1045 double METresY_binSize = 2 * N_METsigma * METresY / NiterMET;
1046
1047 int solution = 0;
1048
1049 std::vector<PtEtaPhiMVector> nu_vec;
1050
1051 m_totalProbSum = 0;
1052 m_mtautauSum = 0;
1053
1054 double metprob = 1.0;
1055 double sign_tmp = 0.0;
1056 double tauprob = 1.0;
1057 double totalProb = 0.0;
1058
1059 m_prob_tmp = 0.0;
1060
1061 double met_smear_x = 0.0;
1062 double met_smear_y = 0.0;
1063 double met_smearL = 0.0;
1064 double met_smearP = 0.0;
1065
1066 double angle1 = 0.0;
1067
1068 if (m_fMfit_all) {
1069 m_fMfit_all->Reset();
1070 }
1071 if (m_fMfit_allNoWeight) {
1072 m_fMfit_allNoWeight->Reset();
1073 }
1074 if (m_fPXfit1) {
1075 m_fPXfit1->Reset();
1076 }
1077 if (m_fPYfit1) {
1078 m_fPYfit1->Reset();
1079 }
1080 if (m_fPZfit1) {
1081 m_fPZfit1->Reset();
1082 }
1083
1084 int iter0 = 0;
1085 m_iter1 = 0;
1086 m_iter2 = 0;
1087 m_iter3 = 0;
1088 m_iter4 = 0;
1089
1090 const double met_coscovphi = cos(preparedInput.m_METcovphi);
1091 const double met_sincovphi = sin(preparedInput.m_METcovphi);
1092
1093 m_iang1low = 0;
1094 m_iang1high = 0;
1095
1096 // double Mvis=(tau_vec1+tau_vec2).M();
1097 // PtEtaPhiMVector met4vec(0.0,0.0,0.0,0.0);
1098 // met4vec.SetPxPyPzE(met_vec.X(),met_vec.Y(),0.0,met_vec.R());
1099 // double Meff=(tau_vec1+tau_vec2+met4vec).M();
1100 // double met_det=met_vec.R();
1101
1102 //---------------------------------------------
1103 if (preparedInput.m_tauTypes == TauTypes::ll) // dilepton case
1104 {
1105 if (preparedInput.m_fUseVerbose == 1) {
1106 Info("DiTauMassTools", "Running in dilepton mode");
1107 }
1108 double input_metX = preparedInput.m_MetVec.X();
1109 double input_metY = preparedInput.m_MetVec.Y();
1110
1111 PtEtaPhiMVector tau_tmp(0.0, 0.0, 0.0, 0.0);
1112 PtEtaPhiMVector lep_tmp(0.0, 0.0, 0.0, 0.0);
1113 int tau_type_tmp;
1114 int tau_ind = 0;
1115
1116 if (preparedInput.m_LFVmode == 1) // muon case: H->mu+tau(->ele) decays
1117 {
1118 if ((preparedInput.m_vistau1.M() > 0.05 &&
1119 preparedInput.m_vistau2.M() < 0.05) != refit) // choosing lepton from Higgs decay
1120 //When the mass calculator is rerun with refit==true the alternative lepton ordering is used
1121 {
1122 tau_tmp = preparedInput.m_vistau2;
1123 lep_tmp = preparedInput.m_vistau1;
1124 tau_type_tmp = preparedInput.m_type_visTau2;
1125 tau_ind = 2;
1126 } else {
1127 tau_tmp = preparedInput.m_vistau1;
1128 lep_tmp = preparedInput.m_vistau2;
1129 tau_type_tmp = preparedInput.m_type_visTau1;
1130 tau_ind = 1;
1131 }
1132 }
1133 if (preparedInput.m_LFVmode == 0) // electron case: H->ele+tau(->mu) decays
1134 {
1135 if ((preparedInput.m_vistau1.M() < 0.05 &&
1136 preparedInput.m_vistau2.M() > 0.05) != refit) // choosing lepton from Higgs decay
1137 //When the mass calculator is rerun with refit=true the alternative lepton ordering is used
1138 {
1139 tau_tmp = preparedInput.m_vistau2;
1140 lep_tmp = preparedInput.m_vistau1;
1141 tau_type_tmp = preparedInput.m_type_visTau2;
1142 tau_ind = 2;
1143 } else {
1144 tau_tmp = preparedInput.m_vistau1;
1145 lep_tmp = preparedInput.m_vistau2;
1146 tau_type_tmp = preparedInput.m_type_visTau1;
1147 tau_ind = 1;
1148 }
1149 }
1150
1151 //------- Settings -------------------------------
1152 double Mlep = tau_tmp.M();
1153 // double dMnu_max=m_MnuScanRange-Mlep;
1154 // double Mnu_binSize=dMnu_max/NiterMnu;
1155 //-------- end of Settings
1156
1157 // double M=Mtau;
1158 double M_nu = 0.0;
1159 double MnuProb = 1.0;
1160 //---------------------------------------------
1161 for (int i3 = 0; i3 < NiterMnu; i3++) //---- loop-3: virtual neutrino mass
1162 {
1163 M_nu = Mnu_binSize * i3;
1164 if (M_nu >= (Mtau - Mlep))
1165 continue;
1166 // M=sqrt(Mtau*Mtau-M_nu*M_nu);
1167 MnuProb = Prob->MnuProbability(preparedInput, M_nu,
1168 Mnu_binSize); // Mnu probability
1169 //---------------------------------------------
1170 for (int i4 = 0; i4 < NiterMET + 1; i4++) // MET_X scan
1171 {
1172 met_smearL = METresX_binSize * i4 - N_METsigma * METresX;
1173 for (int i5 = 0; i5 < NiterMET + 1; i5++) // MET_Y scan
1174 {
1175 met_smearP = METresY_binSize * i5 - N_METsigma * METresY;
1176 if (pow(met_smearL / METresX, 2) + pow(met_smearP / METresY, 2) > pow(N_METsigma, 2))
1177 continue; // use ellipse instead of square
1178 met_smear_x = met_smearL * met_coscovphi - met_smearP * met_sincovphi;
1179 met_smear_y = met_smearL * met_sincovphi + met_smearP * met_coscovphi;
1180 metvec_tmp.SetXY(input_metX + met_smear_x, input_metY + met_smear_y);
1181
1182 solution = NuPsolutionLFV(metvec_tmp, tau_tmp, M_nu, nu_vec);
1183
1184 ++iter0;
1185
1186 if (solution < 1)
1187 continue;
1188 ++m_iter1;
1189
1190 // if fast sin cos, result to not match exactly nupsolutionv2, so skip
1191 // test
1192 // SpeedUp no nested loop to compute individual probability
1193 int ngoodsol1 = 0;
1194
1195 metprob = Prob->MetProbability(preparedInput, met_smearL, met_smearP, METresX, METresY);
1196 if (metprob <= 0)
1197 continue;
1198 for (unsigned int j1 = 0; j1 < nu_vec.size(); j1++) {
1199 if (tau_tmp.E() + nu_vec[j1].E() >= preparedInput.m_beamEnergy)
1200 continue;
1201 const double tau1_tmpp = (tau_tmp + nu_vec[j1]).P();
1202 angle1 = Angle(nu_vec[j1], tau_tmp);
1203
1204 if (angle1 < dTheta3DLimit(tau_type_tmp, 0, tau1_tmpp)) {
1205 ++m_iang1low;
1206 continue;
1207 } // lower 99% bound
1208 if (angle1 > dTheta3DLimit(tau_type_tmp, 1, tau1_tmpp)) {
1209 ++m_iang1high;
1210 continue;
1211 } // upper 99% bound
1212 double tauvecprob1j =
1213 Prob->dTheta3d_probabilityFast(preparedInput, tau_type_tmp, angle1, tau1_tmpp);
1214 if (tauvecprob1j == 0.)
1215 continue;
1216 tauprob = Prob->TauProbabilityLFV(preparedInput, tau_type_tmp, tau_tmp, nu_vec[j1]);
1217 totalProb = tauvecprob1j * metprob * MnuProb * tauprob;
1218
1219 m_tautau_tmp.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);
1220 m_tautau_tmp += tau_tmp;
1221 m_tautau_tmp += lep_tmp;
1222 m_tautau_tmp += nu_vec[j1];
1223
1224 const double mtautau = m_tautau_tmp.M();
1225
1226 m_totalProbSum += totalProb;
1227 m_mtautauSum += mtautau;
1228
1229 fit_code = 1; // at least one solution is found
1230
1231 m_fMfit_all->Fill(mtautau, totalProb);
1232 m_fMfit_allNoWeight->Fill(mtautau, 1.);
1233 //----------------- using P*fit to fill Px,y,z_tau
1234 m_fPXfit1->Fill((tau_tmp + nu_vec[j1]).Px(), totalProb);
1235 m_fPYfit1->Fill((tau_tmp + nu_vec[j1]).Py(), totalProb);
1236 m_fPZfit1->Fill((tau_tmp + nu_vec[j1]).Pz(), totalProb);
1237
1238 if (totalProb > m_prob_tmp) // fill solution with highest probability
1239 {
1240 sign_tmp = -log10(totalProb);
1241 m_prob_tmp = totalProb;
1242 m_fDitauStuffFit.Mditau_best = mtautau;
1243 m_fDitauStuffFit.Sign_best = sign_tmp;
1244 if (tau_ind == 1)
1245 m_fDitauStuffFit.nutau1 = nu_vec[j1];
1246 if (tau_ind == 2)
1247 m_fDitauStuffFit.nutau2 = nu_vec[j1];
1248 }
1249
1250 ++ngoodsol1;
1251 }
1252
1253 if (ngoodsol1 == 0)
1254 continue;
1255 m_iter2 += 1;
1256
1257 m_iter3 += 1;
1258 }
1259 }
1260 }
1261 } else if (preparedInput.m_tauTypes == TauTypes::lh) // lepton+tau case
1262 {
1263 if (preparedInput.m_fUseVerbose == 1) {
1264 Info("DiTauMassTools", "Running in lepton+tau mode");
1265 }
1266 //------- Settings -------------------------------
1267
1268 //----- Stuff below are for Winter 2012 lep-had analysis only; it has to be
1269 // replaced by a more common scheme once other channels are optimized
1270 // XYVector
1271 // mht_vec((tau_vec1+tau_vec2).Px(),(tau_vec1+tau_vec2).Py()); //
1272 // missing Ht vector for Njet25=0 events const double
1273 // mht=mht_vec.R();
1274 double input_metX = preparedInput.m_MetVec.X();
1275 double input_metY = preparedInput.m_MetVec.Y();
1276
1277 // double mht_offset=0.0;
1278 // if(InputInfo.UseHT) // use missing Ht (for 0-jet events only for
1279 // now)
1280 // {
1281 // input_metX=-mht_vec.X();
1282 // input_metY=-mht_vec.Y();
1283 // }
1284 // else // use MET (for 0-jet and 1-jet events)
1285 // {
1286 // input_metX=met_vec.X();
1287 // input_metY=met_vec.Y();
1288 // }
1289
1290 PtEtaPhiMVector tau_tmp(0.0, 0.0, 0.0, 0.0);
1291 PtEtaPhiMVector lep_tmp(0.0, 0.0, 0.0, 0.0);
1292 int tau_type_tmp;
1293 if (preparedInput.m_type_visTau1 == 8) {
1294 tau_tmp = preparedInput.m_vistau2;
1295 lep_tmp = preparedInput.m_vistau1;
1296 tau_type_tmp = preparedInput.m_type_visTau2;
1297 }
1298 if (preparedInput.m_type_visTau2 == 8) {
1299 tau_tmp = preparedInput.m_vistau1;
1300 lep_tmp = preparedInput.m_vistau2;
1301 tau_type_tmp = preparedInput.m_type_visTau1;
1302 }
1303
1304 //---------------------------------------------
1305 for (int i4 = 0; i4 < NiterMET + 1; i4++) // MET_X scan
1306 {
1307 met_smearL = METresX_binSize * i4 - N_METsigma * METresX;
1308 for (int i5 = 0; i5 < NiterMET + 1; i5++) // MET_Y scan
1309 {
1310 met_smearP = METresY_binSize * i5 - N_METsigma * METresY;
1311 if (pow(met_smearL / METresX, 2) + pow(met_smearP / METresY, 2) > pow(N_METsigma, 2))
1312 continue; // use ellipse instead of square
1313 met_smear_x = met_smearL * m_metCovPhiCos - met_smearP * m_metCovPhiSin;
1314 met_smear_y = met_smearL * m_metCovPhiSin + met_smearP * m_metCovPhiCos;
1315 metvec_tmp.SetXY(input_metX + met_smear_x, input_metY + met_smear_y);
1316
1317 solution = NuPsolutionLFV(metvec_tmp, tau_tmp, 0.0, nu_vec);
1318
1319 ++iter0;
1320
1321 if (solution < 1)
1322 continue;
1323 ++m_iter1;
1324
1325 // if fast sin cos, result to not match exactly nupsolutionv2, so skip
1326 // test
1327 // SpeedUp no nested loop to compute individual probability
1328 int ngoodsol1 = 0;
1329
1330 metprob = Prob->MetProbability(preparedInput, met_smearL, met_smearP, METresX, METresY);
1331 if (metprob <= 0)
1332 continue;
1333 for (unsigned int j1 = 0; j1 < nu_vec.size(); j1++) {
1334 if (tau_tmp.E() + nu_vec[j1].E() >= preparedInput.m_beamEnergy)
1335 continue;
1336 const double tau1_tmpp = (tau_tmp + nu_vec[j1]).P();
1337 angle1 = Angle(nu_vec[j1], tau_tmp);
1338
1339 if (angle1 < dTheta3DLimit(tau_type_tmp, 0, tau1_tmpp)) {
1340 ++m_iang1low;
1341 continue;
1342 } // lower 99% bound
1343 if (angle1 > dTheta3DLimit(tau_type_tmp, 1, tau1_tmpp)) {
1344 ++m_iang1high;
1345 continue;
1346 } // upper 99% bound
1347 double tauvecprob1j =
1348 Prob->dTheta3d_probabilityFast(preparedInput, tau_type_tmp, angle1, tau1_tmpp);
1349 if (tauvecprob1j == 0.)
1350 continue;
1351 tauprob = Prob->TauProbabilityLFV(preparedInput, tau_type_tmp, tau_tmp, nu_vec[j1]);
1352 totalProb = tauvecprob1j * metprob * tauprob;
1353
1354 m_tautau_tmp.SetPxPyPzE(0.0, 0.0, 0.0, 0.0);
1355 m_tautau_tmp += tau_tmp;
1356 m_tautau_tmp += lep_tmp;
1357 m_tautau_tmp += nu_vec[j1];
1358
1359 const double mtautau = m_tautau_tmp.M();
1360
1361 m_totalProbSum += totalProb;
1362 m_mtautauSum += mtautau;
1363
1364 fit_code = 1; // at least one solution is found
1365
1366 m_fMfit_all->Fill(mtautau, totalProb);
1367 m_fMfit_allNoWeight->Fill(mtautau, 1.);
1369 // m_fPXfit1->Fill((tau_tmp+nu_vec[j1]).Px(),totalProb);
1370 // m_fPYfit1->Fill((tau_tmp+nu_vec[j1]).Py(),totalProb);
1371 // m_fPZfit1->Fill((tau_tmp+nu_vec[j1]).Pz(),totalProb);
1372
1373 if (totalProb > m_prob_tmp) // fill solution with highest probability
1374 {
1375 sign_tmp = -log10(totalProb);
1376 m_prob_tmp = totalProb;
1377 m_fDitauStuffFit.Mditau_best = mtautau;
1378 m_fDitauStuffFit.Sign_best = sign_tmp;
1379 if (preparedInput.m_type_visTau1 == 8) {
1380 m_fDitauStuffFit.vistau1 = lep_tmp;
1381 m_fDitauStuffFit.vistau2 = tau_tmp;
1382 m_fDitauStuffFit.nutau2 = nu_vec[j1];
1383 } else if (preparedInput.m_type_visTau2 == 8) {
1384 m_fDitauStuffFit.vistau2 = lep_tmp;
1385 m_fDitauStuffFit.vistau1 = tau_tmp;
1386 m_fDitauStuffFit.nutau1 = nu_vec[j1];
1387 }
1388 }
1389
1390 ++ngoodsol1;
1391 }
1392
1393 if (ngoodsol1 == 0)
1394 continue;
1395 m_iter2 += 1;
1396
1397 m_iter3 += 1;
1398 }
1399 }
1400 } else {
1401 Info("DiTauMassTools", "Running in an unknown mode?!?!");
1402 }
1403
1404 OutputInfo.m_NTrials = iter0;
1405 OutputInfo.m_NSuccesses = m_iter3;
1406
1407 if (preparedInput.m_fUseVerbose == 1) {
1408 Info("DiTauMassTools", "%s",
1409 ("SpeedUp niters=" + std::to_string(iter0) + " " + std::to_string(m_iter1) + " " +
1410 std::to_string(m_iter2) + " " + std::to_string(m_iter3) + "skip:" + std::to_string(m_iang1low) +
1411 " " + std::to_string(m_iang1high))
1412 .c_str());
1413 }
1414
1415 if (m_fMfit_all && m_fMfit_all->GetEntries() > 0 && m_iter3 > 0) {
1416#ifdef SMOOTH
1417 m_fMfit_all->Smooth();
1418 m_fMfit_allNoWeight->Smooth();
1419 m_fPXfit1->Smooth();
1420 m_fPYfit1->Smooth();
1421 m_fPZfit1->Smooth();
1422#endif
1423
1424 // default max finding method defined in MissingMassCalculator.h
1425 // note that window defined in terms of number of bin, so depend on binning
1426 std::vector<double> histInfo(HistInfo::MAXHISTINFO);
1427 m_fDitauStuffHisto.Mditau_best = maxFromHist(m_fMfit_all, histInfo);
1428 double prob_hist = histInfo.at(HistInfo::PROB);
1429
1430 if (prob_hist != 0.0)
1431 m_fDitauStuffHisto.Sign_best = -log10(std::abs(prob_hist));
1432 else {
1433 // this mean the histogram is empty.
1434 // possible but very rare if all entries outside histogram range
1435 // fall back to maximum
1436 m_fDitauStuffHisto.Sign_best = -999.;
1437 m_fDitauStuffHisto.Mditau_best = m_fDitauStuffFit.Mditau_best;
1438 }
1439
1440 if (m_fDitauStuffHisto.Mditau_best > 0.0)
1441 m_fDitauStuffHisto.RMSoverMPV = m_fMfit_all->GetRMS() / m_fDitauStuffHisto.Mditau_best;
1442 std::vector<double> histInfoOther(HistInfo::MAXHISTINFO);
1443 //---- getting Nu1
1444 double Px1 = maxFromHist(m_fPXfit1, histInfoOther);
1445 double Py1 = maxFromHist(m_fPYfit1, histInfoOther);
1446 double Pz1 = maxFromHist(m_fPZfit1, histInfoOther);
1447 //---- setting 4-vecs
1448 PxPyPzMVector nu1_tmp(0.0, 0.0, 0.0, 0.0);
1449 PxPyPzMVector nu2_tmp(0.0, 0.0, 0.0, 0.0);
1450 if (preparedInput.m_type_visTau1 == 8) {
1451 nu1_tmp = preparedInput.m_vistau1;
1452 nu2_tmp.SetCoordinates(Px1, Py1, Pz1, ParticleConstants::tauMassInMeV / GEV);
1453 }
1454 if (preparedInput.m_type_visTau2 == 8) {
1455 nu2_tmp = preparedInput.m_vistau2;
1456 nu1_tmp.SetCoordinates(Px1, Py1, Pz1, ParticleConstants::tauMassInMeV / GEV);
1457 }
1458 m_fDitauStuffHisto.nutau1 = nu1_tmp - preparedInput.m_vistau1;
1459 m_fDitauStuffHisto.nutau2 = nu2_tmp - preparedInput.m_vistau2;
1460 }
1461 if (m_lfvLeplepRefit && fit_code==0 && !refit) {
1462 fit_code = DitauMassCalculatorV9lfv(true);
1463 return fit_code;
1464 }
1465
1466
1467
1468 if (preparedInput.m_fUseVerbose == 1) {
1469 if (fit_code == 0) {
1470 Info(
1471 "DiTauMassTools", "%s",
1472 ("!!!----> Warning-3 in MissingMassCalculator::DitauMassCalculatorV9lfv() : fit status=" +
1473 std::to_string(fit_code))
1474 .c_str());
1475 Info("DiTauMassTools", "....... No solution is found. Printing input info .......");
1476
1477 Info("DiTauMassTools", "%s", (" vis Tau-1: Pt="+std::to_string(preparedInput.m_vistau1.Pt())
1478 +" M="+std::to_string(preparedInput.m_vistau1.M())+" eta="+std::to_string(preparedInput.m_vistau1.Eta())
1479 +" phi="+std::to_string(preparedInput.m_vistau1.Phi())
1480 +" type="+std::to_string(preparedInput.m_type_visTau1)).c_str());
1481 Info("DiTauMassTools", "%s", (" vis Tau-2: Pt="+std::to_string(preparedInput.m_vistau2.Pt())
1482 +" M="+std::to_string(preparedInput.m_vistau2.M())+" eta="+std::to_string(preparedInput.m_vistau2.Eta())
1483 +" phi="+std::to_string(preparedInput.m_vistau2.Phi())
1484 +" type="+std::to_string(preparedInput.m_type_visTau2)).c_str());
1485 Info("DiTauMassTools", "%s", (" MET="+std::to_string(preparedInput.m_MetVec.R())+" Met_X="+std::to_string(preparedInput.m_MetVec.X())
1486 +" Met_Y="+std::to_string(preparedInput.m_MetVec.Y())).c_str());
1487 Info("DiTauMassTools", " ---------------------------------------------------------- ");
1488 }
1489 }
1490 return fit_code;
1491}
1492
1493// function to fit maximum
1494Double_t MissingMassCalculator::maxFitting(Double_t *x, Double_t *par)
1495// Double_t maxFitting(Double_t *x, Double_t *par)
1496{
1497 // parabola with parameters max, mean and invwidth
1498 const double mM = x[0];
1499 const double mMax = par[0];
1500 const double mMean = par[1];
1501 const double mInvWidth2 = par[2]; // if param positif distance between intersection of the
1502 // parabola with x axis: 1/Sqrt(mInvWidth2)
1503 const double fitval = mMax * (1 - 4 * mInvWidth2 * std::pow(mM - mMean, 2));
1504 return fitval;
1505}
1506
1507// determine the maximum from the histogram
1508// if input prob not default , compute also some probability
1509// MaxHistStrategy : different method to find maximum
1510// TODO should get the array on work on it
1511// should also find the effective range of the hist
1512
1513double
1514MissingMassCalculator::maxFromHist(TH1F *theHist, std::vector<double> &histInfo,
1515 const MaxHistStrategy::e maxHistStrategy,
1516 const int winHalfWidth, bool debug) {
1517 // namespace HistInfo
1518 // enum e {
1519 // PROB=0,INTEGRAL,CHI2,DISCRI,TANTHETA,TANTHETAW,FITLENGTH,RMS,RMSVSDISCRI,MAXHISTINFO
1520 // };
1521 if (!theHist)[[unlikely]]{
1522 throw std::runtime_error("MissingMassCalculator::maxFromHist: histogram pointer is null.");
1523 }
1524 double maxPos = 0.;
1525 double prob = 0.;
1526
1527 for (std::vector<double>::iterator itr = histInfo.begin(); itr != histInfo.end(); ++itr) {
1528 *itr = -1;
1529 }
1530
1531 histInfo[HistInfo::INTEGRAL] = theHist->Integral();
1532
1533 if (maxHistStrategy == MaxHistStrategy::MAXBIN ||
1534 ((maxHistStrategy == MaxHistStrategy::MAXBINWINDOW ||
1535 maxHistStrategy == MaxHistStrategy::SLIDINGWINDOW) &&
1536 winHalfWidth == 0)) {
1537
1538 // simple max search
1539 // original version, simple bin maximum
1540 int max_bin = theHist->GetMaximumBin();
1541 maxPos = theHist->GetBinCenter(max_bin);
1542
1543 // FIXME GetEntries is unweighted
1544 prob = theHist->GetBinContent(max_bin) / double(theHist->GetEntries());
1545 if (prob > 1.)
1546 prob = 1.;
1547 histInfo[HistInfo::PROB] = prob;
1548 return maxPos;
1549 }
1550
1551 int hNbins = theHist->GetNbinsX();
1552
1553 if (maxHistStrategy == MaxHistStrategy::MAXBINWINDOW) {
1554 // average around maximum bin (nearly useless in fact)
1555 // could be faster
1556 int max_bin = theHist->GetMaximumBin();
1557 int iBinMin = max_bin - winHalfWidth;
1558 if (iBinMin < 0)
1559 iBinMin = 0;
1560 int iBinMax = max_bin + winHalfWidth;
1561 if (iBinMax > hNbins)
1562 iBinMax = hNbins - 1;
1563 double sumw = 0;
1564 double sumx = 0;
1565 for (int iBin = iBinMin; iBin <= iBinMax; ++iBin) {
1566 const double weight = theHist->GetBinContent(iBin);
1567 sumw += weight;
1568 sumx += weight * theHist->GetBinCenter(iBin);
1569 }
1570 maxPos = (sumw != 0.) ? (sumx / sumw) : 0.;
1571
1572 // FIXME GetEntries is unweighted
1573 prob = sumw / theHist->GetEntries();
1574 if (prob > 1.)
1575 prob = 1.;
1576
1577 return maxPos;
1578 }
1579
1580 // now compute sliding window anyway
1581 if (maxHistStrategy != MaxHistStrategy::SLIDINGWINDOW &&
1582 maxHistStrategy != MaxHistStrategy::FIT) {
1583 Error("DiTauMassTools", "%s",
1584 ("ERROR undefined maxHistStrategy:" + std::to_string(maxHistStrategy)).c_str());
1585 return -10.;
1586 }
1587
1588 // first iteration to find the first and last non zero bin, and the histogram
1589 // integral (not same as Entries because of weights)
1590 int lastNonZeroBin = -1;
1591 int firstNonZeroBin = -1;
1592 double totalSumw = 0.;
1593 bool firstNullPart = true;
1594 for (int iBin = 0; iBin < hNbins; ++iBin) {
1595 const double weight = theHist->GetBinContent(iBin);
1596 if (weight > 0) {
1597 totalSumw += weight;
1598 lastNonZeroBin = iBin;
1599 if (firstNullPart) {
1600 firstNullPart = false;
1601 firstNonZeroBin = iBin;
1602 }
1603 }
1604 }
1605
1606 // enlarge first and last non zero bin with window width to avoid side effect
1607 // (maximum close to the edge)
1608 firstNonZeroBin = std::max(0, firstNonZeroBin - winHalfWidth - 1);
1609 lastNonZeroBin = std::min(hNbins - 1, lastNonZeroBin + winHalfWidth + 1);
1610
1611 // if null histogram quit
1612 if (firstNullPart)
1613 return maxPos;
1614
1615 // determine the size of the sliding window in the fit case
1616
1617 // sliding window
1618 const int nwidth = 2 * winHalfWidth + 1;
1619 double winsum = 0.;
1620
1621 for (int ibin = 0; ibin < nwidth; ++ibin) {
1622 winsum += theHist->GetBinContent(ibin);
1623 }
1624 double winmax = winsum;
1625
1626 int max_bin = 0.;
1627 int iBinL = firstNonZeroBin;
1628 int iBinR = iBinL + 2 * winHalfWidth;
1629 bool goingUp = true;
1630
1631 do {
1632 ++iBinL;
1633 ++iBinR;
1634 const double deltawin = theHist->GetBinContent(iBinR) - theHist->GetBinContent(iBinL - 1);
1635
1636 if (deltawin < 0) {
1637 if (goingUp) {
1638 // if were climbing and now loose more on the left
1639 // than win on the right. This was a local maxima
1640 if (winsum > winmax) {
1641 // global maximum one so far
1642 winmax = winsum;
1643 max_bin = (iBinR + iBinL) / 2 - 1;
1644 }
1645 goingUp = false; // now going down
1646 }
1647 } else {
1648 // do not care about minima, simply indicate we are going down
1649 goingUp = true;
1650 }
1651
1652 winsum += deltawin;
1653
1654 } while (iBinR < lastNonZeroBin);
1655
1656 // now compute average
1657 int iBinMin = max_bin - winHalfWidth;
1658 if (iBinMin < 0)
1659 iBinMin = 0;
1660 int iBinMax = max_bin + winHalfWidth;
1661 if (iBinMax >= hNbins)
1662 iBinMax = hNbins - 1;
1663 double sumw = 0;
1664 double sumx = 0;
1665 for (int iBin = iBinMin; iBin <= iBinMax; ++iBin) {
1666 const double weight = theHist->GetBinContent(iBin);
1667 sumw += weight;
1668 sumx += weight * theHist->GetBinCenter(iBin);
1669 }
1670
1671 double maxPosWin = -1.;
1672
1673 if (sumw > 0.) {
1674 maxPosWin = sumx / sumw;
1675 }
1676 // prob if the fraction of events in the window
1677 prob = sumw / totalSumw;
1678
1679 // Definitions of some useful parameters
1680
1681 const double h_rms = theHist->GetRMS(1);
1682 histInfo[HistInfo::RMS] = h_rms;
1683
1684 double num = 0;
1685 double numerator = 0;
1686 double denominator = 0;
1687 bool nullBin = false;
1688
1689 for (int i = iBinMin; i < iBinMax; ++i) {
1690 double binError = theHist->GetBinError(i);
1691 if (binError < 1e-10) {
1692 nullBin = true;
1693 }
1694 double binErrorSquare = std::pow(binError, 2);
1695 num = theHist->GetBinContent(i) / (binErrorSquare);
1696 numerator = numerator + num;
1697 denominator = denominator + (1 / (binErrorSquare));
1698 }
1699 if (numerator < 1e-10 || denominator < 1e-10 || nullBin == true) {
1700 histInfo[HistInfo::MEANBIN] = -1;
1701 } else {
1702 histInfo[HistInfo::MEANBIN] = sqrt(1 / denominator) / (numerator / denominator);
1703 }
1704
1705 // stop here if only looking for sliding window
1706 if (maxHistStrategy == MaxHistStrategy::SLIDINGWINDOW) {
1707 return maxPosWin;
1708 }
1709
1710 maxPos = maxPosWin;
1711 // now FIT maxHistStrategy==MaxHistStrategy::FIT
1712
1713 // now mass fit in range defined by sliding window
1714 // window will be around maxPos
1715 const double binWidth = theHist->GetBinCenter(2) - theHist->GetBinCenter(1);
1716 double fitWidth = (winHalfWidth + 0.5) * binWidth;
1717 // fit range 2 larger than original window range, 3 if less than 20% of the
1718 // histogram in slinding window
1719
1720 if (prob > 0.2) {
1721 fitWidth *= 2.;
1722 } else {
1723 fitWidth *= 3.;
1724 }
1725 // fit option : Q == Quiet, no printout S result of the fit returned in
1726 // TFitResultPtr N do not draw the resulting function
1727
1728 // if debug plot the fitted function
1729 TString fitOption = debug ? "QS" : "QNS";
1730 // root fit
1731 // Sets initial values
1732 m_fFitting->SetParameters(sumw / winHalfWidth, maxPos, 0.0025);
1733 // TFitResultPtr
1734 // fitRes=theHist->Fit("pol2",fitOption,"",maxPos-fitWidth,maxPos+fitWidth);
1735 TFitResultPtr fitRes =
1736 theHist->Fit(m_fFitting, fitOption, "", maxPos - fitWidth, maxPos + fitWidth);
1737
1738 double maxPosFit = -1.;
1739
1740 if (int(fitRes) == 0) {
1741 // root fit
1742 histInfo[HistInfo::CHI2] = fitRes->Chi2();
1743 const double mMax = fitRes->Parameter(0);
1744 const double mMean = fitRes->Parameter(1);
1745 const double mInvWidth2 = fitRes->Parameter(2);
1746 double mMaxError = fitRes->ParError(0);
1747 m_PrintmMaxError = mMaxError;
1748 double mMeanError = fitRes->ParError(1);
1749 m_PrintmMeanError = mMeanError;
1750 double mInvWidth2Error = fitRes->ParError(2);
1751 m_PrintmInvWidth2Error = mInvWidth2Error;
1752 mMeanError = 0.; // avoid warning
1753 mInvWidth2Error = 0.; // avoid warning
1754 const double c = mMax * (1 - 4 * mMean * mMean * mInvWidth2);
1755 const double b = 8 * mMax * mMean * mInvWidth2;
1756 const double a = -4 * mMax * mInvWidth2;
1757 // when built in polynomial fit
1758 // const double c=fitRes->Parameter(0);
1759 // const double b=fitRes->Parameter(1);
1760 // const double a=fitRes->Parameter(2);
1761
1762 const double h_discri = b * b - 4 * a * c;
1763 histInfo[HistInfo::DISCRI] = h_discri;
1764 const double sqrth_discri = sqrt(h_discri);
1765 const double h_fitLength = sqrth_discri / a;
1766 histInfo[HistInfo::FITLENGTH] = h_fitLength;
1767 histInfo[HistInfo::TANTHETA] = 2 * a / sqrth_discri;
1768 histInfo[HistInfo::TANTHETAW] = 2 * a * sumw / sqrth_discri;
1769 histInfo[HistInfo::RMSVSDISCRI] = h_rms / h_fitLength;
1770 // compute maximum position (only if inverted parabola)
1771 if (a < 0)
1772 maxPosFit = -b / (2 * a);
1773 }
1774
1775 // keep fit result only if within 80% of fit window, and fit succeeded
1776 if (maxPosFit >= 0. and std::abs(maxPosFit - maxPosWin) < 0.8 * fitWidth) {
1777 histInfo[HistInfo::PROB] = prob;
1778 return maxPosFit;
1779 } else {
1780 // otherwise keep the weighted average
1781 // negate prob just to flag such event
1782 prob = -prob;
1783 histInfo[HistInfo::PROB] = prob;
1784 return maxPosWin;
1785 }
1786}
1787
1788// compute probability for any input value,can be called from a pure parameter
1789// scan
1790// deltametvec is along phijet
1791// returns number of solution if positive, return code if negative, vector of
1792// probability and mass
1793int MissingMassCalculator::probCalculatorV9fast(const double &phi1, const double &phi2,
1794 const double &M_nu1,
1795 const double &M_nu2) {
1796 // bool debug=true;
1797
1798 int nsol1;
1799 int nsol2;
1800
1801 const int solution = NuPsolutionV3(M_nu1, M_nu2, phi1, phi2, nsol1, nsol2);
1802
1803 if (solution != 1)
1804 return -4;
1805 // refineSolutions ( M_nu1,M_nu2,
1806 // met_smearL,met_smearP,metvec_tmp.R(),
1807 // nsol1, nsol2,m_Mvis,m_Meff);
1808 refineSolutions(M_nu1, M_nu2, nsol1, nsol2, m_Mvis, m_Meff);
1809
1810 if (m_nsol <= 0)
1811 return 0;
1812
1813 // success
1814
1815 return m_nsol; // for backward compatibility
1816}
1817
1818// nuvecsol1 and nuvecsol2 passed by MMC
1819int MissingMassCalculator::refineSolutions(const double &M_nu1, const double &M_nu2,
1820 const int nsol1, const int nsol2,
1821 const double &Mvis, const double &Meff)
1822
1823{
1824 m_nsol = 0;
1825
1826 if (int(m_probFinalSolVec.size()) < m_nsolfinalmax)
1827 Error("DiTauMassTools", "%s",
1828 ("refineSolutions ERROR probFinalSolVec.size() should be " + std::to_string(m_nsolfinalmax))
1829 .c_str());
1830 if (int(m_mtautauFinalSolVec.size()) < m_nsolfinalmax)
1831 Error("DiTauMassTools", "%s",
1832 ("refineSolutions ERROR mtautauSolVec.size() should be " + std::to_string(m_nsolfinalmax))
1833 .c_str());
1834 if (int(m_nu1FinalSolVec.size()) < m_nsolfinalmax)
1835 Error("DiTauMassTools", "%s",
1836 ("refineSolutions ERROR nu1FinalSolVec.size() should be " + std::to_string(m_nsolfinalmax))
1837 .c_str());
1838 if (int(m_nu2FinalSolVec.size()) < m_nsolfinalmax)
1839 Error("DiTauMassTools", "%s",
1840 ("refineSolutions ERROR nu2FinalSolVec.size() should be " + std::to_string(m_nsolfinalmax))
1841 .c_str());
1842 if (nsol1 > int(m_nsolmax))
1843 Error("DiTauMassTools", "%s", ("refineSolutions ERROR nsol1 " + std::to_string(nsol1) +
1844 " > nsolmax !" + std::to_string(m_nsolmax))
1845 .c_str());
1846 if (nsol2 > int(m_nsolmax))
1847 Error("DiTauMassTools", "%s", ("refineSolutions ERROR nsol1 " + std::to_string(nsol2) +
1848 " > nsolmax !" + std::to_string(m_nsolmax))
1849 .c_str());
1850
1851 int ngoodsol1 = 0;
1852 int ngoodsol2 = 0;
1853 double constProb =
1854 Prob->apply(preparedInput, -99, -99, PtEtaPhiMVector(0, 0, 0, 0), PtEtaPhiMVector(0, 0, 0, 0),
1855 PtEtaPhiMVector(0, 0, 0, 0), PtEtaPhiMVector(0, 0, 0, 0), true, false, false);
1856
1857 for (int j1 = 0; j1 < nsol1; ++j1) {
1858 PtEtaPhiMVector &nuvec1_tmpj = m_nuvecsol1[j1];
1859 PtEtaPhiMVector &tauvecsol1j = m_tauvecsol1[j1];
1860 double &tauvecprob1j = m_tauvecprob1[j1];
1861 tauvecprob1j = 0.;
1862 // take first or second solution
1863 // no time to call rndm, switch more or less randomely, according to an
1864 // oscillating switch perturbed by m_phi1
1865 if (nsol1 > 1) {
1866 if (j1 == 0) { // decide at the first solution which one we will take
1867 const int pickInt = std::abs(10000 * m_Phi1);
1868 const int pickDigit = pickInt - 10 * (pickInt / 10);
1869 if (pickDigit < 5)
1870 m_switch1 = !m_switch1;
1871 }
1872 m_switch1 = !m_switch1;
1873 }
1874
1875 if (!m_switch1) {
1876 nuvec1_tmpj.SetCoordinates(nuvec1_tmpj.Pt(), nuvec1_tmpj.Eta(), nuvec1_tmpj.Phi(), M_nu1);
1877 tauvecsol1j.SetPxPyPzE(0., 0., 0., 0.);
1878 tauvecsol1j += nuvec1_tmpj;
1879 tauvecsol1j += m_tauVec1;
1880 if (tauvecsol1j.E() >= preparedInput.m_beamEnergy)
1881 continue;
1882 tauvecprob1j = Prob->apply(preparedInput, preparedInput.m_type_visTau1, -99, m_tauVec1,
1883 PtEtaPhiMVector(0, 0, 0, 0), nuvec1_tmpj,
1884 PtEtaPhiMVector(0, 0, 0, 0), false, true, false);
1885 ++ngoodsol1;
1886 }
1887
1888 for (int j2 = 0; j2 < nsol2; ++j2) {
1889 PtEtaPhiMVector &nuvec2_tmpj = m_nuvecsol2[j2];
1890 PtEtaPhiMVector &tauvecsol2j = m_tauvecsol2[j2];
1891 double &tauvecprob2j = m_tauvecprob2[j2];
1892 if (j1 == 0) {
1893 tauvecprob2j = 0.;
1894 // take first or second solution
1895 // no time to call rndm, switch more or less randomely, according to an
1896 // oscillating switch perturbed by m_phi2
1897 if (nsol2 > 1) {
1898 if (j2 == 0) { // decide at the first solution which one we will take
1899 const int pickInt = std::abs(10000 * m_Phi2);
1900 const int pickDigit = pickInt - 10 * int(pickInt / 10);
1901 if (pickDigit < 5)
1902 m_switch2 = !m_switch2;
1903 }
1904 m_switch2 = !m_switch2;
1905 }
1906
1907 if (!m_switch2) {
1908 nuvec2_tmpj.SetCoordinates(nuvec2_tmpj.Pt(), nuvec2_tmpj.Eta(), nuvec2_tmpj.Phi(), M_nu2);
1909 tauvecsol2j.SetPxPyPzE(0., 0., 0., 0.);
1910 tauvecsol2j += nuvec2_tmpj;
1911 tauvecsol2j += m_tauVec2;
1912 if (tauvecsol2j.E() >= preparedInput.m_beamEnergy)
1913 continue;
1914 tauvecprob2j = Prob->apply(preparedInput, -99, preparedInput.m_type_visTau2,
1915 PtEtaPhiMVector(0, 0, 0, 0), m_tauVec2,
1916 PtEtaPhiMVector(0, 0, 0, 0), nuvec2_tmpj, false, true, false);
1917 ++ngoodsol2;
1918 }
1919 }
1920 if (tauvecprob1j == 0.)
1921 continue;
1922 if (tauvecprob2j == 0.)
1923 continue;
1924
1925 double totalProb = 1.;
1926
1927 m_tautau_tmp.SetPxPyPzE(0., 0., 0., 0.);
1928 m_tautau_tmp += tauvecsol1j;
1929 m_tautau_tmp += tauvecsol2j;
1930 const double mtautau = m_tautau_tmp.M();
1931
1932 if (TailCleanUp(m_tauVec1, nuvec1_tmpj, m_tauVec2, nuvec2_tmpj, mtautau, Mvis, Meff,
1933 preparedInput.m_DelPhiTT) == 0) {
1934 continue;
1935 }
1936
1937 totalProb *=
1938 (constProb * tauvecprob1j * tauvecprob2j *
1939 Prob->apply(preparedInput, preparedInput.m_type_visTau1, preparedInput.m_type_visTau2,
1940 m_tauVec1, m_tauVec2, nuvec1_tmpj, nuvec2_tmpj, false, false, true));
1941
1942 if (totalProb <= 0) {
1943 if (preparedInput.m_fUseVerbose)
1944 Warning("DiTauMassTools", "%s",
1945 ("null proba solution, rejected "+std::to_string(totalProb)).c_str());
1946 } else {
1947 // only count solution with non zero probability
1948 m_totalProbSum += totalProb;
1949 m_mtautauSum += mtautau;
1950
1951 if (m_nsol >= int(m_nsolfinalmax)) {
1952 Error("DiTauMassTools", "%s",
1953 ("refineSolutions ERROR nsol getting larger than nsolfinalmax!!! " +
1954 std::to_string(m_nsol))
1955 .c_str());
1956 Error("DiTauMassTools", "%s",
1957 (" j1 " + std::to_string(j1) + " j2 " + std::to_string(j2) + " nsol1 " +
1958 std::to_string(nsol1) + " nsol2 " + std::to_string(nsol2))
1959 .c_str());
1960 --m_nsol; // overwrite last solution. However this should really never
1961 // happen
1962 }
1963
1964 if ((m_nsol) < 0 or (m_nsol >= m_nsolfinalmax))[[unlikely]]{
1965 throw std::out_of_range("refineSolutions: index m_nsol out of range.");
1966 }
1967 // good solution found, copy in vector
1968 m_mtautauFinalSolVec[m_nsol] = mtautau;
1969 m_probFinalSolVec[m_nsol] = totalProb;
1970
1971 PtEtaPhiMVector &nu1Final = m_nu1FinalSolVec[m_nsol];
1972 PtEtaPhiMVector &nu2Final = m_nu2FinalSolVec[m_nsol];
1973
1974 nu1Final.SetPxPyPzE(nuvec1_tmpj.Px(), nuvec1_tmpj.Py(), nuvec1_tmpj.Pz(), nuvec1_tmpj.E());
1975 nu2Final.SetPxPyPzE(nuvec2_tmpj.Px(), nuvec2_tmpj.Py(), nuvec2_tmpj.Pz(), nuvec2_tmpj.E());
1976
1977 ++m_nsol;
1978 } // else totalProb<=0
1979
1980 } // loop j2
1981 } // loop j1
1982 if (ngoodsol1 == 0) {
1983 return -1;
1984 }
1985 if (ngoodsol2 == 0) {
1986 return -2;
1987 }
1988 return m_nsol;
1989}
1990
1991int MissingMassCalculator::TailCleanUp(const PtEtaPhiMVector &vis1,
1992 const PtEtaPhiMVector &nu1,
1993 const PtEtaPhiMVector &vis2,
1994 const PtEtaPhiMVector &nu2, const double &mmc_mass,
1995 const double &vis_mass, const double &eff_mass,
1996 const double &dphiTT) {
1997
1998 int pass_code = 1;
1999 if (preparedInput.m_fUseTailCleanup == 0)
2000 return pass_code;
2001
2002 // the Clean-up cuts are specifically for rel16 analyses.
2003 // the will change in rel17 analyses and after the MMC is updated
2004
2005 if (preparedInput.m_tauTypes == TauTypes::ll) // lepton-lepton channel
2006 {
2007 const double MrecoMvis = mmc_mass / vis_mass;
2008 if (MrecoMvis > 2.6)
2009 return 0;
2010 const double MrecoMeff = mmc_mass / eff_mass;
2011 if (MrecoMeff > 1.9)
2012 return 0;
2013 const double e1p1 = nu1.E() / vis1.P();
2014 const double e2p2 = nu2.E() / vis2.P();
2015 if ((e1p1 + e2p2) > 4.5)
2016 return 0;
2017 if (e2p2 > 4.0)
2018 return 0;
2019 if (e1p1 > 3.0)
2020 return 0;
2021 }
2022
2023 //-------- these are new cuts for lep-had analysis for Moriond
2024 if (preparedInput.m_tauTypes == TauTypes::lh) // lepton-hadron channel
2025 {
2026
2027 if (m_mmcCalibrationSet == MMCCalibrationSet::MMC2015HIGHMASS ||
2028 m_mmcCalibrationSet == MMCCalibrationSet::MMC2019 ||
2029 m_mmcCalibrationSet == MMCCalibrationSet::MMC2024 ||
2030 m_mmcCalibrationSet == MMCCalibrationSet::UPGRADE)
2031 return pass_code; // don't use TailCleanup for 8 & 13 TeV data
2032
2033 //--------- leave code uncommented to avoid Compilation warnings
2034 if (Prob->GetUseHT()) {
2035 const double MrecoMvis = mmc_mass / vis_mass;
2036 const double MrecoMeff = mmc_mass / eff_mass;
2037 const double x = dphiTT > 1.5 ? dphiTT : 1.5;
2038 if ((MrecoMeff + MrecoMvis) > 5.908 - 1.881 * x + 0.2995 * x * x)
2039 return 0;
2040 }
2041 }
2042 return pass_code;
2043}
2044
2045// note that if MarkovChain the input solutions can be modified
2046void MissingMassCalculator::handleSolutions()
2047
2048{
2049
2050 bool reject = true;
2051 double totalProbSumSol = 0.;
2052 double totalProbSumSolOld = 0.;
2053 bool firstPointWithSol = false;
2054
2055 for (int isol = 0; isol < m_nsol; ++isol) {
2056 totalProbSumSol += m_probFinalSolVec[isol];
2057 }
2058
2059 double uMC = -1.;
2060 bool notSureToKeep = true;
2061 // note : if no solution, the point is treated as having a zero probability
2062 if (m_fullParamSpaceScan) {
2063 reject = false; // accept anyway in this mode
2064 notSureToKeep = false; // do not need to test on prob
2065 if (m_nsol <= 0) {
2066 // if initial full scaning and no sol : continue
2067 m_markovNFullScan += 1;
2068 } else {
2069 // if we were in in full scan mode and we have a solution, switch it off
2070 m_fullParamSpaceScan = false;
2071 firstPointWithSol = true; // as this is the first point without a solution
2072 // there is no old sol
2073 m_iter0 = 0; // reset the counter so that separately the full scan pphase
2074 // and the markov phase use m_niterRandomLocal points
2075 // hack for hh : allow 10 times less iteration for markov than for the
2076 // fullscan phase
2077 if (preparedInput.m_tauTypes == TauTypes::hh) {
2078 m_niterRandomLocal /= 10;
2079 }
2080 }
2081 }
2082
2083 if (notSureToKeep) {
2084 // apply Metropolis algorithm to decide to keep this point.
2085 // compute the probability of the previous point and the current one
2086 for (int isol = 0; isol < m_nsolOld; ++isol) {
2087 totalProbSumSolOld += m_probFinalSolOldVec[isol];
2088 }
2089
2090 // accept anyway if null old probability (should only happen for the very
2091 // first point with a solution)
2092 if (!firstPointWithSol && totalProbSumSolOld <= 0.) {
2093 Error("DiTauMassTools", "%s",
2094 (" ERROR null old probability !!! " + std::to_string(totalProbSumSolOld) + " nsolOld " +
2095 std::to_string(m_nsolOld))
2096 .c_str());
2097 reject = false;
2098 } else if (totalProbSumSol > totalProbSumSolOld) {
2099 // if going up, accept anyway
2100 reject = false;
2101 // else if (totalProbSumSol < 1E-16) { // if null target probability,
2102 // reject anyway
2103 } else if (totalProbSumSol < totalProbSumSolOld * 1E-6) { // if ratio of probability <1e6, point
2104 // will be accepted only every 1E6
2105 // iteration, so can reject anyway
2106 reject = true;
2107 } else if (m_nsol <= 0) { // new parametrisation give prob too small to
2108 // trigger above condition if no solution is found
2109 reject = true;
2110 } else {
2111 // if going down, reject with a probability
2112 // 1-totalProbSum/totalProbSumOld)
2113 uMC = m_randomGen.Rndm();
2114 reject = (uMC > totalProbSumSol / totalProbSumSolOld);
2115 }
2116 } // if reject
2117
2118 // proceed with the handling of the solutions wether the old or the new ones
2119
2120 // optionally fill the vectors with the complete list of points (for all
2121 // walkstrategy)
2122
2123 if (reject) {
2124 // current point reset to the previous one
2125 // Note : only place where m_MEtP etc... are modified outside spacewalkerXYZ
2126 m_MEtP = m_MEtP0;
2127 m_MEtL = m_MEtL0;
2128 m_Phi1 = m_Phi10;
2129 m_Phi2 = m_Phi20;
2130 m_eTau1 = m_eTau10;
2131 m_eTau2 = m_eTau20;
2132 if (m_scanMnu1)
2133 m_Mnu1 = m_Mnu10;
2134 if (m_scanMnu2)
2135 m_Mnu2 = m_Mnu20;
2136 }
2137
2138 // default case : fill the histogram with solution, using current point
2139 bool fillSolution = true;
2140 bool oldToBeUsed = false;
2141
2142 // now handle the reject or accept cases
2143 // the tricky thing is that for markov, we accept the old point as soon as a
2144 // new accepted point is found with a weight equal to one plus the number of
2145 // rejected point inbetween
2146
2147 if (reject) {
2148 fillSolution = false; // do not fill solution, just count number of replication
2149 m_markovCountDuplicate += 1;
2150 if (m_nsol <= 0) {
2151 m_markovNRejectNoSol += 1;
2152 } else {
2153 m_markovNRejectMetropolis += 1;
2154 }
2155
2156 } else {
2157 // if accept, will fill solution (except for very first point) but taking
2158 // the values from the previous point
2159 if (!m_fullParamSpaceScan) {
2160 m_markovNAccept += 1;
2161 }
2162 if (!firstPointWithSol) {
2163 fillSolution = true;
2164 oldToBeUsed = true;
2165 } else {
2166 fillSolution = false;
2167 }
2168 } // else reject
2169
2170 // if do not fill solution exit now
2171 // for the first point with solution we need to copy the new sol into the old
2172 // one before leaving
2173 if (!fillSolution) {
2174 if (firstPointWithSol) {
2175 // current point is the future previous one
2176 m_nsolOld = m_nsol;
2177 for (int isol = 0; isol < m_nsol; ++isol) {
2178 m_probFinalSolOldVec[isol] = m_probFinalSolVec[isol];
2179 m_mtautauFinalSolOldVec[isol] = m_mtautauFinalSolVec[isol];
2180 m_nu1FinalSolOldVec[isol] = m_nu1FinalSolVec[isol];
2181 m_nu2FinalSolOldVec[isol] = m_nu2FinalSolVec[isol];
2182 }
2183 }
2184 return;
2185 }
2186
2187 // compute RMS of the different solutions
2188 double solSum = 0.;
2189 double solSum2 = 0.;
2190
2191 for (int isol = 0; isol < m_nsol; ++isol) {
2192 ++m_iter5;
2193 double totalProb{};
2194 double mtautau{};
2195 const PtEtaPhiMVector *pnuvec1_tmpj;
2196 const PtEtaPhiMVector *pnuvec2_tmpj;
2197
2198 //oldToBeUsed must be true at this point
2199 totalProb = m_probFinalSolOldVec[isol];
2200 mtautau = m_mtautauFinalSolOldVec[isol];
2201 pnuvec1_tmpj = &m_nu1FinalSolOldVec[isol];
2202 pnuvec2_tmpj = &m_nu2FinalSolOldVec[isol];
2203
2204 const PtEtaPhiMVector &nuvec1_tmpj = *pnuvec1_tmpj;
2205 const PtEtaPhiMVector &nuvec2_tmpj = *pnuvec2_tmpj;
2206
2207 solSum += mtautau;
2208 solSum2 += mtautau * mtautau;
2209
2210 double weight;
2211 // MarkovChain : accepted events already distributed according to
2212 // probability distribution, so weight is 1. acutally to have a proper
2213 // estimate of per bin error, instead of putting several time the same point
2214 // when metropolis alg reject one (or no solution), rather put it with the
2215 // multiplicity weight. Should only change the error bars might change if
2216 // weighted markov chain are used there is also an issue with the 4 very
2217 // close nearly identical solution
2218 weight = m_markovCountDuplicate +
2219 1; // incremented only when a point is rejected, hence need to add 1
2220
2221 m_fMfit_all->Fill(mtautau, weight);
2222
2223 if(m_SaveLlhHisto){
2224 m_fMEtP_all->Fill(m_MEtP, weight);
2225 m_fMEtL_all->Fill(m_MEtL, weight);
2226 m_fMnu1_all->Fill(m_Mnu1, weight);
2227 m_fMnu2_all->Fill(m_Mnu2, weight);
2228 m_fPhi1_all->Fill(m_Phi1, weight);
2229 m_fPhi2_all->Fill(m_Phi2, weight);
2230 if (mtautau != 0. && weight != 0.)
2231 m_fMfit_allGraph->SetPoint(m_iter0, mtautau, -TMath::Log(weight));
2232 }
2233
2234 m_fMfit_allNoWeight->Fill(mtautau, 1.);
2235
2236 // m_fPXfit1->Fill(nuvec1_tmpj.Px(),weight);
2237 // m_fPYfit1->Fill(nuvec1_tmpj.Py(),weight);
2238 // m_fPZfit1->Fill(nuvec1_tmpj.Pz(),weight);
2239 // m_fPXfit2->Fill(nuvec2_tmpj.Px(),weight);
2240 // m_fPYfit2->Fill(nuvec2_tmpj.Py(),weight);
2241 // m_fPZfit2->Fill(nuvec2_tmpj.Pz(),weight);
2242
2243 //----------------- using P*fit to fill Px,y,z_tau
2244 // Note that the original vistau are used there deliberately,
2245 // since they will be subtracted after histogram fitting
2246 // DR, kudos Antony Lesage : do not create temporary TLV within each Fill,
2247 // saves 10% CPU
2248 m_fPXfit1->Fill(preparedInput.m_vistau1.Px() + nuvec1_tmpj.Px(), totalProb);
2249 m_fPYfit1->Fill(preparedInput.m_vistau1.Py() + nuvec1_tmpj.Py(), totalProb);
2250 m_fPZfit1->Fill(preparedInput.m_vistau1.Pz() + nuvec1_tmpj.Pz(), totalProb);
2251 m_fPXfit2->Fill(preparedInput.m_vistau2.Px() + nuvec2_tmpj.Px(), totalProb);
2252 m_fPYfit2->Fill(preparedInput.m_vistau2.Py() + nuvec2_tmpj.Py(), totalProb);
2253 m_fPZfit2->Fill(preparedInput.m_vistau2.Pz() + nuvec2_tmpj.Pz(), totalProb);
2254
2255 // fill histograms for floating stopping criterion, split randomly
2256 if (m_fUseFloatStopping) {
2257 if (m_randomGen.Rndm() <= 0.5) {
2258 m_fMmass_split1->Fill(mtautau, weight);
2259 m_fMEtP_split1->Fill(m_MEtP, weight);
2260 m_fMEtL_split1->Fill(m_MEtL, weight);
2261 m_fMnu1_split1->Fill(m_Mnu1, weight);
2262 m_fMnu2_split1->Fill(m_Mnu2, weight);
2263 m_fPhi1_split1->Fill(m_Phi1, weight);
2264 m_fPhi2_split1->Fill(m_Phi2, weight);
2265 } else {
2266 m_fMmass_split2->Fill(mtautau, weight);
2267 m_fMEtP_split2->Fill(m_MEtP, weight);
2268 m_fMEtL_split2->Fill(m_MEtL, weight);
2269 m_fMnu1_split2->Fill(m_Mnu1, weight);
2270 m_fMnu2_split2->Fill(m_Mnu2, weight);
2271 m_fPhi1_split2->Fill(m_Phi1, weight);
2272 m_fPhi2_split2->Fill(m_Phi2, weight);
2273 }
2274 }
2275
2276 if (totalProb > m_prob_tmp) // fill solution with highest probability
2277 {
2278 m_prob_tmp = totalProb;
2279 m_fDitauStuffFit.Mditau_best = mtautau;
2280 m_fDitauStuffFit.Sign_best = -log10(totalProb);
2281 ;
2282 m_fDitauStuffFit.nutau1 = nuvec1_tmpj;
2283 m_fDitauStuffFit.nutau2 = nuvec2_tmpj;
2284 m_fDitauStuffFit.vistau1 = m_tauVec1;
2285 m_fDitauStuffFit.vistau2 = m_tauVec2;
2286 }
2287 } // loop on solutions
2288
2289 m_markovCountDuplicate = 0; // now can reset the duplicate count
2290
2291 if (oldToBeUsed) {
2292 // current point is the future previous one
2293 // TLV copy not super efficient but not dramatic
2294 m_nsolOld = m_nsol;
2295 for (int isol = 0; isol < m_nsol; ++isol) {
2296 m_probFinalSolOldVec[isol] = m_probFinalSolVec[isol];
2297 m_mtautauFinalSolOldVec[isol] = m_mtautauFinalSolVec[isol];
2298 m_nu1FinalSolOldVec[isol] = m_nu1FinalSolVec[isol];
2299 m_nu2FinalSolOldVec[isol] = m_nu2FinalSolVec[isol];
2300 }
2301 }
2302 if (m_nsol == 0) [[unlikely]]{
2303 //throw'ing here causes a ctest to fail; should be investigated
2304 return;
2305 }
2306 // compute rms of solutions
2307 const double solRMS = sqrt(solSum2 / m_nsol - std::pow(solSum / m_nsol, 2));
2308 OutputInfo.m_AveSolRMS += solRMS;
2309
2310 return;
2311}
2312
2313void MissingMassCalculator::SpaceWalkerInit() {
2314 // FIXME could use function pointer to switch between functions
2315 m_nsolOld = 0;
2316
2317 double METresX = preparedInput.m_METsigmaL; // MET resolution in direction parallel to MET
2318 // resolution major axis, for MET scan
2319 double METresY = preparedInput.m_METsigmaP; // MET resolution in direction perpendicular to
2320 // to MET resolution major axis, for MET scan
2321
2322 // precompute some quantities and store in m_ data members
2323 precomputeCache();
2324 if (m_mmcCalibrationSet == MMCCalibrationSet::MMC2019 || m_mmcCalibrationSet == MMCCalibrationSet::MMC2024) {
2325 if (Prob->GetUseMnuProbability() == true && (preparedInput.m_tauTypes == TauTypes::ll || preparedInput.m_tauTypes == TauTypes::lh) ) Prob->setParamNuMass();
2326 Prob->setParamAngle(m_tauVec1, 1, preparedInput.m_type_visTau1);
2327 Prob->setParamAngle(m_tauVec2, 2, preparedInput.m_type_visTau2);
2328 Prob->setParamRatio(1, preparedInput.m_type_visTau1);
2329 Prob->setParamRatio(2, preparedInput.m_type_visTau2);
2330 }
2331
2332 // if m_nsigma_METscan was not set by user, set to default values
2333 if(m_nsigma_METscan == -1){
2334 if (preparedInput.m_tauTypes == TauTypes::ll) // both tau's are leptonic
2335 {
2336 m_nsigma_METscan = m_nsigma_METscan_ll;
2337 } else if (preparedInput.m_tauTypes == TauTypes::lh) // lep had
2338 {
2339 m_nsigma_METscan = m_nsigma_METscan_lh;
2340 } else // hh
2341 {
2342 m_nsigma_METscan = m_nsigma_METscan_hh;
2343 }
2344 }
2345
2346 m_nsigma_METscan2 = std::pow(m_nsigma_METscan, 2);
2347
2348 const double deltaPhi1 = MaxDelPhi(preparedInput.m_type_visTau1, m_tauVec1P, m_dRmax_tau);
2349 const double deltaPhi2 = MaxDelPhi(preparedInput.m_type_visTau2, m_tauVec2P, m_dRmax_tau);
2350
2351 m_walkWeight = 1.;
2352
2353 // dummy initial value to avoid printout with random values
2354 m_Phi10 = 0.;
2355 m_Phi20 = 0.;
2356 m_MEtL0 = 0.;
2357 m_MEtP0 = 0.;
2358 m_Mnu10 = 0.;
2359 m_Mnu20 = 0.;
2360
2361 m_mTau = ParticleConstants::tauMassInMeV / GEV;
2362
2363 // seeds the random generator in a reproducible way from the phi of both tau;
2364 double aux = std::abs(m_tauVec1Phi + double(m_tauVec2Phi) / 100. / TMath::Pi()) * 100;
2365 m_seed = (aux - floor(aux)) * 1E6 * (1 + m_RndmSeedAltering) + 13;
2366
2367 m_randomGen.SetSeed(m_seed);
2368 // int Niter=Niter_fit1; // number of points for each dR loop
2369 // int NiterMET=Niter_fit2; // number of iterations for each MET scan loop
2370 // int NiterMnu=Niter_fit3; // number of iterations for Mnu loop
2371
2372 // approximately compute the number of points from the grid scanning
2373 // divide by abritry number to recover timing with still better results
2374 // m_NiterRandom=(NiterMET+1)*(NiterMET+1)*4*Niter*Niter/10;
2375
2376 m_Phi1Min = m_tauVec1Phi - deltaPhi1;
2377 m_Phi1Max = m_tauVec1Phi + deltaPhi1;
2378 m_Phi1Range = m_Phi1Max - m_Phi1Min;
2379
2380 m_Phi2Min = m_tauVec2Phi - deltaPhi2;
2381 m_Phi2Max = m_tauVec2Phi + deltaPhi2;
2382 m_Phi2Range = m_Phi2Max - m_Phi2Min;
2383
2384 m_Mnu1Min = 0.;
2385 m_scanMnu1 = false;
2386 m_Mnu1 = m_Mnu1Min;
2387
2388 // for markov chain use factor 2
2389 m_niterRandomLocal = m_NiterRandom / 2;
2390
2391 // NiterRandom set by user (default is -1). If negative, defines the default
2392 // here. no more automatic scaling for ll hl hh
2393 if (m_NiterRandom <= 0) {
2394 m_niterRandomLocal = 100000; // number of iterations for Markov for lh
2395 if (preparedInput.m_tauTypes == TauTypes::ll)
2396 m_niterRandomLocal *= 2; // multiplied for ll , unchecked
2397 if (preparedInput.m_tauTypes == TauTypes::hh)
2398 m_niterRandomLocal *= 5; // divided for hh ,checked
2399 } else {
2400 m_niterRandomLocal = m_NiterRandom;
2401 }
2402
2403 if (preparedInput.m_type_visTau1 == 8) {
2404 // m_Mnu1Max=m_mTau-m_tauVec1M;
2405 m_Mnu1Max = m_MnuScanRange - m_tauVec1M;
2406 m_Mnu1Range = m_Mnu1Max - m_Mnu1Min;
2407 m_scanMnu1 = true;
2408 }
2409
2410 m_Mnu2Min = 0.;
2411 m_scanMnu2 = false;
2412 m_Mnu2 = m_Mnu2Min;
2413 if (preparedInput.m_type_visTau2 == 8) {
2414 // m_Mnu2Max=m_mTau-m_tauVec2M;
2415 m_Mnu2Max = m_MnuScanRange - m_tauVec2M;
2416 m_Mnu2Range = m_Mnu2Max - m_Mnu2Min;
2417 m_scanMnu2 = true;
2418 }
2419
2420 m_MEtLMin = -m_nsigma_METscan * METresX;
2421 m_MEtLMax = +m_nsigma_METscan * METresX;
2422 m_MEtLRange = m_MEtLMax - m_MEtLMin;
2423
2424 m_MEtPMin = -m_nsigma_METscan * METresY;
2425 m_MEtPMax = +m_nsigma_METscan * METresY;
2426 m_MEtPRange = m_MEtPMax - m_MEtPMin;
2427
2428 m_eTau1Min = -1;
2429 m_eTau1Max = -1;
2430 m_eTau2Min = -1;
2431 m_eTau2Max = -1;
2432
2433 m_switch1 = true;
2434 m_switch2 = true;
2435
2436 m_nsucStop = m_NsucStop;
2437 m_rmsStop = m_RMSStop;
2438
2439 m_iter0 = -1;
2440 m_iterNuPV3 = 0;
2441 m_testptn1 = 0;
2442 m_testptn2 = 0;
2443 m_testdiscri1 = 0;
2444 m_testdiscri2 = 0;
2445 m_nosol1 = 0;
2446 m_nosol2 = 0;
2447 m_iterNsuc = 0;
2448 if (m_meanbinStop > 0) {
2449 m_meanbinToBeEvaluated = true;
2450 } else {
2451 m_meanbinToBeEvaluated = false;
2452 }
2453
2454 m_markovCountDuplicate = 0;
2455 m_markovNRejectNoSol = 0;
2456 m_markovNRejectMetropolis = 0;
2457 m_markovNAccept = 0;
2458 m_markovNFullScan = 0;
2459 // set full parameter space scannning for the first steps, until a solution is
2460 // found
2461 m_fullParamSpaceScan = true;
2462 // size of step. Needs to be tune. Start with simple heuristic.
2463 if (m_proposalTryMEt < 0) {
2464 m_MEtProposal = m_MEtPRange / 30.;
2465 } else {
2466 m_MEtProposal = m_MEtPRange * m_proposalTryMEt;
2467 }
2468 if (m_ProposalTryPhi < 0) {
2469 m_PhiProposal = 0.04;
2470 } else {
2471 m_PhiProposal = m_ProposalTryPhi;
2472 }
2473 // FIXME if m_Mnu1Range !ne m_Mnu2Range same proposal will be done
2474 if (m_scanMnu1) {
2475 if (m_ProposalTryMnu < 0) {
2476 m_MnuProposal = m_Mnu1Range / 10.;
2477 } else {
2478 m_MnuProposal = m_Mnu1Range * m_ProposalTryMnu;
2479 }
2480 }
2481 if (m_scanMnu2) {
2482 if (m_ProposalTryMnu < 0) {
2483 m_MnuProposal = m_Mnu2Range / 10.;
2484 } else {
2485 m_MnuProposal = m_Mnu2Range * m_ProposalTryMnu;
2486 }
2487 }
2488}
2489
2490// iterator. walk has internal counters, should only be used in a while loop
2491// so far only implement grid strategy
2492// act on MMC data member to be fast
2493bool MissingMassCalculator::SpaceWalkerWalk() {
2494 preparedInput.m_MEtX = -999.;
2495 preparedInput.m_MEtY = -999.;
2496
2497 ++m_iter0;
2498
2499 if (m_meanbinToBeEvaluated && m_iterNsuc == 500) {
2500 Info("DiTauMassTools", " in m_meanbinToBeEvaluated && m_iterNsuc==500 ");
2501 // for markov chain m_iterNsuc is the number of *accepted* points, so there
2502 // can be several iterations without any increment of m_iterNsuc. Hence need
2503 // to make sure meanbin is evaluated only once
2504 m_meanbinToBeEvaluated = false;
2505
2506 // Meanbin stopping criterion
2507 std::vector<double> histInfo(HistInfo::MAXHISTINFO);
2508 // SLIDINGWINDOW strategy to avoid doing the parabola fit now given it will
2509 // not be use
2510 maxFromHist(m_fMfit_all, histInfo, MaxHistStrategy::SLIDINGWINDOW);
2511 double meanbin = histInfo.at(HistInfo::MEANBIN);
2512 if (meanbin < 0) {
2513 m_nsucStop = -1; // no meaningful meanbin switch back to niter criterion
2514 } else {
2515 double stopdouble = 500 * std::pow((meanbin / m_meanbinStop), 2);
2516 int stopint = stopdouble;
2517 m_nsucStop = stopint;
2518 }
2519 if (m_nsucStop < 500)
2520 return false;
2521 }
2522 // should be outside m_meanbinStop test
2523 if (m_iterNsuc == m_nsucStop)
2524 return false; // Critere d'arret pour nombre de succes
2525
2526 if (m_iter0 == m_niterRandomLocal)
2527 return false; // for now simple stopping criterion on number of iteration
2528
2529 // floating stopping criterion, reduces run-time for lh, hh by a factor ~2 and ll by roughly
2530 // factor ~3 check if every scanned variable and resulting mass thermalised after N (default 10k) iterations
2531 // and then every M (default 1k) iterations do this by checking that the means of the split distributions is
2532 // comparable within X% (default 5%) of their sigma
2533 if (m_iter0 >= m_fUseFloatStoppingMinIter && (m_iter0 % m_fUseFloatStoppingCheckFreq) == 0 && m_fUseFloatStopping) {
2534 if (std::abs(m_fMEtP_split1->GetMean() - m_fMEtP_split2->GetMean()) <= m_fUseFloatStoppingComp * m_fMEtP_split1->GetRMS()) {
2535 if (std::abs(m_fMEtL_split1->GetMean() - m_fMEtL_split2->GetMean()) <=
2536 m_fUseFloatStoppingComp * m_fMEtL_split1->GetRMS()) {
2537 if (std::abs(m_fMnu1_split1->GetMean() - m_fMnu1_split2->GetMean()) <=
2538 m_fUseFloatStoppingComp * m_fMnu1_split1->GetRMS()) {
2539 if (std::abs(m_fMnu2_split1->GetMean() - m_fMnu2_split2->GetMean()) <=
2540 m_fUseFloatStoppingComp * m_fMnu2_split1->GetRMS()) {
2541 if (std::abs(m_fPhi1_split1->GetMean() - m_fPhi1_split2->GetMean()) <=
2542 m_fUseFloatStoppingComp * m_fPhi1_split1->GetRMS()) {
2543 if (std::abs(m_fPhi2_split1->GetMean() - m_fPhi2_split2->GetMean()) <=
2544 m_fUseFloatStoppingComp * m_fPhi2_split1->GetRMS()) {
2545 if (std::abs(m_fMmass_split1->GetMean() - m_fMmass_split2->GetMean()) <=
2546 m_fUseFloatStoppingComp * m_fMmass_split1->GetRMS()) {
2547 return false;
2548 }
2549 }
2550 }
2551 }
2552 }
2553 }
2554 }
2555 }
2556
2557 if (m_fullParamSpaceScan) {
2558 // as long as no solution found need to randomise on the full parameter
2559 // space
2560
2561 // cut the corners in MissingET (not optimised at all)
2562 // not needed if distribution is already gaussian
2563 do {
2564 m_MEtP = m_MEtPMin + m_MEtPRange * m_randomGen.Rndm();
2565 m_MEtL = m_MEtLMin + m_MEtLRange * m_randomGen.Rndm();
2566 } while (!checkMEtInRange());
2567
2568 if (m_scanMnu1) {
2569 m_Mnu1 = m_Mnu1Min + m_Mnu1Range * m_randomGen.Rndm();
2570 }
2571
2572 if (m_scanMnu2) {
2573 m_Mnu2 = m_Mnu2Min + m_Mnu2Range * m_randomGen.Rndm();
2574 }
2575
2576 m_Phi1 = m_Phi1Min + m_Phi1Range * m_randomGen.Rndm();
2577 m_Phi2 = m_Phi2Min + m_Phi2Range * m_randomGen.Rndm();
2578
2579 return true;
2580 }
2581
2582 // here the real markov chain takes place : "propose" the new point
2583 // note that if one parameter goes outside range, this should not be fixed
2584 // here but later in handleSolution, otherwise would cause a bias
2585
2586 // m_MEtP0 etc... also store the position of the previous Markov Chain step,
2587 // which is needed by the algorithm
2588 m_MEtP0 = m_MEtP;
2589 m_MEtL0 = m_MEtL;
2590
2591 m_MEtP = m_randomGen.Gaus(m_MEtP0, m_MEtProposal);
2592
2593 m_MEtL = m_randomGen.Gaus(m_MEtL0, m_MEtProposal);
2594
2595 if (m_scanMnu1) {
2596 m_Mnu10 = m_Mnu1;
2597 m_Mnu1 = m_randomGen.Gaus(m_Mnu10, m_MnuProposal);
2598 }
2599
2600 if (m_scanMnu2) {
2601 m_Mnu20 = m_Mnu2;
2602 m_Mnu2 = m_randomGen.Gaus(m_Mnu20, m_MnuProposal);
2603 }
2604
2605 m_Phi10 = m_Phi1;
2606 m_Phi1 = m_randomGen.Gaus(m_Phi10, m_PhiProposal);
2607
2608 m_Phi20 = m_Phi2;
2609
2610 m_Phi2 = m_randomGen.Gaus(m_Phi20, m_PhiProposal);
2611
2612 return true;
2613}
2614
2615// compute cached values (this value do not change within one call of MMC,
2616// except for tau e scanning) return true if cache was already uptodatexs
2617inline bool MissingMassCalculator::precomputeCache() {
2618
2619 // copy tau 4 vect. If tau E scanning, these vectors will be modified
2620 m_tauVec1 = preparedInput.m_vistau1;
2621 m_tauVec2 = preparedInput.m_vistau2;
2622
2623 const XYVector &metVec = preparedInput.m_MetVec;
2624
2625 bool same = true;
2626 same = updateDouble(m_tauVec1.Phi(), m_tauVec1Phi) && same;
2627 same = updateDouble(m_tauVec2.Phi(), m_tauVec2Phi) && same;
2628 same = updateDouble(m_tauVec1.M(), m_tauVec1M) && same;
2629 same = updateDouble(m_tauVec2.M(), m_tauVec2M) && same;
2630 same = updateDouble(m_tauVec1.E(), m_tauVec1E) && same;
2631 same = updateDouble(m_tauVec2.E(), m_tauVec2E) && same;
2632 same = updateDouble(m_tauVec1.Px(), m_tauVec1Px) && same;
2633 same = updateDouble(m_tauVec1.Py(), m_tauVec1Py) && same;
2634 same = updateDouble(m_tauVec1.Pz(), m_tauVec1Pz) && same;
2635 same = updateDouble(m_tauVec2.Px(), m_tauVec2Px) && same;
2636 same = updateDouble(m_tauVec2.Py(), m_tauVec2Py) && same;
2637 same = updateDouble(m_tauVec2.Pz(), m_tauVec2Pz) && same;
2638 same = updateDouble(m_tauVec1.P(), m_tauVec1P) && same;
2639 same = updateDouble(m_tauVec2.P(), m_tauVec2P) && same;
2640
2641 same = updateDouble(ParticleConstants::tauMassInMeV / GEV, m_mTau) && same;
2642 same = updateDouble(std::pow(m_mTau, 2), m_mTau2) && same;
2643 same = updateDouble(cos(preparedInput.m_METcovphi), m_metCovPhiCos) && same;
2644 same = updateDouble(sin(preparedInput.m_METcovphi), m_metCovPhiSin) && same;
2645 same = updateDouble((m_tauVec1 + m_tauVec2).M(), m_Mvis) && same;
2646
2647 PtEtaPhiMVector Met4vec;
2648 Met4vec.SetPxPyPzE(preparedInput.m_MetVec.X(), preparedInput.m_MetVec.Y(), 0.0,
2649 preparedInput.m_MetVec.R());
2650 same = updateDouble((m_tauVec1 + m_tauVec2 + Met4vec).M(), m_Meff) && same;
2651
2652 same = updateDouble(preparedInput.m_HtOffset, preparedInput.m_htOffset) && same;
2653 // note that if useHT met_vec is actually -HT
2654 same = updateDouble(metVec.X(), preparedInput.m_inputMEtX) && same;
2655 same = updateDouble(metVec.Y(), preparedInput.m_inputMEtY) && same;
2656 same = updateDouble(metVec.R(), preparedInput.m_inputMEtT) && same;
2657
2658 return same;
2659}
2660
2661// return true if all parameters are within their domain
2662inline bool MissingMassCalculator::checkAllParamInRange() {
2663
2664 if (m_scanMnu1) {
2665 if (m_Mnu1 < m_Mnu1Min)
2666 return false;
2667 if (m_Mnu1 > m_Mnu1Max)
2668 return false;
2669 if (m_Mnu1 > m_mTau - m_tauVec1M)
2670 return false;
2671 }
2672
2673 if (m_scanMnu2) {
2674 if (m_Mnu2 < m_Mnu2Min)
2675 return false;
2676 if (m_Mnu2 > m_Mnu2Max)
2677 return false;
2678 if (m_Mnu2 > m_mTau - m_tauVec2M)
2679 return false;
2680 }
2681
2682 // FIXME note that since there is a coupling between Met and tau, should
2683 // rigorously test both together however since the 3 sigma range is just a
2684 // hack, it is probably OK
2685
2686 if (m_Phi1 < m_Phi1Min)
2687 return false;
2688 if (m_Phi1 > m_Phi1Max)
2689 return false;
2690
2691 if (m_Phi2 < m_Phi2Min)
2692 return false;
2693 if (m_Phi2 > m_Phi2Max)
2694 return false;
2695
2696 if (!checkMEtInRange())
2697 return false;
2698
2699 return true;
2700}
2701
2702// return true if Met is within disk instead of withing square (cut the corners)
2703inline bool MissingMassCalculator::checkMEtInRange() {
2704 // check MEt is in allowed range
2705 // range is 3sigma disk ("cutting the corners")
2706 if (std::pow(m_MEtL / preparedInput.m_METsigmaL, 2) +
2707 std::pow(m_MEtP / preparedInput.m_METsigmaP, 2) >
2708 m_nsigma_METscan2) {
2709 return false;
2710 } else {
2711 return true;
2712 }
2713}
2714
2715// ----- returns dTheta3D lower and upper boundaries:
2716// limit_code=0: 99% lower limit
2717// limit_code=1; 99% upper limit
2718// limit_code=2; 95% upper limit
2719double MissingMassCalculator::dTheta3DLimit(const int &tau_type, const int &limit_code,
2720 const double &P_tau) {
2721
2722#ifndef WITHDTHETA3DLIM
2723 // make the test ineffective if desired
2724 if (limit_code == 0)
2725 return 0.;
2726 if (limit_code == 1)
2727 return 10.;
2728 if (limit_code == 2)
2729 return 10.;
2730#endif
2731
2732 double limit = 1.0;
2733 // cppcheck-suppress identicalConditionAfterEarlyExit; in #ifdef above
2734 if (limit_code == 0)
2735 limit = 0.0;
2736 double par[3] = {0.0, 0.0, 0.0};
2737 // ---- leptonic tau's
2738 if (tau_type == 8) {
2739 if (limit_code == 0) // lower 99% limit
2740 {
2741 par[0] = 0.3342;
2742 par[1] = -0.3376;
2743 par[2] = -0.001377;
2744 }
2745 if (limit_code == 1) // upper 99% limit
2746 {
2747 par[0] = 3.243;
2748 par[1] = -12.87;
2749 par[2] = 0.009656;
2750 }
2751 if (limit_code == 2) // upper 95% limit
2752 {
2753 par[0] = 2.927;
2754 par[1] = -7.911;
2755 par[2] = 0.007783;
2756 }
2757 }
2758 // ---- 1-prong tau's
2759 if (tau_type >= 0 && tau_type <= 2) {
2760 if (limit_code == 0) // lower 99% limit
2761 {
2762 par[0] = 0.2673;
2763 par[1] = -14.8;
2764 par[2] = -0.0004859;
2765 }
2766 if (limit_code == 1) // upper 99% limit
2767 {
2768 par[0] = 9.341;
2769 par[1] = -15.88;
2770 par[2] = 0.0333;
2771 }
2772 if (limit_code == 2) // upper 95% limit
2773 {
2774 par[0] = 6.535;
2775 par[1] = -8.649;
2776 par[2] = 0.00277;
2777 }
2778 }
2779 // ---- 3-prong tau's
2780 if (tau_type >= 3 && tau_type <= 5) {
2781 if (limit_code == 0) // lower 99% limit
2782 {
2783 par[0] = 0.2308;
2784 par[1] = -15.24;
2785 par[2] = -0.0009458;
2786 }
2787 if (limit_code == 1) // upper 99% limit
2788 {
2789 par[0] = 14.58;
2790 par[1] = -6.043;
2791 par[2] = -0.00928;
2792 }
2793 if (limit_code == 2) // upper 95% limit
2794 {
2795 par[0] = 8.233;
2796 par[1] = -0.3018;
2797 par[2] = -0.009399;
2798 }
2799 }
2800
2801 if (std::abs(P_tau + par[1]) > 0.0)
2802 limit = par[0] / (P_tau + par[1]) + par[2];
2803 if (limit_code == 0) {
2804 if (limit < 0.0) {
2805 limit = 0.0;
2806 } else if (limit > 0.03) {
2807 limit = 0.03;
2808 }
2809 } else {
2810 if (limit < 0.0 || limit > 0.5 * TMath::Pi()) {
2811 limit = 0.5 * TMath::Pi();
2812 } else if (limit < 0.05 && limit > 0.0) {
2813 limit = 0.05; // parameterization only runs up to P~220 GeV in this regime
2814 // will set an upper bound of 0.05
2815 }
2816 }
2817
2818 return limit;
2819}
2820
2821// checks units of input variables, converts into [GeV] if needed, make all
2822// possible corrections DR new : now a second structure preparedInput is derived
2823// from the input one which only has direct user input
2824void MissingMassCalculator::FinalizeSettings(const xAOD::IParticle *part1,
2825 const xAOD::IParticle *part2,
2826 const xAOD::MissingET *met,
2827 const int &njets) {
2828 int mmcType1 = mmcType(part1);
2829 if (mmcType1 < 0)
2830 return; // return CP::CorrectionCode::Error;
2831
2832 int mmcType2 = mmcType(part2);
2833 if (mmcType2 < 0)
2834 return; // return CP::CorrectionCode::Error;
2835
2836 preparedInput.SetLFVmode(-2); // initialise LFV mode value for this event with being *not* LFV
2837 // if(getLFVMode(part1, part2, mmcType1, mmcType2) ==
2838 // CP::CorrectionCode::Error) {
2839 if (m_mmcCalibrationSet == MMCCalibrationSet::LFVMMC2012) {
2840 int LFVMode = getLFVMode(part1, part2, mmcType1, mmcType2);
2841 if (LFVMode == -1) {
2842 return; // return CP::CorrectionCode::Error;
2843 } else if (LFVMode != -2) {
2844 preparedInput.SetLFVmode(LFVMode);
2845 }
2846 }
2847
2848 // this will be in MeV but MMC allows MeV
2849 // assume the mass is correct as well
2850 PtEtaPhiMVector tlvTau1(part1->pt(), part1->eta(), part1->phi(), part1->m());
2851 PtEtaPhiMVector tlvTau2(part2->pt(), part2->eta(), part2->phi(), part2->m());
2852
2853 // Convert to GeV. In principle, MMC should cope with MeV but should check
2854 // thoroughly
2855 PtEtaPhiMVector fixedtau1;
2856 fixedtau1.SetCoordinates(tlvTau1.Pt() / GEV, tlvTau1.Eta(), tlvTau1.Phi(), tlvTau1.M() / GEV);
2857 PtEtaPhiMVector fixedtau2;
2858 fixedtau2.SetCoordinates(tlvTau2.Pt() / GEV, tlvTau2.Eta(), tlvTau2.Phi(), tlvTau2.M() / GEV);
2859
2860 preparedInput.SetVisTauType(0, mmcType1);
2861 preparedInput.SetVisTauType(1, mmcType2);
2862 preparedInput.SetVisTauVec(0, fixedtau1);
2863 preparedInput.SetVisTauVec(1, fixedtau2);
2864
2865 if (mmcType1 == 8 && mmcType2 == 8) {
2866 preparedInput.m_tauTypes = TauTypes::ll;
2867 } else if (mmcType1 >= 0 && mmcType1 <= 5 && mmcType2 >= 0 && mmcType2 <= 5) {
2868 preparedInput.m_tauTypes = TauTypes::hh;
2869 } else {
2870 preparedInput.m_tauTypes = TauTypes::lh;
2871 }
2872 if (preparedInput.m_fUseVerbose)
2873 Info("DiTauMassTools", "%s", ("running for tau types "+std::to_string(preparedInput.m_type_visTau1)+" "+std::to_string(preparedInput.m_type_visTau2)).c_str());
2874 XYVector met_vec(met->mpx() / GEV, met->mpy() / GEV);
2875 preparedInput.SetMetVec(met_vec);
2876 if (preparedInput.m_fUseVerbose)
2877 Info("DiTauMassTools", "%s", ("passing SumEt="+std::to_string(met->sumet() / GEV)).c_str());
2878 preparedInput.SetSumEt(met->sumet() / GEV);
2879 preparedInput.SetNjet25(njets);
2880
2881 // check that the calibration set has been chosen explicitly, otherwise abort
2882 if (m_mmcCalibrationSet == MMCCalibrationSet::MAXMMCCALIBRATIONSET) {
2883 Error("DiTauMassTools", "MMCCalibrationSet has not been set !. Please use "
2884 "fMMC.SetCalibrationSet(MMCCalibrationSet::MMC2019) or fMMC.SetCalibrationSet(MMCCalibrationSet::MMC2024)"
2885 ". Abort now. ");
2886 std::abort();
2887 }
2888 //----------- Re-ordering input info, to make sure there is no dependence of
2889 // results on input order
2890 // this might be needed because a random scan is used
2891 // highest pT tau is always first
2892 preparedInput.m_InputReorder = 0; // set flag to 0 by default, i.e. no re-ordering
2893 if ((preparedInput.m_type_visTau1 >= 0 && preparedInput.m_type_visTau1 <= 5) &&
2894 preparedInput.m_type_visTau2 == 8) // if hadron-lepton, reorder to have lepton first
2895 {
2896 preparedInput.m_InputReorder =
2897 1; // re-order to be done, this flag is to be checked in DoOutputInfo()
2898 } else if (!((preparedInput.m_type_visTau2 >= 0 && preparedInput.m_type_visTau2 <= 5) &&
2899 preparedInput.m_type_visTau1 == 8)) // if not lep-had nor had lep, reorder if tau1 is
2900 // after tau2 clockwise
2901 {
2902 if (fixPhiRange(preparedInput.m_vistau1.Phi() - preparedInput.m_vistau2.Phi()) > 0) {
2903 preparedInput.m_InputReorder = 1; // re-order to be done, this flag is to be
2904 // checked in DoOutputInfo()
2905 }
2906 }
2907
2908 if (preparedInput.m_InputReorder == 1) // copy and re-order
2909 {
2910 std::swap(preparedInput.m_vistau1, preparedInput.m_vistau2);
2911 std::swap(preparedInput.m_type_visTau1, preparedInput.m_type_visTau2);
2912 std::swap(preparedInput.m_Nprong_tau1, preparedInput.m_Nprong_tau2);
2913 }
2914 //--------- re-ordering is done ---------------------------------------
2915
2916 preparedInput.m_DelPhiTT =
2917 std::abs(Phi_mpi_pi(preparedInput.m_vistau1.Phi() - preparedInput.m_vistau2.Phi()));
2918
2919 for (unsigned int i = 0; i < preparedInput.m_jet4vecs.size(); i++) {
2920 // correcting sumEt, give priority to SetMetScanParamsUE()
2921 if (preparedInput.m_METScanScheme == 0) {
2922 if ((preparedInput.m_METsigmaP < 0.1 || preparedInput.m_METsigmaL < 0.1) &&
2923 preparedInput.m_SumEt > preparedInput.m_jet4vecs[i].Pt() &&
2924 preparedInput.m_jet4vecs[i].Pt() > 20.0) {
2925 if (preparedInput.m_fUseVerbose == 1) {
2926 Info("DiTauMassTools", "correcting sumET");
2927 }
2928 preparedInput.m_SumEt -= preparedInput.m_jet4vecs[i].Pt();
2929 }
2930 }
2931 }
2932
2933 // give priority to SetVisTauType, only do this if type_visTau1 and
2934 // type_visTau2 are not set
2935 /*if(type_visTau1<0 && type_visTau2<0 && Nprong_tau1>-1 && Nprong_tau2>-1)
2936 {
2937 if(Nprong_tau1==0) type_visTau1 = 8; // leptonic tau
2938 else if( Nprong_tau1==1) type_visTau1 = 0; // set to 1p0n for now, may use
2939different solution later like explicit integer for this case that pantau info is
2940not available? else if( Nprong_tau1==3) type_visTau1 = 3; // set to 3p0n for
2941now, see above if(Nprong_tau2==0) type_visTau2 = 8; // leptonic tau else if(
2942Nprong_tau2==1) type_visTau2 = 0; // set to 1p0n for now, see above else if(
2943Nprong_tau2==3) type_visTau2=3; // set to 3p0n for now, see above
2944 }
2945 */
2946 // checking input mass of hadronic tau-1
2947 // one prong
2948 // // checking input mass of hadronic tau-1
2949 // DRMERGE LFV addition
2950 if (m_mmcCalibrationSet == MMCCalibrationSet::LFVMMC2012 ) {
2951 if ((preparedInput.m_type_visTau1 >= 0 && preparedInput.m_type_visTau1 <= 2) &&
2952 preparedInput.m_vistau1.M() != 1.1) {
2953 preparedInput.m_vistau1.SetCoordinates(preparedInput.m_vistau1.Pt(), preparedInput.m_vistau1.Eta(),
2954 preparedInput.m_vistau1.Phi(), 1.1);
2955 }
2956 if ((preparedInput.m_type_visTau1 >= 3 && preparedInput.m_type_visTau1 <= 5) &&
2957 preparedInput.m_vistau1.M() != 1.35) {
2958 preparedInput.m_vistau1.SetCoordinates(preparedInput.m_vistau1.Pt(), preparedInput.m_vistau1.Eta(),
2959 preparedInput.m_vistau1.Phi(), 1.35);
2960 }
2961 // checking input mass of hadronic tau-2
2962 if ((preparedInput.m_type_visTau2 >= 0 && preparedInput.m_type_visTau2 <= 2) &&
2963 preparedInput.m_vistau2.M() != 1.1) {
2964 preparedInput.m_vistau2.SetCoordinates(preparedInput.m_vistau2.Pt(), preparedInput.m_vistau2.Eta(),
2965 preparedInput.m_vistau2.Phi(), 1.1);
2966 }
2967 if ((preparedInput.m_type_visTau2 >= 3 && preparedInput.m_type_visTau2 <= 5) &&
2968 preparedInput.m_vistau2.M() != 1.35) {
2969 preparedInput.m_vistau2.SetCoordinates(preparedInput.m_vistau2.Pt(), preparedInput.m_vistau2.Eta(),
2970 preparedInput.m_vistau2.Phi(), 1.35);
2971 }
2972 } else {
2973 // DRMERGE end LFV addition
2974 if ((preparedInput.m_type_visTau1 >= 0 && preparedInput.m_type_visTau1 <= 2) &&
2975 preparedInput.m_vistau1.M() != 0.8) {
2976 preparedInput.m_vistau1.SetCoordinates(preparedInput.m_vistau1.Pt(), preparedInput.m_vistau1.Eta(),
2977 preparedInput.m_vistau1.Phi(), 0.8);
2978 }
2979 // 3 prong
2980 if ((preparedInput.m_type_visTau1 >= 3 && preparedInput.m_type_visTau1 <= 5) &&
2981 preparedInput.m_vistau1.M() != 1.2) {
2982 preparedInput.m_vistau1.SetCoordinates(preparedInput.m_vistau1.Pt(), preparedInput.m_vistau1.Eta(),
2983 preparedInput.m_vistau1.Phi(), 1.2);
2984 }
2985 // checking input mass of hadronic tau-2
2986 // one prong
2987 if ((preparedInput.m_type_visTau2 >= 0 && preparedInput.m_type_visTau2 <= 2) &&
2988 preparedInput.m_vistau2.M() != 0.8) {
2989 preparedInput.m_vistau2.SetCoordinates(preparedInput.m_vistau2.Pt(), preparedInput.m_vistau2.Eta(),
2990 preparedInput.m_vistau2.Phi(), 0.8);
2991 }
2992 // 3 prong
2993 if ((preparedInput.m_type_visTau2 >= 3 && preparedInput.m_type_visTau2 <= 5) &&
2994 preparedInput.m_vistau2.M() != 1.2) {
2995 preparedInput.m_vistau2.SetCoordinates(preparedInput.m_vistau2.Pt(), preparedInput.m_vistau2.Eta(),
2996 preparedInput.m_vistau2.Phi(), 1.2);
2997 }
2998 } // DRDRMERGE LFV else closing
2999
3000 // correcting sumEt for electron pt, give priority to SetMetScanParamsUE()
3001 // DR20150615 in tag 00-00-11 and before. The following was done before the
3002 // mass of the hadronic tau was set which mean that sumEt was wrongly
3003 // corrected for the hadronic tau pt if the hadronic tau mass was set to zero
3004 // Sasha 08/12/15: don't do electron Pt subtraction for high mass studies; in
3005 // the future, need to check if lepton Pt needs to be subtracted for both ele
3006 // and muon
3007 if (preparedInput.m_METsigmaP < 0.1 || preparedInput.m_METsigmaL < 0.1) {
3008
3009 // T. Davidek: hack for lep-lep -- subtract lepton pT both for muon and
3010 // electron
3011 if ((m_mmcCalibrationSet == MMCCalibrationSet::MMC2019 ||
3012 m_mmcCalibrationSet == MMCCalibrationSet::MMC2024) &&
3013 preparedInput.m_vistau1.M() < 0.12 && preparedInput.m_vistau2.M() < 0.12) { // lep-lep channel
3014 if (preparedInput.m_SumEt > preparedInput.m_vistau1.Pt())
3015 preparedInput.m_SumEt -= preparedInput.m_vistau1.Pt();
3016 if (preparedInput.m_SumEt > preparedInput.m_vistau2.Pt())
3017 preparedInput.m_SumEt -= preparedInput.m_vistau2.Pt();
3018 } else {
3019 // continue with the original code
3020 if (preparedInput.m_SumEt > preparedInput.m_vistau1.Pt() && preparedInput.m_vistau1.M() < 0.05 &&
3021 m_mmcCalibrationSet != MMCCalibrationSet::MMC2015HIGHMASS) {
3022 if (preparedInput.m_fUseVerbose == 1) {
3023 Info("DiTauMassTools", "Substracting pt1 from sumEt");
3024 }
3025 preparedInput.m_SumEt -= preparedInput.m_vistau1.Pt();
3026 }
3027 if (preparedInput.m_SumEt > preparedInput.m_vistau2.Pt() && preparedInput.m_vistau2.M() < 0.05 &&
3028 m_mmcCalibrationSet != MMCCalibrationSet::MMC2015HIGHMASS) {
3029 if (preparedInput.m_fUseVerbose == 1) {
3030 Info("DiTauMassTools", "Substracting pt2 from sumEt");
3031 }
3032 preparedInput.m_SumEt -= preparedInput.m_vistau2.Pt();
3033 }
3034 }
3035 }
3036
3037 // controling TauProbability settings for UPGRADE studies
3038 if (m_mmcCalibrationSet == MMCCalibrationSet::UPGRADE &&
3039 preparedInput.m_fUseDefaults == 1) {
3040 if ((preparedInput.m_vistau1.M() < 0.12 && preparedInput.m_vistau2.M() > 0.12) ||
3041 (preparedInput.m_vistau2.M() < 0.12 && preparedInput.m_vistau1.M() > 0.12)) {
3042 Prob->SetUseTauProbability(true); // lep-had case
3043 }
3044 if (preparedInput.m_vistau1.M() > 0.12 && preparedInput.m_vistau2.M() > 0.12) {
3045 Prob->SetUseTauProbability(false); // had-had case
3046 }
3047 }
3048
3049 // change Beam Energy for different running conditions
3050 preparedInput.m_beamEnergy = m_beamEnergy;
3051
3052 //--------------------- pre-set defaults for Run-2. To disable pre-set
3053 // defaults set fUseDefaults=0
3054 if (preparedInput.m_fUseDefaults == 1) {
3055 if (m_mmcCalibrationSet == MMCCalibrationSet::MMC2015HIGHMASS) {
3056 SetNsigmaMETscan_ll(4.0);
3057 SetNsigmaMETscan_lh(4.0);
3058 SetNsigmaMETscan_hh(4.0);
3059 preparedInput.m_fUseTailCleanup = 0;
3060 if ((preparedInput.m_vistau1.M() < 0.12 && preparedInput.m_vistau2.M() > 0.12) ||
3061 (preparedInput.m_vistau2.M() < 0.12 && preparedInput.m_vistau1.M() > 0.12))
3062 Prob->SetUseTauProbability(false); // lep-had
3063 if (preparedInput.m_tauTypes == TauTypes::hh)
3064 Prob->SetUseTauProbability(true); // had-had
3065 Prob->SetUseMnuProbability(false);
3066 }
3067 }
3068
3069 // compute HTOffset if relevant
3070 if (Prob->GetUseHT()) // use missing Ht for Njet25=0 events
3071 {
3072 // dPhi(l-t) dependence of misHt-trueMET
3073 double HtOffset = 0.;
3074 // proper for hh
3075 if (preparedInput.m_tauTypes == TauTypes::hh) {
3076 // hh
3077 double x = preparedInput.m_DelPhiTT;
3078 HtOffset = 87.5 - 27.0 * x;
3079 }
3080
3081 preparedInput.m_HtOffset = HtOffset;
3082
3083 // if use HT, replace MET with HT
3084 preparedInput.m_METsigmaP =
3085 preparedInput.m_MHtSigma2; // sigma of 2nd Gaussian for missing Ht resolution
3086 preparedInput.m_METsigmaL = preparedInput.m_MHtSigma2;
3087
3088 PtEtaPhiMVector tauSum = preparedInput.m_vistau1 + preparedInput.m_vistau2;
3089 preparedInput.m_MetVec.SetXY(-tauSum.Px(), -tauSum.Py()); // WARNING this replace metvec by -mht
3090 }
3091}
3092
3093void MissingMassCalculator::SaveLlhHisto(const bool val){
3094 m_SaveLlhHisto=val;
3095 if(!m_SaveLlhHisto) return;
3096
3097 float hEmax = 3000.0; // maximum energy (GeV)
3098 int hNbins = 1500;
3099 m_fMEtP_all = std::make_shared<TH1F>("MEtP_h1", "M", hNbins, -100.0,
3100 100.); // all solutions
3101 m_fMEtL_all = std::make_shared<TH1F>("MEtL_h1", "M", hNbins, -100.0,
3102 100.); // all solutions
3103 m_fMnu1_all = std::make_shared<TH1F>("Mnu1_h1", "M", hNbins, 0.0,
3104 hEmax); // all solutions
3105 m_fMnu2_all = std::make_shared<TH1F>("Mnu2_h1", "M", hNbins, 0.0,
3106 hEmax); // all solutions
3107 m_fPhi1_all = std::make_shared<TH1F>("Phi1_h1", "M", hNbins, -10.0,
3108 10.); // all solutions
3109 m_fPhi2_all = std::make_shared<TH1F>("Phi2_h1", "M", hNbins, -10.0,
3110 10.); // all solutions
3111 m_fMfit_allGraph = std::make_shared<TGraph>(); // all solutions
3112
3113 m_fMEtP_all->Sumw2();
3114 m_fMEtL_all->Sumw2();
3115 m_fMnu1_all->Sumw2();
3116 m_fMnu2_all->Sumw2();
3117 m_fPhi1_all->Sumw2();
3118 m_fPhi2_all->Sumw2();
3119
3120 m_fMEtP_all->SetDirectory(0);
3121 m_fMEtL_all->SetDirectory(0);
3122 m_fMnu1_all->SetDirectory(0);
3123 m_fMnu2_all->SetDirectory(0);
3124 m_fPhi1_all->SetDirectory(0);
3125 m_fPhi2_all->SetDirectory(0);
3126}
3127
3128void MissingMassCalculator::SetUseFloatStopping(const bool val){
3129 m_fUseFloatStopping=val;
3130 if(!m_fUseFloatStopping) return;
3131
3132 float hEmax = 3000.0; // maximum energy (GeV)
3133 int hNbins = 1500;
3134 m_fMmass_split1 = std::make_shared<TH1F>("mass_h1_1", "M", hNbins, 0.0, hEmax);
3135 m_fMEtP_split1 = std::make_shared<TH1F>("MEtP_h1_1", "M", hNbins, -100.0, 100.0);
3136 m_fMEtL_split1 = std::make_shared<TH1F>("MEtL_h1_1", "M", hNbins, -100.0, 100.0);
3137 m_fMnu1_split1 = std::make_shared<TH1F>("Mnu1_h1_1", "M", hNbins, 0.0, hEmax);
3138 m_fMnu2_split1 = std::make_shared<TH1F>("Mnu2_h1_1", "M", hNbins, 0.0, hEmax);
3139 m_fPhi1_split1 = std::make_shared<TH1F>("Phi1_h1_1", "M", hNbins, -10.0, 10.0);
3140 m_fPhi2_split1 = std::make_shared<TH1F>("Phi2_h1_1", "M", hNbins, -10.0, 10.0);
3141 m_fMmass_split2 = std::make_shared<TH1F>("mass_h1_2", "M", hNbins, 0.0, hEmax);
3142 m_fMEtP_split2 = std::make_shared<TH1F>("MEtP_h1_2", "M", hNbins, -100.0, 100.0);
3143 m_fMEtL_split2 = std::make_shared<TH1F>("MEtL_h1_2", "M", hNbins, -100.0, 100.0);
3144 m_fMnu1_split2 = std::make_shared<TH1F>("Mnu1_h1_2", "M", hNbins, 0.0, hEmax);
3145 m_fMnu2_split2 = std::make_shared<TH1F>("Mnu2_h1_2", "M", hNbins, 0.0, hEmax);
3146 m_fPhi1_split2 = std::make_shared<TH1F>("Phi1_h1_2", "M", hNbins, -10.0, 10.0);
3147 m_fPhi2_split2 = std::make_shared<TH1F>("Phi2_h1_2", "M", hNbins, -10.0, 10.0);
3148
3149 m_fMmass_split1->Sumw2();
3150 m_fMEtP_split1->Sumw2();
3151 m_fMEtL_split1->Sumw2();
3152 m_fMnu1_split1->Sumw2();
3153 m_fMnu2_split1->Sumw2();
3154 m_fPhi1_split1->Sumw2();
3155 m_fPhi2_split1->Sumw2();
3156 m_fMmass_split2->Sumw2();
3157 m_fMEtP_split2->Sumw2();
3158 m_fMEtL_split2->Sumw2();
3159 m_fMnu1_split2->Sumw2();
3160 m_fMnu2_split2->Sumw2();
3161 m_fPhi1_split2->Sumw2();
3162 m_fPhi2_split2->Sumw2();
3163
3164 m_fMmass_split1->SetDirectory(0);
3165 m_fMEtP_split1->SetDirectory(0);
3166 m_fMEtL_split1->SetDirectory(0);
3167 m_fMnu1_split1->SetDirectory(0);
3168 m_fMnu2_split1->SetDirectory(0);
3169 m_fPhi1_split1->SetDirectory(0);
3170 m_fPhi2_split1->SetDirectory(0);
3171 m_fMmass_split2->SetDirectory(0);
3172 m_fMEtP_split2->SetDirectory(0);
3173 m_fMEtL_split2->SetDirectory(0);
3174 m_fMnu1_split2->SetDirectory(0);
3175 m_fMnu2_split2->SetDirectory(0);
3176 m_fPhi1_split2->SetDirectory(0);
3177 m_fPhi2_split2->SetDirectory(0);
3178}
3179
3180// Add CollinearMass calculation
3181bool MissingMassCalculator::MassCollinear(const xAOD::IParticle *p0, const xAOD::IParticle *p1,
3182 const xAOD::MissingET *met, // met
3183 const bool kMMCsynchronize, // mmc sychronization
3184 double &mass, double &xp1, double &xp2) { // result
3185
3186 TLorentzVector k1 = p0->p4();
3187 TLorentzVector k2 = p1->p4();
3188
3190 if (kMMCsynchronize) {
3191 if (p0->type() == xAOD::Type::Tau) {
3192 const xAOD::TauJet *tau0 = static_cast<const xAOD::TauJet *>(p0);
3193 k1.SetPtEtaPhiM(k1.Pt(), k1.Eta(), k1.Phi(),
3194 tau0->nTracks() < 3 ? 800. : 1200.); // MeV
3195 }
3196
3197 if (p1->type() == xAOD::Type::Tau) {
3198 const xAOD::TauJet *tau1 = static_cast<const xAOD::TauJet *>(p1);
3199 k2.SetPtEtaPhiM(k2.Pt(), k2.Eta(), k2.Phi(),
3200 tau1->nTracks() < 3 ? 800. : 1200.); // MeV
3201 }
3202 }
3203
3204 TMatrixD K(2, 2);
3205 K(0, 0) = k1.Px();
3206 K(0, 1) = k2.Px();
3207 K(1, 0) = k1.Py();
3208 K(1, 1) = k2.Py();
3209
3210 if (K.Determinant() == 0)
3211 return false;
3212
3213 TMatrixD M(2, 1);
3214 M(0, 0) = met->mpx();
3215 M(1, 0) = met->mpy();
3216
3217 TMatrixD Kinv = K.Invert();
3218
3219 TMatrixD X(2, 1);
3220 X = Kinv * M;
3221
3222 double X1 = X(0, 0);
3223 double X2 = X(1, 0);
3224 double x1 = 1. / (1. + X1);
3225 double x2 = 1. / (1. + X2);
3226
3227 TLorentzVector par1 = k1 * (1 / x1);
3228 TLorentzVector par2 = k2 * (1 / x2);
3229
3230 double m = (par1 + par2).M();
3231
3232 // return to caller
3233 mass = m;
3234
3235 if (k1.Pt() > k2.Pt()) {
3236 xp1 = x1;
3237 xp2 = x2;
3238 } else {
3239 xp1 = x2;
3240 xp2 = x1;
3241 }
3242
3243 return true;
3244}
3245
3246
3247
Phi_mpi_pi
__HOSTDEV__ double Phi_mpi_pi(double)
Definition GeoRegion.cxx:10
a
static Double_t a
Definition LArPhysWaveHECTool.cxx:38
P
static Double_t P(Double_t *tt, Double_t *par)
Definition LArPhysWaveHECTool.cxx:280
tau0
static Double_t tau0
Definition LArPhysWaveHECTool.cxx:37
debug
const bool debug
Definition MakeUncertaintyPlots.cxx:53
ParticleConstants.h
A number of constexpr particle constants to avoid hardcoding them directly in various places.
GEV
#define GEV
Definition PrintPhotonSF.cxx:25
x
#define x
DiTauMassTools::MissingMassCalculator::TailCleanUp
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)
Definition MissingMassCalculator.cxx:1991
DiTauMassTools::MissingMassCalculator::maxFromHist
double maxFromHist(TH1F *theHist, std::vector< double > &histInfo, const MaxHistStrategy::e maxHistStrategy=MaxHistStrategy::FIT, const int winHalfWidth=2, bool debug=false)
Definition MissingMassCalculator.cxx:1514
DiTauMassTools::MissingMassCalculator::MassCollinear
bool MassCollinear(const xAOD::IParticle *p0, const xAOD::IParticle *p1, const xAOD::MissingET *met, const bool kMMCsynchronize, double &mass, double &xp1, double &xp2)
Definition MissingMassCalculator.cxx:3181
DiTauMassTools::MissingMassCalculator::probCalculatorV9fast
int probCalculatorV9fast(const double &phi1, const double &phi2, const double &M_nu1, const double &M_nu2)
Definition MissingMassCalculator.cxx:1793
DiTauMassTools::MissingMassCalculator::FinalizeSettings
void FinalizeSettings(const xAOD::IParticle *part1, const xAOD::IParticle *part2, const xAOD::MissingET *met, const int &njets)
Definition MissingMassCalculator.cxx:2824
DiTauMassTools::MissingMassCalculator::m_nu2FinalSolOldVec
std::vector< PtEtaPhiMVector > m_nu2FinalSolOldVec
Definition MissingMassCalculator.h:155
DiTauMassTools::MissingMassCalculator::refineSolutions
int refineSolutions(const double &M_nu1, const double &M_nu2, const int nsol1, const int nsol2, const double &Mvis, const double &Meff)
Definition MissingMassCalculator.cxx:1819
DiTauMassTools::MissingMassCalculator::MissingMassCalculator
MissingMassCalculator(MMCCalibrationSet::e aset, std::string paramFilePath)
Definition MissingMassCalculator.cxx:52
DiTauMassTools::MissingMassCalculator::m_nu1FinalSolOldVec
std::vector< PtEtaPhiMVector > m_nu1FinalSolOldVec
Definition MissingMassCalculator.h:154
DiTauMassTools::MissingMassCalculator::m_tauvecsol1
std::vector< PtEtaPhiMVector > m_tauvecsol1
Definition MissingMassCalculator.h:79
DiTauMassTools::MissingMassCalculator::m_nuvecsol1
std::vector< PtEtaPhiMVector > m_nuvecsol1
Definition MissingMassCalculator.h:76
DiTauMassTools::MissingMassCalculator::m_nuvecsol2
std::vector< PtEtaPhiMVector > m_nuvecsol2
Definition MissingMassCalculator.h:77
DiTauMassTools::MissingMassCalculator::m_nu2FinalSolVec
std::vector< PtEtaPhiMVector > m_nu2FinalSolVec
Definition MissingMassCalculator.h:161
DiTauMassTools::MissingMassCalculator::maxFitting
Double_t maxFitting(Double_t *x, Double_t *par)
Definition MissingMassCalculator.cxx:1494
DiTauMassTools::MissingMassCalculator::NuPsolutionV3
int NuPsolutionV3(const double &mNu1, const double &mNu2, const double &phi1, const double &phi2, int &nsol1, int &nsol2)
Definition MissingMassCalculator.cxx:547
DiTauMassTools::MissingMassCalculator::dTheta3DLimit
double dTheta3DLimit(const int &tau_type, const int &limit_code, const double &P_tau)
Definition MissingMassCalculator.cxx:2719
DiTauMassTools::MissingMassCalculator::NuPsolutionLFV
int NuPsolutionLFV(const XYVector &met_vec, const PtEtaPhiMVector &tau, const double &m_nu, std::vector< PtEtaPhiMVector > &nu_vec)
Definition MissingMassCalculator.cxx:752
DiTauMassTools::MissingMassCalculator::m_tauvecsol2
std::vector< PtEtaPhiMVector > m_tauvecsol2
Definition MissingMassCalculator.h:80
DiTauMassTools::MissingMassCalculator::m_nu1FinalSolVec
std::vector< PtEtaPhiMVector > m_nu1FinalSolVec
Definition MissingMassCalculator.h:160
DiTauMassTools::MissingMassCalculator::RunMissingMassCalculator
int RunMissingMassCalculator(const xAOD::IParticle *part1, const xAOD::IParticle *part2, const xAOD::MissingET *met, const int &njets)
Definition MissingMassCalculator.cxx:195
xAOD::IParticle
Class providing the definition of the 4-vector interface.
Definition Event/xAOD/xAODBase/xAODBase/IParticle.h:41
xAOD::TauJet_v3::nTracks
size_t nTracks(TauJetParameters::TauTrackFlag flag=TauJetParameters::TauTrackFlag::classifiedCharged) const
Definition TauJet_v3.cxx:488
binWidth
void binWidth(TH1 *h)
Definition listroot.cxx:80
DiTauMassTools::getLFVMode
int getLFVMode(const xAOD::IParticle *p1, const xAOD::IParticle *p2, int mmcType1, int mmcType2)
Definition PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx:92
DiTauMassTools::MaxDelPhi
double MaxDelPhi(int tau_type, double Pvis, double dRmax_tau)
Definition PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx:24
DiTauMassTools::Angle
double Angle(const VectorType1 &vec1, const VectorType2 &vec2)
Definition PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:60
DiTauMassTools::fastSinCos
void fastSinCos(const double &phi, double &sinPhi, double &cosPhi)
Definition PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx:51
ParticleConstants::PDG2011::tauMassInMeV
constexpr double tauMassInMeV
the mass of the tau (in MeV)
Definition ParticleConstants.h:32
part1
Definition part1.py:1
part2
Definition part2.py:1
std::swap
void swap(ElementLinkVector< DOBJ > &lhs, ElementLinkVector< DOBJ > &rhs)
Definition AthLinks/ElementLinkVector.h:485
xAODType::Tau
@ Tau
The object is a tau (jet).
Definition ObjectType.h:49
xAOD::MissingET
MissingET_v1 MissingET
Version control by type defintion.
Definition Event/xAOD/xAODMissingET/xAODMissingET/MissingET.h:15
xAOD::TauJet
TauJet_v3 TauJet
Definition of the current "tau version".
Definition TauJet.h:17
unlikely
#define unlikely(x)
Definition pythonic_coracool.cxx:9
Generated on for ATLAS Offline Software by doxygen 1.17.0