AbsCor Class Reference

#include <AbsCor.h>

Inheritance diagram for AbsCor:

Public Member Functions
	AbsCor (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode	initialize ()
StatusCode	execute ()
StatusCode	finalize ()
double	corEnergyPi0 (double eg, double theid)

Detailed Description

Definition at line 11 of file AbsCor.h.

Constructor & Destructor Documentation

AbsCor()

AbsCor::AbsCor	(	const std::string &	name,
		ISvcLocator *	pSvcLocator )

Definition at line 54 of file AbsCor.cxx.

    : Algorithm( name, pSvcLocator ) {
 
  // Declare the properties
  ntOut = true;
  declareProperty( "NTupleOut", ntOut = false );
  declareProperty( "McCor", mccor = 0 );
  declareProperty( "EdgeCor", edgecor = 0 );
  declareProperty( "HotCellMask", hotcellmask = 0 );
  declareProperty( "UseTof", usetof = 1 );
  declareProperty( "DoDataCor", dodatacor = 1 );
  declareProperty( "DoPi0Cor", dopi0Cor = 1 );
  declareProperty( "MCUseTof", MCuseTof = 1 );
  declareProperty( "MCCorUseFunction", MCCorUseFunction = 1 );
  declareProperty( "IYear", IYear = 2009 );
  declareProperty( "ifReadDB", m_ifReadDB = false );
  declareProperty( "CorFunparaPath",
                   m_CorFunparaPath = "/afs/ihep.ac.cn/bes3/offline/CalibConst/emc/ShEnCalib/"
                                      "7.0.9/evsetTofCorFunctionPar2021psip.txt" );
  declareProperty(
      "DataPathc3ptof",
      m_DataPathc3ptof =
          "/afs/ihep.ac.cn/bes3/offline/CalibConst/emc/ShEnCalib/7.0.9/c3ptof2021psip.txt" );
 
  runFrom  = 0;
  runTo    = 0;
  m_inFlag = false;
}

Referenced by AbsCor().

Member Function Documentation

corEnergyPi0()

double AbsCor::corEnergyPi0	(	double	eg,
		double	theid )

execute()

StatusCode AbsCor::execute

(

)

to read the parameters of energy correction parameters of single gamma and pi0 calibration//////////

Definition at line 340 of file AbsCor.cxx.

                           {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in execute()" << endmsg;
 
  SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
  int                              runNo = eventHeader->runNumber();
 
  //// by liucx
  int runNum = runNo;
  if ( runNo < 0 ) runNum = -runNo; // for MC
  if ( runNum >= runFrom && runNum <= runTo ) { m_inFlag = true; }
  else { m_inFlag = false; }
 
  if ( m_inFlag == false )
  {
 
    runFrom = m_EmcShEnCalibSvc->getRunFrom();
    runTo   = m_EmcShEnCalibSvc->getRunTo();
 
    ////to read the parameters of energy correction parameters of single gamma and pi0
    /// calibration//////////
    // Read energy correction Function parameters
    cout << "current run=" << runNo << endl;
    cout << "in AbsCor  open getPi0CalibFile()= " << m_EmcShEnCalibSvc->getPi0CalibFile()
         << endl;
    cout << "open getSingleGammaCalibFile()= " << m_EmcShEnCalibSvc->getSingleGammaCalibFile()
         << endl;
    cout << "getRunFrom()=" << runFrom << endl;
    cout << "getRunTo()=" << runTo << endl;
 
    string CorFunparaPath;
    if ( MCuseTof == 1 )
    {
      if ( m_ifReadDB ) { CorFunparaPath = m_EmcShEnCalibSvc->getSingleGammaCalibFile(); }
      else
      {
        CorFunparaPath = m_CorFunparaPath;
        cout << "no read database, but using the set:" << CorFunparaPath << endl;
      }
    }
    if ( MCuseTof == 0 )
    {
      CorFunparaPath = getenv( "ABSCORROOT" );
      CorFunparaPath += "/share/evsetCorFunctionPar.txt";
    }
 
    ifstream in2corfun;
    in2corfun.open( CorFunparaPath.c_str() );
    assert( in2corfun );
 
    for ( int i = 0; i < 28; i++ )
    {
      in2corfun >> m_corFunPar[i][0];
      in2corfun >> m_corFunPar[i][1];
      in2corfun >> m_corFunPar[i][2];
      in2corfun >> m_corFunPar[i][3];
      in2corfun >> m_corFunPar[i][4];
      in2corfun >> m_corFunPar[i][5];
    }
    in2corfun.close();
 
    ////////read the parameter of pi0 calibration for data
    string DataPathc3p;
    DataPathc3p = getenv( "ABSCORROOT" );
    DataPathc3p += "/share/c3p.txt";
 
    string DataPathc3ptof;
    if ( m_ifReadDB ) { DataPathc3ptof = m_EmcShEnCalibSvc->getPi0CalibFile(); }
    else
    {
      DataPathc3ptof = m_DataPathc3ptof;
      cout << "no read database, but using the set:" << DataPathc3ptof << endl;
    }
 
    ifstream inc3p;
    inc3p.exceptions( ifstream::failbit | ifstream::badbit );
    if ( usetof == 0 ) inc3p.open( DataPathc3p.c_str(), ios::in );
    if ( usetof == 1 ) inc3p.open( DataPathc3ptof.c_str(), ios::in );
    for ( int i = 0; i < 4; i++ )
    {
      double am, ame;
      inc3p >> am;
      inc3p >> ame;
      ai[i] = am;
    }
  }
  ////////////////////////end of reading the correction parameters from Svc and file
 
  SmartDataPtr<EvtRecEvent>    evtRecEvent( eventSvc(), EventModel::EvtRec::EvtRecEvent );
  SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), EventModel::EvtRec::EvtRecTrackCol );
  if ( evtRecEvent->totalTracks() > evtRecTrkCol->size() ) return StatusCode::SUCCESS;
  if ( evtRecEvent->totalTracks() > 50 ) return StatusCode::SUCCESS;
 
  IEmcRecGeoSvc* iGeoSvc;
  ISvcLocator*   svcLocator = Gaudi::svcLocator();
  StatusCode     sc         = svcLocator->service( "EmcRecGeoSvc", iGeoSvc );
  if ( sc != StatusCode::SUCCESS ) { cout << "Error: Can't get EmcRecGeoSvc" << endl; }
 
  for ( int i = 0; i < evtRecEvent->totalTracks(); i++ )
  {
    EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
    if ( !( *itTrk )->isEmcShowerValid() ) continue;
    RecEmcShower* emcTrk = ( *itTrk )->emcShower();
    //    if(emcTrk->energy()<0.0005)continue;
 
    /* by liucx
    //the correction is based on e5x5 energy from the version AbsCor-00-00-29.
    //It's impossible to correct twice.
    //So the code is commented out.
 
    int st = emcTrk->status();
    if(st>10000)continue;
     emcTrk->setStatus(st+20000);
    */
 
    //    cout<<"REC after calib EMC status =  "<<emcTrk->status()<<endl;
 
    // if(emcTrk->e5x5()<0.015)continue;
 
    //    if(emcTrk->getTime()<10||emcTrk->getTime()>30)continue;
    /*
        if(emcTrk->e5x5()<0.02){
          emcTrk->setEnergy(emcTrk->e5x5()*1.1);
          continue;
        }
    */
    /*
        int intid = emcTrk->cellId();
        Identifier id=Identifier(intid);
        int getthetaid = EmcID::theta_module(id);
        int getmodule = EmcID::barrel_ec(id);
        if(getmodule==1)getthetaid=getthetaid+6;
        if(getmodule==2)getthetaid=55-getthetaid;
    */
 
    RecEmcID id = Identifier( emcTrk->cellId() );
 
    unsigned int module, ntheta, nphi;
    module = EmcID::barrel_ec( id );
    ntheta = EmcID::theta_module( id );
    nphi   = EmcID::phi_module( id );
 
    // id=EmcID::crystal_id(module,ntheta,nphi);
 
    unsigned int thetaModule = EmcID::theta_module( id );
    unsigned int phiModule   = EmcID::phi_module( id );
 
    double e5x5 = emcTrk->e5x5();
    double etof = 0;
 
    if ( usetof == 1 && ( *itTrk )->isTofTrackValid() )
    {
      SmartRefVector<RecTofTrack> recTofTrackVec = ( *itTrk )->tofTrack();
      if ( !recTofTrackVec.empty() ) etof = recTofTrackVec[0]->energy();
      if ( etof > 100. ) etof = 0;
    }
 
    double energyC = 0;
 
    int thetaId = -99;
    if ( module == 0 || module == 2 ) thetaId = thetaModule;
    if ( module == 1 && thetaModule <= 21 ) thetaId = thetaModule + 6;
    if ( module == 1 && thetaModule > 21 ) thetaId = 43 - thetaModule + 6;
    double DthetaId = double( thetaId );
 
    if ( MCuseTof == 1 )
    {
      if ( thetaId < 6 ) etof = 0.0; // in EMC endcap
      if ( MCCorUseFunction == 1 ) { energyC = ECorrFunctionMC( e5x5 + etof, DthetaId ); }
      else if ( MCCorUseFunction == 0 ) { energyC = ECorrMC( e5x5 + etof, DthetaId ); }
    }
    if ( MCuseTof == 0 )
    {
      if ( MCCorUseFunction == 1 ) { energyC = ECorrFunctionMC( e5x5, DthetaId ); }
      else if ( MCCorUseFunction == 0 ) { energyC = ECorrMC( e5x5, DthetaId ); }
    }
 
    /*
    if(usetof==1){
      energyC=ECorrMC(e5x5+etof,thetaId);
    }
    if (usetof==0) {
      energyC=emcTrk->energy();
    }
    */
    //  double energyC=emcTrk->energy()+etof;
    /*
        if(energyC<=0||energyC>4.99)continue;
        double cor = 1.0;
        if(runNo>0)cor = dt->Eval(energyC);
        if(cor<0.001)continue;
    //  cout<<cor<<endl;
        double energyCC= energyC/cor;
        emcTrk->setEnergy(energyCC);
    */
 
    double lnE = std::log( energyC );
 
    if ( energyC > 1.0 )
      lnE = std::log( 1.0 ); // gamma energy bin from 0.05 to 0.9GeV for pi0 calibration
    if ( energyC < 0.05 )
      lnE = std::log( 0.05 ); // set the range(0.05-1.0GeV) of application for pi0 calibration
                              // function
 
    double lnEcor = 1.0;
    if ( dopi0Cor == 1 )
    {
      if ( runNo > 0 && dodatacor == 1 )
      { lnEcor = ai[0] + ai[1] * lnE + ai[2] * lnE * lnE + ai[3] * lnE * lnE * lnE; }
    }
    if ( lnEcor < 0.5 ) continue;
    //    cout<<lnEcor<<" "<<etof<<endl;
 
    //////////////////////////////////////now  no using in the following
    HepPoint3D pos( emcTrk->position() );
 
    // If it is "hot", return "9999" (Hajime, Jul 2013)
    double enecor = 1.;
    if ( hotcellmask == 1 )
    {
      for ( int ih = 0; ih < 10; ih++ )
      {
        if ( hrunstart[ih] == -1 || hrunend[ih] == -1 || hcell[ih] == -1 ) continue;
        if ( abs( runNo ) < hrunstart[ih] || abs( runNo ) > hrunend[ih] ) continue;
        if ( emcTrk->cellId() == hcell[ih] ) { emcTrk->setStatus( 9999 ); }
      }
    }
 
    if ( edgecor == 1 )
    {
 
      if ( module == 1 )
      { // barrel
        unsigned int theModule;
        if ( thetaModule > 21 )
        {
          theModule = 43 - thetaModule;
          id        = EmcID::crystal_id( module, theModule, phiModule );
          pos.setZ( -pos.z() );
        }
        else { theModule = thetaModule; }
 
        double IRShift;
        if ( theModule == 21 ) { IRShift = 0; }
        else if ( theModule == 20 ) { IRShift = 2.5; }
        else { IRShift = 5.; }
 
        HepPoint3D IR( 0, 0, IRShift );
        HepPoint3D center01, center23; // center of point0,1 and point2,3 in crystal
        center01 =
            ( iGeoSvc->GetCrystalPoint( id, 0 ) + iGeoSvc->GetCrystalPoint( id, 1 ) ) / 2;
        center23 =
            ( iGeoSvc->GetCrystalPoint( id, 2 ) + iGeoSvc->GetCrystalPoint( id, 3 ) ) / 2;
 
        double theta01, theta23, length, d;
        theta01 = ( center01 - IR ).theta();
        theta23 = ( center23 - IR ).theta();
        length  = ( center01 - IR ).mag();
        d       = length * tan( theta01 - theta23 ); // length in x direction
 
        double xIn;
        xIn = length * tan( theta01 - ( pos - IR ).theta() ) - d / 2;
        if ( xIn < ( -d / 2 ) && theModule != 21 )
        {
 
          theModule = theModule + 1;
 
          id = EmcID::crystal_id( module, theModule, phiModule );
          if ( theModule == 21 ) { IRShift = 0; }
          else if ( theModule == 20 ) { IRShift = 2.5; }
          else { IRShift = 5.; }
          HepPoint3D IR1( 0, 0, IRShift );
          IR = IR1;
          center01 =
              ( iGeoSvc->GetCrystalPoint( id, 0 ) + iGeoSvc->GetCrystalPoint( id, 1 ) ) / 2;
          center23 =
              ( iGeoSvc->GetCrystalPoint( id, 2 ) + iGeoSvc->GetCrystalPoint( id, 3 ) ) / 2;
          theta01 = ( center01 - IR ).theta();
          theta23 = ( center23 - IR ).theta();
          length  = ( center01 - IR ).mag();
          d       = length * tan( theta01 - theta23 ); // length in x direction
 
          xIn = length * tan( theta01 - ( pos - IR ).theta() ) - d / 2;
        }
        else if ( xIn > ( d / 2 ) && theModule != 0 )
        {
 
          theModule = theModule - 1;
          id        = EmcID::crystal_id( module, theModule, phiModule );
          if ( theModule == 21 ) { IRShift = 0; }
          else if ( theModule == 20 ) { IRShift = 2.5; }
          else { IRShift = 5.; }
 
          HepPoint3D IR1( 0, 0, IRShift );
          IR = IR1;
          center01 =
              ( iGeoSvc->GetCrystalPoint( id, 0 ) + iGeoSvc->GetCrystalPoint( id, 1 ) ) / 2;
          center23 =
              ( iGeoSvc->GetCrystalPoint( id, 2 ) + iGeoSvc->GetCrystalPoint( id, 3 ) ) / 2;
          theta01 = ( center01 - IR ).theta();
          theta23 = ( center23 - IR ).theta();
          length  = ( center01 - IR ).mag();
          d       = length * tan( theta01 - theta23 ); // length in x direction
 
          xIn = length * tan( theta01 - ( pos - IR ).theta() ) - d / 2;
        }
 
        double yIn;
        //    yIn = pos.phi() < 0 ? pos.phi()*180./CLHEP::pi+360.-phiModule*3.-1.1537
        //      : pos.phi()*180./CLHEP::pi-phiModule*3.-1.1537;
 
        yIn = pos.phi() < 0 ? pos.phi() * 180. / CLHEP::pi + 360. - phiModule * 3. - 1.843
                            : pos.phi() * 180. / CLHEP::pi - phiModule * 3. - 1.843;
 
        if ( nphi == 0 && yIn > 100 ) yIn = yIn - 360;
        if ( nphi == 119 && yIn < -100 ) yIn = yIn + 360;
 
        if ( yIn < -1.5 - 0.15 )
        {
 
          if ( nphi != 0 ) phiModule = phiModule - 1;
          if ( nphi == 0 ) phiModule = 119;
          yIn = pos.phi() < 0 ? pos.phi() * 180. / CLHEP::pi + 360. - phiModule * 3. - 1.843
                              : pos.phi() * 180. / CLHEP::pi - phiModule * 3. - 1.843;
          if ( phiModule == 0 && yIn > 100 ) yIn = yIn - 360;
          if ( phiModule == 119 && yIn < -100 ) yIn = yIn + 360;
        }
 
        if ( yIn > 1.5 - 0.15 )
        {
 
          if ( nphi != 119 ) phiModule = phiModule + 1;
          if ( nphi == 119 ) phiModule = 0;
 
          yIn = pos.phi() < 0 ? pos.phi() * 180. / CLHEP::pi + 360. - phiModule * 3. - 1.843
                              : pos.phi() * 180. / CLHEP::pi - phiModule * 3. - 1.843;
          if ( phiModule == 0 && yIn > 100 ) yIn = yIn - 360;
          if ( phiModule == 119 && yIn < -100 ) yIn = yIn + 360;
        }
 
        enecor = m_par[theModule][0] + m_par[theModule][1] * xIn +
                 m_par[theModule][2] * xIn * xIn + m_par[theModule][3] * xIn * xIn * xIn +
                 m_par[theModule][4] * xIn * xIn * xIn * xIn +
                 m_par[theModule][5] * xIn * xIn * xIn * xIn * xIn +
                 m_par[theModule][6] * xIn * xIn * xIn * xIn * xIn * xIn
                 // +m_par[theModule][7]*xIn*xIn*xIn*xIn*xIn*xIn*xIn
                 // +m_par[theModule][8]*xIn*xIn*xIn*xIn*xIn*xIn*xIn*xIn
                 + m_par[theModule][7] * yIn + m_par[theModule][8] * yIn * yIn +
                 m_par[theModule][9] * yIn * yIn * yIn +
                 m_par[theModule][10] * yIn * yIn * yIn * yIn;
        //////////////////////
      }
      else { enecor = 1; }
    }
    //////////////////////////////////////now  no using in the above
 
    double energyCC = energyC / ( lnEcor * enecor );
    //     cout<<"Trk Level Corr. and Orig. "<<energyCC<<" "<<emcTrk->energy()<<endl;
    emcTrk->setEnergy( energyCC );
 
    if ( ntOut == true )
    {
      m_ef    = energyCC;
      m_e5    = emcTrk->e5x5();
      m_ec    = energyC;
      m_ct    = lnEcor;
      m_cedge = enecor;
      m_tuple1->write();
    }
  }
  //==============Modify Dst===============================================================
  SmartDataPtr<RecEmcShowerCol> recEmcShowerCol( eventSvc(),
                                                 EventModel::Recon::RecEmcShowerCol );
  if ( !recEmcShowerCol ) { return ( StatusCode::SUCCESS ); }
  SmartDataPtr<DstEmcShowerCol> dstEmcShowerCol( eventSvc(),
                                                 EventModel::Dst::DstEmcShowerCol );
  if ( !dstEmcShowerCol ) { return ( StatusCode::SUCCESS ); }
 
  //   cout<<"Rec and Dst Size "<<recEmcShowerCol->size()<<" "<<dstEmcShowerCol->size()<<endl;
  if ( recEmcShowerCol->size() != dstEmcShowerCol->size() ) return StatusCode::SUCCESS;
  for ( int i = 0; i < recEmcShowerCol->size(); i++ )
  {
    RecEmcShowerCol::iterator iter_rec = recEmcShowerCol->begin() + i;
    DstEmcShowerCol::iterator iter_dst = dstEmcShowerCol->begin() + i;
    //     cout<<"Rec and Dst energy       "<<(*iter_rec)->energy()<<"
    //     "<<(*iter_dst)->energy()<<endl;
    ( *iter_dst )->setEnergy( ( *iter_rec )->energy() );
    ( *iter_dst )->setStatus( ( *iter_rec )->status() );
    //     cout<<"DST after calib EMC status =  "<<(*iter_dst)->status()<<endl;
    //     cout<<"Rec == Dst?              "<<(*iter_rec)->energy()<<"
    //     "<<(*iter_dst)->energy()<<endl;
  }
  return ( StatusCode::SUCCESS );
}

finalize()

StatusCode AbsCor::finalize

(

)

Definition at line 736 of file AbsCor.cxx.

                            {
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in finalize()" << endmsg;
  return StatusCode::SUCCESS;
}

initialize()

StatusCode AbsCor::initialize

(

)

Definition at line 85 of file AbsCor.cxx.

                              {
  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 );
        status = m_tuple1->addItem( "cedge", m_cedge );
      }
      else
      {
        log << MSG::ERROR << "    Cannot book N-tuple:" << long( m_tuple1 ) << endmsg;
        return StatusCode::FAILURE;
      }
    }
  }
 
  m_index = new int*[56];
  for ( int j = 0; j < 56; j++ )
  {
    m_index[j] = new int[120];
    for ( int k = 0; k < 120; k++ ) { m_index[j][k] = -1; }
  }
 
  m_par = new double*[22];
  for ( int j = 0; j < 22; j++ )
  {
    m_par[j] = new double[11];
    for ( int k = 0; k < 11; k++ ) { m_par[j][k] = -99; }
  }
 
  m_parphi = new double*[22];
  for ( int j = 0; j < 22; j++ )
  {
    m_parphi[j] = new double[5];
    for ( int k = 0; k < 5; k++ ) { m_parphi[j][k] = -99; }
  }
 
  ifstream EnParIn;
  // EnParIn.exceptions( ifstream::failbit | ifstream::badbit );
  string EnParInFile;
  EnParInFile = getenv( "ABSCORROOT" );
  EnParInFile += "/share/para.dat";
  EnParIn.open( EnParInFile.c_str() );
  for ( int i = 0; i < 22; i++ )
  {
 
    EnParIn >> m_par[i][0] >> m_par[i][1] >> m_par[i][2] >> m_par[i][3] >> m_par[i][4] >>
        m_par[i][5] >> m_par[i][6] >> m_par[i][7] >> m_par[i][8] >> m_par[i][9] >>
        m_par[i][10];
  }
  EnParIn.close();
 
  ifstream EnParInphi;
  // EnParInphi.exceptions( ifstream::failbit | ifstream::badbit );
  string EnParInFilephi = getenv( "ABSCORROOT" );
  EnParInFilephi += "/share/paraphi.dat";
  EnParInphi.open( EnParInFilephi.c_str() );
  for ( int i = 0; i < 22; i++ )
  {
    EnParInphi >> m_parphi[i][0] >> m_parphi[i][1] >> m_parphi[i][2] >> m_parphi[i][3] >>
        m_parphi[i][4];
  }
  EnParInphi.close();
 
  /* liucx
  string DataPathevse;
  DataPathevse = getenv("ABSCORROOT");
  DataPathevse += "/share/barmccorevse10bin.txt";
  ifstream inpi0;
  inpi0.exceptions( ifstream::failbit | ifstream::badbit );
  inpi0.open(DataPathevse.c_str(),ios::in);
 
  double epoint[11],peak[11],peakerr[11];
  for(Int_t i=0;i<11;i++){
    inpi0>>epoint[i];
    inpi0>>peak[i];
    inpi0>>peakerr[i];
  }
  for(Int_t i=0;i<11;i++){
    dt->SetPoint(i, epoint[i],peak[i]);
    dt->SetPointError(i,0,peakerr[i]);
  }
  */
 
  string DataPathcbeam;
  DataPathcbeam = getenv( "ABSCORROOT" );
  DataPathcbeam += "/share/cbeam.txt";
  ifstream incbeam;
  incbeam.exceptions( ifstream::failbit | ifstream::badbit );
  incbeam.open( DataPathcbeam.c_str(), ios::in );
  for ( int i = 0; i < 56; i++ )
  {
    double tid, peak, peake, sig, sige;
    incbeam >> tid;
    incbeam >> peak;
    incbeam >> peake;
    incbeam >> sig;
    incbeam >> sige;
    cbeam[i] = 1.0 / peak;
  }
  /* by liucx
  /////////////////////////////
  string DataPathc3p;
  DataPathc3p = getenv("ABSCORROOT");
  DataPathc3p += "/share/c3p.txt";
 
  string DataPathc3ptof;
  DataPathc3ptof = getenv("ABSCORROOT");
  DataPathc3ptof += "/share/c3ptof.txt";
  cout<<"mccor = "<<mccor<<endl;
 //DataPathc3ptof = m_EmcShEnCalibSvc -> getPi0CalibFile();
 
 
  ifstream inc3p;
  inc3p.exceptions( ifstream::failbit | ifstream::badbit );
  if(usetof==0)inc3p.open(DataPathc3p.c_str(),ios::in);
  if(usetof==1)inc3p.open(DataPathc3ptof.c_str(),ios::in);
  for(int i=0; i<4; i++){
    double am,ame;
    inc3p>>am;
    inc3p>>ame;
    ai[i]=am;
  }
  */
  string paraPath = getenv( "ABSCORROOT" );
  if ( paraPath == "" ) {} // cout<<"AbsCor environment not set!";
 
  // Read energy correction parameters from PhotonCor/McCor
  if ( MCuseTof == 1 ) { paraPath += "/share/evsetTof.txt"; }
  if ( MCuseTof == 0 ) { paraPath += "/share/evset.txt"; }
 
  ifstream in2;
  in2.open( paraPath.c_str() );
  assert( in2 );
  double energy, thetaid, peak1, peakerr1, res, reserr;
  dt    = new TGraph2DErrors();
  dtErr = new TGraph2DErrors();
  // for(int i=0;i<560;i++){
  for ( int i = 0; i < 1484; i++ )
  { // 53*28
    in2 >> energy;
    in2 >> thetaid;
    in2 >> peak1;
    in2 >> peakerr1;
    in2 >> res;
    in2 >> reserr;
    dt->SetPoint( i, energy, thetaid, peak1 );
    dt->SetPointError( i, 0, 0, peakerr1 );
    dtErr->SetPoint( i, energy, thetaid, res );
    dtErr->SetPointError( i, 0, 0, reserr );
    if ( i < 28 ) e25min[int( thetaid )] = energy;
    if ( i >= 1484 - 28 ) e25max[int( thetaid )] = energy;
    // if(i>=560-28) e25max[int(thetaid)]=energy;
  }
  in2.close();
 
  //-------------------------
  // Suppression of hot crystals
  // Reading the map from hotcry.txt (Hajime, Jul 2013)
  for ( int ih = 0; ih < 10; ih++ )
  {
    hrunstart[ih] = -1;
    hrunend[ih]   = -1;
    hcell[ih]     = -1;
  }
  int    numhots = 4; // numbers of hot crystals
  int    dumring, dumphi, dummod, dumid;
  string HotList;
  HotList = getenv( "ABSCORROOT" );
  HotList += "/share/hotcry.txt";
  ifstream hotcrys;
  hotcrys.exceptions( ifstream::failbit | ifstream::badbit );
  hotcrys.open( HotList.c_str(), ios::in );
  for ( int il = 0; il < numhots; il++ )
  {
    hotcrys >> hrunstart[il];
    hotcrys >> hrunend[il];
    hotcrys >> hcell[il];
    hotcrys >> dumring;
    hotcrys >> dumphi;
    hotcrys >> dummod;
    hotcrys >> dumid;
  }
  hotcrys.close();
  //-------------------------
  /* by liucx
  ////////////////////to read the parameters of energy correction function//////////
  string CorFunparaPath;
 
  CorFunparaPath = getenv("ABSCORROOT");
 
  // Read energy correction Function parameters from PhotonCor/McCor
  if(MCuseTof==1){
    if (IYear==2009) CorFunparaPath += "/share/evsetTofCorFunctionPar2009Jpsi.txt";
    if (IYear==2012) CorFunparaPath += "/share/evsetTofCorFunctionPar2012Jpsi.txt";
    if (IYear==2018) CorFunparaPath += "/share/evsetTofCorFunctionPar2018Jpsi.txt";
    if (IYear==2019) CorFunparaPath += "/share/evsetTofCorFunctionPar2019Jpsi.txt";
 
  }
  if(MCuseTof==0){
    CorFunparaPath += "/share/evsetCorFunctionPar.txt";
  }
 
 ifstream in2corfun;
 in2corfun.open(CorFunparaPath.c_str());
 assert(in2corfun);
 
  for(int i=0;i<28;i++){
    in2corfun>>m_corFunPar[i][0];
    in2corfun>>m_corFunPar[i][1];
    in2corfun>>m_corFunPar[i][2];
    in2corfun>>m_corFunPar[i][3];
    in2corfun>>m_corFunPar[i][4];
    in2corfun>>m_corFunPar[i][5];
  }
  in2corfun.close();
  ///////////////////////
 
*/
 
  //////// read the shower correction parameters from database////// by liucx
  ISvcLocator* svcLocator = Gaudi::svcLocator();
  StatusCode   sc         = svcLocator->service( "EmcShEnCalibSvc", m_EmcShEnCalibSvc );
 
  if ( sc != StatusCode::SUCCESS )
  { std::cout << "can not use EmcShEnCalibSvc in AbsCor" << endl; }
  else
  {
    std::cout
        << "getPi0CalibFile() and getSingleGammaCalibFile() still are empty in initialize()"
        << std::endl;
    std::cout << "will be assigned a value in execute()" << std::endl;
    std::cout << "Test EmcShEnCalibSvc in AbsCor initialize getPi0CalibFile()= "
              << m_EmcShEnCalibSvc->getPi0CalibFile()
              << "Test EmcShEnCalibSvc in AbsCor getSingleGammaCalibFile()= "
              << m_EmcShEnCalibSvc->getSingleGammaCalibFile() << std::endl;
  }
 
  log << MSG::INFO << "successfully return from initialize()" << endmsg;
  return StatusCode::SUCCESS;
}

---

The documentation for this class was generated from the following files:

• AbsCor.h
• AbsCor.cxx