DedxCalib.cxx

Go to the documentation of this file.

#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/PropertyMgr.h"
#include "GaudiKernel/StatusCode.h"
 
#include "TFile.h"
#include "TTree.h"
#include <fstream>
#include <string>
 
#include "DedxCalibAlg/DedxCalib.h"
// DECLARE_COMPONENT(DedxCalib)
DedxCalib::DedxCalib( const std::string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator ) {
  declareProperty( "ParticleType", ParticleFlag = -1 );
  declareProperty( "EventType", m_eventType = "isBhabha" );
  declareProperty( "CalibFlag", m_calibflag = 63 );
  declareProperty( "PhShapeFlag", m_phShapeFlag = 0 );
  declareProperty( "TruncRate", truncate = 0.7 );
  declareProperty( "RecFileList", m_recFileList = "no file" );
  declareProperty( "RootFile", m_rootfile = "no file" );
  declareProperty( "CurveFile", m_curvefile = std::string( "no rootfile" ) );
}
 
StatusCode DedxCalib::initialize() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalib initialze()" << endmsg;
 
  StatusCode gesc = service( "MdcGeomSvc", geosvc );
  if ( gesc == StatusCode::SUCCESS ) { log << MSG::INFO << "MdcGeomSvc is running" << endmsg; }
  else
  {
    log << MSG::ERROR << "MdcGeomSvc is failed" << endmsg;
    return StatusCode::SUCCESS;
  }
 
  StatusCode scex = service( "DedxCorrecSvc", exsvc );
  if ( scex == StatusCode::SUCCESS )
  { log << MSG::INFO << "Hi, DedxCorrectSvc is running" << endl; }
  else
  {
    log << MSG::ERROR << "DedxCorrectSvc is failed" << endl;
    return StatusCode::SUCCESS;
  }
  exsvc->set_flag( m_calibflag );
 
  initializing();
  BookHists();
  getCurvePar();
 
  return StatusCode::SUCCESS;
}
 
StatusCode DedxCalib::execute() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalib execute()" << endmsg;
 
  genNtuple();
  FillHists();
  AnalyseHists();
  WriteHists();
 
  return StatusCode::SUCCESS;
}
 
StatusCode DedxCalib::finalize() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalib finalize()" << endmsg;
 
  return StatusCode::SUCCESS;
}
 
void DedxCalib::initializing() {}
 
void DedxCalib::genNtuple() {}
 
void DedxCalib::BookHists() {}
 
void DedxCalib::FillHists() {}
 
void DedxCalib::AnalyseHists() {}
 
void DedxCalib::WriteHists() {}
 
void DedxCalib::ReadRecFileList() {
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalib::ReadRecFileList()" << endmsg;
 
  ifstream filea( m_recFileList.c_str(), ifstream::in );
  cout << m_recFileList.c_str() << endl;
 
  string runlist;
  filea >> runlist;
  do {
    m_recFileVector.push_back( runlist );
    cout << runlist.c_str() << endl;
    filea >> runlist;
  } while ( filea );
}
 
double DedxCalib::cal_dedx_bitrunc( float truncate, std::vector<double> phlist ) {
  sort( phlist.begin(), phlist.end() );
  int      smpl    = (int)( phlist.size() * ( truncate + 0.05 ) );
  int      min_cut = (int)( phlist.size() * 0.05 + 0.5 );
  double   qSum    = 0;
  unsigned i       = 0;
  for ( vector<double>::iterator ql = phlist.begin(); ql != phlist.end(); ql++ )
  {
    i++;
    if ( i <= smpl && i >= min_cut ) qSum += ( *ql );
  }
  double trunc   = -99;
  int    usedhit = smpl - min_cut + 1;
  if ( usedhit > 0 ) trunc = qSum / usedhit;
 
  return trunc;
}
 
double DedxCalib::cal_dedx( float truncate, std::vector<double> phlist ) {
  sort( phlist.begin(), phlist.end() );
  int      smpl    = (int)( phlist.size() * ( truncate + 0.05 ) );
  int      min_cut = 0;
  double   qSum    = 0;
  unsigned i       = 0;
  for ( vector<double>::iterator ql = phlist.begin(); ql != phlist.end(); ql++ )
  {
    i++;
    if ( i <= smpl && i >= min_cut ) qSum += ( *ql );
  }
  double trunc   = -99;
  int    usedhit = smpl - min_cut + 1;
  if ( usedhit > 0 ) trunc = qSum / usedhit;
 
  return trunc;
}
 
void DedxCalib::getCurvePar()
// void DedxCalib::getCurvePar(int runflag)
{
  if ( m_curvefile == "no rootfile" )
  {
    cout << "no curve file! can not calculate chi!!! " << endl;
    return;
  }
 
  double Curve[100];
  double Sigma[100];
  int    CurveSize, SigmaSize;
  TFile* f     = new TFile( m_curvefile.c_str() );
  TTree* curve = (TTree*)f->Get( "curve" );
  TTree* sigma = (TTree*)f->Get( "sigma" );
  curve->SetBranchAddress( "CurveSize", &CurveSize );
  curve->SetBranchAddress( "curve", Curve );
  sigma->SetBranchAddress( "SigmaSize", &SigmaSize );
  sigma->SetBranchAddress( "sigma", Sigma );
  curve->GetEntry( 0 );
  sigma->GetEntry( 0 );
 
  for ( int i = 0; i < CurveSize; i++ ) // get the dedx curve parameters
  {
    cout << Curve[i] << endl;
    if ( i == 0 ) vFlag[0] = (int)Curve[i]; // first element is dedx curve version
    else Curve_Para.push_back( Curve[i] );  // dedx curve parameters
  }
  for ( int i = 0; i < SigmaSize; i++ )
  {
    cout << Sigma[i] << endl;
    if ( i == 0 ) vFlag[1] = (int)Sigma[i]; // dedx sigma version: 2: psip; 3:jpsi
    else Sigma_Para.push_back( Sigma[i] );  // dedx sigma parameters
  }
  // if(runflag>0)   vFlag[2]=0;
  // else            vFlag[2]=1;
}
 
void DedxCalib::set_dEdx( int landau, float dEdx, int trkalg, int runflag, int dedxhit_use,
                          float ptrk, float theta, float pl_rp, int vflag[3], double t0 ) {
  // landau: 1:landau distribution; 0:gaus distribution
  // pl_rp: 1.5, not know what
  // getCurvePar(runflag);
  if ( runflag > 0 ) vFlag[2] = 0;
  else vFlag[2] = 1;
 
  // some old data with old methods
  if ( runflag == 1 || runflag == 2 )
    dedx_pid_exp_old( landau, runflag, (float)dEdx, (int)dedxhit_use, (float)ptrk,
                      (float)theta, (float)t0, (float)pl_rp, m_dedx_exp, m_ex_sigma,
                      m_pid_prob, m_chi_ex );
  // for 2009 psip data and after
  else
    dedx_pid_exp( vflag, (float)dEdx, trkalg, (int)dedxhit_use, (float)ptrk, (float)theta,
                  (float)t0, (float)pl_rp, m_dedx_exp, m_ex_sigma, m_pid_prob, m_chi_ex,
                  Curve_Para, Sigma_Para );
}