AtlfastAlgs: Atlfast::PhotonSmearer Class Reference

Atlfast::PhotonSmearer::PhotonSmearer	(	const int	aseed,
		const int	lumi,
		MsgStream &	log
	)			`[inline]`

Definition at line 44 of file PhotonSmearer.h.

00044                                                                      : 
00045         ISmearer(), DefaultSmearer(aseed), m_lumi(lumi), m_log(log){}

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

Definition at line 46 of file PhotonSmearer.h.

00046 { }

HepLorentzVector Atlfast::PhotonSmearer::smear ( const HepMC::GenParticle & ) [virtual]

Smear methods

Reimplemented from Atlfast::DefaultSmearer.

Definition at line 26 of file PhotonSmearer.cxx.

00026                                                                        {
00027     return smear(particle.momentum());
00028   }

HepLorentzVector Atlfast::PhotonSmearer::smear ( const HepLorentzVector & avec ) [virtual]

Smear method for HepLorentzVector

Reimplemented from Atlfast::DefaultSmearer.

Definition at line 30 of file PhotonSmearer.cxx.

00030                                                                      {
00031     
00032     double rpilup = 0.0;
00033     double aa, bb, sigph1, sigph, sigpu;
00034     double vEta  = fabs(avec.pseudoRapidity());
00035     
00036     //make sure e, momentum in GeV
00037     double pt = avec.perp()/GeV;
00038     double ene = avec.e()/GeV;
00039     double sqrtene = sqrt(ene);
00040     
00041     // The HLV to be returned
00042     HepLorentzVector return_vec(0.,0.,0.,0.);
00043     if (!sqrtene || !ene || !pt) return return_vec;
00044 
00045     if (m_smearParamSchema==1){ // TDR-style smearing
00046       
00047       // do the smearing, copied verbatim (except for ROOT dependencies)
00048       // from photonmaker codein Atlfast++
00049       // The code is very procedural (even more so than Atlfast++!) and should be tidied later!
00050       //
00051       // Parametrisation was by L.Poggioli and F. Gianotti
00052       
00053       // do smearing of theta first
00054       HepLorentzVector bvec(avec);
00055       
00056       //this is equivalent to RESTHE in FORTRAN 
00057       double sigph2 = 0.0;
00058       while (1) {
00059         aa=randGauss()->fire();
00060         if (vEta < 0.8) sigph2                  = aa*m_smearParams[0]/sqrtene;
00061         //if (vEta < 0.8) sigph2                  = aa*0.065/sqrtene;
00062         if (vEta >= 0.8 && vEta < 1.4) sigph2 = aa*m_smearParams[1]/sqrtene;  
00063         //if (vEta >= 0.8 && vEta < 1.4) sigph2 = aa*0.050/sqrtene;
00064         if (vEta >= 1.4 && vEta < 2.5) sigph2 = aa*m_smearParams[2]/sqrtene; 
00065         //if (vEta >= 1.4 && vEta < 2.5) sigph2 = aa*0.040/sqrtene;
00066         if (vEta >= 2.5)                 sigph2 = 0;
00067         if (fabs(bvec.theta()) + sigph2 <= M_PI) break;
00068       }
00069       
00070       // sigph2 is now the offset for theta...
00071       bvec.setPz(bvec.e()*cos(bvec.theta()+sigph2));
00072       
00073       vEta  = fabs(bvec.pseudoRapidity()); //27/02/03 JPC possible bug fix
00074       
00075       //this is RESPHO in FORTRAN
00076       // now do the energy resolution note: use original energy, pt and phi, but new vEta
00077       //
00078       while (1) {
00079         while (1){
00080           aa=randGauss()->fire();
00081           sigph1 = aa*m_smearParams[3]/sqrtene;
00082           //sigph1 = aa*0.10/sqrtene;
00083           if (1.0 + sigph1 > 0.0) break;
00084         }
00085         sigph = sigph1;
00086         if (vEta < 1.4) {
00087           while (1) {
00088             aa=randGauss()->fire();
00089             bb=randGauss()->fire();
00090             sigph1 = aa*m_smearParams[4]/pt + bb*m_smearParams[5];
00091             //sigph1 = aa*0.245/pt + bb*0.007;
00092             if (1.0+sigph1 > 0) break;
00093           }
00094         } else {
00095           while (1) {
00096             aa=randGauss()->fire();
00097             bb=randGauss()->fire();
00098             sigph1 = aa*m_smearParams[6]*((m_smearParams[7]-vEta)+m_smearParams[8])/ene + bb*m_smearParams[9];
00099             // sigph1 = aa*0.306*((2.4-vEta)+0.228)/ene + bb*0.007;
00100             if (1.0+sigph1 > 0) break;
00101           }
00102         }
00103         sigph += sigph1;
00104         if (m_lumi == 2) {
00105           if (vEta < 0.6) rpilup = 0.32;
00106           if (vEta > 0.6 && vEta < 1.4) rpilup = 0.295;
00107           if (vEta > 1.4) rpilup = 0.27;
00108           while (1) {
00109             aa=randGauss()->fire();
00110             sigpu = aa*rpilup/pt;
00111             if (1+sigpu > 0) break;
00112           }
00113           sigph += sigpu;
00114         }
00115         
00116         if (1.0 + sigph > 0) break;
00117       }
00118       
00119       
00120       // Now sigph is the photon "sigma" and we do the following a la Atlfast++ photon maker:
00121       // this seems to be OK for energy resolution. 
00122       //
00123       // It looks to me as though the functional form above for photons is actually the same as
00124       // electrons, just that one magic number is different. Suggest we revisit this and
00125       // exploit commonality.
00126       //
00127       
00128       HepLorentzVector cvec;
00129       cvec.setPx(avec.px()*(1.0+sigph));
00130       cvec.setPy(avec.py()*(1.0+sigph));
00131       cvec.setPz(bvec.pz()*(1.0+sigph));//using recalculated z momentum
00132       
00133       // nb from FORTRAN code, we go back to the original photon energy here
00134       cvec.setE( avec.e() *(1.0+sigph));
00135       //return cvec;
00136       //now do what the fortran does to get massless particle! YUCK!!!
00137       
00138       //if(std::abs (cvec.pz()/cvec.e()) > 1.){
00139       //        theta = acos(cvec.pz()/sqrt(cvec.px()*cvec.px()+cvec.py()*cvec.py()+cvec.pz()*cvec.pz()));
00140       //} else if (cvec.pz()>0) theta = 0.;
00141       //else if (cvec.pz()<0) theta = M_PI;
00142       
00143       double ETA = cvec.pseudoRapidity();//-log(std::abs(tan(.5*theta)));//PPHOT(3)
00144       double PT  = sqrt(cvec.px()*cvec.px()+cvec.py()*cvec.py());//PPHOT(5)
00145       double PHI = cvec.phi();//ANGLE(pxpho,pypho)//PPHOT(4)
00146       
00147       double x = PT*cos(PHI);
00148       double y = PT*sin(PHI);
00149       double z = PT*sinh(ETA);
00150       double t = PT*cosh(ETA);
00151 
00152       return_vec.setPx(x);
00153       return_vec.setPy(y);
00154       return_vec.setPz(z);
00155       return_vec.setE(t);
00156       
00157     } else if (m_smearParamSchema==2){ // Smearing based on CSC sample tuning
00158 
00159       double vE = avec.e();
00160 
00161       double vLimEta=EtaPhoton[NbEtaPhoton-1]+(EtaPhoton[NbEtaPhoton-1]-EtaPhoton[NbEtaPhoton-2])*0.5;
00162     
00163       if( vEta>vLimEta )
00164         {
00165           m_log<<MSG::DEBUG<<"No calibration : "<<vEta<<endreq;
00166           return avec;
00167         }
00168     
00169       int iEnergy=0;
00170       m_log<<MSG::DEBUG<<"Hello";
00171       m_log<<MSG::DEBUG<<vE;
00172       m_log<<MSG::DEBUG<<" > ";
00173       m_log<<MSG::DEBUG<<vEnergiesPhoton[iEnergy];
00174       m_log<<MSG::DEBUG<<" ";
00175       m_log<<MSG::DEBUG<<vE<<" > "<<vEnergiesPhoton[iEnergy]<<" ";
00176       while(iEnergy<NbEnergiesPhoton&&vE>vEnergiesPhoton[iEnergy]*1000)
00177         {
00178           m_log<<MSG::DEBUG<<vEnergiesPhoton[iEnergy]<<" ";
00179           iEnergy++;
00180         }
00181       m_log<<MSG::DEBUG<<endreq;
00182       if(iEnergy==0)iEnergy++;
00183       double ratioEnergy=(vE*0.001-vEnergiesPhoton[iEnergy-1])/(vEnergiesPhoton[iEnergy]-vEnergiesPhoton[iEnergy-1]);
00184 
00185       int iEta=0;
00186       m_log<<MSG::DEBUG<<vEta<<" > "<<EtaPhoton[iEta]<<" ";
00187       while(iEta<NbEtaPhoton&&vEta>EtaPhoton[iEta])
00188         {
00189           m_log<<MSG::DEBUG<<EtaPhoton[iEta]<<" ";
00190           iEta++;
00191         }
00192       m_log<<MSG::DEBUG<<endreq;
00193 
00194       if(iEta==0)iEta++;
00195       double ratioEta=(vEta-EtaPhoton[iEta-1])/(EtaPhoton[iEta]-EtaPhoton[iEta-1]);
00196 
00197       int iPos;
00198       double sigma,sigma1,sigma2;
00199 
00200       iPos=iEnergy-1;
00201       sigma1=CorrFactorPhotonSigma[iPos][iEta-1]+ratioEta*(CorrFactorPhotonSigma[iPos][iEta]-CorrFactorPhotonSigma[iPos][iEta-1]);
00202       iPos=iEnergy;
00203       sigma2=CorrFactorPhotonSigma[iPos][iEta-1]+ratioEta*(CorrFactorPhotonSigma[iPos][iEta]-CorrFactorPhotonSigma[iPos][iEta-1]);
00204 
00205       sigma=sigma1+ratioEnergy*(sigma2-sigma1);
00206 
00207       m_log<<MSG::DEBUG<<" -> interpolation sigma "<<sigma1<<" "<<sigma2<<" "<<sigma<<endreq;
00208 
00209       double alpha=1.0;
00210       double vrandom=alpha;
00211 
00212       int iCmpt=int(randFlat()->fire()*20);
00213       for(int i=0; i<iCmpt; i++)
00214         vrandom=randGauss()->fire(1.0,sigma);
00215     
00216     
00217     
00218       m_log<<MSG::DEBUG<<" -> Energie : "<<vEnergiesPhoton[iEnergy]<<"   Eta : "<<vEnergiesPhoton[iEnergy-1]<<" < "<<vE<<" < "<<vEnergiesPhoton[iEnergy]<<endreq;
00219       m_log<<MSG::DEBUG<<" -> Energie : "<<vEnergiesPhoton[iEnergy]<<"   Eta : "<<EtaPhoton[iEta-1]<<" < "<<vEta<<" < "<<EtaPhoton[iEta]<<endreq;
00220       m_log<<MSG::DEBUG<<"       calib mean-sigma  "<<alpha<<" "<<sigma<<" -> "<<vrandom<<endreq;
00221 
00222       return_vec.set(avec.px()*vrandom, avec.py()*vrandom, avec.pz()*vrandom, vE*vrandom);
00223     
00224       m_log<<MSG::DEBUG<<"-> HepLorentzVector calibration "<<vEta<<" "<<avec<<"  ->  "<<return_vec<<endreq;
00225 
00226     }
00227 
00228     return return_vec;
00229     
00230   }

int Atlfast::PhotonSmearer::setSmearParameters ( const std::vector< double > & smearValues ) [virtual]

and the setSmearParamSchema method

Reimplemented from Atlfast::DefaultSmearer.

Definition at line 232 of file PhotonSmearer.cxx.

00232                                                                             {
00233     m_smearParams=smearValues;
00234     return 0;
00235   }

int Atlfast::PhotonSmearer::setSmearParamSchema ( const int smearSchema ) [virtual]

Sets the smearing schema

Reimplemented from Atlfast::DefaultSmearer.

Definition at line 238 of file PhotonSmearer.cxx.

00238                                                                 {
00239     m_smearParamSchema=smearSchema;
00240     return 0;
00241   }

RandGauss* Atlfast::DefaultSmearer::randGauss ( ) [inline, inherited]

Hook for subclasses etc. to get their hands on a randGauss object

Definition at line 65 of file DefaultSmearer.h.

00065 { return m_randGauss;}

RandFlat* Atlfast::DefaultSmearer::randFlat ( ) [inline, inherited]

Hook for subclasses etc. to get their hands on a randFlat object

Definition at line 71 of file DefaultSmearer.h.

00071 { return m_randFlat;}

int Atlfast::PhotonSmearer::m_lumi [private]

Definition at line 65 of file PhotonSmearer.h.

std::vector<double> Atlfast::PhotonSmearer::m_smearParams [private]

Definition at line 67 of file PhotonSmearer.h.

int Atlfast::PhotonSmearer::m_smearParamSchema [private]

Definition at line 68 of file PhotonSmearer.h.

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

Definition at line 69 of file PhotonSmearer.h.


Public Member Functions
	PhotonSmearer (const int aseed, const int lumi, MsgStream &log)
virtual	~PhotonSmearer ()
virtual HepLorentzVector	smear (const HepMC::GenParticle &)
virtual HepLorentzVector	smear (const HepLorentzVector &avec)
virtual int	setSmearParameters (const std::vector< double > &smearValues)
virtual int	setSmearParamSchema (const int smearSchema)
Private Member Functions
RandGauss *	randGauss ()
RandFlat *	randFlat ()
Private Attributes
int	m_lumi
std::vector< double >	m_smearParams
int	m_smearParamSchema
MsgStream &	m_log

Atlfast::PhotonSmearer Class Reference

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation