00001 #include "MdcCalibAlg/MdcCalib.h"
00002
00003 #include "GaudiKernel/MsgStream.h"
00004 #include "GaudiKernel/IMessageSvc.h"
00005 #include "GaudiKernel/StatusCode.h"
00006 #include "GaudiKernel/ISvcLocator.h"
00007 #include "GaudiKernel/Bootstrap.h"
00008 #include "GaudiKernel/SmartDataPtr.h"
00009 #include "GaudiKernel/IDataProviderSvc.h"
00010 #include "GaudiKernel/PropertyMgr.h"
00011 #include "GaudiKernel/INTupleSvc.h"
00012 #include "GaudiKernel/IService.h"
00013
00014 #include "EventModel/Event.h"
00015 #include "MdcRawEvent/MdcDigi.h"
00016 #include "McTruth/MdcMcHit.h"
00017 #include "Identifier/Identifier.h"
00018 #include "Identifier/MdcID.h"
00019 #include "EventModel/EventHeader.h"
00020 #include "CLHEP/Vector/ThreeVector.h"
00021
00022 #include "TStyle.h"
00023
00024 #include <iostream>
00025 #include <fstream>
00026 #include <iomanip>
00027 #include <string>
00028 #include <cstring>
00029
00030 using namespace Event;
00031 using namespace std;
00032
00033 typedef map<int, int>::value_type valType3;
00034 typedef std::vector<HepLorentzVector> Vp4;
00035
00036 MdcCalib::MdcCalib(){
00037 m_nEvt=0;
00038 m_cut1=0;
00039 m_cut2=0;
00040 m_cut3=0;
00041 m_cut4=0;
00042 m_cut5=0;
00043 m_cut6=0;
00044 m_nGrPoint = 0;
00045 fgReadWireEff = false;
00046
00047 int lay;
00048 for(lay=0; lay<43; lay++){
00049 if(lay < 15) m_nEntr[lay] = 1;
00050 else m_nEntr[lay] = 2;
00051 }
00052 m_dwid = 0.5;
00053 m_fgIni = false;
00054
00055 m_phiWid = PI2 / (double)NPhiBin;
00056 m_theWid = 2.0 / (double)NThetaBin;
00057
00058 m_nEvtNtuple = 0;
00059
00060 for(lay=0; lay<MdcCalNLayer; lay++){
00061 if(lay < 8) m_nBin[lay] = 12;
00062 else m_nBin[lay] = 16;
00063 }
00064
00065
00066 for(lay=0; lay<MdcCalNLayer; lay++){
00067 if((0==lay) || (7==lay) || (8==lay) || (19==lay) || (20==lay) ||
00068 (35==lay) || (36==lay) || (42==lay) ) m_layBound[lay] = true;
00069 else m_layBound[lay] = false;
00070 }
00071 }
00072
00073 MdcCalib::~MdcCalib(){
00074 }
00075
00076 void MdcCalib::clear(){
00077 int lay;
00078 for(lay=0; lay<m_nlayer; lay++){
00079 delete m_htraw[lay];
00080 delete m_htdr[lay];
00081 delete m_hresInc[lay];
00082 delete m_hresExc[lay];
00083 delete m_hresAve[lay];
00084 delete m_hadc[lay];
00085 for (int lr=0; lr<2; lr++){
00086 delete m_htdrlr[lay][lr];
00087 delete m_hreslrInc[lay][lr];
00088 delete m_hreslrExc[lay][lr];
00089 delete m_hreslrAve[lay][lr];
00090 }
00091 }
00092
00093 delete m_effNtrk;
00094 delete m_effNtrkRecHit;
00095 delete m_hresAllInc;
00096 delete m_hresAllExc;
00097 delete m_hresAllAve;
00098 for(int i=0; i<14; i++){
00099 delete m_hresAveAllQ[i];
00100 delete m_hresAveOutQ[i];
00101 }
00102 for(lay=0; lay<43; lay++){
00103 for(int i=0; i<14; i++) delete m_hresAveLayQ[lay][i];
00104 }
00105 delete m_hresInnInc;
00106 delete m_hresInnExc;
00107 delete m_hresStpInc;
00108 delete m_hresStpExc;
00109 delete m_hresOutInc;
00110 delete m_hresOutExc;
00111 for(int iEs=0; iEs<m_param.nEsFlag; iEs++) delete m_hTes[iEs];
00112 delete m_hbbTrkFlg;
00113 delete m_hTesAll;
00114 delete m_hTesGood;
00115 delete m_hTesAllFlag;
00116 delete m_hTesRec;
00117 delete m_hTesCalFlag;
00118 delete m_hTesCalUse;
00119 delete m_hnRawHit;
00120 delete m_hpt;
00121 delete m_hpMax;
00122 delete m_hpMaxCms;
00123 delete m_hptPos;
00124 delete m_hptNeg;
00125 delete m_hp;
00126 delete m_hp_cms;
00127 delete m_hpPos;
00128 delete m_hpNeg;
00129 delete m_hpPoscms;
00130 delete m_hpNegcms;
00131 delete m_hp_cut;
00132 delete m_hchisq;
00133 delete m_hnTrk;
00134 delete m_hnTrkCal;
00135 delete m_hnhitsRec;
00136 delete m_hnhitsRecInn;
00137 delete m_hnhitsRecStp;
00138 delete m_hnhitsRecOut;
00139 delete m_hnhitsCal;
00140 delete m_hnhitsCalInn;
00141 delete m_hnhitsCalStp;
00142 delete m_hnhitsCalOut;
00143 delete m_wirehitmap;
00144 delete m_layerhitmap;
00145 delete m_hnoisephi;
00146 delete m_hnoiselay;
00147 delete m_hnoisenhits;
00148 delete m_hratio;
00149 delete m_hdr;
00150 delete m_hphi0;
00151 delete m_hkap;
00152 delete m_hdz;
00153 delete m_htanl;
00154 delete m_hcosthe;
00155 delete m_hcostheNeg;
00156 delete m_hcosthePos;
00157 delete m_hx0;
00158 delete m_hy0;
00159 delete m_hdelZ0;
00160 delete m_grX0Y0;
00161 delete m_hitEffAll;
00162 delete m_hitEffRaw;
00163 delete m_hitEffRec;
00164 int bin;
00165 int thbin;
00166 for(bin=0; bin<NPhiBin; bin++){
00167 delete m_ppPhi[bin];
00168 delete m_pnPhi[bin];
00169 delete m_ppPhiCms[bin];
00170 delete m_pnPhiCms[bin];
00171 for(thbin=0; thbin<NThetaBin; thbin++){
00172 delete m_ppThePhi[thbin][bin];
00173 delete m_pnThePhi[thbin][bin];
00174 delete m_ppThePhiCms[thbin][bin];
00175 delete m_pnThePhiCms[thbin][bin];
00176 }
00177 }
00178 for(thbin=0; thbin<NThetaBin; thbin++){
00179 delete m_ppThe[thbin];
00180 delete m_pnThe[thbin];
00181 delete m_ppTheCms[thbin];
00182 delete m_pnTheCms[thbin];
00183 }
00184
00185 for(unsigned i=0; i<m_hr2dInc.size(); i++){
00186 delete m_hr2dInc[i];
00187 delete m_hr2dExc[i];
00188 }
00189 m_hr2dInc.clear();
00190 m_hr2dExc.clear();
00191 m_mapr2d.clear();
00192
00193 delete m_fdTime;
00194 delete m_fdAdc;
00195 delete m_fdres;
00196 delete m_fdresAve;
00197 delete m_fdres2d;
00198 delete m_fdcom;
00199 delete m_fdResQ;
00200 }
00201
00202 void MdcCalib::initialize(TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
00203 IMdcCalibFunSvc* mdcFunSvc, IMdcUtilitySvc* mdcUtilitySvc) {
00204 IMessageSvc* msgSvc;
00205 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00206 MsgStream log(msgSvc, "MdcCalib");
00207 log << MSG::INFO << "MdcCalib::initialize()" << endreq;
00208
00209 m_hlist = hlist;
00210 m_mdcGeomSvc = mdcGeomSvc;
00211 m_mdcFunSvc = mdcFunSvc;
00212 m_mdcUtilitySvc = mdcUtilitySvc;
00213
00214 int lay;
00215 int iEntr;
00216 int lr;
00217 int bin;
00218 char hname[200];
00219
00220 m_nlayer = m_mdcGeomSvc -> getLayerSize();
00221
00222 for(lay=0; lay<m_nlayer; lay++){
00223 m_radii[lay] = m_mdcGeomSvc->Layer(lay)->Radius();
00224 }
00225 ofstream fwpc("wirelog.txt");
00226 for(int wir=0; wir<MdcCalTotCell; wir++){
00227 m_xe[wir] = m_mdcGeomSvc->Wire(wir)->Backward().x();
00228 m_ye[wir] = m_mdcGeomSvc->Wire(wir)->Backward().y();
00229 m_ze[wir] = m_mdcGeomSvc->Wire(wir)->Backward().z();
00230 m_xw[wir] = m_mdcGeomSvc->Wire(wir)->Forward().x();
00231 m_yw[wir] = m_mdcGeomSvc->Wire(wir)->Forward().y();
00232 m_zw[wir] = m_mdcGeomSvc->Wire(wir)->Forward().z();
00233 fwpc << setw(6) << wir << setw(15) << m_xe[wir] << setw(15) << m_ye[wir]
00234 << setw(15) << m_ze[wir] << setw(15) << m_xw[wir]
00235 << setw(15) << m_yw[wir] << setw(15) << m_zw[wir] << endl;
00236 }
00237 fwpc.close();
00238
00239 m_fdcom = new TFolder("common", "common");
00240 m_hlist -> Add(m_fdcom);
00241
00242 m_effNtrk = new TH1F("effNtrk", "", 43, -0.5, 42.5);
00243 m_fdcom->Add(m_effNtrk);
00244
00245 m_effNtrkRecHit = new TH1F("effNtrkRecHit", "", 43, -0.5, 42.5);
00246 m_fdcom->Add(m_effNtrkRecHit);
00247
00248 m_hresAllInc = new TH1F("HResAllInc", "", 200, -1.0, 1.0);
00249 m_fdcom -> Add(m_hresAllInc);
00250
00251 m_hresAllExc = new TH1F("HResAllExc", "", 200, -1.0, 1.0);
00252 m_fdcom -> Add(m_hresAllExc);
00253
00254 m_hresAllAve = new TH1F("HResAllAve", "", 200, -1.0, 1.0);
00255 m_fdcom -> Add(m_hresAllAve);
00256
00257 m_hresInnInc = new TH1F("HResInnInc", "", 200, -1.0, 1.0);
00258 m_fdcom -> Add(m_hresInnInc);
00259
00260 m_hresInnExc = new TH1F("HResInnExc", "", 200, -1.0, 1.0);
00261 m_fdcom -> Add(m_hresInnExc);
00262
00263 m_hresStpInc = new TH1F("HResStpInc", "", 200, -1.0, 1.0);
00264 m_fdcom -> Add(m_hresStpInc);
00265
00266 m_hresStpExc = new TH1F("HResStpExc", "", 200, -1.0, 1.0);
00267 m_fdcom -> Add(m_hresStpExc);
00268
00269 m_hresOutInc = new TH1F("HResOutInc", "", 200, -1.0, 1.0);
00270 m_fdcom -> Add(m_hresOutInc);
00271
00272 m_hresOutExc = new TH1F("HResOutExc", "", 200, -1.0, 1.0);
00273 m_fdcom -> Add(m_hresOutExc);
00274
00275 m_fdResQ = new TFolder("ResQ", "ResQ");
00276 m_hlist->Add(m_fdResQ);
00277 for(int i=0; i<14; i++){
00278 sprintf(hname, "resoAll_qbin%02d", i);
00279 m_hresAveAllQ[i] = new TH1F(hname, "", 200, -1, 1);
00280 m_fdResQ->Add(m_hresAveAllQ[i]);
00281
00282 sprintf(hname, "resoOut_qbin%02d", i);
00283 m_hresAveOutQ[i] = new TH1F(hname, "", 200, -1, 1);
00284 m_fdResQ->Add(m_hresAveOutQ[i]);
00285 }
00286 for(lay=0; lay<43; lay++){
00287 for(int i=0; i<14; i++){
00288 sprintf(hname, "resoLay%02d_qbin%02d", lay, i);
00289 m_hresAveLayQ[lay][i] = new TH1F(hname, "", 200, -1, 1);
00290 m_fdResQ->Add(m_hresAveLayQ[lay][i]);
00291 }
00292 }
00293
00294 for(int iEs=0; iEs<m_param.nEsFlag; iEs++){
00295 sprintf(hname, "Tes_%d", m_param.esFlag[iEs]);
00296 m_hTes[iEs] = new TH1F(hname, "", 750, 0, 1500);
00297 m_fdcom->Add(m_hTes[iEs]);
00298 }
00299
00300 m_hbbTrkFlg = new TH1F("BbTrkFlg", "", 100, 0, 6);
00301 m_fdcom -> Add(m_hbbTrkFlg);
00302
00303 m_hTesAll = new TH1F("TesAll", "", 1000, 0, 2000);
00304 m_fdcom -> Add(m_hTesAll);
00305
00306 m_hTesGood = new TH1F("TesGood", "", 1000, 0, 2000);
00307 m_fdcom -> Add(m_hTesGood);
00308
00309 m_hTesAllFlag = new TH1F("TesAllFlag", "", 300, -0.5, 299.5);
00310 m_fdcom -> Add(m_hTesAllFlag);
00311
00312 m_hTesRec = new TH1F("TesRec", "", 1000, 0, 2000);
00313 m_fdcom -> Add(m_hTesRec);
00314
00315 m_hTesCalFlag = new TH1F("TesCalFlag", "", 1000, 0, 2000);
00316 m_fdcom -> Add(m_hTesCalFlag);
00317
00318 m_hTesCalUse = new TH1F("TesCalUse", "", 1000, 0, 2000);
00319 m_fdcom -> Add(m_hTesCalUse);
00320
00321 m_hnRawHit = new TH1F("nRawHit", "", 6797, -0.5, 6796.5);
00322 m_fdcom -> Add(m_hnRawHit);
00323
00324 m_hpt = new TH1F("HPt", "", 800, 0, 3);
00325 m_fdcom -> Add(m_hpt);
00326
00327 m_hptPos = new TH1F("HPtPos", "", 800, 0, 3);
00328 m_fdcom -> Add(m_hptPos);
00329
00330 m_hptNeg = new TH1F("HPtNeg", "", 800, 0, 3);
00331 m_fdcom -> Add(m_hptNeg);
00332
00333 m_hp = new TH1F("HP", "", 800, 0, 3);
00334 m_fdcom -> Add(m_hp);
00335
00336 m_hp_cms = new TH1F("HPCMS", "", 800, 0, 3);
00337 m_fdcom -> Add(m_hp_cms);
00338
00339 m_hpMax = new TH1F("HPMax", "", 800, 0, 3);
00340 m_fdcom -> Add(m_hpMax);
00341
00342 m_hpMaxCms = new TH1F("HPMax_Cms", "", 800, 0, 3);
00343 m_fdcom -> Add(m_hpMaxCms);
00344
00345 m_hpPos = new TH1F("HP_Pos", "", 800, 0, 3);
00346 m_fdcom -> Add(m_hpPos);
00347
00348 m_hpNeg = new TH1F("HP_Neg", "", 800, 0, 3);
00349 m_fdcom -> Add(m_hpNeg);
00350
00351 m_hpPoscms = new TH1F("HP_Pos_cms", "", 800, 0, 3);
00352 m_fdcom -> Add(m_hpPoscms);
00353
00354 m_hpNegcms = new TH1F("HP_Neg_cms", "", 800, 0, 3);
00355 m_fdcom -> Add(m_hpNegcms);
00356
00357 m_hp_cut = new TH1F("HPCut", "", 800, 0, 3);
00358 m_fdcom -> Add(m_hp_cut);
00359
00360 m_hchisq = new TH1F("Chisq", "", 10, 0, 100);
00361 m_fdcom -> Add(m_hchisq);
00362
00363 m_hnTrk = new TH1F("HNtrack", "HNtrack", 10, -0.5, 9.5);
00364 m_fdcom -> Add(m_hnTrk);
00365
00366 m_hnTrkCal = new TH1F("HNtrackCal", "HNtrackCal", 10, -0.5, 9.5);
00367 m_fdcom -> Add(m_hnTrkCal);
00368
00369 m_hnhitsRec = new TH1F("HNhitsRec", "", 100, -0.5, 99.5);
00370 m_fdcom -> Add(m_hnhitsRec);
00371
00372 m_hnhitsRecInn = new TH1F("HNhitsInnRec", "", 60, 0.5, 60.5);
00373 m_fdcom -> Add(m_hnhitsRecInn);
00374
00375 m_hnhitsRecStp = new TH1F("HNhitsStpRec", "", 60, 0.5, 60.5);
00376 m_fdcom -> Add(m_hnhitsRecStp);
00377
00378 m_hnhitsRecOut = new TH1F("HNhitsOutRec", "", 60, 0.5, 60.5);
00379 m_fdcom -> Add(m_hnhitsRecOut);
00380
00381 m_hnhitsCal = new TH1F("HNhitsCal", "", 100, -0.5, 99.5);
00382 m_fdcom -> Add(m_hnhitsCal);
00383
00384 m_hnhitsCalInn = new TH1F("HNhitsCalInn", "", 60, 0.5, 60.5);
00385 m_fdcom -> Add(m_hnhitsCalInn);
00386
00387 m_hnhitsCalStp = new TH1F("HNhitsCalStp", "", 60, 0.5, 60.5);
00388 m_fdcom -> Add(m_hnhitsCalStp);
00389
00390 m_hnhitsCalOut = new TH1F("HNhitsCalOut", "", 60, 0.5, 60.5);
00391 m_fdcom -> Add(m_hnhitsCalOut);
00392
00393 m_wirehitmap = new TH1F("Wire_HitMap", "Wire_HitMap", 6796, -0.5, 6795.5);
00394 m_fdcom -> Add(m_wirehitmap);
00395
00396 m_layerhitmap = new TH1F("Layer_HitMap", "Layer_HitMap", 43, -0.5, 42.5);
00397 m_fdcom -> Add(m_layerhitmap);
00398
00399 m_hnoisephi = new TH1F("phi_noise", "", 100, 0, 6.284);
00400 m_fdcom -> Add(m_hnoisephi);
00401
00402 m_hnoiselay = new TH1F("Layer_noise", "Layer_noise", 43, -0.5, 42.5);
00403 m_fdcom -> Add(m_hnoiselay);
00404
00405 m_hnoisenhits = new TH1F("nhits_noise", "nhits_noise", 6796, -0.5, 6795.5);
00406 m_fdcom -> Add(m_hnoisenhits);
00407
00408 m_hratio = new TH1F("ratio", "", 100, 0, 1);
00409 m_fdcom -> Add(m_hratio);
00410
00411 m_hdr = new TH1F("dr", "", 500, -500, 500);
00412 m_fdcom -> Add(m_hdr);
00413
00414 m_hphi0 = new TH1F("phi0", "", 100, 0, 6.284);
00415 m_fdcom -> Add(m_hphi0);
00416
00417 m_hkap = new TH1F("kappa", "", 400, -50, 50);
00418 m_fdcom -> Add(m_hkap);
00419
00420 m_hdz = new TH1F("dz", "", 500, -1000, 1000);
00421 m_fdcom -> Add(m_hdz);
00422
00423 m_htanl = new TH1F("tanl", "", 200, -5, 5);
00424 m_fdcom -> Add(m_htanl);
00425
00426 m_hcosthe = new TH1F("costheta", "", 200, -1, 1);
00427 m_fdcom -> Add(m_hcosthe);
00428
00429 m_hcostheNeg = new TH1F("costhetaNeg", "", 200, -1, 1);
00430 m_fdcom -> Add(m_hcostheNeg);
00431
00432 m_hcosthePos = new TH1F("costhetaPos", "", 200, -1, 1);
00433 m_fdcom -> Add(m_hcosthePos);
00434
00435 m_hx0 = new TH1F("x0", "", 100, -10, 10);
00436 m_fdcom -> Add(m_hx0);
00437
00438 m_hy0 = new TH1F("y0", "", 100, -10, 10);
00439 m_fdcom -> Add(m_hy0);
00440
00441 m_hdelZ0 = new TH1F("delta_z0", "", 100, -50, 50);
00442 m_fdcom -> Add(m_hdelZ0);
00443
00444 m_grX0Y0 = new TGraph();
00445 m_grX0Y0->SetName("x0y0");
00446 m_fdcom -> Add(m_grX0Y0);
00447
00448 m_hitEffAll = new TH1F("hitEffAll", "", 6800, -0.5, 6799.5);
00449 m_fdcom -> Add(m_hitEffAll);
00450
00451 m_hitEffRaw = new TH1F("hitEffRaw", "", 6800, -0.5, 6799.5);
00452 m_fdcom -> Add(m_hitEffRaw);
00453
00454 m_hitEffRec = new TH1F("hitEffRec", "", 6800, -0.5, 6799.5);
00455 m_fdcom -> Add(m_hitEffRec);
00456
00457
00458 m_fdTime = new TFolder("time", "time");
00459 m_hlist -> Add(m_fdTime);
00460
00461 for(lay=0; lay<m_nlayer; lay++){
00462 sprintf(hname, "Traw%02d", lay);
00463 m_htraw[lay] = new TH1F(hname, "", 1000, 0, 2000);
00464 m_fdTime -> Add(m_htraw[lay]);
00465
00466 sprintf(hname, "Tdr%02d", lay);
00467 m_htdr[lay] = new TH1F(hname, "", 510, -10, 500);
00468 m_fdTime -> Add(m_htdr[lay]);
00469
00470 for (lr=0; lr<2; lr++){
00471 sprintf(hname, "Tdr%02d_lr%01d", lay, lr);
00472 m_htdrlr[lay][lr] = new TH1F(hname, "", 510, -10, 500);
00473 m_fdTime -> Add(m_htdrlr[lay][lr]);
00474 }
00475 }
00476
00477
00478 m_fdAdc = new TFolder("adc", "adc");
00479 m_hlist -> Add(m_fdAdc);
00480
00481 for(lay=0; lay<m_nlayer; lay++){
00482 sprintf(hname, "adc%02d", lay);
00483 m_hadc[lay] = new TH1F(hname, "", 1500, 0, 3000);
00484 m_fdAdc -> Add(m_hadc[lay]);
00485 }
00486
00487 m_fdres = new TFolder("resolution", "resolution");
00488 m_hlist -> Add(m_fdres);
00489
00490 m_fdresAve = new TFolder("resAve", "resAve");
00491 m_hlist -> Add(m_fdresAve);
00492
00493 for(lay=0; lay<m_nlayer; lay++){
00494 sprintf(hname, "Reso%02dInc", lay);
00495 m_hresInc[lay] = new TH1F(hname, "", 1000, -5, 5);
00496 m_fdres -> Add(m_hresInc[lay]);
00497
00498 sprintf(hname, "Reso%02dExc", lay);
00499 m_hresExc[lay] = new TH1F(hname, "", 1000, -5, 5);
00500 m_fdres -> Add(m_hresExc[lay]);
00501
00502 sprintf(hname, "Reso%02d", lay);
00503 m_hresAve[lay] = new TH1F(hname, "", 1000, -5, 5);
00504 m_fdresAve -> Add(m_hresAve[lay]);
00505
00506 for (lr=0; lr<2; lr++){
00507 sprintf(hname, "Reso%02dInc_lr%01d", lay, lr);
00508
00509 m_hreslrInc[lay][lr] = new TH1F(hname, "", 1000, -5, 5);
00510 m_fdres->Add(m_hreslrInc[lay][lr]);
00511
00512 sprintf(hname, "Reso%02dExc_lr%01d", lay, lr);
00513
00514 m_hreslrExc[lay][lr] = new TH1F(hname, "", 1000, -5, 5);
00515 m_fdres->Add(m_hreslrExc[lay][lr]);
00516
00517 sprintf(hname, "Reso%02d_lr%01d", lay, lr);
00518
00519 m_hreslrAve[lay][lr] = new TH1F(hname, "", 1000, -5, 5);
00520 m_fdresAve->Add(m_hreslrAve[lay][lr]);
00521 }
00522 for(int phi=0; phi<20; phi++){
00523 sprintf(hname, "ResoPhi%02d_phi%02d", lay, phi);
00524 m_hresphi[lay][phi] = new TH1F(hname, "", 200, -1, 1);
00525 m_fdres->Add(m_hresphi[lay][phi]);
00526 }
00527 }
00528
00529
00530 m_fdmomPhi = new TFolder("momPhi", "momPhi");
00531 m_hlist -> Add(m_fdmomPhi);
00532
00533 int thbin;
00534 for(bin=0; bin<NPhiBin; bin++){
00535 sprintf(hname, "hPpos_phi%02d", bin);
00536 m_ppPhi[bin] = new TH1F(hname, "", 400, 1.0, 2.5);
00537 m_fdmomPhi->Add(m_ppPhi[bin]);
00538
00539 sprintf(hname, "hPneg_phi%02d", bin);
00540 m_pnPhi[bin] = new TH1F(hname, "", 400, 1.0, 2.5);
00541 m_fdmomPhi->Add(m_pnPhi[bin]);
00542
00543 sprintf(hname, "hPpos_phi_cms%02d", bin);
00544 m_ppPhiCms[bin] = new TH1F(hname, "", 400, 1.0, 2.5);
00545 m_fdmomPhi->Add(m_ppPhiCms[bin]);
00546
00547 sprintf(hname, "hPneg_phi_cms%02d", bin);
00548 m_pnPhiCms[bin] = new TH1F(hname, "", 400, 1.0, 2.5);
00549 m_fdmomPhi->Add(m_pnPhiCms[bin]);
00550
00551 for(thbin=0; thbin<NThetaBin; thbin++){
00552 sprintf(hname, "hPpos_theta%02d_phi%02d", thbin, bin);
00553 m_ppThePhi[thbin][bin] = new TH1F(hname, "", 400, 1.0, 2.5);
00554 m_fdmomPhi->Add(m_ppThePhi[thbin][bin]);
00555
00556 sprintf(hname, "hPneg_theta%02d_phi%02d", thbin, bin);
00557 m_pnThePhi[thbin][bin] = new TH1F(hname, "", 400, 1.0, 2.5);
00558 m_fdmomPhi->Add(m_pnThePhi[thbin][bin]);
00559
00560 sprintf(hname, "hPposCms_theta%02d_phi%02d", thbin, bin);
00561 m_ppThePhiCms[thbin][bin] = new TH1F(hname, "", 400, 1.0, 2.5);
00562 m_fdmomPhi->Add(m_ppThePhiCms[thbin][bin]);
00563
00564 sprintf(hname, "hPnegCms_theta%02d_phi%02d", thbin, bin);
00565 m_pnThePhiCms[thbin][bin] = new TH1F(hname, "", 400, 1.0, 2.5);
00566 m_fdmomPhi->Add(m_pnThePhiCms[thbin][bin]);
00567 }
00568 }
00569 for(thbin=0; thbin<NThetaBin; thbin++){
00570 sprintf(hname, "hPpos_the%02d", thbin);
00571 m_ppThe[thbin] = new TH1F(hname, "", 400, 1.0, 2.5);
00572 m_fdmomPhi->Add(m_ppThe[thbin]);
00573
00574 sprintf(hname, "hPneg_the%02d", thbin);
00575 m_pnThe[thbin] = new TH1F(hname, "", 400, 1.0, 2.5);
00576 m_fdmomPhi->Add(m_pnThe[thbin]);
00577
00578 sprintf(hname, "hPposCms_the%02d", thbin);
00579 m_ppTheCms[thbin] = new TH1F(hname, "", 400, 1.0, 2.5);
00580 m_fdmomPhi->Add(m_ppTheCms[thbin]);
00581
00582 sprintf(hname, "hPnegCms_the%02d", thbin);
00583 m_pnTheCms[thbin] = new TH1F(hname, "", 400, 1.0, 2.5);
00584 m_fdmomPhi->Add(m_pnTheCms[thbin]);
00585 }
00586
00587
00588 m_fdres2d = new TFolder("res2d", "res2d");
00589 m_hlist -> Add(m_fdres2d);
00590
00591 int hid = 0;
00592 int key;
00593 TH1F* hist;
00594 for(lay=0; lay<m_nlayer; lay++){
00595 for(iEntr=0; iEntr<MdcCalNENTRSD; iEntr++){
00596 for(lr=0; lr<2; lr++){
00597 for(bin=0; bin<MdcCalSdNBIN; bin++){
00598 sprintf(hname, "r2d%02d_%02d_%01d_%02dInc", lay, iEntr, lr, bin);
00599 hist = new TH1F(hname, "", 200, -1, 1);
00600 m_hr2dInc.push_back(hist);
00601 m_fdres2d -> Add(hist);
00602
00603 sprintf(hname, "r2d%02d_%02d_%01d_%02dExc", lay, iEntr, lr, bin);
00604 hist = new TH1F(hname, "", 200, -1, 1);
00605 m_hr2dExc.push_back(hist);
00606 m_fdres2d -> Add(hist);
00607
00608 key = getHresId(lay, iEntr, lr, bin);
00609 m_mapr2d.insert( valType3(key, hid) );
00610 hid++;
00611 }
00612 }
00613 }
00614 }
00615
00616 m_fdres2t = new TFolder("res2t", "res2t");
00617
00618
00619 for(lay=0; lay<m_nlayer; lay++){
00620 for(iEntr=0; iEntr<MdcCalNENTRXT; iEntr++){
00621 for(lr=0; lr<2; lr++){
00622 for(bin=0; bin<45; bin++){
00623 sprintf(hname, "r2t%02d_%02d_%01d_%02d", lay, iEntr, lr, bin);
00624 m_hr2t[lay][iEntr][lr][bin] = new TH1F(hname, "", 600, -3, 3);
00625 m_fdres2t -> Add(m_hr2t[lay][iEntr][lr][bin]);
00626 }
00627 }
00628 }
00629 }
00630
00631 INTupleSvc* ntupleSvc;
00632 Gaudi::svcLocator() -> service("NTupleSvc", ntupleSvc);
00633 for(lay=0; lay<m_nlayer; lay++){
00634 sprintf(hname, "FILE136/xt%02d", lay);
00635 NTuplePtr nt(ntupleSvc, hname);
00636 if ( nt ) m_xtTuple[lay] = nt;
00637 else{
00638 m_xtTuple[lay] = ntupleSvc->book(hname, CLID_ColumnWiseTuple, "MdcXtNtuple");
00639 if( m_xtTuple[lay] ){
00640 m_xtTuple[lay]->addItem("cel", m_cel[lay]);
00641 m_xtTuple[lay]->addItem("lr", m_lr[lay]);
00642 m_xtTuple[lay]->addItem("run", m_run[lay]);
00643 m_xtTuple[lay]->addItem("evt", m_evt[lay]);
00644 m_xtTuple[lay]->addItem("doca", m_doca[lay]);
00645 m_xtTuple[lay]->addItem("dm", m_dm[lay]);
00646 m_xtTuple[lay]->addItem("tdr", m_tdr[lay]);
00647 m_xtTuple[lay]->addItem("tdc", m_tdc[lay]);
00648 m_xtTuple[lay]->addItem("entr", m_entr[lay]);
00649 m_xtTuple[lay]->addItem("zhit", m_zhit[lay]);
00650 m_xtTuple[lay]->addItem("qhit", m_qhit[lay]);
00651 m_xtTuple[lay]->addItem("p", m_p[lay]);
00652 m_xtTuple[lay]->addItem("pt", m_pt[lay]);
00653 m_xtTuple[lay]->addItem("phi0", m_phi0[lay]);
00654 m_xtTuple[lay]->addItem("tanl", m_tanl[lay]);
00655 m_xtTuple[lay]->addItem("hitphi", m_hitphi[lay]);
00656 } else{
00657 log << MSG::FATAL << "Cannot book Xt N-tuple:"
00658 << long(m_xtTuple[lay]) << endreq;
00659 }
00660 }
00661 }
00662
00663 if(5==m_param.particle){
00664 sprintf(hname, "FILE136/cosmic");
00665 NTuplePtr nt(ntupleSvc, hname);
00666 if ( nt ) m_cosmic = nt;
00667 else{
00668 m_cosmic = ntupleSvc->book(hname, CLID_ColumnWiseTuple, "MdcXtNtuple");
00669 if( m_cosmic ){
00670 m_cosmic->addItem("pUp", m_pUpcos);
00671 m_cosmic->addItem("pDw", m_pDwcos);
00672 m_cosmic->addItem("ptUp", m_ptUpcos);
00673 m_cosmic->addItem("ptDw", m_ptDwcos);
00674 m_cosmic->addItem("phiUp", m_phiUpcos);
00675 m_cosmic->addItem("phiDw", m_phiDwcos);
00676 m_cosmic->addItem("drUp", m_drUpcos);
00677 m_cosmic->addItem("drDw", m_drDwcos);
00678 m_cosmic->addItem("dzUp", m_dzUpcos);
00679 m_cosmic->addItem("dzDw", m_dzDwcos);
00680 m_cosmic->addItem("ctheUp", m_ctheUpcos);
00681 m_cosmic->addItem("ctheDw", m_ctheDwcos);
00682 m_cosmic->addItem("nhitUp", m_nhitUpcos);
00683 m_cosmic->addItem("nhitDw", m_nhitDwcos);
00684 m_cosmic->addItem("char", m_chargecos);
00685 m_cosmic->addItem("tesfg", m_tesFlagcos);
00686 m_cosmic->addItem("tes", m_tescos);
00687 }
00688 }
00689 }
00690 }
00691
00692 int MdcCalib::fillHist(MdcCalEvent* event){
00693 IMessageSvc* msgSvc;
00694 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00695 MsgStream log(msgSvc, "MdcCalib");
00696 log << MSG::DEBUG << "MdcCalib::fillHist()" << endreq;
00697
00698 int i;
00699 int k;
00700 int lay;
00701 int cel;
00702 int wir;
00703 int lr;
00704 int stat;
00705
00706 int hid;
00707 int key;
00708 int iEntr;
00709 int bin;
00710
00711 int phiBin;
00712 int phiBinCms;
00713 int theBin;
00714 int theBinCms;
00715 double lamda;
00716 double theta;
00717
00718 double qhit;
00719 double traw;
00720 double tdr;
00721 double doca;
00722 double resiInc;
00723 double resiExc;
00724 double entr;
00725 double pt;
00726 double p;
00727 double p_cms;
00728 double chisq;
00729 double ecm = m_param.ecm;
00730 double xboost = m_param.boostPar[0] * ecm;
00731 double yboost = m_param.boostPar[1];
00732 double zboost = m_param.boostPar[2];
00733 HepLorentzVector p4;
00734
00735 double dr;
00736 double phi0;
00737 double dz;
00738 double kap;
00739 double tanl;
00740
00741 double x0;
00742 double y0;
00743 double zminus = 9999.0;
00744 double zplus = -9999.0;
00745
00746 double hitphi;
00747 double philab;
00748 double phicms;
00749 double thetacms;
00750 double costheCMS;
00751
00752 double dphi;
00753 double wphi;
00754 double xx;
00755 double yy;
00756 double rr;
00757
00758 int nhitlay;
00759 bool fgHitLay[MdcCalNLayer];
00760 bool fgTrk;
00761
00762 int ntrk = event -> getNTrk();
00763 int nhitRec;
00764 int nhitRecInn;
00765 int nhitRecStp;
00766 int nhitRecOut;
00767 int nhitCal;
00768 int nhitCalInn;
00769 int nhitCalStp;
00770 int nhitCalOut;
00771 MdcCalRecTrk* rectrk;
00772 MdcCalRecHit* rechit;
00773
00774 int fgAdd[43];
00775
00776
00777 if(!fgReadWireEff){
00778 for(lay=0; lay<MdcCalNLayer; lay++){
00779 int ncel = m_mdcGeomSvc->Layer(lay)->NCell();
00780 for(cel=0; cel<ncel; cel++){
00781 double eff = m_mdcFunSvc->getWireEff(lay, cel);
00782 if(eff > 0.5) m_fgGoodWire[lay][cel] = true;
00783 else m_fgGoodWire[lay][cel] = false;
00784 if(eff<0.9) cout << "dead channel: " << setw(5) << lay << setw(5) << cel << endl;
00785 }
00786 }
00787 fgReadWireEff = true;
00788 }
00789
00790 int nRawHit = event->getNRawHitTQ();
00791 m_hnRawHit->Fill(nRawHit);
00792
00793 IDataProviderSvc* eventSvc = NULL;
00794 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
00795
00796 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc,"/Event/EventHeader");
00797 if (!eventHeader) {
00798 log << MSG::FATAL << "Could not find Event Header" << endreq;
00799 return( StatusCode::FAILURE);
00800 }
00801 int iRun = eventHeader->runNumber();
00802 int iEvt = eventHeader->eventNumber();
00803
00804 int esTimeflag = event->getEsFlag();
00805 double tes = event->getTes();
00806 bool esCutFg = event->getEsCutFlag();
00807 int iEs = event->getNesCutFlag();
00808
00809 if (-1 != esTimeflag) {
00810 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc, "/Event/Recon/RecMdcTrackCol");
00811 if(!newtrkCol){
00812 log << MSG::ERROR << "Could not find RecMdcTrackCol" << endreq;
00813 return ( StatusCode::FAILURE );
00814 }
00815 int nGoodTrk = 0;
00816 int icharge = 0;
00817 Vp4 p4p;
00818 Vp4 p4m;
00819 RecMdcTrackCol::iterator it_trk = newtrkCol->begin();
00820 for(; it_trk != newtrkCol->end(); it_trk++){
00821 double mass = 0.000511;
00822 HepLorentzVector p4 = (*it_trk)->p4(mass);
00823 icharge = (*it_trk)->charge();
00824 if (icharge > 0) p4p.push_back(p4);
00825 else p4m.push_back(p4);
00826 }
00827 if (1 == p4p.size() * p4m.size()){
00828 double dang = p4p[0].vect().angle(p4m[0].vect());
00829 m_hbbTrkFlg->Fill(1);
00830 if (dang > 2.94) {
00831 m_hbbTrkFlg->Fill(2);
00832 }
00833 }
00834
00835 }
00836 m_hTesAll->Fill(tes);
00837
00838 if (-1 != esTimeflag) m_hTesGood->Fill(tes);
00839 m_hTesAllFlag->Fill(esTimeflag);
00840 if(ntrk > 0) m_hTesRec->Fill(tes);
00841 if( (iEs >= 0) && (iEs < m_param.nEsFlag) ) m_hTes[iEs]->Fill(tes);
00842 if( esCutFg ) m_hTesCalFlag->Fill(tes);
00843 else return -1;
00844
00845 if(! m_fgIni){
00846 for(lay=0; lay<MdcCalNLayer; lay++){
00847 if(lay < 8) m_docaMax[lay] = m_param.maxDocaInner;
00848 else m_docaMax[lay] = m_param.maxDocaOuter;
00849 }
00850 m_fgIni = true;
00851 }
00852
00853 bool trkFlag[200];
00854 for(i=0; i<200; i++) trkFlag[i] = false;
00855
00856 int ntrkCal = 0;
00857 double pTrk[2];
00858 double pTrkcms[2];
00859 double hitphiplus = 9999.0;
00860 double hitthetaplus = 9999.0;
00861 double hitphiminus = -9999.0;
00862 double hitthetaminus = -9999.0;
00863 Vp4 pp;
00864 Vp4 pm;
00865 pp.clear();
00866 pm.clear();
00867 int phibinp;
00868 int phibinm;
00869
00870 m_hnTrk->Fill(ntrk);
00871 if((ntrk < m_param.nTrkCut[0]) || (ntrk > m_param.nTrkCut[1])){
00872 m_cut1++;
00873 return -1;
00874 }
00875
00876 for(i=0; i<ntrk; i++){
00877 fgTrk = true;
00878 rectrk = event -> getRecTrk(i);
00879 nhitRec = rectrk -> getNHits();
00880 m_hnhitsRec -> Fill( nhitRec );
00881
00882 for(lay=0; lay<MdcCalNLayer; lay++){
00883 fgHitLay[lay] = false;
00884 }
00885
00886 nhitRecInn = 0;
00887 nhitRecStp = 0;
00888 nhitRecOut = 0;
00889 for(k=0; k<nhitRec; k++){
00890 rechit = rectrk -> getRecHit(k);
00891 lay = rechit -> getLayid();
00892 doca = rechit -> getDocaExc();
00893 resiExc = rechit -> getResiExc();
00894 fgHitLay[lay] = true;
00895
00896 if(lay < 8) nhitRecInn++;
00897 else if(lay < 20) nhitRecStp++;
00898 else nhitRecOut++;
00899 }
00900 m_hnhitsRecInn->Fill(nhitRecInn);
00901 m_hnhitsRecStp->Fill(nhitRecStp);
00902 m_hnhitsRecOut->Fill(nhitRecOut);
00903
00904
00905 pt = rectrk -> getPt();
00906 p = rectrk -> getP();
00907
00908
00909 p4 = rectrk->getP4();
00910 HepLorentzVector psip(xboost, yboost, zboost, ecm);
00911 Hep3Vector boostv = psip.boostVector();
00912 p4.boost(- boostv);
00913 p_cms = p4.rho();
00914 phicms = p4.phi();
00915 if(phicms < 0) phicms += PI2;
00916 thetacms = p4.theta();
00917 costheCMS = cos(thetacms);
00918 if (pt < 0) p_cms *= -1.0;
00919 p4.boost(boostv);
00920
00921
00922
00923 if( (costheCMS < m_param.costheCut[0]) || (costheCMS > m_param.costheCut[1]) ){
00924 m_cut2++;
00925 continue;
00926 }
00927
00928
00929 dr = rectrk->getDr();
00930 if(fabs(dr) > m_param.drCut){
00931 m_cut3++;
00932 continue;
00933 }
00934
00935
00936 dz = rectrk->getDz();
00937 if(fabs(dz) > m_param.dzCut){
00938 m_cut4++;
00939 continue;
00940 }
00941
00942
00943
00944
00945 nhitlay = 0;
00946 for(lay=0; lay<MdcCalNLayer; lay++){
00947 if(fgHitLay[lay]) nhitlay++;
00948 }
00949 if(nhitlay < m_param.nHitLayCut){
00950 m_cut5++;
00951 continue;
00952 }
00953
00954
00955 if(nhitRec < m_param.nHitCut){
00956 m_cut6++;
00957 continue;
00958 }
00959
00960
00961
00962
00963
00964
00965
00966
00967
00968
00969 int cellTrkPass[43];
00970 bool fgPass = getCellTrkPass(event, i, cellTrkPass);
00971 for(lay=0; lay<m_nlayer; lay++){
00972 fgAdd[lay] = 0;
00973
00974 if((15==lay) || (16==lay) || (18==lay) || (19==lay) || (40==lay) || (41==lay) || (42==lay)){
00975 int iCell = cellTrkPass[lay];
00976 if(fgPass && (iCell >= 0) && m_fgGoodWire[lay][iCell]) m_effNtrk->Fill(lay);
00977 else fgAdd[lay] = 1;
00978 } else{
00979 m_effNtrk->Fill(lay);
00980 }
00981 }
00982
00983 chisq = rectrk -> getChisq();
00984 m_hchisq -> Fill( chisq );
00985
00986 if(pt > 0){
00987 m_hpt -> Fill(pt);
00988 m_hptPos -> Fill(pt);
00989 m_hp -> Fill(p);
00990 m_hp_cms -> Fill(p_cms);
00991 m_hpPos -> Fill(p);
00992 m_hpPoscms -> Fill(p_cms);
00993 } else{
00994 m_hpt -> Fill(-1.0*pt);
00995 m_hptNeg -> Fill(-1.0*pt);
00996 m_hp -> Fill(-1.0*p);
00997 m_hp_cms -> Fill(-1.0*p_cms);
00998 m_hpNeg -> Fill(-1.0*p);
00999 m_hpNegcms -> Fill(-1.0*p_cms);
01000 }
01001 if(2 == ntrk){
01002 pTrk[i] = fabs(p);
01003 pTrkcms[i] = fabs(p_cms);
01004 }
01005
01006 dr = rectrk -> getDr();
01007 phi0 = rectrk -> getPhi0();
01008 kap = rectrk -> getKappa();
01009 dz = rectrk -> getDz();
01010 tanl = rectrk -> getTanLamda();
01011 lamda = atan(tanl);
01012 theta = HFPI - lamda;
01013
01014 m_hdr -> Fill(dr);
01015 m_hphi0 -> Fill(phi0);
01016 m_hkap -> Fill(kap);
01017 m_hdz -> Fill(dz);
01018 m_htanl -> Fill(tanl);
01019 m_hcosthe -> Fill(cos(theta));
01020 if(pt > 0) m_hcosthePos->Fill(cos(theta));
01021 else m_hcostheNeg->Fill(cos(theta));
01022
01023 philab = phi0 + HFPI;
01024 if(philab > PI2) philab -= PI2;
01025
01026
01027 phiBin = (int)(philab / m_phiWid);
01028 phiBinCms = (int)(phicms / m_phiWid);
01029 theBin = (int)((cos(theta) + 1.0) / m_theWid);
01030 theBinCms = (int)((cos(thetacms) + 1.0) / m_theWid);
01031 if(phiBin < 0) phiBin = 0;
01032 if(phiBin >= NPhiBin) phiBin = NPhiBin-1;
01033 if(phiBinCms < 0) phiBinCms = 0;
01034 if(phiBinCms >= NPhiBin) phiBinCms = NPhiBin-1;
01035 if(theBin < 0) theBin = 0;
01036 if(theBin >= NThetaBin) theBin = NThetaBin-1;
01037 if(theBinCms < 0) theBinCms = 0;
01038 if(theBinCms >= NThetaBin) theBinCms = NThetaBin-1;
01039
01040 if(pt > 0){
01041 m_ppPhi[phiBin]->Fill(p);
01042 m_ppPhiCms[phiBinCms]->Fill(p_cms);
01043 m_ppThe[theBin]->Fill(p);
01044 m_ppTheCms[theBinCms]->Fill(p_cms);
01045 m_ppThePhi[theBin][phiBin]->Fill(p);
01046 m_ppThePhiCms[theBinCms][phiBinCms]->Fill(p_cms);
01047 } else{
01048 m_pnPhi[phiBin]->Fill(-1.0*p);
01049 m_pnPhiCms[phiBinCms]->Fill(-1.0*p_cms);
01050 m_pnThe[theBin]->Fill(-1.0*p);
01051 m_pnTheCms[theBinCms]->Fill(-1.0*p_cms);
01052 m_pnThePhi[theBin][phiBin]->Fill(-1.0*p);
01053 m_pnThePhiCms[theBinCms][phiBinCms]->Fill(-1.0*p_cms);
01054 }
01055
01056 x0 = dr * cos(phi0);
01057 y0 = dr * sin(phi0);
01058 m_hx0 -> Fill(x0);
01059 m_hy0 -> Fill(y0);
01060 if(m_nGrPoint < 10000){
01061 m_grX0Y0->SetPoint(m_nGrPoint, x0, y0);
01062 m_nGrPoint++;
01063 }
01064
01065 if(kap < 0) {
01066 zminus = dz;
01067 pm.push_back(p4);
01068 phibinm = phiBinCms;
01069 } else {
01070 zplus = dz;
01071 pp.push_back(p4);
01072 phibinp = phiBinCms;
01073 }
01074
01075
01076 ntrkCal++;
01077 trkFlag[i] = true;
01078 nhitCal = 0;
01079 nhitCalInn = 0;
01080 nhitCalStp = 0;
01081 nhitCalOut = 0;
01082 for(k=0; k<nhitRec; k++){
01083 rechit = rectrk -> getRecHit(k);
01084
01085 lay = rechit -> getLayid();
01086 cel = rechit -> getCellid();
01087 lr = rechit -> getLR();
01088 stat = rechit -> getStat();
01089 doca = rechit -> getDocaExc();
01090 resiInc = rechit -> getResiIncLR();
01091 resiExc = rechit -> getResiExcLR();
01092 entr = rechit -> getEntra();
01093 tdr = rechit -> getTdrift();
01094 traw = (rechit -> getTdc()) * MdcCalTdcCnv;
01095 wir = m_mdcGeomSvc -> Wire(lay, cel) -> Id();
01096
01097 m_cel[lay] = (long)cel;
01098 m_lr[lay] = (long)lr;
01099 m_run[lay] = iRun;
01100 m_evt[lay] = iEvt;
01101 m_doca[lay] = doca;
01102 m_dm[lay] = rechit -> getDmeas();
01103 m_tdr[lay] = tdr;
01104 m_tdc[lay] = traw;
01105 m_entr[lay] = entr*180.0/3.14;
01106 m_zhit[lay] = rechit -> getZhit();
01107 m_qhit[lay] = rechit -> getQhit();
01108 m_p[lay] = p;
01109 m_pt[lay] = pt;
01110 m_phi0[lay] = phi0;
01111 m_tanl[lay] = tanl;
01112 qhit = rechit -> getQhit();
01113
01114
01115 xx = (m_zhit[lay] - m_zw[wir]) * (m_xe[wir] - m_xw[wir]) /
01116 (m_ze[wir] - m_zw[wir]) + m_xw[wir];
01117 yy = (m_zhit[lay] - m_zw[wir]) * (m_ye[wir] - m_yw[wir]) /
01118 (m_ze[wir] - m_zw[wir]) + m_yw[wir];
01119 rr = sqrt( (xx * xx) + (yy * yy) );
01120 dphi = fabs(doca) / m_radii[lay];
01121
01122 if( yy >= 0 ) wphi = acos(xx / rr);
01123 else wphi = PI2 - acos(xx / rr);
01124 if(1 == lr) hitphi = wphi + dphi;
01125 else hitphi = wphi - dphi;
01126 if(hitphi < 0) hitphi += PI2;
01127 else if(hitphi > PI2) hitphi -= PI2;
01128
01129 m_hitphi[lay] = hitphi;
01130
01131 if( (fabs(doca) > m_docaMax[lay]) ||
01132 (fabs(resiExc) > m_param.resiCut[lay]) ){
01133 continue;
01134 }
01135
01136 if(m_param.fgAdjacLayerCut){
01137 if(0 == lay){
01138 if( ! fgHitLay[1] ) continue;
01139 } else if(42 == lay){
01140 if( ! fgHitLay[41] ) continue;
01141 } else{
01142 if( (!fgHitLay[lay-1]) && (!fgHitLay[lay+1]) ) continue;
01143
01144
01145 if( m_param.fgBoundLayerCut && m_layBound[lay] &&
01146 ((!fgHitLay[lay-1]) || (!fgHitLay[lay+1])) ) continue;
01147 }
01148 }
01149
01150 if((1 == m_param.hitStatCut) && (1 != stat)) continue;
01151
01152
01153 if((1 == m_param.fillNtuple) && (m_nEvtNtuple < m_param.nEvtNtuple)){
01154 m_xtTuple[lay] -> write();
01155 }
01156
01157 if(1 == m_param.hitStatCut){
01158 if( (0 == fgAdd[lay]) && (1 == stat) ){
01159 m_effNtrkRecHit->Fill(lay);
01160 fgAdd[lay] = 1;
01161 }
01162 } else{
01163 if(0 == fgAdd[lay]){
01164 m_effNtrkRecHit->Fill(lay);
01165 fgAdd[lay] = 1;
01166 }
01167 }
01168
01169 nhitCal++;
01170 if(lay < 8) nhitCalInn++;
01171 else if(lay < 20) nhitCalStp++;
01172 else nhitCalOut++;
01173
01174 m_wirehitmap -> Fill(wir);
01175 m_layerhitmap -> Fill(lay);
01176
01177 m_htraw[lay] -> Fill(traw);
01178 m_htdr[lay] -> Fill(tdr);
01179 m_htdrlr[lay][lr]->Fill(tdr);
01180 m_hadc[lay] -> Fill(qhit);
01181
01182 m_hresAllInc -> Fill(resiInc);
01183 m_hresAllExc -> Fill(resiExc);
01184 double resiAve = 0.5 * (resiInc + resiExc);
01185 m_hresAllAve -> Fill(resiAve);
01186
01187 if(lay < 8){
01188 m_hresInnInc -> Fill(resiInc);
01189 m_hresInnExc -> Fill(resiExc);
01190 } else if(lay < 20){
01191 m_hresStpInc -> Fill(resiInc);
01192 m_hresStpExc -> Fill(resiExc);
01193 } else{
01194 m_hresOutInc -> Fill(resiInc);
01195 m_hresOutExc -> Fill(resiExc);
01196 }
01197
01198 int qbin = (int)((qhit-100.0)/100.0);
01199 if(qbin>=0 && qbin<14){
01200 m_hresAveAllQ[qbin]->Fill(resiAve);
01201 m_hresAveLayQ[lay][qbin]->Fill(resiAve);
01202 if(lay > 7) m_hresAveOutQ[qbin]->Fill(resiAve);
01203 }
01204
01205 m_hresInc[lay] -> Fill(resiInc);
01206 m_hreslrInc[lay][lr]->Fill(resiInc);
01207 m_hresExc[lay] -> Fill(resiExc);
01208 m_hreslrExc[lay][lr]->Fill(resiExc);
01209 m_hresAve[lay] -> Fill(resiAve);
01210 m_hreslrAve[lay][lr]->Fill(resiAve);
01211
01212 int iPhi = (int)(hitphi*20.0/PI2);
01213 if(iPhi>=20) iPhi = 19;
01214 m_hresphi[lay][iPhi]->Fill((resiInc+resiExc)*0.5);
01215
01216
01217 bin = (int)(fabs(rechit->getDmeas()) / m_dwid);
01218 iEntr = m_mdcFunSvc -> getSdEntrIndex(entr);
01219 if(1 == m_nEntr[lay]){
01220 iEntr = 0;
01221 } else if(2 == m_nEntr[lay]){
01222 if(entr > 0.0) iEntr = 1;
01223 else iEntr = 0;
01224 }
01225 if((iEntr < MdcCalNENTRSD) && (bin < MdcCalSdNBIN)){
01226 key = getHresId(lay, iEntr, lr, bin);
01227 if( 1 == m_mapr2d.count(key) ){
01228 hid = m_mapr2d[key];
01229 m_hr2dInc[hid] -> Fill(resiInc);
01230 m_hr2dExc[hid] -> Fill(resiExc);
01231 }
01232 }
01233
01234 if((tdr>0) && (tdr<750)){
01235 if(tdr<300) bin = (int)(tdr/10.0);
01236 else bin = (int)((tdr-300.0)/30.0) + 29;
01237 m_hr2t[lay][iEntr][lr][bin]->Fill(resiExc);
01238 }
01239 }
01240 m_nEvtNtuple++;
01241 m_hnhitsCal->Fill(nhitCal);
01242 m_hnhitsCalInn->Fill(nhitCalInn);
01243 m_hnhitsCalStp->Fill(nhitCalStp);
01244 m_hnhitsCalOut->Fill(nhitCalOut);
01245 }
01246 m_hnTrkCal->Fill(ntrkCal);
01247 if(2 == ntrkCal){
01248 if(pTrk[0] > pTrk[1]) m_hpMax->Fill(pTrk[0]);
01249 else m_hpMax->Fill(pTrk[1]);
01250
01251 if(pTrkcms[0] > pTrkcms[1]) m_hpMaxCms->Fill(pTrkcms[0]);
01252 else m_hpMaxCms->Fill(pTrkcms[1]);
01253 }
01254 if(ntrkCal > 0) m_hTesCalUse->Fill(tes);
01255
01256 double delZ0;
01257 if((fabs(zminus) < 9000.0) && (fabs(zplus) < 9000.0)) delZ0 = zplus - zminus;
01258 m_hdelZ0 -> Fill(delZ0);
01259
01260 if (1 == pp.size() * pm.size()){
01261 HepLorentzVector ptot = pp[0] + pm[0];
01262 bool fourmomcut = false;
01263
01264
01265 fourmomcut = (fabs(ptot.x()-0.04)<0.026) && (fabs(ptot.y()) < 0.026)
01266 && (fabs(ptot.z()-0.005)<0.036) && (fabs(ptot.e()-ecm)<0.058);
01267
01268 if (fourmomcut) {
01269 HepLorentzVector psip(xboost, yboost, zboost, ecm);
01270 Hep3Vector boostv = psip.boostVector();
01271 pp[0].boost(- boostv);
01272 pm[0].boost(- boostv);
01273 m_hp_cut->Fill(pp[0].rho());
01274 m_hp_cut->Fill(pm[0].rho());
01275 }
01276 }
01277
01278 if(2==ntrk) for(i=0; i<ntrk; i++) pTrk[i] = (event -> getRecTrk(i)) -> getP();
01279 if((5==m_param.particle) && (2==ntrk) && (fabs(pTrk[0])<5) && (fabs(pTrk[1])<5)){
01280
01281
01282 m_tescos = tes;
01283 m_tesFlagcos = esTimeflag;
01284 for(i=0; i<ntrk; i++){
01285 rectrk = event -> getRecTrk(i);
01286 phi0 = rectrk -> getPhi0();
01287 phi0 = ((phi0+HFPI) > PI2) ? (phi0+HFPI-PI2) : (phi0+HFPI);
01288
01289 tanl = rectrk -> getTanLamda();
01290 lamda = atan(tanl);
01291 theta = HFPI - lamda;
01292
01293 if(phi0 < (2.0*HFPI)){
01294 m_nhitUpcos = rectrk -> getNHits();
01295 m_pUpcos = rectrk -> getP();
01296 m_ptUpcos = rectrk -> getPt();
01297 m_phiUpcos = phi0;
01298 m_drUpcos = rectrk->getDr();
01299 m_dzUpcos = rectrk->getDz();
01300 m_ctheUpcos = cos(theta);
01301 } else{
01302 m_nhitDwcos = rectrk -> getNHits();
01303 m_pDwcos = rectrk -> getP();
01304 m_ptDwcos = rectrk -> getPt();
01305 m_phiDwcos = phi0;
01306 m_drDwcos = rectrk->getDr();
01307 m_dzDwcos = rectrk->getDz();
01308 m_ctheDwcos = cos(theta);
01309
01310 if(m_pDwcos > 0) m_chargecos = 1;
01311 else m_chargecos = 0;
01312 }
01313 }
01314 m_cosmic->write();
01315 }
01316
01317 if(1 == m_param.fgCalDetEffi) calDetEffi();
01318
01319 return 1;
01320 }
01321
01322 int MdcCalib::updateConst(MdcCalibConst* calconst){
01323 IMessageSvc* msgSvc;
01324 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
01325 MsgStream log(msgSvc, "MdcCalib");
01326 log << MSG::DEBUG << "MdcCalib::updateConst()" << endreq;
01327
01328 cout << "Tot " << m_hnTrk->GetEntries()
01329 << ", nTrkCut " << m_cut1 << ", cos(theta)_cut " << m_cut2 << ", drCut " << m_cut3
01330 << ", dzCut " << m_cut4 << ", nHitLayer_cut " << m_cut5 << ", nHit_cut " << m_cut6 << endl;
01331
01332 int lay;
01333 double effi;
01334 double effErr;
01335
01336 int nGoodAll = 0;
01337 int nGoodInn = 0;
01338 int nGoodStp = 0;
01339 int nGoodOut = 0;
01340 int nTotAll = 0;
01341 int nTotInn = 0;
01342 int nTotStp = 0;
01343 int nTotOut = 0;
01344 ofstream feffi("MdcLayerEffi.dat");
01345 for(lay=0; lay<m_nlayer; lay++){
01346 double effNtrk = m_effNtrk->GetBinContent(lay+1);
01347 double effGoodHit = m_effNtrkRecHit->GetBinContent(lay+1);
01348 nGoodAll += effGoodHit;
01349 if(lay < 8) nGoodInn += effGoodHit;
01350 else if (lay < 20) nGoodStp += effGoodHit;
01351 else nGoodOut += effGoodHit;
01352
01353 nTotAll += effNtrk;
01354 if(lay < 8) nTotInn += effNtrk;
01355 else if (lay < 20) nTotStp += effNtrk;
01356 else nTotOut += effNtrk;
01357
01358 effi = (double)effGoodHit / (double)effNtrk;
01359 effErr = sqrt(effi * (1-effi) / (double)effNtrk);
01360 feffi << setw(5) << lay << setw(15) << effi << setw(15) << effErr
01361 << setw(15) << effGoodHit << setw(15) << effNtrk << endl;
01362 }
01363 double effiAll = (double)nGoodAll / (double)(nTotAll);
01364 double errAll = sqrt(effiAll * (1-effiAll) / (double)(nTotAll));
01365 double effiInn = (double)nGoodInn / (double)(nTotInn);
01366 double errInn = sqrt(effiInn * (1-effiInn) / (double)(nTotInn));
01367 double effiStp = (double)nGoodStp / (double)(nTotStp);
01368 double errStp = sqrt(effiStp * (1-effiStp) / (double)(nTotStp));
01369 double effiOut = (double)nGoodOut / (double)(nTotOut);
01370 double errOut = sqrt(effiOut * (1-effiOut) / (double)(nTotOut));
01371 feffi << endl << "EffiAll: " << setw(15) << effiAll << setw(15) << errAll
01372 << setw(15) << nGoodAll << setw(15) << nTotAll << endl;
01373 feffi << endl << "EffiInn: " << setw(15) << effiInn << setw(15) << errInn
01374 << setw(15) << nGoodInn << setw(15) << nTotInn << endl;
01375 feffi << endl << "EffiStp: " << setw(15) << effiStp << setw(15) << errStp
01376 << setw(15) << nGoodStp << setw(15) << nTotStp << endl;
01377 feffi << endl << "EffiOut: " << setw(15) << effiOut << setw(15) << errOut
01378 << setw(15) << nGoodOut << setw(15) << nTotOut << endl;
01379 feffi.close();
01380
01381
01382 if(0 != m_param.fgCalDetEffi){
01383 int nHitAll[] = {0, 0};
01384 int nHitInn[] = {0, 0};
01385 int nHitStp[] = {0, 0};
01386 int nHitOut[] = {0, 0};
01387 ofstream feff2("MdcHitEffi.dat");
01388 for(lay=0; lay<m_nlayer; lay++){
01389 nHitAll[0] += m_hitNum[lay][0];
01390 nHitAll[1] += m_hitNum[lay][1];
01391 if(lay < 8){
01392 nHitInn[0] += m_hitNum[lay][0];
01393 nHitInn[1] += m_hitNum[lay][1];
01394 } else if (lay < 20){
01395 nHitStp[0] += m_hitNum[lay][0];
01396 nHitStp[1] += m_hitNum[lay][1];
01397 } else{
01398 nHitOut[0] += m_hitNum[lay][0];
01399 nHitOut[1] += m_hitNum[lay][1];
01400 }
01401
01402 effi = (double)m_hitNum[lay][1] / (double)m_hitNum[lay][0];
01403 effErr = sqrt(effi * (1-effi) / (double)m_hitNum[lay][0]);
01404 feff2 << setw(5) << lay << setw(15) << effi << setw(15) << effErr
01405 << setw(15) << m_hitNum[lay][1] << setw(15) << m_hitNum[lay][0] << endl;
01406 }
01407 effiAll = (double)nHitAll[1] / (double)nHitAll[0];
01408 errAll = sqrt(effiAll * (1-effiAll)) / (double)nHitAll[0];
01409 effiInn = (double)nHitInn[1] / (double)nHitInn[0];
01410 errInn = sqrt(effiInn * (1-effiInn)) / (double)nHitInn[0];
01411 effiStp = (double)nHitStp[1] / (double)nHitStp[0];
01412 errStp = sqrt(effiStp * (1-effiStp)) / (double)nHitStp[0];
01413 effiOut = (double)nHitOut[1] / (double)nHitOut[0];
01414 errOut = sqrt(effiOut * (1-effiOut)) / (double)nHitOut[0];
01415 feff2 << endl << "EffiAll: " << setw(15) << effiAll << setw(15) << errAll
01416 << setw(15) << nHitAll[1] << setw(15) << nHitAll[0] << endl;
01417 feff2 << endl << "EffiInn: " << setw(15) << effiInn << setw(15) << errInn
01418 << setw(15) << nHitInn[1] << setw(15) << nHitInn[0] << endl;
01419 feff2 << endl << "EffiStp: " << setw(15) << effiStp << setw(15) << errStp
01420 << setw(15) << nHitStp[1] << setw(15) << nHitStp[0] << endl;
01421 feff2 << endl << "EffiOut: " << setw(15) << effiOut << setw(15) << errOut
01422 << setw(15) << nHitOut[1] << setw(15) << nHitOut[0] << endl;
01423 feff2.close();
01424 }
01425
01426
01427 int i;
01428 int iEntr;
01429 int lr;
01430 int bin;
01431 int key;
01432 int hid;
01433
01434 Stat_t entry;
01435 double sigm[MdcCalSdNBIN];
01436 if(m_param.calSigma){
01437 ofstream fr2d("logr2d.dat");
01438 for(lay=0; lay<m_nlayer; lay++){
01439 for(iEntr=0; iEntr<m_nEntr[lay]; iEntr++){
01440 for(lr=0; lr<2; lr++){
01441 fr2d << setw(3) << lay << setw(3) << iEntr << setw(3) << lr << endl;
01442 for(bin=0; bin<m_nBin[lay]; bin++){
01443 key = getHresId(lay, iEntr, lr, bin);
01444 hid = m_mapr2d[key];
01445
01446 if(1 == m_param.resiType){
01447 entry = m_hr2dExc[hid] -> GetEntries();
01448 if(entry > 500){
01449 m_hr2dExc[hid] -> Fit("gaus", "Q");
01450 sigm[bin] = m_hr2dExc[hid]->GetFunction("gaus")->GetParameter(2);
01451 } else if(entry > 100){
01452 sigm[bin] = m_hr2dExc[hid]->GetRMS();
01453 } else{
01454 sigm[bin] = 0.2;
01455 }
01456 } else{
01457 entry = m_hr2dInc[hid] -> GetEntries();
01458 if(entry > 500){
01459 m_hr2dInc[hid] -> Fit("gaus", "Q");
01460 sigm[bin] = m_hr2dInc[hid]->GetFunction("gaus")->GetParameter(2);
01461 } else if(entry > 100){
01462 sigm[bin] = m_hr2dInc[hid]->GetRMS();
01463 } else{
01464 sigm[bin] = 0.2;
01465 }
01466 }
01467 if(sigm[bin] < 0.05) sigm[bin] = 0.05;
01468 }
01469
01470 for(bin=m_nBin[lay]; bin<MdcCalSdNBIN; bin++){
01471 sigm[bin] = sigm[m_nBin[lay]-1];
01472 }
01473
01474 for(bin=0; bin<MdcCalSdNBIN; bin++){
01475 if(1 == m_param.fgCalib[lay]){
01476
01477 if(1 == m_nEntr[lay]){
01478 for(i=0; i<6; i++) calconst -> resetSdpar(lay, i, lr, bin, sigm[bin]);
01479 } else if(2 == m_nEntr[lay]){
01480 if(0 == iEntr){
01481 for(i=0; i<3; i++){
01482 calconst -> resetSdpar(lay, i, lr, bin, sigm[bin]);
01483 }
01484 } else{
01485 for(i=3; i<6; i++){
01486 calconst -> resetSdpar(lay, i, lr, bin, sigm[bin]);
01487 }
01488 }
01489 }
01490 } else{
01491 sigm[bin] = calconst->getSdpar(lay, iEntr, lr, bin);
01492 }
01493 fr2d << setw(5) << bin << setw(15) << sigm[bin] << endl;
01494 }
01495 }
01496 }
01497 }
01498 fr2d.close();
01499 }
01500
01501 return 1;
01502 }
01503
01504 int MdcCalib::getHresId(int lay, int entr, int lr, int bin) const{
01505 int index = ( (lay << HRESLAYER_INDEX) & HRESLAYER_MASK ) |
01506 ( (entr << HRESENTRA_INDEX) & HRESENTRA_MASK ) |
01507 ( (lr << HRESLR_INDEX) & HRESLR_MASK ) |
01508 ( (bin << HRESBIN_INDEX) & HRESBIN_MASK );
01509 return index;
01510 }
01511
01512 int MdcCalib::calDetEffi(){
01513 IMessageSvc* msgSvc;
01514 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
01515 MsgStream log(msgSvc, "MdcCalib");
01516
01517 IDataProviderSvc* eventSvc = NULL;
01518 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
01519 SmartDataPtr<MdcDigiCol> mdcDigiCol(eventSvc,"/Event/Digi/MdcDigiCol");
01520 if (!mdcDigiCol) {
01521 log << MSG::FATAL << "Could not find event" << endreq;
01522 }
01523
01524 int lay;
01525 int cel;
01526 bool hitCel[43][288];
01527 int hitCel2[43][288];
01528 for(lay=0; lay<43; lay++){
01529 for(cel=0; cel<288; cel++){
01530 hitCel[lay][cel] = false;
01531 hitCel2[lay][cel] = 0;
01532 }
01533 }
01534
01535 MdcDigiCol::iterator iter = mdcDigiCol->begin();
01536 unsigned fgOverFlow;
01537 int digiId = 0;
01538 Identifier id;
01539 for(; iter != mdcDigiCol->end(); iter++, digiId++) {
01540 MdcDigi *aDigi = (*iter);
01541 id = (aDigi)->identify();
01542
01543 lay = MdcID::layer(id);
01544 cel = MdcID::wire(id);
01545 fgOverFlow = (aDigi) -> getOverflow();
01546
01547
01548
01549
01550 if ( ((fgOverFlow & 3) !=0 ) || ((fgOverFlow & 12) != 0) ||
01551 (aDigi->getTimeChannel() == 0x7FFFFFFF) ||
01552 (aDigi->getChargeChannel() == 0x7FFFFFFF) ) continue;
01553 hitCel[lay][cel] = true;
01554 hitCel2[lay][cel] = 1;
01555 }
01556
01557 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc, "/Event/Recon/RecMdcTrackCol");
01558 if(!newtrkCol){
01559 log << MSG::ERROR << "Could not find RecMdcTrackCol" << endreq;
01560 return -1;
01561 }
01562
01563 double dphi = 1.0;
01564 Identifier identifier;
01565 MdcID mdcid;
01566 RecMdcTrackCol::iterator it_trk = newtrkCol->begin();
01567 for(it_trk = newtrkCol->begin(); it_trk != newtrkCol->end(); it_trk++){
01568 HitRefVec gothits = (*it_trk) -> getVecHits();
01569 HitRefVec::iterator it_hit = gothits.begin();
01570 for(; it_hit != gothits.end(); it_hit++){
01571 identifier = (*it_hit)->getMdcId();
01572 lay = mdcid.layer(identifier);
01573 cel = mdcid.wire(identifier);
01574 hitCel2[lay][cel] = 2;
01575 }
01576 }
01577 for(it_trk = newtrkCol->begin(); it_trk != newtrkCol->end(); it_trk++){
01578 HepVector helix = (*it_trk)->helix();
01579 HepSymMatrix helixErr = (*it_trk)->err();
01580 double phi0 = (*it_trk)->helix(1);
01581 double phiTrk = phi0 + HFPI;
01582 if(phiTrk > PI2) phiTrk -= PI2;
01583
01584 for(lay=0; lay<43; lay++){
01585 double docamin = 0.9;
01586 if(lay<8) docamin = 0.6;
01587 int celmin = -1;
01588 int ncel = m_mdcGeomSvc->Layer(lay)->NCell();
01589 for(cel=0; cel<ncel; cel++){
01590 double wphi;
01591 const MdcGeoWire* pWire = m_mdcGeomSvc -> Wire(lay, cel);
01592 double xx = pWire->Backward().x();
01593 double yy = pWire->Backward().y();
01594 double rr = sqrt( (xx * xx) + (yy * yy) );
01595 if( yy >= 0 ) wphi = acos(xx / rr);
01596 else wphi = CLHEP::twopi - acos(xx / rr);
01597
01598 if( !( (fabs(wphi-phiTrk) < dphi) || (fabs(wphi+PI2-phiTrk) < dphi) ||
01599 (fabs(wphi-PI2-phiTrk) < dphi) ) ){
01600 continue;
01601 }
01602
01603 double doca = m_mdcUtilitySvc->doca(lay, cel, helix, helixErr);
01604
01605 if(fabs(doca) < fabs(docamin)){
01606 docamin = doca;
01607 celmin = cel;
01608 }
01609 }
01610 if(celmin > -1){
01611 int wir = m_mdcGeomSvc -> Wire(lay, celmin) -> Id();
01612 m_hitEffAll->Fill(wir);
01613 m_hitEffAll->Fill(6799);
01614 if(lay<8) m_hitEffAll->Fill(6796);
01615 else if(lay<20) m_hitEffAll->Fill(6797);
01616 else m_hitEffAll->Fill(6798);
01617
01618 if(hitCel[lay][celmin]){
01619 m_hitEffRaw->Fill(wir);
01620 m_hitEffRaw->Fill(6799);
01621 if(lay<8) m_hitEffRaw->Fill(6796);
01622 else if(lay<20) m_hitEffRaw->Fill(6797);
01623 else m_hitEffRaw->Fill(6798);
01624 }
01625 if(2==hitCel2[lay][celmin]){
01626 m_hitEffRec->Fill(wir);
01627 m_hitEffRec->Fill(6799);
01628 if(lay<8) m_hitEffRec->Fill(6796);
01629 else if(lay<20) m_hitEffRec->Fill(6797);
01630 else m_hitEffRec->Fill(6798);
01631 }
01632 }
01633 }
01634 }
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656
01657
01658
01659
01660
01661
01662
01663
01664
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674
01675
01676
01677
01678
01679
01680
01681
01682
01683
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707 return 1;
01708 }
01709
01710 bool MdcCalib::getCellTrkPass(MdcCalEvent* event, int iTrk, int cellTrkPass[]){
01711 MdcCalRecTrk* rectrk = event -> getRecTrk(iTrk);
01712 int nHits = rectrk -> getNHits();
01713 int hitInTrk[100];
01714 for(int k=0; k<nHits; k++){
01715 MdcCalRecHit* rechit = rectrk->getRecHit(k);
01716 int lay = rechit->getLayid();
01717 int cel = rechit->getCellid();
01718 int wir = m_mdcGeomSvc -> Wire(lay, cel) -> Id();
01719 hitInTrk[k] = wir;
01720 }
01721
01722 IDataProviderSvc* eventSvc = NULL;
01723 Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
01724
01725 SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc, "/Event/Recon/RecMdcTrackCol");
01726 if(!newtrkCol){
01727
01728 return false;
01729 }
01730 MdcID mdcid;
01731 Identifier identifier;
01732 double dphi = 1.0;
01733 RecMdcTrackCol::iterator it_trk = newtrkCol->begin();
01734 for(it_trk = newtrkCol->begin(); it_trk != newtrkCol->end(); it_trk++){
01735 int nRecHits = (*it_trk)->getNhits();
01736 if(nRecHits < nHits) continue;
01737
01738 int hitInRecTrk[100];
01739 int iRecHit = 0;
01740 HitRefVec gothits = (*it_trk) -> getVecHits();
01741 HitRefVec::iterator it_hit = gothits.begin();
01742 for(; it_hit != gothits.end(); it_hit++){
01743 identifier = (*it_hit)->getMdcId();
01744 int lay = mdcid.layer(identifier);
01745 int cel = mdcid.wire(identifier);
01746 int wir = m_mdcGeomSvc -> Wire(lay, cel) -> Id();
01747 hitInRecTrk[iRecHit] = wir;
01748 iRecHit++;
01749 }
01750
01751
01752 bool matchSuccess = true;
01753 for(int i=0; i<nHits; i++){
01754 bool findHit = false;
01755 for(int k=0; k<nRecHits; k++){
01756 if(hitInTrk[i] == hitInRecTrk[k]){
01757 findHit = true;
01758 break;
01759 }
01760 }
01761 if(!findHit){
01762 matchSuccess = false;
01763 break;
01764 }
01765 }
01766 if(!matchSuccess) continue;
01767
01768 HepVector helix = (*it_trk)->helix();
01769 HepSymMatrix helixErr = (*it_trk)->err();
01770 double phi0 = (*it_trk)->helix(1);
01771 double phiTrk = phi0 + HFPI;
01772 if(phiTrk > PI2) phiTrk -= PI2;
01773
01774 for(int lay=0; lay<43; lay++){
01775 double docamin = 0.9;
01776 if(lay<8) docamin = 0.6;
01777 int celmin = -1;
01778 int ncel = m_mdcGeomSvc->Layer(lay)->NCell();
01779 for(int cel=0; cel<ncel; cel++){
01780 double wphi;
01781 const MdcGeoWire* pWire = m_mdcGeomSvc -> Wire(lay, cel);
01782 double xx = pWire->Backward().x();
01783 double yy = pWire->Backward().y();
01784 double rr = sqrt( (xx * xx) + (yy * yy) );
01785 if( yy >= 0 ) wphi = acos(xx / rr);
01786 else wphi = CLHEP::twopi - acos(xx / rr);
01787
01788 if( !( (fabs(wphi-phiTrk) < dphi) || (fabs(wphi+PI2-phiTrk) < dphi) ||
01789 (fabs(wphi-PI2-phiTrk) < dphi) ) ){
01790 continue;
01791 }
01792
01793 double doca = m_mdcUtilitySvc->doca(lay, cel, helix, helixErr);
01794
01795 if(fabs(doca) < fabs(docamin)){
01796 docamin = doca;
01797 celmin = cel;
01798 }
01799 }
01800 if(celmin > -1){
01801 cellTrkPass[lay] = celmin;
01802 } else{
01803 cellTrkPass[lay] = -1;
01804 }
01805 }
01806 return true;
01807 }
01808 return false;
01809 }