00001 #include "GaudiKernel/MsgStream.h" 00002 00003 #include <sstream> 00004 #include <string> 00005 #include "TFile.h" 00006 #include "TF1.h" 00007 #include "TH1F.h" 00008 #include "TTree.h" 00009 00010 #include "DedxCalibAlg/DedxCalibRunByRun.h" 00011 00012 using namespace std; 00013 00014 int runNo = 0; 00015 00016 DedxCalibRunByRun::DedxCalibRunByRun(const std::string& name, ISvcLocator* pSvcLocator): DedxCalib(name, pSvcLocator){} 00017 00018 00019 void DedxCalibRunByRun::BookHists() 00020 { 00021 MsgStream log(msgSvc(), name()); 00022 log<<MSG::INFO<<"DedxCalibRunByRun::BookHists()"<<endreq; 00023 00024 ReadRecFileList(); 00025 runNo = m_recFileVector.size(); 00026 00027 m_hist = new TH1F("dEdx_allRun","dEdx_allRun",NumHistBins,MinHistValue1,MaxHistValue1); 00028 m_hist->GetXaxis()->SetTitle("dE/dx"); 00029 m_hist->GetYaxis()->SetTitle("entries"); 00030 m_rungain = new TH1F*[runNo]; 00031 stringstream hlabel; 00032 for(int i=0; i<runNo; i++) 00033 { 00034 hlabel.str(""); 00035 hlabel<<"dEdx_Run_"<<i; 00036 m_rungain[i] = new TH1F(hlabel.str().c_str(), hlabel.str().c_str(),NumHistBins,MinHistValue1,MaxHistValue1); 00037 } 00038 } 00039 00040 void DedxCalibRunByRun::FillHists() 00041 { 00042 MsgStream log(msgSvc(), name()); 00043 log<<MSG::INFO<<"DedxCalibRunByRun::FillHists()"<<endreq; 00044 00045 TFile* f; 00046 TTree* n103; 00047 string runlist; 00048 00049 int ndedxhit=0; 00050 double dEdx[100]={0},pathlength[100]={0},wid[100]={0},layid[100]={0},dd_in[100]={0},eangle[100]={0},zhit[100]={0}; 00051 double dedx=0; 00052 float runNO=0,runFlag=0,costheta=0,tes=0; 00053 vector<double> phlist; 00054 for(int i=0; i<runNo; i++) 00055 { 00056 runlist = m_recFileVector[i]; 00057 f = new TFile(runlist.c_str()); 00058 n103 = (TTree*)f->Get("n103"); 00059 n103-> SetBranchAddress("ndedxhit", &ndedxhit); 00060 n103-> SetBranchAddress("dEdx_hit",dEdx); 00061 n103-> SetBranchAddress("pathlength_hit",pathlength); 00062 n103-> SetBranchAddress("wid_hit",wid); 00063 n103-> SetBranchAddress("layid_hit",layid); 00064 n103-> SetBranchAddress("dd_in_hit",dd_in); 00065 n103-> SetBranchAddress("eangle_hit",eangle); 00066 n103-> SetBranchAddress("zhit_hit",zhit); 00067 n103-> SetBranchAddress("runNO",&runNO); 00068 n103-> SetBranchAddress("runFlag",&runFlag); 00069 n103-> SetBranchAddress("costheta",&costheta); 00070 n103-> SetBranchAddress("t0",&tes); 00071 for(int j=0;j<n103->GetEntries();j++) 00072 { 00073 phlist.clear(); 00074 n103->GetEntry(j); 00075 for(int k=0;k<ndedxhit;k++) 00076 { 00077 dEdx[k] = exsvc->StandardHitCorrec(0,runFlag,2,runNO,pathlength[k],wid[k],layid[k],dEdx[k],dd_in[k],eangle[k],costheta); 00078 phlist.push_back(dEdx[k]); 00079 } 00080 dedx = cal_dedx_bitrunc(truncate, phlist); 00081 dedx = exsvc->StandardTrackCorrec(0, runFlag, 2, runNO, dedx, costheta, tes); 00082 m_rungain[i]->Fill(dedx); 00083 m_hist->Fill(dedx); 00084 } 00085 cout<<"runNO = "<<runNO<<endl; 00086 m_runNOVector.push_back(runNO); 00087 } 00088 } 00089 00090 void DedxCalibRunByRun::AnalyseHists() 00091 { 00092 MsgStream log(msgSvc(), name()); 00093 log<<MSG::INFO<<"DedxCalibRunByRun::AnalyseHists()"<<endreq; 00094 00095 m_hist->Fit("gaus", "Q" ); 00096 for(int i=0; i<runNo; i++) 00097 { 00098 m_rungain[i]->Fit("gaus", "Q" ); 00099 } 00100 } 00101 00102 void DedxCalibRunByRun::WriteHists() 00103 { 00104 MsgStream log(msgSvc(), name()); 00105 log<<MSG::INFO<<"DedxCalibRunByRun::WriteHists()"<<endreq; 00106 00107 TFile* f = new TFile(m_rootfile.c_str(), "recreate"); 00108 00109 Double_t gainpar=0, resolpar=0; 00110 gainpar = m_hist -> GetFunction("gaus") -> GetParameter(1); 00111 resolpar = m_hist -> GetFunction("gaus") -> GetParameter(2); 00112 00113 TTree* gain = new TTree("gaincalib", "gaincalib"); 00114 gain -> Branch("gain", &gainpar, "gain[1]/D"); 00115 gain->Fill(); 00116 TTree* resol = new TTree("resolcalib", "resolcalib"); 00117 resol -> Branch("resol", &resolpar, "resol[1]/D"); 00118 resol->Fill(); 00119 00120 Int_t runno=0; 00121 Double_t runmean=0, rungain=0, runresol=0; 00122 TTree* runbyrun = new TTree("runcalib", "runcalib"); 00123 runbyrun -> Branch("runno", &runno, "runno/I"); 00124 runbyrun -> Branch("runmean", &runmean, "runmean/D"); 00125 runbyrun -> Branch("rungain", &rungain, "rungain/D"); 00126 runbyrun -> Branch("runresol", &runresol, "runresol/D"); 00127 for(int i=0; i<runNo; i++) 00128 { 00129 runno = m_runNOVector[i]; 00130 runmean = m_rungain[i] -> GetFunction("gaus") -> GetParameter(1); 00131 runresol = m_rungain[i] -> GetFunction("gaus") -> GetParameter(2); 00132 rungain = runmean/NormalMean; 00133 cout<<"i = "<<i<<" runno = "<<runno <<" @ runmean = "<<runmean<<" @ rungain = "<<rungain<<" @ runresol = "<<runresol<<endl; 00134 runbyrun -> Fill(); 00135 } 00136 00137 m_hist->Write(); 00138 for(int i=0; i<runNo; i++) m_rungain[i]->Write(); 00139 gain->Write(); 00140 resol -> Write(); 00141 runbyrun -> Write(); 00142 00143 delete m_hist; 00144 for(int i=0; i<runNo; i++) delete m_rungain[i]; 00145 delete gain; 00146 delete resol; 00147 delete runbyrun; 00148 f->Close(); 00149 }