00001 #include "GaudiKernel/MsgStream.h"
00002
00003 #include <sstream>
00004 #include <string>
00005 #include "TTree.h"
00006 #include "TH1F.h"
00007 #include "TF1.h"
00008 #include "TFile.h"
00009 #include "TStyle.h"
00010 #include "TCanvas.h"
00011
00012 #include "DedxCalibAlg/DedxCalibLayerGain.h"
00013
00014 using namespace std;
00015
00016 DedxCalibLayerGain::DedxCalibLayerGain(const std::string& name, ISvcLocator* pSvcLocator): DedxCalib(name, pSvcLocator){}
00017
00018 const int layNo = 43;
00019 double m_mean,m_rms;
00020 void calculate(vector<double> phlist)
00021 {
00022 double mean=0,rms=0;
00023 for(int i=0;i<phlist.size();i++)
00024 {
00025 mean+=phlist[i];
00026 }
00027 mean/=phlist.size();
00028 for(int i=0;i<phlist.size();i++)
00029 {
00030 rms += pow((phlist[i]-mean),2);
00031 }
00032
00033 rms = sqrt(rms/phlist.size());
00034
00035 m_mean = mean;
00036 m_rms = rms;
00037 }
00038 double getMean() {return m_mean;}
00039 double getRms() {return m_rms;}
00040
00041 void DedxCalibLayerGain::BookHists()
00042 {
00043 MsgStream log(msgSvc(), name());
00044 log<<MSG::INFO<<"DedxCalibLayerGain::BookHists()"<<endreq;
00045
00046 ReadRecFileList();
00047
00048 m_laygain = new TH1F*[layNo];
00049 m_laygain_gaus = new TH1F*[layNo];
00050 stringstream hlabel;
00051 for(int i=0; i<layNo; i++)
00052 {
00053 hlabel.str("");
00054 hlabel<<"dEdx_Lay_"<<i;
00055 m_laygain[i] = new TH1F(hlabel.str().c_str(), hlabel.str().c_str(), NumHistBins, MinHistValue, MaxHistValue);
00056 hlabel.str("");
00057 hlabel<<"dEdx_gaus_Lay_"<<i;
00058 m_laygain_gaus[i] = new TH1F(hlabel.str().c_str(), hlabel.str().c_str(), NumHistBins, MinHistValue1, MaxHistValue1);
00059 }
00060 }
00061
00062 void DedxCalibLayerGain::FillHists()
00063 {
00064 MsgStream log(msgSvc(), name());
00065 log<<MSG::INFO<<"DedxCalibLayerGain::FillHists()"<<endreq;
00066
00067 TFile* f;
00068 TTree* n102;
00069 string runlist;
00070
00071 double dedx=0;
00072 float runNO=0,runFlag=0,pathlength=0,wid=0,layid=0,dd_in=0,driftdist=0,eangle=0,zhit=0,costheta=0,tes=0,ptrk=0;
00073 cout<<"m_recFileVector.size()= "<<m_recFileVector.size()<<endl;
00074 for(int i=0; i<m_recFileVector.size(); i++)
00075 {
00076 runlist = m_recFileVector[i];
00077 cout<<"runlist: "<<runlist.c_str()<<endl;
00078 f = new TFile(runlist.c_str());
00079 n102 = (TTree*)f->Get("n102");
00080 n102-> SetBranchAddress("adc_raw",&dedx);
00081 n102-> SetBranchAddress("path_rphi",&pathlength);
00082 n102-> SetBranchAddress("wire",&wid);
00083 n102-> SetBranchAddress("layer",&layid);
00084 n102-> SetBranchAddress("doca_in",&dd_in);
00085 n102-> SetBranchAddress("driftdist",&driftdist);
00086 n102-> SetBranchAddress("eangle",&eangle);
00087 n102-> SetBranchAddress("zhit",&zhit);
00088 n102-> SetBranchAddress("runNO",&runNO);
00089 n102-> SetBranchAddress("runFlag",&runFlag);
00090 n102-> SetBranchAddress("costheta1",&costheta);
00091 n102-> SetBranchAddress("t01",&tes);
00092 n102->SetBranchAddress("ptrk1",&ptrk);
00093
00094 for(int j=0;j<n102->GetEntries();j++)
00095 {
00096 n102->GetEntry(j);
00097 if(tes>1400) continue;
00098 if (ptrk<0.1) continue;
00099
00100 if(layid <8)
00101 {if(dedx<MinADCValuecut || dedx>MaxADCValuecut || pathlength>RPhi_PathMaxCut || pathlength<Iner_RPhi_PathMinCut || driftdist>Iner_DriftDistCut) continue;}
00102 else
00103 {if(dedx<MinADCValuecut || dedx>MaxADCValuecut || pathlength>RPhi_PathMaxCut || pathlength<Out_RPhi_PathMinCut || driftdist>Out_DriftDistCut) continue;}
00104 dedx = exsvc->StandardHitCorrec(0,runFlag,2,runNO,pathlength,wid,layid,dedx,dd_in,eangle,costheta);
00105 dedx = exsvc->StandardTrackCorrec(0, runFlag, 2, runNO, dedx, costheta, tes);
00106 m_laygain[(int)layid]->Fill(dedx);
00107
00108 m_vector[(int)layid].push_back(dedx);
00109 }
00110 }
00111
00112 }
00113
00114 void DedxCalibLayerGain::AnalyseHists()
00115 {
00116 MsgStream log(msgSvc(), name());
00117 log<<MSG::INFO<<"DedxCalibLayerGain::AnalyseHists()"<<endreq;
00118
00119 TF1* func;
00120 Double_t entry=0,mean=0,rms=0;
00121 Double_t binmax=0;
00122 Double_t lp[3]={0};
00123 gStyle->SetOptFit(1111);
00124
00125 stringstream hlabel;
00126 double dedxt;
00127 vector<double> phlist;
00128 vector<double> phlist_gaus;
00129 for(int i=0; i<layNo; i++)
00130 {
00131 entry = m_laygain[i]->GetEntries();
00132 if(entry<10) continue;
00133 mean = m_laygain[i]->GetMean();
00134 rms = m_laygain[i]->GetRMS();
00135 binmax = m_laygain[i]->GetMaximumBin();
00136 lp[1] = m_laygain[i]->GetBinCenter(binmax);
00137 lp[2] = 200;
00138 lp[0] = (m_laygain[i]->GetBinContent(binmax)+m_laygain[i]->GetBinContent(binmax-1)+m_laygain[i]->GetBinContent(binmax+1))/3;
00139
00140 if(m_phShapeFlag==0)
00141 {
00142 func = new TF1("mylan",mylan,10,3000,4);
00143 func->SetParameters(entry, mean, rms, -0.9);
00144 }
00145 else if(m_phShapeFlag==1)
00146 {
00147 func = new TF1("landaun",landaun,10,3000,3);
00148 func->SetParameters(lp[0], lp[1], lp[2] );
00149 }
00150 else if(m_phShapeFlag==2)
00151 {
00152 func = new TF1("Landau",Landau,10,3000,2);
00153 func->SetParameters(0.7*mean, rms );
00154 }
00155 else if(m_phShapeFlag==3)
00156 {
00157 func = new TF1("Vavilov",Vavilov,10,3000,2);
00158 func->SetParameters(0.05, 0.6);
00159 }
00160 else if(m_phShapeFlag==4)
00161 {
00162 func = new TF1("AsymGauss",AsymGauss,10,3000,4);
00163 func->SetParameters(entry, mean, rms, 1.0 );
00164 }
00165 func->SetLineWidth(2.1);
00166 func->SetLineColor(2);
00167
00168 m_laygain[i]->Fit(func, "r" );
00169
00170
00171
00172 for(int j=0;j<m_vector[i].size();j+=100)
00173 {
00174 for(int k=0;k<100;k++) phlist.push_back(m_vector[i][j+k]);
00175 dedxt = cal_dedx_bitrunc(truncate, phlist);
00176 phlist_gaus.push_back(dedxt);
00177 phlist.clear();
00178 }
00179
00180 hlabel.str("");
00181 hlabel<<"dEdx_gaus_Lay_"<<i;
00182 calculate(phlist_gaus);
00183
00184 m_laygain_gaus[i] = new TH1F(hlabel.str().c_str(), hlabel.str().c_str(), NumHistBins, getMean()-10*getRms(), getMean()+10*getRms());
00185 for(int j=0;j<phlist_gaus.size();j++) m_laygain_gaus[i]->Fill(phlist_gaus[j]);
00186 phlist_gaus.clear();
00187 m_laygain_gaus[i]->Fit("gaus","Q");
00188
00189
00190 }
00191 }
00192
00193 void DedxCalibLayerGain::WriteHists()
00194 {
00195 MsgStream log(msgSvc(), name());
00196 log<<MSG::INFO<<"DedxCalibLayerGain::WriteHists()"<<endreq;
00197
00198
00199 double entryNo[layNo]={0},mean[layNo]={0},sigma[layNo]={0},proper[layNo]={0},fitmean[layNo]={0},fitmeanerr[layNo]={0},fitsigma[layNo]={0},layerg[layNo]={0},layer[layNo]={0},chisq[layNo]={0};
00200 double fitmean_gaus[layNo]={0},fitmeanerr_gaus[layNo]={0},fitsigma_gaus[layNo]={0},layerg_gaus[layNo]={0};
00201 for(int i=0; i<layNo; i++)
00202 {
00203 fitmean_gaus[i] = m_laygain_gaus[i]->GetFunction("gaus")->GetParameter(1);
00204 fitmeanerr_gaus[i] = m_laygain_gaus[i]->GetFunction("gaus")->GetParError(1);
00205 fitsigma_gaus[i] = m_laygain_gaus[i]->GetFunction("gaus")->GetParameter(2);
00206
00207 entryNo[i] = m_laygain[i]->GetEntries();
00208 mean[i] = m_laygain[i]->GetMean();
00209 sigma[i] = m_laygain[i]->GetRMS();
00210 proper[i] = m_laygain[i]->GetBinCenter(m_laygain[i]->GetMaximumBin());
00211 if(entryNo[i]<10) continue;
00212
00213 if(m_phShapeFlag==0)
00214 {
00215 fitmean[i] = m_laygain[i]->GetFunction("mylan")->GetParameter(1);
00216 fitmeanerr[i] = m_laygain[i]->GetFunction("mylan")->GetParError(1);
00217 fitsigma[i] = m_laygain[i]->GetFunction("mylan")->GetParameter(2);
00218 chisq[i] = (m_laygain[i]->GetFunction("mylan")-> GetChisquare())/(m_laygain[i]->GetFunction("mylan")-> GetNDF());
00219 }
00220 else if(m_phShapeFlag==1)
00221 {
00222 fitmean[i] = m_laygain[i]->GetFunction("landaun")->GetParameter(1);
00223 fitmeanerr[i] = m_laygain[i]->GetFunction("landaun")->GetParError(1);
00224 fitsigma[i] = m_laygain[i]->GetFunction("landaun")->GetParameter(2);
00225 chisq[i] = (m_laygain[i]->GetFunction("landaun")-> GetChisquare())/(m_laygain[i]->GetFunction("landaun")-> GetNDF());
00226 }
00227 else if(m_phShapeFlag==2)
00228 {
00229 fitmean[i] = m_laygain[i]->GetFunction("Landau")->GetParameter(1);
00230 fitmeanerr[i] = m_laygain[i]->GetFunction("Landau")->GetParError(1);
00231 fitsigma[i] = m_laygain[i]->GetFunction("Landau")->GetParameter(2);
00232 chisq[i] = (m_laygain[i]->GetFunction("Landau")-> GetChisquare())/(m_laygain[i]->GetFunction("Landau")-> GetNDF());
00233 }
00234 else if(m_phShapeFlag==3)
00235 {
00236 fitmean[i] = m_laygain[i]->GetFunction("Vavilov")->GetParameter(1);
00237 fitmeanerr[i] = m_laygain[i]->GetFunction("Vavilov")->GetParError(1);
00238 fitsigma[i] = m_laygain[i]->GetFunction("Vavilov")->GetParameter(2);
00239 chisq[i] = (m_laygain[i]->GetFunction("Vavilov")-> GetChisquare())/(m_laygain[i]->GetFunction("Vavilov")-> GetNDF());
00240 }
00241 else if(m_phShapeFlag==4)
00242 {
00243 fitmean[i] = m_laygain[i]->GetFunction("AsymGauss")->GetParameter(1);
00244 fitmeanerr[i] = m_laygain[i]->GetFunction("AsymGauss")->GetParError(1);
00245 fitsigma[i] = m_laygain[i]->GetFunction("AsymGauss")->GetParameter(2);
00246 chisq[i] = (m_laygain[i]->GetFunction("AsymGauss")-> GetChisquare())/(m_laygain[i]->GetFunction("AsymGauss")-> GetNDF());
00247 }
00248
00249 }
00250
00251 double sum=0, sum_gaus=0, n=0;
00252 for(int i=4; i<layNo; i++)
00253 {
00254 if(fitmean[i]>0 && fitmeanerr[i]>0 && fitmeanerr[i]<10 && fitsigma[i]>0 && fitsigma[i]<200 && entryNo[i]>1500)
00255 {
00256 sum += fitmean[i];
00257 sum_gaus+= fitmean_gaus[i];
00258 n++;
00259 }
00260 }
00261 sum/=n;
00262 sum_gaus/=n;
00263 cout<<"average value of good layer: "<<sum<<endl;
00264
00265 for(int i=0;i<layNo;i++)
00266 {
00267 layerg[i] = fitmean[i]/sum;
00268 layerg_gaus[i] = fitmean_gaus[i]/sum_gaus;
00269 layer[i] = i;
00270 }
00271
00272
00273 log<<MSG::INFO<<"begin generating root file!!! "<<endreq;
00274 TFile* f = new TFile(m_rootfile.c_str(), "recreate");
00275 for(int i=0;i<layNo;i++)
00276 {
00277 m_laygain[i]->Write();
00278 m_laygain_gaus[i]->Write();
00279 }
00280
00281 TTree* layergcalib = new TTree("layergcalib", "layergcalib");
00282 layergcalib->Branch("layerg_gaus", layerg_gaus, "layerg_gaus[43]/D");
00283 layergcalib->Branch("layerg", layerg, "layerg[43]/D");
00284 layergcalib->Branch("layer", layer, "layer[43]/D");
00285 layergcalib->Branch("entryNo", entryNo, "entryNo[43]/D");
00286 layergcalib->Branch("mean", mean, "mean[43]/D");
00287 layergcalib->Branch("sigma", sigma, "sigma[43]/D");
00288 layergcalib->Branch("fitmean", fitmean, "fitmean[43]/D");
00289 layergcalib->Branch("fitmeanerr", fitmeanerr, "fitmeanerr[43]/D");
00290 layergcalib->Branch("fitsigma", fitsigma, "fitsigma[43]/D");
00291 layergcalib->Branch("chisq", chisq, "chisq[43]/D");
00292 layergcalib->Branch("fitmean_gaus", fitmean_gaus, "fitmean_gaus[43]/D");
00293 layergcalib->Branch("fitmeanerr_gaus", fitmeanerr_gaus, "fitmeanerr_gaus[43]/D");
00294 layergcalib->Branch("fitsigma_gaus", fitsigma_gaus, "fitsigma_gaus[43]/D");
00295 layergcalib->Fill();
00296 layergcalib->Write();
00297
00298 f->Close();
00299
00300 TCanvas c1("c1", "canvas", 500, 400);
00301 c1.Print((m_rootfile+".ps[").c_str());
00302 for(int i=0;i<layNo;i++)
00303 {
00304 m_laygain[i]->Draw();
00305 c1.Print((m_rootfile+".ps").c_str());
00306 }
00307 c1.Print((m_rootfile+".ps]").c_str());
00308
00309 for(int i=0;i<layNo;i++)
00310 {
00311 delete m_laygain[i];
00312 delete m_laygain_gaus[i];
00313 }
00314 }