ATLAS Offline Software
Loading...
Searching...
No Matches
muCombUtil.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
3*/
4
5// *********************************************************************
6//
7// NAME: muCombUtil.cxx
8// PACKAGE: Trigger/TrigAlgorithms/TrigmuComb
9//
10// AUTHOR: S. Giagu <stefano.giagu@cern.ch>
11//
12// PURPOSE: Utility namespace for LVL2 Combined Muon Reco FEX Algorithm
13// *********************************************************************
14#include <iostream>
15#include <math.h>
16#include <utility>
17
18#include "muCombUtil.h"
19#include "TrigMuonEvent/MuonFeature.h"
20#include "TrigInDetEvent/TrigInDetTrack.h"
21
22namespace muCombUtil {
23
24 void setMuFastRes(std::vector<double>& vec, double p1,double p2,
25 double p3,double p4,double p5,double p6) {
26 vec.clear();
27 vec.push_back(p1);
28 vec.push_back(p2);
29 vec.push_back(p3);
30 vec.push_back(p4);
31 vec.push_back(p5);
32 vec.push_back(p6);
33 }
34
35 void setIDSCANRes(std::vector<double>& vec, double p1,double p2) {
36 vec.clear();
37 vec.push_back(p1);
38 vec.push_back(p2);
39 }
40
41
42 double getMuFastRes(const std::vector<double>& vec, const MuonFeature* feature) {
43
44 double ptGev = feature->pt();
45 double pt = ptGev*1000.; //muFast Pt(MeV)
46
47 if (pt == 0) return 1.0e33;
48
49 double AbsPtInv = fabs(1./pt); //muFast 1/Pt (1/MeV)
50 int add = feature->saddress();
51 double AbsEta = fabs(feature->eta());
52 double phi = feature->phi();
53
54 if ( add != -1) {
55 if (AbsPtInv < 0.000186) {
56 return vec[0]*AbsPtInv + vec[1]/1000.;
57 }
58 else {
59 double AbsPtInv3 = AbsPtInv*AbsPtInv*AbsPtInv;
60 double AbsPtInv2 = AbsPtInv*AbsPtInv;
61 return vec[2]*AbsPtInv3/(1000.*1000.) +
62 vec[3]*AbsPtInv2/(1000.) +
63 vec[4]*AbsPtInv +
64 vec[5]/1000.;
65 }
66 }
67 else {//Takuya/Kunihiro updated numbers
68
69 const int N_PARAMS = 5;
70 const double vparEC1[N_PARAMS] = {0.291483, -6.11348, 65.1099, -285.664, 440.041};
71 const double vparEC2[N_PARAMS] = {0.286307, -4.6759, 43.2815, -163.185, 210.786};
72 const double vparEC3[N_PARAMS] = {0.330699, -6.70755, 70.4725, -291.85, 408.739};
73 const double vparEC4[N_PARAMS] = {0.261738, -4.69971, 47.4762, -183.98, 236.813};
74 const double vparEC5[N_PARAMS] = {0.196301, -3.57276, 38.3744, -159.808, 228.256};
75 const double vparEC6[N_PARAMS] = {0.172939, -3.10788, 33.3823, -142.996, 212.957};
76 const double vparEC7[N_PARAMS] = {0.233017, -4.377, 42.5691, -171.752, 245.702};
77 const double vparEC8[N_PARAMS] = {0.22389, -4.16259, 40.1369, -162.824, 236.39};
78 const double vparEC9[N_PARAMS] = {0.197992, -3.52117, 33.5997, -136.014, 197.474};
79 const double vparECA[N_PARAMS] = {0.417289, -0.852254,-31.9257, 308.873, -719.591};
80 const double vparECB[N_PARAMS] = {0.526612, -8.04087, 82.1906, -336.87, 462.973};
81
82 double AbsPtInvGeV = AbsPtInv * 1000;
83
84 const double AbsPtInvGeVMin = 5e-3; // 200 GeV
85 const double AbsPtInvGeVMax = 0.25; // 4 GeV
86 if( AbsPtInvGeV < AbsPtInvGeVMin ) AbsPtInvGeV = AbsPtInvGeVMin;
87 if( AbsPtInvGeV > AbsPtInvGeVMax ) AbsPtInvGeV = AbsPtInvGeVMax;
88
89 //cout << "...AbsPtInvGeV=" << AbsPtInv << endl;
90
91 const double* vpar;
92 int spReg = whichECRegion(AbsEta,phi);
93 //cout << "...spReg=" << spReg << endl;
94 if ( spReg==1 ) { vpar = vparECA; }
95 else if( spReg==2 ) { vpar = vparECB; }
96 else {
97 if ( AbsEta < 1.20) { vpar = vparEC1; }
98 else if( AbsEta < 1.35) { vpar = vparEC2; }
99 else if( AbsEta < 1.50) { vpar = vparEC3; }
100 else if( AbsEta < 1.65) { vpar = vparEC4; }
101 else if( AbsEta < 1.80) { vpar = vparEC5; }
102 else if( AbsEta < 1.95) { vpar = vparEC6; }
103 else if( AbsEta < 2.10) { vpar = vparEC7; }
104 else if( AbsEta < 2.35) { vpar = vparEC8; }
105 else { vpar = vparEC9; }
106 }
107
108 double fracRes = vpar[0] + vpar[1]*AbsPtInvGeV
109 + vpar[2]*AbsPtInvGeV*AbsPtInvGeV
110 + vpar[3]*AbsPtInvGeV*AbsPtInvGeV*AbsPtInvGeV
111 + vpar[4]*AbsPtInvGeV*AbsPtInvGeV*AbsPtInvGeV*AbsPtInvGeV;
112
113 return fabs(fracRes * AbsPtInv);
114 }
115}
116
117int whichECRegion( const float eta, const float phi ) {
118 // 0: bulk
119 // 1: WeakBfield A
120 // 2: WeakBfield B
121
122 float absEta = fabs(eta);
123
124 if( ( 1.3 <= absEta && absEta < 1.45) &&
125 ( (0 <= fabs(phi) && fabs(phi) < M_PI/48. ) ||
126 (M_PI*11./48. <= fabs(phi) && fabs(phi) < M_PI*13./48. ) ||
127 (M_PI*23./48. <= fabs(phi) && fabs(phi) < M_PI*25./48. ) ||
128 (M_PI*35./48. <= fabs(phi) && fabs(phi) < M_PI*37./48. ) ||
129 (M_PI*47./48. <= fabs(phi) && fabs(phi) < M_PI )
130 )
131 ) return 1;
132
133 else if( ( 1.5 <= absEta && absEta < 1.65 ) &&
134 ( (M_PI*3./32. <= fabs(phi) && fabs(phi) < M_PI*5./32. ) ||
135 (M_PI*11./32. <= fabs(phi) && fabs(phi) < M_PI*13./32.) ||
136 (M_PI*19./32. <= fabs(phi) && fabs(phi) < M_PI*21./32.) ||
137 (M_PI*27./32. <= fabs(phi) && fabs(phi) < M_PI*29./32.)
138 )
139 ) return 2;
140
141 else return 0;
142}
143
144/*
145 double getMuFastRes(std::vector<double> barrelvec, std::vector<double> ec1vec,
146 std::vector<double> ec2vec, std::vector<double> ec3vec, std::vector<double> ec4vec,
147 const MuonFeature* feature) {
148
149 double pt = feature->pt()*1000.; //muFast 1/Pt (1/MeV)
150 if (pt == 0) return 1.0e33;
151
152 double AbsPtInv = fabs(1./pt); //muFast 1/Pt (1/MeV)
153 int add = feature->saddress();
154 double AbsEta = fabs(feature->eta());
155
156 std::vector<double> vec;
157 if ( add != -1) {
158 vec = barrelvec;
159 if (AbsPtInv < 0.000186) {
160 return vec[0]*AbsPtInv + vec[1]/1000.;
161 } else {
162 double AbsPtInv3 = AbsPtInv*AbsPtInv*AbsPtInv;
163 double AbsPtInv2 = AbsPtInv*AbsPtInv;
164 return vec[2]*AbsPtInv3/(1000.*1000.) +
165 vec[3]*AbsPtInv2/(1000.) +
166 vec[4]*AbsPtInv +
167 vec[5]/1000.;
168 }
169 } else {
170 if (AbsEta < 1.35) vec = ec1vec;
171 else if (AbsEta>=1.35 && AbsEta<1.65) vec = ec2vec;
172 else if (AbsEta>=1.65 && AbsEta<2.0) vec = ec3vec;
173 else vec = ec4vec;
174
175 if (AbsPtInv < 0.000191) {
176 return vec[0]*AbsPtInv + vec[1]/1000.;
177 } else {
178 double AbsPtInv3 = AbsPtInv*AbsPtInv*AbsPtInv;
179 double AbsPtInv2 = AbsPtInv*AbsPtInv;
180 return vec[2]*AbsPtInv3/(1000.*1000.) +
181 vec[3]*AbsPtInv2/(1000.) +
182 vec[4]*AbsPtInv +
183 vec[5]/1000.;
184 }
185 }
186 return 9999.;
187 }
188*/
189
190 double getIDSCANRes(std::vector<double> barrelvec, std::vector<double> ec1vec,
191 std::vector<double> ec2vec, std::vector<double> ec3vec, std::vector<double> ec4vec,
192 double pt_id, double eta_id) {
193
194 double pt = pt_id;
195 if (pt == 0) return 1.0e33;
196
197 double AbsPtInv = fabs(1./pt);
198 double AbsEta = fabs(eta_id);
199
200 std::vector<double> vec;
201 if (AbsEta < 1.05) vec = std::move(barrelvec);
202 else if (AbsEta>=1.05 && AbsEta<1.35) vec = std::move(ec1vec);
203 else if (AbsEta>=1.35 && AbsEta<1.65) vec = std::move(ec2vec);
204 else if (AbsEta>=1.65 && AbsEta<2.0) vec = std::move(ec3vec);
205 else vec = std::move(ec4vec);
206
207 return vec[0]*AbsPtInv+vec[1]/1000.;
208 }
209
210
211 double getMuFastEtaRes(const MuonFeature* feature) {
212
213 double pt = feature->pt();
214 double eta = feature->eta();
215 if (pt < 4. ) pt = 4.;
216 if (pt > 40.) pt = 40.;
217 bool ts = false;
218 if (feature->radius() <= 10.) ts = true;
219
220 if (fabs(eta) < 1.) {//barrel
221 return 7.75e-2/pt + 8.1e-3;
222 }
223 else {
224 if (ts) {//trigger station
225 if (fabs(eta) >= 1. && fabs(eta) < 1.2) {//ec0
226 return 2.0e-1/pt + 1.5e-3;
227 }
228 else if (fabs(eta) >= 1.2 && fabs(eta) < 1.5) {//ec1
229 return 2.5e-1/pt + 0.1e-3;
230 }
231 else if (fabs(eta) >= 1.5 && fabs(eta) < 2.0) {//ec2
232 return 1.3e-1/pt + 6.0e-3;
233 }
234 else {//ec3
235 return 1.3e-1/pt + 8.0e-3;
236 }
237 }
238 else {
239 if (fabs(eta) >= 1. && fabs(eta) < 1.6) {//ec0/1
240 return 0.071/pt + 0.0055;
241 }
242 else {//ec2/3
243 return 0.055/pt + 0.0063;
244 }
245 }
246 }
247 }
248
249
250 double getMuFastPhiRes(const MuonFeature* feature) {
251
252 double pt = feature->pt();
253 double eta = feature->eta();
254 if (pt < 4. ) pt = 4.;
255 if (pt > 40.) pt = 40.;
256 bool ts = false;
257 if (feature->radius() <= 10.) ts = true;
258
259 if (fabs(eta) < 1.) {//barrel
260 return 1.9e-1/pt - 3.4e-4;
261 }
262 else {
263 if (ts) {//trigger station
264 if (fabs(eta) >= 1. && fabs(eta) < 1.2) {//ec0
265 return 2.0e-1/pt + 2.0e-3;
266 }
267 else if (fabs(eta) >= 1.2 && fabs(eta) < 1.5) {//ec1
268 return 1.8e-1/pt + 3.0e-3;
269 }
270 else if (fabs(eta) >= 1.5 && fabs(eta) < 2.0) {//ec2
271 return 1.5e-1/pt + 8.0e-3;
272 }
273 else {//ec3
274 return 1.0e-1/pt + 1.8e-2;
275 }
276 }
277 else {
278 //if (fabs(eta) >= 1. && fabs(eta) < 1.6) {//ec0/1
279 return 0.08/pt + 0.0025;
280 //}
281 //else {//ec2/3
282 // return 0.08/pt + 0.0025;
283 //}
284 }
285 }
286 }
287
288 double getG4ExtEtaRes(double pt, double eta) {
289
290 //double pt = feature->pt();
291 //double eta = feature->eta();
292 if (pt < 4. ) pt = 4.;
293 if (pt > 40.) pt = 40.;
294 // bool ts = false;
295 // if (feature->radius() <= 10.) ts = true;
296
297 //G4 extrapolator eta resolution parametrized as A/pt^2 + B/pt + C + D*pt
298 // Resolution evaluated in regions of eta: 0.0,0.5,1.0,1.5,2.0,above...
299 if (fabs(eta) < 0.5) {
300 return 6.0e-3/(pt*pt) + 4.4e-2/pt + 1.2e-2 + (-3.2e-5*pt);
301 }
302 else if (fabs(eta) >= 0.5 && fabs(eta) < 1.0) {
303 return 5.6e-3/(pt*pt) + 3.9e-2/pt + 1.3e-2 + (-3.7e-5*pt);
304 }
305 else if (fabs(eta) >= 1.0 && fabs(eta) < 1.5) {
306 return 2.1e-1/(pt*pt) + 8.8e-2/pt + 7.1e-4 + (1.7e-4*pt);
307 }
308 else if (fabs(eta) >= 1.5 && fabs(eta) < 2.0) {
309 return 4.5e-3/(pt*pt) + 1.4e-2/pt + 1.1e-2 + (-5.7e-5*pt);
310 }
311 else { // if (fabs(eta) >= 2.0) {
312 return 3.1e-2/(pt*pt) + 4.9e-1/pt + (-2.6e-3) + (3.2e-5*pt);
313 }
314 }
315
316
317 double getG4ExtPhiRes(double pt, double eta) {
318
319 //double pt = feature->pt();
320 //double eta = feature->eta();
321 if (pt < 4. ) pt = 4.;
322 if (pt > 40.) pt = 40.;
323 // bool ts = false;
324 // if (feature->radius() <= 10.) ts = true;
325
326 //G4 extrapolator phi resolution parametrized as A/pt^2 + B/pt + C + D*pt
327 // Resolution evaluated in regions of eta: 0.0,0.5,1.0,1.5,2.0,above...
328 if (fabs(eta) < 0.5) {
329 return 4.1e-2/(pt*pt) + 1.5e-1/pt + 1.4e-3 + (5.8e-5*pt);
330 }
331 else if (fabs(eta) >= 0.5 && fabs(eta) < 1.0) {
332 return 9.4e-3/(pt*pt) + 8.8e-2/pt + 1.1e-2 + (-1.e-4*pt);
333 }
334 else if (fabs(eta) >= 1.0 && fabs(eta) < 1.5) {
335 return 3.5e-1/(pt*pt) + 9.4e-2/pt + 8.2e-3 + (-3.3e-7*pt);
336 }
337 else if (fabs(eta) >= 1.5 && fabs(eta) < 2.0) {
338 return 2.6e-2/(pt*pt) + 1.2e-1/pt + 3.8e-3 + (2.2e-5*pt);
339 }
340 else { // if (fabs(eta) >= 2.0) {
341 return 2.7e-1/(pt*pt) + 3.2e-2/pt + (9.0e-3) + (-2.6e-5*pt);
342 }
343 }
344
345 double getDeltaPhi(double phi1, double phi2) {
346 double dphi = phi1 - phi2;
347 if (dphi > M_PI) dphi -= 2*M_PI;
348 if (dphi < -M_PI) dphi += 2*M_PI;
349 return fabs(dphi);
350 }
351
352 double getDeltaEta(double eta1, double eta2) {
353 return fabs(eta1-eta2);
354 }
355
356 double getDeltaR(double eta1, double phi1, double eta2, double phi2) {
357 double deta = getDeltaEta(eta1,eta2);
358 double dphi = getDeltaPhi(phi1,phi2);
359 return sqrt(deta*deta + dphi*dphi);
360 }
361
362 double getCombinedAverage(double p1, double sp1, double p2, double sp2) {
363 if (sp1 != 0 && sp2 != 0) return (sp2*sp2*fabs(p1) + sp1*sp1*fabs(p2))/(sp1*sp1 + sp2*sp2);
364 else if (sp1 != 0 && sp2 == 0) return fabs(p2);
365 else if (sp1 == 0 && sp2 != 0) return fabs(p1);
366 else return (fabs(p1)+fabs(p2))*0.5;
367 }
368
369 double getCombinedAverageSigma(double sp1, double sp2) {
370 if (sp1 != 0 && sp2 != 0) return sqrt((sp1*sp1*sp2*sp2)/(sp1*sp1 + sp2*sp2));
371 else return 0.0;
372 }
373
374 double getChi2(int& ndof, double ipt,
375 double eta1, double seta1, double phi1, double sphi1, double ipt1, double sipt1,
376 double eta2, double seta2, double phi2, double sphi2, double ipt2, double sipt2, bool useAbsPt) {
377
378 double deta = getDeltaEta(eta1,eta2);
379 double sdeta = seta1*seta1+seta2*seta2;
380 double dphi = getDeltaPhi(phi1,phi2);
381 double sdphi = sphi1*sphi1+sphi2*sphi2;
382 double dipt_1 = ipt - ipt1;
383 if (useAbsPt) dipt_1 = fabs(ipt) - fabs(ipt1);
384 double sdipt_1 = sipt1*sipt1;
385 double dipt_2 = ipt - ipt2;
386 if (useAbsPt) dipt_2 = fabs(ipt) - fabs(ipt2);
387 double sdipt_2 = sipt2*sipt2;
388
389 double chi2 = 0.0;
390 ndof = 0;
391 if (sdeta != 0) { chi2 += deta*deta/sdeta; ndof++; }
392 if (sdphi != 0) { chi2 += dphi*dphi/sdphi; ndof++; }
393 if (sdipt_1 != 0) { chi2 += dipt_1*dipt_1/sdipt_1; ndof++; }
394 if (sdipt_2 != 0) { chi2 += dipt_2*dipt_2/sdipt_2; ndof++; }
395
396 if (ndof == 0) return 1.0e30;
397 else return chi2;
398 }
399
400
401 double getStdChi2(int& ndof,
402 double eta1, double seta1, double phi1, double sphi1, double qOvpt1, double sqOvpt1,
403 double eta2, double seta2, double phi2, double sphi2, double qOvpt2, double sqOvpt2, bool useAbsPt) {
404
405 double deta = getDeltaEta(eta1,eta2);
406 double sdeta = seta1*seta1+seta2*seta2;
407 double dphi = getDeltaPhi(phi1,phi2);
408 double sdphi = sphi1*sphi1+sphi2*sphi2;
409 double dipt = qOvpt1 - qOvpt2;
410 if (useAbsPt) dipt = fabs(qOvpt1) - fabs(qOvpt2);
411 double sdipt = sqOvpt1*sqOvpt1+sqOvpt2*sqOvpt2;
412
413 double chi2 = 0.0;
414 ndof = 0;
415 if (sdeta != 0) { chi2 += deta*deta/sdeta; ndof++; }
416 if (sdphi != 0) { chi2 += dphi*dphi/sdphi; ndof++; }
417 if (sdipt != 0) { chi2 += dipt*dipt/sdipt; ndof++; }
418
419 if (ndof == 0) return 1.0e30;
420 else return chi2;
421 }
422
423
424
425
426}//muCombUtil
M_PI
#define M_PI
Definition ActiveFraction.h:11
eta
Scalar eta() const
pseudorapidity method
Definition AmgMatrixBasePlugin.h:83
phi
Scalar phi() const
phi method
Definition AmgMatrixBasePlugin.h:67
vec
std::vector< size_t > vec
Definition CombinationsGeneratorTest.cxx:9
MuonFeature::saddress
int saddress(void) const
Definition MuonFeature.h:47
MuonFeature::pt
float pt(void) const
Definition MuonFeature.h:48
MuonFeature::phi
float phi(void) const
Definition MuonFeature.h:51
MuonFeature::radius
float radius(void) const
Definition MuonFeature.h:49
MuonFeature::eta
float eta(void) const
Definition MuonFeature.h:50
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition comparitor.cxx:525
add
bool add(const std::string &hname, TKey *tobj)
Definition fastadd.cxx:55
ts
int ts
Definition globals.cxx:24
muCombUtil::whichECRegion
int whichECRegion(const float eta, const float phi)
utility function for getMuFastRes
Definition muCombUtil.cxx:117
muCombUtil::getIDSCANRes
double getIDSCANRes(std::vector< double > barrelvec, std::vector< double > ec1vec, std::vector< double > ec2vec, std::vector< double > ec3vec, std::vector< double > ec4vec, double pt_id, double eta_id)
Get parametrized IDSCAN 1/pt resolution.
Definition muCombUtil.cxx:190
muCombUtil::getChi2
double getChi2(int &ndof, double ipt, double eta1, double seta1, double phi1, double sphi1, double ipt1, double sipt1, double eta2, double seta2, double phi2, double sphi2, double ipt2, double sipt2, bool useAbsPt)
Get OLD style (i.e. muFast time) Chi2.
Definition muCombUtil.cxx:374
muCombUtil::getCombinedAverage
double getCombinedAverage(double p1, double sp1, double p2, double sp2)
Get weighted mean.
Definition muCombUtil.cxx:362
muCombUtil::getG4ExtEtaRes
double getG4ExtEtaRes(double pt, double eta)
Get parametrized Geant4 Eta resolution (extrapolated)
Definition muCombUtil.cxx:288
muCombUtil::getStdChi2
double getStdChi2(int &ndof, double eta1, double seta1, double phi1, double sphi1, double qOvpt1, double sqOvpt1, double eta2, double seta2, double phi2, double sphi2, double qOvpt2, double sqOvpt2, bool useAbsPt)
Get Std Chi2.
Definition muCombUtil.cxx:401
muCombUtil::getCombinedAverageSigma
double getCombinedAverageSigma(double sp1, double sp2)
Get sigma of weighted mean.
Definition muCombUtil.cxx:369
muCombUtil::getDeltaEta
double getDeltaEta(double eta1, double eta2)
Get DeltaEta.
Definition muCombUtil.cxx:352
muCombUtil::getMuFastPhiRes
double getMuFastPhiRes(const MuonFeature *feature)
Get parametrized muFast Phi resolution (extrapolated)
Definition muCombUtil.cxx:250
muCombUtil::getMuFastRes
double getMuFastRes(const std::vector< double > &vec, const MuonFeature *feature)
Definition muCombUtil.cxx:42
muCombUtil::getDeltaR
double getDeltaR(double eta1, double phi1, double eta2, double phi2)
Get DeltaR.
Definition muCombUtil.cxx:356
muCombUtil::setMuFastRes
void setMuFastRes(std::vector< double > &vec, double p1, double p2, double p3, double p4, double p5, double p6)
utility vector set
Definition muCombUtil.cxx:24
muCombUtil::getMuFastEtaRes
double getMuFastEtaRes(const MuonFeature *feature)
Get parametrized muFast Eta resolution (extrapolated)
Definition muCombUtil.cxx:211
muCombUtil::setIDSCANRes
void setIDSCANRes(std::vector< double > &vec, double p1, double p2)
utility vector set
Definition muCombUtil.cxx:35
muCombUtil::getDeltaPhi
double getDeltaPhi(double phi1, double phi2)
Get DeltaPhi.
Definition muCombUtil.cxx:345
muCombUtil::getG4ExtPhiRes
double getG4ExtPhiRes(double pt, double eta)
Get parametrized Geant4 Phi resolution (extrapolated)
Definition muCombUtil.cxx:317
Generated on for ATLAS Offline Software by doxygen 1.14.0