DedxCalibTzero.cxx

Go to the documentation of this file.

#include "GaudiKernel/MsgStream.h"
 
#include <sstream>
#include <string>
#include "TFile.h"
#include "TF1.h"
#include "TH1F.h"
#include "TTree.h"
#include "TStyle.h"
#include "TCanvas.h"
 
#include "DedxCalibAlg/DedxCalibTzero.h"
 
#define Size  700000
#define numcut 50
 
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()"<<endreq;
 
    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()"<<endreq;
 
    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()"<<endreq;
 
    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()"<<endreq;
 
    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!!! "<<endreq;
    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;
}