BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
DedxCalibMomentum.cxx
Go to the documentation of this file.
1#include "GaudiKernel/MsgStream.h"
2
3#include <sstream>
4#include <string>
5#include "cmath"
6#include "TTree.h"
7#include "TH1F.h"
8#include "TF1.h"
9#include "TFile.h"
10#include "TStyle.h"
11#include "TCanvas.h"
12
13#include "DedxCalibAlg/DedxCalibMomentum.h"
14
15#define Size 700000
16DECLARE_COMPONENT(DedxCalibMomentum)
17
18using namespace std;
19const int nbins=50;
20const double chihist_min=-10;
21const double chihist_max=10;
22const double dedxhist_min=350;
23const double dedxhist_max=750;
24const double hits_min=0;
25const double hits_max=35;
26
27DedxCalibMomentum::DedxCalibMomentum(const std::string& name, ISvcLocator* pSvcLocator): DedxCalib(name, pSvcLocator){
28 declareProperty("PMin", pMin=0.2);
29 declareProperty("PMax", pMax=1.0);
30}
31
32void DedxCalibMomentum::BookHists()
33{
34 MsgStream log(msgSvc(), name());
35 log<<MSG::INFO<<"DedxCalibMomentum::BookHists()"<<endreq;
36 bin_step = (pMax-pMin)/nbins;
37 ReadRecFileList();
38
39 m_chi = new TH1F*[nbins];
40 m_pos_chi = new TH1F*[nbins];
41 m_neg_chi = new TH1F*[nbins];
42 m_dedx = new TH1F*[nbins];
43 m_pos_dedx = new TH1F*[nbins];
44 m_neg_dedx = new TH1F*[nbins];
45 m_hits = new TH1F*[nbins];
46
47 stringstream hlabel;
48 for(int i=0;i<nbins;i++)
49 {
50 hlabel.str("");
51 hlabel<<"chi"<<i;
52 m_chi[i] = new TH1F(hlabel.str().c_str(), hlabel.str().c_str(), nbins, chihist_min, chihist_max);
53 hlabel.str("");
54 hlabel<<"pos_chi"<<i;
55 m_pos_chi[i] = new TH1F(hlabel.str().c_str(), hlabel.str().c_str(), nbins, chihist_min, chihist_max);
56 hlabel.str("");
57 hlabel<<"neg_chi"<<i;
58 m_neg_chi[i] = new TH1F(hlabel.str().c_str(), hlabel.str().c_str(), nbins, chihist_min, chihist_max);
59 hlabel.str("");
60 hlabel<<"dedx"<<i;
61 m_dedx[i] = new TH1F(hlabel.str().c_str(), hlabel.str().c_str(), nbins, dedxhist_min, dedxhist_max);
62 hlabel.str("");
63 hlabel<<"pos_dedx"<<i;
64 m_pos_dedx[i] = new TH1F(hlabel.str().c_str(), hlabel.str().c_str(), nbins, dedxhist_min, dedxhist_max);
65 hlabel.str("");
66 hlabel<<"neg_dedx"<<i;
67 m_neg_dedx[i] = new TH1F(hlabel.str().c_str(), hlabel.str().c_str(), nbins, dedxhist_min, dedxhist_max);
68 hlabel.str("");
69 hlabel<<"hits"<<i;
70 m_hits[i] = new TH1F(hlabel.str().c_str(), hlabel.str().c_str(), nbins, hits_min, hits_max);
71 }
72
73 Vec_dedx.clear();
74 Vec_ptrk.clear();
75}
76
77void DedxCalibMomentum::FillHists()
78{
79 MsgStream log(msgSvc(), name());
80 log<<MSG::INFO<<"DedxCalibMomentum::FillHists()"<<endreq;
81
82 TFile* f;
83 TTree* n103;
84 string runlist;
85
86 int ndedxhit=0;
87 double dEdx[100]={0},pathlength[100]={0},wid[100]={0},layid[100]={0},dd_in[100]={0},eangle[100]={0},zhit[100]={0};
88 double dedx=0;
89 float runNO=0,evtNO=0,runFlag=0,costheta=0,tes=0,charge=0,recalg=0,ptrk=0,chidedx=0;
90 int usedhit=0;
91 vector<double> phlist;
92 cout<<"m_recFileVector.size()= "<<m_recFileVector.size()<<endl;
93 for(int i=0; i<m_recFileVector.size(); i++)
94 {
95 runlist = m_recFileVector[i];
96 f = new TFile(runlist.c_str());
97 n103 = (TTree*)f->Get("n103");
98 n103-> SetBranchAddress("ndedxhit", &ndedxhit);
99 n103-> SetBranchAddress("dEdx_hit",dEdx);
100 n103-> SetBranchAddress("pathlength_hit",pathlength);
101 n103-> SetBranchAddress("wid_hit",wid);
102 n103-> SetBranchAddress("layid_hit",layid);
103 n103-> SetBranchAddress("dd_in_hit",dd_in);
104 n103-> SetBranchAddress("eangle_hit",eangle);
105 n103-> SetBranchAddress("zhit_hit",zhit);
106 n103-> SetBranchAddress("runNO",&runNO);
107 n103-> SetBranchAddress("evtNO",&evtNO);
108 n103-> SetBranchAddress("runFlag",&runFlag);
109 n103-> SetBranchAddress("costheta",&costheta);
110 n103-> SetBranchAddress("t0",&tes);
111 n103-> SetBranchAddress("charge",&charge);
112 n103-> SetBranchAddress("recalg",&recalg);
113 n103-> SetBranchAddress("ndedxhit",&ndedxhit);
114 n103-> SetBranchAddress("ptrk",&ptrk);
115 n103-> SetBranchAddress("chidedx_e",&chidedx);
116 for(int j=0;j<n103->GetEntries();j++)
117 {
118 phlist.clear();
119 n103->GetEntry(j);
120 if(ptrk>pMax || ptrk<pMin) continue;
121 if(tes>1400) continue;
122 for(int k=0;k<ndedxhit;k++)
123 {
124 dEdx[k] = exsvc->StandardHitCorrec(0,runFlag,2,runNO,evtNO,pathlength[k],wid[k],layid[k],dEdx[k],dd_in[k],eangle[k],costheta);
125 phlist.push_back(dEdx[k]);
126 }
127 dedx = cal_dedx_bitrunc(truncate, phlist);
128 dedx = exsvc->StandardTrackCorrec(0, runFlag, 2, runNO, evtNO, dedx, costheta, tes);
129 int ipBin = (Int_t)floor((ptrk-pMin)/bin_step);
130 m_dedx[ipBin]->Fill(dedx);
131 if(charge>0) m_pos_dedx[ipBin]->Fill(dedx);
132 if(charge<0) m_neg_dedx[ipBin]->Fill(dedx);
133
134 usedhit = ndedxhit*truncate;
135 set_dEdx(1,dedx,recalg,runFlag,usedhit,ptrk,acos(costheta),1.5,vFlag,tes);
136 double chi = m_chi_ex[ParticleFlag];
137 //cout<<"chidedx= "<<chidedx<<" particleType= "<<ParticleFlag<<" new chi= "<<chi<<endl;
138 m_chi[ipBin]->Fill(chi);
139 if(charge>0) m_pos_chi[ipBin]->Fill(chi);
140 if(charge<0) m_neg_chi[ipBin]->Fill(chi);
141 m_hits[ipBin]->Fill(usedhit);
142
143 Vec_dedx.push_back(dedx);
144 Vec_ptrk.push_back(ptrk);
145 }
146 }
147}
148
149void DedxCalibMomentum::AnalyseHists()
150{
151 MsgStream log(msgSvc(), name());
152 log<<MSG::INFO<<"DedxCalibMomentum::AnalyseHists()"<<endreq;
153
154 gStyle->SetOptFit(1111);
155 for(int i=0; i<nbins; i++)
156 {
157 m_dedx[i]->Fit("gaus", "Q" );
158 m_pos_dedx[i]->Fit("gaus", "Q" );
159 m_neg_dedx[i]->Fit("gaus", "Q" );
160 m_chi[i]->Fit("gaus", "Q" );
161 m_pos_chi[i]->Fit("gaus", "Q" );
162 m_neg_chi[i]->Fit("gaus", "Q" );
163 }
164}
165
166void DedxCalibMomentum::WriteHists()
167{
168 MsgStream log(msgSvc(), name());
169 log<<MSG::INFO<<"DedxCalibMomentum::WriteHists()"<< endreq;
170
171 double chientryNo[nbins]={0},chimean[nbins]={0},chimeanerr[nbins]={0},chisigma[nbins]={0},chisq[nbins]={0};
172 double pos_chientryNo[nbins]={0},pos_chimean[nbins]={0},pos_chimeanerr[nbins]={0},pos_chisigma[nbins]={0},pos_chisq[nbins]={0};
173 double neg_chientryNo[nbins]={0},neg_chimean[nbins]={0},neg_chimeanerr[nbins]={0},neg_chisigma[nbins]={0},neg_chisq[nbins]={0};
174 double fitentryNo[nbins]={0},fitmean[nbins]={0},fitmeanerr[nbins]={0},fitsigma[nbins]={0},fitchisq[nbins]={0};
175 double pos_fitentryNo[nbins]={0},pos_fitmean[nbins]={0},pos_fitmeanerr[nbins]={0},pos_fitsigma[nbins]={0},pos_fitchisq[nbins]={0};
176 double neg_fitentryNo[nbins]={0},neg_fitmean[nbins]={0},neg_fitmeanerr[nbins]={0},neg_fitsigma[nbins]={0},neg_fitchisq[nbins]={0};
177 double hits_mean[nbins]={0},hits_sigma[nbins]={0};
178 double pBin[nbins]={0};
179
180 for(int i=0;i<nbins;i++)
181 {
182 pBin[i] = (i+0.5)*bin_step + pMin;
183
184 chientryNo[i] = m_chi[i]->Integral();
185 if(chientryNo[i]>100)
186 {
187 chimean[i] = m_chi[i]->GetFunction("gaus")->GetParameter(1);
188 chimeanerr[i] = m_chi[i]->GetFunction("gaus")->GetParError(1);
189 chisigma[i] = m_chi[i]->GetFunction("gaus")->GetParameter(2);
190 chisq[i] = (m_chi[i]->GetFunction("gaus")->GetChisquare())/(m_chi[i]->GetFunction("gaus")->GetNDF());
191 }
192 pos_chientryNo[i] = m_pos_chi[i]->Integral();
193 if(pos_chientryNo[i]>100)
194 {
195 pos_chimean[i] = m_pos_chi[i]->GetFunction("gaus")->GetParameter(1);
196 pos_chimeanerr[i] = m_pos_chi[i]->GetFunction("gaus")->GetParError(1);
197 pos_chisigma[i] = m_pos_chi[i]->GetFunction("gaus")->GetParameter(2);
198 pos_chisq[i] = (m_pos_chi[i]->GetFunction("gaus")->GetChisquare())/(m_pos_chi[i]->GetFunction("gaus")->GetNDF());
199 }
200 neg_chientryNo[i] = m_neg_chi[i]->Integral();
201 if(neg_chientryNo[i]>100)
202 {
203 neg_chimean[i] = m_neg_chi[i]->GetFunction("gaus")->GetParameter(1);
204 neg_chimeanerr[i] = m_neg_chi[i]->GetFunction("gaus")->GetParError(1);
205 neg_chisigma[i] = m_neg_chi[i]->GetFunction("gaus")->GetParameter(2);
206 neg_chisq[i] = (m_neg_chi[i]->GetFunction("gaus")->GetChisquare())/(m_neg_chi[i]->GetFunction("gaus")->GetNDF());
207 }
208
209 fitentryNo[i] = m_dedx[i]->Integral();
210 if(fitentryNo[i]>100)
211 {
212 fitmean[i] = m_dedx[i]->GetFunction("gaus")->GetParameter(1);
213 fitmeanerr[i] = m_dedx[i]->GetFunction("gaus")->GetParError(1);
214 fitsigma[i] = m_dedx[i]->GetFunction("gaus")->GetParameter(2);
215 fitchisq[i] = (m_dedx[i]->GetFunction("gaus")->GetChisquare())/(m_dedx[i]->GetFunction("gaus")->GetNDF());
216 }
217 pos_fitentryNo[i] = m_pos_dedx[i]->Integral();
218 if(pos_fitentryNo[i]>100)
219 {
220 pos_fitmean[i] = m_pos_dedx[i]->GetFunction("gaus")->GetParameter(1);
221 pos_fitmeanerr[i] = m_pos_dedx[i]->GetFunction("gaus")->GetParError(1);
222 pos_fitsigma[i] = m_pos_dedx[i]->GetFunction("gaus")->GetParameter(2);
223 pos_fitchisq[i] = (m_pos_dedx[i]->GetFunction("gaus")->GetChisquare())/(m_pos_dedx[i]->GetFunction("gaus")->GetNDF());
224 }
225 neg_fitentryNo[i] = m_neg_dedx[i]->Integral();
226 if(neg_fitentryNo[i]>100)
227 {
228 neg_fitmean[i] = m_neg_dedx[i]->GetFunction("gaus")->GetParameter(1);
229 neg_fitmeanerr[i] = m_neg_dedx[i]->GetFunction("gaus")->GetParError(1);
230 neg_fitsigma[i] = m_neg_dedx[i]->GetFunction("gaus")->GetParameter(2);
231 neg_fitchisq[i] = (m_neg_dedx[i]->GetFunction("gaus")->GetChisquare())/(m_neg_dedx[i]->GetFunction("gaus")->GetNDF());
232 }
233
234 hits_mean[i] = m_hits[i]->GetMean();
235 hits_sigma[i] = m_hits[i]->GetRMS();
236 }
237
238 double dedx1[Size] = {0};
239 double ptrk1[Size] = {0};
240 cout << "Vec_dedx.size() = " << Vec_dedx.size() << endl;
241 for(int i=0;i<Size && i< Vec_dedx.size();i++)
242 {
243 dedx1[i] = Vec_dedx[i];
244 ptrk1[i] = Vec_ptrk[i];
245 //cout<<"i= "<<i<<"dedx= "<<dedx1[i]<<" ptrk= "<<ptrk1[i]<<endl;
246 }
247
248 log<<MSG::INFO<<"begin generating root file!!! "<<endreq;
249 TFile* f = new TFile(m_rootfile.c_str(), "RECREATE");
250 for(int i=0;i<nbins;i++)
251 {
252 m_chi[i]->Write();
253 m_pos_chi[i]->Write();
254 m_neg_chi[i]->Write();
255 m_dedx[i]->Write();
256 m_pos_dedx[i]->Write();
257 m_neg_dedx[i]->Write();
258 m_hits[i]->Write();
259 }
260
261 TTree* momInfor = new TTree("momInfor","momInfor");
262 momInfor->Branch("chientryNo",chientryNo,"chientryNo[50]/D");
263 momInfor->Branch("chimean",chimean,"chimean[50]/D");
264 momInfor->Branch("chimeanerr",chimeanerr,"chimeanerr[50]/D");
265 momInfor->Branch("chisigma",chisigma,"chisigma[50]/D");
266 momInfor->Branch("chisq",chisq,"chisq[50]/D");
267 momInfor->Branch("pos_chientryNo",pos_chientryNo,"pos_chientryNo[50]/D");
268 momInfor->Branch("pos_chimean",pos_chimean,"pos_chimean[50]/D");
269 momInfor->Branch("pos_chimeanerr",pos_chimeanerr,"pos_chimeanerr[50]/D");
270 momInfor->Branch("pos_chisigma",pos_chisigma,"pos_chisigma[50]/D");
271 momInfor->Branch("pos_chisq",pos_chisq,"pos_chisq[50]/D");
272 momInfor->Branch("neg_chientryNo",neg_chientryNo,"neg_chientryNo[50]/D");
273 momInfor->Branch("neg_chimean",neg_chimean,"neg_chimean[50]/D");
274 momInfor->Branch("neg_chimeanerr",neg_chimeanerr,"neg_chimeanerr[50]/D");
275 momInfor->Branch("neg_chisigma",neg_chisigma,"neg_chisigma[50]/D");
276 momInfor->Branch("neg_chisq",neg_chisq,"neg_chisq[50]/D");
277 momInfor->Branch("fitentryNo",fitentryNo,"fitentryNo[50]/D");
278 momInfor->Branch("fitmean",fitmean,"fitmean[50]/D");
279 momInfor->Branch("fitmeanerr",fitmeanerr,"fitmeanerr[50]/D");
280 momInfor->Branch("fitsigma",fitsigma,"fitsigma[50]/D");
281 momInfor->Branch("fitchisq",fitchisq,"fitchisq[50]/D");
282 momInfor->Branch("pos_fitentryNo",pos_fitentryNo,"pos_fitentryNo[50]/D");
283 momInfor->Branch("pos_fitmean",pos_fitmean,"pos_fitmean[50]/D");
284 momInfor->Branch("pos_fitmeanerr",pos_fitmeanerr,"pos_fitmeanerr[50]/D");
285 momInfor->Branch("pos_fitsigma",pos_fitsigma,"pos_fitsigma[50]/D");
286 momInfor->Branch("pos_fitchisq",pos_fitchisq,"pos_fitchisq[50]/D");
287 momInfor->Branch("neg_fitentryNo",neg_fitentryNo,"neg_fitentryNo[50]/D");
288 momInfor->Branch("neg_fitmean",neg_fitmean,"neg_fitmean[50]/D");
289 momInfor->Branch("neg_fitmeanerr",neg_fitmeanerr,"neg_fitmeanerr[50]/D");
290 momInfor->Branch("neg_fitsigma",neg_fitsigma,"neg_fitsigma[50]/D");
291 momInfor->Branch("neg_fitchisq",neg_fitchisq,"neg_fitchisq[50]/D");
292 momInfor->Branch("hits_mean",hits_mean,"hits_mean[50]/D");
293 momInfor->Branch("hits_sigma",hits_sigma,"hits_sigma[50]/D");
294 momInfor->Branch("pBin",pBin,"pBin[50]/D");
295 momInfor-> Branch("ptrk1",ptrk1,"ptrk1[700000]/D");
296 momInfor-> Branch("dedx1",dedx1,"dedx1[700000]/D");
297 momInfor->Fill();
298 momInfor->Write();
299
300 TCanvas c1("c1", "canvas", 500, 400);
301 c1.Print((m_rootfile+".ps[").c_str());
302 momInfor -> Draw("dedx1:ptrk1","dedx1>200 && dedx1<1200");
303 c1.Print((m_rootfile+".ps").c_str());
304 for(Int_t i=0;i<nbins;i++)
305 {
306 m_chi[i]->Draw();
307 c1.Print((m_rootfile+".ps").c_str());
308 }
309 for(Int_t i=0;i<nbins;i++)
310 {
311 m_pos_chi[i]->Draw();
312 c1.Print((m_rootfile+".ps").c_str());
313 }
314 for(Int_t i=0;i<nbins;i++)
315 {
316 m_neg_chi[i]->Draw();
317 c1.Print((m_rootfile+".ps").c_str());
318 }
319 for(Int_t i=0;i<nbins;i++)
320 {
321 m_dedx[i]->Draw();
322 c1.Print((m_rootfile+".ps").c_str());
323 }
324 for(Int_t i=0;i<nbins;i++)
325 {
326 m_pos_dedx[i]->Draw();
327 c1.Print((m_rootfile+".ps").c_str());
328 }
329 for(Int_t i=0;i<nbins;i++)
330 {
331 m_neg_dedx[i]->Draw();
332 c1.Print((m_rootfile+".ps").c_str());
333 }
334 c1.Print((m_rootfile+".ps]").c_str());
335 f->Close();
336
337 for(int i=0;i<nbins;i++)
338 {
339 delete m_chi[i];
340 delete m_pos_chi[i];
341 delete m_neg_chi[i];
342 delete m_dedx[i];
343 delete m_pos_dedx[i];
344 delete m_neg_dedx[i];
345 delete m_hits[i];
346 }
347}
nbins
const int nbins
Definition: CheckRung/check/dEdxVsCos.cxx:18
Branch
curve Branch("CurveSize",&CurveSize,"CurveSize/I")
Draw
legend Draw()
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
SetBranchAddress
data SetBranchAddress("time",&time)
Size
#define Size
Definition: DedxCalibCostheta.cxx:15
chihist_min
const double chihist_min
Definition: DedxCalibMomentum.cxx:20
nbins
const int nbins
Definition: DedxCalibMomentum.cxx:19
dedxhist_max
const double dedxhist_max
Definition: DedxCalibMomentum.cxx:23
hits_max
const double hits_max
Definition: DedxCalibMomentum.cxx:25
chihist_max
const double chihist_max
Definition: DedxCalibMomentum.cxx:21
dedxhist_min
const double dedxhist_min
Definition: DedxCalibMomentum.cxx:22
hits_min
const double hits_min
Definition: DedxCalibMomentum.cxx:24
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
DedxCalibMomentum::DedxCalibMomentum
DedxCalibMomentum(const std::string &name, ISvcLocator *pSvcLocator)
Definition: DedxCalibMomentum.cxx:27
DedxCalib::m_chi_ex
double m_chi_ex[5]
Definition: DedxCalib.h:46
DedxCalib::set_dEdx
void set_dEdx(int landau, float dEdx, int trkalg, int runflag, int dedxhit_use, float ptrk, float theta, float pl_rp, int vflag[3], double t0)
Definition: DedxCalib.cxx:184
DedxCalib::ParticleFlag
int ParticleFlag
Definition: DedxCalib.h:50
DedxCalib::ReadRecFileList
void ReadRecFileList()
Definition: DedxCalib.cxx:89
DedxCalib::m_rootfile
std::string m_rootfile
Definition: DedxCalib.h:55
DedxCalib::vFlag
int vFlag[3]
Definition: DedxCalib.h:41
DedxCalib::exsvc
IDedxCorrecSvc * exsvc
Definition: DedxCalib.h:30
DedxCalib::cal_dedx_bitrunc
double cal_dedx_bitrunc(float truncate, std::vector< double > phlist)
Definition: DedxCalib.cxx:108
DedxCalib::m_recFileVector
vector< string > m_recFileVector
Definition: DedxCalib.h:48
DedxCalib::truncate
float truncate
Definition: DedxCalib.h:33
IDedxCorrecSvc::StandardHitCorrec
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
IDedxCorrecSvc::StandardTrackCorrec
virtual double StandardTrackCorrec(int calib_rec_Flag, int typFlag, int ntpFlag, int runNO, int evtNO, double ex, double costheta, double t0) const =0
costheta
float costheta
Definition: prepare_data.cxx:28
charge
float charge
Definition: prepare_data.cxx:28
ptrk
float ptrk
Definition: prepare_data.cxx:28