Athena-Atlfast

ATLFAST SMEARERS

The electron smearer is parameterised as follows:
// 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

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

while ( 1 ) {
    aa = randGauss( )->fire( );
    //bb = randGauss( )->fire( );
    sigph1 = aa*0.12/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/avec.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 )/avec.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/avec.perp( );
            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:
HepLorentzVector bvec( avec.px( )*( 1.0+sigph ), avec.py( )*( 1.0+sigph ), avec.pz( )*( 1.0+sigph ), avec.e( ) *( 1.0+sigph ) );

return bvec;
The photon smearer is parameterised as follows:
// do the smearing, copied verbatim ( except for ROOT dependencies )
// from photonmaker codein Atlfast++
// 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.
//

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;
The muon smearer is parameterised as follows:
// ATLFAST fortran routines
// RESMUO( KEYMUO, KEYFUN, PT, ETA, PHI )
// RESOMU ( this had been reduced to merely converting radians to degrees )
// have been converted to C++ here, after restructuring the code.
// RESOLUMU ( in Atlfast/mu_atlfast/resmuo.F ) is left in fortran.

float eta = avec.pseudoRapidity( );
float pt = avec.perp( );
float phi = avec.phi( )*deg;
float sigmams;
float sigmamuon;
float sigmaid;
float sigmatrack;
float sigma = 0.;
pair<float, float> sigmapr;
switch ( m_keymuo ) {
    case 1:
        while( ( ( sigmapr = resolms( pt, eta, phi ) ).first ) <-1. ) {}
        sigma = sigmapr.first;
    case 2:
        while( ( sigma = resolid1( pt, eta ) ) <-1. ) {}
    case 3:
        while( ( ( sigmapr = resolms( pt, eta, phi ) ).first ) <-1. ) {}
        sigmams = sigmapr.first;
        sigmamuon = sigmapr.second;
        while( ( ( sigmapr = resolid2( pt, eta ) ).first ) <-1. ) {}
        sigmaid = sigmapr.first;
        sigmatrack = sigmapr.second;

        sigma = resol1( sigmams, sigmamuon, sigmaid, sigmatrack );
    case 4:
        while( ( ( sigmapr = resolms( pt, eta, phi ) ).first ) <-1. ) {}
        sigmams = sigmapr.first;
        sigmamuon = sigmapr.second;
        while( ( ( sigmapr = resolid2( pt, eta ) ).first ) <-1. ) {}
        sigmaid = sigmapr.first;
        sigmatrack = sigmapr.second;

        sigma = resol2( sigmamuon, sigmatrack );
    case 5:
        sigma = resol3( );
    case 6:
        sigma = resol4( );

    }

HepLorentzVector bvec( avec.px( )/( 1.0+sigma ), avec.py( )/( 1.0+sigma ),
avec.pz( )/( 1.0+sigma ), avec.e( ) /( 1.0+sigma ) );
return bvec;
}

pair<float, float> MuonSmearer::resolms ( float pt, float eta, float phi ) {
    float cpt = pt;
    float ceta = eta;
    float cphi = phi;
    float resol;
    resolumu_( &cpt, &ceta, &cphi, &resol );
    double aa = randGauss( )->fire( );
    float sigmamuon = m_percent*resol;
    float sigmams = aa*sigmamuon;
    pair<float, float> sigmapr( sigmams, sigmamuon );

    return sigmapr;
    }

float MuonSmearer::resolid1 ( float pt, float eta ) {
    double aa = randGauss( )->fire( );
    double bb = randGauss( )->fire( );
    float factor = 1. + pow( fabs( eta ), m_magic2 );
    float atrack = m_magic1*factor;
    float sigma = atrack*pt*aa + m_magic3*bb;
    sigma = sigma*factor;

    return sigma;
    }

pair<float, float> MuonSmearer::resolid2 ( float pt, float eta ) {
    double aa = randGauss( )->fire( );
    double bb = randGauss( )->fire( );
    double factor = 1. + pow( fabs( eta ), m_magic2 );
    float atrack = m_magic1*factor;
    float sigmatr = atrack*pt*aa + m_magic3*bb;
    float sigmatrack = sqrt( ( atrack*pt )*( atrack*pt ) + m_magic3*m_magic3 );
    pair<float, float> sigmapr( sigmatr, sigmatrack );

    return sigmapr;
    }

float MuonSmearer::resol1 ( float sigmams, float sigmamuon,
float sigmaid, float sigmatrack ) {
    double wmuon = 1./( sigmamuon*sigmamuon );
    double wtrack = 1./( sigmatrack*sigmatrack );
    double wtot = wmuon+wtrack;
    double corr = ( wmuon*( 1.0+sigmams )+wtrack*( 1.0+sigmaid ) )/wtot;
    float sigma = corr-1.0;

    return sigma;
    }

float MuonSmearer::resol2 ( float sigmamuon, float sigmatrack ) {
    double wmuon = 1./( sigmamuon*sigmamuon );
    double wtrack = 1./( sigmatrack*sigmatrack );
    double wtot = wmuon+wtrack;
    double aa = randGauss( )->fire( );
    float sigma = aa/sqrt( wtot );

    return sigma;
    }

float MuonSmearer::resol3 ( ) {
    float sigma = m_magic4*randGauss( )->fire( );

    return sigma;
    }

float MuonSmearer::resol4 ( ) {
    float sigmams = m_magic4;
    float sigmaid = m_magic4;
    float wmuon = 1.0/( sigmams*sigmams );
    float wtrak = 1.0/( sigmaid*sigmaid );
    float wtot = wmuon + wtrak;
    double aa = randGauss( )->fire( );
    float sigma = aa/sqrt( wtot );

    return sigma;