DedxCalibTzero.cxx

Go to the documentation of this file.

#include "GaudiKernel/MsgStream.h"
 
#include "TCanvas.h"
#include "TF1.h"
#include "TFile.h"
#include "TH1F.h"
#include "TStyle.h"
#include "TTree.h"
#include <sstream>
#include <string>
 
#include "DedxCalibAlg/DedxCalibTzero.h"
 
#define Size 700000
#define numcut 50
DECLARE_COMPONENT( DedxCalibTzero )
using namespace std;
 
DedxCalibTzero::DedxCalibTzero( const std::string& name, ISvcLocator* pSvcLocator )
    : DedxCalib( name, pSvcLocator ) {
  declareProperty( "Debug", m_debug = false );
}
 
void DedxCalibTzero::BookHists() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalibTzero::BookHists()" << endmsg;
 
  ReadRecFileList();
 
  m_tzero     = new TH1F*[NumBinCalT0];
  m_pos_tzero = new TH1F*[NumBinCalT0];
  m_neg_tzero = new TH1F*[NumBinCalT0];
  m_chi       = new TH1F*[NumBinCalT0];
  m_pos_chi   = new TH1F*[NumBinCalT0];
  m_neg_chi   = new TH1F*[NumBinCalT0];
  stringstream hlabel;
  for ( int i = 0; i < NumBinCalT0; i++ )
  {
    hlabel.str( "" );
    hlabel << "dEdx_tzero" << i;
    m_tzero[i] = new TH1F( hlabel.str().c_str(), hlabel.str().c_str(), NumHistBins,
                           MinHistValue, MaxHistValue );
    hlabel.str( "" );
    hlabel << "pos_dEdx_tzero" << i;
    m_pos_tzero[i] = new TH1F( hlabel.str().c_str(), hlabel.str().c_str(), NumHistBins,
                               MinHistValue, MaxHistValue );
    hlabel.str( "" );
    hlabel << "neg_dEdx_tzero" << i;
    m_neg_tzero[i] = new TH1F( hlabel.str().c_str(), hlabel.str().c_str(), NumHistBins,
                               MinHistValue, MaxHistValue );
    hlabel.str( "" );
    hlabel << "chi_tzero" << i;
    m_chi[i] = new TH1F( hlabel.str().c_str(), hlabel.str().c_str(), NumHistBins, MinChiValue,
                         MaxChiValue );
    hlabel.str( "" );
    hlabel << "pos_chi_tzero" << i;
    m_pos_chi[i] = new TH1F( hlabel.str().c_str(), hlabel.str().c_str(), NumHistBins,
                             MinChiValue, MaxChiValue );
    hlabel.str( "" );
    hlabel << "neg_chi_tzero" << i;
    m_neg_chi[i] = new TH1F( hlabel.str().c_str(), hlabel.str().c_str(), NumHistBins,
                             MinChiValue, MaxChiValue );
  }
  hlabel.str( "" );
  hlabel << "dEdxVsTzero";
  m_dEdxVsTzero = new TH2F( hlabel.str().c_str(), "dEdx vs t0", NumBinTzero, TzeroMin,
                            TzeroMax, NumHistBins, MinHistValue, MaxHistValue );
  m_dEdxVsTzero->GetXaxis()->SetTitle( "t0" );
  m_dEdxVsTzero->GetYaxis()->SetTitle( "dE/dx" );
  hlabel.str( "" );
  hlabel << "pos_dEdxVsTzero";
  m_pos_dEdxVsTzero = new TH2F( hlabel.str().c_str(), "positron dEdx vs t0", NumBinTzero,
                                TzeroMin, TzeroMax, NumHistBins, MinHistValue, MaxHistValue );
  m_pos_dEdxVsTzero->GetXaxis()->SetTitle( "pos t0" );
  m_pos_dEdxVsTzero->GetYaxis()->SetTitle( "dE/dx" );
  hlabel.str( "" );
  hlabel << "neg_dEdxVsTzero";
  m_neg_dEdxVsTzero = new TH2F( hlabel.str().c_str(), "electron dEdx vs tzero", NumBinTzero,
                                TzeroMin, TzeroMax, NumHistBins, MinHistValue, MaxHistValue );
  m_neg_dEdxVsTzero->GetXaxis()->SetTitle( "neg t0" );
  m_neg_dEdxVsTzero->GetYaxis()->SetTitle( "dE/dx" );
 
  Vec_dedx.clear();
  Vec_tzero.clear();
}
 
void DedxCalibTzero::FillHists() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalibTzero::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 },
         tzero_bin_low[NumBinCalT0] = { 0 }, tzero_bin_high[NumBinCalT0] = { 0 };
  double dedx  = 0;
  float  runNO = 0, evtNO = 0, runFlag = 0, costheta = 0, tes = 0, charge = 0, recalg = 0,
        ptrk = 0, chidedx = 0;
  int usedhit = 0;
  for ( int i = 0; i < NumBinCalT0; i++ )
  {
    tzero_bin_low[i]  = ( i == 0 ) ? 0 : ( tzero_bin_value[i - 1] + tzero_bin_value[i] ) / 2;
    tzero_bin_high[i] = ( i == NumBinCalT0 - 1 )
                            ? TzeroMax
                            : ( tzero_bin_value[i] + tzero_bin_value[i + 1] ) / 2;
  }
  vector<double> phlist;
  cout << "m_recFileVector.size()= " << m_recFileVector.size() << endl;
  for ( unsigned int i = 0; i < m_recFileVector.size(); 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 );
    n103->SetBranchAddress( "charge", &charge );
    n103->SetBranchAddress( "recalg", &recalg );
    n103->SetBranchAddress( "ndedxhit", &ndedxhit );
    n103->SetBranchAddress( "ptrk", &ptrk );
    n103->SetBranchAddress( "chidedx_e", &chidedx );
    for ( int j = 0; j < n103->GetEntries(); j++ )
    {
      phlist.clear();
      n103->GetEntry( j );
      if ( costheta > CosthetaMax || costheta < CosthetaMin ) continue;
      if ( tes > TzeroMax ) continue;
      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 );
      int iTes( 0 );
      for ( ; iTes < NumBinCalT0; iTes++ )
        if ( tes > tzero_bin_low[iTes] && tes <= tzero_bin_high[iTes] ) break;
      if ( m_debug && i < 10 && j < 100 )
        cout << "t0 " << tes << "  iTes " << iTes << "  dedx " << dedx << endl;
      // if(m_debug) cout << "before cor, dedx " << pre_dedx << "  with t0 " << tes << endl;
      dedx = exsvc->StandardTrackCorrec( 0, runFlag, 2, runNO, evtNO, dedx, costheta, tes );
      // if(m_debug) cout << "after cor, dedx " << dedx << " with gain " << pre_dedx/dedx <<
      // endl;
      m_tzero[iTes]->Fill( dedx );
      m_dEdxVsTzero->Fill( tes, dedx );
      if ( charge > 0 )
      {
        m_pos_tzero[iTes]->Fill( dedx );
        m_pos_dEdxVsTzero->Fill( tes, dedx );
      }
      if ( charge < 0 )
      {
        m_neg_tzero[iTes]->Fill( dedx );
        m_neg_dEdxVsTzero->Fill( tes, dedx );
      }
 
      usedhit = ndedxhit * truncate;
      set_dEdx( 1, dedx, recalg, runFlag, usedhit, ptrk, acos( costheta ), 1.5, vFlag, tes );
      double chi = m_chi_ex[ParticleFlag];
      // cout<<"chidedx= "<<chidedx<<" particleType= "<<ParticleFlag<<"  new chi= "<<chi<<endl;
      m_chi[iTes]->Fill( chi );
      if ( charge > 0 ) m_pos_chi[iTes]->Fill( chi );
      if ( charge < 0 ) m_neg_chi[iTes]->Fill( chi );
 
      Vec_dedx.push_back( dedx );
      Vec_tzero.push_back( tes );
    }
  }
}
 
void DedxCalibTzero::AnalyseHists() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalibTzero::AnalyseHists()" << endmsg;
 
  gStyle->SetOptFit( 1111 );
  for ( int i = 0; i < NumBinCalT0; i++ )
  {
    if ( m_debug ) cout << "num of bin " << i << endl;
    if ( m_tzero[i]->Integral() > numcut ) m_tzero[i]->Fit( "gaus", "Q" );
    if ( m_pos_tzero[i]->Integral() > numcut ) m_pos_tzero[i]->Fit( "gaus", "Q" );
    if ( m_neg_tzero[i]->Integral() > numcut ) m_neg_tzero[i]->Fit( "gaus", "Q" );
    if ( m_chi[i]->Integral() > numcut ) m_chi[i]->Fit( "gaus", "Q" );
    if ( m_pos_chi[i]->Integral() > numcut ) m_pos_chi[i]->Fit( "gaus", "Q" );
    if ( m_neg_chi[i]->Integral() > numcut ) m_neg_chi[i]->Fit( "gaus", "Q" );
  }
}
 
void DedxCalibTzero::WriteHists() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalibTzero::WriteHists()" << endmsg;
 
  double fitentryNo[NumBinCalT0] = { 0 }, fitmean[NumBinCalT0] = { 0 },
         fitmeanerr[NumBinCalT0] = { 0 }, fitsigma[NumBinCalT0] = { 0 },
         gain[NumBinCalT0] = { 0 }, fitchisq[NumBinCalT0] = { 0 }, t0[NumBinCalT0] = { 0 },
         fakegain[NumBinCalT0]       = { 0 };
  double pos_fitentryNo[NumBinCalT0] = { 0 }, pos_fitmean[NumBinCalT0] = { 0 },
         pos_fitmeanerr[NumBinCalT0] = { 0 }, pos_fitsigma[NumBinCalT0] = { 0 },
         pos_gain[NumBinCalT0] = { 0 }, pos_fitchisq[NumBinCalT0] = { 0 };
  double neg_fitentryNo[NumBinCalT0] = { 0 }, neg_fitmean[NumBinCalT0] = { 0 },
         neg_fitmeanerr[NumBinCalT0] = { 0 }, neg_fitsigma[NumBinCalT0] = { 0 },
         neg_gain[NumBinCalT0] = { 0 }, neg_fitchisq[NumBinCalT0] = { 0 };
 
  //    double avegain(m_dEdxVsTzero->ProjectionX()->GetMean());
  double avegain( 550. );
 
  for ( int i = 0; i < NumBinCalT0; i++ )
  {
    t0[i]         = tzero_bin_value[i];
    fitentryNo[i] = m_tzero[i]->Integral();
    if ( m_debug ) cout << "num of bin " << i << "  with events " << fitentryNo[i] << endl;
    if ( fitentryNo[i] <= numcut ) continue;
    if ( m_debug )
      cout << "get results at " << i << " bin with fit entries " << fitentryNo[i] << endl;
    if ( m_tzero[i]->GetFunction( "gaus" ) )
    {
      fitmean[i]    = m_tzero[i]->GetFunction( "gaus" )->GetParameter( 1 );
      fitmeanerr[i] = m_tzero[i]->GetFunction( "gaus" )->GetParError( 1 );
      fitsigma[i]   = m_tzero[i]->GetFunction( "gaus" )->GetParameter( 2 );
      gain[i]       = fitmean[i] / avegain;
      fakegain[i]   = 1.0;
      fitchisq[i]   = ( m_tzero[i]->GetFunction( "gaus" )->GetChisquare() ) /
                    ( m_tzero[i]->GetFunction( "gaus" )->GetNDF() );
    }
    pos_fitentryNo[i] = m_pos_tzero[i]->Integral();
    if ( pos_fitentryNo[i] <= numcut ) continue;
    if ( m_debug )
      cout << "get results at " << i << " bin with pos_fit entries " << pos_fitentryNo[i]
           << endl;
    if ( m_pos_tzero[i]->GetFunction( "gaus" ) )
    {
      pos_fitmean[i]    = m_pos_tzero[i]->GetFunction( "gaus" )->GetParameter( 1 );
      pos_fitmeanerr[i] = m_pos_tzero[i]->GetFunction( "gaus" )->GetParError( 1 );
      pos_fitsigma[i]   = m_pos_tzero[i]->GetFunction( "gaus" )->GetParameter( 2 );
      pos_gain[i]       = pos_fitmean[i] / avegain;
      pos_fitchisq[i]   = ( m_pos_tzero[i]->GetFunction( "gaus" )->GetChisquare() ) /
                        ( m_pos_tzero[i]->GetFunction( "gaus" )->GetNDF() );
    }
    neg_fitentryNo[i] = m_neg_tzero[i]->Integral();
    if ( neg_fitentryNo[i] <= numcut ) continue;
    if ( m_debug )
      cout << "get results at " << i << " bin with neg_fit entries " << neg_fitentryNo[i]
           << endl;
    if ( m_neg_tzero[i]->GetFunction( "gaus" ) )
    {
      neg_fitmean[i]    = m_neg_tzero[i]->GetFunction( "gaus" )->GetParameter( 1 );
      neg_fitmeanerr[i] = m_neg_tzero[i]->GetFunction( "gaus" )->GetParError( 1 );
      neg_fitsigma[i]   = m_neg_tzero[i]->GetFunction( "gaus" )->GetParameter( 2 );
      neg_gain[i]       = neg_fitmean[i] / avegain;
      neg_fitchisq[i]   = ( m_neg_tzero[i]->GetFunction( "gaus" )->GetChisquare() ) /
                        ( m_neg_tzero[i]->GetFunction( "gaus" )->GetNDF() );
    }
  }
 
  log << MSG::INFO << "begin generating root file!!! " << endmsg;
  TFile* f = new TFile( m_rootfile.c_str(), "RECREATE" );
  for ( int i = 0; i < NumBinCalT0; i++ )
  {
    m_chi[i]->Write();
    m_pos_chi[i]->Write();
    m_neg_chi[i]->Write();
    m_tzero[i]->Write();
    m_pos_tzero[i]->Write();
    m_neg_tzero[i]->Write();
  }
  m_dEdxVsTzero->Write();
  m_pos_dEdxVsTzero->Write();
  m_neg_dEdxVsTzero->Write();
 
  TTree* tzerocalib = new TTree( "gaincalib", "gaincalib" );
  tzerocalib->Branch( "gain", &avegain, "gain/D" );
  tzerocalib->Branch( "t0", t0, "t0[35]/D" );
  tzerocalib->Branch( "dedx", gain, "dedx[35]/D" );
  tzerocalib->Branch( "fakededx", fakegain, "fakededx[35]/D" );
  tzerocalib->Branch( "fitmean", fitmean, "fitmean[35]/D" );
  tzerocalib->Branch( "fitmeanerr", fitmeanerr, "fitmeanerr[35]/D" );
  tzerocalib->Branch( "fitsigma", fitsigma, "fitsigma[35]/D" );
  tzerocalib->Branch( "fitchisq", fitchisq, "fitchisq[35]/D" );
  tzerocalib->Fill();
  tzerocalib->Write();
 
  TCanvas c1( "c1", "canvas", 500, 400 );
  c1.Print( ( m_rootfile + ".ps[" ).c_str() );
  tzerocalib->Draw( "dedx:t0" );
  f->Close(); // delay this close for drawing plots
 
  c1.Print( ( m_rootfile + ".ps" ).c_str() );
  m_dEdxVsTzero->Draw();
  c1.Print( ( m_rootfile + ".ps" ).c_str() );
  m_pos_dEdxVsTzero->Draw();
  c1.Print( ( m_rootfile + ".ps" ).c_str() );
  m_neg_dEdxVsTzero->Draw();
  c1.Print( ( m_rootfile + ".ps" ).c_str() );
  for ( Int_t i = 0; i < NumBinCalT0; i++ )
  {
    m_tzero[i]->Draw();
    c1.Print( ( m_rootfile + ".ps" ).c_str() );
  }
  for ( Int_t i = 0; i < NumBinCalT0; i++ )
  {
    m_pos_tzero[i]->Draw();
    c1.Print( ( m_rootfile + ".ps" ).c_str() );
  }
  for ( Int_t i = 0; i < NumBinCalT0; i++ )
  {
    m_neg_tzero[i]->Draw();
    c1.Print( ( m_rootfile + ".ps" ).c_str() );
  }
  c1.Print( ( m_rootfile + ".ps]" ).c_str() );
  if ( m_debug ) cout << "before delete" << endl;
  for ( Int_t i = 0; i < NumBinCalT0; i++ )
  {
    delete m_chi[i];
    delete m_pos_chi[i];
    delete m_neg_chi[i];
    delete m_tzero[i];
    delete m_pos_tzero[i];
    delete m_neg_tzero[i];
  }
  delete m_dEdxVsTzero;
  delete m_pos_dEdxVsTzero;
  delete m_neg_dEdxVsTzero;
}