#include <DedxCalibSaturation.h>
Inheritance diagram for DedxCalibSaturation:
Public Member Functions | |
void | AnalyseHists () |
void | AnalyseHists () |
void | BookHists () |
void | BookHists () |
void | checkSelections () |
void | checkSelections () |
DedxCalibSaturation (const std::string &name, ISvcLocator *pSvcLocator) | |
DedxCalibSaturation (const std::string &name, ISvcLocator *pSvcLocator) | |
virtual StatusCode | execute () |
virtual StatusCode | execute () |
void | FillHists () |
void | FillHists () |
void | FillHistsFromTree () |
void | FillHistsFromTree () |
void | FillTestHists () |
void | FillTestHists () |
void | FillTree () |
void | FillTree () |
virtual StatusCode | finalize () |
virtual StatusCode | finalize () |
virtual Float_t | GetChargeOffCorr (Int_t il, Int_t iw, Float_t ChargeIn, Float_t path, Float_t doca, Float_t EAngle, Float_t dipAngle, Float_t z) |
virtual Float_t | GetChargeOffCorr (Int_t il, Int_t iw, Float_t ChargeIn, Float_t path, Float_t doca, Float_t EAngle, Float_t dipAngle, Float_t z) |
Float_t | GetChargeOnCorr (Int_t il, Int_t iw, Float_t ChargeIn, Float_t path, Float_t doca, Float_t EAngle, Float_t dipAngle, Float_t z) |
Float_t | GetChargeOnCorr (Int_t il, Int_t iw, Float_t ChargeIn, Float_t path, Float_t doca, Float_t EAngle, Float_t dipAngle, Float_t z) |
virtual StatusCode | initialize () |
virtual StatusCode | initialize () |
void | ReadBetheBlochParameters () |
void | ReadBetheBlochParameters () |
void | ReadCalibdEdxParameters () |
void | ReadCalibdEdxParameters () |
void | WriteCalibdEdxParameters () |
void | WriteCalibdEdxParameters () |
void | WriteHists () |
void | WriteHists () |
Protected Attributes | |
int | calib_flag |
bool | ddgflag |
IDedxCorrecSvc * | exsvc |
IDedxCorrecSvc * | exsvc |
IMdcGeomSvc * | geosvc |
IMdcGeomSvc * | geosvc |
bool | ggsflag |
bool | layergflag |
MsgStream | log |
std::string | m_constrootfile |
int | m_correc_flag |
std::string | m_inputfile |
int | m_par_flag |
int | m_phshape_flag |
std::string | m_rootfile |
bool | wiregflag |
bool | zdepflag |
Private Attributes | |
TObjArray * | m_caliblist |
TObjArray * | m_caliblist |
vector< TH1F * > | m_hist |
vector< TH1F * > | m_hist |
vector< TH1F * > | m_histcalib |
vector< TH1F * > | m_histcalib |
TRandom * | m_rand |
TRandom * | m_rand |
TObjArray * | m_satlist |
TObjArray * | m_satlist |
|
00026 : DedxCalib(name, pSvcLocator) 00027 { 00028 log<<MSG::INFO<<"DedxCalibSaturation::DedxCalibSaturation( )..."<<endreq; 00029 }
|
|
|
|
Implements DedxCalib. |
|
Implements DedxCalib. 00201 { 00202 log<<MSG::INFO<<"DedxCalibSaturation::AnalyseHists( )..."<<endreq; 00203 00204 TF1* func; 00205 if(m_phshape_flag==0) { 00206 func = new TF1("mylan",mylan,1000,4000,4); 00207 func->SetParameters(10000.0,2010.1,300,-1.0); 00208 } else if(m_phshape_flag==1) { 00209 func = new TF1("landaun",landaun,1000,4000,3); 00210 func->SetParameters(2000.0,300.0); 00211 } else if(m_phshape_flag==2) { 00212 func = new TF1("land",Landau,1000,4000,2); 00213 func->SetParameters(2010.1,300); 00214 } else if(m_phshape_flag==3) { 00215 func = new TF1("vavi",Vavilov,1000,4000,2); 00216 func->SetParameters(3.1,0.4); 00217 } else if(m_phshape_flag==4) { 00218 func = new TF1("asym",AsymGauss,1000,4000,4); 00219 func->SetParameters(10000.0,2000.0,300.0,100.); 00220 } 00221 00222 func->SetLineWidth(2.1); 00223 func->SetLineColor(2); 00224 Double_t entry, mean, rms; 00225 Double_t binmax; 00226 Double_t lp[3]; 00227 for(int hid=0; hid<43*NumBinCostheta; hid++){ 00228 entry = m_hist[hid]->GetEntries(); 00229 mean = m_hist[hid]->GetMean(); 00230 rms = m_hist[hid]->GetRMS(); 00231 binmax = m_hist[hid]->GetMaximumBin(); 00232 lp[1] = m_hist[hid]->GetBinCenter(binmax); 00233 lp[2] = 200; 00234 lp[0] = (m_hist[hid]->GetBinContent(binmax)+m_hist[hid]->GetBinContent(binmax-1)+m_hist[hid]->GetBinContent(binmax+1))/3; 00235 00236 if(m_phshape_flag==0) { 00237 func->SetParameters(entry, mean, rms, -0.8); 00238 } else if(m_phshape_flag==1) { 00239 //func->SetParameters(entry, 0.7*mean, rms ); 00240 func->SetParameters(lp[0], lp[1], lp[2] ); 00241 } else if(m_phshape_flag==2) { 00242 func->SetParameters(0.7*mean, rms ); 00243 } else if(m_phshape_flag==3) { 00244 func->SetParameters(0.05, 0.6); 00245 } else if(m_phshape_flag==4) { 00246 func->SetParameters(entry, mean, rms, 1.0 ); 00247 } 00248 m_hist[hid]->Fit(func, "Q" ); 00249 } 00250 00251 double fitmean[NumBinCostheta], mean0; 00252 if(m_phshape_flag==0) { 00253 mean0 = m_hist[CosthetaStandard] -> GetFunction("mylan") -> GetParameter(1); 00254 } else if(m_phshape_flag==1) { 00255 mean0 = m_hist[CosthetaStandard] -> GetFunction("landaun") -> GetParameter(1); 00256 } else if(m_phshape_flag==2) { 00257 mean0 = m_hist[CosthetaStandard] -> GetFunction("land") -> GetParameter(1); 00258 } else if(m_phshape_flag==3) { 00259 mean0 = m_hist[CosthetaStandard] -> GetFunction("vavi") -> GetParameter(1); 00260 } else if(m_phshape_flag==4) { 00261 mean0 = m_hist[CosthetaStandard] -> GetFunction("asym") -> GetParameter(1); 00262 } 00263 00264 for(int lay=0; lay<43; lay++) { 00265 for(int bin=0; bin<NumBinCostheta; bin++) { 00266 int hid=lay*NumBinCostheta+bin; 00267 if(m_phshape_flag==0) { 00268 fitmean[bin] = m_hist[hid] -> GetFunction("mylan") -> GetParameter(1); 00269 } else if(m_phshape_flag==1) { 00270 fitmean[bin] = m_hist[hid] -> GetFunction("landaun") -> GetParameter(1); 00271 } else if(m_phshape_flag==2) { 00272 fitmean[bin] = m_hist[hid] -> GetFunction("land") -> GetParameter(1); 00273 } else if(m_phshape_flag==3) { 00274 fitmean[bin] = m_hist[hid] -> GetFunction("vavi") -> GetParameter(1); 00275 } else if(m_phshape_flag==4) { 00276 fitmean[bin] = m_hist[hid] -> GetFunction("asym") -> GetParameter(1); 00277 } 00278 00279 /*if(mean0) { 00280 fitmean[bin] = fitmean[bin]/mean0; 00281 } else { 00282 fitmean[bin] = fitmean[bin]/2000.0; 00283 }*/ 00284 fitmean[bin] = fitmean[bin]/NormalMean; 00285 if(fitmean[bin]>1.12) fitmean[bin] =1; 00286 m_histcalib[lay]->Fill(Costheta_step*bin+CosthetaMin+Costheta_step/2,fitmean[bin]); 00287 } 00288 } 00289 00290 TF1* parfunc; 00291 parfunc = new TF1("cheby3",ChebyP3,0.0,15.0,4); 00292 parfunc->SetParameters(10.0,1.0,3.1,0.4); 00293 parfunc->SetLineWidth(2.1); 00294 parfunc->SetLineColor(2); 00295 //cout<<"par_flag = "<<m_par_flag<<endl; 00296 for(int lay=0; lay<43; lay++) { 00297 if(m_par_flag==0) m_histcalib[lay]->Fit(parfunc,"Q"); 00298 if(m_par_flag==1) m_histcalib[lay]->Fit("pol3","Q"); 00299 } 00300 }
|
|
Implements DedxCalib. |
|
Implements DedxCalib. 00081 { 00082 log<<MSG::INFO<<"DedxCalibSaturation::BookHists( )..."<<endreq; 00083 log<<MSG::DEBUG<<" bookhist...1" <<endreq; 00084 m_satlist = new TObjArray(0); 00085 m_caliblist = new TObjArray(0); 00086 log<<MSG::DEBUG<<" bookhist...2" <<endreq; 00087 m_rand = new TRandom(); 00088 std::string hlabel; 00089 log<<MSG::DEBUG<<" bookhist...3" <<endreq; 00090 00091 00092 for(int jj=0; jj<43; jj++) { 00093 for(int mm=0; mm<NumBinCostheta; mm++) { 00094 std::ostringstream strout; 00095 strout << "Gain_Layer_" << jj <<"_Costheta_"<<mm; 00096 hlabel = strout.str(); 00097 TH1F *sat = new TH1F(hlabel.c_str(), hlabel.c_str(), NumHistBins,MinHistValue,MaxHistValue); 00098 // wireg->Print("base"); 00099 sat->GetXaxis()->SetTitle("dE/dx(eV/1.5cm)"); 00100 sat->GetYaxis()->SetTitle("entries/(40eV/1.5cm)"); 00101 m_hist.push_back(sat); 00102 m_satlist -> Add(sat); 00103 } 00104 std::ostringstream strout2; 00105 strout2 << "Saturation_Layer_" << jj ; 00106 hlabel = strout2.str(); 00107 TH1F *calib = new TH1F(hlabel.c_str(), hlabel.c_str(), NumBinCostheta, CosthetaMin, CosthetaMax); 00108 // wireg->Print("base"); 00109 calib->GetXaxis()->SetTitle("cos#theta"); 00110 calib->GetYaxis()->SetTitle("normalized dE/dx"); 00111 m_histcalib.push_back(calib); 00112 m_caliblist -> Add(calib); 00113 00114 } 00115 log<<MSG::DEBUG<<" m_hist size..." <<m_hist.size()<<endreq; 00116 }
|
|
|
|
00121 { 00122 log<<MSG::INFO<<"DedxCalib::checkSelections()...."<<endreq; 00123 }
|
|
|
|
00075 { 00076 00077 this->FillTree(); 00078 00079 this->FillHists(); 00080 return StatusCode::SUCCESS; 00081 }
|
|
Implements DedxCalib. |
|
Implements DedxCalib. 00119 { 00120 log<<MSG::INFO<<"DedxCalibSaturation::FillHists( )..."<<endreq; 00121 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader"); 00122 if (!eventHeader) { 00123 log << MSG::FATAL << "Could not find Event Header" << endreq; 00124 return; 00125 } 00126 int event = eventHeader->eventNumber(); 00127 int runNO = eventHeader->runNumber(); 00128 int typFlag; //data type flag: 0:MC data, 1: 2150, 2:2200 00129 if(runNO<0) typFlag =0; 00130 else if(runNO>=0 && runNO<5459) typFlag =1; 00131 else if(runNO>=5459 && runNO<8093) typFlag =2; 00132 else if(runNO>=8093) typFlag =3; 00133 00134 Identifier mdcid; 00135 SmartDataPtr<RecMdcTrackCol> trkCol(eventSvc(),"/Event/Recon/RecMdcTrackCol"); 00136 if (!trkCol) { 00137 log << MSG::WARNING << "Could not find RecMdcTrackCol" << endreq; 00138 return; 00139 } 00140 int lr, layid, hid,localwid; 00141 double adc_raw, costhe, zhit, ddl, ddr, dd, ph, sinenta, pathlength, sintheta, rphi_path; 00142 double m_driftd, driftd; 00143 log << MSG::INFO << "Begin to check RecMdcTrackCol"<<endreq; 00144 RecMdcTrackCol::iterator iter_trk = trkCol->begin(); 00145 for( ; iter_trk != trkCol->end(); iter_trk++) { 00146 HepVector a = (*iter_trk)->helix(); 00147 HepSymMatrix tkErrM = (*iter_trk)->err(); 00148 double m_d0E = tkErrM[0][0]; 00149 double m_phi0E = tkErrM[1][1]; 00150 double m_cpaE = tkErrM[2][2]; 00151 double m_z0E = tkErrM[3][3]; 00152 HepPoint3D IP(0,0,0); 00153 Dedx_Helix exhel(IP, a); 00154 float m_Pt = 1.0/fabs( a(3) ); 00155 float m_P = m_Pt*sqrt(1 + a(5)*a(5)); 00156 if(m_P>1.95||m_P<1.8) continue; 00157 00158 log << MSG::DEBUG << "hitList of this track:" << endreq; 00159 HitRefVec gothits = (*iter_trk)->getVecHits(); 00160 HitRefVec::iterator it_gothit = gothits.begin(); 00161 for( ; it_gothit != gothits.end(); it_gothit++) { 00162 mdcid = (*it_gothit)->getMdcId(); 00163 layid = MdcID::layer(mdcid); 00164 localwid = MdcID::wire(mdcid); 00165 adc_raw = (*it_gothit)->getAdc()*1000000; 00166 zhit = (*it_gothit)->getZhit(); 00167 lr = (*it_gothit)->getFlagLR(); 00168 if(lr == 0 || lr == 2) m_driftd = (*it_gothit)->getDriftDistLeft(); 00169 if(lr == 1 ) m_driftd = (*it_gothit)->getDriftDistRight(); 00170 driftd = abs(m_driftd); 00171 dd = fabs((*it_gothit)->getDoca()); 00172 sinenta = sin((*it_gothit)->getEntra()); 00173 //costhe = sin(atan(a[4])); 00174 costhe = abs(cos(M_PI_2-atan(a[4]))); 00175 sintheta = sin(M_PI_2-atan(a[4])); 00176 int ntpFlag = 1; 00177 ph = exsvc->StandardCorrec(typFlag, ntpFlag, runNO, exhel, mdcid, adc_raw, dd, sinenta, zhit, costhe ); 00178 pathlength = exsvc->PathL(ntpFlag, exhel, layid, localwid, zhit); 00179 00180 //ph = exsvc->StandardCorrec( exhel, mdcid, adc_raw, dd, sinenta, zhit, costhe ); 00181 //pathlength = exsvc->PathLCosmic(exhel, layid, localwid, zhit,m_z0E); 00182 rphi_path = pathlength * sintheta; 00183 if(runNO<=5911 && runNO>=5854 && layid ==14) continue; 00184 if(layid <8) 00185 {if( adc_raw<MinADCValuecut || adc_raw>MaxADCValuecut || rphi_path>RPhi_PathMaxCut || rphi_path<Iner_RPhi_PathMinCut || driftd>Iner_DriftDistCut) continue;} 00186 else if(layid >= 8) 00187 {if(adc_raw<MinADCValuecut || adc_raw>MaxADCValuecut || rphi_path>RPhi_PathMaxCut || rphi_path<Out_RPhi_PathMinCut || driftd>Out_DriftDistCut) continue;} 00188 if(costhe > CosthetaMax) continue; 00189 for(int kk=0; kk<NumBinCostheta; kk++) { 00190 if((costhe>=CosthetaMin+Costheta_step*kk)&&(costhe<CosthetaMin+Costheta_step*(kk+1))) { 00191 hid=layid*NumBinCostheta+kk; 00192 } 00193 } 00194 m_hist[hid]->Fill(ph); 00195 } 00196 } 00197 }
|
|
Implements DedxCalib. |
|
Implements DedxCalib. 00303 { 00304 log<<MSG::INFO<<"DedxCalibSaturation::FillHistsFromTree( )..."<<endreq; 00305 }
|
|
|
|
00117 { 00118 log<<MSG::INFO<<"DedxCalib::FillTestHists()...."<<endreq; 00119 }
|
|
|
|
00105 { 00106 log<<MSG::INFO<<"DedxCalib::FillTree()...."<<endreq; 00107 }
|
|
|
|
00084 { 00085 00086 log << MSG::INFO << "DedxCalib finalize() ..." << endreq; 00087 this->AnalyseHists(); 00088 00089 this->WriteCalibdEdxParameters(); 00090 this->WriteHists(); 00091 00092 this->FillTestHists(); 00093 00094 00095 00096 return StatusCode::SUCCESS; 00097 }
|
|
Implements DedxCalib. |
|
Implements DedxCalib. 00034 { 00035 log<<MSG::INFO<<"DedxCalibSaturation::GetChargeOffCorr( )..."<<endreq; 00036 return 1.0; 00037 }
|
|
|
|
00127 { 00128 00129 }
|
|
|
|
00035 { 00036 log << MSG::INFO << "DedxCalib initialze() ..." << endreq; 00037 /*StatusCode scint = Service::initialize(); 00038 if( scint.isFailure() ) return scint; 00039 IIncidentSvc* incsvc; 00040 scint = service("IncidentSvc", incsvc); 00041 int priority = 100; 00042 if( sc.isSuccess() ){ 00043 incsvc -> addListener(this, "BeginEvent", priority); 00044 //incsvc -> addListener(this, "NewRun", priority); 00045 }*/ 00046 00047 StatusCode sc = service("MdcGeomSvc", geosvc); 00048 if (sc == StatusCode::SUCCESS) { 00049 log << MSG::INFO <<"MdcGeomSvc is running"<<endl; 00050 } else { 00051 log << MSG::ERROR <<"MdcGeomSvc is failed"<<endl; 00052 return StatusCode::SUCCESS; 00053 } 00054 00055 StatusCode scex = service("DedxCorrecSvc", exsvc); 00056 if (scex == StatusCode::SUCCESS) { 00057 log << MSG::INFO <<"Hi, DedxCorrectSvc is running"<<endl; 00058 } else { 00059 log << MSG::ERROR <<"DedxCorrectSvc is failed"<<endl; 00060 return StatusCode::SUCCESS; 00061 } 00062 00063 exsvc->set_flag( calib_flag ); 00064 this->checkSelections(); 00065 log << MSG::INFO <<"DedxCalib: read correction parameters"<<endreq; 00066 this->ReadCalibdEdxParameters(); 00067 this->FillHistsFromTree(); 00068 00069 this->BookHists(); 00070 00071 return StatusCode::SUCCESS; 00072 }
|
|
|
|
|
|
|
|
00109 { 00110 log<<MSG::INFO<<"DedxCalib::ReadCalibdEdxParameters()...."<<endreq; 00111 }
|
|
Reimplemented from DedxCalib. |
|
Reimplemented from DedxCalib. 00308 { 00309 log<<MSG::INFO<<"DedxCalibSaturation::WriteCalibdEdxParameters( )..."<<endreq; 00310 }
|
|
Implements DedxCalib. |
|
Implements DedxCalib. 00041 { 00042 log<<MSG::INFO<<"DedxCalibSaturation::WriteHists( )..."<<endreq; 00043 double mean[4][43]; 00044 for(int lay=0; lay<43; lay++) { 00045 if(m_par_flag==0) { 00046 mean[0][lay] = m_histcalib[lay] -> GetFunction("cheby3") -> GetParameter(0); 00047 mean[1][lay] = m_histcalib[lay] -> GetFunction("cheby3") -> GetParameter(1); 00048 mean[2][lay] = m_histcalib[lay] -> GetFunction("cheby3") -> GetParameter(2); 00049 mean[3][lay] = m_histcalib[lay] -> GetFunction("cheby3") -> GetParameter(3); 00050 } 00051 if(m_par_flag==1) { 00052 mean[0][lay] = m_histcalib[lay] -> GetFunction("pol3") -> GetParameter(0); 00053 mean[1][lay] = m_histcalib[lay] -> GetFunction("pol3") -> GetParameter(1); 00054 mean[2][lay] = m_histcalib[lay] -> GetFunction("pol3") -> GetParameter(2); 00055 mean[3][lay] = m_histcalib[lay] -> GetFunction("pol3") -> GetParameter(3); 00056 } 00057 } 00058 00059 TTree* ggs = new TTree("ggscalib", "ggscalib"); 00060 ggs -> Branch("ggs0", &mean[0], "ggs[43]/D"); 00061 ggs -> Branch("ggs1", &mean[1], "ggs[43]/D"); 00062 ggs -> Branch("ggs2", &mean[2], "ggs[43]/D"); 00063 ggs -> Branch("ggs3", &mean[3], "ggs[43]/D"); 00064 ggs -> Fill(); 00065 //TFile fhist("/ihepbatch/d09/xcao/calibconst_634/ggscalib.root", "recreate"); 00066 TFile fhist("/ihepbatch/besd09/xcao/calibconst_temp/ggscalib.root", "recreate"); 00067 00068 ggs -> Write(); 00069 m_satlist -> Write(); 00070 m_caliblist -> Write(); 00071 fhist.Close(); 00072 for(int i=0; i<43*NumBinCostheta; i++) delete m_hist[i]; 00073 for(int i=0; i<43; i++) delete m_histcalib[i]; 00074 delete m_satlist; 00075 delete m_caliblist; 00076 delete ggs; 00077 delete m_rand; 00078 }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|