00001 #include "GaudiKernel/MsgStream.h"
00002
00003 #include <sstream>
00004 #include <string>
00005 #include "TTree.h"
00006 #include "TH1F.h"
00007 #include "TH2F.h"
00008 #include "TF1.h"
00009 #include "TFile.h"
00010 #include "TStyle.h"
00011 #include "TCanvas.h"
00012
00013 #include "DedxCalibAlg/DedxCalibDocaEAng.h"
00014
00015 using namespace std;
00016 const double pMin=0.2;
00017 const double pMax=2.1;
00018 const double StepDoca = (DocaMax-DocaMin)/NumSlicesDoca;
00019
00020 DedxCalibDocaEAng::DedxCalibDocaEAng(const std::string& name, ISvcLocator* pSvcLocator): DedxCalib(name, pSvcLocator){
00021 declareProperty("Debug", m_debug=false);
00022 declareProperty("Idoca", m_idoca_test=1);
00023 declareProperty("Ieang", m_ieang_test=1);
00024 }
00025
00026 void DedxCalibDocaEAng::BookHists()
00027 {
00028 MsgStream log(msgSvc(), name());
00029 log<<MSG::INFO<<"DedxCalibDocaEAng::BookHists()"<<endreq;
00030
00031 ReadRecFileList();
00032
00033 m_docaeangle = new TH1F**[NumSlicesDoca];
00034 stringstream hlabel;
00035 for(int i=0;i<NumSlicesDoca;i++)
00036 {
00037 m_docaeangle[i] = new TH1F*[NumSlicesEAng];
00038 for(int j=0;j<NumSlicesEAng;j++)
00039 {
00040 hlabel.str("");
00041 hlabel<<"dEdx_doca"<<i<<"_eang"<<j;
00042 m_docaeangle[i][j] = new TH1F(hlabel.str().c_str(), hlabel.str().c_str(), NumHistBins, MinHistValue, MaxHistValue);
00043 }
00044 }
00045 hlabel.str("");
00046 hlabel<<"dEdxVsDocaEAng";
00047 m_DocaEAngAverdE = new TH2F(hlabel.str().c_str(),"dEdx vs Doca & EAng",NumSlicesDoca,DocaMin,DocaMax,NumSlicesEAng,PhiMin,PhiMax);
00048 }
00049
00050 void DedxCalibDocaEAng::FillHists()
00051 {
00052 MsgStream log(msgSvc(), name());
00053 log<<MSG::INFO<<"DedxCalibDocaEAng::FillHists()"<<endreq;
00054
00055 TFile* f, *f_test(0);
00056 TTree* n103, *t_test(0);
00057 string runlist;
00058
00059 double doca_test(0), eang_test(0);
00060 int ndedxhit=0;
00061 double dEdx[100]={0},pathlength[100]={0},wid[100]={0},layid[100]={0},dd_in[100]={0},eangle[100]={0},zhit[100]={0};
00062 float runNO=0,runFlag=0,costheta=0,tes=0,charge=0,recalg=0,ptrk=0,chidedx=0;
00063 cout<<"m_recFileVector.size()= "<<m_recFileVector.size()<<endl;
00064 if(m_debug){ f_test = new TFile("doca_eangle_test.root","recreate");
00065 t_test = new TTree("inftest", "information for test doca eangle calibration");
00066 t_test->Branch("idoca", &m_idoca_test, "idoca/I");
00067 t_test->Branch("ieang", &m_ieang_test, "ieang/I");
00068 t_test->Branch("doca", &doca_test, "doca/D");
00069 t_test->Branch("eang", &eang_test, "eang/D");
00070 }
00071
00072 for(unsigned int i=0; i<m_recFileVector.size(); i++){
00073 runlist = m_recFileVector[i];
00074 cout<<"runlist: "<<runlist.c_str()<<endl;
00075 f = new TFile(runlist.c_str());
00076
00077 n103 = (TTree*)f->Get("n103");
00078 n103-> SetBranchAddress("ndedxhit", &ndedxhit);
00079 n103-> SetBranchAddress("dEdx_hit",dEdx);
00080 n103-> SetBranchAddress("pathlength_hit",pathlength);
00081 n103-> SetBranchAddress("wid_hit",wid);
00082 n103-> SetBranchAddress("layid_hit",layid);
00083 n103-> SetBranchAddress("dd_in_hit",dd_in);
00084 n103-> SetBranchAddress("eangle_hit",eangle);
00085 n103-> SetBranchAddress("zhit_hit",zhit);
00086 n103-> SetBranchAddress("runNO",&runNO);
00087 n103-> SetBranchAddress("runFlag",&runFlag);
00088 n103-> SetBranchAddress("costheta",&costheta);
00089 n103-> SetBranchAddress("t0",&tes);
00090 n103-> SetBranchAddress("charge",&charge);
00091 n103-> SetBranchAddress("recalg",&recalg);
00092 n103-> SetBranchAddress("ptrk",&ptrk);
00093 n103-> SetBranchAddress("chidedx_e",&chidedx);
00094 for(int j=0;j<n103->GetEntries();j++){
00095 n103->GetEntry(j);
00096 if(ptrk>pMax || ptrk<pMin) continue;
00097 if(tes>1400) continue;
00098 for(int k=0;k<ndedxhit;k++){
00099 dd_in[k] = dd_in[k]/doca_norm[int(layid[k])];
00100 if(eangle[k] >0.25) eangle[k] = eangle[k] -0.5;
00101 else if(eangle[k] <-0.25) eangle[k] = eangle[k] +0.5;
00102 if(eangle[k]>PhiMin && eangle[k]<PhiMax && dd_in[k]>DocaMin && dd_in[k]<DocaMax)
00103 {
00104
00105
00106
00107
00108
00109 int iDoca =(Int_t)floor((dd_in[k]-DocaMin)/StepDoca);
00110 int iEAng = 0;
00111 for(int i =0; i<14; i++)
00112 {
00113 if(eangle[k] < Eangle_value_cut[i] || eangle[k] >= Eangle_value_cut[i+1]) continue;
00114 else if(eangle[k] >= Eangle_value_cut[i] && eangle[k] < Eangle_value_cut[i+1])
00115 {
00116 for(int m =0; m<i; m++)
00117 {
00118 iEAng += Eangle_cut_bin[m];
00119 }
00120 double eangle_bin_cut_step = (Eangle_value_cut[i+1] - Eangle_value_cut[i])/Eangle_cut_bin[i];
00121 int temp_bin = int((eangle[k]-Eangle_value_cut[i])/eangle_bin_cut_step);
00122 iEAng += temp_bin;
00123 }
00124 }
00125 if(m_debug && (iDoca==m_idoca_test) && (iEAng==m_ieang_test)){
00126 doca_test = dd_in[k];
00127 eang_test = eangle[k];
00128 t_test->Fill();
00129 cout << "dedx before cor: " << dEdx[k]*1.5/pathlength[k] << endl;
00130 }
00131 dEdx[k] = exsvc->StandardHitCorrec(0,runFlag,2,runNO,pathlength[k],wid[k],layid[k],dEdx[k],(dd_in[k])*(doca_norm[int(layid[k])]),eangle[k],costheta);
00132 dEdx[k] = exsvc->StandardTrackCorrec(0, runFlag, 2, runNO, dEdx[k], costheta, tes);
00133 m_docaeangle[iDoca][iEAng]->Fill(dEdx[k]);
00134 if(m_debug && (iDoca==m_idoca_test) && (iEAng==m_ieang_test)){
00135 cout << "dedx after cor: " << dEdx[k] << endl << endl;
00136 }
00137 }
00138 }
00139 }
00140 }
00141 if(m_debug){ f_test->cd(); t_test->Write(); f_test->Close(); }
00142 }
00143
00144 void DedxCalibDocaEAng::AnalyseHists()
00145 {
00146 MsgStream log(msgSvc(), name());
00147 log<<MSG::INFO<<"DedxCalibDocaEAng::AnalyseHists()"<<endreq;
00148
00149 TF1* func;
00150 Double_t entry=0,mean=0,rms=0;
00151 Double_t binmax=0;
00152 Double_t lp[3]={0};
00153 gStyle->SetOptFit(1111);
00154 for(Int_t id=0;id<NumSlicesDoca;id++)
00155 {
00156 for(Int_t ip=0;ip<NumSlicesEAng;ip++)
00157 {
00158 entry = m_docaeangle[id][ip]->GetEntries();
00159 if(entry<10) continue;
00160 mean = m_docaeangle[id][ip]->GetMean();
00161 rms = m_docaeangle[id][ip]->GetRMS();
00162 binmax = m_docaeangle[id][ip]->GetMaximumBin();
00163 lp[1] = m_docaeangle[id][ip]->GetBinCenter(binmax);
00164 lp[2] = 200;
00165 lp[0] = (m_docaeangle[id][ip]->GetBinContent(binmax)+m_docaeangle[id][ip]->GetBinContent(binmax-1)+m_docaeangle[id][ip]->GetBinContent(binmax+1))/3;
00166
00167 if(m_phShapeFlag==0)
00168 {
00169 func = new TF1("mylan",mylan,10,3000,4);
00170 func->SetParameters(entry, mean, rms, -0.8);
00171 func->SetParLimits(0, 0, entry);
00172 func->SetParLimits(1, 0, mean);
00173 func->SetParLimits(2, 10, rms);
00174 func->SetParLimits(3, -1, 0);
00175 func->SetParameters(entry/36., mean-200, rms/2.3, -0.5);
00176 }
00177 else if(m_phShapeFlag==1)
00178 {
00179 func = new TF1("landaun",landaun,10,3000,3);
00180 func->SetParameters(lp[0], lp[1], lp[2] );
00181 }
00182 else if(m_phShapeFlag==2)
00183 {
00184 func = new TF1("Landau",Landau,10,3000,2);
00185 func->SetParameters(0.7*mean, rms );
00186 }
00187 else if(m_phShapeFlag==3)
00188 {
00189 func = new TF1("Vavilov",Vavilov,10,3000,2);
00190 func->SetParameters(0.05, 0.6);
00191 }
00192 else if(m_phShapeFlag==4)
00193 {
00194 func = new TF1("AsymGauss",AsymGauss,10,3000,4);
00195 func->SetParameters(entry, mean, rms, 1.0 );
00196 }
00197 func->SetLineWidth(2.1);
00198 func->SetLineColor(2);
00199
00200 m_docaeangle[id][ip]->Fit(func, "r");
00201 }
00202 }
00203 }
00204
00205 void DedxCalibDocaEAng::WriteHists()
00206 {
00207 MsgStream log(msgSvc(), name());
00208 log<<MSG::INFO<<"DedxCalibDocaEAng::WriteHists()"<<endreq;
00209
00210 const int totalNo = NumSlicesDoca*NumSlicesEAng;
00211 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};
00212 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};
00213 double Id_doca[totalNo]={0},Ip_eangle[totalNo]={0},doca[totalNo]={0},eangle[totalNo]={0},chisq[totalNo]={0};
00214
00215 double Norm(0), count(0);
00216 for(Int_t id=0;id<NumSlicesDoca;id++)
00217 {
00218 for(Int_t ip=0;ip<NumSlicesEAng;ip++)
00219 {
00220 Id_doca[id*NumSlicesEAng+ip] = id;
00221 Ip_eangle[id*NumSlicesEAng+ip] = ip;
00222 doca[id*NumSlicesEAng+ip] = (id+0.5)*StepDoca + DocaMin;
00223 int iE=0,i=0;
00224 for(;i<14;i++)
00225 {
00226 if(ip>=iE && ip<iE+Eangle_cut_bin[i]) break;
00227 iE=iE+Eangle_cut_bin[i];
00228 }
00229
00230 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];
00231
00232
00233 Out_entryNo[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetEntries();
00234 Out_mean[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetMean();
00235 Out_sigma[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetRMS();
00236
00237 if(Out_entryNo[id*NumSlicesEAng+ip]<10) continue;
00238
00239 if(m_phShapeFlag==0)
00240 {
00241 Out_fitmean[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("mylan")->GetParameter(1);
00242 Out_fitmeanerr[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("mylan")->GetParError(1);
00243 Out_fitsigma[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("mylan")->GetParameter(2);
00244
00245 chisq[id*NumSlicesEAng+ip] = (m_docaeangle[id][ip]->GetFunction("mylan")-> GetChisquare())/(m_docaeangle[id][ip]->GetFunction("mylan")-> GetNDF());
00246 }
00247 if(m_phShapeFlag==1)
00248 {
00249 Out_fitmean[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("landaun")->GetParameter(1);
00250 Out_fitmeanerr[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("landaun")->GetParError(1);
00251 Out_fitsigma[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("landaun")->GetParameter(2);
00252 chisq[id*NumSlicesEAng+ip] = (m_docaeangle[id][ip]->GetFunction("landaun")-> GetChisquare())/(m_docaeangle[id][ip]->GetFunction("landaun")-> GetNDF());
00253 }
00254 if(m_phShapeFlag==2)
00255 {
00256 Out_fitmean[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("Landau")->GetParameter(1);
00257 Out_fitmeanerr[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("Landau")->GetParError(1);
00258 Out_fitsigma[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("Landau")->GetParameter(2);
00259 chisq[id*NumSlicesEAng+ip] = (m_docaeangle[id][ip]->GetFunction("Landau")-> GetChisquare())/(m_docaeangle[id][ip]->GetFunction("Landau")-> GetNDF());
00260 }
00261 if(m_phShapeFlag==3)
00262 {
00263 Out_fitmean[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("Vavilov")->GetParameter(1);
00264 Out_fitmeanerr[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("Vavilov")->GetParError(1);
00265 Out_fitsigma[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("Vavilov")->GetParameter(2);
00266 chisq[id*NumSlicesEAng+ip] = (m_docaeangle[id][ip]->GetFunction("Vavilov")-> GetChisquare())/(m_docaeangle[id][ip]->GetFunction("Vavilov")-> GetNDF());
00267 }
00268 if(m_phShapeFlag==4)
00269 {
00270 Out_fitmean[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("AsymGauss")->GetParameter(1);
00271 Out_fitmeanerr[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("AsymGauss")->GetParError(1);
00272 Out_fitsigma[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("AsymGauss")->GetParameter(2);
00273 chisq[id*NumSlicesEAng+ip] = (m_docaeangle[id][ip]->GetFunction("AsymGauss")-> GetChisquare())/(m_docaeangle[id][ip]->GetFunction("AsymGauss")-> GetNDF());
00274 }
00275
00276 m_DocaEAngAverdE -> SetBinContent(id+1, ip+1, Out_fitmean[id*NumSlicesEAng+ip]);
00277 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){
00278 Norm += (Out_fitmean[id*NumSlicesEAng+ip])*(Out_entryNo[id*NumSlicesEAng+ip]);
00279 count += Out_entryNo[id*NumSlicesEAng+ip];
00280 }
00281 }
00282 }
00283 Norm = Norm/count;
00284 cout<<"count= "<<count<<" average value of doca eangle bins: "<<Norm<<endl;
00285 for(Int_t i=0;i<totalNo;i++)
00286 {
00287 if(Out_entryNo[i]>10 && Out_fitmean[i]>0 && Out_fitsigma[i]>0) Out_gain[i]=Out_fitmean[i]/Norm;
00288 else Out_gain[i] = 0;
00289 }
00290
00291 log<<MSG::INFO<<"begin generating root file!!! "<<endreq;
00292 TFile* f = new TFile(m_rootfile.c_str(), "recreate");
00293 for(Int_t id=0;id<NumSlicesDoca;id++)
00294 {
00295 for(Int_t ip=0;ip<NumSlicesEAng;ip++)
00296 {
00297 m_docaeangle[id][ip]->Write();
00298 }
00299 }
00300 m_DocaEAngAverdE->Write();
00301
00302 TTree* ddgcalib = new TTree("ddgcalib", "ddgcalib");
00303 ddgcalib -> Branch("Out_hits", Out_entryNo, "Out_entryNo[1600]/D");
00304 ddgcalib -> Branch("Out_mean", Out_mean, "Out_mean[1600]/D");
00305 ddgcalib -> Branch("Out_sigma", Out_sigma, "Out_sigma[1600]/D");
00306 ddgcalib -> Branch("Out_fitmean", Out_fitmean, "Out_fitmean[1600]/D");
00307 ddgcalib -> Branch("Out_fitmeanerr", Out_fitmeanerr, "Out_fitmeanerr[1600]/D");
00308 ddgcalib -> Branch("Out_chi", Out_fitsigma, "Out_fitsigma[1600]/D");
00309 ddgcalib -> Branch("Out_gain", Out_gain, "Out_gain[1600]/D");
00310 ddgcalib -> Branch("Iner_hits", Iner_entryNo, "Iner_entryNo[1600]/D");
00311 ddgcalib -> Branch("Iner_mean", Iner_mean, "Iner_mean[1600]/D");
00312 ddgcalib -> Branch("Iner_sigma", Iner_sigma, "Iner_sigma[1600]/D");
00313 ddgcalib -> Branch("Iner_fitmean", Iner_fitmean, "Iner_fitmean[1600]/D");
00314 ddgcalib -> Branch("Iner_fitmeanerr", Iner_fitmeanerr, "Iner_fitmeanerr[1600]/D");
00315 ddgcalib -> Branch("Iner_chi", Iner_fitsigma, "Iner_fitsigma[1600]/D");
00316 ddgcalib -> Branch("Iner_gain", Iner_gain, "Iner_gain[1600]/D");
00317 ddgcalib -> Branch("Id_doca", Id_doca, "Id_doca[1600]/D");
00318 ddgcalib -> Branch("Ip_eangle", Ip_eangle, "Ip_eangle[1600]/D");
00319 ddgcalib -> Branch("chisq", chisq, "chisq[1600]/D");
00320 ddgcalib -> Branch("doca", doca, "doca[1600]/D");
00321 ddgcalib -> Branch("eangle", eangle, "eangle[1600]/D");
00322 ddgcalib -> Fill();
00323 ddgcalib -> Write();
00324 f->Close();
00325
00326 TCanvas c1("c1", "canvas", 500, 400);
00327 c1.Print((m_rootfile+".ps[").c_str());
00328 for(Int_t id=0;id<NumSlicesDoca;id++)
00329 {
00330 for(Int_t ip=0;ip<NumSlicesEAng;ip++)
00331 {
00332 m_docaeangle[id][ip]->Draw();
00333 c1.Print((m_rootfile+".ps").c_str());
00334 }
00335 }
00336 c1.Print((m_rootfile+".ps]").c_str());
00337
00338 for(Int_t id=0;id<NumSlicesDoca;id++)
00339 {
00340 for(Int_t ip=0;ip<NumSlicesEAng;ip++)
00341 {
00342 delete m_docaeangle[id][ip];
00343 }
00344 }
00345 delete m_DocaEAngAverdE;
00346 }