DedxCalibRunByRun.cxx

Go to the documentation of this file.

#include "GaudiKernel/MsgStream.h"
 
#include "TF1.h"
#include "TFile.h"
#include "TH1F.h"
#include "TTree.h"
#include <sstream>
#include <string>
 
#include "DedxCalibAlg/DedxCalibRunByRun.h"
DECLARE_COMPONENT( DedxCalibRunByRun )
using namespace std;
 
int runNo = 0;
 
DedxCalibRunByRun::DedxCalibRunByRun( const std::string& name, ISvcLocator* pSvcLocator )
    : DedxCalib( name, pSvcLocator ) {}
 
void DedxCalibRunByRun::BookHists() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalibRunByRun::BookHists()" << endmsg;
 
  ReadRecFileList();
  runNo = m_recFileVector.size();
 
  m_hist = new TH1F( "dEdx_allRun", "dEdx_allRun", NumHistBins, MinHistValue1, MaxHistValue1 );
  m_hist->GetXaxis()->SetTitle( "dE/dx" );
  m_hist->GetYaxis()->SetTitle( "entries" );
  m_rungain = new TH1F*[runNo * NumEvtBlocksInRun];
  stringstream hlabel;
  for ( int i = 0; i < runNo; i++ )
  {
    for ( int j = 0; j < NumEvtBlocksInRun; j++ )
    {
      hlabel.str( "" );
      hlabel << "dEdx_Run_" << i << "_sub_" << j;
      m_rungain[i * NumEvtBlocksInRun + j] =
          new TH1F( hlabel.str().c_str(), hlabel.str().c_str(), NumHistBins, MinHistValue1,
                    MaxHistValue1 );
    }
  }
}
 
void DedxCalibRunByRun::FillHists() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalibRunByRun::FillHists()" << endmsg;
 
  TFile* f;
  TTree* n103;
  string runlist;
 
  int    ndedxhit  = 0;
  double dEdx[100] = { 0 }, pathlength[100] = { 0 }, wid[100] = { 0 }, layid[100] = { 0 },
         dd_in[100] = { 0 }, eangle[100] = { 0 }, zhit[100] = { 0 };
  double         dedx  = 0;
  float          runNO = 0, evtNO = 0, runFlag = 0, costheta = 0, tes = 0;
  vector<double> phlist;
  for ( int i = 0; i < runNo; i++ )
  {
    runlist = m_recFileVector[i];
    f       = new TFile( runlist.c_str() );
    n103    = (TTree*)f->Get( "n103" );
    n103->SetBranchAddress( "ndedxhit", &ndedxhit );
    n103->SetBranchAddress( "dEdx_hit", dEdx );
    n103->SetBranchAddress( "pathlength_hit", pathlength );
    n103->SetBranchAddress( "wid_hit", wid );
    n103->SetBranchAddress( "layid_hit", layid );
    n103->SetBranchAddress( "dd_in_hit", dd_in );
    n103->SetBranchAddress( "eangle_hit", eangle );
    n103->SetBranchAddress( "zhit_hit", zhit );
    n103->SetBranchAddress( "runNO", &runNO );
    n103->SetBranchAddress( "evtNO", &evtNO );
    n103->SetBranchAddress( "runFlag", &runFlag );
    n103->SetBranchAddress( "costheta", &costheta );
    n103->SetBranchAddress( "t0", &tes );
    int           m( 0 ); // count for events in a block
    int           index( 0 );
    float         evtfrom( 0 );
    struct runevt m_runevt;
    for ( int j = 0; j < n103->GetEntries(); j++ )
    {
      phlist.clear();
      n103->GetEntry( j );
      for ( int k = 0; k < ndedxhit; k++ )
      {
        dEdx[k] = exsvc->StandardHitCorrec( 0, runFlag, 2, runNO, evtNO, pathlength[k], wid[k],
                                            layid[k], dEdx[k], dd_in[k], eangle[k], costheta );
        phlist.push_back( dEdx[k] );
      }
      dedx = cal_dedx_bitrunc( truncate, phlist );
      dedx = exsvc->StandardTrackCorrec( 0, runFlag, 2, runNO, evtNO, dedx, costheta, tes );
      m_rungain[i * NumEvtBlocksInRun + index]->Fill( dedx );
      m_hist->Fill( dedx );
      m++;
      if ( m == NumEvtInBlocks )
      {
        m_runevt.runno        = runNO;
        m_runevt.evtfrom      = evtfrom;
        m_runevt.evtto        = evtNO;
        m_runevt.global_index = i * NumEvtBlocksInRun + index;
        m_runevt.local_index  = index;
        m_evtNOVector.push_back( m_runevt );
        evtfrom = evtNO + 1;
        index++;
        m = 0;
      }
      if ( index == NumEvtBlocksInRun )
      {
        log << MSG::FATAL << "DedxCalibRunByRun: Too many events in a run!!!" << endmsg;
        break;
      }
    }
    if ( index ) m_evtNOVector.pop_back();
    else
    {
      m_runevt.runno        = runNO;
      m_runevt.evtfrom      = 0;
      m_runevt.global_index = i * NumEvtBlocksInRun;
      m_runevt.local_index  = 0;
    }
    m_runevt.evtto = MaxNumEvt;
    m_evtNOVector.push_back( m_runevt );
    cout << "runNO = " << m_runevt.runno << endl;
  }
}
 
void DedxCalibRunByRun::AnalyseHists() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalibRunByRun::AnalyseHists()" << endmsg;
 
  m_hist->Fit( "gaus", "Q" );
  for ( int i = 0; i < runNo; i++ )
  {
    for ( int j = 0; j < NumEvtBlocksInRun; j++ )
    { m_rungain[i * NumEvtBlocksInRun + j]->Fit( "gaus", "Q" ); }
  }
}
 
void DedxCalibRunByRun::WriteHists() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalibRunByRun::WriteHists()" << endmsg;
 
  TFile* f = new TFile( m_rootfile.c_str(), "recreate" );
 
  Double_t gainpar = 0, resolpar = 0;
  gainpar  = m_hist->GetFunction( "gaus" )->GetParameter( 1 );
  resolpar = m_hist->GetFunction( "gaus" )->GetParameter( 2 );
 
  TTree* gain = new TTree( "gaincalib", "gaincalib" );
  gain->Branch( "gain", &gainpar, "gain[1]/D" );
  gain->Fill();
  TTree* resol = new TTree( "resolcalib", "resolcalib" );
  resol->Branch( "resol", &resolpar, "resol[1]/D" );
  resol->Fill();
 
  Int_t    runno( 0 ), evtfrom( 0 ), evtto( 0 ), index( 0 ), global_index( 0 );
  Double_t runmean = 0, rungain = 0, runresol = 0;
  TTree*   runbyrun = new TTree( "runcalib", "runcalib" );
  runbyrun->Branch( "runno", &runno, "runno/I" );
  runbyrun->Branch( "index", &index, "index/I" );
  runbyrun->Branch( "evtfrom", &evtfrom, "evtfrom/I" );
  runbyrun->Branch( "evtto", &evtto, "evtto/I" );
  runbyrun->Branch( "runmean", &runmean, "runmean/D" );
  runbyrun->Branch( "rungain", &rungain, "rungain/D" );
  runbyrun->Branch( "runresol", &runresol, "runresol/D" );
  vector<struct runevt>::iterator p;
  p = m_evtNOVector.begin();
  while ( p != m_evtNOVector.end() )
  {
    runno        = ( *p ).runno;
    evtfrom      = ( *p ).evtfrom;
    evtto        = ( *p ).evtto;
    index        = ( *p ).local_index;
    global_index = ( *p ).global_index;
    runmean      = m_rungain[global_index]->GetFunction( "gaus" )->GetParameter( 1 );
    runresol     = m_rungain[global_index]->GetFunction( "gaus" )->GetParameter( 2 );
    rungain      = runmean / NormalMean;
    cout << "   runno = " << runno << "   @  evtfrom = " << evtfrom << "  @ evtto " << evtto
         << "  @ index " << index << "  @  runmean = " << runmean
         << "  @  rungain = " << rungain << "  @  runresol =  " << runresol << endl;
    runbyrun->Fill();
    m_rungain[global_index]->Write();
    p++;
  }
 
  m_hist->Write();
  gain->Write();
  resol->Write();
  runbyrun->Write();
 
  delete m_hist;
  for ( int i = 0; i < runNo; i++ )
  {
    for ( int j = 0; j < NumEvtBlocksInRun; j++ ) delete m_rungain[i * NumEvtBlocksInRun + j];
  }
  delete gain;
  delete resol;
  delete runbyrun;
  f->Close();
}