BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
DedxCalibEAng.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 <sstream>
10#include <string>
11
12#include "DedxCalibAlg/DedxCalibEAng.h"
13DECLARE_COMPONENT( DedxCalibEAng )
14
15using namespace std;
16
17const double StepEAng = ( PhiMax - PhiMin ) / NumSlices;
18
19DedxCalibEAng::DedxCalibEAng( const std::string& name, ISvcLocator* pSvcLocator )
20 : DedxCalib( name, pSvcLocator ) {
21 declareProperty( "Truncate", m_truncate = 0.7 );
22 declareProperty( "PMin", pMin = 0.2 );
23 declareProperty( "PMax", pMax = 2.3 );
24}
25
26void DedxCalibEAng::BookHists() {
27 MsgStream log( msgSvc(), name() );
28 log << MSG::INFO << "DedxCalibEAng::BookHists()" << endmsg;
29
30 ReadRecFileList();
31
32 m_eangle = new TH1F*[NumSlices];
33 m_trunc_eangle = new TH1F*[NumSlices];
34 m_pos_eangle = new TH1F*[NumSlices];
35 m_neg_eangle = new TH1F*[NumSlices];
36 stringstream hlabel;
37 for ( int i = 0; i < NumSlices; i++ )
38 {
39 hlabel.str( "" );
40 hlabel << "dEdx"
41 << "_eang" << i;
42 m_eangle[i] = new TH1F( hlabel.str().c_str(), hlabel.str().c_str(), NumHistBins,
43 MinHistValue, MaxHistValue );
44 hlabel.str( "" );
45 hlabel << "trunc_dEdx"
46 << "_eang" << i;
47 m_trunc_eangle[i] = new TH1F( hlabel.str().c_str(), hlabel.str().c_str(), NumHistBins,
48 MinHistValue, MaxHistValue );
49 hlabel.str( "" );
50 hlabel << "pos_dEdx"
51 << "_eang" << i;
52 m_pos_eangle[i] = new TH1F( hlabel.str().c_str(), hlabel.str().c_str(), NumHistBins,
53 MinHistValue, MaxHistValue );
54 hlabel.str( "" );
55 hlabel << "neg_dEdx"
56 << "_eang" << i;
57 m_neg_eangle[i] = new TH1F( hlabel.str().c_str(), hlabel.str().c_str(), NumHistBins,
58 MinHistValue, MaxHistValue );
59 }
60 hlabel.str( "" );
61 hlabel << "dEdxVsEAng";
62 m_EAngAverdE = new TH1F( hlabel.str().c_str(), "dEdx & EAng", NumSlices, PhiMin, PhiMax );
63 hlabel.str( "" );
64 hlabel << "trunc_dEdxVsEAng";
65 m_trunc_EAngAverdE =
66 new TH1F( hlabel.str().c_str(), "trunc_dEdx & EAng", NumSlices, PhiMin, PhiMax );
67 hlabel.str( "" );
68 hlabel << "pos_dEdxVsEAng";
69 m_pos_EAngAverdE =
70 new TH1F( hlabel.str().c_str(), "pos_dEdx & EAng", NumSlices, PhiMin, PhiMax );
71 hlabel.str( "" );
72 hlabel << "neg_dEdxVsEAng";
73 m_neg_EAngAverdE =
74 new TH1F( hlabel.str().c_str(), "neg_dEdx & EAng", NumSlices, PhiMin, PhiMax );
75}
76
77void DedxCalibEAng::FillHists() {
78 MsgStream log( msgSvc(), name() );
79 log << MSG::INFO << "DedxCalibEAng::FillHists()" << endmsg;
80
81 TFile* f;
82 TTree* n102;
83 string runlist;
84
85 double dedx = 0;
86 float runNO = 0, evtNO = 0, runFlag = 0, pathlength = 0, wid = 0, layid = 0, dd_in = 0,
87 driftdist = 0, eangle = 0, zhit = 0, costheta = 0, tes = 0, ptrk = 0, charge = 0;
88 cout << "m_recFileVector.size()= " << m_recFileVector.size() << endl;
89
90 vector<double> m_phlist( 0 );
91 vector<double> m_phlist_test( 0 );
92 vector<int> m_iEAng_list( 0 );
93 double old_costheta( 0 ), old_ptrk( 0 );
94
95 for ( unsigned int i = 0; i < m_recFileVector.size(); i++ )
96 {
97 runlist = m_recFileVector[i];
98 cout << "runlist: " << runlist.c_str() << endl;
99 f = new TFile( runlist.c_str() );
100 n102 = (TTree*)f->Get( "n102" );
101 n102->SetBranchAddress( "adc_raw", &dedx );
102 n102->SetBranchAddress( "path_rphi", &pathlength );
103 n102->SetBranchAddress( "wire", &wid );
104 n102->SetBranchAddress( "layer", &layid );
105 n102->SetBranchAddress( "doca_in", &dd_in );
106 n102->SetBranchAddress( "driftdist", &driftdist );
107 n102->SetBranchAddress( "eangle", &eangle );
108 n102->SetBranchAddress( "zhit", &zhit );
109 n102->SetBranchAddress( "runNO", &runNO );
110 n102->SetBranchAddress( "evtNO", &evtNO );
111 n102->SetBranchAddress( "runFlag", &runFlag );
112 n102->SetBranchAddress( "costheta1", &costheta );
113 n102->SetBranchAddress( "t01", &tes );
114 n102->SetBranchAddress( "ptrk1", &ptrk );
115 n102->SetBranchAddress( "charge", &charge );
116 for ( int j = 0; j < n102->GetEntries(); j++ )
117 {
118 n102->GetEntry( j );
119 // begin of save truncated hits information for EACH track
120 if ( ( fabs( costheta - old_costheta ) > 0.001 || fabs( ptrk - old_ptrk > 0.001 ) ) &&
121 m_phlist.size() >= 2 )
122 {
123 sort( m_phlist_test.begin(), m_phlist_test.end() );
124 unsigned int nsample = (unsigned int)( m_phlist_test.size() * m_truncate );
125 double m_trunc_value = m_phlist_test[nsample];
126 for ( unsigned int i = 0; i < m_phlist.size(); i++ )
127 {
128 if ( m_phlist[i] > m_trunc_value ) continue;
129 m_trunc_eangle[m_iEAng_list[i]]->Fill( m_phlist[i] );
130 }
131 m_phlist.clear();
132 m_phlist_test.clear();
133 m_iEAng_list.clear();
134 }
135 old_costheta = costheta;
136 old_ptrk = ptrk;
137 // end of save truncated information
138 if ( eangle > 0.25 ) eangle = eangle - 0.5;
139 else if ( eangle < -0.25 ) eangle = eangle + 0.5;
140 if ( eangle > PhiMin && eangle < PhiMax )
141 {
142 if ( ptrk > pMax || ptrk < pMin ) continue;
143 if ( tes > 1400 ) continue;
144 // if (ptrk<0.1) continue;
145 if ( layid < 4 ) continue;
146 else if ( layid < 8 )
147 {
148 if ( dedx < MinADCValuecut || dedx > MaxADCValuecut ||
149 pathlength > RPhi_PathMaxCut || pathlength < Iner_RPhi_PathMinCut ||
150 driftdist > Iner_DriftDistCut )
151 continue;
152 }
153 else
154 {
155 if ( dedx < MinADCValuecut || dedx > MaxADCValuecut ||
156 pathlength > RPhi_PathMaxCut || pathlength < Out_RPhi_PathMinCut ||
157 driftdist > Out_DriftDistCut )
158 continue;
159 }
160 int iEAng = 0;
161 iEAng = (int)( ( eangle - PhiMin ) / StepEAng );
162 dedx = exsvc->StandardHitCorrec( 0, runFlag, 2, runNO, evtNO, pathlength, wid, layid,
163 dedx, dd_in, eangle, costheta );
164 dedx = exsvc->StandardTrackCorrec( 0, runFlag, 2, runNO, evtNO, dedx, costheta, tes );
165 m_eangle[iEAng]->Fill( dedx );
166 m_phlist.push_back( dedx );
167 m_phlist_test.push_back( dedx );
168 m_iEAng_list.push_back( iEAng );
169 if ( charge > 0 ) m_pos_eangle[iEAng]->Fill( dedx );
170 if ( charge < 0 ) m_neg_eangle[iEAng]->Fill( dedx );
171 }
172 }
173 }
174}
175
176void DedxCalibEAng::AnalyseHists() // fit histograms
177{
178 MsgStream log( msgSvc(), name() );
179 log << MSG::INFO << "DedxCalibEAng::AnalyseHists()" << endmsg;
180
181 TF1* func;
182 Double_t entry = 0, mean = 0, rms = 0;
183 Double_t binmax = 0;
184 Double_t lp[3] = { 0 };
185 gStyle->SetOptFit( 1111 );
186 for ( Int_t ip = 0; ip < NumSlices; ip++ )
187 {
188 entry = m_eangle[ip]->Integral();
189 if ( entry < 10 ) continue;
190 mean = m_eangle[ip]->GetMean();
191 rms = m_eangle[ip]->GetRMS();
192 binmax = m_eangle[ip]->GetMaximumBin();
193 lp[1] = m_eangle[ip]->GetBinCenter( binmax );
194 lp[2] = 200;
195 lp[0] =
196 ( m_eangle[ip]->GetBinContent( binmax ) + m_eangle[ip]->GetBinContent( binmax - 1 ) +
197 m_eangle[ip]->GetBinContent( binmax + 1 ) ) /
198 3;
199
200 if ( m_phShapeFlag == 0 )
201 {
202 func = new TF1( "mylan", mylan, 10, 3000, 4 );
203 func->SetParLimits( 0, 0, entry );
204 func->SetParLimits( 1, 0, mean );
205 func->SetParLimits( 2, 10, rms );
206 func->SetParLimits( 3, -1, 0 );
207 func->SetParameters( entry / 36., mean - 200, rms / 2.3, -0.5 );
208 }
209 else if ( m_phShapeFlag == 1 )
210 {
211 func = new TF1( "landaun", landaun, 10, 3000, 3 );
212 func->SetParameters( lp[0], lp[1], lp[2] );
213 }
214 else if ( m_phShapeFlag == 2 )
215 {
216 func = new TF1( "Landau", Landau, 10, 3000, 2 );
217 func->SetParameters( 0.7 * mean, rms );
218 }
219 else if ( m_phShapeFlag == 3 )
220 {
221 func = new TF1( "Vavilov", Vavilov, 10, 3000, 2 );
222 func->SetParameters( 0.05, 0.6 );
223 }
224 else if ( m_phShapeFlag == 4 )
225 {
226 func = new TF1( "AsymGauss", AsymGauss, 10, 3000, 4 );
227 func->SetParameters( entry, mean, rms, 1.0 );
228 }
229 func->SetLineWidth( 2.1 );
230 func->SetLineColor( 2 );
231
232 m_eangle[ip]->Fit( func, "rq" ); // , "", mean-rms, mean+rms/2);
233 m_pos_eangle[ip]->Fit( func, "rq" ); // , "", mean-rms, mean+rms/2);
234 m_neg_eangle[ip]->Fit( func, "rq" ); // , "", mean-rms, mean+rms/2);
235 }
236}
237
238void DedxCalibEAng::WriteHists() {
239 MsgStream log( msgSvc(), name() );
240 log << MSG::INFO << "DedxCalibEAng::WriteHists()" << endmsg;
241
242 double entryNo[NumSlices] = { 0 }, trunc_mean[NumSlices] = { 0 },
243 trunc_sigma[NumSlices] = { 0 }, fitmean[NumSlices] = { 0 },
244 fitmeanerr[NumSlices] = { 0 }, fitsigma[NumSlices] = { 0 }, gain[NumSlices] = { 0 },
245 chisq[NumSlices] = { 0 };
246 double trunc_gain[NumSlices] = { 0 };
247 double pos_entryNo[NumSlices] = { 0 }, pos_mean[NumSlices] = { 0 },
248 pos_sigma[NumSlices] = { 0 }, pos_fitmean[NumSlices] = { 0 },
249 pos_fitmeanerr[NumSlices] = { 0 }, pos_fitsigma[NumSlices] = { 0 },
250 pos_gain[NumSlices] = { 0 }, pos_chisq[NumSlices] = { 0 };
251 double neg_entryNo[NumSlices] = { 0 }, neg_mean[NumSlices] = { 0 },
252 neg_sigma[NumSlices] = { 0 }, neg_fitmean[NumSlices] = { 0 },
253 neg_fitmeanerr[NumSlices] = { 0 }, neg_fitsigma[NumSlices] = { 0 },
254 neg_gain[NumSlices] = { 0 }, neg_chisq[NumSlices] = { 0 };
255 double Ip_eangle[NumSlices] = { 0 }, eangle[NumSlices] = { 0 };
256
257 double Norm( 0 ), trunc_Norm( 0 );
258 int count( 0 );
259 for ( Int_t ip = 0; ip < NumSlices; ip++ )
260 {
261 Ip_eangle[ip] = ip;
262 eangle[ip] = ( ip + 0.5 ) * StepEAng + PhiMin;
263
264 entryNo[ip] = m_eangle[ip]->Integral();
265 trunc_mean[ip] = m_trunc_eangle[ip]->GetMean();
266 trunc_sigma[ip] = m_trunc_eangle[ip]->GetRMS();
267 pos_entryNo[ip] = m_pos_eangle[ip]->Integral();
268 pos_mean[ip] = m_pos_eangle[ip]->GetMean();
269 pos_sigma[ip] = m_pos_eangle[ip]->GetRMS();
270 neg_entryNo[ip] = m_neg_eangle[ip]->Integral();
271 neg_mean[ip] = m_neg_eangle[ip]->GetMean();
272 neg_sigma[ip] = m_neg_eangle[ip]->GetRMS();
273 if ( entryNo[ip] < 10 || pos_entryNo[ip] < 10 || neg_entryNo[ip] < 10 ) continue;
274
275 if ( m_phShapeFlag == 0 )
276 {
277 fitmean[ip] = m_eangle[ip]->GetFunction( "mylan" )->GetParameter( 1 );
278 fitmeanerr[ip] = m_eangle[ip]->GetFunction( "mylan" )->GetParError( 1 );
279 fitsigma[ip] = m_eangle[ip]->GetFunction( "mylan" )->GetParameter( 2 );
280 chisq[ip] = ( m_eangle[ip]->GetFunction( "mylan" )->GetChisquare() ) /
281 ( m_eangle[ip]->GetFunction( "mylan" )->GetNDF() );
282
283 pos_fitmean[ip] = m_pos_eangle[ip]->GetFunction( "mylan" )->GetParameter( 1 );
284 pos_fitmeanerr[ip] = m_pos_eangle[ip]->GetFunction( "mylan" )->GetParError( 1 );
285 pos_fitsigma[ip] = m_pos_eangle[ip]->GetFunction( "mylan" )->GetParameter( 2 );
286 pos_chisq[ip] = ( m_pos_eangle[ip]->GetFunction( "mylan" )->GetChisquare() ) /
287 ( m_pos_eangle[ip]->GetFunction( "mylan" )->GetNDF() );
288
289 neg_fitmean[ip] = m_neg_eangle[ip]->GetFunction( "mylan" )->GetParameter( 1 );
290 neg_fitmeanerr[ip] = m_neg_eangle[ip]->GetFunction( "mylan" )->GetParError( 1 );
291 neg_fitsigma[ip] = m_neg_eangle[ip]->GetFunction( "mylan" )->GetParameter( 2 );
292 neg_chisq[ip] = ( m_neg_eangle[ip]->GetFunction( "mylan" )->GetChisquare() ) /
293 ( m_neg_eangle[ip]->GetFunction( "mylan" )->GetNDF() );
294 }
295 if ( m_phShapeFlag == 1 )
296 {
297 fitmean[ip] = m_eangle[ip]->GetFunction( "landaun" )->GetParameter( 1 );
298 fitmeanerr[ip] = m_eangle[ip]->GetFunction( "landaun" )->GetParError( 1 );
299 fitsigma[ip] = m_eangle[ip]->GetFunction( "landaun" )->GetParameter( 2 );
300 chisq[ip] = ( m_eangle[ip]->GetFunction( "landaun" )->GetChisquare() ) /
301 ( m_eangle[ip]->GetFunction( "mylan" )->GetNDF() );
302
303 pos_fitmean[ip] = m_pos_eangle[ip]->GetFunction( "landaun" )->GetParameter( 1 );
304 pos_fitmeanerr[ip] = m_pos_eangle[ip]->GetFunction( "landaun" )->GetParError( 1 );
305 pos_fitsigma[ip] = m_pos_eangle[ip]->GetFunction( "landaun" )->GetParameter( 2 );
306 pos_chisq[ip] = ( m_pos_eangle[ip]->GetFunction( "landaun" )->GetChisquare() ) /
307 ( m_pos_eangle[ip]->GetFunction( "landaun" )->GetNDF() );
308
309 neg_fitmean[ip] = m_neg_eangle[ip]->GetFunction( "landaun" )->GetParameter( 1 );
310 neg_fitmeanerr[ip] = m_neg_eangle[ip]->GetFunction( "landaun" )->GetParError( 1 );
311 neg_fitsigma[ip] = m_neg_eangle[ip]->GetFunction( "landaun" )->GetParameter( 2 );
312 neg_chisq[ip] = ( m_neg_eangle[ip]->GetFunction( "landaun" )->GetChisquare() ) /
313 ( m_neg_eangle[ip]->GetFunction( "landaun" )->GetNDF() );
314 }
315 if ( m_phShapeFlag == 2 )
316 {
317 fitmean[ip] = m_eangle[ip]->GetFunction( "Landau" )->GetParameter( 1 );
318 fitmeanerr[ip] = m_eangle[ip]->GetFunction( "Landau" )->GetParError( 1 );
319 fitsigma[ip] = m_eangle[ip]->GetFunction( "Landau" )->GetParameter( 2 );
320 chisq[ip] = ( m_eangle[ip]->GetFunction( "Landau" )->GetChisquare() ) /
321 ( m_eangle[ip]->GetFunction( "Landau" )->GetNDF() );
322
323 pos_fitmean[ip] = m_pos_eangle[ip]->GetFunction( "Landau" )->GetParameter( 1 );
324 pos_fitmeanerr[ip] = m_pos_eangle[ip]->GetFunction( "Landau" )->GetParError( 1 );
325 pos_fitsigma[ip] = m_pos_eangle[ip]->GetFunction( "Landau" )->GetParameter( 2 );
326 pos_chisq[ip] = ( m_pos_eangle[ip]->GetFunction( "Landau" )->GetChisquare() ) /
327 ( m_pos_eangle[ip]->GetFunction( "Landau" )->GetNDF() );
328
329 neg_fitmean[ip] = m_neg_eangle[ip]->GetFunction( "Landau" )->GetParameter( 1 );
330 neg_fitmeanerr[ip] = m_neg_eangle[ip]->GetFunction( "Landau" )->GetParError( 1 );
331 neg_fitsigma[ip] = m_neg_eangle[ip]->GetFunction( "Landau" )->GetParameter( 2 );
332 neg_chisq[ip] = ( m_neg_eangle[ip]->GetFunction( "Landau" )->GetChisquare() ) /
333 ( m_neg_eangle[ip]->GetFunction( "Landau" )->GetNDF() );
334 }
335 if ( m_phShapeFlag == 3 )
336 {
337 fitmean[ip] = m_eangle[ip]->GetFunction( "Vavilov" )->GetParameter( 1 );
338 fitmeanerr[ip] = m_eangle[ip]->GetFunction( "Vavilov" )->GetParError( 1 );
339 fitsigma[ip] = m_eangle[ip]->GetFunction( "Vavilov" )->GetParameter( 2 );
340 chisq[ip] = ( m_eangle[ip]->GetFunction( "Vavilov" )->GetChisquare() ) /
341 ( m_eangle[ip]->GetFunction( "Vavilov" )->GetNDF() );
342
343 pos_fitmean[ip] = m_pos_eangle[ip]->GetFunction( "Vavilov" )->GetParameter( 1 );
344 pos_fitmeanerr[ip] = m_pos_eangle[ip]->GetFunction( "Vavilov" )->GetParError( 1 );
345 pos_fitsigma[ip] = m_pos_eangle[ip]->GetFunction( "Vavilov" )->GetParameter( 2 );
346 pos_chisq[ip] = ( m_pos_eangle[ip]->GetFunction( "Vavilov" )->GetChisquare() ) /
347 ( m_pos_eangle[ip]->GetFunction( "Vavilov" )->GetNDF() );
348
349 neg_fitmean[ip] = m_neg_eangle[ip]->GetFunction( "Vavilov" )->GetParameter( 1 );
350 neg_fitmeanerr[ip] = m_neg_eangle[ip]->GetFunction( "Vavilov" )->GetParError( 1 );
351 neg_fitsigma[ip] = m_neg_eangle[ip]->GetFunction( "Vavilov" )->GetParameter( 2 );
352 neg_chisq[ip] = ( m_neg_eangle[ip]->GetFunction( "Vavilov" )->GetChisquare() ) /
353 ( m_neg_eangle[ip]->GetFunction( "Vavilov" )->GetNDF() );
354 }
355 if ( m_phShapeFlag == 4 )
356 {
357 fitmean[ip] = m_eangle[ip]->GetFunction( "AsymGauss" )->GetParameter( 1 );
358 fitmeanerr[ip] = m_eangle[ip]->GetFunction( "AsymGauss" )->GetParError( 1 );
359 fitsigma[ip] = m_eangle[ip]->GetFunction( "AsymGauss" )->GetParameter( 2 );
360 chisq[ip] = ( m_eangle[ip]->GetFunction( "AsymGauss" )->GetChisquare() ) /
361 ( m_eangle[ip]->GetFunction( "AsymGauss" )->GetNDF() );
362
363 pos_fitmean[ip] = m_pos_eangle[ip]->GetFunction( "AsymGauss" )->GetParameter( 1 );
364 pos_fitmeanerr[ip] = m_pos_eangle[ip]->GetFunction( "AsymGauss" )->GetParError( 1 );
365 pos_fitsigma[ip] = m_pos_eangle[ip]->GetFunction( "AsymGauss" )->GetParameter( 2 );
366 pos_chisq[ip] = ( m_pos_eangle[ip]->GetFunction( "AsymGauss" )->GetChisquare() ) /
367 ( m_pos_eangle[ip]->GetFunction( "AsymGauss" )->GetNDF() );
368
369 neg_fitmean[ip] = m_neg_eangle[ip]->GetFunction( "AsymGauss" )->GetParameter( 1 );
370 neg_fitmeanerr[ip] = m_neg_eangle[ip]->GetFunction( "AsymGauss" )->GetParError( 1 );
371 neg_fitsigma[ip] = m_neg_eangle[ip]->GetFunction( "AsymGauss" )->GetParameter( 2 );
372 neg_chisq[ip] = ( m_neg_eangle[ip]->GetFunction( "AsymGauss" )->GetChisquare() ) /
373 ( m_neg_eangle[ip]->GetFunction( "AsymGauss" )->GetNDF() );
374 }
375 if ( entryNo[ip] > 1500 && fitmean[ip] > 0 && fitmeanerr[ip] > 0 && fitmeanerr[ip] < 10 &&
376 fitsigma[ip] > 0 && fitsigma[ip] < 200 )
377 {
378 Norm += fitmean[ip];
379 trunc_Norm += trunc_mean[ip];
380 count++;
381 }
382 }
383 Norm = Norm / count;
384 trunc_Norm = trunc_Norm / count;
385 // Norm = 550; // use an absolute normalization
386 cout << "count= " << count << " average value: " << Norm << " truncted: " << trunc_Norm
387 << endl;
388 for ( Int_t i = 0; i < NumSlices; i++ )
389 { // the difference between gain and mean is the normalization
390 if ( entryNo[i] > 10 && fitmean[i] > 0 && fitsigma[i] > 0 ) gain[i] = fitmean[i] / Norm;
391 else gain[i] = 1;
392 if ( entryNo[i] > 10 ) trunc_gain[i] = trunc_mean[i] / trunc_Norm;
393 else trunc_mean[i] = 1;
394 if ( pos_entryNo[i] > 10 && pos_fitmean[i] > 0 && pos_fitsigma[i] > 0 )
395 pos_gain[i] = pos_fitmean[i] / Norm;
396 else pos_gain[i] = 1;
397 if ( neg_entryNo[i] > 10 && neg_fitmean[i] > 0 && neg_fitsigma[i] > 0 )
398 neg_gain[i] = neg_fitmean[i] / Norm;
399 else neg_gain[i] = 1;
400 if ( gain[i] != 1 )
401 {
402 m_EAngAverdE->SetBinContent( i + 1, gain[i] );
403 m_EAngAverdE->SetBinError( i + 1, fitmeanerr[i] / Norm );
404 }
405 if ( trunc_gain[i] != 1 )
406 {
407 m_trunc_EAngAverdE->SetBinContent( i + 1, trunc_gain[i] );
408 m_trunc_EAngAverdE->SetBinError( i + 1, trunc_sigma[i] / trunc_Norm );
409 }
410 if ( pos_gain[i] != 1 )
411 {
412 m_pos_EAngAverdE->SetBinContent( i + 1, pos_gain[i] );
413 m_pos_EAngAverdE->SetBinError( i + 1, pos_fitmeanerr[i] / Norm );
414 }
415 if ( neg_gain[i] != 1 )
416 {
417 m_neg_EAngAverdE->SetBinContent( i + 1, neg_gain[i] );
418 m_neg_EAngAverdE->SetBinError( i + 1, neg_fitmeanerr[i] / Norm );
419 }
420 }
421
422 log << MSG::INFO << "begin getting calibration constant!!! " << endmsg;
423 double denangle[NumSlices] = { 0 };
424 int denangle_entry[1] = { 100 };
425 m_EAngAverdE->Fit( "pol1", "r", "", -0.25, 0.25 );
426 double par0 = m_EAngAverdE->GetFunction( "pol1" )->GetParameter( 0 );
427 double par1 = m_EAngAverdE->GetFunction( "pol1" )->GetParameter( 1 );
428 // double par2 = m_EAngAverdE->GetFunction("pol2")->GetParameter(2);
429 // par1 = -par1; // a trick is used here, because we don't really understand
430 // for(int i=0;i<NumSlices;i++) denangle[i] =
431 // (par0+par1*eangle[i])/par0;//+par2*pow(eangle[i],2))/par0;
432 for ( int i = 0; i < NumSlices; i++ ) denangle[i] = gain[i]; // use raw values instead of fit
433
434 log << MSG::INFO << "begin generating root file!!! " << endmsg;
435 TFile* f = new TFile( m_rootfile.c_str(), "recreate" );
436 for ( Int_t ip = 0; ip < NumSlices; ip++ )
437 {
438 m_eangle[ip]->Write();
439 m_trunc_eangle[ip]->Write();
440 m_pos_eangle[ip]->Write();
441 m_neg_eangle[ip]->Write();
442 }
443 m_EAngAverdE->Write();
444 m_trunc_EAngAverdE->Write();
445 m_pos_EAngAverdE->Write();
446 m_neg_EAngAverdE->Write();
447
448 TTree* entracalib = new TTree( "entracalib", "entracalib" );
449 entracalib->Branch( "1denangle", denangle, "denangle[100]/D" );
450 entracalib->Branch( "1denangle_entry", denangle_entry, "denangle_entry[1]/I" );
451 entracalib->Fill();
452 entracalib->Write();
453
454 TTree* ddgcalib = new TTree( "ddgcalib", "ddgcalib" );
455 ddgcalib->Branch( "hits", entryNo, "entryNo[100]/D" );
456 ddgcalib->Branch( "truncmean", trunc_mean, "trunc_mean[100]/D" );
457 ddgcalib->Branch( "truncsigma", trunc_sigma, "trunc_sigma[100]/D" );
458 ddgcalib->Branch( "fitmean", fitmean, "fitmean[100]/D" );
459 ddgcalib->Branch( "fitmeanerr", fitmeanerr, "fitmeanerr[100]/D" );
460 ddgcalib->Branch( "chi", fitsigma, "fitsigma[100]/D" );
461 ddgcalib->Branch( "gain", gain, "gain[100]/D" );
462 ddgcalib->Branch( "truncgain", trunc_gain, "trunc_gain[100]/D" );
463 ddgcalib->Branch( "chisq", chisq, "chisq[100]/D" );
464 ddgcalib->Branch( "pos_hits", pos_entryNo, "pos_entryNo[100]/D" );
465 ddgcalib->Branch( "pos_mean", pos_mean, "pos_mean[100]/D" );
466 ddgcalib->Branch( "pos_sigma", pos_sigma, "pos_sigma[100]/D" );
467 ddgcalib->Branch( "pos_fitmean", pos_fitmean, "pos_fitmean[100]/D" );
468 ddgcalib->Branch( "pos_fitmeanerr", pos_fitmeanerr, "pos_fitmeanerr[100]/D" );
469 ddgcalib->Branch( "pos_chi", pos_fitsigma, "pos_fitsigma[100]/D" );
470 ddgcalib->Branch( "pos_gain", pos_gain, "pos_gain[100]/D" );
471 ddgcalib->Branch( "pos_chisq", pos_chisq, "pos_chisq[100]/D" );
472 ddgcalib->Branch( "neg_hits", neg_entryNo, "neg_entryNo[100]/D" );
473 ddgcalib->Branch( "neg_mean", neg_mean, "neg_mean[100]/D" );
474 ddgcalib->Branch( "neg_sigma", neg_sigma, "neg_sigma[100]/D" );
475 ddgcalib->Branch( "neg_fitmean", neg_fitmean, "neg_fitmean[100]/D" );
476 ddgcalib->Branch( "neg_fitmeanerr", neg_fitmeanerr, "neg_fitmeanerr[100]/D" );
477 ddgcalib->Branch( "neg_chi", neg_fitsigma, "neg_fitsigma[100]/D" );
478 ddgcalib->Branch( "neg_gain", neg_gain, "neg_gain[100]/D" );
479 ddgcalib->Branch( "neg_chisq", neg_chisq, "neg_chisq[100]/D" );
480 ddgcalib->Branch( "Ip_eangle", Ip_eangle, "Ip_eangle[100]/D" );
481 ddgcalib->Branch( "eangle", eangle, "eangle[100]/D" );
482 ddgcalib->Fill();
483 ddgcalib->Write();
484 f->Close();
485
486 TCanvas c1( "c1", "canvas", 500, 400 );
487 c1.Print( ( m_rootfile + ".ps[" ).c_str() );
488 m_EAngAverdE->Draw();
489 c1.Print( ( m_rootfile + ".ps" ).c_str() );
490 m_trunc_EAngAverdE->Draw();
491 c1.Print( ( m_rootfile + ".ps" ).c_str() );
492 m_pos_EAngAverdE->Draw();
493 c1.Print( ( m_rootfile + ".ps" ).c_str() );
494 m_neg_EAngAverdE->Draw();
495 c1.Print( ( m_rootfile + ".ps" ).c_str() );
496 for ( Int_t ip = 0; ip < NumSlices; ip++ )
497 {
498 m_eangle[ip]->Draw();
499 c1.Print( ( m_rootfile + ".ps" ).c_str() );
500 }
501 for ( Int_t ip = 0; ip < NumSlices; ip++ )
502 {
503 m_trunc_eangle[ip]->Draw();
504 c1.Print( ( m_rootfile + ".ps" ).c_str() );
505 }
506 for ( Int_t ip = 0; ip < NumSlices; ip++ )
507 {
508 m_pos_eangle[ip]->Draw();
509 c1.Print( ( m_rootfile + ".ps" ).c_str() );
510 }
511 for ( Int_t ip = 0; ip < NumSlices; ip++ )
512 {
513 m_neg_eangle[ip]->Draw();
514 c1.Print( ( m_rootfile + ".ps" ).c_str() );
515 }
516 c1.Print( ( m_rootfile + ".ps]" ).c_str() );
517
518 for ( Int_t ip = 0; ip < NumSlices; ip++ )
519 {
520 delete m_eangle[ip];
521 delete m_trunc_eangle[ip];
522 delete m_pos_eangle[ip];
523 delete m_neg_eangle[ip];
524 }
525 delete m_EAngAverdE;
526 delete m_trunc_EAngAverdE;
527 delete m_pos_EAngAverdE;
528 delete m_neg_EAngAverdE;
529}
DECLARE_COMPONENT
DECLARE_COMPONENT(BesBdkRc)
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
StepEAng
const double StepEAng
Definition DedxCalibEAng.cxx:17
Landau
double Landau(double *x, double *par)
Definition DedxCalibParameters.cxx:30
Vavilov
double Vavilov(double *x, double *par)
Definition DedxCalibParameters.cxx:35
NumSlices
#define NumSlices
Definition DedxCalibParameters.h:61
landaun
double landaun(double *x, double *par)
Definition DedxCalibParameters.cxx:22
AsymGauss
double AsymGauss(double *x, double *par)
Definition DedxCalibParameters.cxx:43
PhiMax
#define PhiMax
Definition DedxCalibParameters.h:67
PhiMin
#define PhiMin
Definition DedxCalibParameters.h:66
mylan
double mylan(double *x, double *par)
Definition DedxCalibParameters.cxx:15
count
DOUBLE_PRECISION count[3]
Definition EvtPhokharaDef.hh:63
MaxADCValuecut
#define MaxADCValuecut
Definition InstallArea/x86_64-el9-gcc13-dbg/include/DedxCorrecSvc/DedxCorrParameters.h:34
Out_RPhi_PathMinCut
#define Out_RPhi_PathMinCut
Definition InstallArea/x86_64-el9-gcc13-dbg/include/DedxCorrecSvc/DedxCorrParameters.h:30
Iner_DriftDistCut
#define Iner_DriftDistCut
Definition InstallArea/x86_64-el9-gcc13-dbg/include/DedxCorrecSvc/DedxCorrParameters.h:31
RPhi_PathMaxCut
#define RPhi_PathMaxCut
Definition InstallArea/x86_64-el9-gcc13-dbg/include/DedxCorrecSvc/DedxCorrParameters.h:28
MinHistValue
#define MinHistValue
Definition InstallArea/x86_64-el9-gcc13-dbg/include/DedxCorrecSvc/DedxCorrParameters.h:23
Out_DriftDistCut
#define Out_DriftDistCut
Definition InstallArea/x86_64-el9-gcc13-dbg/include/DedxCorrecSvc/DedxCorrParameters.h:32
Iner_RPhi_PathMinCut
#define Iner_RPhi_PathMinCut
Definition InstallArea/x86_64-el9-gcc13-dbg/include/DedxCorrecSvc/DedxCorrParameters.h:29
NumHistBins
#define NumHistBins
Definition InstallArea/x86_64-el9-gcc13-dbg/include/DedxCorrecSvc/DedxCorrParameters.h:22
MaxHistValue
#define MaxHistValue
Definition InstallArea/x86_64-el9-gcc13-dbg/include/DedxCorrecSvc/DedxCorrParameters.h:24
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
DedxCalibEAng::DedxCalibEAng
DedxCalibEAng(const std::string &name, ISvcLocator *pSvcLocator)
Definition DedxCalibEAng.cxx:19
DedxCalib::DedxCalib
DedxCalib(const std::string &name, ISvcLocator *pSvcLocator)
Definition DedxCalib.cxx:12
DedxCalib::ReadRecFileList
void ReadRecFileList()
Definition DedxCalib.cxx:84
DedxCalib::m_phShapeFlag
int m_phShapeFlag
Definition DedxCalib.h:54
DedxCalib::m_rootfile
std::string m_rootfile
Definition DedxCalib.h:57
DedxCalib::exsvc
IDedxCorrecSvc * exsvc
Definition DedxCalib.h:29
DedxCalib::m_recFileVector
vector< string > m_recFileVector
Definition DedxCalib.h:50