#include <DedxCalibZpos.h>
Inheritance diagram for DedxCalibZpos:
Public Member Functions | |
void | AnalyseHists () |
void | AnalyseHists () |
void | BookHists () |
void | BookHists () |
void | checkSelections () |
void | checkSelections () |
DedxCalibZpos (const std::string &name, ISvcLocator *pSvcLocator) | |
DedxCalibZpos (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 |
TObjArray * | m_list |
TObjArray * | m_list |
TRandom * | m_rand |
TRandom * | m_rand |
|
00029 : DedxCalib(name, pSvcLocator) 00030 { 00031 log<<MSG::INFO<<"DedxCalibZpos::DedxCalibZpos( )..."<<endreq; 00032 }
|
|
|
|
Implements DedxCalib. |
|
Implements DedxCalib. 00201 { 00202 log<<MSG::INFO<<"DedxCalibZpos::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*NumBinZhit; 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[NumBinZhit], mean0; 00252 if(m_phshape_flag==0) { 00253 mean0 = m_hist[748] -> GetFunction("mylan") -> GetParameter(1); 00254 } else if(m_phshape_flag==1) { 00255 mean0 = m_hist[748] -> GetFunction("landaun") -> GetParameter(1); 00256 } else if(m_phshape_flag==2) { 00257 mean0 = m_hist[748] -> GetFunction("land") -> GetParameter(1); 00258 } else if(m_phshape_flag==3) { 00259 mean0 = m_hist[748] -> GetFunction("vavi") -> GetParameter(1); 00260 } else if(m_phshape_flag==4) { 00261 mean0 = m_hist[748] -> GetFunction("asym") -> GetParameter(1); 00262 } 00263 00264 for(int lay=0; lay<43; lay++) { 00265 for(int bin=0; bin<NumBinZhit; bin++) { 00266 int hid=lay*NumBinZhit+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 m_histcalib[lay]->Fill(step_lay[lay]*bin+ZhitMin+step_lay[lay]/2,fitmean[bin]); 00285 } 00286 } 00287 00288 TF1* parfunc; 00289 parfunc = new TF1("cheby3",ChebyP3,0.0,15.0,4); 00290 parfunc->SetParameters(10.0,1.0,3.1,0.4); 00291 parfunc->SetLineWidth(2.1); 00292 parfunc->SetLineColor(2); 00293 for(int lay=0; lay<43; lay++) { 00294 if(m_par_flag==0) m_histcalib[lay]->Fit(parfunc,"Q"); 00295 if(m_par_flag==1) m_histcalib[lay]->Fit("pol3","Q"); 00296 } 00297 }
|
|
Implements DedxCalib. |
|
Implements DedxCalib. 00082 { 00083 log<<MSG::INFO<<"DedxCalibZpos::BookHists( )..."<<endreq; 00084 log<<MSG::DEBUG<<" bookhist...1" <<endreq; 00085 m_list = new TObjArray(0); 00086 m_caliblist = new TObjArray(0); 00087 log<<MSG::DEBUG<<" bookhist...2" <<endreq; 00088 m_rand = new TRandom(); 00089 std::string hlabel; 00090 log<<MSG::DEBUG<<" bookhist...3" <<endreq; 00091 00092 00093 for(int jj=0; jj<43; jj++) { 00094 for(int mm=0; mm<NumBinZhit; mm++) { 00095 std::ostringstream strout; 00096 strout << "Gain_Layer_" << jj <<"_Zpos_"<<mm; 00097 hlabel = strout.str(); 00098 TH1F *zdep = new TH1F(hlabel.c_str(), hlabel.c_str(), NumHistBins,MinHistValue,MaxHistValue); 00099 // wireg->Print("base"); 00100 zdep->GetXaxis()->SetTitle("dE/dx(eV/1.5cm)"); 00101 zdep->GetYaxis()->SetTitle("entries/(40eV/1.5cm)"); 00102 m_hist.push_back(zdep); 00103 m_list -> Add(zdep); 00104 } 00105 std::ostringstream strout2; 00106 strout2 << "Zpos_Layer_" << jj ; 00107 hlabel = strout2.str(); 00108 TH1F *calib = new TH1F(hlabel.c_str(), hlabel.c_str(), NumBinZhit, ZhitMin, Zhit_Max[jj]); 00109 // wireg->Print("base"); 00110 calib->GetXaxis()->SetTitle("z position(cm)"); 00111 calib->GetYaxis()->SetTitle("normalized dE/dx"); 00112 m_histcalib.push_back(calib); 00113 m_caliblist -> Add(calib); 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<<"DedxCalibZpos::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 for(int i=0; i<43; i++){ 00141 step_lay[i] = (Zhit_Max[i] - ZhitMin)/NumBinZhit; 00142 } 00143 int lr, layid, hid,localwid; 00144 double adc_raw, costhe, zhit, ddl, ddr, dd, ph, sinenta, pathlength, sintheta, rphi_path; 00145 double m_driftd, driftd; 00146 hid=0; 00147 log << MSG::INFO << "Begin to check RecMdcTrackCol"<<endreq; 00148 RecMdcTrackCol::iterator iter_trk = trkCol->begin(); 00149 for( ; iter_trk != trkCol->end(); iter_trk++) { 00150 HepVector a = (*iter_trk)->helix(); 00151 HepSymMatrix tkErrM = (*iter_trk)->err(); 00152 double m_d0E = tkErrM[0][0]; 00153 double m_phi0E = tkErrM[1][1]; 00154 double m_cpaE = tkErrM[2][2]; 00155 double m_z0E = tkErrM[3][3]; 00156 00157 HepPoint3D IP(0,0,0); 00158 Dedx_Helix exhel(IP, a); 00159 log << MSG::DEBUG << "hitList of this track:" << endreq; 00160 HitRefVec gothits = (*iter_trk)->getVecHits(); 00161 HitRefVec::iterator it_gothit = gothits.begin(); 00162 for( ; it_gothit != gothits.end(); it_gothit++) { 00163 mdcid = (*it_gothit)->getMdcId(); 00164 layid = MdcID::layer(mdcid); 00165 localwid = MdcID::wire(mdcid); 00166 adc_raw = (*it_gothit)->getAdc()*1000000; 00167 zhit = fabs((*it_gothit)->getZhit()); 00168 lr = (*it_gothit)->getFlagLR(); 00169 if(lr == 0 || lr == 2) m_driftd = (*it_gothit)->getDriftDistLeft(); 00170 if(lr == 1 ) m_driftd = (*it_gothit)->getDriftDistRight(); 00171 driftd = abs(m_driftd); 00172 dd = fabs((*it_gothit)->getDoca()); 00173 sinenta = sin((*it_gothit)->getEntra()); 00174 costhe = 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->PathL(exhel, layid, localwid, zhit); 00182 // pathlength = exsvc->PathLCosmic(exhel, layid, localwid, zhit,m_z0E); 00183 rphi_path = pathlength * sintheta; 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(zhit >Zhit_Max[layid]) continue; 00189 for(int kk=0; kk<NumBinZhit; kk++) { 00190 if((zhit>step_lay[layid]*kk)&&(zhit<step_lay[layid]*(kk+1))) { 00191 hid=layid*NumBinZhit+kk; 00192 } 00193 m_hist[hid]->Fill(ph); 00194 } 00195 } 00196 } 00197 }
|
|
Implements DedxCalib. |
|
Implements DedxCalib. 00301 { 00302 log<<MSG::INFO<<"DedxCalibZpos::FillHistsFromTree( )..."<<endreq; 00303 }
|
|
|
|
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. 00037 { 00038 log<<MSG::INFO<<"DedxCalibZpos::GetChargeOffCorr( )..."<<endreq; 00039 return 1.0; 00040 }
|
|
|
|
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. 00306 { 00307 log<<MSG::INFO<<"DedxCalibZpos::WriteCalibdEdxParameters( )..."<<endreq; 00308 }
|
|
Implements DedxCalib. |
|
Implements DedxCalib. 00044 { 00045 log<<MSG::INFO<<"DedxCalibZpos::WriteHists( )..."<<endreq; 00046 double mean[4][43]; 00047 for(int lay=0; lay<43; lay++) { 00048 if(m_par_flag==0) { 00049 mean[0][lay] = m_histcalib[lay] -> GetFunction("cheby3") -> GetParameter(0); 00050 mean[1][lay] = m_histcalib[lay] -> GetFunction("cheby3") -> GetParameter(1); 00051 mean[2][lay] = m_histcalib[lay] -> GetFunction("cheby3") -> GetParameter(2); 00052 mean[3][lay] = m_histcalib[lay] -> GetFunction("cheby3") -> GetParameter(3); 00053 } 00054 if(m_par_flag==1) { 00055 mean[0][lay] = m_histcalib[lay] -> GetFunction("pol3") -> GetParameter(0); 00056 mean[1][lay] = m_histcalib[lay] -> GetFunction("pol3") -> GetParameter(1); 00057 mean[2][lay] = m_histcalib[lay] -> GetFunction("pol3") -> GetParameter(2); 00058 mean[3][lay] = m_histcalib[lay] -> GetFunction("pol3") -> GetParameter(3); 00059 } 00060 } 00061 TTree* zdep = new TTree("zdepcalib", "zdepcalib"); 00062 zdep -> Branch("zdep0", &mean[0], "zdep[43]/D"); 00063 zdep -> Branch("zdep1", &mean[1], "zdep[43]/D"); 00064 zdep -> Branch("zdep2", &mean[2], "zdep[43]/D"); 00065 zdep -> Branch("zdep3", &mean[3], "zdep[43]/D"); 00066 zdep -> Fill(); 00067 //TFile fhist("/ihepbatch/d09/xcao/calibconst_634/zdepcalib.root", "recreate"); 00068 TFile fhist("/ihepbatch/besd09/xcao/calibconst_temp/zdepcalib.root", "recreate"); 00069 zdep -> Write(); 00070 m_list -> Write(); 00071 m_caliblist -> Write(); 00072 fhist.Close(); 00073 for(int i=0; i<860; i++) delete m_hist[i]; 00074 for(int i=0; i<43; i++) delete m_histcalib[i]; 00075 delete m_list; 00076 delete m_caliblist; 00077 delete zdep; 00078 delete m_rand; 00079 }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|