---

PhotonSmearer.cxx

Go to the documentation of this file.

// PhotonSmearer.cxx
//
// Implementation of the Photon Smearer class
//
// Only the smear() method is overridden
//
//
// Authors: H.T. Phillips, P.Clarke, E.Richter-Was, P.Sherwood, R.Steward
//
//

#include "AtlfastAlgs/PhotonSmearer.h"

#include <cmath>

#include "CLHEP/Vector/LorentzVector.h"

#include "CLHEP/Random/JamesRandom.h"
#include "CLHEP/Random/RandGauss.h"

namespace Atlfast {
  using std::abs;

  HepLorentzVector PhotonSmearer::smear (const HepLorentzVector& avec) {
    // do the smearing, copied verbatim (except for ROOT dependencies)
    // from photonmaker codein Atlfast++
    // The code is very procedural (even more so than Atlfast++!) and should be tidied later!
    //
    // Parametrisation was by L.Poggioli and F. Gianotti

    double rpilup = 0.0;
    double aa, bb, sigph1, sigph, sigpu;
    double sqrtene = sqrt(avec.e());
    double abseta  = fabs(avec.pseudoRapidity());

    // do smearing of theta first
    HepLorentzVector bvec(avec);
    
    double sigph2 = 0.0;
    while (1) {
      aa=randGauss()->fire();
      if (abseta < 0.8) sigph2                  = aa*0.065/sqrtene;
      if (abseta >= 0.8 && abseta < 1.4) sigph2 = aa*0.050/sqrtene;
      if (abseta >= 1.4 && abseta < 2.5) sigph2 = aa*0.040/sqrtene;
      if (abseta >= 2.5)                 sigph2 = 0;
      if (fabs(bvec.theta()) + sigph2 <= M_PI) break;
    }

    // sigph2 is now the offset for theta...
    bvec.setPz(bvec.e()*cos(bvec.theta()+sigph2));


    // now do the energy resolution
    //
    while (1) {
      aa=randGauss()->fire();
      //bb=randGauss()->fire();
      sigph1 = aa*0.10/sqrtene;
      if (1.0 + sigph1 <= 0.0) continue;
      sigph = sigph1;
      if (abseta < 1.4) {
        while (1) {
          aa=randGauss()->fire();
          bb=randGauss()->fire();
          sigph1 = aa*0.245/bvec.perp() + bb*0.007;
          if (1.0+sigph1 > 0) break;
        }
      } else {
        while (1) {
          aa=randGauss()->fire();
          bb=randGauss()->fire();
          sigph1 = aa*0.306*((2.4-abseta)+0.228)/bvec.e() + bb*0.007;
          if (1.0+sigph1 > 0) break;
        }
      }
      sigph += sigph1;
      if (m_lumi == 2) {
        if (abseta < 0.6) rpilup = 0.32;
        if (abseta > 0.6 && abseta < 1.4) rpilup = 0.295;
        if (abseta > 1.4) rpilup = 0.27;
        while (1) {
          aa=randGauss()->fire();
          //bb=randGauss()->fire();
          sigpu = aa*rpilup/bvec.perp();
          if (1+sigpu > 0) break;
        }
        sigph += sigpu;
      }
      if (1.0 + sigph > 0) break;
    }


    // Now sigph is the photon "sigma" and we do the following a la Atlfast++ photon maker:
    // this seems to be OK for energy resolution. 
    //
    // It looks to me as though the functional form above for photons is actually the same as
    // electrons, just that one magic number is different. Suggest we revisit this and
    // exploit commonality.
    //

    bvec.setPx(bvec.px()*(1.0+sigph));
    bvec.setPy(bvec.py()*(1.0+sigph));
    bvec.setPz(bvec.pz()*(1.0+sigph));

    // nb from Atlfast++ code, we go back to the original photon energy here
    bvec.setE( avec.e() *(1.0+sigph));
  
    return bvec;
  }

} // end of namespace bracket

---