DedxCalibWireGain.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/DedxCalibWireGain.h"
DECLARE_COMPONENT( DedxCalibWireGain )
using namespace std;
 
DedxCalibWireGain::DedxCalibWireGain( const std::string& name, ISvcLocator* pSvcLocator )
    : DedxCalib( name, pSvcLocator ) {
  declareProperty( "PMin", pMin = 0.2 );
  declareProperty( "PMax", pMax = 2.3 );
}
 
const int wireNo = 6796;
 
void DedxCalibWireGain::BookHists() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalibWireGain::BookHists()" << endmsg;
 
  ReadRecFileList();
 
  m_wiregain = new TH1F*[wireNo];
  stringstream hlabel;
  for ( int i = 0; i < wireNo; i++ )
  {
    hlabel.str( "" );
    hlabel << "dEdx_Wire_" << i;
    m_wiregain[i] = new TH1F( hlabel.str().c_str(), hlabel.str().c_str(), NumHistBins,
                              MinHistValue, MaxHistValue );
  }
}
 
void DedxCalibWireGain::FillHists() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalibWireGain::FillHists()" << endmsg;
 
  string runlist;
 
  double dedx  = 0;
  float  runNO = 0, evtNO = 0, runFlag = 0, pathlength = 0, wid = 0, layid = 0, dd_in = 0,
        driftdist = 0, eangle = 0, zhit = 0, costheta = 0, tes = 0, ptrk = 0;
  cout << "m_recFileVector.size()= " << m_recFileVector.size() << endl;
  for ( int i = 0; i < (int)m_recFileVector.size(); i++ )
  {
    runlist = m_recFileVector[i];
    cout << "runlist: " << runlist.c_str() << endl;
    TFile* f;
    TTree* n102;
    f    = new TFile( runlist.c_str() );
    n102 = (TTree*)f->Get( "n102" );
    n102->SetBranchAddress( "adc_raw", &dedx );
    n102->SetBranchAddress( "path_rphi", &pathlength );
    n102->SetBranchAddress( "wire", &wid );
    n102->SetBranchAddress( "layer", &layid );
    n102->SetBranchAddress( "doca_in", &dd_in );
    n102->SetBranchAddress( "driftdist", &driftdist );
    n102->SetBranchAddress( "eangle", &eangle );
    n102->SetBranchAddress( "zhit", &zhit );
    n102->SetBranchAddress( "runNO", &runNO );
    n102->SetBranchAddress( "evtNO", &evtNO );
    n102->SetBranchAddress( "runFlag", &runFlag );
    n102->SetBranchAddress( "costheta1", &costheta );
    n102->SetBranchAddress( "t01", &tes );
    n102->SetBranchAddress( "ptrk1", &ptrk );
 
    cout << "entries in this file " << n102->GetEntries() << endl;
 
    for ( int j = 0; j < n102->GetEntries(); j++ )
    {
      n102->GetEntry( j );
      if ( tes > 1400 ) continue;
      if ( ptrk > pMax || ptrk < pMin ) continue;
      if ( layid < 8 )
      {
        if ( dedx < MinADCValuecut || dedx > MaxADCValuecut || pathlength > RPhi_PathMaxCut ||
             pathlength < Iner_RPhi_PathMinCut || driftdist > Iner_DriftDistCut )
          continue;
      }
      else
      {
        if ( dedx < MinADCValuecut || dedx > MaxADCValuecut || pathlength > RPhi_PathMaxCut ||
             pathlength < Out_RPhi_PathMinCut || driftdist > Out_DriftDistCut )
          continue;
      }
      dedx = exsvc->StandardHitCorrec( 0, runFlag, 2, runNO, evtNO, pathlength, wid, layid,
                                       dedx, dd_in, eangle, costheta );
      dedx = exsvc->StandardTrackCorrec( 0, runFlag, 2, runNO, evtNO, dedx, costheta, tes );
      m_wiregain[(int)wid]->Fill( dedx );
    }
    f->Delete();
  }
}
 
void DedxCalibWireGain::AnalyseHists() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalibWireGain::AnalyseHists()" << endmsg;
 
  TF1*     func  = new TF1();
  Double_t entry = 0, mean = 0, rms = 0;
  Double_t binmax = 0;
  Double_t lp[3]  = { 0 };
  gStyle->SetOptFit( 1111 );
  for ( int wire = 0; wire < wireNo; wire++ )
  {
    entry = m_wiregain[wire]->Integral();
    if ( entry < 10 ) continue;
    mean   = m_wiregain[wire]->GetMean();
    rms    = m_wiregain[wire]->GetRMS();
    binmax = m_wiregain[wire]->GetMaximumBin();
    lp[1]  = m_wiregain[wire]->GetBinCenter( binmax );
    lp[2]  = 200;
    lp[0]  = ( m_wiregain[wire]->GetBinContent( binmax ) +
              m_wiregain[wire]->GetBinContent( binmax - 1 ) +
              m_wiregain[wire]->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, -3, 3 );
    }
    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_wiregain[wire]->Fit( func, "rq", "", mean - rms, mean + rms / 2 );
  }
}
 
void DedxCalibWireGain::WriteHists() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalibWireGain::WriteHists()" << endmsg;
 
  double entryNo[wireNo] = { 0 }, mean[wireNo] = { 0 }, sigma[wireNo] = { 0 },
         proper[wireNo] = { 0 }, fitmean[wireNo] = { 0 }, fitmeanerr[wireNo] = { 0 },
         fitsigma[wireNo] = { 0 }, wireg[wireNo] = { 0 }, wire[wireNo] = { 0 },
         chisq[wireNo] = { 0 };
  for ( int i = 0; i < wireNo; i++ )
  {
    entryNo[i] = m_wiregain[i]->Integral();
    mean[i]    = m_wiregain[i]->GetMean();
    sigma[i]   = m_wiregain[i]->GetRMS();
    proper[i]  = m_wiregain[i]->GetBinCenter( m_wiregain[i]->GetMaximumBin() );
    if ( entryNo[i] < 10 ) continue;
 
    if ( m_phShapeFlag == 0 )
    {
      fitmean[i]    = m_wiregain[i]->GetFunction( "mylan" )->GetParameter( 1 );
      fitmeanerr[i] = m_wiregain[i]->GetFunction( "mylan" )->GetParError( 1 );
      fitsigma[i]   = m_wiregain[i]->GetFunction( "mylan" )->GetParameter( 2 );
      chisq[i]      = ( m_wiregain[i]->GetFunction( "mylan" )->GetChisquare() ) /
                 ( m_wiregain[i]->GetFunction( "mylan" )->GetNDF() );
    }
    else if ( m_phShapeFlag == 1 )
    {
      fitmean[i]    = m_wiregain[i]->GetFunction( "landaun" )->GetParameter( 1 );
      fitmeanerr[i] = m_wiregain[i]->GetFunction( "landaun" )->GetParError( 1 );
      fitsigma[i]   = m_wiregain[i]->GetFunction( "landaun" )->GetParameter( 2 );
      chisq[i]      = ( m_wiregain[i]->GetFunction( "landaun" )->GetChisquare() ) /
                 ( m_wiregain[i]->GetFunction( "landaun" )->GetNDF() );
    }
    else if ( m_phShapeFlag == 2 )
    {
      fitmean[i]    = m_wiregain[i]->GetFunction( "Landau" )->GetParameter( 1 );
      fitmeanerr[i] = m_wiregain[i]->GetFunction( "Landau" )->GetParError( 1 );
      fitsigma[i]   = m_wiregain[i]->GetFunction( "Landau" )->GetParameter( 2 );
      chisq[i]      = ( m_wiregain[i]->GetFunction( "Landau" )->GetChisquare() ) /
                 ( m_wiregain[i]->GetFunction( "Landau" )->GetNDF() );
    }
    else if ( m_phShapeFlag == 3 )
    {
      fitmean[i]    = m_wiregain[i]->GetFunction( "Vavilov" )->GetParameter( 1 );
      fitmeanerr[i] = m_wiregain[i]->GetFunction( "Vavilov" )->GetParError( 1 );
      fitsigma[i]   = m_wiregain[i]->GetFunction( "Vavilov" )->GetParameter( 2 );
      chisq[i]      = ( m_wiregain[i]->GetFunction( "Vavilov" )->GetChisquare() ) /
                 ( m_wiregain[i]->GetFunction( "Vavilov" )->GetNDF() );
    }
    else if ( m_phShapeFlag == 4 )
    {
      fitmean[i]    = m_wiregain[i]->GetFunction( "AsymGauss" )->GetParameter( 1 );
      fitmeanerr[i] = m_wiregain[i]->GetFunction( "AsymGauss" )->GetParError( 1 );
      fitsigma[i]   = m_wiregain[i]->GetFunction( "AsymGauss" )->GetParameter( 2 );
      chisq[i]      = ( m_wiregain[i]->GetFunction( "AsymGauss" )->GetChisquare() ) /
                 ( m_wiregain[i]->GetFunction( "AsymGauss" )->GetNDF() );
    }
    // cout<<"fitmean[i]= "<<fitmean[i]<<" fitmeanerr[i]= "<<fitmeanerr[i]<<" fitsigma[i]=
    // "<<fitsigma[i]<<endl;
  }
 
  double Norm = 0, count = 0;
  for ( int i = 188; i < wireNo; i++ )
  {
    if ( fitmean[i] > 0 && fitmeanerr[i] > 0 && fitmeanerr[i] < 10 && fitsigma[i] > 0 &&
         fitsigma[i] < 200 && entryNo[i] > 1500 )
    {
      Norm += fitmean[i];
      count++;
    }
  }
  Norm /= count;
  cout << "count=  " << count << "   average value of good wire: " << Norm << endl;
  //  double Norm(550); // use 550 instead of average because we will handle wireg later
  //  outside of the program
 
  for ( int i = 0; i < wireNo; i++ )
  {
    wireg[i] = fitmean[i] / Norm;
    wire[i]  = i;
  }
 
  log << MSG::INFO << "begin generating root file!!! " << endmsg;
  TFile* f = new TFile( m_rootfile.c_str(), "recreate" );
  for ( int i = 0; i < wireNo; i++ ) m_wiregain[i]->Write();
 
  TTree* wiregcalib = new TTree( "wiregcalib", "wiregcalib" );
  wiregcalib->Branch( "wireg", wireg, "wireg[6796]/D" );
  wiregcalib->Branch( "wire", wire, "wire[6796]/D" );
  wiregcalib->Branch( "entryNo", entryNo, "entryNo[6796]/D" );
  wiregcalib->Branch( "proper", proper, "proper[6796]/D" );
  wiregcalib->Branch( "mean", mean, "mean[6796]/D" );
  wiregcalib->Branch( "sigma", sigma, "sigma[6796]/D" );
  wiregcalib->Branch( "fitmean", fitmean, "fitmean[6796]/D" );
  wiregcalib->Branch( "fitmeanerr", fitmeanerr, "fitmeanerr[6796]/D" );
  wiregcalib->Branch( "fitsigma", fitsigma, "fitsigma[6796]/D" );
  wiregcalib->Branch( "chisq", chisq, "chisq[6796]/D" );
  wiregcalib->Fill();
  wiregcalib->Write();
 
  f->Close();
 
  TCanvas c1( "c1", "canvas", 500, 400 );
  c1.Print( ( m_rootfile + ".ps[" ).c_str() );
  for ( int i = 0; i < wireNo; i++ )
  {
    m_wiregain[i]->Draw();
    c1.Print( ( m_rootfile + ".ps" ).c_str() );
  }
  c1.Print( ( m_rootfile + ".ps]" ).c_str() );
 
  for ( int i = 0; i < wireNo; i++ ) delete m_wiregain[i];
}