ATLAS Offline Software
PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 // vim: ts=2 sw=2
6 
7 // if included, do not use fast calculation of sin/cos but built-in functions
8 //#define NOFASTSINCOS
9 
10 // local include(s)
12 
13 using namespace DiTauMassTools;
14 
15 double DiTauMassTools::MaxDelPhi(int tau_type, double Pvis, double dRmax_tau)
16 {
17  return dRmax_tau+0*Pvis*tau_type; // hack to avoid warning
18 }
19 
20 //put back phi within -pi, +pi
21 double DiTauMassTools::fixPhiRange (const double & phi)
22 {
23  double phiOut=phi;
24 
25  if (phiOut>0){
26  while (phiOut>TMath::Pi()) {
27  phiOut-=TMath::TwoPi();
28  }
29  }
30 
31  else{
32  while (phiOut<-TMath::Pi()) {
33  phiOut+=TMath::TwoPi();
34  }
35  }
36  return phiOut;
37 }
38 
39 // fast approximate calculation of sin and cos
40 // approximation good to 1 per mill. c^2+s^2=1 strictly exact though
41 // it is like using slightly different values of phi1 and phi2
42 void DiTauMassTools::fastSinCos (const double & phiInput, double & sinPhi, double & cosPhi)
43 {
44  const double fastB=4/TMath::Pi();
45  const double fastC=-4/(TMath::Pi()*TMath::Pi());
46  const double fastP=9./40.;
47  const double fastQ=31./40.;
48  // use normal sin cos if switch off
49 #ifdef NOFASTSINCOS
50  sinPhi=sin(phiInput);
51  cosPhi=cos(phiInput);
52 #else
53 
54 
55  double phi = fixPhiRange(phiInput);
56 
57  // http://devmaster.net/forums/topic/4648-fast-and-accurate-sinecosine/
58  // accurate to 1 per mille
59 
60  // even faster
61  // const double y=fastB*phi+fastC*phi*std::abs(phi);
62  // sinPhi=fastP*(y*std::abs(y)-y)+y;
63  const double y=phi*(fastB+fastC*std::abs(phi));
64  sinPhi=y*(fastP*std::abs(y)+fastQ);
65 
66 
67  //note that one could use cos(phi)=sin(phi+pi/2), however then one would not have c^2+s^2=1 (would get it only within 1 per mille)
68  // the choice here is to keep c^2+s^2=1 so everything is as one would compute c and s from a slightly (1 per mille) different angle
69  cosPhi=sqrt(1-std::pow(sinPhi,2));
70  if (std::abs(phi)>TMath::PiOver2()) cosPhi=-cosPhi;
71 #endif
72 }
73 
74 // return true if cache was uptodate
75 bool DiTauMassTools::updateDouble (const double in, double & out)
76 {
77  if (out==in) return true;
78  out=in;
79  return false;
80 }
81 
82 //CP::CorrectionCode MissingMassTool::getLFVMode( const xAOD::IParticle* p1, const xAOD::IParticle* p2, int mmcType1, int mmcType2) {
83 int DiTauMassTools::getLFVMode( const xAOD::IParticle* p1, const xAOD::IParticle* p2, int mmcType1, int mmcType2) {
84 
85  // Check if particles pointers are null
86  if(p1 == nullptr || p2 == nullptr) {
87  Info("DiTauMassTools", "MissingMassTool::getLFVMode() got a nullptr - returning -1");
88  return -1;
89  }
90 
91  // In case we're in LFV calibration, pass the LFV mode to the tool
92  if (mmcType1 == -1) mmcType1 = mmcType(p1);
93  if (mmcType1 < 0) return -1;//return CP::CorrectionCode::Error;
94 
95  if (mmcType2 == -1) mmcType2 = mmcType(p2);
96  if (mmcType2 < 0) return -1;//return CP::CorrectionCode::Error;
97 
98  // we don't use mmcType as it's 0 for both leptons
99  const xAOD::IParticle *p;
100  if(mmcType1 == 8 && mmcType2 == 8) {
101  // both leptonic; find whichever has the highest pT
102  if(p1->pt() > p2->pt() ) {
103  p = p1;
104  } else {
105  p = p2;
106  }
107 
108  } else {
109  // one of them is a lepton, find it
110  if(mmcType1 == 8) {
111  p = p1;
112  } else if(mmcType2 == 8) {
113  p = p2;
114  } else {
115  // if you're here, you've passed 2 taus to the LFV mode.
116  Warning("DiTauMassTools", "Trying to set LFV mode for 2 taus!");
117  return -1;//return CP::CorrectionCode::Error;
118  }
119  }
120 
121  if(!p) return -1;//return CP::CorrectionCode::Error;
122 
123  int LFVMode = -1;
124  if(p->type() == xAOD::Type::Muon) {
125  // mu+tau mode
126  LFVMode = 1;
127  } else {
128  // e+tau mode
129  LFVMode = 0;
130  }
131 
132  return LFVMode;
133  //return CP::CorrectionCode::Ok;
134 }
135 
137 {
138  if(part == nullptr) {
139  // idea for fix: pass logger object to this function or whole algorithm that inherited from athena as pointer
140  Info("DiTauMassTools", "MissingMassTool::mmcType() got a nullptr - returning -1");
141  return -1;
142  }
143 
144  xAOD::Type::ObjectType partType=part->type();
145  int aType=-1;
146  if (partType==xAOD::Type::Electron || partType==xAOD::Type::Muon)
147  {
148  aType=8; // PanTau_DecayMode for leptonic taus
149  }
150  else if (partType==xAOD::Type::Tau)
151  {
152  const xAOD::TauJet * aTauJet=dynamic_cast<const xAOD::TauJet*>(part);
153  if (aTauJet==0)
154  {
155  Warning("DiTauMassTools", "MissingMassTool::mmcType() dynamic_cast of TauJet failed");
156  aType=-2;
157  }
158  else
159  {
160  aTauJet->panTauDetail(xAOD::TauJetParameters::PanTau_DecayMode, aType); // 0: 1p0n, 1: 1p1n, 2: 1pXn, 3: 3p0n, 4: 3pXn, 5: Other (2p, 4p, 5p), 6: not set (0p, >=6p), 7: error, 8: leptonic
161  }
162  }
163  else
164  {
165  Warning("DiTauMassTools", "MissingMassTool::mmcType() unrecognised particle type! Only Electron, Muon, TauJet allowed. If no mistake, please call MissingMassTool::calculate() directly.");
166  aType=-1;
167  }
168 
169  return aType;
170 }
171 //________________________________________________________________________
172 
173 void DiTauMassTools::readInParams(TDirectory* dir, MMCCalibrationSet::e aset, std::vector<TF1*>& lep_numass, std::vector<TF1*>& lep_angle, std::vector<TF1*>& lep_ratio, std::vector<TF1*>& had_angle, std::vector<TF1*>& had_ratio) {
174  std::string paramcode;
175  if (aset == MMCCalibrationSet::MMC2019) paramcode = "MMC2019MC16";
176  else if (aset == MMCCalibrationSet::MMC2024) paramcode = "MMC2024MC23";
177  else {
178  Info("DiTauMassTools", "The specified calibration version does not support root file parametrisations");
179  return;
180  }
181  TIter next(dir->GetListOfKeys());
182  TKey* key;
183  while ((key = (TKey*)next())) {
184  TClass *cl = gROOT->GetClass(key->GetClassName());
185  // if there is another subdirectory, go into that dir
186  if (cl->InheritsFrom("TDirectory")) {
187  dir->cd(key->GetName());
188  TDirectory *subdir = gDirectory;
189  readInParams(subdir, aset, lep_numass, lep_angle, lep_ratio, had_angle, had_ratio);
190  dir->cd();
191  }
192  else if (cl->InheritsFrom("TF1") || cl->InheritsFrom("TGraph")) {
193  // get parametrisations and sort them into their corresponding vectors
194  std::string total_path = dir->GetPath();
195  if (total_path.find(paramcode) == std::string::npos) continue;
196  TF1* func = (TF1*)dir->Get( (const char*) key->GetName() );
197  TF1* f = new TF1(*func);
198  if (total_path.find("lep") != std::string::npos)
199  {
200  if (total_path.find("Angle") != std::string::npos){
201  lep_angle.push_back(f);
202  }
203  else if (total_path.find("Ratio") != std::string::npos)
204  {
205  lep_ratio.push_back(f);
206  }
207  else if (total_path.find("Mass") != std::string::npos)
208  {
209  lep_numass.push_back(f);
210  }
211  else
212  {
213  Warning("DiTauMassTools", "Undefined leptonic PDF term in input file.");
214  }
215  }
216  else if (total_path.find("had") != std::string::npos)
217  {
218  if (total_path.find("Angle") != std::string::npos){
219  had_angle.push_back(f);
220  }
221  else if (total_path.find("Ratio") != std::string::npos)
222  {
223  had_ratio.push_back(f);
224  }
225  else
226  {
227  Warning("DiTauMassTools", "Undefined hadronic PDF term in input file.");
228  }
229  }
230  else
231  {
232  Warning("DiTauMassTools", "Undefined decay channel in input file.");
233  }
234  }
235  else {
236  Warning("DiTauMassTools", "Class in input file not recognized.");
237  }
238  }
239 }
DiTauMassTools::MMCCalibrationSet::MMC2024
@ MMC2024
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:40
LArG4FSStartPointFilter.part
part
Definition: LArG4FSStartPointFilter.py:21
xAOD::TauJetParameters::PanTau_DecayMode
@ PanTau_DecayMode
Definition: TauDefs.h:360
xAOD::Electron
Electron_v1 Electron
Definition of the current "egamma version".
Definition: Event/xAOD/xAODEgamma/xAODEgamma/Electron.h:17
ObjectType
ObjectType
Definition: BaseObject.h:11
TRTCalib_cfilter.p1
p1
Definition: TRTCalib_cfilter.py:130
python.AthDsoLogger.out
out
Definition: AthDsoLogger.py:71
HelperFunctions.h
xAOD::TauJet_v3::panTauDetail
bool panTauDetail(TauJetParameters::PanTauDetails panTauDetail, int &value) const
Get and set values of pantau details variables via enum.
Definition: TauJet_v3.cxx:367
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
xAOD::IParticle
Class providing the definition of the 4-vector interface.
Definition: Event/xAOD/xAODBase/xAODBase/IParticle.h:41
DiTauMassTools::fastSinCos
void fastSinCos(const double &phi, double &sinPhi, double &cosPhi)
Definition: PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx:42
DiTauMassTools::updateDouble
bool updateDouble(const double in, double &out)
Definition: PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx:75
DiTauMassTools
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:24
TRTCalib_cfilter.p2
p2
Definition: TRTCalib_cfilter.py:131
DiTauMassTools::MaxDelPhi
double MaxDelPhi(int tau_type, double Pvis, double dRmax_tau)
Definition: PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx:15
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
fillPileUpNoiseLumi.next
next
Definition: fillPileUpNoiseLumi.py:52
ParseInputs.gDirectory
gDirectory
Definition: Final2012/ParseInputs.py:133
DiTauMassTools::readInParams
void readInParams(TDirectory *dir, MMCCalibrationSet::e aset, std::vector< TF1 * > &lep_numass, std::vector< TF1 * > &lep_angle, std::vector< TF1 * > &lep_ratio, std::vector< TF1 * > &had_angle, std::vector< TF1 * > &had_ratio)
Definition: PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx:173
xAOD::TauJet_v3
Class describing a tau jet.
Definition: TauJet_v3.h:41
DiTauMassTools::getLFVMode
int getLFVMode(const xAOD::IParticle *p1, const xAOD::IParticle *p2, int mmcType1, int mmcType2)
Definition: PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx:83
DiTauMassTools::MMCCalibrationSet::MMC2019
@ MMC2019
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:40
hist_file_dump.f
f
Definition: hist_file_dump.py:141
DiTauMassTools::MMCCalibrationSet::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:40
DiTauMassTools::mmcType
int mmcType(const xAOD::IParticle *part)
Definition: PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx:136
beamspotman.dir
string dir
Definition: beamspotman.py:623
DiTauMassTools::fixPhiRange
double fixPhiRange(const double &phi)
Definition: PhysicsAnalysis/TauID/DiTauMassTools/Root/HelperFunctions.cxx:21
Muon
struct TBPatternUnitContext Muon
y
#define y
xAODType::Tau
@ Tau
The object is a tau (jet)
Definition: ObjectType.h:49
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
python.LumiCalcRecover.subdir
subdir
Definition: LumiCalcRecover.py:25
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
dq_make_web_display.cl
cl
print [x.__class__ for x in toList(dqregion.getSubRegions()) ]
Definition: dq_make_web_display.py:26
mapkey::key
key
Definition: TElectronEfficiencyCorrectionTool.cxx:37