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