BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
DedxCalibDocaEAng.cxx
Go to the documentation of this file.
1#include "GaudiKernel/MsgStream.h"
2
3#include <sstream>
4#include <string>
5#include "TTree.h"
6#include "TH1F.h"
7#include "TH2F.h"
8#include "TF1.h"
9#include "TFile.h"
10#include "TStyle.h"
11#include "TCanvas.h"
12
14DECLARE_COMPONENT(DedxCalibDocaEAng)
15
16using namespace std;
18
19DedxCalibDocaEAng::DedxCalibDocaEAng(const std::string& name, ISvcLocator* pSvcLocator): DedxCalib(name, pSvcLocator){
20 declareProperty("Debug", m_debug=false);
21 declareProperty("Idoca", m_idoca_test=1);
22 declareProperty("Ieang", m_ieang_test=1);
23 declareProperty("PMin", pMin=0.2);
24 declareProperty("PMax", pMax=2.3);
25 declareProperty("Truncate", m_truncate=0.7);
26}
27
29{
30 MsgStream log(msgSvc(), name());
31 log<<MSG::INFO<<"DedxCalibDocaEAng::BookHists()"<<endreq;
32
34
35 m_docaeangle = new TH1F**[NumSlicesDoca];
36 stringstream hlabel;
37 for(int i=0;i<NumSlicesDoca;i++)
38 {
39 m_docaeangle[i] = new TH1F*[NumSlicesEAng];
40 for(int j=0;j<NumSlicesEAng;j++)
41 {
42 hlabel.str("");
43 hlabel<<"dEdx_doca"<<i<<"_eang"<<j;
44 m_docaeangle[i][j] = new TH1F(hlabel.str().c_str(), hlabel.str().c_str(), NumHistBins, MinHistValue, MaxHistValue);
45 }
46 }
47 hlabel.str("");
48 hlabel<<"dEdxVsDocaEAng";
49 m_DocaEAngAverdE = new TH2F(hlabel.str().c_str(),"dEdx vs Doca & EAng",NumSlicesDoca,DocaMin,DocaMax,NumSlicesEAng,PhiMin,PhiMax);
50}
51
53{
54 MsgStream log(msgSvc(), name());
55 log<<MSG::INFO<<"DedxCalibDocaEAng::FillHists()"<<endreq;
56
57 TFile* f, *f_test(0);
58 TTree* n103, *t_test(0);
59 string runlist;
60
61 double doca_test(0), eang_test(0);
62 int ndedxhit=0;
63 double dEdx[100]={0},pathlength[100]={0},wid[100]={0},layid[100]={0},dd_in[100]={0},eangle[100]={0},zhit[100]={0};
64 float runNO=0,evtNO=0,runFlag=0,costheta=0,tes=0,charge=0,recalg=0,ptrk=0,chidedx=0;
65 cout<<"m_recFileVector.size()= "<<m_recFileVector.size()<<endl;
66 if(m_debug){ f_test = new TFile("doca_eangle_test.root","recreate");
67 t_test = new TTree("inftest", "information for test doca eangle calibration");
68 t_test->Branch("idoca", &m_idoca_test, "idoca/I");
69 t_test->Branch("ieang", &m_ieang_test, "ieang/I");
70 t_test->Branch("doca", &doca_test, "doca/D");
71 t_test->Branch("eang", &eang_test, "eang/D");
72 }
73
74 for(unsigned int i=0; i<m_recFileVector.size(); i++){
75 runlist = m_recFileVector[i];
76 cout<<"runlist: "<<runlist.c_str()<<endl;
77 f = new TFile(runlist.c_str());
78 n103 = (TTree*)f->Get("n103");
79 if(f==0){ cout << "the file is empty" << endl; continue;}
80 if(n103==0){ cout << "the tree is empty" << endl; continue;}
81 n103-> SetBranchAddress("ndedxhit", &ndedxhit);
82 n103-> SetBranchAddress("dEdx_hit",dEdx);
83 n103-> SetBranchAddress("pathlength_hit",pathlength);
84 n103-> SetBranchAddress("wid_hit",wid);
85 n103-> SetBranchAddress("layid_hit",layid);
86 n103-> SetBranchAddress("dd_in_hit",dd_in);
87 n103-> SetBranchAddress("eangle_hit",eangle);
88 n103-> SetBranchAddress("zhit_hit",zhit);
89 n103-> SetBranchAddress("runNO",&runNO);
90 n103-> SetBranchAddress("evtNO",&evtNO);
91 n103-> SetBranchAddress("runFlag",&runFlag);
92 n103-> SetBranchAddress("costheta",&costheta);
93 n103-> SetBranchAddress("t0",&tes);
94 n103-> SetBranchAddress("charge",&charge);
95 n103-> SetBranchAddress("recalg",&recalg);
96 n103-> SetBranchAddress("ptrk",&ptrk);
97 n103-> SetBranchAddress("chidedx_e",&chidedx);
98
99 vector<double> m_phlist_test(0);
100 for(int j=0;j<n103->GetEntries();j++){
101 n103->GetEntry(j);
102 if(ptrk>pMax || ptrk<pMin) continue;
103 if(tes>1400) continue;
104 // begin of save truncated hits information for EACH track
105 // m_phlist_test.clear();
106 // for(int ktest=0; ktest<ndedxhit; ktest++){
107 // if(layid[ktest]<4) continue;
108 // m_phlist_test.push_back(dEdx[ktest]);
109 // }
110 // sort(m_phlist_test.begin(), m_phlist_test.end());
111 // if(m_debug){
112 // for(unsigned int i=0; i<m_phlist_test.size(); i++)
113 // cout << "i " << i << " " << m_phlist_test[i] << endl;
114 // }
115 // unsigned int nsample = (unsigned int)(m_phlist_test.size()*m_truncate);
116 // double m_trunc_value = m_phlist_test[nsample];
117 // if(m_debug) cout << "trunctate value " << m_trunc_value << endl;
118 for(int k=0;k<ndedxhit;k++){
119 // if(dEdx[k]>m_trunc_value) continue;
120 dd_in[k] = dd_in[k]/doca_norm[int(layid[k])];
121 if(eangle[k] >0.25) eangle[k] = eangle[k] -0.5;
122 else if(eangle[k] <-0.25) eangle[k] = eangle[k] +0.5;
123 if(eangle[k]>PhiMin && eangle[k]<PhiMax && dd_in[k]>DocaMin && dd_in[k]<DocaMax)
124 {
125 if(layid[k] <4) continue;
126// else if(layid[k] <8)
127// {if(dEdx[k]<MinADCValuecut || dEdx[k]>MaxADCValuecut || pathlength[k]>RPhi_PathMaxCut || pathlength[k]<Iner_RPhi_PathMinCut || driftdist[k]>Iner_DriftDistCut) continue;}
128// else
129// {if(dEdx[k]<MinADCValuecut || dEdx[k]>MaxADCValuecut || pathlength[k]>RPhi_PathMaxCut || pathlength[k]<Out_RPhi_PathMinCut || driftdist[k]>Out_DriftDistCut) continue;}
130 int iDoca =(Int_t)floor((dd_in[k]-DocaMin)/StepDoca);
131 int iEAng = 0;
132 for(int i =0; i<14; i++)
133 {
134 if(eangle[k] < Eangle_value_cut[i] || eangle[k] >= Eangle_value_cut[i+1]) continue;
135 else if(eangle[k] >= Eangle_value_cut[i] && eangle[k] < Eangle_value_cut[i+1])
136 {
137 for(int m =0; m<i; m++)
138 {
139 iEAng += Eangle_cut_bin[m];
140 }
141 double eangle_bin_cut_step = (Eangle_value_cut[i+1] - Eangle_value_cut[i])/Eangle_cut_bin[i];
142 int temp_bin = int((eangle[k]-Eangle_value_cut[i])/eangle_bin_cut_step);
143 iEAng += temp_bin;
144 }
145 }
146 if(m_debug && (iDoca==m_idoca_test) && (iEAng==m_ieang_test)){
147 doca_test = dd_in[k];
148 eang_test = eangle[k];
149 t_test->Fill();
150 cout << "dedx before cor: " << dEdx[k]*1.5/pathlength[k] << endl;
151 }
152 dEdx[k] = exsvc->StandardHitCorrec(0,runFlag,2,runNO,evtNO,pathlength[k],wid[k],layid[k],dEdx[k],(dd_in[k])*(doca_norm[int(layid[k])]),eangle[k],costheta);
153 dEdx[k] = exsvc->StandardTrackCorrec(0, runFlag, 2, runNO, evtNO, dEdx[k], costheta, tes);
154 m_docaeangle[iDoca][iEAng]->Fill(dEdx[k]);
155 if(m_debug && (iDoca==m_idoca_test) && (iEAng==m_ieang_test)){
156 cout << "dedx after cor: " << dEdx[k] << endl << endl;
157 }
158 } // end of if(eangle, doca) are normal
159 } // end of loop in k, hit list
160 } // end of loop in j, entry list
161 f->Close(); f = 0; n103 = 0;
162 }// end of loop in i, filelist
163 if(m_debug){ f_test->cd(); t_test->Write(); f_test->Close(); }
164} // end of FillHists()
165
167{
168 MsgStream log(msgSvc(), name());
169 log<<MSG::INFO<<"DedxCalibDocaEAng::AnalyseHists()"<<endreq;
170
171 TF1* func;
172 Double_t entry=0,mean=0,rms=0;
173 Double_t binmax=0;
174 Double_t lp[3]={0};
175 gStyle->SetOptFit(1111);
176 for(Int_t id=0;id<NumSlicesDoca;id++)
177 {
178 for(Int_t ip=0;ip<NumSlicesEAng;ip++)
179 {
180 entry = m_docaeangle[id][ip]->GetEntries();
181 if(entry<10) continue;
182 mean = m_docaeangle[id][ip]->GetMean();
183 rms = m_docaeangle[id][ip]->GetRMS();
184 binmax = m_docaeangle[id][ip]->GetMaximumBin();
185 lp[1] = m_docaeangle[id][ip]->GetBinCenter(binmax);
186 lp[2] = 200;
187 lp[0] = (m_docaeangle[id][ip]->GetBinContent(binmax)+m_docaeangle[id][ip]->GetBinContent(binmax-1)+m_docaeangle[id][ip]->GetBinContent(binmax+1))/3;
188
189 if(m_phShapeFlag==0)
190 {
191 func = new TF1("mylan",mylan,10,3000,4);
192 func->SetParameters(entry, mean, rms, -0.8);
193 func->SetParLimits(0, 0, entry);
194 func->SetParLimits(1, 0, mean);
195 func->SetParLimits(2, 10, rms);
196 func->SetParLimits(3, -1, 0);
197 func->SetParameters(entry/36., mean-200, rms/2.3, -0.5);
198 }
199 else if(m_phShapeFlag==1)
200 {
201 func = new TF1("landaun",landaun,10,3000,3);
202 func->SetParameters(lp[0], lp[1], lp[2] );
203 }
204 else if(m_phShapeFlag==2)
205 {
206 func = new TF1("Landau",Landau,10,3000,2);
207 func->SetParameters(0.7*mean, rms );
208 }
209 else if(m_phShapeFlag==3)
210 {
211 func = new TF1("Vavilov",Vavilov,10,3000,2);
212 func->SetParameters(0.05, 0.6);
213 }
214 else if(m_phShapeFlag==4)
215 {
216 func = new TF1("AsymGauss",AsymGauss,10,3000,4);
217 func->SetParameters(entry, mean, rms, 1.0 );
218 }
219 func->SetLineWidth(2.1);
220 func->SetLineColor(2);
221
222 m_docaeangle[id][ip]->Fit(func, "rq");
223 }
224 }
225}
226
228{
229 MsgStream log(msgSvc(), name());
230 log<<MSG::INFO<<"DedxCalibDocaEAng::WriteHists()"<<endreq;
231
232 const int totalNo = NumSlicesDoca*NumSlicesEAng;
233 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};
234 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};
235 double Id_doca[totalNo]={0},Ip_eangle[totalNo]={0},doca[totalNo]={0},eangle[totalNo]={0},chisq[totalNo]={0};
236
237 double Norm(0), count(0);
238 for(Int_t id=0;id<NumSlicesDoca;id++)
239 {
240 for(Int_t ip=0;ip<NumSlicesEAng;ip++)
241 {
242 Id_doca[id*NumSlicesEAng+ip] = id;
243 Ip_eangle[id*NumSlicesEAng+ip] = ip;
244 doca[id*NumSlicesEAng+ip] = (id+0.5)*StepDoca + DocaMin;
245 int iE=0,i=0;
246 for(;i<14;i++)
247 {
248 if(ip>=iE && ip<iE+Eangle_cut_bin[i]) break;
249 iE=iE+Eangle_cut_bin[i];
250 }
251 //cout<<"iE= "<<iE<<" Eangle_cut_bin["<<i<<"]= "<<Eangle_cut_bin[i]<<endl;
252 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];
253 //cout<<"eangle["<<id<<"]["<<ip<<"]= "<<eangle[id*NumSlicesEAng+ip]<<endl;
254
255 Out_entryNo[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetEntries();
256 Out_mean[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetMean();
257 Out_sigma[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetRMS();
258 //cout<<"Out_entryNo["<<id<<"]["<<ip<<"]= "<<Out_entryNo[id*NumSlicesEAng+ip]<<" "<<Out_mean[id*NumSlicesEAng+ip]<<" "<<Out_sigma[id*NumSlicesEAng+ip]<<endl;
259 if(Out_entryNo[id*NumSlicesEAng+ip]<10) continue;
260
261 if(m_phShapeFlag==0)
262 {
263 Out_fitmean[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("mylan")->GetParameter(1);
264 Out_fitmeanerr[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("mylan")->GetParError(1);
265 Out_fitsigma[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("mylan")->GetParameter(2);
266 //cout<<Out_fitmean[id*NumSlicesEAng+ip]<<" "<<Out_fitmeanerr[id*NumSlicesEAng+ip]<<" "<<Out_fitsigma[id*NumSlicesEAng+ip]<<endl;
267 chisq[id*NumSlicesEAng+ip] = (m_docaeangle[id][ip]->GetFunction("mylan")-> GetChisquare())/(m_docaeangle[id][ip]->GetFunction("mylan")-> GetNDF());
268 }
269 if(m_phShapeFlag==1)
270 {
271 Out_fitmean[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("landaun")->GetParameter(1);
272 Out_fitmeanerr[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("landaun")->GetParError(1);
273 Out_fitsigma[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("landaun")->GetParameter(2);
274 chisq[id*NumSlicesEAng+ip] = (m_docaeangle[id][ip]->GetFunction("landaun")-> GetChisquare())/(m_docaeangle[id][ip]->GetFunction("landaun")-> GetNDF());
275 }
276 if(m_phShapeFlag==2)
277 {
278 Out_fitmean[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("Landau")->GetParameter(1);
279 Out_fitmeanerr[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("Landau")->GetParError(1);
280 Out_fitsigma[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("Landau")->GetParameter(2);
281 chisq[id*NumSlicesEAng+ip] = (m_docaeangle[id][ip]->GetFunction("Landau")-> GetChisquare())/(m_docaeangle[id][ip]->GetFunction("Landau")-> GetNDF());
282 }
283 if(m_phShapeFlag==3)
284 {
285 Out_fitmean[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("Vavilov")->GetParameter(1);
286 Out_fitmeanerr[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("Vavilov")->GetParError(1);
287 Out_fitsigma[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("Vavilov")->GetParameter(2);
288 chisq[id*NumSlicesEAng+ip] = (m_docaeangle[id][ip]->GetFunction("Vavilov")-> GetChisquare())/(m_docaeangle[id][ip]->GetFunction("Vavilov")-> GetNDF());
289 }
290 if(m_phShapeFlag==4)
291 {
292 Out_fitmean[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("AsymGauss")->GetParameter(1);
293 Out_fitmeanerr[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("AsymGauss")->GetParError(1);
294 Out_fitsigma[id*NumSlicesEAng+ip] = m_docaeangle[id][ip]->GetFunction("AsymGauss")->GetParameter(2);
295 chisq[id*NumSlicesEAng+ip] = (m_docaeangle[id][ip]->GetFunction("AsymGauss")-> GetChisquare())/(m_docaeangle[id][ip]->GetFunction("AsymGauss")-> GetNDF());
296 }
297
298 m_DocaEAngAverdE -> SetBinContent(id+1, ip+1, Out_fitmean[id*NumSlicesEAng+ip]);
299 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){
300 Norm += (Out_fitmean[id*NumSlicesEAng+ip])*(Out_entryNo[id*NumSlicesEAng+ip]);
301 count += Out_entryNo[id*NumSlicesEAng+ip];
302 }
303 } // end of loop in ip
304 } // end of loop in id
305 Norm = Norm/count;
306 cout<<"count= "<<count<<" average value of doca eangle bins: "<<Norm<<endl;
307 for(Int_t i=0;i<totalNo;i++)
308 {
309 if(Out_entryNo[i]>10 && Out_fitmean[i]>0 && Out_fitsigma[i]>0) Out_gain[i]=Out_fitmean[i]/Norm;
310 else Out_gain[i] = 0;
311 }
312
313 log<<MSG::INFO<<"begin generating root file!!! "<<endreq;
314 TFile* f = new TFile(m_rootfile.c_str(), "recreate");
315 for(Int_t id=0;id<NumSlicesDoca;id++)
316 {
317 for(Int_t ip=0;ip<NumSlicesEAng;ip++)
318 {
319 m_docaeangle[id][ip]->Write();
320 }
321 }
322 m_DocaEAngAverdE->Write();
323
324 TTree* ddgcalib = new TTree("ddgcalib", "ddgcalib");
325 ddgcalib -> Branch("Out_hits", Out_entryNo, "Out_entryNo[1600]/D");
326 ddgcalib -> Branch("Out_mean", Out_mean, "Out_mean[1600]/D");
327 ddgcalib -> Branch("Out_sigma", Out_sigma, "Out_sigma[1600]/D");
328 ddgcalib -> Branch("Out_fitmean", Out_fitmean, "Out_fitmean[1600]/D");
329 ddgcalib -> Branch("Out_fitmeanerr", Out_fitmeanerr, "Out_fitmeanerr[1600]/D");
330 ddgcalib -> Branch("Out_chi", Out_fitsigma, "Out_fitsigma[1600]/D");
331 ddgcalib -> Branch("Out_gain", Out_gain, "Out_gain[1600]/D");
332 ddgcalib -> Branch("Iner_hits", Iner_entryNo, "Iner_entryNo[1600]/D");
333 ddgcalib -> Branch("Iner_mean", Iner_mean, "Iner_mean[1600]/D");
334 ddgcalib -> Branch("Iner_sigma", Iner_sigma, "Iner_sigma[1600]/D");
335 ddgcalib -> Branch("Iner_fitmean", Iner_fitmean, "Iner_fitmean[1600]/D");
336 ddgcalib -> Branch("Iner_fitmeanerr", Iner_fitmeanerr, "Iner_fitmeanerr[1600]/D");
337 ddgcalib -> Branch("Iner_chi", Iner_fitsigma, "Iner_fitsigma[1600]/D");
338 ddgcalib -> Branch("Iner_gain", Iner_gain, "Iner_gain[1600]/D");
339 ddgcalib -> Branch("Id_doca", Id_doca, "Id_doca[1600]/D");
340 ddgcalib -> Branch("Ip_eangle", Ip_eangle, "Ip_eangle[1600]/D");
341 ddgcalib -> Branch("chisq", chisq, "chisq[1600]/D");
342 ddgcalib -> Branch("doca", doca, "doca[1600]/D");
343 ddgcalib -> Branch("eangle", eangle, "eangle[1600]/D");
344 ddgcalib -> Fill();
345 ddgcalib -> Write();
346 f->Close();
347
348 TCanvas c1("c1", "canvas", 500, 400);
349 c1.Print((m_rootfile+".ps[").c_str());
350 for(Int_t id=0;id<NumSlicesDoca;id++)
351 {
352 for(Int_t ip=0;ip<NumSlicesEAng;ip++)
353 {
354 m_docaeangle[id][ip]->Draw();
355 c1.Print((m_rootfile+".ps").c_str());
356 }
357 }
358 c1.Print((m_rootfile+".ps]").c_str());
359
360 for(Int_t id=0;id<NumSlicesDoca;id++)
361 {
362 for(Int_t ip=0;ip<NumSlicesEAng;ip++)
363 {
364 delete m_docaeangle[id][ip];
365 }
366 }
367 delete m_DocaEAngAverdE;
368}
curve Branch("CurveSize",&CurveSize,"CurveSize/I")
curve Fill()
curve Write()
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
data SetBranchAddress("time",&time)
const double StepDoca
double Landau(double *x, double *par)
double Vavilov(double *x, double *par)
#define DocaMin
#define NumSlicesEAng
#define DocaMax
#define MinHistValue
double landaun(double *x, double *par)
double AsymGauss(double *x, double *par)
#define PhiMax
#define NumSlicesDoca
#define PhiMin
#define NumHistBins
#define MaxHistValue
double mylan(double *x, double *par)
DOUBLE_PRECISION count[3]
IMessageSvc * msgSvc()
DedxCalibDocaEAng(const std::string &name, ISvcLocator *pSvcLocator)
void ReadRecFileList()
Definition: DedxCalib.cxx:89
int m_phShapeFlag
Definition: DedxCalib.h:52
std::string m_rootfile
Definition: DedxCalib.h:55
IDedxCorrecSvc * exsvc
Definition: DedxCalib.h:30
vector< string > m_recFileVector
Definition: DedxCalib.h:48
virtual double StandardHitCorrec(int calib_rec_Flag, int runFlag, int ntpFlag, int runNO, int evtNO, double pathl, int wid, int layid, double adc, double dd, double eangle, double costheta) const =0
virtual double StandardTrackCorrec(int calib_rec_Flag, int typFlag, int ntpFlag, int runNO, int evtNO, double ex, double costheta, double t0) const =0
float costheta
float charge
float ptrk