001
002
003
004
005
006
007
008
009
010
011
012 #include "AtlfastAlgs/MuonAcceptor.h"
013
014 #include "CLHEP/Random/JamesRandom.h"
015 #include "CLHEP/Random/RandGauss.h"
016 #include "CLHEP/Units/SystemOfUnits.h"
017
018 #include <iostream>
019 #include <fstream>
020 #include <iostream>
021 #include <sstream>
022 #include <cmath>
023
024 #include "PathResolver/PathResolver.h"
025
026
027
028 namespace Atlfast {
029
030 using std::abs;
031 using std::pair;
032 using std::make_pair;
033
034
035
036
037
038
039
040
041
042
043
044 muEffdata::muEffdata(){}
045 muEffdata::~muEffdata(){}
046
047
048 void muEffdata::Print(std::ostream* out){
049 *out <<std::setw(12)<<std::setprecision(3)<<GetPt()
050 <<std::setw(12)<<std::setprecision(3)<<GetEtaMin()
051 <<std::setw(12)<<std::setprecision(3)<<GetEtaMax()
052 << std::endl;
053 }
054
055 void muEffdata::SetPt (double pt) { m_pt = pt ;}
056
057 void muEffdata::SetEtaMin (double etamin) { m_etamin = etamin ;}
058
059 void muEffdata::SetEtaMax (double etamax) { m_etamax = etamax ;}
060
061 void muEffdata::SetPtEtaNum(std::vector<int> ptEtaNum){m_ptEtaNum=ptEtaNum;}
062
063 void muEffdata::SetPtEtaNumeDeno(std::vector<std::pair<int, int> > VecPaire)
064 {m_VecPaire = VecPaire;}
065
066 double muEffdata::GetPt() {return m_pt ;}
067
068 double muEffdata::GetEtaMin(){return m_etamin;}
069
070 double muEffdata::GetEtaMax(){return m_etamax;}
071
072 int muEffdata::GetPtEtaNum(std::vector<int> &ptEtaNum){
073
074 int SizeOfVectorOfmuEffdata = m_ptEtaNum.size();
075 ptEtaNum = m_ptEtaNum;
076 return SizeOfVectorOfmuEffdata;
077
078 }
079
080 void muEffdata::GetPtEtaNumeDeno(double &pt,double &etamin,double &etamax,
081 std::vector<std::pair<int, int> > &VecPaire){
082
083 pt = m_pt ;
084 etamin = m_etamin ;
085 etamax = m_etamax ;
086 VecPaire = m_VecPaire;
087
088 }
089
090 muEffdata::muEffdata(double &pt,double &etamin,double &etamax,
091 std::vector<std::pair<int, int> > &VecPaire){
092
093 pt = m_pt ;
094 etamin = m_etamin ;
095 etamax = m_etamax ;
096 VecPaire = m_VecPaire;
097
098 }
099
100
101
102
103
104
105
106 MuonAcceptor::MuonAcceptor(int muSmearKey, MsgStream& log) : m_log(log) {
107
108 m_log << MSG::DEBUG << "MuonAcceptor::MuonAcceptor" << endreq;
109
110 m_randomEngine = new HepJamesRandom(12345);
111
112
113
114 std::string filename = muSmearKey == 1 ?
115 "atlfastDatafiles/MuonSpectrometerEff.txt" :
116 "atlfastDatafiles/MuonCombinedEff.txt";
117 filename = PathResolver::find_file (filename, "DATAPATH");
118 m_log << MSG::DEBUG << "Opening file: " << filename << endreq;
119 std::ifstream InputFile;
120 InputFile.open(filename.c_str());
121 if(!InputFile) {
122 m_log << MSG::ERROR << "Error: file could not be opened" << endreq;
123 exit(1);
124 }
125
126 std::vector<muEffdata>VectorOfmuEffdata;
127 std::vector<std::pair<int, int> > ptEtaNumeDeno;
128 VectorOfmuEffdata.clear();
129 ptEtaNumeDeno.clear();
130
131 m_muEffdataObject.clear();
132 m_ptValues.clear();
133 m_etaminValues.clear();
134 m_etamaxValues.clear();
135
136
137
138 int Istatus = 1;
139 while ( Istatus == 1 ) {
140
141
142
143
144 Istatus = muReadEff(InputFile,VectorOfmuEffdata);
145
146 if (Istatus == 1){
147
148 std::vector<int> ptEtaNum;
149 int SizeOfVectorOfmuEffdata = VectorOfmuEffdata.size() ;
150 for (int Item=0; Item <SizeOfVectorOfmuEffdata ; Item++){
151 double pt;double etamin;double etamax;
152 VectorOfmuEffdata[Item].GetPtEtaNumeDeno(pt,etamin,etamax,ptEtaNumeDeno);
153 m_muEffdataObject.push_back(VectorOfmuEffdata[Item]);
154 m_ptValues.push_back(pt);
155 m_etaminValues.push_back(etamin);
156 m_etamaxValues.push_back(etamax);
157 }
158 }
159 }
160
161 InputFile.close();
162 m_log << MSG::DEBUG << "Done with MuonAcceptor constructor" << endreq;
163 }
164
165 bool MuonAcceptor::accept( const ReconstructedParticle &particle, MsgStream &log ){
166
167 log << MSG::DEBUG << "MuonAcceptor::accept" << endreq;
168
169
170 double pTinGeV = particle.pT() / GeV;
171
172 double muEff = muEfficiency(pTinGeV,
173 particle.eta(),
174 particle.phi()).first;
175 log << MSG::DEBUG << "pT,eta,phi,efficiency = " << pTinGeV
176 << "," << particle.eta() << ","
177 << particle.phi() << "," << muEff << endreq;
178 double random = m_randomEngine->flat();
179 log << MSG::DEBUG << "random = " << random << endreq;
180
181 return ( random < muEff) ? true : false;
182 }
183
184
185
186 int MuonAcceptor::muReadEff(std::ifstream& InputFile,
187 std::vector<muEffdata>& VectorOfmuEffdata){
188
189 VectorOfmuEffdata.clear();
190
191 std::vector<int> ptEtaNume;
192 std::vector<int> ptEtaDeno;
193
194 std::vector<int> ptEtaNume1;
195 std::vector<int> ptEtaDeno1;
196
197 std::vector<int> ptEtaNume2;
198 std::vector<int> ptEtaDeno2;
199
200 std::vector<int> ptEtaNume3;
201 std::vector<int> ptEtaDeno3;
202
203 std::vector<int> ptEtaNume4;
204 std::vector<int> ptEtaDeno4;
205
206 std::vector<std::pair<int, int> > VecPaire;
207 std::pair<int, int> mapaire;
208
209 std::vector<std::pair<int, int> > VecPaire1;
210 std::pair<int, int> mapaire1;
211
212 std::vector<std::pair<int, int> > VecPaire2;
213 std::pair<int, int> mapaire2;
214
215 std::vector<std::pair<int, int> > VecPaire3;
216 std::pair<int, int> mapaire3;
217
218 std::vector<std::pair<int, int> > VecPaire4;
219 std::pair<int, int> mapaire4;
220
221 ptEtaNume.clear();
222 ptEtaDeno.clear();
223 ptEtaNume1.clear();
224 ptEtaDeno1.clear();
225 ptEtaNume2.clear();
226 ptEtaDeno2.clear();
227 ptEtaNume3.clear();
228 ptEtaDeno3.clear();
229 ptEtaNume4.clear();
230 ptEtaDeno4.clear();
231
232 muEffdata pmuEffdata;
233
234 if(!InputFile) return 0;
235 if(InputFile.eof()) return 0;
236
237 double pt;
238 double etamin,etamax;
239 int nBin;
240 int nnum,ndeno;
241
242 if(InputFile.eof()) return 0;
243 InputFile >> pt;
244 if(InputFile.eof()) return 0;
245 InputFile >>etamin>>etamax>>nBin;
246 if(InputFile.eof()) return 0;
247
248 pmuEffdata.SetPt(pt);
249 pmuEffdata.SetEtaMin(etamin );
250 pmuEffdata.SetEtaMax(etamax);
251
252 for (int ibin=0; ibin < nBin ; ibin++){
253 if(InputFile.eof()) return 0;
254 InputFile >> nnum ;
255 ptEtaNume.push_back(nnum);
256 }
257 if(InputFile.eof()) return 0;
258 for (int ibin=0; ibin < nBin ; ibin++){
259 if(InputFile.eof()) return 0;
260 InputFile >> ndeno ;
261 ptEtaDeno.push_back(ndeno);
262 }
263 pmuEffdata.SetPtEtaNum(ptEtaNume);
264 pmuEffdata.SetPtEtaNum(ptEtaDeno);
265
266 int SizeOfptEtaNume =ptEtaNume.size();
267 for(int i=0;i<SizeOfptEtaNume;i++){
268 mapaire.first=ptEtaNume[i];
269 mapaire.second = ptEtaDeno[i];
270 VecPaire.push_back(mapaire);
271 }
272 pmuEffdata.SetPtEtaNumeDeno(VecPaire);
273
274
275 if(InputFile.eof()) return 0;
276
277 VectorOfmuEffdata.push_back(pmuEffdata);
278
279 int nBin1 =0;
280 double etamin1=0;
281 double etamax1=0;
282 int ndeno1=0;
283 int nnum1=0;
284
285 InputFile >>etamin1>>etamax1>>nBin1;
286 if(InputFile.eof()) return 0;
287
288 pmuEffdata.SetEtaMin(etamin1 );
289 pmuEffdata.SetEtaMax(etamax1);
290
291 for (int Item=0; Item < nBin1 ; Item++){
292 if(InputFile.eof()) return 0;
293 InputFile >> nnum1 ;
294 ptEtaNume1.push_back(nnum1);
295 }
296 if(InputFile.eof()) return 0;
297 for (int Item=0; Item < nBin1 ; Item++){
298 if(InputFile.eof()) return 0;
299 InputFile >> ndeno1 ;
300 ptEtaDeno1.push_back(ndeno1);
301 }
302 pmuEffdata.SetPtEtaNum(ptEtaNume1);
303 pmuEffdata.SetPtEtaNum(ptEtaDeno1);
304
305
306 if(InputFile.eof()) return 0;
307
308 int SizeOfptEtaNume1 =ptEtaNume1.size();
309 for(int i=0;i<SizeOfptEtaNume1;i++){
310 mapaire1.first=ptEtaNume1[i];
311 mapaire1.second = ptEtaDeno1[i];
312 VecPaire1.push_back(mapaire1);
313 }
314 pmuEffdata.SetPtEtaNumeDeno(VecPaire1);
315 VectorOfmuEffdata.push_back(pmuEffdata);
316
317
318 int nBin2 =0;
319 double etamin2=0;
320 double etamax2=0;
321 int ndeno2=0;
322 int nnum2=0;
323
324 InputFile >>etamin2>>etamax2>>nBin2;
325 if(InputFile.eof()) return 0;
326
327 pmuEffdata.SetEtaMin(etamin2 );
328 pmuEffdata.SetEtaMax(etamax2);
329
330 for (int Item=0; Item < nBin2 ; Item++){
331 if(InputFile.eof()) return 0;
332 InputFile >> nnum2 ;
333 ptEtaNume2.push_back(nnum2);
334 }
335 if(InputFile.eof()) return 0;
336 for (int Item=0; Item < nBin2 ; Item++){
337 if(InputFile.eof()) return 0;
338 InputFile >> ndeno2 ;
339 ptEtaDeno2.push_back(ndeno2);
340 }
341 pmuEffdata.SetPtEtaNum(ptEtaNume2);
342 pmuEffdata.SetPtEtaNum(ptEtaDeno2);
343
344
345 if(InputFile.eof()) return 0;
346
347 int SizeOfptEtaNume2 =ptEtaNume2.size();
348 for(int i=0;i<SizeOfptEtaNume2;i++){
349 mapaire2.first=ptEtaNume2[i];
350 mapaire2.second = ptEtaDeno2[i];
351 VecPaire2.push_back(mapaire2);
352 }
353 pmuEffdata.SetPtEtaNumeDeno(VecPaire2);
354 VectorOfmuEffdata.push_back(pmuEffdata);
355
356
357
358 int nBin3 =0;
359 double etamin3=0;
360 double etamax3=0;
361 int ndeno3=0;
362 int nnum3=0;
363
364 InputFile >>etamin3>>etamax3>>nBin3;
365 if(InputFile.eof()) return 0;
366
367 pmuEffdata.SetEtaMin(etamin3 );
368 pmuEffdata.SetEtaMax(etamax3);
369
370 for (int Item=0; Item < nBin3 ; Item++){
371 if(InputFile.eof()) return 0;
372 InputFile >> nnum3 ;
373 ptEtaNume3.push_back(nnum3);
374 }
375 if(InputFile.eof()) return 0;
376 for (int Item=0; Item < nBin3 ; Item++){
377 if(InputFile.eof()) return 0;
378 InputFile >> ndeno3 ;
379 ptEtaDeno3.push_back(ndeno3);
380 }
381 pmuEffdata.SetPtEtaNum(ptEtaNume3);
382 pmuEffdata.SetPtEtaNum(ptEtaDeno3);
383
384
385 if(InputFile.eof()) return 0;
386
387 int SizeOfptEtaNume3 =ptEtaNume3.size();
388 for(int i=0;i<SizeOfptEtaNume3;i++){
389 mapaire3.first=ptEtaNume3[i];
390 mapaire3.second = ptEtaDeno3[i];
391 VecPaire3.push_back(mapaire3);
392 }
393 pmuEffdata.SetPtEtaNumeDeno(VecPaire3);
394 VectorOfmuEffdata.push_back(pmuEffdata);
395
396
397
398
399 int nBin4 =0;
400 double etamin4=0;
401 double etamax4=0;
402 int ndeno4=0;
403 int nnum4=0;
404
405 InputFile >>etamin4>>etamax4>>nBin4;
406 if(InputFile.eof()) return 0;
407
408 pmuEffdata.SetEtaMin(etamin4 );
409 pmuEffdata.SetEtaMax(etamax4);
410
411 for (int Item=0; Item < nBin4 ; Item++){
412 if(InputFile.eof()) return 0;
413 InputFile >> nnum4 ;
414 ptEtaNume4.push_back(nnum4);
415 }
416 if(InputFile.eof()) return 0;
417 for (int Item=0; Item < nBin4 ; Item++){
418 if(InputFile.eof()) return 0;
419 InputFile >> ndeno4 ;
420 ptEtaDeno4.push_back(ndeno4);
421 }
422 pmuEffdata.SetPtEtaNum(ptEtaNume4);
423 pmuEffdata.SetPtEtaNum(ptEtaDeno4);
424
425
426 if(InputFile.eof()) return 0;
427
428 int SizeOfptEtaNume4 =ptEtaNume4.size();
429 for(int i=0;i<SizeOfptEtaNume4;i++){
430 mapaire4.first=ptEtaNume4[i];
431 mapaire4.second = ptEtaDeno4[i];
432 VecPaire4.push_back(mapaire4);
433 }
434 pmuEffdata.SetPtEtaNumeDeno(VecPaire4);
435 VectorOfmuEffdata.push_back(pmuEffdata);
436
437
438 return 1;
439 }
440
441
442
443
444
445
446 std::pair<double,double>
447 MuonAcceptor::muEfficiency(double ptIn,double etaIn, double phi)
448 {
449
450 double efficiency =0;
451 double errEfficiency =0;
452
453 double eta = fabs(etaIn);
454 double pt = fabs(ptIn);
455
456
457 double pt1 = 0.;
458 double pt2 = 0.;
459
460
461 std::pair<double,double>ptLowHigh;
462 ptLowHigh= muGetPtLowHigh(pt);
463 pt1 = ptLowHigh.first;
464 pt2 = ptLowHigh.second;
465
466 if(pt >= pt1 && pt < pt2){
467 std::pair<double,double> eff1_errff1;
468 std::pair<double,double> eff2_errff2;
469
470 eff1_errff1 = muEff(pt1,eta,phi);
471 eff2_errff2 = muEff(pt2,eta,phi);
472
473 double eff1 = eff1_errff1.first ;
474 double eeff1 = eff1_errff1.second;
475
476 double eff2 = eff2_errff2.first;
477 double eeff2 = eff2_errff2.second;
478
479 efficiency = eff2 *log(pt/pt1)/log(pt2/pt1);
480 efficiency = efficiency + eff1 *log(pt2/pt)/log(pt2/pt1);
481 errEfficiency = pow((eeff2*log(pt/pt1)/log(pt2/pt1)),2);
482 errEfficiency = errEfficiency +pow((eeff1*log(pt2/pt)/log(pt2/pt1)),2);
483 errEfficiency = sqrt(errEfficiency);
484 }
485
486 if(efficiency < 0){
487 efficiency =0.;
488 errEfficiency=0.;
489 }
490 if(efficiency > 1){
491 efficiency =1.;
492 errEfficiency=0.;
493 }
494
495 return make_pair(efficiency,errEfficiency);
496 }
497
498
499
500
501
502
503
504
505 std::pair<double,double>
506 MuonAcceptor::muEff(double ptIn,double etaIn, double phi)
507 {
508 double eta = fabs(etaIn);
509 double pt = fabs(ptIn);
510
511 double xEtaMinOut = 0.;
512 double xEtaMaxOut = 2.5;
513
514 int nNumOut = 0;
515 int nDenoOut = 0;
516
517 int nNumTotOut = 0;
518 int nDenoTotOut= 0;
519
520 int iEtaBin = muGetEtaBin(eta);
521
522 std::vector<std::pair<int, int> > ptEtaNumeDeno;
523 ptEtaNumeDeno.clear();
524 ptEtaNumeDeno= muEffGetPaire(pt,xEtaMinOut,xEtaMaxOut);
525
526 nNumOut=ptEtaNumeDeno[iEtaBin].first;
527 nDenoOut=ptEtaNumeDeno[iEtaBin].second;
528
529 int nBinOut = ptEtaNumeDeno.size();
530
531 for(unsigned int i=0 ; i<ptEtaNumeDeno.size() ; i++){
532 nNumTotOut = nNumTotOut+ptEtaNumeDeno[i].first;
533 nDenoTotOut = nDenoTotOut+ptEtaNumeDeno[i].second;
534 }
535
536 double eff1 = double(nNumOut)/double(nDenoOut);
537 double eeff1 = sqrt(eff1*(1.-eff1)/double(nDenoOut));
538 double delta = (xEtaMaxOut-xEtaMinOut)/double(nBinOut) ;
539 double xcbin1 = xEtaMinOut+delta/2.+delta*double(iEtaBin) ;
540
541 muPhiEfficiency(pt,xcbin1,phi,eff1,eeff1);
542
543 return make_pair(eff1,eeff1);
544 }
545
546
547 int MuonAcceptor::muGetEtaBin(double etaIn){
548
549 int iEtaBin = -1;
550 double eta = fabs(etaIn);
551
552 double xMaxOut =2.5;
553 double xMinOut = 0.;
554 double pt=50.;
555
556 std::vector<std::pair<int, int> > ptEtaNumeDeno;
557 ptEtaNumeDeno.clear();
558 ptEtaNumeDeno= muEffGetPaire(pt,xMinOut,xMaxOut);
559
560 int NbinOut = ptEtaNumeDeno.size();
561
562 double delta=double(xMaxOut-xMinOut)/double(NbinOut);
563 iEtaBin = int(floor(eta/delta));
564 if( iEtaBin >= NbinOut) iEtaBin= NbinOut-1;
565
566 return iEtaBin;
567 }
568
569
570
571 void MuonAcceptor::muPhiEfficiency(double pt,double etaIn,double phi,double &eff,double &eeff){
572
573 double pi = 3.1415927;
574 double phig = phi;
575 if (phi < 0.) {
576 phig = pi*2 + phi;
577 }
578
579 double effin = eff;
580
581
582 double eta = fabs(etaIn);
583
584
585
586
587
588 double etamin1 = 0.;
589 double etamax1 = 0.15;
590 if (eta >= etamin1 && eta < etamax1) {
591 double pion4 = pi/4.;
592 double pion8 = pi/8.;
593 double phig2 = fmod(phig,pion4);
594 if (phig2 >= pion8) phig2 = pion4-phig2;
595
596 std::pair<double,double> eff_erreff;
597 eff_erreff = getPhiEff(pt,effin,phig2,etamin1,etamax1);
598 eff = eff_erreff.first ;
599 eeff = eff_erreff.second;
600 }
601
602
603
604
605
606 double etamin2 = 0.3;
607 double etamax2 = 0.4;
608 double etamin22 = 0.8;
609 double etamax22 = 0.9;
610 if ((eta >= etamin2 && eta < etamax2) ||
611 (eta >= etamin22 && eta < etamax22) ) {
612 std::pair<double,double> eff_erreff;
613 eff_erreff = getPhiEff(pt,effin,phig,etamin2,etamax2);
614 eff = eff_erreff.first ;
615 eeff = eff_erreff.second;
616 }
617
618
619
620 double etamin3 =1.05;
621 double etamax3 =1.15;
622 if (eta >= etamin3 && eta < etamax3) {
623
624 std::pair<double,double> eff_erreff;
625 eff_erreff = getPhiEff(pt,effin,phig,etamin3,etamax3);
626 eff = eff_erreff.first ;
627 eeff = eff_erreff.second;
628 }
629
630
631
632
633 double etamin4 =1.20;
634 double etamax4 =1.25;
635
636 if (eta >= etamin4 && eta < etamax4) {
637
638 std::pair<double,double> eff_erreff;
639 eff_erreff = getPhiEff(pt,effin,phig,etamin4,etamax4);
640 eff = eff_erreff.first ;
641 eeff = eff_erreff.second;
642 }
643
644 }
645
646
647
648 int MuonAcceptor::muGetPhiBin(double phi){
649
650 double pi = 3.1415927;
651 int iPhiBin = -1;
652
653 double xMaxOut =360;
654 double xMinOut = 0.;
655 int NbinOut = 20;
656
657 double phig = phi;
658 if (phi < 0.) {
659 phig = pi*2 + phi;
660 }
661 double phideg = phig*180./ pi;
662
663 double delta=double(xMaxOut-xMinOut)/double(NbinOut);
664 iPhiBin = int(floor(phideg/delta));
665 if( iPhiBin >= NbinOut) iPhiBin= NbinOut-1;
666
667 return iPhiBin;
668 }
669
670 int MuonAcceptor::muGetPhiBin1(double phi){
671
672 double pi = 3.1415927;
673 int iPhiBin = -1;
674
675 double xMaxOut = 22.5;
676 double xMinOut = 0.;
677 int NbinOut = 20;
678
679 double phig = phi;
680 if (phi < 0.) {
681 phig = pi*2 + phi;
682 }
683 double phideg = phig*180./ pi;
684
685 double delta=double(xMaxOut-xMinOut)/double(NbinOut);
686 iPhiBin = int(floor(phideg/delta));
687 if( iPhiBin >= NbinOut) iPhiBin= NbinOut-1;
688
689 return iPhiBin;
690 }
691
692
693 std::pair<double,double>
694 MuonAcceptor::getPhiEff(double pt,double effin,double phig,double etamin,double etamax)
695 {
696 double eff = 1.;
697 double eeff =0.;
698
699 double phivalue[21];
700 double efficiency[21];
701 double errefficiency[21];
702
703 for (int j = 0; j < 21; j++) {
704 efficiency[j] = 0.;
705 errefficiency[j] = 0.;
706 phivalue[j] = 0.;
707 }
708
709 double pi = 3.1415927;
710 double degrad = pi/180.;
711
712 for (int j = 0; j < 21; j++) {
713 phivalue[j] = j*18.*degrad;
714 if(etamin ==0. && etamax ==0.15)phivalue[j] = j*1.125*degrad;
715 }
716
717
718
719 std::pair<double,double> eff_errff;
720 double phiLow = phivalue[0];
721 double phiHigh = phivalue[20];
722 eff_errff = muPhiEff(pt,etamin,etamax,phiLow,phiHigh);
723 efficiency[20] = eff_errff.first ;
724 errefficiency[20] = eff_errff.second;
725
726
727 double delta2 = 8.9*degrad;
728 if(etamin ==0. && etamax ==0.15)delta2 = 0.55*degrad;
729 for (int j = 0; j < 20; j++) {
730 double phicenter = (phivalue[j+1]+phivalue[j])/2.;
731 double phimin = phicenter - delta2;
732 double phimax = phicenter + delta2;
733 eff_errff = muPhiEff(pt,etamin,etamax,phimin,phimax);
734 efficiency[j] = eff_errff.first ;
735 errefficiency[j] = eff_errff.second;
736 }
737
738 double averaefficiency = effin/efficiency[20];
739 double factornorm = averaefficiency;
740 if(factornorm > 1. || factornorm < 0.9)factornorm =1;
741
742 if (phig >= phivalue[0] && phig <= phivalue[1]) {
743 eff = efficiency[0] * factornorm;
744 eeff = errefficiency[0];
745 }
746 for (int j = 1; j < 20; j++) {
747 if (phig > phivalue[j] && phig <= phivalue[j+1]) {
748 eff = efficiency[j] * factornorm;
749 eeff= errefficiency[j];
750 }
751 }
752 if (phig > phivalue[20] && phig <= phivalue[21]) {
753 eff = efficiency[19] * factornorm;
754 eeff= errefficiency[19];
755 }
756
757 if (eff < 0.) {
758 eff = 0.;
759 eeff = 0.;
760 }
761 if (eff > 1.) {
762 eff = 1.;
763 eeff =0.;
764 }
765
766 return make_pair(eff,eeff);
767 }
768
769
770 std::pair<double,double>
771 MuonAcceptor::muPhiEff(double pt,double etamin,double etamax,double phimin,double phimax){
772
773 int NnumTot=0 ;
774 int NdenoTot=0;
775
776 double eff =0;
777 double eeff=0;
778
779 int ibinmin;
780 int ibinmax;
781
782 if(etamin==0. && etamax==0.15){
783 ibinmin = muGetPhiBin1(phimin) ;
784 ibinmax = muGetPhiBin1(phimax) ;
785 }
786 else{
787 ibinmin = muGetPhiBin(phimin) ;
788 ibinmax = muGetPhiBin(phimax) ;
789 }
790 std::vector<std::pair<int, int> > ptEtaNumeDeno;
791 ptEtaNumeDeno.clear();
792 ptEtaNumeDeno= muEffGetPaire(pt,etamin,etamax);
793
794 for (int ibin=ibinmin;ibin<=ibinmax;ibin++){
795 NnumTot =NnumTot+ptEtaNumeDeno[ibin].first;
796 NdenoTot=NdenoTot+ptEtaNumeDeno[ibin].second;
797 }
798
799 eff =double(NnumTot)/double(NdenoTot);
800 eeff=sqrt(eff*(1.-eff)/double(NdenoTot));
801
802 return make_pair(eff,eeff);
803 }
804
805
806 std::vector<std::pair<int, int> >
807 MuonAcceptor::muEffGetPaire(double pt,double etamin,double etamax){
808
809 std::vector<std::pair<int, int> > vecPaireNumeDeno;
810 vecPaireNumeDeno.clear();
811
812 std::vector<std::pair<int, int> > ptEtaNumeDeno;
813 ptEtaNumeDeno.clear();
814
815 int iposition = muGetPaireNumeDenoPosi(pt,etamin,etamax);
816
817 if(iposition != -1){
818 int i=iposition;
819 m_muEffdataObject[i].GetPtEtaNumeDeno(pt,etamin,etamax,ptEtaNumeDeno);
820 vecPaireNumeDeno=ptEtaNumeDeno;
821 }
822
823 return vecPaireNumeDeno;
824 }
825
826 int MuonAcceptor::muGetPaireNumeDenoPosi(double pt,double etamin,double etamax){
827
828 int iposition =-1;
829 const double epsilon=1.e-6;
830 for(unsigned int i= 0;i<m_ptValues.size();i++){
831 if( fabs(pt-m_ptValues[i]) < epsilon &&
832 fabs(etamin - m_etaminValues[i]) < epsilon &&
833 fabs( etamax - m_etamaxValues[i]) < epsilon
834 )
835 iposition=i;
836 }
837 return iposition;
838 }
839
840 std::pair<double,double>
841 MuonAcceptor::muGetPtLowHigh(double pt){
842
843 std::pair<double,double> ptLowHigh;
844
845 int sizeOfptvalues = m_ptValues.size();
846 double ptmin= m_ptValues.at(0);
847 double ptmax= m_ptValues.at(sizeOfptvalues-1);
848
849 for (int i=0; i <sizeOfptvalues-1 ; i++){
850 ptmin= m_ptValues[i];
851 ptmax= m_ptValues[i+1];
852 if(pt >=ptmin && ptmin<ptmax){
853 double ptlow =m_ptValues[i] ;
854 double pthigh =m_ptValues[i+1] ;
855 ptLowHigh=make_pair(ptlow,pthigh);
856 }
857 }
858 if(pt>=m_ptValues.at(sizeOfptvalues-1)){
859 double ptlow = m_ptValues.at(sizeOfptvalues-2);
860 double pthigh = m_ptValues.at(sizeOfptvalues-1);
861 ptLowHigh=make_pair(ptlow,pthigh);
862 }
863 if(pt <=m_ptValues.at(0)){
864 double ptlow = m_ptValues.at(0);
865 double pthigh = m_ptValues.at(0) ;
866 ptLowHigh=make_pair(ptlow,pthigh);
867 }
868 return ptLowHigh ;
869 }
870
871 }
872
| Due to the LXR bug, the updates fail sometimes to remove references to deleted files. The Saturday's full rebuilds fix these problems |
|
This page was automatically generated by the
LXR engine.
|
|