BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
MdcCalib.cxx
Go to the documentation of this file.
1#include "MdcCalibAlg/MdcCalib.h"
2
3#include "GaudiKernel/Bootstrap.h"
4#include "GaudiKernel/IDataProviderSvc.h"
5#include "GaudiKernel/IMessageSvc.h"
6#include "GaudiKernel/INTupleSvc.h"
7#include "GaudiKernel/IService.h"
8#include "GaudiKernel/ISvcLocator.h"
9#include "GaudiKernel/MsgStream.h"
10#include "GaudiKernel/PropertyMgr.h"
11#include "GaudiKernel/SmartDataPtr.h"
12#include "GaudiKernel/StatusCode.h"
13
14#include "CLHEP/Vector/ThreeVector.h"
15#include "EventModel/Event.h"
16#include "EventModel/EventHeader.h"
17#include "Identifier/Identifier.h"
18#include "Identifier/MdcID.h"
19#include "McTruth/MdcMcHit.h"
20#include "MdcRawEvent/MdcDigi.h"
21
22#include "TStyle.h"
23
24#include <cstring>
25#include <fstream>
26#include <iomanip>
27#include <iostream>
28#include <string>
29
30using namespace Event;
31using namespace std;
32
33typedef map<int, int>::value_type valType3;
34typedef std::vector<HepLorentzVector> Vp4;
35
36MdcCalib::MdcCalib() {
37 m_nEvt = 0;
38 m_cut1 = 0;
39 m_nTrkAfTrkCut = 0;
40 m_cut2 = 0;
41 m_cut3 = 0;
42 m_cut4 = 0;
43 m_cut5 = 0;
44 m_cut6 = 0;
45 m_cut7 = 0;
46 m_nTrkCal = 0;
47 m_nGrPoint = 0;
48 fgReadWireEff = false;
49
50 int lay;
51 for ( lay = 0; lay < 43; lay++ )
52 {
53 if ( lay < 15 ) m_nEntr[lay] = 1;
54 else m_nEntr[lay] = 2;
55 // m_nEntr[lay] = 3;
56 }
57 m_dwid = 0.5; // mm
58 m_fgIni = false;
59
60 m_phiWid = PI2 / (double)NPhiBin;
61 m_theWid = 2.0 / (double)NThetaBin;
62
63 m_nEvtNtuple = 0;
64
65 for ( lay = 0; lay < MdcCalNLayer; lay++ )
66 {
67 if ( lay < 8 ) m_nBin[lay] = 12;
68 else m_nBin[lay] = 16;
69 }
70
71 // setting boundary layer flags
72 for ( lay = 0; lay < MdcCalNLayer; lay++ )
73 {
74 if ( ( 0 == lay ) || ( 7 == lay ) || ( 8 == lay ) || ( 19 == lay ) || ( 20 == lay ) ||
75 ( 35 == lay ) || ( 36 == lay ) || ( 42 == lay ) )
76 m_layBound[lay] = true;
77 else m_layBound[lay] = false;
78 }
79}
80
81MdcCalib::~MdcCalib() {}
82
83void MdcCalib::clear() {
84 int lay;
85 for ( lay = 0; lay < m_nlayer; lay++ )
86 {
87 delete m_htraw[lay];
88 delete m_htdr[lay];
89 delete m_hresInc[lay];
90 delete m_hresExc[lay];
91 delete m_hresAve[lay];
92 delete m_hadc[lay];
93 for ( int lr = 0; lr < 2; lr++ )
94 {
95 delete m_htdrlr[lay][lr];
96 delete m_hreslrInc[lay][lr];
97 delete m_hreslrExc[lay][lr];
98 delete m_hreslrAve[lay][lr];
99 }
100 }
101
102 delete m_effNtrk;
103 delete m_effNtrkRecHit;
104 delete m_hresAllInc;
105 delete m_hresAllExc;
106 delete m_hresAllAve;
107 for ( int i = 0; i < 14; i++ )
108 {
109 delete m_hresAveAllQ[i];
110 delete m_hresAveOutQ[i];
111 }
112 for ( lay = 0; lay < 43; lay++ )
113 {
114 for ( int i = 0; i < 14; i++ ) delete m_hresAveLayQ[lay][i];
115 }
116 delete m_hresInnInc;
117 delete m_hresInnExc;
118 delete m_hresStpInc;
119 delete m_hresStpExc;
120 delete m_hresOutInc;
121 delete m_hresOutExc;
122 delete m_hresAllExcUp;
123 delete m_hresAllExcDw;
124 delete m_hresInnExcUp;
125 delete m_hresInnExcDw;
126 delete m_hresStpExcUp;
127 delete m_hresStpExcDw;
128 delete m_hresOutExcUp;
129 delete m_hresOutExcDw;
130 for ( int iEs = 0; iEs < m_param.nEsFlag; iEs++ ) delete m_hTes[iEs];
131 delete m_hbbTrkFlg;
132 delete m_hTesAll;
133 delete m_hTesGood;
134 delete m_hTesAllFlag;
135 delete m_hTesRec;
136 delete m_hTesCalFlag;
137 delete m_hTesCalUse;
138 delete m_hnRawHit;
139 delete m_hpt;
140 delete m_hpMax;
141 delete m_hpMaxCms;
142 delete m_hptPos;
143 delete m_hptNeg;
144 delete m_hp;
145 delete m_hp_cms;
146 delete m_hpPos;
147 delete m_hpNeg;
148 delete m_hpPoscms;
149 delete m_hpNegcms;
150 delete m_hp_cut;
151 delete m_hchisq;
152 delete m_hnTrk;
153 delete m_hnTrkCal;
154 delete m_hnhitsRec;
155 delete m_hnhitsRecInn;
156 delete m_hnhitsRecStp;
157 delete m_hnhitsRecOut;
158 delete m_hnhitsCal;
159 delete m_hnhitsCalInn;
160 delete m_hnhitsCalStp;
161 delete m_hnhitsCalOut;
162 delete m_wirehitmap;
163 delete m_layerhitmap;
164 delete m_hnoisephi;
165 delete m_hnoiselay;
166 delete m_hnoisenhits;
167 delete m_hratio;
168 delete m_hdr;
169 delete m_hphi0;
170 delete m_hkap;
171 delete m_hdz;
172 delete m_htanl;
173 delete m_hcosthe;
174 delete m_hcostheNeg;
175 delete m_hcosthePos;
176 delete m_hx0;
177 delete m_hy0;
178 delete m_hdelZ0;
179 delete m_grX0Y0;
180 delete m_hitEffAll;
181 delete m_hitEffRaw;
182 delete m_hitEffRec;
183 int bin;
184 int thbin;
185 for ( bin = 0; bin < NPhiBin; bin++ )
186 {
187 delete m_ppPhi[bin];
188 delete m_pnPhi[bin];
189 delete m_ppPhiCms[bin];
190 delete m_pnPhiCms[bin];
191 for ( thbin = 0; thbin < NThetaBin; thbin++ )
192 {
193 delete m_ppThePhi[thbin][bin];
194 delete m_pnThePhi[thbin][bin];
195 delete m_ppThePhiCms[thbin][bin];
196 delete m_pnThePhiCms[thbin][bin];
197 }
198 }
199 for ( thbin = 0; thbin < NThetaBin; thbin++ )
200 {
201 delete m_ppThe[thbin];
202 delete m_pnThe[thbin];
203 delete m_ppTheCms[thbin];
204 delete m_pnTheCms[thbin];
205 }
206
207 for ( unsigned i = 0; i < m_hr2dInc.size(); i++ )
208 {
209 delete m_hr2dInc[i];
210 delete m_hr2dExc[i];
211 }
212 m_hr2dInc.clear();
213 m_hr2dExc.clear();
214 m_mapr2d.clear();
215
216 delete m_fdTime;
217 delete m_fdAdc;
218 delete m_fdres;
219 delete m_fdresAve;
220 delete m_fdres2d;
221 delete m_fdcom;
222 delete m_fdResQ;
223}
224
225void MdcCalib::initialize( TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
226 IMdcCalibFunSvc* mdcFunSvc, IMdcUtilitySvc* mdcUtilitySvc ) {
227 IMessageSvc* msgSvc;
228 Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
229 MsgStream log( msgSvc, "MdcCalib" );
230 log << MSG::INFO << "MdcCalib::initialize()" << endmsg;
231
232 m_hlist = hlist;
233 m_mdcGeomSvc = mdcGeomSvc;
234 m_mdcFunSvc = mdcFunSvc;
235 m_mdcUtilitySvc = mdcUtilitySvc;
236
237 int lay;
238 int iEntr;
239 int lr;
240 int bin;
241 char hname[200];
242
243 int nbinMom = 2400;
244
245 m_nlayer = m_mdcGeomSvc->getLayerSize();
246
247 for ( lay = 0; lay < m_nlayer; lay++ )
248 { m_radii[lay] = m_mdcGeomSvc->Layer( lay )->Radius(); }
249 ofstream fwpc( "wirelog.txt" );
250 for ( int wir = 0; wir < MdcCalTotCell; wir++ )
251 {
252 m_xe[wir] = m_mdcGeomSvc->Wire( wir )->Backward().x();
253 m_ye[wir] = m_mdcGeomSvc->Wire( wir )->Backward().y();
254 m_ze[wir] = m_mdcGeomSvc->Wire( wir )->Backward().z();
255 m_xw[wir] = m_mdcGeomSvc->Wire( wir )->Forward().x();
256 m_yw[wir] = m_mdcGeomSvc->Wire( wir )->Forward().y();
257 m_zw[wir] = m_mdcGeomSvc->Wire( wir )->Forward().z();
258 fwpc << setw( 6 ) << wir << setw( 15 ) << m_xe[wir] << setw( 15 ) << m_ye[wir]
259 << setw( 15 ) << m_ze[wir] << setw( 15 ) << m_xw[wir] << setw( 15 ) << m_yw[wir]
260 << setw( 15 ) << m_zw[wir] << endl;
261 }
262 fwpc.close();
263
264 m_fdcom = new TFolder( "common", "common" );
265 m_hlist->Add( m_fdcom );
266
267 m_effNtrk = new TH1F( "effNtrk", "", 43, -0.5, 42.5 );
268 m_fdcom->Add( m_effNtrk );
269
270 m_effNtrkRecHit = new TH1F( "effNtrkRecHit", "", 43, -0.5, 42.5 );
271 m_fdcom->Add( m_effNtrkRecHit );
272
273 m_hresAllInc = new TH1F( "HResAllInc", "", 200, -1.0, 1.0 );
274 m_fdcom->Add( m_hresAllInc );
275
276 m_hresAllExc = new TH1F( "HResAllExc", "", 200, -1.0, 1.0 );
277 m_fdcom->Add( m_hresAllExc );
278
279 m_hresAllAve = new TH1F( "HResAllAve", "", 200, -1.0, 1.0 );
280 m_fdcom->Add( m_hresAllAve );
281
282 m_hresInnInc = new TH1F( "HResInnInc", "", 200, -1.0, 1.0 );
283 m_fdcom->Add( m_hresInnInc );
284
285 m_hresInnExc = new TH1F( "HResInnExc", "", 200, -1.0, 1.0 );
286 m_fdcom->Add( m_hresInnExc );
287
288 m_hresStpInc = new TH1F( "HResStpInc", "", 200, -1.0, 1.0 );
289 m_fdcom->Add( m_hresStpInc );
290
291 m_hresStpExc = new TH1F( "HResStpExc", "", 200, -1.0, 1.0 );
292 m_fdcom->Add( m_hresStpExc );
293
294 m_hresOutInc = new TH1F( "HResOutInc", "", 200, -1.0, 1.0 );
295 m_fdcom->Add( m_hresOutInc );
296
297 m_hresOutExc = new TH1F( "HResOutExc", "", 200, -1.0, 1.0 );
298 m_fdcom->Add( m_hresOutExc );
299
300 m_hresAllExcUp = new TH1F( "HResAllExcUp", "", 200, -1.0, 1.0 );
301 m_fdcom->Add( m_hresAllExcUp );
302
303 m_hresAllExcDw = new TH1F( "HResAllExcDw", "", 200, -1.0, 1.0 );
304 m_fdcom->Add( m_hresAllExcDw );
305
306 m_hresInnExcUp = new TH1F( "HResInnExcUp", "", 200, -1.0, 1.0 );
307 m_fdcom->Add( m_hresInnExcUp );
308
309 m_hresInnExcDw = new TH1F( "HResInnExcDw", "", 200, -1.0, 1.0 );
310 m_fdcom->Add( m_hresInnExcDw );
311
312 m_hresStpExcUp = new TH1F( "HResStpExcUp", "", 200, -1.0, 1.0 );
313 m_fdcom->Add( m_hresStpExcUp );
314
315 m_hresStpExcDw = new TH1F( "HResStpExcDw", "", 200, -1.0, 1.0 );
316 m_fdcom->Add( m_hresStpExcDw );
317
318 m_hresOutExcUp = new TH1F( "HResOutExcUp", "", 200, -1.0, 1.0 );
319 m_fdcom->Add( m_hresOutExcUp );
320
321 m_hresOutExcDw = new TH1F( "HResOutExcDw", "", 200, -1.0, 1.0 );
322 m_fdcom->Add( m_hresOutExcDw );
323
324 m_fdResQ = new TFolder( "ResQ", "ResQ" );
325 m_hlist->Add( m_fdResQ );
326 for ( int i = 0; i < 14; i++ )
327 {
328 sprintf( hname, "resoAll_qbin%02d", i );
329 m_hresAveAllQ[i] = new TH1F( hname, "", 200, -1, 1 );
330 m_fdResQ->Add( m_hresAveAllQ[i] );
331
332 sprintf( hname, "resoOut_qbin%02d", i );
333 m_hresAveOutQ[i] = new TH1F( hname, "", 200, -1, 1 );
334 m_fdResQ->Add( m_hresAveOutQ[i] );
335 }
336 for ( lay = 0; lay < 43; lay++ )
337 {
338 for ( int i = 0; i < 14; i++ )
339 {
340 sprintf( hname, "resoLay%02d_qbin%02d", lay, i );
341 m_hresAveLayQ[lay][i] = new TH1F( hname, "", 200, -1, 1 );
342 m_fdResQ->Add( m_hresAveLayQ[lay][i] );
343 }
344 }
345
346 for ( int iEs = 0; iEs < m_param.nEsFlag; iEs++ )
347 {
348 sprintf( hname, "Tes_%d", m_param.esFlag[iEs] );
349 m_hTes[iEs] = new TH1F( hname, "", 750, 0, 1500 );
350 m_fdcom->Add( m_hTes[iEs] );
351 }
352
353 m_hbbTrkFlg = new TH1F( "BbTrkFlg", "", 100, 0, 6 );
354 m_fdcom->Add( m_hbbTrkFlg );
355
356 m_hTesAll = new TH1F( "TesAll", "", 2000, 0, 2000 );
357 m_fdcom->Add( m_hTesAll );
358
359 m_hTesGood = new TH1F( "TesGood", "", 2000, 0, 2000 );
360 m_fdcom->Add( m_hTesGood );
361
362 m_hTesAllFlag = new TH1F( "TesAllFlag", "", 300, -0.5, 299.5 );
363 m_fdcom->Add( m_hTesAllFlag );
364
365 m_hTesRec = new TH1F( "TesRec", "", 2000, 0, 2000 );
366 m_fdcom->Add( m_hTesRec );
367
368 m_hTesCalFlag = new TH1F( "TesCalFlag", "", 2000, 0, 2000 );
369 m_fdcom->Add( m_hTesCalFlag );
370
371 m_hTesCalUse = new TH1F( "TesCalUse", "", 2000, 0, 2000 );
372 m_fdcom->Add( m_hTesCalUse );
373
374 m_hnRawHit = new TH1F( "nRawHit", "", 6797, -0.5, 6796.5 );
375 m_fdcom->Add( m_hnRawHit );
376
377 m_hpt = new TH1F( "HPt", "", nbinMom, 0, 3 );
378 m_fdcom->Add( m_hpt );
379
380 m_hptPos = new TH1F( "HPtPos", "", nbinMom, 0, 3 );
381 m_fdcom->Add( m_hptPos );
382
383 m_hptNeg = new TH1F( "HPtNeg", "", nbinMom, 0, 3 );
384 m_fdcom->Add( m_hptNeg );
385
386 m_hp = new TH1F( "HP", "", nbinMom, 0, 3 );
387 m_fdcom->Add( m_hp );
388
389 m_hp_cms = new TH1F( "HPCMS", "", nbinMom, 0, 3 );
390 m_fdcom->Add( m_hp_cms );
391
392 m_hpMax = new TH1F( "HPMax", "", nbinMom, 0, 3 );
393 m_fdcom->Add( m_hpMax );
394
395 m_hpMaxCms = new TH1F( "HPMax_Cms", "", nbinMom, 0, 3 );
396 m_fdcom->Add( m_hpMaxCms );
397
398 m_hpPos = new TH1F( "HP_Pos", "", nbinMom, 0, 3 );
399 m_fdcom->Add( m_hpPos );
400
401 m_hpNeg = new TH1F( "HP_Neg", "", nbinMom, 0, 3 );
402 m_fdcom->Add( m_hpNeg );
403
404 m_hpPoscms = new TH1F( "HP_Pos_cms", "", nbinMom, 0, 3 );
405 m_fdcom->Add( m_hpPoscms );
406
407 m_hpNegcms = new TH1F( "HP_Neg_cms", "", nbinMom, 0, 3 );
408 m_fdcom->Add( m_hpNegcms );
409
410 m_hp_cut = new TH1F( "HPCut", "", nbinMom, 0, 3 );
411 m_fdcom->Add( m_hp_cut );
412
413 m_hchisq = new TH1F( "Chisq", "", 10, 0, 100 );
414 m_fdcom->Add( m_hchisq );
415
416 m_hnTrk = new TH1F( "HNtrack", "HNtrack", 10, -0.5, 9.5 );
417 m_fdcom->Add( m_hnTrk );
418
419 m_hnTrkCal = new TH1F( "HNtrackCal", "HNtrackCal", 10, -0.5, 9.5 );
420 m_fdcom->Add( m_hnTrkCal );
421
422 m_hnhitsRec = new TH1F( "HNhitsRec", "", 100, -0.5, 99.5 );
423 m_fdcom->Add( m_hnhitsRec );
424
425 m_hnhitsRecInn = new TH1F( "HNhitsInnRec", "", 60, 0.5, 60.5 );
426 m_fdcom->Add( m_hnhitsRecInn );
427
428 m_hnhitsRecStp = new TH1F( "HNhitsStpRec", "", 60, 0.5, 60.5 );
429 m_fdcom->Add( m_hnhitsRecStp );
430
431 m_hnhitsRecOut = new TH1F( "HNhitsOutRec", "", 60, 0.5, 60.5 );
432 m_fdcom->Add( m_hnhitsRecOut );
433
434 m_hnhitsCal = new TH1F( "HNhitsCal", "", 100, -0.5, 99.5 );
435 m_fdcom->Add( m_hnhitsCal );
436
437 m_hnhitsCalInn = new TH1F( "HNhitsCalInn", "", 60, 0.5, 60.5 );
438 m_fdcom->Add( m_hnhitsCalInn );
439
440 m_hnhitsCalStp = new TH1F( "HNhitsCalStp", "", 60, 0.5, 60.5 );
441 m_fdcom->Add( m_hnhitsCalStp );
442
443 m_hnhitsCalOut = new TH1F( "HNhitsCalOut", "", 60, 0.5, 60.5 );
444 m_fdcom->Add( m_hnhitsCalOut );
445
446 m_wirehitmap = new TH1F( "Wire_HitMap", "Wire_HitMap", 6796, -0.5, 6795.5 );
447 m_fdcom->Add( m_wirehitmap );
448
449 m_layerhitmap = new TH1F( "Layer_HitMap", "Layer_HitMap", 43, -0.5, 42.5 );
450 m_fdcom->Add( m_layerhitmap );
451
452 m_hnoisephi = new TH1F( "phi_noise", "", 100, 0, 6.284 );
453 m_fdcom->Add( m_hnoisephi );
454
455 m_hnoiselay = new TH1F( "Layer_noise", "Layer_noise", 43, -0.5, 42.5 );
456 m_fdcom->Add( m_hnoiselay );
457
458 m_hnoisenhits = new TH1F( "nhits_noise", "nhits_noise", 6796, -0.5, 6795.5 );
459 m_fdcom->Add( m_hnoisenhits );
460
461 m_hratio = new TH1F( "ratio", "", 100, 0, 1 );
462 m_fdcom->Add( m_hratio );
463
464 m_hdr = new TH1F( "dr", "", 2000, -100, 100 );
465 m_fdcom->Add( m_hdr );
466
467 m_hphi0 = new TH1F( "phi0", "", 100, 0, 6.284 );
468 m_fdcom->Add( m_hphi0 );
469
470 m_hkap = new TH1F( "kappa", "", 400, -50, 50 );
471 m_fdcom->Add( m_hkap );
472
473 m_hdz = new TH1F( "dz", "", 2000, -500, 500 );
474 m_fdcom->Add( m_hdz );
475
476 m_htanl = new TH1F( "tanl", "", 200, -5, 5 );
477 m_fdcom->Add( m_htanl );
478
479 m_hcosthe = new TH1F( "costheta", "", 200, -1, 1 );
480 m_fdcom->Add( m_hcosthe );
481
482 m_hcostheNeg = new TH1F( "costhetaNeg", "", 200, -1, 1 );
483 m_fdcom->Add( m_hcostheNeg );
484
485 m_hcosthePos = new TH1F( "costhetaPos", "", 200, -1, 1 );
486 m_fdcom->Add( m_hcosthePos );
487
488 m_hx0 = new TH1F( "x0", "", 100, -10, 10 );
489 m_fdcom->Add( m_hx0 );
490
491 m_hy0 = new TH1F( "y0", "", 100, -10, 10 );
492 m_fdcom->Add( m_hy0 );
493
494 m_hdelZ0 = new TH1F( "delta_z0", "", 100, -50, 50 );
495 m_fdcom->Add( m_hdelZ0 );
496
497 m_grX0Y0 = new TGraph();
498 m_grX0Y0->SetName( "x0y0" );
499 m_fdcom->Add( m_grX0Y0 );
500
501 m_hitEffAll = new TH1F( "hitEffAll", "", 6800, -0.5, 6799.5 );
502 m_fdcom->Add( m_hitEffAll );
503
504 m_hitEffRaw = new TH1F( "hitEffRaw", "", 6800, -0.5, 6799.5 );
505 m_fdcom->Add( m_hitEffRaw );
506
507 m_hitEffRec = new TH1F( "hitEffRec", "", 6800, -0.5, 6799.5 );
508 m_fdcom->Add( m_hitEffRec );
509
510 // histograms for drift time
511 m_fdTime = new TFolder( "time", "time" );
512 m_hlist->Add( m_fdTime );
513
514 for ( lay = 0; lay < m_nlayer; lay++ )
515 {
516 sprintf( hname, "Traw%02d", lay );
517 m_htraw[lay] = new TH1F( hname, "", 1000, 0, 2000 );
518 m_fdTime->Add( m_htraw[lay] );
519
520 sprintf( hname, "Tdr%02d", lay );
521 m_htdr[lay] = new TH1F( hname, "", 510, -10, 500 );
522 m_fdTime->Add( m_htdr[lay] );
523
524 for ( lr = 0; lr < 2; lr++ )
525 {
526 sprintf( hname, "Tdr%02d_lr%01d", lay, lr );
527 m_htdrlr[lay][lr] = new TH1F( hname, "", 510, -10, 500 );
528 m_fdTime->Add( m_htdrlr[lay][lr] );
529 }
530 }
531
532 // histograms of adc
533 m_fdAdc = new TFolder( "adc", "adc" );
534 m_hlist->Add( m_fdAdc );
535
536 for ( lay = 0; lay < m_nlayer; lay++ )
537 {
538 sprintf( hname, "adc%02d", lay );
539 m_hadc[lay] = new TH1F( hname, "", 1500, 0, 3000 );
540 m_fdAdc->Add( m_hadc[lay] );
541 }
542 // histograms for resolution
543 m_fdres = new TFolder( "resolution", "resolution" );
544 m_hlist->Add( m_fdres );
545
546 m_fdresAve = new TFolder( "resAve", "resAve" );
547 m_hlist->Add( m_fdresAve );
548
549 for ( lay = 0; lay < m_nlayer; lay++ )
550 {
551 sprintf( hname, "Reso%02dInc", lay );
552 m_hresInc[lay] = new TH1F( hname, "", 1000, -5, 5 );
553 m_fdres->Add( m_hresInc[lay] );
554
555 sprintf( hname, "Reso%02dExc", lay );
556 m_hresExc[lay] = new TH1F( hname, "", 1000, -5, 5 );
557 m_fdres->Add( m_hresExc[lay] );
558
559 sprintf( hname, "Reso%02d", lay );
560 m_hresAve[lay] = new TH1F( hname, "", 1000, -5, 5 );
561 m_fdresAve->Add( m_hresAve[lay] );
562
563 for ( lr = 0; lr < 2; lr++ )
564 {
565 sprintf( hname, "Reso%02dInc_lr%01d", lay, lr );
566 // m_hreslrInc[lay][lr] = new TH1F(hname, "", 200, -1, 1);
567 m_hreslrInc[lay][lr] = new TH1F( hname, "", 1000, -5, 5 );
568 m_fdres->Add( m_hreslrInc[lay][lr] );
569
570 sprintf( hname, "Reso%02dExc_lr%01d", lay, lr );
571 // m_hreslrExc[lay][lr] = new TH1F(hname, "", 200, -1, 1);
572 m_hreslrExc[lay][lr] = new TH1F( hname, "", 1000, -5, 5 );
573 m_fdres->Add( m_hreslrExc[lay][lr] );
574
575 sprintf( hname, "Reso%02d_lr%01d", lay, lr );
576 // m_hreslrAve[lay][lr] = new TH1F(hname, "", 200, -1, 1);
577 m_hreslrAve[lay][lr] = new TH1F( hname, "", 1000, -5, 5 );
578 m_fdresAve->Add( m_hreslrAve[lay][lr] );
579 }
580 for ( int phi = 0; phi < 20; phi++ )
581 {
582 sprintf( hname, "ResoPhi%02d_phi%02d", lay, phi );
583 m_hresphi[lay][phi] = new TH1F( hname, "", 200, -1, 1 );
584 m_fdres->Add( m_hresphi[lay][phi] );
585 }
586 }
587
588 /* histograms for momentum vs phi */
589 m_fdmomPhi = new TFolder( "momPhi", "momPhi" );
590 m_hlist->Add( m_fdmomPhi );
591
592 int thbin;
593 for ( bin = 0; bin < NPhiBin; bin++ )
594 {
595 sprintf( hname, "hPpos_phi%02d", bin );
596 m_ppPhi[bin] = new TH1F( hname, "", nbinMom, 0.0, 3.0 );
597 m_fdmomPhi->Add( m_ppPhi[bin] );
598
599 sprintf( hname, "hPneg_phi%02d", bin );
600 m_pnPhi[bin] = new TH1F( hname, "", nbinMom, 0.0, 3.0 );
601 m_fdmomPhi->Add( m_pnPhi[bin] );
602
603 sprintf( hname, "hPpos_phi_cms%02d", bin );
604 m_ppPhiCms[bin] = new TH1F( hname, "", nbinMom, 0.0, 3.0 );
605 m_fdmomPhi->Add( m_ppPhiCms[bin] );
606
607 sprintf( hname, "hPneg_phi_cms%02d", bin );
608 m_pnPhiCms[bin] = new TH1F( hname, "", nbinMom, 0.0, 3.0 );
609 m_fdmomPhi->Add( m_pnPhiCms[bin] );
610
611 for ( thbin = 0; thbin < NThetaBin; thbin++ )
612 {
613 sprintf( hname, "hPpos_theta%02d_phi%02d", thbin, bin );
614 m_ppThePhi[thbin][bin] = new TH1F( hname, "", nbinMom, 0.0, 3.0 );
615 m_fdmomPhi->Add( m_ppThePhi[thbin][bin] );
616
617 sprintf( hname, "hPneg_theta%02d_phi%02d", thbin, bin );
618 m_pnThePhi[thbin][bin] = new TH1F( hname, "", nbinMom, 0.0, 3.0 );
619 m_fdmomPhi->Add( m_pnThePhi[thbin][bin] );
620
621 sprintf( hname, "hPposCms_theta%02d_phi%02d", thbin, bin );
622 m_ppThePhiCms[thbin][bin] = new TH1F( hname, "", nbinMom, 0.0, 3.0 );
623 m_fdmomPhi->Add( m_ppThePhiCms[thbin][bin] );
624
625 sprintf( hname, "hPnegCms_theta%02d_phi%02d", thbin, bin );
626 m_pnThePhiCms[thbin][bin] = new TH1F( hname, "", nbinMom, 0.0, 3.0 );
627 m_fdmomPhi->Add( m_pnThePhiCms[thbin][bin] );
628 }
629 }
630 for ( thbin = 0; thbin < NThetaBin; thbin++ )
631 {
632 sprintf( hname, "hPpos_the%02d", thbin );
633 m_ppThe[thbin] = new TH1F( hname, "", nbinMom, 0.0, 3.0 );
634 m_fdmomPhi->Add( m_ppThe[thbin] );
635
636 sprintf( hname, "hPneg_the%02d", thbin );
637 m_pnThe[thbin] = new TH1F( hname, "", nbinMom, 0.0, 3.0 );
638 m_fdmomPhi->Add( m_pnThe[thbin] );
639
640 sprintf( hname, "hPposCms_the%02d", thbin );
641 m_ppTheCms[thbin] = new TH1F( hname, "", nbinMom, 0.0, 3.0 );
642 m_fdmomPhi->Add( m_ppTheCms[thbin] );
643
644 sprintf( hname, "hPnegCms_the%02d", thbin );
645 m_pnTheCms[thbin] = new TH1F( hname, "", nbinMom, 0.0, 3.0 );
646 m_fdmomPhi->Add( m_pnTheCms[thbin] );
647 }
648
649 // histograms for resolution vs distance
650 m_fdres2d = new TFolder( "res2d", "res2d" );
651 m_hlist->Add( m_fdres2d );
652
653 int hid = 0;
654 int key;
655 TH1F* hist;
656 for ( lay = 0; lay < m_nlayer; lay++ )
657 {
658 for ( iEntr = 0; iEntr < MdcCalNENTRSD; iEntr++ )
659 {
660 for ( lr = 0; lr < 2; lr++ )
661 {
662 for ( bin = 0; bin < MdcCalSdNBIN; bin++ )
663 {
664 sprintf( hname, "r2d%02d_%02d_%01d_%02dInc", lay, iEntr, lr, bin );
665 hist = new TH1F( hname, "", 200, -1, 1 );
666 m_hr2dInc.push_back( hist );
667 m_fdres2d->Add( hist );
668
669 sprintf( hname, "r2d%02d_%02d_%01d_%02dExc", lay, iEntr, lr, bin );
670 hist = new TH1F( hname, "", 200, -1, 1 );
671 m_hr2dExc.push_back( hist );
672 m_fdres2d->Add( hist );
673
674 key = getHresId( lay, iEntr, lr, bin );
675 m_mapr2d.insert( valType3( key, hid ) );
676 hid++;
677 }
678 }
679 }
680 } // end of layer loop
681
682 m_fdres2t = new TFolder( "res2t", "res2t" );
683 // m_hlist -> Add(m_fdres2t);
684
685 for ( lay = 0; lay < m_nlayer; lay++ )
686 {
687 for ( iEntr = 0; iEntr < MdcCalNENTRXT; iEntr++ )
688 {
689 for ( lr = 0; lr < 2; lr++ )
690 {
691 for ( bin = 0; bin < 45; bin++ )
692 {
693 sprintf( hname, "r2t%02d_%02d_%01d_%02d", lay, iEntr, lr, bin );
694 m_hr2t[lay][iEntr][lr][bin] = new TH1F( hname, "", 600, -3, 3 );
695 m_fdres2t->Add( m_hr2t[lay][iEntr][lr][bin] );
696 }
697 }
698 }
699 }
700
701 INTupleSvc* ntupleSvc;
702 Gaudi::svcLocator()->service( "NTupleSvc", ntupleSvc );
703 for ( lay = 0; lay < m_nlayer; lay++ )
704 {
705 sprintf( hname, "FILE136/xt%02d", lay );
706 NTuplePtr nt( ntupleSvc, hname );
707 if ( nt ) m_xtTuple[lay] = nt;
708 else
709 {
710 m_xtTuple[lay] = ntupleSvc->book( hname, CLID_ColumnWiseTuple, "MdcXtNtuple" );
711 if ( m_xtTuple[lay] )
712 {
713 m_xtTuple[lay]->addItem( "cel", m_cel[lay] );
714 m_xtTuple[lay]->addItem( "lr", m_lr[lay] );
715 m_xtTuple[lay]->addItem( "run", m_run[lay] );
716 m_xtTuple[lay]->addItem( "evt", m_evt[lay] );
717 m_xtTuple[lay]->addItem( "doca", m_doca[lay] );
718 m_xtTuple[lay]->addItem( "dm", m_dm[lay] );
719 m_xtTuple[lay]->addItem( "tdr", m_tdr[lay] );
720 m_xtTuple[lay]->addItem( "tdc", m_tdc[lay] );
721 m_xtTuple[lay]->addItem( "entr", m_entr[lay] );
722 m_xtTuple[lay]->addItem( "zhit", m_zhit[lay] );
723 m_xtTuple[lay]->addItem( "qhit", m_qhit[lay] );
724 m_xtTuple[lay]->addItem( "p", m_p[lay] );
725 m_xtTuple[lay]->addItem( "pt", m_pt[lay] );
726 m_xtTuple[lay]->addItem( "phi0", m_phi0[lay] );
727 m_xtTuple[lay]->addItem( "tanl", m_tanl[lay] );
728 m_xtTuple[lay]->addItem( "hitphi", m_hitphi[lay] );
729 }
730 else
731 { log << MSG::FATAL << "Cannot book Xt N-tuple:" << long( m_xtTuple[lay] ) << endmsg; }
732 }
733 }
734
735 if ( 5 == m_param.particle )
736 {
737 sprintf( hname, "FILE136/cosmic" );
738 NTuplePtr nt( ntupleSvc, hname );
739 if ( nt ) m_cosmic = nt;
740 else
741 {
742 m_cosmic = ntupleSvc->book( hname, CLID_ColumnWiseTuple, "MdcXtNtuple" );
743 if ( m_cosmic )
744 {
745 m_cosmic->addItem( "pUp", m_pUpcos );
746 m_cosmic->addItem( "pDw", m_pDwcos );
747 m_cosmic->addItem( "ptUp", m_ptUpcos );
748 m_cosmic->addItem( "ptDw", m_ptDwcos );
749 m_cosmic->addItem( "phiUp", m_phiUpcos );
750 m_cosmic->addItem( "phiDw", m_phiDwcos );
751 m_cosmic->addItem( "drUp", m_drUpcos );
752 m_cosmic->addItem( "drDw", m_drDwcos );
753 m_cosmic->addItem( "dzUp", m_dzUpcos );
754 m_cosmic->addItem( "dzDw", m_dzDwcos );
755 m_cosmic->addItem( "ctheUp", m_ctheUpcos );
756 m_cosmic->addItem( "ctheDw", m_ctheDwcos );
757 m_cosmic->addItem( "nhitUp", m_nhitUpcos );
758 m_cosmic->addItem( "nhitDw", m_nhitDwcos );
759 m_cosmic->addItem( "char", m_chargecos );
760 m_cosmic->addItem( "tesfg", m_tesFlagcos );
761 m_cosmic->addItem( "tes", m_tescos );
762 }
763 }
764 }
765}
766
767int MdcCalib::fillHist( MdcCalEvent* event ) {
768 IMessageSvc* msgSvc;
769 Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
770 MsgStream log( msgSvc, "MdcCalib" );
771 log << MSG::DEBUG << "MdcCalib::fillHist()" << endmsg;
772
773 int i;
774 int k;
775 int lay;
776 int cel;
777 int wir;
778 int lr;
779 int stat;
780
781 int hid;
782 int key;
783 int iEntr;
784 int bin;
785
786 int phiBin;
787 int phiBinCms;
788 int theBin;
789 int theBinCms;
790 double lamda;
791 double theta;
792
793 double qhit;
794 double traw;
795 double tdr;
796 double doca;
797 double resiInc;
798 double resiExc;
799 double entr;
800 double pt;
801 double p;
802 double p_cms;
803 double chisq;
804 double ecm = m_param.ecm;
805 double xboost = m_param.boostPar[0] * ecm;
806 double yboost = m_param.boostPar[1];
807 double zboost = m_param.boostPar[2];
808 HepLorentzVector p4;
809
810 double dr;
811 double phi0;
812 double dz;
813 double kap;
814 double tanl;
815
816 double x0;
817 double y0;
818 double zminus = 9999.0;
819 double zplus = -9999.0;
820
821 double hitphi;
822 double philab;
823 double phicms;
824 double thetacms;
825 double costheCMS;
826
827 double dphi;
828 double wphi;
829 double xx;
830 double yy;
831 double rr;
832
833 int nhitlay;
834 bool fgHitLay[MdcCalNLayer];
835 bool fgTrk;
836
837 int ntrk = event->getNTrk();
838 int nhitRec;
839 int nhitRecInn;
840 int nhitRecStp;
841 int nhitRecOut;
842 int nhitCal;
843 int nhitCalInn;
844 int nhitCalStp;
845 int nhitCalOut;
846 MdcCalRecTrk* rectrk;
847 MdcCalRecHit* rechit;
848
849 int fgAdd[43]; // for calculating layer efficiency
850
851 // read dead wire
852 if ( !fgReadWireEff )
853 {
854 for ( lay = 0; lay < MdcCalNLayer; lay++ )
855 {
856 int ncel = m_mdcGeomSvc->Layer( lay )->NCell();
857 for ( cel = 0; cel < ncel; cel++ )
858 {
859 double eff = m_mdcFunSvc->getWireEff( lay, cel );
860 if ( eff > 0.5 ) m_fgGoodWire[lay][cel] = true;
861 else m_fgGoodWire[lay][cel] = false;
862 int wireid = m_mdcGeomSvc->Wire( lay, cel )->Id();
863 if ( eff < 0.9 )
864 cout << "dead channel: " << setw( 5 ) << lay << setw( 5 ) << cel << setw( 8 )
865 << wireid << endl;
866 }
867 }
868 fgReadWireEff = true;
869
870 ofstream fwpc( "wirelog_align.txt" );
871 for ( wir = 0; wir < MdcCalTotCell; wir++ )
872 {
873 m_xe[wir] = m_mdcGeomSvc->Wire( wir )->Backward().x();
874 m_ye[wir] = m_mdcGeomSvc->Wire( wir )->Backward().y();
875 m_ze[wir] = m_mdcGeomSvc->Wire( wir )->Backward().z();
876 m_xw[wir] = m_mdcGeomSvc->Wire( wir )->Forward().x();
877 m_yw[wir] = m_mdcGeomSvc->Wire( wir )->Forward().y();
878 m_zw[wir] = m_mdcGeomSvc->Wire( wir )->Forward().z();
879 fwpc << setw( 6 ) << wir << setw( 15 ) << m_xe[wir] << setw( 15 ) << m_ye[wir]
880 << setw( 15 ) << m_ze[wir] << setw( 15 ) << m_xw[wir] << setw( 15 ) << m_yw[wir]
881 << setw( 15 ) << m_zw[wir] << endl;
882 }
883 fwpc.close();
884 }
885
886 int nRawHit = event->getNRawHitTQ();
887 m_hnRawHit->Fill( nRawHit );
888
889 IDataProviderSvc* eventSvc = NULL;
890 Gaudi::svcLocator()->service( "EventDataSvc", eventSvc );
891
892 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc, "/Event/EventHeader" );
893 if ( !eventHeader )
894 {
895 log << MSG::FATAL << "Could not find Event Header" << endmsg;
896 return ( StatusCode::FAILURE );
897 }
898 int iRun = eventHeader->runNumber();
899 int iEvt = eventHeader->eventNumber();
900
901 int esTimeflag = event->getEsFlag();
902 double tes = event->getTes();
903 bool esCutFg = event->getEsCutFlag();
904 int iEs = event->getNesCutFlag();
905 // calculate the efficiency of Bhabha event
906 if ( -1 != esTimeflag )
907 {
908 SmartDataPtr<RecMdcTrackCol> newtrkCol( eventSvc, "/Event/Recon/RecMdcTrackCol" );
909 if ( !newtrkCol )
910 {
911 log << MSG::ERROR << "Could not find RecMdcTrackCol" << endmsg;
912 return ( StatusCode::FAILURE );
913 }
914 int nGoodTrk = 0;
915 int icharge = 0;
916 Vp4 p4p;
917 Vp4 p4m;
918 RecMdcTrackCol::iterator it_trk = newtrkCol->begin();
919 for ( ; it_trk != newtrkCol->end(); it_trk++ )
920 {
921 double mass = 0.000511; // only for eletron
922 HepLorentzVector p4 = ( *it_trk )->p4( mass );
923 icharge = ( *it_trk )->charge();
924 if ( icharge > 0 ) p4p.push_back( p4 );
925 else p4m.push_back( p4 );
926 }
927 if ( 1 == p4p.size() * p4m.size() )
928 {
929 double dang = p4p[0].vect().angle( p4m[0].vect() );
930 m_hbbTrkFlg->Fill( 1 );
931 if ( dang > 2.94 ) { m_hbbTrkFlg->Fill( 2 ); }
932 }
933 }
934 m_hTesAll->Fill( tes );
935 // cout << "tes " << tes << endl;
936 // if (-1 != esTimeflag) m_hTesGood->Fill(tes);
937 if ( -1 == esTimeflag ) return -1;
938 m_hTesGood->Fill( tes );
939 m_hTesAllFlag->Fill( esTimeflag );
940 if ( ntrk > 0 ) m_hTesRec->Fill( tes );
941 if ( ( iEs >= 0 ) && ( iEs < m_param.nEsFlag ) ) m_hTes[iEs]->Fill( tes );
942
943 if ( esCutFg ) m_hTesCalFlag->Fill( tes );
944 else return -1;
945
946 if ( !m_fgIni )
947 {
948 for ( lay = 0; lay < MdcCalNLayer; lay++ )
949 {
950 if ( lay < 8 ) m_docaMax[lay] = m_param.maxDocaInner;
951 else m_docaMax[lay] = m_param.maxDocaOuter;
952 }
953 m_fgIni = true;
954 }
955
956 bool trkFlag[200]; // for calculating hit efficiency without impact of track fitting
957 for ( i = 0; i < 200; i++ ) trkFlag[i] = false;
958
959 int ntrkCal = 0;
960 double pTrk[2];
961 double pTrkcms[2];
962 double hitphiplus = 9999.0;
963 double hitthetaplus = 9999.0;
964 double hitphiminus = -9999.0;
965 double hitthetaminus = -9999.0;
966 Vp4 pp;
967 Vp4 pm;
968 pp.clear();
969 pm.clear();
970 int phibinp;
971 int phibinm;
972
973 m_hnTrk->Fill( ntrk );
974 if ( ( ntrk < m_param.nTrkCut[0] ) || ( ntrk > m_param.nTrkCut[1] ) )
975 {
976 m_cut1++;
977 return -1;
978 }
979 // if(ntrk > 2) return -1;
980 for ( i = 0; i < ntrk; i++ )
981 {
982 m_nTrkAfTrkCut++;
983 fgTrk = true;
984 rectrk = event->getRecTrk( i );
985 nhitRec = rectrk->getNHits();
986 m_hnhitsRec->Fill( nhitRec );
987
988 for ( lay = 0; lay < MdcCalNLayer; lay++ ) { fgHitLay[lay] = false; }
989
990 nhitRecInn = 0;
991 nhitRecStp = 0;
992 nhitRecOut = 0;
993 for ( k = 0; k < nhitRec; k++ )
994 {
995 rechit = rectrk->getRecHit( k );
996 lay = rechit->getLayid();
997 doca = rechit->getDocaExc();
998 resiExc = rechit->getResiExc();
999 fgHitLay[lay] = true;
1000
1001 if ( lay < 8 ) nhitRecInn++;
1002 else if ( lay < 20 ) nhitRecStp++;
1003 else nhitRecOut++;
1004 }
1005 m_hnhitsRecInn->Fill( nhitRecInn );
1006 m_hnhitsRecStp->Fill( nhitRecStp );
1007 m_hnhitsRecOut->Fill( nhitRecOut );
1008
1009 // get momentum
1010 pt = rectrk->getPt();
1011 p = rectrk->getP();
1012
1013 // boost P to the cms
1014 p4 = rectrk->getP4();
1015 HepLorentzVector psip( xboost, yboost, zboost, ecm );
1016 Hep3Vector boostv = psip.boostVector();
1017 p4.boost( -boostv );
1018 p_cms = p4.rho();
1019 phicms = p4.phi();
1020 if ( phicms < 0 ) phicms += PI2;
1021 thetacms = p4.theta();
1022 costheCMS = cos( thetacms );
1023 if ( pt < 0 ) p_cms *= -1.0;
1024 p4.boost( boostv );
1025 // cout << setw(15) << p << setw(15) << p_cms << setw(15) << costheCMS << endl;
1026
1027 // cos(theta) cut
1028 if ( ( costheCMS < m_param.costheCut[0] ) || ( costheCMS > m_param.costheCut[1] ) )
1029 {
1030 m_cut2++;
1031 continue;
1032 }
1033
1034 // dr cut
1035 dr = rectrk->getDr();
1036 if ( fabs( dr ) > m_param.drCut )
1037 {
1038 m_cut3++;
1039 continue;
1040 }
1041
1042 // dz cut
1043 dz = rectrk->getDz();
1044 if ( fabs( dz ) > m_param.dzCut )
1045 {
1046 m_cut4++;
1047 continue;
1048 }
1049
1050 // momentum cut
1051 if ( ( fabs( p ) < m_param.pCut[0] ) || ( fabs( p ) > m_param.pCut[1] ) )
1052 {
1053 m_cut5++;
1054 continue;
1055 }
1056
1057 // charge cut, to check +/- asymmetry
1058 int charge = 1;
1059 if ( p < 0 ) charge = -1;
1060 if ( ( abs( m_param.charge ) == 1 ) && ( charge != m_param.charge ) ) { continue; }
1061
1062 // if(! fgTrk) continue;
1063
1064 // hit layer cut
1065 nhitlay = 0;
1066 for ( lay = 0; lay < MdcCalNLayer; lay++ )
1067 {
1068 if ( fgHitLay[lay] ) nhitlay++;
1069 }
1070 if ( nhitlay < m_param.nHitLayCut )
1071 {
1072 m_cut6++;
1073 continue;
1074 }
1075
1076 // nhit cut
1077 if ( nhitRec < m_param.nHitCut )
1078 {
1079 m_cut7++;
1080 continue;
1081 }
1082
1083 // bool fgNoise = rectrk->getFgNoiseRatio();
1084 // if(m_param.noiseCut && (!fgNoise)) continue;
1085 // cout << setw(10) << iRun << setw(15) << iEvt << setw(5) << fgNoise << endl;
1086
1087 // if(! ((fgHitLay[0]||fgHitLay[1]) && (fgHitLay[41]||fgHitLay[42])) ){
1088 // continue;
1089 // }
1090
1091 // calculate cell on the track
1092 int cellTrkPass[43];
1093 bool fgPass = getCellTrkPass( event, i, cellTrkPass );
1094 for ( lay = 0; lay < m_nlayer; lay++ )
1095 {
1096 fgAdd[lay] = 0;
1097 // if((16==lay) || (18==lay) || (19==lay) || (41==lay)){ // hv2200 2009-3
1098 if ( ( 15 == lay ) || ( 16 == lay ) || ( 18 == lay ) || ( 19 == lay ) || ( 40 == lay ) ||
1099 ( 41 == lay ) || ( 42 == lay ) )
1100 {
1101 int iCell = cellTrkPass[lay];
1102 if ( fgPass && ( iCell >= 0 ) && m_fgGoodWire[lay][iCell] ) m_effNtrk->Fill( lay );
1103 else fgAdd[lay] = 1;
1104 }
1105 else { m_effNtrk->Fill( lay ); }
1106 }
1107
1108 m_nTrkCal++;
1109 chisq = rectrk->getChisq();
1110 m_hchisq->Fill( chisq );
1111
1112 if ( pt > 0 )
1113 {
1114 m_hpt->Fill( pt );
1115 m_hptPos->Fill( pt );
1116 m_hp->Fill( p );
1117 m_hp_cms->Fill( p_cms );
1118 m_hpPos->Fill( p );
1119 m_hpPoscms->Fill( p_cms );
1120 }
1121 else
1122 {
1123 m_hpt->Fill( -1.0 * pt );
1124 m_hptNeg->Fill( -1.0 * pt );
1125 m_hp->Fill( -1.0 * p );
1126 m_hp_cms->Fill( -1.0 * p_cms );
1127 m_hpNeg->Fill( -1.0 * p );
1128 m_hpNegcms->Fill( -1.0 * p_cms );
1129 }
1130 if ( 2 == ntrk )
1131 {
1132 pTrk[i] = fabs( p );
1133 pTrkcms[i] = fabs( p_cms );
1134 }
1135
1136 dr = rectrk->getDr();
1137 phi0 = rectrk->getPhi0();
1138 kap = rectrk->getKappa();
1139 dz = rectrk->getDz();
1140 tanl = rectrk->getTanLamda();
1141 lamda = atan( tanl );
1142 theta = HFPI - lamda;
1143
1144 m_hdr->Fill( dr );
1145 m_hphi0->Fill( phi0 );
1146 m_hkap->Fill( kap );
1147 m_hdz->Fill( dz );
1148 m_htanl->Fill( tanl );
1149 m_hcosthe->Fill( cos( theta ) );
1150 if ( pt > 0 ) m_hcosthePos->Fill( cos( theta ) );
1151 else m_hcostheNeg->Fill( cos( theta ) );
1152
1153 philab = phi0 + HFPI;
1154 if ( philab > PI2 ) philab -= PI2;
1155 // cout << setw(15) << phi0 << setw(15) << philab << setw(15) << phicms << endl;
1156
1157 phiBin = (int)( philab / m_phiWid );
1158 phiBinCms = (int)( phicms / m_phiWid );
1159 theBin = (int)( ( cos( theta ) + 1.0 ) / m_theWid );
1160 theBinCms = (int)( ( cos( thetacms ) + 1.0 ) / m_theWid );
1161 if ( phiBin < 0 ) phiBin = 0;
1162 if ( phiBin >= NPhiBin ) phiBin = NPhiBin - 1;
1163 if ( phiBinCms < 0 ) phiBinCms = 0;
1164 if ( phiBinCms >= NPhiBin ) phiBinCms = NPhiBin - 1;
1165 if ( theBin < 0 ) theBin = 0;
1166 if ( theBin >= NThetaBin ) theBin = NThetaBin - 1;
1167 if ( theBinCms < 0 ) theBinCms = 0;
1168 if ( theBinCms >= NThetaBin ) theBinCms = NThetaBin - 1;
1169
1170 if ( pt > 0 )
1171 {
1172 m_ppPhi[phiBin]->Fill( p );
1173 m_ppPhiCms[phiBinCms]->Fill( p_cms );
1174 m_ppThe[theBin]->Fill( p );
1175 m_ppTheCms[theBinCms]->Fill( p_cms );
1176 m_ppThePhi[theBin][phiBin]->Fill( p );
1177 m_ppThePhiCms[theBinCms][phiBinCms]->Fill( p_cms );
1178 }
1179 else
1180 {
1181 m_pnPhi[phiBin]->Fill( -1.0 * p );
1182 m_pnPhiCms[phiBinCms]->Fill( -1.0 * p_cms );
1183 m_pnThe[theBin]->Fill( -1.0 * p );
1184 m_pnTheCms[theBinCms]->Fill( -1.0 * p_cms );
1185 m_pnThePhi[theBin][phiBin]->Fill( -1.0 * p );
1186 m_pnThePhiCms[theBinCms][phiBinCms]->Fill( -1.0 * p_cms );
1187 }
1188
1189 x0 = dr * cos( phi0 );
1190 y0 = dr * sin( phi0 );
1191 m_hx0->Fill( x0 );
1192 m_hy0->Fill( y0 );
1193 if ( m_nGrPoint < 10000 )
1194 {
1195 m_grX0Y0->SetPoint( m_nGrPoint, x0, y0 );
1196 m_nGrPoint++;
1197 }
1198
1199 if ( kap < 0 )
1200 {
1201 zminus = dz;
1202 pm.push_back( p4 );
1203 phibinm = phiBinCms;
1204 }
1205 else
1206 {
1207 zplus = dz;
1208 pp.push_back( p4 );
1209 phibinp = phiBinCms;
1210 }
1211
1212 // cout << "phi = " << setw(15) << philab << setw(15) << philab*180./3.14159 << setw(15)
1213 // << p << endl;
1214 ntrkCal++;
1215 trkFlag[i] = true;
1216 nhitCal = 0;
1217 nhitCalInn = 0;
1218 nhitCalStp = 0;
1219 nhitCalOut = 0;
1220 for ( k = 0; k < nhitRec; k++ )
1221 {
1222 rechit = rectrk->getRecHit( k );
1223
1224 lay = rechit->getLayid();
1225 cel = rechit->getCellid();
1226 lr = rechit->getLR();
1227 stat = rechit->getStat();
1228 doca = rechit->getDocaExc();
1229 resiInc = rechit->getResiIncLR();
1230 resiExc = rechit->getResiExcLR();
1231 entr = rechit->getEntra();
1232 tdr = rechit->getTdrift();
1233 traw = ( rechit->getTdc() ) * MdcCalTdcCnv;
1234 wir = m_mdcGeomSvc->Wire( lay, cel )->Id();
1235
1236 m_cel[lay] = (long)cel;
1237 m_lr[lay] = (long)lr;
1238 m_run[lay] = iRun;
1239 m_evt[lay] = iEvt;
1240 m_doca[lay] = doca;
1241 m_dm[lay] = rechit->getDmeas();
1242 m_tdr[lay] = tdr;
1243 m_tdc[lay] = traw;
1244 m_entr[lay] = entr * 180.0 / 3.14;
1245 m_zhit[lay] = rechit->getZhit();
1246 m_qhit[lay] = rechit->getQhit();
1247 m_p[lay] = p_cms;
1248 m_pt[lay] = pt;
1249 m_phi0[lay] = phi0;
1250 m_tanl[lay] = tanl;
1251 qhit = rechit->getQhit();
1252 // cout << "layer " << setw(5) << lay << setw(5) << cel << setw(5) << lr <<
1253 // endl; calculating hitphi
1254 xx =
1255 ( m_zhit[lay] - m_zw[wir] ) * ( m_xe[wir] - m_xw[wir] ) / ( m_ze[wir] - m_zw[wir] ) +
1256 m_xw[wir];
1257 yy =
1258 ( m_zhit[lay] - m_zw[wir] ) * ( m_ye[wir] - m_yw[wir] ) / ( m_ze[wir] - m_zw[wir] ) +
1259 m_yw[wir];
1260 rr = sqrt( ( xx * xx ) + ( yy * yy ) );
1261 dphi = fabs( doca ) / m_radii[lay];
1262
1263 if ( yy >= 0 ) wphi = acos( xx / rr );
1264 else wphi = PI2 - acos( xx / rr );
1265 if ( 1 == lr ) hitphi = wphi + dphi; // mention
1266 else hitphi = wphi - dphi;
1267 if ( hitphi < 0 ) hitphi += PI2;
1268 else if ( hitphi > PI2 ) hitphi -= PI2;
1269
1270 m_hitphi[lay] = hitphi;
1271
1272 if ( ( fabs( doca ) > m_docaMax[lay] ) || ( fabs( resiExc ) > m_param.resiCut[lay] ) )
1273 { continue; }
1274
1275 if ( m_param.fgAdjacLayerCut )
1276 {
1277 if ( 0 == lay )
1278 {
1279 if ( !fgHitLay[1] ) continue;
1280 }
1281 else if ( 42 == lay )
1282 {
1283 if ( !fgHitLay[41] ) continue;
1284 }
1285 else
1286 {
1287 if ( ( !fgHitLay[lay - 1] ) && ( !fgHitLay[lay + 1] ) ) continue;
1288
1289 // for boundary layers
1290 if ( m_param.fgBoundLayerCut && m_layBound[lay] &&
1291 ( ( !fgHitLay[lay - 1] ) || ( !fgHitLay[lay + 1] ) ) )
1292 continue;
1293 }
1294 }
1295
1296 if ( ( 1 == m_param.hitStatCut ) && ( 1 != stat ) ) continue;
1297
1298 // fill xtplot tree
1299 if ( ( 1 == m_param.fillNtuple ) && ( m_nEvtNtuple < m_param.nEvtNtuple ) )
1300 { m_xtTuple[lay]->write(); }
1301
1302 if ( 1 == m_param.hitStatCut )
1303 {
1304 if ( ( 0 == fgAdd[lay] ) && ( 1 == stat ) )
1305 {
1306 m_effNtrkRecHit->Fill( lay );
1307 fgAdd[lay] = 1;
1308 }
1309 }
1310 else
1311 {
1312 if ( 0 == fgAdd[lay] )
1313 {
1314 m_effNtrkRecHit->Fill( lay );
1315 fgAdd[lay] = 1;
1316 }
1317 }
1318
1319 nhitCal++;
1320 if ( lay < 8 ) nhitCalInn++;
1321 else if ( lay < 20 ) nhitCalStp++;
1322 else nhitCalOut++;
1323
1324 m_wirehitmap->Fill( wir );
1325 m_layerhitmap->Fill( lay );
1326
1327 m_htraw[lay]->Fill( traw );
1328 m_htdr[lay]->Fill( tdr );
1329 m_htdrlr[lay][lr]->Fill( tdr );
1330 m_hadc[lay]->Fill( qhit );
1331
1332 m_hresAllInc->Fill( resiInc );
1333 m_hresAllExc->Fill( resiExc );
1334 double resiAve = 0.5 * ( resiInc + resiExc );
1335 m_hresAllAve->Fill( resiAve );
1336 if ( yy > 0 ) m_hresAllExcUp->Fill( resiExc );
1337 else m_hresAllExcDw->Fill( resiExc );
1338
1339 if ( lay < 8 )
1340 {
1341 m_hresInnInc->Fill( resiInc );
1342 m_hresInnExc->Fill( resiExc );
1343 if ( yy > 0 ) m_hresInnExcUp->Fill( resiExc );
1344 else m_hresInnExcDw->Fill( resiExc );
1345 }
1346 else if ( lay < 20 )
1347 {
1348 m_hresStpInc->Fill( resiInc );
1349 m_hresStpExc->Fill( resiExc );
1350 if ( yy > 0 ) m_hresStpExcUp->Fill( resiExc );
1351 else m_hresStpExcDw->Fill( resiExc );
1352 }
1353 else
1354 {
1355 m_hresOutInc->Fill( resiInc );
1356 m_hresOutExc->Fill( resiExc );
1357 if ( yy > 0 ) m_hresOutExcUp->Fill( resiExc );
1358 else m_hresOutExcDw->Fill( resiExc );
1359 }
1360
1361 int qbin = (int)( ( qhit - 100.0 ) / 100.0 );
1362 if ( qbin >= 0 && qbin < 14 )
1363 {
1364 m_hresAveAllQ[qbin]->Fill( resiAve );
1365 m_hresAveLayQ[lay][qbin]->Fill( resiAve );
1366 if ( lay > 7 ) m_hresAveOutQ[qbin]->Fill( resiAve );
1367 }
1368
1369 if ( ( !m_param.cosmicDwTrk ) || ( m_param.cosmicDwTrk && ( yy < 0 ) ) )
1370 {
1371 m_hresInc[lay]->Fill( resiInc );
1372 m_hreslrInc[lay][lr]->Fill( resiInc );
1373 m_hresExc[lay]->Fill( resiExc );
1374 m_hreslrExc[lay][lr]->Fill( resiExc );
1375 m_hresAve[lay]->Fill( resiAve );
1376 m_hreslrAve[lay][lr]->Fill( resiAve );
1377 }
1378
1379 int iPhi = (int)( hitphi * 20.0 / PI2 );
1380 if ( iPhi >= 20 ) iPhi = 19;
1381 m_hresphi[lay][iPhi]->Fill( ( resiInc + resiExc ) * 0.5 );
1382
1383 // bin = (int)(fabs(doca) / m_dwid);
1384 bin = (int)( fabs( rechit->getDmeas() ) / m_dwid );
1385 iEntr = m_mdcFunSvc->getSdEntrIndex( entr );
1386 if ( 1 == m_nEntr[lay] ) { iEntr = 0; }
1387 else if ( 2 == m_nEntr[lay] )
1388 {
1389 if ( entr > 0.0 ) iEntr = 1;
1390 else iEntr = 0;
1391 }
1392 else if ( 3 == m_nEntr[lay] )
1393 {
1394 double entrBinWid = 0.174532925; // 10 degree
1395 iEntr = (int)( fabs( entr ) / entrBinWid );
1396 if ( iEntr > 2 ) iEntr = 2;
1397 }
1398 if ( ( iEntr < MdcCalNENTRSD ) && ( bin < MdcCalSdNBIN ) )
1399 {
1400 key = getHresId( lay, iEntr, lr, bin );
1401 if ( 1 == m_mapr2d.count( key ) )
1402 {
1403 hid = m_mapr2d[key];
1404 m_hr2dInc[hid]->Fill( resiInc );
1405 m_hr2dExc[hid]->Fill( resiExc );
1406 }
1407 }
1408
1409 if ( ( tdr > 0 ) && ( tdr < 750 ) )
1410 {
1411 if ( tdr < 300 ) bin = (int)( tdr / 10.0 );
1412 else bin = (int)( ( tdr - 300.0 ) / 30.0 ) + 29;
1413 m_hr2t[lay][iEntr][lr][bin]->Fill( resiExc );
1414 }
1415 } // loop of nhits
1416 m_nEvtNtuple++;
1417 m_hnhitsCal->Fill( nhitCal );
1418 m_hnhitsCalInn->Fill( nhitCalInn );
1419 m_hnhitsCalStp->Fill( nhitCalStp );
1420 m_hnhitsCalOut->Fill( nhitCalOut );
1421 } // end of track loop
1422 m_hnTrkCal->Fill( ntrkCal );
1423 if ( 2 == ntrkCal )
1424 {
1425 if ( pTrk[0] > pTrk[1] ) m_hpMax->Fill( pTrk[0] );
1426 else m_hpMax->Fill( pTrk[1] );
1427
1428 if ( pTrkcms[0] > pTrkcms[1] ) m_hpMaxCms->Fill( pTrkcms[0] );
1429 else m_hpMaxCms->Fill( pTrkcms[1] );
1430 }
1431 if ( ntrkCal > 0 ) m_hTesCalUse->Fill( tes );
1432
1433 double delZ0;
1434 if ( ( fabs( zminus ) < 9000.0 ) && ( fabs( zplus ) < 9000.0 ) ) delZ0 = zplus - zminus;
1435 m_hdelZ0->Fill( delZ0 );
1436
1437 if ( 1 == pp.size() * pm.size() )
1438 {
1439 HepLorentzVector ptot = pp[0] + pm[0];
1440 bool fourmomcut = false;
1441 // fourmomcut = (ptot.x()>0.02 && ptot.x()<0.06) && (fabs(ptot.y()) < 0.02)
1442 // && (ptot.z()>-0.01 && ptot.z()<0.03) && (ptot.e()>3.4 && ptot.e()<4.0);
1443 fourmomcut = ( fabs( ptot.x() - 0.04 ) < 0.026 ) && ( fabs( ptot.y() ) < 0.026 ) &&
1444 ( fabs( ptot.z() - 0.005 ) < 0.036 ) && ( fabs( ptot.e() - ecm ) < 0.058 );
1445 // cout << "x = " << ptot.x() << ", y = " << ptot.y() << ", z = " << ptot.z() << ", e = "
1446 // << ptot.e() << endl;
1447 if ( fourmomcut )
1448 {
1449 HepLorentzVector psip( xboost, yboost, zboost, ecm );
1450 Hep3Vector boostv = psip.boostVector();
1451 pp[0].boost( -boostv );
1452 pm[0].boost( -boostv );
1453 m_hp_cut->Fill( pp[0].rho() );
1454 m_hp_cut->Fill( pm[0].rho() );
1455 }
1456 }
1457
1458 if ( 2 == ntrk )
1459 for ( i = 0; i < ntrk; i++ ) pTrk[i] = ( event->getRecTrk( i ) )->getP();
1460 if ( ( 5 == m_param.particle ) && ( 2 == ntrk ) && ( fabs( pTrk[0] ) < 5 ) &&
1461 ( fabs( pTrk[1] ) < 5 ) )
1462 {
1463 // if(1==ntrk) p = (event->getRecTrk(0)) -> getP();
1464 // if((5==m_param.particle) && (1==ntrk) && (fabs(p)<5)){
1465 m_tescos = tes;
1466 m_tesFlagcos = esTimeflag;
1467 for ( i = 0; i < ntrk; i++ )
1468 {
1469 rectrk = event->getRecTrk( i );
1470 phi0 = rectrk->getPhi0();
1471 phi0 = ( ( phi0 + HFPI ) > PI2 ) ? ( phi0 + HFPI - PI2 ) : ( phi0 + HFPI );
1472
1473 tanl = rectrk->getTanLamda();
1474 lamda = atan( tanl );
1475 theta = HFPI - lamda;
1476
1477 if ( phi0 < ( 2.0 * HFPI ) )
1478 {
1479 m_nhitUpcos = rectrk->getNHits();
1480 m_pUpcos = rectrk->getP();
1481 m_ptUpcos = rectrk->getPt();
1482 m_phiUpcos = phi0;
1483 m_drUpcos = rectrk->getDr();
1484 m_dzUpcos = rectrk->getDz();
1485 m_ctheUpcos = cos( theta );
1486 }
1487 else
1488 {
1489 m_nhitDwcos = rectrk->getNHits();
1490 m_pDwcos = rectrk->getP();
1491 m_ptDwcos = rectrk->getPt();
1492 m_phiDwcos = phi0;
1493 m_drDwcos = rectrk->getDr();
1494 m_dzDwcos = rectrk->getDz();
1495 m_ctheDwcos = cos( theta );
1496
1497 if ( m_pDwcos > 0 ) m_chargecos = 1;
1498 else m_chargecos = 0;
1499 }
1500 }
1501 m_cosmic->write();
1502 }
1503
1504 if ( 1 == m_param.fgCalDetEffi ) calDetEffi();
1505
1506 return 1;
1507}
1508
1509void MdcCalib::printCut() const {
1510 cout << "Events " << m_hTesAll->GetEntries() << ", Tes cut " << m_hTesCalFlag->GetEntries()
1511 << ", nTrkCut " << m_cut1 << endl;
1512 cout << "Tot Tracks " << m_nTrkAfTrkCut << ", cos(theta)_cut " << m_cut2 << ", drCut "
1513 << m_cut3 << ", dzCut " << m_cut4 << ", momCut" << m_cut5 << ", nHitLayer_cut "
1514 << m_cut6 << ", nHit_cut " << m_cut7 << endl;
1515 cout << "nTrack after all cuts: " << m_nTrkCal << endl;
1516}
1517
1518int MdcCalib::updateConst( MdcCalibConst* calconst ) {
1519 IMessageSvc* msgSvc;
1520 Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
1521 MsgStream log( msgSvc, "MdcCalib" );
1522 log << MSG::DEBUG << "MdcCalib::updateConst()" << endmsg;
1523
1524 int lay;
1525 double effi;
1526 double effErr;
1527
1528 int nGoodAll = 0;
1529 int nGoodInn = 0;
1530 int nGoodStp = 0;
1531 int nGoodOut = 0;
1532 int nTotAll = 0;
1533 int nTotInn = 0;
1534 int nTotStp = 0;
1535 int nTotOut = 0;
1536 ofstream feffi( "MdcLayerEffi.dat" );
1537 for ( lay = 0; lay < m_nlayer; lay++ )
1538 {
1539 double effNtrk = m_effNtrk->GetBinContent( lay + 1 );
1540 double effGoodHit = m_effNtrkRecHit->GetBinContent( lay + 1 );
1541 nGoodAll += effGoodHit;
1542 if ( lay < 8 ) nGoodInn += effGoodHit;
1543 else if ( lay < 20 ) nGoodStp += effGoodHit;
1544 else nGoodOut += effGoodHit;
1545
1546 nTotAll += effNtrk;
1547 if ( lay < 8 ) nTotInn += effNtrk;
1548 else if ( lay < 20 ) nTotStp += effNtrk;
1549 else nTotOut += effNtrk;
1550
1551 effi = (double)effGoodHit / (double)effNtrk;
1552 effErr = sqrt( effi * ( 1 - effi ) / (double)effNtrk );
1553 feffi << setw( 5 ) << lay << setw( 15 ) << effi << setw( 15 ) << effErr << setw( 15 )
1554 << effGoodHit << setw( 15 ) << effNtrk << endl;
1555 }
1556 double effiAll = (double)nGoodAll / (double)( nTotAll );
1557 double errAll = sqrt( effiAll * ( 1 - effiAll ) / (double)( nTotAll ) );
1558 double effiInn = (double)nGoodInn / (double)( nTotInn );
1559 double errInn = sqrt( effiInn * ( 1 - effiInn ) / (double)( nTotInn ) );
1560 double effiStp = (double)nGoodStp / (double)( nTotStp );
1561 double errStp = sqrt( effiStp * ( 1 - effiStp ) / (double)( nTotStp ) );
1562 double effiOut = (double)nGoodOut / (double)( nTotOut );
1563 double errOut = sqrt( effiOut * ( 1 - effiOut ) / (double)( nTotOut ) );
1564 feffi << endl
1565 << "EffiAll: " << setw( 15 ) << effiAll << setw( 15 ) << errAll << setw( 15 )
1566 << nGoodAll << setw( 15 ) << nTotAll << endl;
1567 feffi << endl
1568 << "EffiInn: " << setw( 15 ) << effiInn << setw( 15 ) << errInn << setw( 15 )
1569 << nGoodInn << setw( 15 ) << nTotInn << endl;
1570 feffi << endl
1571 << "EffiStp: " << setw( 15 ) << effiStp << setw( 15 ) << errStp << setw( 15 )
1572 << nGoodStp << setw( 15 ) << nTotStp << endl;
1573 feffi << endl
1574 << "EffiOut: " << setw( 15 ) << effiOut << setw( 15 ) << errOut << setw( 15 )
1575 << nGoodOut << setw( 15 ) << nTotOut << endl;
1576 feffi.close();
1577
1578 // calculate efficiency without the impact of track fitting
1579 if ( 0 != m_param.fgCalDetEffi )
1580 {
1581 int nHitAll[] = { 0, 0 };
1582 int nHitInn[] = { 0, 0 };
1583 int nHitStp[] = { 0, 0 };
1584 int nHitOut[] = { 0, 0 };
1585 ofstream feff2( "MdcHitEffi.dat" );
1586 for ( lay = 0; lay < m_nlayer; lay++ )
1587 {
1588 nHitAll[0] += m_hitNum[lay][0];
1589 nHitAll[1] += m_hitNum[lay][1];
1590 if ( lay < 8 )
1591 {
1592 nHitInn[0] += m_hitNum[lay][0];
1593 nHitInn[1] += m_hitNum[lay][1];
1594 }
1595 else if ( lay < 20 )
1596 {
1597 nHitStp[0] += m_hitNum[lay][0];
1598 nHitStp[1] += m_hitNum[lay][1];
1599 }
1600 else
1601 {
1602 nHitOut[0] += m_hitNum[lay][0];
1603 nHitOut[1] += m_hitNum[lay][1];
1604 }
1605
1606 effi = (double)m_hitNum[lay][1] / (double)m_hitNum[lay][0];
1607 effErr = sqrt( effi * ( 1 - effi ) / (double)m_hitNum[lay][0] );
1608 feff2 << setw( 5 ) << lay << setw( 15 ) << effi << setw( 15 ) << effErr << setw( 15 )
1609 << m_hitNum[lay][1] << setw( 15 ) << m_hitNum[lay][0] << endl;
1610 }
1611 effiAll = (double)nHitAll[1] / (double)nHitAll[0];
1612 errAll = sqrt( effiAll * ( 1 - effiAll ) ) / (double)nHitAll[0];
1613 effiInn = (double)nHitInn[1] / (double)nHitInn[0];
1614 errInn = sqrt( effiInn * ( 1 - effiInn ) ) / (double)nHitInn[0];
1615 effiStp = (double)nHitStp[1] / (double)nHitStp[0];
1616 errStp = sqrt( effiStp * ( 1 - effiStp ) ) / (double)nHitStp[0];
1617 effiOut = (double)nHitOut[1] / (double)nHitOut[0];
1618 errOut = sqrt( effiOut * ( 1 - effiOut ) ) / (double)nHitOut[0];
1619 feff2 << endl
1620 << "EffiAll: " << setw( 15 ) << effiAll << setw( 15 ) << errAll << setw( 15 )
1621 << nHitAll[1] << setw( 15 ) << nHitAll[0] << endl;
1622 feff2 << endl
1623 << "EffiInn: " << setw( 15 ) << effiInn << setw( 15 ) << errInn << setw( 15 )
1624 << nHitInn[1] << setw( 15 ) << nHitInn[0] << endl;
1625 feff2 << endl
1626 << "EffiStp: " << setw( 15 ) << effiStp << setw( 15 ) << errStp << setw( 15 )
1627 << nHitStp[1] << setw( 15 ) << nHitStp[0] << endl;
1628 feff2 << endl
1629 << "EffiOut: " << setw( 15 ) << effiOut << setw( 15 ) << errOut << setw( 15 )
1630 << nHitOut[1] << setw( 15 ) << nHitOut[0] << endl;
1631 feff2.close();
1632 }
1633
1634 // get resolution
1635 int i;
1636 int iEntr;
1637 int lr;
1638 int bin;
1639 int key;
1640 int hid;
1641
1642 Stat_t entry;
1643 double sigm[MdcCalSdNBIN];
1644 if ( m_param.calSigma )
1645 {
1646 ofstream fr2d( "logr2d.dat" );
1647 for ( lay = 0; lay < m_nlayer; lay++ )
1648 {
1649 for ( iEntr = 0; iEntr < m_nEntr[lay]; iEntr++ )
1650 {
1651 for ( lr = 0; lr < 2; lr++ )
1652 {
1653 fr2d << setw( 3 ) << lay << setw( 3 ) << iEntr << setw( 3 ) << lr << endl;
1654 for ( bin = 0; bin < m_nBin[lay]; bin++ )
1655 {
1656 key = getHresId( lay, iEntr, lr, bin );
1657 hid = m_mapr2d[key];
1658
1659 if ( 1 == m_param.resiType )
1660 {
1661 entry = m_hr2dExc[hid]->GetEntries();
1662 if ( entry > 500 )
1663 {
1664 m_hr2dExc[hid]->Fit( "gaus", "Q" );
1665 sigm[bin] = m_hr2dExc[hid]->GetFunction( "gaus" )->GetParameter( 2 );
1666 }
1667 else if ( entry > 100 ) { sigm[bin] = m_hr2dExc[hid]->GetRMS(); }
1668 else { sigm[bin] = 0.2; }
1669 }
1670 else
1671 {
1672 entry = m_hr2dInc[hid]->GetEntries();
1673 if ( entry > 500 )
1674 {
1675 m_hr2dInc[hid]->Fit( "gaus", "Q" );
1676 sigm[bin] = m_hr2dInc[hid]->GetFunction( "gaus" )->GetParameter( 2 );
1677 }
1678 else if ( entry > 100 ) { sigm[bin] = m_hr2dInc[hid]->GetRMS(); }
1679 else { sigm[bin] = 0.2; }
1680 }
1681 if ( sigm[bin] < 0.05 ) sigm[bin] = 0.05; // for boundary layers
1682 } // end of bin loop
1683
1684 for ( bin = m_nBin[lay]; bin < MdcCalSdNBIN; bin++ )
1685 { sigm[bin] = sigm[m_nBin[lay] - 1]; }
1686
1687 for ( bin = 0; bin < MdcCalSdNBIN; bin++ )
1688 {
1689 if ( 1 == m_param.fgCalib[lay] )
1690 {
1691 // calconst -> resetSdpar(lay, iEntr, lr, bin,
1692 // sigm[bin]);
1693 if ( 1 == m_nEntr[lay] )
1694 {
1695 for ( i = 0; i < 6; i++ ) calconst->resetSdpar( lay, i, lr, bin, sigm[bin] );
1696 }
1697 else if ( 2 == m_nEntr[lay] )
1698 {
1699 if ( 0 == iEntr )
1700 {
1701 for ( i = 0; i < 3; i++ )
1702 { // entr<0
1703 calconst->resetSdpar( lay, i, lr, bin, sigm[bin] );
1704 }
1705 }
1706 else
1707 {
1708 for ( i = 3; i < 6; i++ )
1709 { // entr>0
1710 calconst->resetSdpar( lay, i, lr, bin, sigm[bin] );
1711 }
1712 }
1713 }
1714 }
1715 else { sigm[bin] = calconst->getSdpar( lay, iEntr, lr, bin ); }
1716 fr2d << setw( 5 ) << bin << setw( 15 ) << sigm[bin] << endl;
1717 } // end of bin loop
1718 }
1719 } // end of entr loop
1720 }
1721 fr2d.close();
1722 }
1723
1724 return 1;
1725}
1726
1727int MdcCalib::getHresId( int lay, int entr, int lr, int bin ) const {
1728 int index = ( ( lay << HRESLAYER_INDEX ) & HRESLAYER_MASK ) |
1729 ( ( entr << HRESENTRA_INDEX ) & HRESENTRA_MASK ) |
1730 ( ( lr << HRESLR_INDEX ) & HRESLR_MASK ) |
1731 ( ( bin << HRESBIN_INDEX ) & HRESBIN_MASK );
1732 return index;
1733}
1734
1735int MdcCalib::calDetEffi() {
1736 IMessageSvc* msgSvc;
1737 Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
1738 MsgStream log( msgSvc, "MdcCalib" );
1739
1740 IDataProviderSvc* eventSvc = NULL;
1741 Gaudi::svcLocator()->service( "EventDataSvc", eventSvc );
1742 SmartDataPtr<MdcDigiCol> mdcDigiCol( eventSvc, "/Event/Digi/MdcDigiCol" );
1743 if ( !mdcDigiCol ) { log << MSG::FATAL << "Could not find event" << endmsg; }
1744
1745 int lay;
1746 int cel;
1747 bool hitCel[43][288];
1748 int hitCel2[43][288];
1749 for ( lay = 0; lay < 43; lay++ )
1750 {
1751 for ( cel = 0; cel < 288; cel++ )
1752 {
1753 hitCel[lay][cel] = false;
1754 hitCel2[lay][cel] = 0;
1755 }
1756 }
1757
1758 MdcDigiCol::iterator iter = mdcDigiCol->begin();
1759 unsigned fgOverFlow;
1760 int digiId = 0;
1761 Identifier id;
1762 for ( ; iter != mdcDigiCol->end(); iter++, digiId++ )
1763 {
1764 MdcDigi* aDigi = ( *iter );
1765 id = ( aDigi )->identify();
1766
1767 lay = MdcID::layer( id );
1768 cel = MdcID::wire( id );
1769 fgOverFlow = ( aDigi )->getOverflow();
1770
1771 // if ( !((fgOverFlow == 0)||(fgOverFlow ==12)) ||
1772 // (aDigi->getTimeChannel() == 0x7FFFFFFF) ||
1773 // (aDigi->getChargeChannel() == 0x7FFFFFFF) ) continue;
1774 if ( ( ( fgOverFlow & 3 ) != 0 ) || ( ( fgOverFlow & 12 ) != 0 ) ||
1775 ( aDigi->getTimeChannel() == 0x7FFFFFFF ) ||
1776 ( aDigi->getChargeChannel() == 0x7FFFFFFF ) )
1777 continue;
1778 hitCel[lay][cel] = true;
1779 hitCel2[lay][cel] = 1;
1780 }
1781
1782 SmartDataPtr<RecMdcTrackCol> newtrkCol( eventSvc, "/Event/Recon/RecMdcTrackCol" );
1783 if ( !newtrkCol )
1784 {
1785 log << MSG::ERROR << "Could not find RecMdcTrackCol" << endmsg;
1786 return -1;
1787 }
1788
1789 double dphi = 1.0;
1790 Identifier identifier;
1791 MdcID mdcid;
1792 RecMdcTrackCol::iterator it_trk = newtrkCol->begin();
1793 for ( it_trk = newtrkCol->begin(); it_trk != newtrkCol->end(); it_trk++ )
1794 {
1795 HitRefVec gothits = ( *it_trk )->getVecHits();
1796 HitRefVec::iterator it_hit = gothits.begin();
1797 for ( ; it_hit != gothits.end(); it_hit++ )
1798 {
1799 identifier = ( *it_hit )->getMdcId();
1800 lay = mdcid.layer( identifier );
1801 cel = mdcid.wire( identifier );
1802 hitCel2[lay][cel] = 2;
1803 }
1804 }
1805 for ( it_trk = newtrkCol->begin(); it_trk != newtrkCol->end(); it_trk++ )
1806 {
1807 HepVector helix = ( *it_trk )->helix();
1808 HepSymMatrix helixErr = ( *it_trk )->err();
1809 double phi0 = ( *it_trk )->helix( 1 );
1810 double phiTrk = phi0 + HFPI;
1811 if ( phiTrk > PI2 ) phiTrk -= PI2;
1812
1813 for ( lay = 0; lay < 43; lay++ )
1814 {
1815 double docamin = 0.9; // cm
1816 if ( lay < 8 ) docamin = 0.6; // cm
1817 int celmin = -1;
1818 int ncel = m_mdcGeomSvc->Layer( lay )->NCell();
1819 for ( cel = 0; cel < ncel; cel++ )
1820 {
1821 double wphi;
1822 const MdcGeoWire* pWire = m_mdcGeomSvc->Wire( lay, cel );
1823 double xx = pWire->Backward().x();
1824 double yy = pWire->Backward().y();
1825 double rr = sqrt( ( xx * xx ) + ( yy * yy ) );
1826 if ( yy >= 0 ) wphi = acos( xx / rr );
1827 else wphi = CLHEP::twopi - acos( xx / rr );
1828
1829 if ( !( ( fabs( wphi - phiTrk ) < dphi ) || ( fabs( wphi + PI2 - phiTrk ) < dphi ) ||
1830 ( fabs( wphi - PI2 - phiTrk ) < dphi ) ) )
1831 { continue; }
1832
1833 double doca = m_mdcUtilitySvc->doca( lay, cel, helix, helixErr );
1834 // cout << setw(5) << lay << setw(5) << cel << setw(15) << doca <<
1835 // endl;
1836 if ( fabs( doca ) < fabs( docamin ) )
1837 {
1838 docamin = doca;
1839 celmin = cel;
1840 }
1841 }
1842
1843 if ( ( celmin > -1 ) && ( m_fgGoodWire[lay][celmin] ) )
1844 {
1845 int wir = m_mdcGeomSvc->Wire( lay, celmin )->Id();
1846 m_hitEffAll->Fill( wir );
1847 m_hitEffAll->Fill( 6799 );
1848 if ( lay < 8 ) m_hitEffAll->Fill( 6796 );
1849 else if ( lay < 20 ) m_hitEffAll->Fill( 6797 );
1850 else m_hitEffAll->Fill( 6798 );
1851 // double doca = m_mdcUtilitySvc->doca(lay, celmin, helix, helixErr,
1852 // true, true); // test doca cout << "662p01 " << setw(5) << lay << setw(5) << celmin
1853 // << setw(15) << doca << endl;
1854
1855 if ( hitCel[lay][celmin] )
1856 {
1857 m_hitEffRaw->Fill( wir );
1858 m_hitEffRaw->Fill( 6799 );
1859 if ( lay < 8 ) m_hitEffRaw->Fill( 6796 );
1860 else if ( lay < 20 ) m_hitEffRaw->Fill( 6797 );
1861 else m_hitEffRaw->Fill( 6798 );
1862 }
1863 if ( 2 == hitCel2[lay][celmin] )
1864 {
1865 m_hitEffRec->Fill( wir );
1866 m_hitEffRec->Fill( 6799 );
1867 if ( lay < 8 ) m_hitEffRec->Fill( 6796 );
1868 else if ( lay < 20 ) m_hitEffRec->Fill( 6797 );
1869 else m_hitEffRec->Fill( 6798 );
1870 }
1871 }
1872 }
1873 }
1874
1875 // int ncel;
1876 // int fgLayHit[43];
1877 // int celHit[43];
1878 // int celDocaMin;
1879 // double docaMin;
1880 // double trkpar[5];
1881
1882 // for(i=0; i<ntrk; i++){
1883 // if(!trkFlag[i]) continue;
1884 // rectrk = event -> getRecTrk(i);
1885 // nhitRec = rectrk -> getNHits();
1886
1887 // HepVector helix = rectrk->getHelix();
1888 // HepSymMatrix helixErr = rectrk->getHelixErr();
1889
1890 // for(lay=0; lay<43; lay++){
1891 // // m_hitNum[lay][0]++;
1892 // fgLayHit[lay] = false;
1893 // celHit[lay] = -1;
1894 // }
1895 // for(k=0; k<nhitRec; k++){
1896 // rechit = rectrk -> getRecHit(k);
1897 // lay = rechit -> getLayid();
1898 // cel = rechit -> getCellid();
1899 // fgLayHit[lay] = true;
1900 // celHit[lay] = cel;
1901 // hitCel2[lay][cel] = 2;
1902
1903 // // HepVector helix = rechit->getHelix();
1904 // // HepSymMatrix helixErr = rechit->getHelixErr();
1905 // double doca_rec = m_mdcUtilitySvc->doca(lay, cel, helix, helixErr);
1906 // cout << setw(15) << doca_rec*10. << setw(15) << rechit -> getDocaInc() << endl;
1907 // }
1908
1909 // // for(lay=0; lay<43; lay++){
1910 // // if(fgLayHit[lay]){
1911 // // m_hitNum[lay][0]++;
1912 // // m_hitNum[lay][1]++;
1913 // // } else{
1914 // // if(lay<8) docaMin = 8.0;
1915 // // else docaMin = 10.5;
1916 // // celDocaMin = -1;
1917
1918 // // ncel = m_mdcGeomSvc->Layer(lay)->NCell();
1919 // // for(cel=0; cel<ncel; cel++){
1920 // // double doca_rec = m_mdcUtilitySvc->doca(lay, cel, helix, helixErr);
1921 // // if(fabs(docaCal) < docaMin){
1922 // // docaMin = fabs(docaCal);
1923 // // celDocaMin = cel;
1924 // // }
1925 // // }
1926 // // if(celDocaMin > -1){
1927 // // m_hitNum[lay][0]++;
1928 // // if(hitCel[lay][celDocaMin]) m_hitNum[lay][1]++;
1929 // // }
1930 // // }
1931 // // }
1932 // }
1933
1934 // int nraw = 0;
1935 // int nrec = 0;
1936 // for(lay=0; lay<43; lay++){
1937 // for(cel=0; cel<288; cel++){
1938 // if(hitCel2[lay][cel] > 0) nraw++;
1939 // if(hitCel2[lay][cel] > 1) nrec++;
1940 // }
1941 // }
1942 // double ratio = (double)nrec / (double)nraw;
1943 // m_hratio->Fill(ratio);
1944 return 1;
1945}
1946
1947bool MdcCalib::getCellTrkPass( MdcCalEvent* event, int iTrk, int cellTrkPass[] ) {
1948 MdcCalRecTrk* rectrk = event->getRecTrk( iTrk );
1949 int nHits = rectrk->getNHits();
1950 int hitInTrk[100];
1951 for ( int k = 0; k < nHits; k++ )
1952 {
1953 MdcCalRecHit* rechit = rectrk->getRecHit( k );
1954 int lay = rechit->getLayid();
1955 int cel = rechit->getCellid();
1956 int wir = m_mdcGeomSvc->Wire( lay, cel )->Id();
1957 hitInTrk[k] = wir;
1958 }
1959
1960 IDataProviderSvc* eventSvc = NULL;
1961 Gaudi::svcLocator()->service( "EventDataSvc", eventSvc );
1962
1963 SmartDataPtr<RecMdcTrackCol> newtrkCol( eventSvc, "/Event/Recon/RecMdcTrackCol" );
1964 if ( !newtrkCol )
1965 {
1966 // log << MSG::ERROR << "Could not find RecMdcTrackCol" << endmsg;
1967 return false;
1968 }
1969 MdcID mdcid;
1970 Identifier identifier;
1971 double dphi = 1.0;
1972 RecMdcTrackCol::iterator it_trk = newtrkCol->begin();
1973 for ( it_trk = newtrkCol->begin(); it_trk != newtrkCol->end(); it_trk++ )
1974 {
1975 int nRecHits = ( *it_trk )->getNhits();
1976 if ( nRecHits < nHits ) continue;
1977
1978 int hitInRecTrk[100];
1979 int iRecHit = 0;
1980 HitRefVec gothits = ( *it_trk )->getVecHits();
1981 HitRefVec::iterator it_hit = gothits.begin();
1982 for ( ; it_hit != gothits.end(); it_hit++ )
1983 {
1984 identifier = ( *it_hit )->getMdcId();
1985 int lay = mdcid.layer( identifier );
1986 int cel = mdcid.wire( identifier );
1987 int wir = m_mdcGeomSvc->Wire( lay, cel )->Id();
1988 hitInRecTrk[iRecHit] = wir;
1989 iRecHit++;
1990 }
1991
1992 // match the track
1993 bool matchSuccess = true;
1994 for ( int i = 0; i < nHits; i++ )
1995 {
1996 bool findHit = false;
1997 for ( int k = 0; k < nRecHits; k++ )
1998 {
1999 if ( hitInTrk[i] == hitInRecTrk[k] )
2000 {
2001 findHit = true;
2002 break;
2003 }
2004 }
2005 if ( !findHit )
2006 {
2007 matchSuccess = false;
2008 break;
2009 }
2010 }
2011 if ( !matchSuccess ) continue;
2012
2013 HepVector helix = ( *it_trk )->helix();
2014 HepSymMatrix helixErr = ( *it_trk )->err();
2015 double phi0 = ( *it_trk )->helix( 1 );
2016 double phiTrk = phi0 + HFPI;
2017 if ( phiTrk > PI2 ) phiTrk -= PI2;
2018
2019 for ( int lay = 0; lay < 43; lay++ )
2020 {
2021 double docamin = 0.9; // cm
2022 if ( lay < 8 ) docamin = 0.6; // cm
2023 int celmin = -1;
2024 int ncel = m_mdcGeomSvc->Layer( lay )->NCell();
2025 for ( int cel = 0; cel < ncel; cel++ )
2026 {
2027 double wphi;
2028 const MdcGeoWire* pWire = m_mdcGeomSvc->Wire( lay, cel );
2029 double xx = pWire->Backward().x();
2030 double yy = pWire->Backward().y();
2031 double rr = sqrt( ( xx * xx ) + ( yy * yy ) );
2032 if ( yy >= 0 ) wphi = acos( xx / rr );
2033 else wphi = CLHEP::twopi - acos( xx / rr );
2034
2035 if ( !( ( fabs( wphi - phiTrk ) < dphi ) || ( fabs( wphi + PI2 - phiTrk ) < dphi ) ||
2036 ( fabs( wphi - PI2 - phiTrk ) < dphi ) ) )
2037 { continue; }
2038
2039 double doca = m_mdcUtilitySvc->doca( lay, cel, helix, helixErr );
2040 // cout << setw(5) << lay << setw(5) << cel << setw(15) << doca <<
2041 // endl;
2042 if ( fabs( doca ) < fabs( docamin ) )
2043 {
2044 docamin = doca;
2045 celmin = cel;
2046 }
2047 }
2048 if ( celmin > -1 ) { cellTrkPass[lay] = celmin; }
2049 else { cellTrkPass[lay] = -1; }
2050 }
2051 return true;
2052 }
2053 return false;
2054}
mass
double mass
Definition CosmicGenerator.cxx:128
p_cms
const HepLorentzVector p_cms(0.034067, 0.0, 0.0, 3.097)
sprintf
sprintf(cut, "kal_costheta0_em>-0.93&&kal_costheta0_em<0.93&&kal_pxy0_em>=0.05+%d*0.1&&kal_" "pxy0_em<0.15+%d*0.1&&NGch>=2", j, j)
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:83
bin
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per bin
Definition FoamA.h:85
Vp4
std::vector< HepLorentzVector > Vp4
Definition Gam4pikp.cxx:38
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
HitRefVec
SmartRefVector< RecMdcHit > HitRefVec
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcRecEvent/RecMdcTrack.h:22
MdcCalNLayer
const int MdcCalNLayer
Definition MdcCalParams.h:6
MdcCalNENTRSD
const int MdcCalNENTRSD
Definition MdcCalParams.h:19
MdcCalTdcCnv
const double MdcCalTdcCnv
Definition MdcCalParams.h:24
MdcCalSdNBIN
const int MdcCalSdNBIN
Definition MdcCalParams.h:20
MdcCalNENTRXT
const int MdcCalNENTRXT
Definition MdcCalParams.h:12
MdcCalTotCell
const int MdcCalTotCell
Definition MdcCalParams.h:9
valType3
map< int, int >::value_type valType3
Definition MdcCalib.cxx:33
ntupleSvc
INTupleSvc * ntupleSvc()
Definition ServiceAccessor.h:51
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:12
key
*************DOUBLE PRECISION m_pi *DOUBLE PRECISION m_HvecTau2 DOUBLE PRECISION m_HvClone2 DOUBLE PRECISION m_gamma1 DOUBLE PRECISION m_gamma2 DOUBLE PRECISION m_thet1 DOUBLE PRECISION m_thet2 INTEGER m_IFPHOT *COMMON c_Taupair $ !Spin Polarimeter vector first Tau $ !Spin Polarimeter vector second Tau $ !Clone Spin Polarimeter vector first Tau $ !Clone Spin Polarimeter vector second Tau $ !Random Euler angle for cloning st tau $ !Random Euler angle for cloning st tau $ !Random Euler angle for cloning st tau $ !Random Euler angle for cloning nd tau $ !Random Euler angle for cloning nd tau $ !Random Euler angle for cloning nd tau $ !phi of HvecTau1 $ !theta of HvecTau1 $ !phi of HvecTau2 $ !theta of HvecTau2 $ !super key
Definition Taupair.h:42
MdcCalEvent::getRecTrk
MdcCalRecTrk * getRecTrk(int index) const
Definition MdcCalEvent.h:34
MdcCalRecHit::getDocaExc
double getDocaExc() const
Definition MdcCalRecHit.h:28
MdcCalRecHit::getDmeas
double getDmeas() const
Definition MdcCalRecHit.h:30
MdcCalRecHit::getResiExcLR
double getResiExcLR() const
Definition MdcCalRecHit.h:37
MdcCalRecHit::getCellid
int getCellid() const
Definition MdcCalRecHit.h:23
MdcCalRecHit::getTdrift
double getTdrift() const
Definition MdcCalRecHit.h:40
MdcCalRecHit::getQhit
double getQhit() const
Definition MdcCalRecHit.h:41
MdcCalRecHit::getEntra
double getEntra() const
Definition MdcCalRecHit.h:43
MdcCalRecHit::getLayid
int getLayid() const
Definition MdcCalRecHit.h:22
MdcCalRecHit::getResiExc
double getResiExc() const
Definition MdcCalRecHit.h:36
MdcCalRecHit::getResiIncLR
double getResiIncLR() const
Definition MdcCalRecHit.h:34
MdcCalRecHit::getZhit
double getZhit() const
Definition MdcCalRecHit.h:42
MdcCalRecHit::getLR
int getLR() const
Definition MdcCalRecHit.h:24
MdcCalRecHit::getTdc
double getTdc() const
Definition MdcCalRecHit.h:39
MdcCalRecHit::getStat
int getStat() const
Definition MdcCalRecHit.h:25
MdcCalRecTrk::getRecHit
MdcCalRecHit * getRecHit(int index) const
Definition MdcCalRecTrk.h:40
MdcCalRecTrk::getP
double getP() const
Definition MdcCalRecTrk.h:36
MdcCalRecTrk::getTanLamda
double getTanLamda() const
Definition MdcCalRecTrk.h:34
MdcCalRecTrk::getDz
double getDz() const
Definition MdcCalRecTrk.h:33
MdcCalRecTrk::getNHits
int getNHits() const
Definition MdcCalRecTrk.h:38
MdcCalRecTrk::getPhi0
double getPhi0() const
Definition MdcCalRecTrk.h:31
MdcCalRecTrk::getP4
HepLorentzVector getP4() const
Definition MdcCalRecTrk.h:39
MdcCalRecTrk::getKappa
double getKappa() const
Definition MdcCalRecTrk.h:32
MdcCalRecTrk::getChisq
double getChisq() const
Definition MdcCalRecTrk.h:37
MdcCalRecTrk::getPt
double getPt() const
Definition MdcCalRecTrk.h:35
MdcCalRecTrk::getDr
double getDr() const
Definition MdcCalRecTrk.h:30
MdcCalibConst::resetSdpar
void resetSdpar(int lay, int entr, int lr, int bin, double val)
Definition MdcCalibConst.cpp:84
MdcCalibConst::getSdpar
double getSdpar(int lay, int entr, int lr, int bin)
Definition MdcCalibConst.cpp:89
MdcCalib::printCut
virtual void printCut() const =0
Definition MdcCalib.cxx:1509
MdcCalib::clear
virtual void clear()=0
Definition MdcCalib.cxx:83
MdcCalib::~MdcCalib
virtual ~MdcCalib()
Definition MdcCalib.cxx:81
MdcCalib::initialize
virtual void initialize(TObjArray *hlist, IMdcGeomSvc *mdcGeomSvc, IMdcCalibFunSvc *mdcFunSvc, IMdcUtilitySvc *mdcUtilitySvc)=0
Definition MdcCalib.cxx:225
MdcCalib::updateConst
virtual int updateConst(MdcCalibConst *calconst)=0
Definition MdcCalib.cxx:1518
MdcCalib::fillHist
virtual int fillHist(MdcCalEvent *event)=0
Definition MdcCalib.cxx:767
MdcID::layer
static int layer(const Identifier &id)
Values of different levels (failure returns 0).
Definition MdcID.cxx:47
MdcID::wire
static int wire(const Identifier &id)
Definition MdcID.cxx:52
RawData::getChargeChannel
unsigned int getChargeChannel() const
Definition RawData.cxx:35
RawData::getTimeChannel
unsigned int getTimeChannel() const
Definition RawData.cxx:32