DedxCalibRunByRun.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 "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()"<<endreq;
 
    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()"<<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};
    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!!!"<<endreq;
        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()"<<endreq;
 
    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()"<<endreq;
 
    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(); 
}