DedxCalibDocaEAng.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 "TH2F.h"
#include "TStyle.h"
#include "TTree.h"
#include <sstream>
#include <string>
 
#include "DedxCalibAlg/DedxCalibDocaEAng.h"
DECLARE_COMPONENT( DedxCalibDocaEAng )
 
using namespace std;
const double StepDoca = ( DocaMax - DocaMin ) / NumSlicesDoca;
 
DedxCalibDocaEAng::DedxCalibDocaEAng( const std::string& name, ISvcLocator* pSvcLocator )
    : DedxCalib( name, pSvcLocator ) {
  declareProperty( "Debug", m_debug = false );
  declareProperty( "Idoca", m_idoca_test = 1 );
  declareProperty( "Ieang", m_ieang_test = 1 );
  declareProperty( "PMin", pMin = 0.2 );
  declareProperty( "PMax", pMax = 2.3 );
  declareProperty( "Truncate", m_truncate = 0.7 );
}
 
void DedxCalibDocaEAng::BookHists() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalibDocaEAng::BookHists()" << endmsg;
 
  ReadRecFileList();
 
  m_docaeangle = new TH1F**[NumSlicesDoca];
  stringstream hlabel;
  for ( int i = 0; i < NumSlicesDoca; i++ )
  {
    m_docaeangle[i] = new TH1F*[NumSlicesEAng];
    for ( int j = 0; j < NumSlicesEAng; j++ )
    {
      hlabel.str( "" );
      hlabel << "dEdx_doca" << i << "_eang" << j;
      m_docaeangle[i][j] = new TH1F( hlabel.str().c_str(), hlabel.str().c_str(), NumHistBins,
                                     MinHistValue, MaxHistValue );
    }
  }
  hlabel.str( "" );
  hlabel << "dEdxVsDocaEAng";
  m_DocaEAngAverdE = new TH2F( hlabel.str().c_str(), "dEdx vs Doca & EAng", NumSlicesDoca,
                               DocaMin, DocaMax, NumSlicesEAng, PhiMin, PhiMax );
}
 
void DedxCalibDocaEAng::FillHists() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalibDocaEAng::FillHists()" << endmsg;
 
  TFile *f, *f_test( 0 );
  TTree *n103, *t_test( 0 );
  string runlist;
 
  double doca_test( 0 ), eang_test( 0 );
  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 };
  float runNO = 0, evtNO = 0, runFlag = 0, costheta = 0, tes = 0, charge = 0, recalg = 0,
        ptrk = 0, chidedx = 0;
  cout << "m_recFileVector.size()= " << m_recFileVector.size() << endl;
  if ( m_debug )
  {
    f_test = new TFile( "doca_eangle_test.root", "recreate" );
    t_test = new TTree( "inftest", "information for test doca eangle calibration" );
    t_test->Branch( "idoca", &m_idoca_test, "idoca/I" );
    t_test->Branch( "ieang", &m_ieang_test, "ieang/I" );
    t_test->Branch( "doca", &doca_test, "doca/D" );
    t_test->Branch( "eang", &eang_test, "eang/D" );
  }
 
  for ( unsigned int i = 0; i < m_recFileVector.size(); i++ )
  {
    runlist = m_recFileVector[i];
    cout << "runlist: " << runlist.c_str() << endl;
    f    = new TFile( runlist.c_str() );
    n103 = (TTree*)f->Get( "n103" );
    if ( f == 0 )
    {
      cout << "the file is empty" << endl;
      continue;
    }
    if ( n103 == 0 )
    {
      cout << "the tree is empty" << endl;
      continue;
    }
    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( "ptrk", &ptrk );
    n103->SetBranchAddress( "chidedx_e", &chidedx );
 
    vector<double> m_phlist_test( 0 );
    for ( int j = 0; j < n103->GetEntries(); j++ )
    {
      n103->GetEntry( j );
      if ( ptrk > pMax || ptrk < pMin ) continue;
      if ( tes > 1400 ) continue;
      // begin of save truncated hits information for EACH track
      // m_phlist_test.clear();
      // for(int ktest=0; ktest<ndedxhit; ktest++){
      //   if(layid[ktest]<4) continue;
      //   m_phlist_test.push_back(dEdx[ktest]);
      // }
      // sort(m_phlist_test.begin(), m_phlist_test.end());
      // if(m_debug){
      //   for(unsigned int i=0; i<m_phlist_test.size(); i++)
      //     cout << "i " << i << " " << m_phlist_test[i] << endl;
      // }
      // unsigned int nsample = (unsigned int)(m_phlist_test.size()*m_truncate);
      // double m_trunc_value = m_phlist_test[nsample];
      // if(m_debug) cout << "trunctate value " << m_trunc_value << endl;
      for ( int k = 0; k < ndedxhit; k++ )
      {
        //  if(dEdx[k]>m_trunc_value) continue;
        dd_in[k] = dd_in[k] / doca_norm[int( layid[k] )];
        if ( eangle[k] > 0.25 ) eangle[k] = eangle[k] - 0.5;
        else if ( eangle[k] < -0.25 ) eangle[k] = eangle[k] + 0.5;
        if ( eangle[k] > PhiMin && eangle[k] < PhiMax && dd_in[k] > DocaMin &&
             dd_in[k] < DocaMax )
        {
          if ( layid[k] < 4 ) continue;
          //          else if(layid[k] <8)
          //        {if(dEdx[k]<MinADCValuecut || dEdx[k]>MaxADCValuecut ||
          //pathlength[k]>RPhi_PathMaxCut || pathlength[k]<Iner_RPhi_PathMinCut ||
          //driftdist[k]>Iner_DriftDistCut) continue;}        else      {if(dEdx[k]<MinADCValuecut ||
          //dEdx[k]>MaxADCValuecut || pathlength[k]>RPhi_PathMaxCut ||
          //pathlength[k]<Out_RPhi_PathMinCut || driftdist[k]>Out_DriftDistCut) continue;}
          int iDoca = (Int_t)floor( ( dd_in[k] - DocaMin ) / StepDoca );
          int iEAng = 0;
          for ( int i = 0; i < 14; i++ )
          {
            if ( eangle[k] < Eangle_value_cut[i] || eangle[k] >= Eangle_value_cut[i + 1] )
              continue;
            else if ( eangle[k] >= Eangle_value_cut[i] && eangle[k] < Eangle_value_cut[i + 1] )
            {
              for ( int m = 0; m < i; m++ ) { iEAng += Eangle_cut_bin[m]; }
              double eangle_bin_cut_step =
                  ( Eangle_value_cut[i + 1] - Eangle_value_cut[i] ) / Eangle_cut_bin[i];
              int temp_bin = int( ( eangle[k] - Eangle_value_cut[i] ) / eangle_bin_cut_step );
              iEAng += temp_bin;
            }
          }
          if ( m_debug && ( iDoca == m_idoca_test ) && ( iEAng == m_ieang_test ) )
          {
            doca_test = dd_in[k];
            eang_test = eangle[k];
            t_test->Fill();
            cout << "dedx before cor: " << dEdx[k] * 1.5 / pathlength[k] << endl;
          }
          dEdx[k] = exsvc->StandardHitCorrec(
              0, runFlag, 2, runNO, evtNO, pathlength[k], wid[k], layid[k], dEdx[k],
              ( dd_in[k] ) * ( doca_norm[int( layid[k] )] ), eangle[k], costheta );
          dEdx[k] = exsvc->StandardTrackCorrec( 0, runFlag, 2, runNO, evtNO, dEdx[k], costheta,
                                                tes );
          m_docaeangle[iDoca][iEAng]->Fill( dEdx[k] );
          if ( m_debug && ( iDoca == m_idoca_test ) && ( iEAng == m_ieang_test ) )
          { cout << "dedx after cor: " << dEdx[k] << endl << endl; }
        } // end of if(eangle, doca) are normal
      }   // end of loop in k, hit list
    }     // end of loop in j, entry list
    f->Close();
    f    = 0;
    n103 = 0;
  } // end of loop in i, filelist
  if ( m_debug )
  {
    f_test->cd();
    t_test->Write();
    f_test->Close();
  }
} // end of FillHists()
 
void DedxCalibDocaEAng::AnalyseHists() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalibDocaEAng::AnalyseHists()" << endmsg;
 
  TF1*     func;
  Double_t entry = 0, mean = 0, rms = 0;
  Double_t binmax = 0;
  Double_t lp[3]  = { 0 };
  gStyle->SetOptFit( 1111 );
  for ( Int_t id = 0; id < NumSlicesDoca; id++ )
  {
    for ( Int_t ip = 0; ip < NumSlicesEAng; ip++ )
    {
      entry = m_docaeangle[id][ip]->GetEntries();
      if ( entry < 10 ) continue;
      mean   = m_docaeangle[id][ip]->GetMean();
      rms    = m_docaeangle[id][ip]->GetRMS();
      binmax = m_docaeangle[id][ip]->GetMaximumBin();
      lp[1]  = m_docaeangle[id][ip]->GetBinCenter( binmax );
      lp[2]  = 200;
      lp[0]  = ( m_docaeangle[id][ip]->GetBinContent( binmax ) +
                m_docaeangle[id][ip]->GetBinContent( binmax - 1 ) +
                m_docaeangle[id][ip]->GetBinContent( binmax + 1 ) ) /
              3;
 
      if ( m_phShapeFlag == 0 )
      {
        func = new TF1( "mylan", mylan, 10, 3000, 4 );
        func->SetParameters( entry, mean, rms, -0.8 );
        func->SetParLimits( 0, 0, entry );
        func->SetParLimits( 1, 0, mean );
        func->SetParLimits( 2, 10, rms );
        func->SetParLimits( 3, -1, 0 );
        func->SetParameters( entry / 36., mean - 200, rms / 2.3, -0.5 );
      }
      else if ( m_phShapeFlag == 1 )
      {
        func = new TF1( "landaun", landaun, 10, 3000, 3 );
        func->SetParameters( lp[0], lp[1], lp[2] );
      }
      else if ( m_phShapeFlag == 2 )
      {
        func = new TF1( "Landau", Landau, 10, 3000, 2 );
        func->SetParameters( 0.7 * mean, rms );
      }
      else if ( m_phShapeFlag == 3 )
      {
        func = new TF1( "Vavilov", Vavilov, 10, 3000, 2 );
        func->SetParameters( 0.05, 0.6 );
      }
      else if ( m_phShapeFlag == 4 )
      {
        func = new TF1( "AsymGauss", AsymGauss, 10, 3000, 4 );
        func->SetParameters( entry, mean, rms, 1.0 );
      }
      func->SetLineWidth( 2.1 );
      func->SetLineColor( 2 );
 
      m_docaeangle[id][ip]->Fit( func, "rq" );
    }
  }
}
 
void DedxCalibDocaEAng::WriteHists() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalibDocaEAng::WriteHists()" << endmsg;
 
  const int totalNo           = NumSlicesDoca * NumSlicesEAng;
  double Out_entryNo[totalNo] = { 0 }, Out_mean[totalNo] = { 0 }, Out_sigma[totalNo] = { 0 },
         Out_fitmean[totalNo] = { 0 }, Out_fitmeanerr[totalNo] = { 0 },
         Out_fitsigma[totalNo] = { 0 }, Out_gain[totalNo] = { 0 };
  double Iner_entryNo[totalNo] = { 0 }, Iner_mean[totalNo] = { 0 },
         Iner_sigma[totalNo] = { 0 }, Iner_fitmean[totalNo] = { 0 },
         Iner_fitmeanerr[totalNo] = { 0 }, Iner_fitsigma[totalNo] = { 0 },
         Iner_gain[totalNo] = { 0 };
  double Id_doca[totalNo] = { 0 }, Ip_eangle[totalNo] = { 0 }, doca[totalNo] = { 0 },
         eangle[totalNo] = { 0 }, chisq[totalNo] = { 0 };
 
  double Norm( 0 ), count( 0 );
  for ( Int_t id = 0; id < NumSlicesDoca; id++ )
  {
    for ( Int_t ip = 0; ip < NumSlicesEAng; ip++ )
    {
      Id_doca[id * NumSlicesEAng + ip]   = id;
      Ip_eangle[id * NumSlicesEAng + ip] = ip;
      doca[id * NumSlicesEAng + ip]      = ( id + 0.5 ) * StepDoca + DocaMin;
      int iE = 0, i = 0;
      for ( ; i < 14; i++ )
      {
        if ( ip >= iE && ip < iE + Eangle_cut_bin[i] ) break;
        iE = iE + Eangle_cut_bin[i];
      }
      // cout<<"iE= "<<iE<<" Eangle_cut_bin["<<i<<"]= "<<Eangle_cut_bin[i]<<endl;
      eangle[id * NumSlicesEAng + ip] =
          Eangle_value_cut[i] + ( ip - iE + 0.5 ) *
                                    ( Eangle_value_cut[i + 1] - Eangle_value_cut[i] ) /
                                    Eangle_cut_bin[i];
      // cout<<"eangle["<<id<<"]["<<ip<<"]=  "<<eangle[id*NumSlicesEAng+ip]<<endl;
 
      Out_entryNo[id * NumSlicesEAng + ip] = m_docaeangle[id][ip]->GetEntries();
      Out_mean[id * NumSlicesEAng + ip]    = m_docaeangle[id][ip]->GetMean();
      Out_sigma[id * NumSlicesEAng + ip]   = m_docaeangle[id][ip]->GetRMS();
      // cout<<"Out_entryNo["<<id<<"]["<<ip<<"]= "<<Out_entryNo[id*NumSlicesEAng+ip]<<"
      // "<<Out_mean[id*NumSlicesEAng+ip]<<"  "<<Out_sigma[id*NumSlicesEAng+ip]<<endl;
      if ( Out_entryNo[id * NumSlicesEAng + ip] < 10 ) continue;
 
      if ( m_phShapeFlag == 0 )
      {
        Out_fitmean[id * NumSlicesEAng + ip] =
            m_docaeangle[id][ip]->GetFunction( "mylan" )->GetParameter( 1 );
        Out_fitmeanerr[id * NumSlicesEAng + ip] =
            m_docaeangle[id][ip]->GetFunction( "mylan" )->GetParError( 1 );
        Out_fitsigma[id * NumSlicesEAng + ip] =
            m_docaeangle[id][ip]->GetFunction( "mylan" )->GetParameter( 2 );
        // cout<<Out_fitmean[id*NumSlicesEAng+ip]<<"  "<<Out_fitmeanerr[id*NumSlicesEAng+ip]<<"
        // "<<Out_fitsigma[id*NumSlicesEAng+ip]<<endl;
        chisq[id * NumSlicesEAng + ip] =
            ( m_docaeangle[id][ip]->GetFunction( "mylan" )->GetChisquare() ) /
            ( m_docaeangle[id][ip]->GetFunction( "mylan" )->GetNDF() );
      }
      if ( m_phShapeFlag == 1 )
      {
        Out_fitmean[id * NumSlicesEAng + ip] =
            m_docaeangle[id][ip]->GetFunction( "landaun" )->GetParameter( 1 );
        Out_fitmeanerr[id * NumSlicesEAng + ip] =
            m_docaeangle[id][ip]->GetFunction( "landaun" )->GetParError( 1 );
        Out_fitsigma[id * NumSlicesEAng + ip] =
            m_docaeangle[id][ip]->GetFunction( "landaun" )->GetParameter( 2 );
        chisq[id * NumSlicesEAng + ip] =
            ( m_docaeangle[id][ip]->GetFunction( "landaun" )->GetChisquare() ) /
            ( m_docaeangle[id][ip]->GetFunction( "landaun" )->GetNDF() );
      }
      if ( m_phShapeFlag == 2 )
      {
        Out_fitmean[id * NumSlicesEAng + ip] =
            m_docaeangle[id][ip]->GetFunction( "Landau" )->GetParameter( 1 );
        Out_fitmeanerr[id * NumSlicesEAng + ip] =
            m_docaeangle[id][ip]->GetFunction( "Landau" )->GetParError( 1 );
        Out_fitsigma[id * NumSlicesEAng + ip] =
            m_docaeangle[id][ip]->GetFunction( "Landau" )->GetParameter( 2 );
        chisq[id * NumSlicesEAng + ip] =
            ( m_docaeangle[id][ip]->GetFunction( "Landau" )->GetChisquare() ) /
            ( m_docaeangle[id][ip]->GetFunction( "Landau" )->GetNDF() );
      }
      if ( m_phShapeFlag == 3 )
      {
        Out_fitmean[id * NumSlicesEAng + ip] =
            m_docaeangle[id][ip]->GetFunction( "Vavilov" )->GetParameter( 1 );
        Out_fitmeanerr[id * NumSlicesEAng + ip] =
            m_docaeangle[id][ip]->GetFunction( "Vavilov" )->GetParError( 1 );
        Out_fitsigma[id * NumSlicesEAng + ip] =
            m_docaeangle[id][ip]->GetFunction( "Vavilov" )->GetParameter( 2 );
        chisq[id * NumSlicesEAng + ip] =
            ( m_docaeangle[id][ip]->GetFunction( "Vavilov" )->GetChisquare() ) /
            ( m_docaeangle[id][ip]->GetFunction( "Vavilov" )->GetNDF() );
      }
      if ( m_phShapeFlag == 4 )
      {
        Out_fitmean[id * NumSlicesEAng + ip] =
            m_docaeangle[id][ip]->GetFunction( "AsymGauss" )->GetParameter( 1 );
        Out_fitmeanerr[id * NumSlicesEAng + ip] =
            m_docaeangle[id][ip]->GetFunction( "AsymGauss" )->GetParError( 1 );
        Out_fitsigma[id * NumSlicesEAng + ip] =
            m_docaeangle[id][ip]->GetFunction( "AsymGauss" )->GetParameter( 2 );
        chisq[id * NumSlicesEAng + ip] =
            ( m_docaeangle[id][ip]->GetFunction( "AsymGauss" )->GetChisquare() ) /
            ( m_docaeangle[id][ip]->GetFunction( "AsymGauss" )->GetNDF() );
      }
 
      m_DocaEAngAverdE->SetBinContent( id + 1, ip + 1, Out_fitmean[id * NumSlicesEAng + ip] );
      if ( Out_fitmean[id * NumSlicesEAng + ip] > 0 &&
           Out_fitmeanerr[id * NumSlicesEAng + ip] > 0 &&
           Out_fitmeanerr[id * NumSlicesEAng + ip] < 10 &&
           Out_fitsigma[id * NumSlicesEAng + ip] > 0 &&
           Out_fitsigma[id * NumSlicesEAng + ip] < 200 )
      {
        Norm += ( Out_fitmean[id * NumSlicesEAng + ip] ) *
                ( Out_entryNo[id * NumSlicesEAng + ip] );
        count += Out_entryNo[id * NumSlicesEAng + ip];
      }
    } // end of loop in ip
  }   // end of loop in id
  Norm = Norm / count;
  cout << "count=  " << count << "  average value of doca eangle bins: " << Norm << endl;
  for ( Int_t i = 0; i < totalNo; i++ )
  {
    if ( Out_entryNo[i] > 10 && Out_fitmean[i] > 0 && Out_fitsigma[i] > 0 )
      Out_gain[i] = Out_fitmean[i] / Norm;
    else Out_gain[i] = 0;
  }
 
  log << MSG::INFO << "begin generating root file!!! " << endmsg;
  TFile* f = new TFile( m_rootfile.c_str(), "recreate" );
  for ( Int_t id = 0; id < NumSlicesDoca; id++ )
  {
    for ( Int_t ip = 0; ip < NumSlicesEAng; ip++ ) { m_docaeangle[id][ip]->Write(); }
  }
  m_DocaEAngAverdE->Write();
 
  TTree* ddgcalib = new TTree( "ddgcalib", "ddgcalib" );
  ddgcalib->Branch( "Out_hits", Out_entryNo, "Out_entryNo[1600]/D" );
  ddgcalib->Branch( "Out_mean", Out_mean, "Out_mean[1600]/D" );
  ddgcalib->Branch( "Out_sigma", Out_sigma, "Out_sigma[1600]/D" );
  ddgcalib->Branch( "Out_fitmean", Out_fitmean, "Out_fitmean[1600]/D" );
  ddgcalib->Branch( "Out_fitmeanerr", Out_fitmeanerr, "Out_fitmeanerr[1600]/D" );
  ddgcalib->Branch( "Out_chi", Out_fitsigma, "Out_fitsigma[1600]/D" );
  ddgcalib->Branch( "Out_gain", Out_gain, "Out_gain[1600]/D" );
  ddgcalib->Branch( "Iner_hits", Iner_entryNo, "Iner_entryNo[1600]/D" );
  ddgcalib->Branch( "Iner_mean", Iner_mean, "Iner_mean[1600]/D" );
  ddgcalib->Branch( "Iner_sigma", Iner_sigma, "Iner_sigma[1600]/D" );
  ddgcalib->Branch( "Iner_fitmean", Iner_fitmean, "Iner_fitmean[1600]/D" );
  ddgcalib->Branch( "Iner_fitmeanerr", Iner_fitmeanerr, "Iner_fitmeanerr[1600]/D" );
  ddgcalib->Branch( "Iner_chi", Iner_fitsigma, "Iner_fitsigma[1600]/D" );
  ddgcalib->Branch( "Iner_gain", Iner_gain, "Iner_gain[1600]/D" );
  ddgcalib->Branch( "Id_doca", Id_doca, "Id_doca[1600]/D" );
  ddgcalib->Branch( "Ip_eangle", Ip_eangle, "Ip_eangle[1600]/D" );
  ddgcalib->Branch( "chisq", chisq, "chisq[1600]/D" );
  ddgcalib->Branch( "doca", doca, "doca[1600]/D" );
  ddgcalib->Branch( "eangle", eangle, "eangle[1600]/D" );
  ddgcalib->Fill();
  ddgcalib->Write();
  f->Close();
 
  TCanvas c1( "c1", "canvas", 500, 400 );
  c1.Print( ( m_rootfile + ".ps[" ).c_str() );
  for ( Int_t id = 0; id < NumSlicesDoca; id++ )
  {
    for ( Int_t ip = 0; ip < NumSlicesEAng; ip++ )
    {
      m_docaeangle[id][ip]->Draw();
      c1.Print( ( m_rootfile + ".ps" ).c_str() );
    }
  }
  c1.Print( ( m_rootfile + ".ps]" ).c_str() );
 
  for ( Int_t id = 0; id < NumSlicesDoca; id++ )
  {
    for ( Int_t ip = 0; ip < NumSlicesEAng; ip++ ) { delete m_docaeangle[id][ip]; }
  }
  delete m_DocaEAngAverdE;
}