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/DedxCalibWireGain.h"
00013
00014 using namespace std;
00015
00016 DedxCalibWireGain::DedxCalibWireGain(const std::string& name, ISvcLocator* pSvcLocator): DedxCalib(name, pSvcLocator){}
00017
00018 const int wireNo = 6796;
00019
00020 void DedxCalibWireGain::BookHists()
00021 {
00022 MsgStream log(msgSvc(), name());
00023 log<<MSG::INFO<<"DedxCalibWireGain::BookHists()"<<endreq;
00024
00025 ReadRecFileList();
00026
00027 m_wiregain = new TH1F*[wireNo];
00028 stringstream hlabel;
00029 for(int i=0; i<wireNo; i++)
00030 {
00031 hlabel.str("");
00032 hlabel<<"dEdx_Wire_"<<i;
00033 m_wiregain[i] = new TH1F(hlabel.str().c_str(), hlabel.str().c_str(), NumHistBins, MinHistValue, MaxHistValue);
00034 }
00035 }
00036
00037 void DedxCalibWireGain::FillHists()
00038 {
00039 MsgStream log(msgSvc(), name());
00040 log<<MSG::INFO<<"DedxCalibWireGain::FillHists()"<<endreq;
00041
00042 TFile* f;
00043 TTree* n102;
00044 string runlist;
00045
00046 double dedx=0;
00047 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;
00048 cout<<"m_recFileVector.size()= "<<m_recFileVector.size()<<endl;
00049 for(int i=0; i<m_recFileVector.size(); i++)
00050 {
00051 runlist = m_recFileVector[i];
00052 cout<<"runlist: "<<runlist.c_str()<<endl;
00053 f = new TFile(runlist.c_str());
00054 n102 = (TTree*)f->Get("n102");
00055 n102-> SetBranchAddress("adc_raw",&dedx);
00056 n102-> SetBranchAddress("path_rphi",&pathlength);
00057 n102-> SetBranchAddress("wire",&wid);
00058 n102-> SetBranchAddress("layer",&layid);
00059 n102-> SetBranchAddress("doca_in",&dd_in);
00060 n102-> SetBranchAddress("driftdist",&driftdist);
00061 n102-> SetBranchAddress("eangle",&eangle);
00062 n102-> SetBranchAddress("zhit",&zhit);
00063 n102-> SetBranchAddress("runNO",&runNO);
00064 n102-> SetBranchAddress("runFlag",&runFlag);
00065 n102-> SetBranchAddress("costheta1",&costheta);
00066 n102-> SetBranchAddress("t01",&tes);
00067 n102->SetBranchAddress("ptrk1",&ptrk);
00068
00069 for(int j=0;j<n102->GetEntries();j++)
00070 {
00071 n102->GetEntry(j);
00072 if(tes>1400) continue;
00073 if (ptrk<0.1) continue;
00074 if(layid <8)
00075 {if(dedx<MinADCValuecut || dedx>MaxADCValuecut || pathlength>RPhi_PathMaxCut || pathlength<Iner_RPhi_PathMinCut || driftdist>Iner_DriftDistCut) continue;}
00076 else
00077 {if(dedx<MinADCValuecut || dedx>MaxADCValuecut || pathlength>RPhi_PathMaxCut || pathlength<Out_RPhi_PathMinCut || driftdist>Out_DriftDistCut) continue;}
00078 dedx = exsvc->StandardHitCorrec(0,runFlag,2,runNO,pathlength,wid,layid,dedx,dd_in,eangle,costheta);
00079 dedx = exsvc->StandardTrackCorrec(0, runFlag, 2, runNO, dedx, costheta, tes);
00080 m_wiregain[(int)wid]->Fill(dedx);
00081 }
00082 }
00083 }
00084
00085 void DedxCalibWireGain::AnalyseHists()
00086 {
00087 MsgStream log(msgSvc(), name());
00088 log<<MSG::INFO<<"DedxCalibWireGain::AnalyseHists()"<<endreq;
00089
00090 TF1* func;
00091 Double_t entry=0,mean=0,rms=0;
00092 Double_t binmax=0;
00093 Double_t lp[3]={0};
00094 gStyle->SetOptFit(1111);
00095 for(int wire=0; wire<wireNo; wire++)
00096 {
00097 entry = m_wiregain[wire]->Integral();
00098 if(entry<10) continue;
00099 mean = m_wiregain[wire]->GetMean();
00100 rms = m_wiregain[wire]->GetRMS();
00101 binmax = m_wiregain[wire]->GetMaximumBin();
00102 lp[1] = m_wiregain[wire]->GetBinCenter(binmax);
00103 lp[2] = 200;
00104 lp[0] = (m_wiregain[wire]->GetBinContent(binmax)+m_wiregain[wire]->GetBinContent(binmax-1)+m_wiregain[wire]->GetBinContent(binmax+1))/3;
00105
00106 if(m_phShapeFlag==0)
00107 {
00108 func = new TF1("mylan",mylan,10,3000,4);
00109 func->SetParameters(entry, mean, rms, -0.8);
00110 func->SetParLimits(0, 0, entry);
00111 func->SetParLimits(1, 0, mean);
00112 func->SetParLimits(2, 10, rms);
00113 func->SetParLimits(3, -3, 3);
00114 }
00115 else if(m_phShapeFlag==1)
00116 {
00117 func = new TF1("landaun",landaun,10,3000,3);
00118 func->SetParameters(lp[0], lp[1], lp[2] );
00119 }
00120 else if(m_phShapeFlag==2)
00121 {
00122 func = new TF1("Landau",Landau,10,3000,2);
00123 func->SetParameters(0.7*mean, rms );
00124 }
00125 else if(m_phShapeFlag==3)
00126 {
00127 func = new TF1("Vavilov",Vavilov,10,3000,2);
00128 func->SetParameters(0.05, 0.6);
00129 }
00130 else if(m_phShapeFlag==4)
00131 {
00132 func = new TF1("AsymGauss",AsymGauss,10,3000,4);
00133 func->SetParameters(entry, mean, rms, 1.0 );
00134 }
00135 func->SetLineWidth(2.1);
00136 func->SetLineColor(2);
00137
00138 m_wiregain[wire]->Fit(func, "rq", "", mean-rms, mean+rms/2 );
00139 }
00140 }
00141
00142 void DedxCalibWireGain::WriteHists()
00143 {
00144 MsgStream log(msgSvc(), name());
00145 log<<MSG::INFO<<"DedxCalibWireGain::WriteHists()"<<endreq;
00146
00147
00148 double entryNo[wireNo]={0},mean[wireNo]={0},sigma[wireNo]={0},proper[wireNo]={0},fitmean[wireNo]={0},fitmeanerr[wireNo]={0},fitsigma[wireNo]={0},wireg[wireNo]={0},wire[wireNo]={0},chisq[wireNo]={0};
00149 for(int i=0; i<wireNo; i++)
00150 {
00151 entryNo[i] = m_wiregain[i]->Integral();
00152 mean[i] = m_wiregain[i]->GetMean();
00153 sigma[i] = m_wiregain[i]->GetRMS();
00154 proper[i] = m_wiregain[i]->GetBinCenter(m_wiregain[i]->GetMaximumBin());
00155 if(entryNo[i]<10) continue;
00156
00157 if(m_phShapeFlag==0)
00158 {
00159 fitmean[i] = m_wiregain[i]->GetFunction("mylan")->GetParameter(1);
00160 fitmeanerr[i] = m_wiregain[i]->GetFunction("mylan")->GetParError(1);
00161 fitsigma[i] = m_wiregain[i]->GetFunction("mylan")->GetParameter(2);
00162 chisq[i] = (m_wiregain[i]->GetFunction("mylan")-> GetChisquare())/(m_wiregain[i]->GetFunction("mylan")-> GetNDF());
00163 }
00164 else if(m_phShapeFlag==1)
00165 {
00166 fitmean[i] = m_wiregain[i]->GetFunction("landaun")->GetParameter(1);
00167 fitmeanerr[i] = m_wiregain[i]->GetFunction("landaun")->GetParError(1);
00168 fitsigma[i] = m_wiregain[i]->GetFunction("landaun")->GetParameter(2);
00169 chisq[i] = (m_wiregain[i]->GetFunction("landaun")-> GetChisquare())/(m_wiregain[i]->GetFunction("landaun")-> GetNDF());
00170 }
00171 else if(m_phShapeFlag==2)
00172 {
00173 fitmean[i] = m_wiregain[i]->GetFunction("Landau")->GetParameter(1);
00174 fitmeanerr[i] = m_wiregain[i]->GetFunction("Landau")->GetParError(1);
00175 fitsigma[i] = m_wiregain[i]->GetFunction("Landau")->GetParameter(2);
00176 chisq[i] = (m_wiregain[i]->GetFunction("Landau")-> GetChisquare())/(m_wiregain[i]->GetFunction("Landau")-> GetNDF());
00177 }
00178 else if(m_phShapeFlag==3)
00179 {
00180 fitmean[i] = m_wiregain[i]->GetFunction("Vavilov")->GetParameter(1);
00181 fitmeanerr[i] = m_wiregain[i]->GetFunction("Vavilov")->GetParError(1);
00182 fitsigma[i] = m_wiregain[i]->GetFunction("Vavilov")->GetParameter(2);
00183 chisq[i] = (m_wiregain[i]->GetFunction("Vavilov")-> GetChisquare())/(m_wiregain[i]->GetFunction("Vavilov")-> GetNDF());
00184 }
00185 else if(m_phShapeFlag==4)
00186 {
00187 fitmean[i] = m_wiregain[i]->GetFunction("AsymGauss")->GetParameter(1);
00188 fitmeanerr[i] = m_wiregain[i]->GetFunction("AsymGauss")->GetParError(1);
00189 fitsigma[i] = m_wiregain[i]->GetFunction("AsymGauss")->GetParameter(2);
00190 chisq[i] = (m_wiregain[i]->GetFunction("AsymGauss")-> GetChisquare())/(m_wiregain[i]->GetFunction("AsymGauss")-> GetNDF());
00191 }
00192
00193 }
00194
00195 double Norm=0, count=0;
00196 for(int i=188; i<wireNo; i++)
00197 {
00198 if(fitmean[i]>0 && fitmeanerr[i]>0 && fitmeanerr[i]<10 && fitsigma[i]>0 && fitsigma[i]<200 && entryNo[i]>1500)
00199 {
00200 Norm += fitmean[i];
00201 count++;
00202 }
00203 }
00204 Norm/=count;
00205 cout<<"count= "<<count<<" average value of good wire: "<<Norm<<endl;
00206
00207
00208 for(int i=0;i<wireNo;i++)
00209 {
00210 wireg[i] = fitmean[i]/Norm;
00211 wire[i] = i;
00212 }
00213
00214
00215 log<<MSG::INFO<<"begin generating root file!!! "<<endreq;
00216 TFile* f = new TFile(m_rootfile.c_str(), "recreate");
00217 for(int i=0;i<wireNo;i++) m_wiregain[i]->Write();
00218
00219 TTree* wiregcalib = new TTree("wiregcalib", "wiregcalib");
00220 wiregcalib->Branch("wireg", wireg, "wireg[6796]/D");
00221 wiregcalib->Branch("wire", wire, "wire[6796]/D");
00222 wiregcalib->Branch("entryNo", entryNo, "entryNo[6796]/D");
00223 wiregcalib->Branch("proper", proper, "proper[6796]/D");
00224 wiregcalib->Branch("mean", mean, "mean[6796]/D");
00225 wiregcalib->Branch("sigma", sigma, "sigma[6796]/D");
00226 wiregcalib->Branch("fitmean", fitmean, "fitmean[6796]/D");
00227 wiregcalib->Branch("fitmeanerr", fitmeanerr, "fitmeanerr[6796]/D");
00228 wiregcalib->Branch("fitsigma", fitsigma, "fitsigma[6796]/D");
00229 wiregcalib->Branch("chisq", chisq, "chisq[6796]/D");
00230 wiregcalib->Fill();
00231 wiregcalib->Write();
00232
00233 f->Close();
00234
00235 TCanvas c1("c1", "canvas", 500, 400);
00236 c1.Print((m_rootfile+".ps[").c_str());
00237 for(int i=0;i<wireNo;i++)
00238 {
00239 m_wiregain[i]->Draw();
00240 c1.Print((m_rootfile+".ps").c_str());
00241 }
00242 c1.Print((m_rootfile+".ps]").c_str());
00243
00244 for(int i=0;i<wireNo;i++) delete m_wiregain[i];
00245
00246 }