Atlfast::MuonAcceptor Class Reference

#include <MuonAcceptor.h>

Inheritance diagram for Atlfast::MuonAcceptor:

Inheritance graph
[legend]
Collaboration diagram for Atlfast::MuonAcceptor:

Collaboration graph
[legend]
List of all members.

Public Member Functions

 MuonAcceptor (int muSmearKey, MsgStream &log)
virtual ~MuonAcceptor ()
virtual bool accept (const ReconstructedParticle &particle, MsgStream &log)
int muReadEff (std::ifstream &InputFile, vector< muEffdata > &VectorOfmuEffdata)
vector< pair< int, int > > muEffGetPaire (double pt, double etamin, double etamax)
pair< double, double > muGetPtLowHigh (double pt)
int muGetPaireNumeDenoPosi (double pt, double etamin, double etamax)
pair< double, double > muEfficiency (double ptIn, double etaIn, double phi)
pair< double, double > muEff (double ptIn, double etaIn, double phi)
int muGetEtaBin (double etaIn)
int muGetPhiBin (double phi)
int muGetPhiBin1 (double phi)
pair< double, double > muPhiEff (double pt, double etamin, double etamax, double phimin, double phimax)
pair< double, double > getPhiEff (double pt, double effin, double phig, double etamin, double etamax)
void muPhiEfficiency (double pt, double etaIn, double phi, double &eff, double &eeff)

Private Attributes

vector< muEffdatam_muEffdataObject
vector< double > m_ptValues
vector< double > m_etaminValues
vector< double > m_etamaxValues
MsgStream & m_log
HepRandomEngine * m_randomEngine

Detailed Description

Definition at line 102 of file MuonAcceptor.h.


Constructor & Destructor Documentation

Atlfast::MuonAcceptor::MuonAcceptor ( int  muSmearKey,
MsgStream &  log 
)

Definition at line 92 of file MuonAcceptor.cxx.

00092                                                            : m_log(log) {
00093     
00094     m_log << MSG::DEBUG << "MuonAcceptor::MuonAcceptor" << endreq;
00095     // Initialize random number generator
00096     m_randomEngine = new HepJamesRandom(12345);
00097     
00098     //Open input file
00099 
00100     std::string filename = muSmearKey == 1 ? 
00101       "atlfastDatafiles/MuonSpectrometerEff.txt" : 
00102       "atlfastDatafiles/MuonCombinedEff.txt";
00103     filename = PathResolver::find_file (filename, "DATAPATH");
00104     m_log << MSG::DEBUG << "Opening file: " << filename << endreq;
00105     std::ifstream InputFile;
00106     InputFile.open(filename.c_str());
00107     if(!InputFile) { 
00108       m_log << MSG::ERROR << "Error: file could not be opened" << endreq;
00109       exit(1);
00110     }
00111     // 
00112     std::vector<muEffdata>VectorOfmuEffdata;
00113     std::vector<std::pair<int, int> > ptEtaNumeDeno;
00114     VectorOfmuEffdata.clear();
00115     ptEtaNumeDeno.clear();
00116     //
00117     m_muEffdataObject.clear();
00118     m_ptValues.clear();
00119     m_etaminValues.clear();
00120     m_etamaxValues.clear();
00121 
00122 
00123     //Loop on evts
00124     int Istatus = 1;
00125     while ( Istatus == 1 ) {
00126   
00127       //  Get evt data
00128       //
00129       Istatus = muReadEff(InputFile,VectorOfmuEffdata);
00130       //  Do something
00131       if (Istatus == 1){
00132 
00133         std::vector<int> ptEtaNum;
00134         int SizeOfVectorOfmuEffdata = VectorOfmuEffdata.size() ;
00135         for (int Item=0; Item <SizeOfVectorOfmuEffdata  ; Item++){
00136           double pt;double etamin;double etamax;
00137           VectorOfmuEffdata[Item].GetPtEtaNumeDeno(pt,etamin,etamax,ptEtaNumeDeno);
00138           m_muEffdataObject.push_back(VectorOfmuEffdata[Item]);
00139           m_ptValues.push_back(pt);
00140           m_etaminValues.push_back(etamin);
00141           m_etamaxValues.push_back(etamax);
00142         }
00143       }
00144     }
00145     //Close input file
00146     InputFile.close();
00147     m_log << MSG::DEBUG << "Done with MuonAcceptor constructor" << endreq;
00148   }

virtual Atlfast::MuonAcceptor::~MuonAcceptor (  )  [inline, virtual]

Definition at line 106 of file MuonAcceptor.h.

00106 { };


Member Function Documentation

bool Atlfast::MuonAcceptor::accept ( const ReconstructedParticle particle,
MsgStream &  log 
) [virtual]

Accept function, must be overridden in derived class

Implements Atlfast::IAcceptor.

Definition at line 150 of file MuonAcceptor.cxx.

00150                                                                                   {
00151     
00152     log << MSG::DEBUG << "MuonAcceptor::accept" << endreq;
00153 
00154     // The efficiencies are a function of pT in GeV
00155     double pTinGeV = particle.pT() / GeV;
00156  
00157     double muEff = muEfficiency(pTinGeV,
00158                                 particle.eta(),
00159                                 particle.phi()).first;
00160     log << MSG::DEBUG << "pT,eta,phi,efficiency = " << pTinGeV
00161         << "," << particle.eta() << "," 
00162         << particle.phi() << "," << muEff << endreq;
00163     double random = m_randomEngine->flat();
00164     log << MSG::DEBUG << "random = " << random << endreq;
00165     return ( random < muEff) ? true : false;
00166   }

int Atlfast::MuonAcceptor::muReadEff ( std::ifstream &  InputFile,
vector< muEffdata > &  VectorOfmuEffdata 
)

std::vector< std::pair< int, int > > Atlfast::MuonAcceptor::muEffGetPaire ( double  pt,
double  etamin,
double  etamax 
)

Definition at line 776 of file MuonAcceptor.cxx.

00776                                                                   {
00777     
00778     std::vector<std::pair<int, int> > vecPaireNumeDeno;
00779     vecPaireNumeDeno.clear();
00780     //
00781     std::vector<std::pair<int, int> > ptEtaNumeDeno;
00782     ptEtaNumeDeno.clear();
00783     // 
00784     int iposition = muGetPaireNumeDenoPosi(pt,etamin,etamax);
00785     
00786     if(iposition != -1){
00787       int i=iposition;
00788       m_muEffdataObject[i].GetPtEtaNumeDeno(pt,etamin,etamax,ptEtaNumeDeno);
00789       vecPaireNumeDeno=ptEtaNumeDeno;
00790     }
00791     //
00792     return vecPaireNumeDeno;
00793   }

std::pair< double, double > Atlfast::MuonAcceptor::muGetPtLowHigh ( double  pt  ) 

Definition at line 810 of file MuonAcceptor.cxx.

00810                                        {
00811 
00812     std::pair<double,double> ptLowHigh;
00813 
00814     int sizeOfptvalues   = m_ptValues.size();
00815     double ptmin= m_ptValues.at(0);
00816     double ptmax= m_ptValues.at(sizeOfptvalues-1);
00817 
00818     for (int i=0; i <sizeOfptvalues-1 ; i++){     
00819       ptmin= m_ptValues[i];
00820       ptmax= m_ptValues[i+1];
00821       if(pt >=ptmin && ptmin<ptmax){
00822         double ptlow  =m_ptValues[i] ;
00823         double pthigh =m_ptValues[i+1] ;
00824         ptLowHigh=make_pair(ptlow,pthigh);     
00825       } 
00826     }  
00827     if(pt>=m_ptValues.at(sizeOfptvalues-1)){
00828       double ptlow  = m_ptValues.at(sizeOfptvalues-2);
00829       double pthigh = m_ptValues.at(sizeOfptvalues-1);
00830       ptLowHigh=make_pair(ptlow,pthigh);
00831     }
00832     if(pt <=m_ptValues.at(0)){
00833       double ptlow  = m_ptValues.at(0);
00834       double pthigh = m_ptValues.at(0) ;
00835       ptLowHigh=make_pair(ptlow,pthigh);
00836     }
00837     return ptLowHigh ;
00838   }

int Atlfast::MuonAcceptor::muGetPaireNumeDenoPosi ( double  pt,
double  etamin,
double  etamax 
)

Definition at line 795 of file MuonAcceptor.cxx.

00795                                                                                {
00796     
00797     int iposition =-1;
00798     const double epsilon=1.e-6;
00799     for(unsigned int i= 0;i<m_ptValues.size();i++){
00800       if( fabs(pt-m_ptValues[i]) < epsilon &&
00801           fabs(etamin - m_etaminValues[i]) < epsilon &&
00802           fabs( etamax - m_etamaxValues[i]) < epsilon 
00803           )     
00804         iposition=i;
00805     }
00806     return iposition;
00807   }

std::pair< double, double > Atlfast::MuonAcceptor::muEfficiency ( double  ptIn,
double  etaIn,
double  phi 
)

Definition at line 427 of file MuonAcceptor.cxx.

00428   {
00429 
00430     double efficiency    =0;
00431     double errEfficiency =0; 
00432    
00433     double eta = fabs(etaIn);
00434     double pt  = fabs(ptIn);
00435 
00436     // Interpolate between pt points
00437     double pt1 = 0.;
00438     double pt2 = 0.;
00439 
00440     // Get the two pt points
00441     std::pair<double,double>ptLowHigh;
00442     ptLowHigh= muGetPtLowHigh(pt);
00443     pt1 = ptLowHigh.first;
00444     pt2 = ptLowHigh.second;
00445 
00446     if(pt >= pt1 && pt < pt2){    
00447       std::pair<double,double> eff1_errff1;
00448       std::pair<double,double> eff2_errff2;
00449 
00450       eff1_errff1 = muEff(pt1,eta,phi);
00451       eff2_errff2 = muEff(pt2,eta,phi); 
00452 
00453       double eff1  =  eff1_errff1.first ;
00454       double eeff1 =  eff1_errff1.second;
00455 
00456       double eff2  =  eff2_errff2.first;
00457       double eeff2 =  eff2_errff2.second;
00458 
00459       efficiency    =              eff2 *log(pt/pt1)/log(pt2/pt1);     
00460       efficiency    = efficiency + eff1 *log(pt2/pt)/log(pt2/pt1);     
00461       errEfficiency = pow((eeff2*log(pt/pt1)/log(pt2/pt1)),2); 
00462       errEfficiency = errEfficiency +pow((eeff1*log(pt2/pt)/log(pt2/pt1)),2); 
00463       errEfficiency = sqrt(errEfficiency); 
00464     }
00465 
00466     if(efficiency < 0){
00467       efficiency   =0.;
00468       errEfficiency=0.;
00469     }
00470     if(efficiency > 1){
00471       efficiency   =1.;
00472       errEfficiency=0.;
00473     }
00474 
00475     return make_pair(efficiency,errEfficiency);
00476   }

std::pair< double, double > Atlfast::MuonAcceptor::muEff ( double  ptIn,
double  etaIn,
double  phi 
)

Definition at line 479 of file MuonAcceptor.cxx.

00480   {
00481     double eta = fabs(etaIn);
00482     double pt  = fabs(ptIn);
00483  
00484     double xEtaMinOut = 0.;
00485     double xEtaMaxOut = 2.5;
00486 
00487     int nNumOut  = 0;            
00488     int nDenoOut = 0;  
00489 
00490     int nNumTotOut = 0;
00491     int nDenoTotOut= 0;
00492     //
00493     int iEtaBin = muGetEtaBin(eta);
00494     //
00495     std::vector<std::pair<int, int> > ptEtaNumeDeno;
00496     ptEtaNumeDeno.clear();
00497     ptEtaNumeDeno= muEffGetPaire(pt,xEtaMinOut,xEtaMaxOut);
00498 
00499     nNumOut=ptEtaNumeDeno[iEtaBin].first;
00500     nDenoOut=ptEtaNumeDeno[iEtaBin].second;
00501 
00502     int nBinOut = ptEtaNumeDeno.size();
00503 
00504     for(unsigned int i=0 ; i<ptEtaNumeDeno.size() ; i++){
00505       nNumTotOut = nNumTotOut+ptEtaNumeDeno[i].first;
00506       nDenoTotOut = nDenoTotOut+ptEtaNumeDeno[i].second;
00507     }
00508 
00509     double eff1    = double(nNumOut)/double(nDenoOut);                     
00510     double eeff1   = sqrt(eff1*(1.-eff1)/double(nDenoOut));                
00511     double delta   = (xEtaMaxOut-xEtaMinOut)/double(nBinOut) ;                 
00512     double xcbin1  = xEtaMinOut+delta/2.+delta*double(iEtaBin) ; 
00513     //
00514     muPhiEfficiency(pt,xcbin1,phi,eff1,eeff1);
00515     //
00516     return make_pair(eff1,eeff1);
00517   }

int Atlfast::MuonAcceptor::muGetEtaBin ( double  etaIn  ) 

Definition at line 519 of file MuonAcceptor.cxx.

00519                                            {
00520 
00521     int iEtaBin = -1;
00522     double eta = fabs(etaIn);
00523 
00524     double xMaxOut =2.5;
00525     double xMinOut = 0.;
00526     double pt=50.;
00527     //
00528     std::vector<std::pair<int, int> > ptEtaNumeDeno;
00529     ptEtaNumeDeno.clear();
00530     ptEtaNumeDeno= muEffGetPaire(pt,xMinOut,xMaxOut);
00531     //
00532     int  NbinOut   = ptEtaNumeDeno.size();
00533 
00534     double delta=double(xMaxOut-xMinOut)/double(NbinOut); 
00535     iEtaBin = int(floor(eta/delta));                         
00536     if( iEtaBin >= NbinOut) iEtaBin= NbinOut-1;  
00537 
00538     return iEtaBin;
00539   }

int Atlfast::MuonAcceptor::muGetPhiBin ( double  phi  ) 

Definition at line 617 of file MuonAcceptor.cxx.

00617                                          {
00618 
00619     double pi = 3.1415927;
00620     int iPhiBin = -1;
00621 
00622     double xMaxOut =360;
00623     double xMinOut = 0.;
00624     int  NbinOut   = 20;
00625 
00626     double phig = phi;
00627     if (phi < 0.) {
00628       phig = pi*2 + phi;
00629     }
00630     double phideg = phig*180./ pi;
00631 
00632     double delta=double(xMaxOut-xMinOut)/double(NbinOut); 
00633     iPhiBin = int(floor(phideg/delta));                         
00634     if( iPhiBin >= NbinOut) iPhiBin= NbinOut-1;
00635 
00636     return iPhiBin;
00637   }

int Atlfast::MuonAcceptor::muGetPhiBin1 ( double  phi  ) 

Definition at line 639 of file MuonAcceptor.cxx.

00639                                           {
00640 
00641     double pi = 3.1415927;
00642     int iPhiBin = -1;
00643 
00644     double xMaxOut = 22.5;
00645     double xMinOut = 0.;
00646     int  NbinOut   = 20;
00647 
00648     double phig = phi;
00649     if (phi < 0.) {
00650       phig = pi*2 + phi;
00651     }
00652     double phideg = phig*180./ pi;
00653 
00654     double delta=double(xMaxOut-xMinOut)/double(NbinOut); 
00655     iPhiBin = int(floor(phideg/delta)); 
00656     if( iPhiBin >= NbinOut) iPhiBin= NbinOut-1;               
00657  
00658     return iPhiBin;
00659   }

std::pair< double, double > Atlfast::MuonAcceptor::muPhiEff ( double  pt,
double  etamin,
double  etamax,
double  phimin,
double  phimax 
)

Definition at line 740 of file MuonAcceptor.cxx.

00740                                                                                          {
00741 
00742     int NnumTot=0 ;                                  
00743     int NdenoTot=0;  
00744 
00745     double eff =0;
00746     double eeff=0;  
00747 
00748     int ibinmin;
00749     int ibinmax;
00750 
00751     if(etamin==0. && etamax==0.15){
00752       ibinmin = muGetPhiBin1(phimin) ;
00753       ibinmax = muGetPhiBin1(phimax) ; 
00754     }
00755     else{
00756       ibinmin = muGetPhiBin(phimin) ;
00757       ibinmax = muGetPhiBin(phimax) ; 
00758     }
00759     std::vector<std::pair<int, int> > ptEtaNumeDeno;
00760     ptEtaNumeDeno.clear();
00761     ptEtaNumeDeno= muEffGetPaire(pt,etamin,etamax);
00762 
00763     for (int ibin=ibinmin;ibin<=ibinmax;ibin++){
00764       NnumTot =NnumTot+ptEtaNumeDeno[ibin].first;
00765       NdenoTot=NdenoTot+ptEtaNumeDeno[ibin].second;
00766     }
00767  
00768     eff =double(NnumTot)/double(NdenoTot); 
00769     eeff=sqrt(eff*(1.-eff)/double(NdenoTot)); 
00770     //
00771     return make_pair(eff,eeff);
00772   }

std::pair< double, double > Atlfast::MuonAcceptor::getPhiEff ( double  pt,
double  effin,
double  phig,
double  etamin,
double  etamax 
)

Definition at line 663 of file MuonAcceptor.cxx.

00664   {
00665     double eff = 1.;
00666     double eeff =0.;
00667 
00668     double phivalue[21];
00669     double efficiency[21];
00670     double errefficiency[21];
00671 
00672     for (int j = 0; j < 21; j++) {
00673       efficiency[j]    = 0.;
00674       errefficiency[j] = 0.;
00675       phivalue[j]      = 0.;
00676     }
00677     //
00678     double pi = 3.1415927;
00679     double degrad  = pi/180.;
00680     //
00681     for (int j = 0; j < 21; j++) {
00682       phivalue[j] = j*18.*degrad;
00683       if(etamin ==0. && etamax ==0.15)phivalue[j] = j*1.125*degrad;
00684     }
00685     //
00686     // Mean Efficiency in phi
00687 
00688     std::pair<double,double> eff_errff;
00689     double phiLow  = phivalue[0];
00690     double phiHigh = phivalue[20];
00691     eff_errff         = muPhiEff(pt,etamin,etamax,phiLow,phiHigh);
00692     efficiency[20]    = eff_errff.first ;
00693     errefficiency[20] = eff_errff.second;
00694 
00695     // The other bins
00696     double delta2 = 8.9*degrad;
00697     if(etamin ==0. && etamax ==0.15)delta2 = 0.55*degrad;
00698     for (int j = 0; j < 20; j++) {
00699       double phicenter = (phivalue[j+1]+phivalue[j])/2.;
00700       double phimin    =  phicenter - delta2;
00701       double phimax    =  phicenter + delta2;
00702       eff_errff = muPhiEff(pt,etamin,etamax,phimin,phimax);
00703       efficiency[j] =  eff_errff.first ;
00704       errefficiency[j] = eff_errff.second;
00705     }
00706     //
00707     double averaefficiency = effin/efficiency[20];
00708     double factornorm = averaefficiency;
00709     if(factornorm > 1. || factornorm < 0.9)factornorm =1;
00710     //
00711     if (phig >= phivalue[0] && phig <= phivalue[1]) {
00712       eff  = efficiency[0] * factornorm;
00713       eeff = errefficiency[0];
00714     }
00715     for (int j = 1; j < 20; j++) {
00716       if (phig > phivalue[j] && phig <= phivalue[j+1]) {
00717         eff = efficiency[j] * factornorm;
00718         eeff= errefficiency[j];
00719       }
00720     }
00721     if (phig > phivalue[20] && phig <= phivalue[21]) {
00722       eff = efficiency[19] * factornorm;
00723       eeff= errefficiency[19];
00724     }
00725     //
00726     if (eff < 0.) {
00727       eff = 0.;
00728       eeff = 0.;
00729     }
00730     if (eff > 1.) {
00731       eff = 1.;
00732       eeff =0.;
00733     }
00734     //
00735     return make_pair(eff,eeff);
00736   }

void Atlfast::MuonAcceptor::muPhiEfficiency ( double  pt,
double  etaIn,
double  phi,
double &  eff,
double &  eeff 
)

Definition at line 541 of file MuonAcceptor.cxx.

00541                                                                                               {
00542 
00543     double pi = 3.1415927;
00544     double phig = phi;
00545     if (phi < 0.) {
00546       phig = pi*2 + phi;
00547     }
00548 
00549     double effin  = eff;
00550     //double eeffin = eeff;
00551 
00552     double eta = fabs(etaIn);
00553 
00554     //---------------------------
00555     //Eta region 0. <|eta|<0.15
00556     //---------------------------
00557 
00558     double etamin1 = 0.;
00559     double etamax1 = 0.15;
00560     if (eta >= etamin1 && eta < etamax1) {
00561       double pion4  = pi/4.;
00562       double pion8  = pi/8.;
00563       double phig2  = fmod(phig,pion4);
00564       if (phig2 >= pion8) phig2 = pion4-phig2;
00565       //
00566       std::pair<double,double> eff_erreff;
00567       eff_erreff = getPhiEff(pt,effin,phig2,etamin1,etamax1);
00568       eff  =  eff_erreff.first ;
00569       eeff =  eff_erreff.second;
00570     }
00571 
00572     //------------------------------------------------------
00573     //Feets: Eta region 0.3 <|eta|<0.4 && 0.8 <|eta|<0.9
00574     //------------------------------------------------------
00575 
00576     double etamin2 = 0.3;
00577     double etamax2 = 0.4;
00578     double etamin22 = 0.8;
00579     double etamax22 = 0.9;
00580     if ((eta >= etamin2 && eta < etamax2) ||
00581         (eta >= etamin22 && eta < etamax22) ) {
00582       std::pair<double,double> eff_erreff;
00583       eff_erreff = getPhiEff(pt,effin,phig,etamin2,etamax2);
00584       eff  =  eff_erreff.first ;
00585       eeff =  eff_erreff.second;
00586     }
00587     //---------------------------
00588     // Eta region 1.05 <|eta|<1.15
00589     //---------------------------
00590     double etamin3 =1.05;
00591     double etamax3 =1.15;
00592     if (eta >= etamin3 && eta < etamax3) {
00593 
00594       std::pair<double,double> eff_erreff;
00595       eff_erreff = getPhiEff(pt,effin,phig,etamin3,etamax3);
00596       eff  =  eff_erreff.first ;
00597       eeff =  eff_erreff.second;
00598     }
00599     //
00600     //---------------------------
00601     // Eta region 1.20 <|eta|<1.25
00602     //---------------------------
00603     double etamin4 =1.20;
00604     double etamax4 =1.25;
00605 
00606     if (eta >= etamin4 && eta < etamax4) {
00607 
00608       std::pair<double,double> eff_erreff;
00609       eff_erreff = getPhiEff(pt,effin,phig,etamin4,etamax4);
00610       eff  =  eff_erreff.first ;
00611       eeff =  eff_erreff.second;
00612     }
00613 
00614   }


Member Data Documentation

vector<muEffdata> Atlfast::MuonAcceptor::m_muEffdataObject [private]

Definition at line 136 of file MuonAcceptor.h.

vector<double> Atlfast::MuonAcceptor::m_ptValues [private]

Definition at line 137 of file MuonAcceptor.h.

vector<double> Atlfast::MuonAcceptor::m_etaminValues [private]

Definition at line 138 of file MuonAcceptor.h.

vector<double> Atlfast::MuonAcceptor::m_etamaxValues [private]

Definition at line 139 of file MuonAcceptor.h.

MsgStream& Atlfast::MuonAcceptor::m_log [private]

Definition at line 140 of file MuonAcceptor.h.

HepRandomEngine* Atlfast::MuonAcceptor::m_randomEngine [private]

Definition at line 141 of file MuonAcceptor.h.


The documentation for this class was generated from the following files:
Generated on Mon Sep 24 14:19:41 2007 for AtlfastAlgs by  doxygen 1.5.1