Report problems to ATLAS LXR Team (with time and IP address indicated)

The LXR Cross Referencer
[ source navigation ]
[ diff markup ]
[ identifier search ]
[ general search ]
 
 
Architecture: [ linux ]
Version: [ head ] [ nightly ] [ GaudiDev ]
  Links to LXR source navigation pages for stable releases [ 12.*.* ]   [ 13.*.* ]   [ 14.*.* ]   [ 15.*.* ]  
001 // PhotonReconstructor.cxx
002 //
003 // Implementation of the Photon Smearer class
004 //
005 // Only the reconstruct() method is overridden
006 //
007 //
008 // Authors: H.T. Phillips, P.Clarke, E.Richter-Was, P.Sherwood, R.Steward
009 // Update:  Sabine Elles, Simon Dean
010 //
011 // Smearer - Reconstructor Conversion: Alexander Richards (UCL)
012 //
013 
014 #include "AtlfastAlgs/PhotonReconstructor.h"
015 #include "AtlfastAlgs/CorrectionFactorPhoton.h"
016 
017 #include <cmath>
018 
019 #include "CLHEP/Vector/LorentzVector.h"
020 
021 #include "CLHEP/Random/JamesRandom.h"
022 #include "CLHEP/Random/RandGauss.h"
023 #include "CLHEP/Units/SystemOfUnits.h"
024 
025 namespace Atlfast {
026   using std::abs;
027   
028   ReconstructedParticle PhotonReconstructor::reconstruct( const ReconstructedParticle& particle )const{
029     
030     double rpilup = 0.0;
031     double aa, bb, sigph1, sigph, sigpu;
032     double vEta  = fabs(particle.eta());
033     
034     //make sure e, momentum in GeV
035     double pt = particle.pt()/GeV;
036     double ene = particle.e()/GeV;
037     double sqrtene = sqrt(ene);
038     
039     // The HLV to be returned
040     HepLorentzVector return_vec(0.,0.,0.,0.);
041     if (!sqrtene || !ene || !pt) {
042       ReconstructedParticle return_particle(particle);
043       return_particle.set_momentum(return_vec);
044       return return_particle;
045     }
046     if (m_smearParamSchema==1){ // TDR-style smearing
047       
048       // do the smearing, copied verbatim (except for ROOT dependencies)
049       // from photonmaker codein Atlfast++
050       // The code is very procedural (even more so than Atlfast++!) and should be tidied later!
051       //
052       // Parametrisation was by L.Poggioli and F. Gianotti
053       
054       // do smearing of theta first
055       HepLorentzVector bvec(particle.momentum());
056       
057       //this is equivalent to RESTHE in FORTRAN 
058       double sigph2 = 0.0;
059       while (1) {
060         aa=m_randGauss->fire();
061         if (vEta < 0.8) sigph2                  = aa*m_smearParams[0]/sqrtene;
062         //if (vEta < 0.8) sigph2                  = aa*0.065/sqrtene;
063         if (vEta >= 0.8 && vEta < 1.4) sigph2 = aa*m_smearParams[1]/sqrtene;  
064         //if (vEta >= 0.8 && vEta < 1.4) sigph2 = aa*0.050/sqrtene;
065         if (vEta >= 1.4 && vEta < 2.5) sigph2 = aa*m_smearParams[2]/sqrtene; 
066         //if (vEta >= 1.4 && vEta < 2.5) sigph2 = aa*0.040/sqrtene;
067         if (vEta >= 2.5)                 sigph2 = 0;
068         if (fabs(bvec.theta()) + sigph2 <= M_PI) break;
069       }
070       
071       // sigph2 is now the offset for theta...
072       bvec.setPz(bvec.e()*cos(bvec.theta()+sigph2));
073       
074       vEta  = fabs(bvec.pseudoRapidity()); //27/02/03 JPC possible bug fix
075       
076       //this is RESPHO in FORTRAN
077       // now do the energy resolution note: use original energy, pt and phi, but new vEta
078       //
079       while (1) {
080         while (1){
081           aa=m_randGauss->fire();
082           sigph1 = aa*m_smearParams[3]/sqrtene;
083           //sigph1 = aa*0.10/sqrtene;
084           if (1.0 + sigph1 > 0.0) break;
085         }
086         sigph = sigph1;
087         if (vEta < 1.4) {
088           while (1) {
089             aa=m_randGauss->fire();
090             bb=m_randGauss->fire();
091             sigph1 = aa*m_smearParams[4]/pt + bb*m_smearParams[5];
092             //sigph1 = aa*0.245/pt + bb*0.007;
093             if (1.0+sigph1 > 0) break;
094           }
095         } else {
096           while (1) {
097             aa=m_randGauss->fire();
098             bb=m_randGauss->fire();
099             sigph1 = aa*m_smearParams[6]*((m_smearParams[7]-vEta)+m_smearParams[8])/ene + bb*m_smearParams[9];
100             // sigph1 = aa*0.306*((2.4-vEta)+0.228)/ene + bb*0.007;
101             if (1.0+sigph1 > 0) break;
102           }
103         }
104         sigph += sigph1;
105         if (m_lumi == 2) {
106           if (vEta < 0.6) rpilup = 0.32;
107           if (vEta > 0.6 && vEta < 1.4) rpilup = 0.295;
108           if (vEta > 1.4) rpilup = 0.27;
109           while (1) {
110             aa=m_randGauss->fire();
111             sigpu = aa*rpilup/pt;
112             if (1+sigpu > 0) break;
113           }
114           sigph += sigpu;
115         }
116         
117         if (1.0 + sigph > 0) break;
118       }
119       
120       
121       // Now sigph is the photon "sigma" and we do the following a la Atlfast++ photon maker:
122       // this seems to be OK for energy resolution. 
123       //
124       // It looks to me as though the functional form above for photons is actually the same as
125       // electrons, just that one magic number is different. Suggest we revisit this and
126       // exploit commonality.
127       //
128       
129       HepLorentzVector cvec;
130       cvec.setPx(particle.px()*(1.0+sigph));
131       cvec.setPy(particle.py()*(1.0+sigph));
132       cvec.setPz(bvec.pz()*(1.0+sigph));//using recalculated z momentum
133       
134       // nb from FORTRAN code, we go back to the original photon energy here
135       cvec.setE( particle.e() *(1.0+sigph));
136       //return cvec;
137       //now do what the fortran does to get massless particle! YUCK!!!
138       
139       //if(std::abs (cvec.pz()/cvec.e()) > 1.){
140       //        theta = acos(cvec.pz()/sqrt(cvec.px()*cvec.px()+cvec.py()*cvec.py()+cvec.pz()*cvec.pz()));
141       //} else if (cvec.pz()>0) theta = 0.;
142       //else if (cvec.pz()<0) theta = M_PI;
143       
144       double ETA = cvec.pseudoRapidity();//-log(std::abs(tan(.5*theta)));//PPHOT(3)
145       double PT  = sqrt(cvec.px()*cvec.px()+cvec.py()*cvec.py());//PPHOT(5)
146       double PHI = cvec.phi();//ANGLE(pxpho,pypho)//PPHOT(4)
147       
148       double x = PT*cos(PHI);
149       double y = PT*sin(PHI);
150       double z = PT*sinh(ETA);
151       double t = PT*cosh(ETA);
152 
153       return_vec.setPx(x);
154       return_vec.setPy(y);
155       return_vec.setPz(z);
156       return_vec.setE(t);
157       
158     } else if (m_smearParamSchema==2){ // Smearing based on CSC sample tuning - - data in CorrectionFactorPhoton.h
159 
160       double vE = particle.e();
161 
162       double vLimEta=EtaPhoton[NbEtaPhoton-1]+(EtaPhoton[NbEtaPhoton-1]-EtaPhoton[NbEtaPhoton-2])*0.5;
163     
164       if( vEta>vLimEta )
165         {
166           m_log<<MSG::DEBUG<<"No calibration : "<<vEta<<endreq;
167           return particle;
168         }
169     
170       int iEnergy=0;
171       m_log<<MSG::DEBUG<<vE<<" > "<<vEnergiesPhoton[iEnergy]<<" ";
172       while(iEnergy<NbEnergiesPhoton&&vE>vEnergiesPhoton[iEnergy]*1000)
173         {
174           m_log<<MSG::DEBUG<<vEnergiesPhoton[iEnergy]<<" ";
175           iEnergy++;
176         }
177       m_log<<MSG::DEBUG<<endreq;
178       if(iEnergy==0)iEnergy++;
179       double ratioEnergy=(vE*0.001-vEnergiesPhoton[iEnergy-1])/(vEnergiesPhoton[iEnergy]-vEnergiesPhoton[iEnergy-1]);
180 
181       int iEta=0;
182       m_log<<MSG::DEBUG<<vEta<<" > "<<EtaPhoton[iEta]<<" ";
183       while(iEta<NbEtaPhoton&&vEta>EtaPhoton[iEta])
184         {
185           m_log<<MSG::DEBUG<<EtaPhoton[iEta]<<" ";
186           iEta++;
187         }
188       m_log<<MSG::DEBUG<<endreq;
189 
190       if(iEta==0)iEta++;
191       double ratioEta=(vEta-EtaPhoton[iEta-1])/(EtaPhoton[iEta]-EtaPhoton[iEta-1]);
192 
193       int iPos;
194       double sigma,sigma1,sigma2;
195 
196       iPos=iEnergy-1;
197       sigma1=CorrFactorPhotonSigma[iPos][iEta-1]+ratioEta*(CorrFactorPhotonSigma[iPos][iEta]-CorrFactorPhotonSigma[iPos][iEta-1]);
198       iPos=iEnergy;
199       sigma2=CorrFactorPhotonSigma[iPos][iEta-1]+ratioEta*(CorrFactorPhotonSigma[iPos][iEta]-CorrFactorPhotonSigma[iPos][iEta-1]);
200 
201       sigma=sigma1+ratioEnergy*(sigma2-sigma1);
202 
203       m_log<<MSG::DEBUG<<" -> interpolation sigma "<<sigma1<<" "<<sigma2<<" "<<sigma<<endreq;
204 
205       double alpha=1.0;
206       double vrandom=alpha;
207 
208       int iCmpt=int(m_randFlat->fire()*20);
209       for(int i=0; i<iCmpt; i++)
210         vrandom=m_randGauss->fire(1.0,sigma);
211     
212     
213     
214       m_log<<MSG::DEBUG<<" -> Energie : "<<vEnergiesPhoton[iEnergy]<<"   Eta : "<<vEnergiesPhoton[iEnergy-1]<<" < "<<vE<<" < "<<vEnergiesPhoton[iEnergy]<<endreq;
215       m_log<<MSG::DEBUG<<" -> Energie : "<<vEnergiesPhoton[iEnergy]<<"   Eta : "<<EtaPhoton[iEta-1]<<" < "<<vEta<<" < "<<EtaPhoton[iEta]<<endreq;
216       m_log<<MSG::DEBUG<<"       calib mean-sigma  "<<alpha<<" "<<sigma<<" -> "<<vrandom<<endreq;
217 
218       return_vec.set(particle.px()*vrandom, particle.py()*vrandom, particle.pz()*vrandom, vE*vrandom);
219     
220       m_log<<MSG::DEBUG<<"-> HepLorentzVector calibration "<<vEta<<" "<<particle<<"  ->  "<<return_vec<<endreq;
221 
222     }
223 
224     ReconstructedParticle return_particle(particle);
225     return_particle.set_momentum(return_vec);
226     return return_particle;
227     
228   }
229   //cct: implement the setSmearParameters method
230   int PhotonReconstructor::setSmearParameters (const std::vector<double>& smearValues){
231     m_smearParams=smearValues;
232     return 0;
233   }
234   
235   //cct: implement the setSmearParamSchema method
236   int PhotonReconstructor::setSmearParamSchema ( const int smearSchema) {
237     m_smearParamSchema=smearSchema;
238     return 0;
239   }
240   
241 } // end of namespace bracket
242
[ source navigation ] [ diff markup ] [ identifier search ] [ general search ]
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. Valid HTML 4.01!