5
5
6
6
#include < iostream>
7
7
#include < string>
8
- #include < assert.h>
8
+ #include < assert.h>
9
9
#include < math.h>
10
10
#include < cmath>
11
11
@@ -30,7 +30,7 @@ HadTauTFCrystalBall::HadTauTFCrystalBall(const std::string& inputFilePath)
30
30
pp_(0 ),
31
31
genPt_cache_(-1 .),
32
32
genEta_cache_(0 .)
33
- {
33
+ {
34
34
// std::cout << "<HadTauTFCrystalBall::HadTauTFCrystalBall>:" << std::endl;
35
35
36
36
xx_ = new double [1 ];
@@ -45,7 +45,7 @@ HadTauTFCrystalBall::HadTauTFCrystalBall(const std::string& inputFilePath)
45
45
parDir[5 ] = inputFilePath + " /par5/" ;
46
46
parDir[6 ] = inputFilePath + " /par6/" ;
47
47
parDir[7 ] = inputFilePath + " /par7/" ;
48
-
48
+
49
49
std::map<int , std::string> fileL2;
50
50
fileL2[kAll ] = " TauJec11V1_L2Relative_AK5tauHPSlooseCombDBcorrAll.txt" ;
51
51
fileL2[kOneProng0Pi0 ] = " TauJec11V1_L2Relative_AK5tauHPSlooseCombDBcorrOneProng0Pi0.txt" ;
@@ -85,19 +85,19 @@ HadTauTFCrystalBall::~HadTauTFCrystalBall()
85
85
delete mapPar_[decayMode][par];
86
86
}
87
87
}
88
- }
88
+ }
89
89
90
90
void HadTauTFCrystalBall::setDecayMode (int decayMode) const
91
- {
92
- decayMode_ = decayMode;
91
+ {
92
+ decayMode_ = decayMode;
93
93
94
- int idxDecayMode = -1 ;
94
+ int idxDecayMode = -1 ;
95
95
if ( decayMode_ == reco::PFTau::kOneProng0PiZero ) idxDecayMode = kOneProng0Pi0 ;
96
96
else if ( decayMode_ == reco::PFTau::kOneProng1PiZero ) idxDecayMode = kOneProng1Pi0 ;
97
97
else if ( decayMode_ == reco::PFTau::kOneProng2PiZero ) idxDecayMode = kOneProng2Pi0 ;
98
98
else if ( decayMode_ == reco::PFTau::kThreeProng0PiZero ) idxDecayMode = kThreeProng0Pi0 ;
99
99
else {
100
- std::cerr << " Warning: No transfer function defined for decay mode = " << decayMode_ << " !!" ;
100
+ std::cerr << " Warning: No transfer function defined for decay mode = " << decayMode_ << " !!" ;
101
101
idxDecayMode = kAll ;
102
102
}
103
103
assert (mapPar_.find (idxDecayMode) != mapPar_.end ());
@@ -110,8 +110,8 @@ void HadTauTFCrystalBall::setDecayMode(int decayMode) const
110
110
111
111
namespace
112
112
{
113
- double fnc_dscb (double * xx, double * pp)
114
- {
113
+ double fnc_dscb (double * xx, double * pp)
114
+ {
115
115
// std::cout << "<fnc_dscb>:" << std::endl;
116
116
double x = xx[0 ];
117
117
// gaussian core
@@ -154,8 +154,8 @@ namespace
154
154
// std::cout << "A2 = " << A2 << ", a2 + B2 = " << a2_plus_B2 << ", 2. - mu + B2*sig = " << (2. - mu + B2*sig) << ", -1 + p2 = " << (-1. + p2) << std::endl;
155
155
term3 = A2*TMath::Power (a2_plus_B2, -p2)*(a2_plus_B2*sig - TMath::Power (a2_plus_B2*sig, p2)*TMath::Power (TMath::Min (2 ., 2 . - mu + B2*sig), 1 . - p2))/((-1 . + p2)*sig);
156
156
if ( !std::isfinite (term3) ) term3 = 0 .;
157
- }
158
- // std::cout << "term1 = " << term1 << ", term2 = " << term2 << ", term3 = " << term3 << ", sig = " << sig << std::endl;
157
+ }
158
+ // std::cout << "term1 = " << term1 << ", term2 = " << term2 << ", term3 = " << term3 << ", sig = " << sig << std::endl;
159
159
double one_over_N = term1 + term2 + term3;
160
160
one_over_N *= sig; // CV: multiply result obtained from Mathematica by Jacobi factor dx/du = sig for variable transformation from x to u
161
161
// std::cout << "1/N = " << one_over_N << std::endl;
@@ -167,7 +167,7 @@ namespace
167
167
result = N*A1*TMath::Power (B1 - u, -p1);
168
168
} else if ( u < a2 ) {
169
169
double arg = -0.5 *u*u;
170
- // std::cout << "N = " << N << ", u = " << u << ", arg = " << arg << std::endl;
170
+ // std::cout << "N = " << N << ", u = " << u << ", arg = " << arg << std::endl;
171
171
result = N*TMath::Exp (arg);
172
172
} else {
173
173
// std::cout << "N = " << N << ", A2 = " << A2 << ", B2 + u = " << (B2 + u) << std::endl;
@@ -179,26 +179,26 @@ namespace
179
179
}
180
180
181
181
double HadTauTFCrystalBall::operator ()(double recPt, double genPt, double genEta) const
182
- {
182
+ {
183
183
// std::cout << "<HadTauTFCrystalBall::operator()>:" << std::endl;
184
184
// std::cout << " pT: rec = " << recPt << ", gen = " << genPt << std::endl;
185
185
// std::cout << " eta = " << genEta << std::endl;
186
186
187
- xx_[0 ] = recPt/genPt;
187
+ xx_[0 ] = recPt/genPt;
188
188
// std::cout << " xx = " << xx_[0] << std::endl;
189
189
190
190
if ( genPt != genPt_cache_ || genEta != genEta_cache_ ) {
191
191
for ( int iPar = 1 ; iPar < cbParSize_; ++iPar ) {
192
192
pp_[iPar] = (*thePar_[iPar])(genPt, genEta);
193
- pp_[iPar] = 1 ./(*thePar_[iPar])(genPt, genEta); // CV: temporary fix
193
+ pp_[iPar] = 1 ./(*thePar_[iPar])(genPt, genEta); // CV: temporary fix
194
194
// std::cout << " pp(" << iPar << ") = " << pp_[iPar] << std::endl;
195
- }
195
+ }
196
196
genPt_cache_ = genPt;
197
197
genEta_cache_ = genEta;
198
198
}
199
199
200
200
double retVal = fnc_dscb (xx_, pp_);
201
- // std::cout << "--> returning retVal = " << retVal << std::endl;
201
+ // std::cout << "--> returning retVal = " << retVal << std::endl;
202
202
return retVal;
203
203
}
204
204
@@ -207,23 +207,23 @@ double HadTauTFCrystalBall::integral(double recPt_low, double recPt_up, double g
207
207
// std::cout << "<HadTauTFCrystalBall::integral()>:" << std::endl;
208
208
// std::cout << " pT: rec = " << recPt_low << ".." << recPt_up << ", gen = " << genPt << std::endl;
209
209
// std::cout << " eta = " << genEta << std::endl;
210
-
210
+
211
211
if ( !(recPt_low < recPt_up) ) return 0 .;
212
212
213
- double x_low = recPt_low/genPt;
214
- double x_up = recPt_up/genPt;
213
+ double x_low = recPt_low/genPt;
214
+ double x_up = recPt_up/genPt;
215
215
// std::cout << " x = " << x_low << ".." << x_up << std::endl;
216
-
216
+
217
217
if ( genPt != genPt_cache_ || genEta != genEta_cache_ ) {
218
218
for ( int iPar = 1 ; iPar < cbParSize_; ++iPar ) {
219
219
pp_[iPar] = (*thePar_[iPar])(genPt, genEta);
220
- pp_[iPar] = 1 ./(*thePar_[iPar])(genPt, genEta); // CV: temporary fix
220
+ pp_[iPar] = 1 ./(*thePar_[iPar])(genPt, genEta); // CV: temporary fix
221
221
// std::cout << " pp(" << iPar << ") = " << pp_[iPar] << std::endl;
222
- }
222
+ }
223
223
genPt_cache_ = genPt;
224
224
genEta_cache_ = genEta;
225
225
}
226
-
226
+
227
227
// -----------------------------------------------------------------------------
228
228
// Warning: this code needs to be identical to code in fnc_dscb function !!
229
229
double mu = pp_[1 ]; // mean
@@ -281,7 +281,7 @@ double HadTauTFCrystalBall::integral(double recPt_low, double recPt_up, double g
281
281
double term3 = A2*TMath::Power (B2_plus_u5*B2_plus_u6, -p2)*(B2_plus_u5*TMath::Power (B2_plus_u6, p2) - TMath::Power (B2_plus_u5, p2)*B2_plus_u6)/(-one_minus_p2);
282
282
if ( !std::isfinite (term3) ) term3 = 0 .;
283
283
integral += term3;
284
- }
284
+ }
285
285
return integral;
286
286
}
287
287
@@ -298,10 +298,10 @@ HadTauTFCrystalBall* HadTauTFCrystalBall::Clone(const std::string& label) const
298
298
clone->pp_ [iPar] = pp_[iPar];
299
299
}
300
300
for ( std::map<int , vHadTauTFCrystalBallParPtr>::const_iterator entryPar = mapPar_.begin ();
301
- entryPar != mapPar_.end (); ++entryPar ) {
301
+ entryPar != mapPar_.end (); ++entryPar ) {
302
302
const vHadTauTFCrystalBallParPtr& cbPars = entryPar->second ;
303
303
vHadTauTFCrystalBallParPtr cbPars_cloned;
304
- for ( int iPar = 1 ; iPar < cbParSize_; ++iPar ) {
304
+ for ( int iPar = 1 ; iPar < cbParSize_; ++iPar ) {
305
305
const HadTauTFCrystalBallPar* cbPar = cbPars[iPar];
306
306
cbPars_cloned.push_back (new HadTauTFCrystalBallPar (*cbPar));
307
307
}
0 commit comments