// ElectronReconstructor.cxx
 //
 // Implementation of the Electron Reconstructor class
 //
 // Only the reconstruct() method is overridden
 //
 //
 // Authors: H.T. Phillips, P.Clarke, E.Richter-Was, P.Sherwood, R.Steward
 // Update:  Sabine Elles, Simon Dean
 //
 // Smearer - Reconstructor Conversion: Alexander Richards (UCL)
 //
 #include "AtlfastAlgs/ElectronReconstructor.h"
 #include "AtlfastAlgs/CorrectionFactorElectron.h"
 
 #include <cmath>
 
 #include "CLHEP/Vector/LorentzVector.h"
 
 
 #include "CLHEP/Random/JamesRandom.h"
 #include "CLHEP/Random/RandGauss.h"
 #include "CLHEP/Units/SystemOfUnits.h"
 
 /** @brief Electron implementation of the IDefaultReconstructor interface.
 *
 * Privately inherits from the DefaultReconstructor class to use its implementation of
 * random number services etc.
 * ElectronReconstructor can run with TDR parameterisation (m_smearParamSchema==1)
 * or with CSC parameterisation (m_smearParamSchema==2)
 */
 
 namespace Atlfast {
 
   ReconstructedParticle ElectronReconstructor::reconstruct( const ReconstructedParticle &particle )const{
     double rpilup = 0.0;
     double aa, bb, sigph1, sigph, sigpu;
     //Make sure momentum, aenergy units are Gev
     double ene     = particle.e()/GeV ;
     double sqrtene = sqrt(ene);
     double pt      = particle.pt()/GeV ;
     double vEta    = fabs(particle.eta());
     HepLorentzVector bvec(0,0,0,0);
     
     if(m_smearParamSchema==1){ // TDR-style smearing
       
       // do the smearing, copied verbatim (except for ROOT dependencies)
       // from electronmaker codein Atlfast++
       //
       // This code has not otherwise been altered which is why it is all very
       // procedural
       //
       // Parametrisation was by L.Poggioli
       
       while (1) {
         aa=m_randGauss->fire();
         //bb=m_randGauss-->fire();
         sigph1 = aa*m_smearParams[0]/sqrtene;
         if (1.0 + sigph1 <= 0.0) continue;
         sigph = sigph1;
         if (vEta < 1.4) {
           while (1) {
             aa=m_randGauss->fire();
             bb=m_randGauss->fire();
             sigph1 = aa*m_smearParams[1]/pt + bb*m_smearParams[2];
             if (1.0+sigph1 > 0) break;
           }
         } else {
           while (1) {
             aa=m_randGauss->fire();
             bb=m_randGauss->fire();
             sigph1 = aa*m_smearParams[3]*((m_smearParams[4]-vEta)+m_smearParams[5])/ene + bb*m_smearParams[6];
             //    sigph1 = aa*0.306*((2.4-vEta)+0.228)/ene + bb*0.007;
             if (1.0+sigph1 > 0) break;
           }
         }
         sigph += sigph1;
         if (m_lumi == 2) {
           if (vEta < 0.6) rpilup = 0.32;
           if (vEta > 0.6 && vEta < 1.4) rpilup = 0.295;
           if (vEta > 1.4) rpilup = 0.27;
           while (1) {
             aa=m_randGauss->fire();
             //bb=m_randGauss-->fire();
             sigpu = aa*rpilup/pt;
             if (1.0+sigpu > 0) break;
           }
           sigph += sigpu;
         }
         if (1.0 + sigph > 0) break;
       }
       
       //now sigph is the electron "sigma" and we do the following a la Atlfast++ electron maker:
       bvec.set(particle.px()*(1.0+sigph), particle.py()*(1.0+sigph), particle.pz()*(1.0+sigph), particle.e() *(1.0+sigph));
       
     } else if (m_smearParamSchema==2){ // Smearing based on CSC sample tuning - data in CorrectionFactorElectron.h
 
         
     
       double vE = particle.e();
 
       double vLimEta=EtaElectron[NbEtaElectron-1]+(EtaElectron[NbEtaElectron-1]-EtaElectron[NbEtaElectron-2])*0.5;
     
       if( vEta>vLimEta )
         {
           m_log<<MSG::DEBUG<<"No calibration : "<<vEta<<endreq;
           return particle;
         }
     
       int iEnergy=0;
       m_log<<MSG::DEBUG<<vE<<" > "<<vEnergiesElectron[iEnergy]<<" ";
       while(iEnergy<NbEnergiesElectron&&vE>vEnergiesElectron[iEnergy]*1000)
         {
           m_log<<MSG::DEBUG<<vEnergiesElectron[iEnergy]<<" ";
           iEnergy++;
         }
       m_log<<MSG::DEBUG<<endreq;
       if(iEnergy==0)iEnergy++;
       double ratioEnergy=(vE*0.001-vEnergiesElectron[iEnergy-1])/
              (vEnergiesElectron[iEnergy]-vEnergiesElectron[iEnergy-1]);
 
       int iEta=0;
       m_log<<MSG::DEBUG<<vEta<<" > "<<EtaElectron[iEta]<<" ";
       while(iEta<NbEtaElectron&&vEta>EtaElectron[iEta])
         {
           m_log<<MSG::DEBUG<<EtaElectron[iEta]<<" ";
           iEta++;
         }
       m_log<<MSG::DEBUG<<endreq;
 
       if(iEta==0)iEta++;
       double ratioEta=(vEta-EtaElectron[iEta-1])/(EtaElectron[iEta]-EtaElectron[iEta-1]);
 
       int iPos;
       double sigma,sigma1,sigma2;
 
       iPos=iEnergy-1;
       sigma1=CorrFactorElectronSigma[iPos][iEta-1]+ratioEta*
   (CorrFactorElectronSigma[iPos][iEta]-CorrFactorElectronSigma[iPos][iEta-1]);
       iPos=iEnergy;
       sigma2=CorrFactorElectronSigma[iPos][iEta-1]+ratioEta*
   (CorrFactorElectronSigma[iPos][iEta]-CorrFactorElectronSigma[iPos][iEta-1]);
 
       sigma=sigma1+ratioEnergy*(sigma2-sigma1);
 
       m_log<<MSG::DEBUG<<" -> interpolation sigma "<<sigma1<<" "<<sigma2<<" "<<sigma<<endreq;
 
       double alpha=1.0;
       double vrandom=alpha;
 
       int iCmpt=int(m_randFlat->fire()*20);
       for(int i=0; i<iCmpt; i++)
         vrandom=m_randGauss->fire(1.0,sigma);
     
       m_log<<MSG::DEBUG<<" -> Energie : "<<vEnergiesElectron[iEnergy]<<
         "   Eta : "<<vEnergiesElectron[iEnergy-1]<<" < "<<vE<<" < "<<
          vEnergiesElectron[iEnergy]<<endreq;
       m_log<<MSG::DEBUG<<" -> Energie : "<<vEnergiesElectron[iEnergy]<<
         "   Eta : "<<EtaElectron[iEta-1]<<" < "<<vEta<<" < "<< 
          EtaElectron[iEta]<<endreq;
       m_log<<MSG::DEBUG<<"       calib mean-sigma  "<<alpha<<" "<<sigma<<" -> "<<vrandom<<endreq;
 
       bvec.set(particle.px()*vrandom, particle.py()*vrandom, particle.pz()*vrandom, vE*vrandom);
     
       m_log<<MSG::DEBUG<<"-> HepLorentzVector calibration "<<vEta<<" "<<particle<<"  ->  "<<bvec<<endreq;
 
     }
 
     ReconstructedParticle return_particle(particle);
     return_particle.set_momentum(bvec);
     return return_particle;
   } //end of smear method
   
   //cct: implement the setSmearParameters method
   int ElectronReconstructor::setSmearParameters (const std::vector<double>& smearValues){
     m_smearParams=smearValues;
     return 0;
   }
   
   //cct: implement the setSmearParamSchema method
   int ElectronReconstructor::setSmearParamSchema ( const int smearSchema){
     m_smearParamSchema=smearSchema;
     return 0;
   }
   
 } // end of namespace bracket
 

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.