EmcRecShowerEnergy.cxx

Go to the documentation of this file.

 
#include "EmcRec/EmcRecShowerEnergy.h"
#include "EmcRec/EmcRecNeighbor.h"
#include "EmcRec/EmcRecParameter.h"
 
void EmcRecShowerEnergy::Energy( RecEmcShower& aShower ) {
  RecEmcFractionMap::const_iterator cit;
  RecEmcEnergy                      e1      = 0;
  RecEmcEnergy                      e9      = 0;
  RecEmcEnergy                      e25     = 0;
  RecEmcEnergy                      elepton = 0;
  RecEmcEnergy                      eall    = 0;
 
  EmcRecNeighbor nhb;
 
  RecEmcID                             CellId       = aShower.getShowerId();
  int                                  module       = EmcID::barrel_ec( CellId );
  RecEmcIDVector                       NearCell     = nhb.GetNeighbors( CellId );
  RecEmcIDVector                       NextNearCell = nhb.GetNextNeighbors( CellId );
  RecEmcIDVector                       tmpNearCell;
  RecEmcIDVector                       tmpNextNearCell;
  i_RecEmcIDVector                     pNearCell;
  i_RecEmcIDVector                     pNextNearCell;
  vector<RecEmcEnergy>                 eVec;
  vector<RecEmcEnergy>::const_iterator ciVec;
 
  tmpNearCell.push_back( CellId );
  tmpNextNearCell.push_back( CellId );
 
  cit = aShower.Find( CellId );
  // int time_seed = cit->second.getTime();
  e1 = ( cit->second.getEnergy() ) * ( cit->second.getFraction() );
  e9 += ( cit->second.getEnergy() ) * ( cit->second.getFraction() );
  e25 += ( cit->second.getEnergy() ) * ( cit->second.getFraction() );
 
  // e3x3
  for ( pNearCell = NearCell.begin(); pNearCell != NearCell.end(); pNearCell++ )
  {
    cit = aShower.Find( *pNearCell );
    if ( cit != aShower.End() )
    {
      tmpNearCell.push_back( *pNearCell );
      tmpNextNearCell.push_back( *pNearCell );
      e9 += cit->second.getEnergy() * cit->second.getFraction();
      e25 += cit->second.getEnergy() * cit->second.getFraction();
    }
  }
 
  // e5x5
  for ( pNextNearCell = NextNearCell.begin(); pNextNearCell != NextNearCell.end();
        pNextNearCell++ )
  {
    cit = aShower.Find( *pNextNearCell );
    if ( cit != aShower.End() )
    {
      tmpNextNearCell.push_back( *pNextNearCell );
      e25 += cit->second.getEnergy() * cit->second.getFraction();
    }
  }
 
  // eall
  for ( cit = aShower.Begin(); cit != aShower.End(); ++cit )
  {
    eall += ( cit->second.getEnergy() ) * ( cit->second.getFraction() );
    eVec.push_back( cit->second.getEnergy() * cit->second.getFraction() );
  }
 
  // calculate number of hits from MC
  int              nHit, n;
  EmcRecParameter& Para = EmcRecParameter::GetInstance();
  nHit = (int)( Para.HitNb( 0 ) * log( Para.HitNb( 1 ) * e9 + Para.HitNb( 2 ) ) );
  n    = 0;
 
  // sort by energy
  sort( eVec.begin(), eVec.end(), greater<RecEmcEnergy>() );
 
  for ( ciVec = eVec.begin(); ciVec != eVec.end(); ciVec++ )
  {
    if ( n < nHit )
    {
      elepton += *ciVec;
      n++;
    }
  }
 
  // energy correction
  // RecEmcEnergy eCorr=ECorrTheta(e25,CellId);
  // RecEmcEnergy eCorr=ECorrection(e25);
  int getthetaid = EmcID::theta_module( CellId );
  int getmodule  = EmcID::barrel_ec( CellId );
  if ( getthetaid > 21 ) getthetaid = 43 - getthetaid;
  if ( getmodule == 1 ) getthetaid = getthetaid + 6;
  double dthetaid = double( getthetaid );
  double eCorr    = Para.ECorrMC( e25, dthetaid );
 
  // energy error
  RecEmcEnergy de, de1, de2, de3;
  de1 = Para.SigE( 0 ) / eCorr;
  de2 = Para.SigE( 1 ) / pow( eCorr, 0.25 );
  de3 = Para.SigE( 2 );
  de  = sqrt( de1 * de1 + de2 * de2 + de3 * de3 ) * perCent * eCorr;
 
  double err = Para.ErrMC( e25, dthetaid );
  if ( err > 0 ) de = err * e25;
 
  aShower.setTrackId( -1 );
  aShower.setCellId( CellId );
  aShower.setModule( module );
  aShower.setNumHits( aShower.getSize() );
  aShower.setDE( de );
  aShower.CellId3x3( tmpNearCell );
  aShower.CellId5x5( tmpNextNearCell );
  aShower.setESeed( e1 );
  aShower.setE3x3( e9 );
  aShower.setE5x5( e25 );
  aShower.ELepton( elepton );
  aShower.EAll( eall );
  aShower.setEnergy( eCorr );
 
  // cout<<"e1="<<aShower.eSeed()
  //   <<"\te9="<<aShower.e3x3()
  //   <<"\te25="<<aShower.e5x5()
  //   <<"\telepton="<<aShower.getELepton()
  //   <<"\teall="<<aShower.getEAll()
  //   <<"\tenergy="<<aShower.energy()<<endl;
}
 
RecEmcEnergy EmcRecShowerEnergy::ECorrection( const RecEmcEnergy eIn ) {
  if ( eIn > 3. ) return eIn;
 
  EmcRecParameter& Para = EmcRecParameter::GetInstance();
 
  RecEmcEnergy eOut = 0;
  double       par[4];
  for ( int i = 0; i < 4; i++ ) { par[i] = Para.ECorr( i ); }
 
  eOut = eIn / ( par[0] + par[1] * eIn + par[2] * eIn * eIn + par[3] * eIn * eIn * eIn );
  return eOut;
}
 
RecEmcEnergy EmcRecShowerEnergy::ECorrTheta( const RecEmcEnergy eIn, const RecEmcID& id ) {
  EmcRecParameter& Para = EmcRecParameter::GetInstance();
  RecEmcEnergy     eOut = eIn;
 
  unsigned int npart  = EmcID::barrel_ec( id );
  unsigned int ntheta = EmcID::theta_module( id );
 
  if ( npart == 1 ) { eOut *= 1.843 / Para.Peak( ntheta ); }
  else if ( npart == 0 ) { eOut *= 1.843 / Para.Peak( ntheta + 44 ); }
  else if ( npart == 2 ) { eOut *= 1.843 / Para.Peak( ntheta + 50 ); }
 
  return eOut;
}