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