McCor.cxx

Go to the documentation of this file.

#include "EventModel/Event.h"
#include "EventModel/EventModel.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/PropertyMgr.h"
#include "GaudiKernel/SmartDataPtr.h"
 
#include "EvtRecEvent/EvtRecEvent.h"
#include "EvtRecEvent/EvtRecTrack.h"
 
#include "CLHEP/Vector/LorentzVector.h"
#include "CLHEP/Vector/ThreeVector.h"
#include "CLHEP/Vector/TwoVector.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IHistogramSvc.h"
#include "GaudiKernel/INTupleSvc.h"
#include "GaudiKernel/NTuple.h"
#include "TMath.h"
using CLHEP::Hep2Vector;
using CLHEP::Hep3Vector;
using CLHEP::HepLorentzVector;
#include "CLHEP/Geometry/Point3D.h"
#ifndef ENABLE_BACKWARDS_COMPATIBILITY
typedef HepGeom::Point3D<double> HepPoint3D;
#endif
#include "McCor.h"
 
#include "Identifier/Identifier.h"
 
#include "TGraph2D.h"
#include "TGraph2DErrors.h"
#include "TGraphErrors.h"
#include <fstream>
#include <string>
#include <strstream>
#include <vector>
 
/////////////////////////////////////////////////////////////////////////////
//
TGraph2DErrors* dt = new TGraph2DErrors();
 
DECLARE_COMPONENT( McCor )
McCor::McCor( const std::string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator ) {
 
  // Declare the properties
  ntOut = true;
  declareProperty( "NTupleOut", ntOut );
 
  /*
    declareProperty("b0",m_b[0] = 0.9976);
    declareProperty("b1",m_b[1] = -0.01718);
    declareProperty("b2",m_b[2] = 0.01634);
  */
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
StatusCode McCor::initialize() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in initialize()" << endmsg;
  StatusCode status;
  if ( ntOut == true )
  {
    NTuplePtr nt1( ntupleSvc(), "FILE1/ec" );
    if ( nt1 ) m_tuple1 = nt1;
    else
    {
      m_tuple1 = ntupleSvc()->book( "FILE1/ec", CLID_ColumnWiseTuple, "ks N-Tuple example" );
      if ( m_tuple1 )
      {
        status = m_tuple1->addItem( "ef", m_ef );
        status = m_tuple1->addItem( "e5", m_e5 );
        status = m_tuple1->addItem( "ec", m_ec );
        status = m_tuple1->addItem( "ct", m_ct );
      }
      else
      {
        log << MSG::ERROR << "    Cannot book N-tuple:" << long( m_tuple1 ) << endmsg;
        return StatusCode::FAILURE;
      }
    }
  }
  /*
    string CoeffPath=getenv("MCCORROOT");
    CoeffPath +="/share/McCorCoeff.txt";
    ifstream in;
    in.open(CoeffPath.c_str(),ios::in);
    for(int i=0;i<4;i++){
      in>>m_a[i];
      in>>m_ae[i];
    }
  */
  double energy, thetaid, peak, peakerr, res, reserr;
  string DataPath;
  DataPath = getenv( "MCCORROOT" );
  DataPath += "/share/evset.txt";
  ifstream in1;
  in1.open( DataPath.c_str(), ios::in );
  //  in.open("$MCCORROOT/share/evsc.txt");
  double ep[18] = { 0.03, 0.04, 0.05, 0.075, 0.1, 0.125, 0.15, 0.2,  0.25,
                    0.3,  0.4,  0.5,  0.75,  1.0, 1.25,  1.5,  1.75, 2.0 };
  for ( int i = 0; i < 504; i++ )
  {
    in1 >> energy;
    in1 >> thetaid;
    in1 >> peak;
    in1 >> peakerr;
    in1 >> res;
    in1 >> reserr;
    int j = i / 28;
    dt->SetPoint( i, energy, thetaid, peak );
    dt->SetPointError( i, 0, 0, peakerr );
  }
  in1.close();
  log << MSG::INFO << "successfully return from initialize()" << endmsg;
  return StatusCode::SUCCESS;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
StatusCode McCor::execute() {
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in execute()" << endmsg;
 
  SmartDataPtr<Event::EventH> eventHeader( eventSvc(), "/Event" );
  SmartDataPtr<EvtRecEvent>   evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
  //  log << MSG::INFO << "get event tag OK" << endmsg;
  log << MSG::DEBUG << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
      << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
 
  SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
 
  for ( int i = 0; i < evtRecEvent->totalTracks(); i++ )
  {
 
    EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
    if ( !( *itTrk )->isEmcShowerValid() ) continue;
    RecEmcShower* emcTrk = ( *itTrk )->emcShower();
 
    int        intid      = emcTrk->cellId();
    Identifier id         = Identifier( intid );
    int        getthetaid = EmcID::theta_module( id );
    int        getmodule  = EmcID::barrel_ec( id );
    if ( getthetaid > 21 ) getthetaid = 43 - getthetaid;
    if ( getmodule == 1 ) getthetaid = getthetaid + 6;
    double energyF  = emcTrk->energy();
    double e5x5     = emcTrk->e5x5();
    double costheta = cos( emcTrk->theta() );
    double dthetaid = double( getthetaid );
    double energyC  = McCor::corEnergyMc( e5x5, dthetaid );
 
    if ( ntOut == true )
    {
      m_ct = costheta;
      m_ef = energyF;
      m_e5 = e5x5;
      m_ec = energyC;
      m_tuple1->write();
    }
    emcTrk->setEnergy( energyC );
  }
 
  return StatusCode::SUCCESS;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
StatusCode McCor::finalize() {
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in finalize()" << endmsg;
  return StatusCode::SUCCESS;
}
 
double McCor::corEnergyMc( double eg, double theid ) {
 
  double Energy5x5 = eg;
  if ( eg < 0.029 ) eg = 0.029;
  if ( eg > 1.76 ) eg = 1.76;
  if ( theid <= 0 ) theid = 0.001;
  if ( theid >= 27 ) theid = 26.999;
  Float_t einter = eg + 0.00001;
  Float_t tinter = theid + 0.0001;
  double  ecor   = dt->Interpolate( einter, tinter );
  if ( !( ecor ) ) return Energy5x5;
  if ( ecor < 0.5 ) return Energy5x5;
  double EnergyCor = Energy5x5 / ecor;
  return EnergyCor;
}