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