AtlfastB.cxx Source File

00001 // ================================================
00002 // AtlfastB class Implementation
00003 // ================================================
00004 //
00005 // THIS TEXT TO BE REPLACED BY ATLAS STANDARD FORMAT
00006 //
00007 // Namespace ATLF::
00008 //
00009 
00010 #include <fstream>
00011 #include <cmath>
00012 
00013 #include "AtlfastAlgs/AtlfastB.h"
00014 #include "AtlfastAlgs/GlobalEventData.h"
00015 
00016 #include "AtlfastEvent/Jet.h"
00017 #include "AtlfastEvent/CollectionDefs.h"
00018 
00019 #include "AtlfastUtils/TesIO.h"
00020 #include "AtlfastUtils/HeaderPrinter.h"
00021 
00022 // Gaudi includes
00023 #include "GaudiKernel/DataSvc.h"
00024 #include "StoreGate/DataHandle.h"
00025 #include "GaudiKernel/ISvcLocator.h"
00026 #include "GaudiKernel/MsgStream.h"
00027 
00028 //Other Packages used by this class
00029 // Other Packages used by this class:-
00030 
00031 #include "CLHEP/Random/Ranlux64Engine.h"
00032 #include "CLHEP/Random/RandFlat.h"
00033 
00034 
00035 //--------------------------------
00036 // Constructors and destructors
00037 //--------------------------------
00038 
00039 namespace Atlfast {
00040   
00041   using std::abs;
00042   AtlfastB::AtlfastB( const std::string& name, ISvcLocator* pSvcLocator ) 
00043     : Algorithm( name, pSvcLocator ){
00044      //Default private parameters
00045     
00046     m_AtlfBJetSwitch    =  true;
00047     m_AtlfCalSwitch     =  true;
00048     m_AtlfTauSwitch     =  true;
00049     m_AtlfTauVetoSwitch =  false;
00050     m_AtlfTrigMuoSwitch =  false;
00051     m_atlfBNSet         =  1;
00052     m_epsitau           =  0.5;
00053     m_indtauveto        =  1;
00054     m_corrjfile         =  "AtlfastBjet.dat";
00055     m_corrcfile         =  "AtlfastBcjet.dat";
00056   
00057     
00058     // Default paths for entities in the TES
00059     m_inputLocation     = "/Event/AtlfastJets";
00060     m_outputLocation    = "/Event/AtlfastCalJets";
00061     
00062     // This is how you declare the paramemters to Gaudi so that
00063     // they can be over-written via the job options file
00064     
00065     
00066     
00067     
00068     declareProperty( "InputLocation",           m_inputLocation ) ;
00069     declareProperty( "OutputLocation",          m_outputLocation ) ;
00070     
00071     declareProperty( "AtlfBjetSwitch",           m_AtlfBJetSwitch ) ;
00072     declareProperty( "AtlfCalSwitch",            m_AtlfCalSwitch ) ;
00073     declareProperty( "AtlfTauSwitch",            m_AtlfTauSwitch );
00074     declareProperty( "AtlfTauVetoSwitch",        m_AtlfTauVetoSwitch );
00075     declareProperty( "AtlfTrigMuoSwitch",        m_AtlfTrigMuoSwitch );
00076     
00077     
00078     declareProperty( "AtlfBNSet", m_atlfBNSet ) ;
00079     declareProperty( "TauEff", m_epsitau ) ;
00080     declareProperty( "TauVetoOption", m_indtauveto ) ;
00081     
00082     declareProperty( "JetCorrFile", m_corrjfile ) ;
00083     declareProperty( "CJetCorrFile", m_corrcfile ) ;
00084 
00085   }
00086   
00087   AtlfastB::~AtlfastB() {
00088     
00089     MsgStream log( messageService(), name() ) ;
00090     log << MSG::DEBUG << "Destructor" << endreq;
00091     if( m_tesIO) {
00092       delete m_tesIO;
00093     }
00094     if( m_pRandomEngine) {
00095       delete m_pRandomEngine;
00096     }
00097     if( m_pRandFlatGenerator) {
00098       delete m_pRandFlatGenerator;
00099     }
00100 
00101   }
00102   
00103   
00104   
00105   //---------------------------------
00106   // initialise() 
00107   //---------------------------------
00108   
00109   StatusCode AtlfastB::initialize(){
00110     MsgStream log( messageService(), name() ) ;
00111     log << MSG::DEBUG << "instantiating a AtlfastB" << endreq;
00112     
00113     
00114     m_tesIO= new TesIO();
00115     
00116     //open files for jet correction factors
00117     
00118     std::ifstream inputjetfile;
00119     inputjetfile.open(m_corrjfile.c_str());
00120     
00121     if(inputjetfile){
00122       
00123       log << MSG::INFO 
00124           << "Pt jet correction factors file " 
00125           << m_corrjfile
00126           <<" open."
00127           <<endreq;
00128       int nrow;
00129       int ncolumn;
00130       inputjetfile>>nrow;
00131       inputjetfile>>ncolumn;
00132       if(nrow != 5||ncolumn != 8) {
00133         log << MSG::ERROR
00134             <<"no. of input rows,columns: "
00135             <<nrow<<" "
00136             <<ncolumn<<" "
00137             <<"expected 5, 8 "
00138             <<endreq;  
00139         return StatusCode::FAILURE ;
00140       }
00141 
00142       for(int i=0; i<nrow;i++){
00143         for(int j=0; j<ncolumn; j++){
00144           
00145           int nset;
00146           switch (j){
00147             
00148           case 0:
00149             nset=1;
00150             break;
00151           case 1:
00152             nset=2;
00153             break;
00154           case 2:
00155             nset=3;
00156             break;
00157           case 3:
00158             nset=5;
00159             break;
00160           case 4:
00161             nset=11;
00162             break;
00163           case 5:
00164             nset=12;
00165             break;
00166           case 6:
00167             nset=13;
00168             break;
00169           case 7:
00170             nset=14;
00171             break;
00172           default:
00173             log << MSG::ERROR<<"j="<<j<<" Nset= "<<nset<<"  value not correct. Nset must b one of those: 1,2,3,5,11,12,13,14 "<<endreq;
00174             return StatusCode::FAILURE;
00175             break;
00176             
00177 
00178           }
00179           inputjetfile>>m_corrj[i][nset];
00180           
00181         }
00182       }
00183       
00184     }else{
00185       log << MSG::ERROR << "Pt jet correction factors file not found" <<endreq;
00186       return StatusCode::FAILURE ;
00187     }
00188     inputjetfile.close();
00189     
00190     std::ifstream inputcjetfile;
00191     inputcjetfile.open(m_corrcfile.c_str());
00192     
00193     if(inputcjetfile){
00194       
00195       log << MSG::INFO 
00196           << "Pt c jet correction factors file " 
00197           << m_corrcfile
00198           <<" open."
00199           <<endreq;
00200       int nrow;
00201       int ncolumn;
00202       inputcjetfile>>nrow;
00203       inputcjetfile>>ncolumn;
00204        if(nrow != 5||ncolumn != 8) {
00205         log << MSG::ERROR
00206             <<"no. of input rows,columns: "
00207             <<nrow<<" "
00208             <<ncolumn<<" "
00209             <<"expected 5, 8 "
00210             <<endreq;  
00211         return StatusCode::FAILURE ;
00212       }
00213       for(int i=0; i<nrow;i++){
00214         for(int j=0; j<ncolumn; j++){
00215           
00216           int nset;
00217           switch (j){
00218             
00219           case 0:
00220             nset=1;
00221             break;
00222           case 1:
00223             nset=2;
00224             break;
00225           case 2:
00226             nset=3;
00227             break;
00228           case 3:
00229             nset=5;
00230             break;
00231           case 4:
00232             nset=11;
00233             break;
00234           case 5:
00235             nset=12;
00236             break;
00237           case 6:
00238             nset=13;
00239             break;
00240           case 7:
00241             nset=14;
00242             break;
00243           default:
00244             log << MSG::ERROR<<"Nset value not correct. Nset must b one of those: 1,2,3,5,11,12,13,14 "<<endreq;
00245             return StatusCode::FAILURE;
00246             break;
00247           }
00248           inputcjetfile>>m_corrc[i][nset];
00249         }
00250       }
00251       
00252     }else{
00253       log << MSG::ERROR 
00254           << "Pt c jet correction factors file not found" 
00255           <<endreq;
00256       return StatusCode::FAILURE ;
00257     }
00258     inputcjetfile.close();
00259 
00260     //Flat random number generator
00261     //get the Global Event Data using singleton pattern
00262     
00263     GlobalEventData* ged = GlobalEventData::Instance();
00264     int randSeed = ged->randSeed() ;
00265     
00266     m_pRandomEngine = new Ranlux64Engine(randSeed);
00267     m_pRandFlatGenerator=new RandFlat(*m_pRandomEngine);
00268     
00269     
00270     HeaderPrinter hp("AtlfastB:", log);
00271     
00272     hp.add("TES Locations:              ");
00273     hp.add(" Jets from                     ", m_inputLocation); 
00274     hp.add(" Calibrated and tagged Jets to ", m_outputLocation);
00275     hp.add("AtlfBje switch:                ", m_AtlfBJetSwitch);
00276     hp.add("AtlfBje NSet:                  ", m_atlfBNSet);
00277     hp.add("AtlfCal switch:                ", m_AtlfCalSwitch);
00278     hp.add("AtlfTau switch:                ", m_AtlfTauSwitch);    
00279     hp.add("AtlfTauVeto switch:            ", m_AtlfTauVetoSwitch);
00280     hp.add("AtlfTrigMuo switch:            ", m_AtlfTrigMuoSwitch); 
00281     hp.print();
00282     
00283     return StatusCode::SUCCESS ;
00284   }
00285   
00286   //---------------------------------
00287   // finalise() 
00288   //---------------------------------
00289   
00290   StatusCode AtlfastB::finalize(){
00291     
00292     MsgStream log( messageService(), name() ) ;
00293     log << MSG::INFO << "finalizing" << endreq;    
00294     return StatusCode::SUCCESS ;
00295   }
00296   
00297   
00298   //----------------------------------------------
00299   // execute() method called once per event
00300   //----------------------------------------------
00301   
00302   StatusCode AtlfastB::execute(){
00303     
00304     StatusCode sc;
00305     MsgStream log( messageService(), name() ) ;  
00306     
00307     log << MSG::DEBUG<<"In execute"<<endreq;
00308 
00309     
00310     std::vector<Jet*>  originalJets;
00311     if( ! m_tesIO->copy<JetCollection>( originalJets, m_inputLocation ) ) {
00312       log << MSG::INFO << "No Jets in TES " << endreq;
00313       return StatusCode::FAILURE;
00314       
00315     }  
00316     
00317     
00318     
00319     log << MSG::DEBUG<<"Found Jets in TES"<<endreq;
00320     
00321     
00322     m_Jets.clear();
00323     for(m_it=originalJets.begin(); m_it<originalJets.end(); ++m_it){
00324       m_Jets.push_back(new Jet(**m_it));
00325       log << MSG::DEBUG<<"Jet "<<**m_it<<endreq;
00326     }
00327     
00328     
00329     if(m_AtlfBJetSwitch){
00330       //Apply efficiency for b-jet identification
00331       
00332 
00333       
00334       sc=atlfBje(m_atlfBNSet);
00335       
00336       if(!sc){
00337         
00338         log << MSG::ERROR<<"Error in atlfBje"<<endreq;
00339         return StatusCode::FAILURE;
00340       }
00341     }
00342    
00343 
00344 
00345 
00346     if(m_AtlfTauSwitch){
00347       
00348       //Apply efficiency for tau-jet identification
00349       
00350       sc=atlfTau(m_epsitau);
00351       
00352       if(!sc){
00353         
00354         log << MSG::ERROR<<"Error in atlfTau"<<endreq;
00355         return StatusCode::FAILURE;
00356       }
00357     }
00358     
00359     if(m_AtlfTauVetoSwitch){
00360       
00361       //Apply veto for tau-jet
00362       
00363       sc=atlfTauVeto(m_indtauveto);
00364       
00365       if(!sc){
00366         
00367         log << MSG::ERROR<<"Error in atlfTauVeto"<<endreq;
00368         return StatusCode::FAILURE;
00369       }
00370     }
00371     
00372  
00373     if(m_AtlfTrigMuoSwitch){
00374       
00375       //Apply efficiency for muon trigger. Not implemented yet!
00376       
00377       sc=atlfTrigMuo();
00378       
00379       if(!sc){
00380         
00381         log << MSG::ERROR<<"Error in atlfTrigMuo"<<endreq;
00382         return StatusCode::FAILURE;
00383       }
00384     }   
00385     
00386     
00387     if(m_AtlfCalSwitch){
00388       
00389       //Recalibrate jet energies
00390       
00391       sc=atlfCal();
00392       
00393       if(!sc){
00394         
00395         log << MSG::ERROR<<"Error in atlfCal"<<endreq;
00396         return StatusCode::FAILURE;
00397       }
00398     }
00399     
00400     
00401    
00402    
00403     
00404     
00405     
00406     //......................................
00407     // Fill jets that have been successfully tagged and calibrated in JetCollection 
00408     //and register them in the transient event store. 
00409     
00410     
00411     JetCollection* myCalJets=new JetCollection;
00412     
00413     for(m_it=m_Jets.begin();m_it<m_Jets.end();++m_it){
00414       myCalJets->push_back(*m_it);
00415     }
00416 
00417     Atlfast::JetCollection::const_iterator jetcolit = myCalJets->begin();
00418     m_it = originalJets.begin();
00419     
00420 
00421     for( ;jetcolit != myCalJets->end(); ++jetcolit ){
00422 
00423        log << MSG::DEBUG
00424         <<"original jets: PDGID= "<<(*m_it)->pdg_id()
00425           <<" B tag= "    <<(*m_it)->isBTag()
00426           <<" C tag= "    <<(*m_it)->isCTag()
00427           <<" Light tag= "<<(*m_it)->isLightTag()
00428           <<" Tau tag= "  <<(*m_it)->isTauTag()<<endreq;
00429       log << MSG::DEBUG
00430           <<"Jets After AtlfastB: PDGID= "<<(*jetcolit)->pdg_id()
00431           <<" B tag= "    <<(*jetcolit)->isBTag()
00432           <<" C tag= "    <<(*jetcolit)->isCTag()
00433           <<" Light tag= "<<(*jetcolit)->isLightTag()
00434           <<" Tau tag= "  <<(*jetcolit)->isTauTag()<<endreq;
00435       
00436       //   log << MSG::DEBUG
00437       //          <<"calibrated: PDGID= "<<(*jetcolit)->pdg_id()<<" pt= "<<(*jetcolit)->pT()<<endreq;
00438 
00439       ++m_it;
00440     }
00441     TesIoStat stat;
00442     
00443     stat = m_tesIO->store( myCalJets, m_outputLocation ) ;
00444     if(!stat){
00445       log<<MSG::ERROR<<"Could not store jets"<<endreq;
00446       return stat;
00447     }
00448     
00449     return StatusCode::SUCCESS; 
00450   }
00451   
00452   //-------------------------------
00453   // helper methods implementation
00454   //-------------------------------
00455   
00456   StatusCode AtlfastB::atlfBje(int nset){
00457     
00458     MsgStream log( messageService(), name() ) ; 
00459     
00460 
00461     switch(nset){
00462       
00463     case 1:
00464       m_epsib=0.5;
00465       m_epsic=1./10.9;
00466       m_epsij=1./231.;
00467       break;
00468       
00469     case 2:
00470       m_epsib=0.60;
00471       m_epsic=1./6.7;
00472       m_epsij=1./93.;
00473       break;
00474       
00475     case 3:
00476       m_epsib=0.70;
00477       m_epsic=1./4.3;
00478       m_epsij=1./34.1;
00479       break;
00480       
00481     case 5:
00482       m_epsib=0.60;
00483       m_epsic=1./10.;
00484       m_epsij=1./100.;
00485       break;
00486       
00487     case 11:
00488       m_epsib=0.33;
00489       m_epsic=1./22.9;
00490       m_epsij=1./1381.;
00491       break;
00492       
00493     case 12:
00494       m_epsib=0.43;
00495       m_epsic=1./10.8;
00496       m_epsij=1./219.;
00497       break;
00498       
00499     case 13:
00500       m_epsib=0.53;
00501       m_epsic=1./6.7;
00502       m_epsij=1./91.;
00503       break;
00504       
00505     case 14:
00506       m_epsib=0.624;
00507       m_epsic=1./6.7;
00508       m_epsij=1./91.;
00509       break;
00510       
00511     default:
00512       log << MSG::ERROR<<"No b jet tagging efficiencies applied"<<endreq;
00513       return StatusCode::FAILURE;
00514       break;
00515     }
00516     
00517     
00518     //Apply randomized B tagging
00519     
00520     log<<MSG::DEBUG<<"Applying randomized B tagging...."<<endreq;
00521     
00522     for(m_it=m_Jets.begin();m_it<m_Jets.end();++m_it){
00523       
00524       double randnum;
00525       double corpt;    
00526       double corr;
00527       
00528       double ptjet=(*m_it)->pT();
00529       double etajet=(*m_it)->eta();
00530       int pdgid=(*m_it)->pdg_id();
00531       
00532       //apply b tagging as for recalibrated jets
00533       
00534       
00535       corr=fitcoreb(ptjet);
00536 
00537       ptjet=corr*ptjet;
00538 
00539       randnum=m_pRandFlatGenerator->shoot(m_pRandomEngine);
00540       
00541       corpt=0.;
00542       
00543       if(abs(pdgid)==5){
00544         if((abs(etajet)<=2.5)&&(randnum<m_epsib))
00545           {
00546             
00547             (*m_it)->setBTag();
00548           } else{
00549 
00550             (*m_it)->setLightTag();
00551 
00552 
00553           }
00554       }else if(abs(pdgid)==4){
00555 
00556 
00557         if(ptjet<=30.) corpt=m_corrc[0][m_atlfBNSet];
00558         if((ptjet>30.)&&(ptjet<=45.)) corpt=m_corrc[1][m_atlfBNSet];
00559         if((ptjet>45.)&&(ptjet<=60.)) corpt=m_corrc[2][m_atlfBNSet];    
00560         if((ptjet>60.)&&(ptjet<=100.)) corpt=m_corrc[3][m_atlfBNSet];
00561         if(ptjet>100.) corpt=m_corrc[4][m_atlfBNSet];
00562         if((abs(etajet)<=2.5)&&(randnum<(m_epsic/corpt))) {
00563           
00564           (*m_it)->setBTag();
00565           //log << MSG::DEBUG << "misstaged u jet. ptjet=" <<ptjet<< endreq;
00566 
00567         }else{
00568           (*m_it)->setLightTag();
00569 
00570         }
00571     
00572 
00573     }else if(((abs(pdgid)!=4)&&(abs(pdgid)!=5))||abs(pdgid)==15){
00574 
00575 
00576         if(ptjet<=30.) corpt=m_corrj[0][m_atlfBNSet];
00577         if((ptjet>30.)&&(ptjet<=45.)) corpt=m_corrj[1][m_atlfBNSet];
00578         if((ptjet>45.)&&(ptjet<=60.)) corpt=m_corrj[2][m_atlfBNSet];    
00579         if((ptjet>60.)&&(ptjet<=100.)) corpt=m_corrj[3][m_atlfBNSet];
00580         if(ptjet>100.) corpt=m_corrj[4][m_atlfBNSet];
00581         if((abs(etajet)<=2.5)&&(randnum<(m_epsij/corpt))) {
00582           
00583           (*m_it)->setBTag();
00584           //      log << MSG::DEBUG << "misstaged tau jet. ptjet=" <<ptjet<< endreq;
00585         }
00586       }
00587       
00588     } 
00589     return StatusCode::SUCCESS;
00590     
00591   }
00592   
00593   
00594   StatusCode AtlfastB::atlfCal(){
00595     
00596     //It recalibrates the energy of all jets using parameterisation 
00597     //obtained with 'Z+jet' equivalent method
00598     
00599     MsgStream log( messageService(), name() ) ; 
00600     
00601     
00602     
00603     for(m_it=m_Jets.begin();m_it<m_Jets.end();++m_it){
00604       
00605       
00606       double ptjet=(*m_it)->pT();
00607       
00608       
00609       double corr;
00610       
00611       if((*m_it)->isBTag()){
00612         
00613         corr=fitcoreb(ptjet);
00614       }else{
00615         corr=fitcoreu(ptjet);
00616       }
00617       
00618       (*m_it)->setMomentum(corr*((*m_it)->momentum()));
00619       
00620       
00621       
00622     }
00623     return StatusCode::SUCCESS;
00624   }
00625   
00626   
00627   StatusCode AtlfastB::atlfTau(double epsitau){
00628     
00629     MsgStream log( messageService(), name() ) ; 
00630     StatusCode sc;
00631     
00632     for(m_it=m_Jets.begin();m_it<m_Jets.end();++m_it){
00633       
00634       double randnum;
00635       
00636       
00637       
00638       
00639       double etajet=(*m_it)->eta();      
00640       double ptjet=(*m_it)->pT();
00641       int pdgid=(*m_it)->pdg_id();
00642       
00643       randnum=m_pRandFlatGenerator->shoot(m_pRandomEngine);
00644       
00645       int iflag=0;
00646 
00647       
00648       
00649       if(abs(pdgid)==15){
00650         
00651         if(randnum<epsitau){
00652 
00653 
00654           (*m_it)->setTauTag(15);
00655         }else{
00656           (*m_it)->setLightTag();
00657         }
00658       }else if(abs(etajet)<=2.5){
00659         
00660         double rjet;
00661         //Tau tag it computes the rjet tau rejection value
00662 
00663         sc=tautag(ptjet,etajet,epsitau*100,rjet,iflag);
00664                   
00665           if((iflag==0)&&(randnum<(1./rjet))){
00666 
00667             (*m_it)->setTauTag(15);
00668           }
00669           
00670       }
00671       
00672     }
00673     return StatusCode::SUCCESS;
00674     
00675   }
00676   
00677   
00678   StatusCode AtlfastB::tautag(double pt, 
00679                               double eta, 
00680                               double efftau, 
00681                               double &rjet, int &iflag){
00682     
00683     //It computes the corresponding jet rejection rjet
00684     //inputs: pt and eta of cluster, efftau-> eff(in %)  for tau identification
00685     
00686     MsgStream log( messageService(), name() ) ; 
00687 
00688     iflag=0;    
00689     
00690     if((pt<15.)||(pt>150.)){
00691       
00692       
00693       log << MSG::DEBUG<<"Tau pt value out of range (15<pt<150): pt= "<<pt<<" No tau tagging efficiencies will be applied"<<endreq;
00694 
00695       iflag=1;
00696       
00697     }
00698     if(abs(eta)>2.5){
00699       
00700       log << MSG::DEBUG<<"Tau eta value out of range (eta<2.5): eta= "<<eta<<" No tau tagging efficiencies will be applied"<<endreq;
00701       
00702 
00703       iflag=1;
00704       
00705     }
00706     
00707 
00708     //eta dependence
00709     double eff=efftau;
00710     double coefe1;
00711     double coefe3;
00712     double coef1;
00713     double coef2;
00714     
00715     if(abs(eta)<0.7){
00716       coefe1=1.35-0.0035*efftau;
00717       eff=efftau/coefe1;
00718     }
00719     if(abs(eta)>1.5){
00720       coefe3=0.70+0.0030*efftau;
00721       eff=efftau/coefe3;
00722     }
00723     if(eff>100){
00724       
00725       
00726       log << MSG::WARNING<<"Tau efficiency value > 100 ! eff= "<<eff<<endreq;
00727       iflag=1;
00728     }
00729     
00730     //pt dependence
00731     
00732     coef1=0.027+0.00024*pt;
00733     coef2=2.28+0.027*pt;
00734     rjet=pow(10,(-coef1*eff+coef2));
00735     
00736   
00737     
00738     //     log << MSG::DEBUG<<"Tau rejection value computed in tautag= "<<rjet<<endreq;
00739     
00740     
00741     return StatusCode::SUCCESS;
00742     
00743   }
00744   
00745   
00746   StatusCode  AtlfastB::atlfTauVeto(int ind){
00747     
00748     MsgStream log( messageService(), name() ) ; 
00749     
00750     
00751     StatusCode sc;
00752     
00753     for(m_it=m_Jets.begin();m_it<m_Jets.end();++m_it){
00754       
00755       double randnum;    
00756       double efftau;
00757       double effjet;
00758       
00759       double ptjet=(*m_it)->pT();
00760       int pdgid=(*m_it)->pdg_id();
00761 
00762       sc=tauveto(ptjet,ind,efftau,effjet);
00763       
00764       if(!sc){
00765         log << MSG::ERROR<<"Cannot compute tau jet rejection"<<endreq;
00766         return StatusCode::FAILURE;
00767       }
00768       
00769       randnum=m_pRandFlatGenerator->shoot(m_pRandomEngine);
00770       
00771 
00772       if(abs(pdgid)==15) {
00773         if(randnum<efftau/100.){
00774           (*m_it)->setLightTag();
00775         }
00776       }else{
00777         if(randnum>effjet/100.){
00778 
00779           (*m_it)->setTauTag(15);
00780         }
00781       }
00782     }
00783     
00784     return StatusCode::SUCCESS;
00785   }
00786   
00787   
00788   
00789   StatusCode AtlfastB::tauveto(double pt, int ind, double &efftau, double &effjet){ 
00790     MsgStream log( messageService(), name() ) ; 
00791     
00792     
00793 
00794     //input: pt of cluster
00795     //ind=1->fix efftau to 5% and gives effjet
00796     //ind=2->fix effjet to 90% and gives efftau
00797     //output: efftau tau efficiency in % or effjet jets efficiency in %
00798     
00799     
00800     //using calo criteria + tracks
00801     
00802     
00803     
00804 
00805     if(ind==1){
00806      
00807       efftau=5.;
00808       effjet=36.+1.5*pt-0.01*pow(pt,2);
00809       if(pt>60) effjet=90.;
00810     }
00811     else if(ind==2){
00812       effjet=90.;
00813       efftau=27.-0.35*pt;
00814       if(pt>60.)efftau=5.;
00815 
00816     }else{
00817       
00818       log << MSG::ERROR<<"Wrong ind values for tauveto method. ind= "<<ind<<"ind=1 or ind=2"<<endreq;
00819       
00820       return StatusCode::FAILURE;
00821     }
00822     
00823     
00824     
00825     return StatusCode::SUCCESS;
00826   }
00827   
00828   StatusCode AtlfastB::atlfTrigMuo(){
00829     MsgStream log( messageService(), name() ) ; 
00830     
00831     log << MSG::INFO<<"atlfTrigMuo is not yet implemented!!!!!!"<<endreq;
00832     
00833     
00834     return StatusCode::SUCCESS;
00835   }
00836   double AtlfastB::fitcoreb(double pt){
00837     
00838     
00839     MsgStream log( messageService(), name() ) ; 
00840     
00841     double core;
00842     double a0,a1,a2,a3,a4,a5;
00843     double xc;
00844     
00845     
00846     
00847     if(pt<10.){
00848       
00849       core=0.;
00850     }else if(pt<55.){
00851       
00852       a0=1.2715;
00853       a1=0.12241;
00854       a2=-0.10480e-01;
00855       a3=0.33310e-03;
00856       a4=-0.47454e-05;
00857       a5=0.25436e-07;
00858       core=a0+a1*pt+a2*pow(pt,2)+a3*pow(pt,3)+a4*pow(pt,4)+a5*pow(pt,5);
00859       core=core*1.006;
00860     }else if(pt<200.){
00861       
00862       a0=1.18;
00863       a1=-0.16672e-02;
00864       a2=0.44414e-05;
00865       core=a0+a1*pt+a2*pow(pt,2);
00866       
00867     }else{
00868       
00869       xc=200.;
00870       a0=1.18;
00871       a1=-0.16672e-02;
00872       a2=0.44414e-05;
00873       core=a0+a1*xc+a2*pow(xc,2);
00874     }
00875     
00876     
00877     return core;
00878     
00879   }
00880   
00881   double AtlfastB::fitcoreu(double pt){
00882     
00883     MsgStream log( messageService(), name() ) ;     
00884     
00885     double core;
00886     double a0,a1,a2,a3,a4,a5;
00887     double xc;
00888     
00889     
00890     
00891     if(pt<10.){
00892       
00893       core=0.;
00894       
00895       
00896       
00897       
00898       
00899     }else if((pt<45.)&&(pt>10.)){
00900       
00901       a0=1.5085;
00902       a1=0.31468e-01;
00903       a2=-0.36973e-02;
00904       a3=0.11220e-03;
00905       a4=-0.13921e-05;
00906       a5=0.61538e-08;
00907       core=a0+a1*pt+a2*pow(pt,2)+a3*pow(pt,3)+a4*pow(pt,4)+a5*pow(pt,5);
00908       
00909     }else if(pt<200.){
00910       
00911       a0=1.18;
00912       a1=-0.16672e-02;
00913       a2=0.44414e-05;
00914       core=a0+a1*pt+a2*pow(pt,2);
00915       core=core/1.025;
00916     }else{
00917       
00918       xc=200.;
00919       a0=1.18;
00920       a1=-0.16672e-02;
00921       a2=0.44414e-05;
00922       core=a0+a1*xc+a2*pow(xc,2);
00923       core=core/1.025;
00924     }
00925     
00926     
00927     return core;
00928     
00929   }
00930   
00931   
00932 } //end namespace bracket