CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
DotsConnection.cxx
Go to the documentation of this file.
1
2#include "DotsConnection/DotsConnection.h"
3
4#include "GaudiKernel/MsgStream.h"
5//#include "GaudiKernel/ISvcLocator.h"
6#include "GaudiKernel/SmartDataPtr.h"
7#include "GaudiKernel/IPartPropSvc.h"
8#include "GaudiKernel/INTupleSvc.h"
9//#include "GaudiKernel/IDataManagerSvc.h"
10//#include "GaudiKernel/IDataProviderSvc.h"
11#include "HepPDT/ParticleDataTable.hh"
12#include "McTruth/McParticle.h"
13#include "McTruth/MdcMcHit.h"
14#include "EventModel/EventHeader.h"
15#include "EvTimeEvent/RecEsTime.h"
16#include "Identifier/MdcID.h"
17#include "RawEvent/RawDataUtil.h"
18// #include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
19
20using namespace Event; // for McParticleCol
21//using namespace EventModel; // for EventHeader
22
23DotsConnection::DotsConnection(const std::string& name, ISvcLocator* pSvcLocator) :
24 Algorithm(name, pSvcLocator)
25{
26 myMdcCalibFunSvc=NULL;
27
28 // Part 1: Declare the properties
29 //declareProperty("MyInt", m_myInt);
30 declareProperty("Ntuple", myNtProd=0);
31 declareProperty("driftTimeUpLimit", myDriftTimeUpLimit = 400);
32 declareProperty("MdcHitChi2Cut", myMdcHitChi2Cut = 100);
33 declareProperty("ChiCut_circle", myChiCut_circle=5);
34 declareProperty("NmaxDeact_circle", myNmaxDeact_circle=1);
35 declareProperty("ChiCut_helix", myChiCut_helix=5);
36 declareProperty("NmaxDeact_helix", myNmaxDeact_helix=1);
37 declareProperty("Debug", myDebug=0);
38 declareProperty("Chi2CutDiverge", myChi2CutDiverge=99999999);
39 //declareProperty("MCparID", myMCparID=1);//1: e, 2: mu, 3: pi, 4: K, 5: p
40}
41
42// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
43
44StatusCode DotsConnection::initialize(){
45
46 // Part 1: Get the messaging service, print where you are
47 MsgStream log(msgSvc(), name());
48 log << MSG::INFO << " DotsConnection initialize()" << endreq;
49
50 // Part 2: Print out the property values
51 // log << MSG::INFO << " MyInt = " << m_myInt << endreq;
52
53 // --- Get MdcGeomSvc ---
54 IMdcGeomSvc* imdcGeomSvc;
55 //StatusCode sc1 =Gaudi::svcLocator()->service("MdcGeomSvc", imdcGeomSvc);
56 StatusCode sc = service ("MdcGeomSvc", imdcGeomSvc);
57 myMdcGeomSvc = dynamic_cast<MdcGeomSvc*> (imdcGeomSvc);
58 if (sc.isFailure()) {
59 return( StatusCode::FAILURE);
60 }
61 int nWire = myMdcGeomSvc->getWireSize();
62 int nLayer= myMdcGeomSvc->getLayerSize();
63 int nSuperLayer= myMdcGeomSvc->getSuperLayerSize();
64 int nSeg= myMdcGeomSvc->getSegmentNo();
65 /*cout<<"nWire = "<<nWire
66 <<" nLayer="<<nLayer
67 <<" nSuperLayer="<<nSuperLayer
68 <<" nSeg="<<nSeg
69 <<endl;*/
70 int nCellTot=0;
71 for(int i=0; i<nLayer; i++)
72 {
73 const MdcGeoLayer* aLayer = myMdcGeomSvc->Layer(i);
74 myNWire[i]=aLayer->NCell();
75 myRLayer[i]=aLayer->Radius();
76 //cout<<"layer "<<i<<", "<<aLayer->NCell()<<" cells, slant "<<aLayer->Slant()<<", R="<<aLayer->Radius()<<endl;
77
78 double nomShift = aLayer->nomShift();
79 if(nomShift>0) myWireFlag[i]=1;
80 else if(nomShift<0) myWireFlag[i]=-1;
81 else myWireFlag[i]=0;
82
83 //nCellTot+=aLayer->NCell();
84 nCellTot+=myNWire[i];
85 for(int j=0; j<myNWire[i]; j++)
86 {
87 const MdcGeoWire* aWire = myMdcGeomSvc->Wire(i,j);
88 //cout<<" wire "<<j<<", BWpos ="<<aWire->Backward()
89 // <<", "<<aWire->BWirePos()
90 // <<", FWpos ="<<aWire->Forward()
91 // <<", "<<aWire->FWirePos()
92 // <<endl;
93 myWirePhi[i][j]=aWire->Forward().phi();
94 int iInnerLayer = i-1;
95 if(iInnerLayer>=0) {
96 for(int k=0; k<myNWire[iInnerLayer]; k++)
97 {
98 int k_last = k-1;
99 if(k_last<0) k_last=myNWire[iInnerLayer]-1;
100 double dphi_last = dPhi(myWirePhi[iInnerLayer][k_last], myWirePhi[i][j]);
101 double dphi = dPhi(myWirePhi[iInnerLayer][k], myWirePhi[i][j]);
102 if(dphi_last<0&&dphi>0) {
103 myInnerWire[i][j][0]=k_last;
104 myInnerWire[i][j][1]=k;
105 //cout<<"k_last, k ="<<k_last<<", "<<k<<endl;
106 break;
107 }
108 }
109 }
110 }
111 }
112 for(int i=0; i<nLayer; i++)
113 {
114 for(int j=0; j<myNWire[i]; j++)
115 {
116 int iOuterLayer = i+1;
117 if(iOuterLayer<nLayer) {
118 for(int k=0; k<myNWire[iOuterLayer]; k++)
119 {
120 int k_last = k-1;
121 if(k_last<0) k_last=myNWire[iOuterLayer]-1;
122 double dphi_last = dPhi(myWirePhi[iOuterLayer][k_last], myWirePhi[i][j]);
123 double dphi = dPhi(myWirePhi[iOuterLayer][k], myWirePhi[i][j]);
124 if(dphi_last<0&&dphi>0) {
125 myOuterWire[i][j][0]=k_last;
126 myOuterWire[i][j][1]=k;
127 //cout<<"k_last, k ="<<k_last<<", "<<k<<endl;
128 //cout<<"phi="<<myWirePhi[i][j]<<", outer phi = "<<myWirePhi[iOuterLayer][k_last]<<", "<<myWirePhi[iOuterLayer][k]<<endl;
129 break;
130 }
131 }
132 }
133 }
134 }
135 //cout<<"Total "<<nCellTot<<" cells"<<endl;
136 // ----------------------
137
138
139 // --- Initialize RawDataProviderSvc ---
140 IRawDataProviderSvc* irawDataProviderSvc;
141 sc = service ("RawDataProviderSvc", irawDataProviderSvc);
142 myRawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
143 if ( sc.isFailure() ){
144 log << MSG::FATAL << name()<<" Could not load RawDataProviderSvc!" << endreq;
145 return StatusCode::FAILURE;
146 }
147
148 // --- MdcCalibFunSvc ---
149 IMdcCalibFunSvc* imdcCalibSvc;
150 sc = service("MdcCalibFunSvc", imdcCalibSvc);
151 if ( sc.isSuccess() ){
152 myMdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
153 }
154 else {
155 cout<<"DotsConnection::initialize(): can not get MdcCalibFunSvc"<<endl;
156 }
157 // --- initialize DotsHelixFitter ---
158 myDotsHelixFitter.initialize();
159 myDotsHelixFitter.setChi2Diverge(myChi2CutDiverge);
160 if(myNtProd&1)
161 {
162 NTuplePtr nt_ptr(ntupleSvc(),"TestDotsHelixFitter/trkPar");
163 if( nt_ptr ) myNtHelixFitter = nt_ptr;
164 else
165 {
166 myNtHelixFitter = ntupleSvc()->book("TestDotsHelixFitter/trkPar",CLID_ColumnWiseTuple,"trkPar");
167 if( myNtHelixFitter ) {
168 myNtHelixFitter->addItem ("run", myRUN);
169 myNtHelixFitter->addItem ("evt", myEVT);
170 myNtHelixFitter->addItem ("pid", myPID);
171 myNtHelixFitter->addItem ("Npar", myNPar, 0,5);
172 myNtHelixFitter->addIndexedItem("HelixMC", myNPar, myArrayHelixMC);
173 myNtHelixFitter->addIndexedItem("HelixFitted", myNPar, myArrayHelixFitted);
174 myNtHelixFitter->addItem ("NhitCircle", myNHitsCircle, 0,100);
175 myNtHelixFitter->addIndexedItem("LayerHitsCircle", myNHitsCircle, myLayerHitsCircle);
176 myNtHelixFitter->addIndexedItem("ChiHitsCircle", myNHitsCircle, myChiHitsCircle);
177 myNtHelixFitter->addItem ("Nhit", myNHits, 0,100);
178 myNtHelixFitter->addIndexedItem("LayerHits", myNHits, myLayerHits);
179 myNtHelixFitter->addIndexedItem("ChiHits", myNHits, myChiHits);
180 myNtHelixFitter->addItem ("NXhit", myNXHits, 0,100);
181 myNtHelixFitter->addItem ("NVhit", myNVHits, 0,100);
182 myNtHelixFitter->addItem ("NXCluster", myNXCluster, 0,100);
183 myNtHelixFitter->addItem ("NVCluster", myNVCluster, 0,100);
184 myNtHelixFitter->addItem ("TrkIdBest", myTrkIdBest);
185 myNtHelixFitter->addItem ("NHitsBestTrk", myNHitsBestTrk);
186 myNtHelixFitter->addItem ("NSameHitsBestTrk", myNSameHitsBestTrk);
187 //cout<<"myNtHelixFitter added !"<<endl;
188 }
189 }
190 }
191
192 // --- get MdcUtilitySvc
193 /*
194 sc = service("MdcUtilitySvc", myMdcUtilitySvc);
195 if( sc != StatusCode::SUCCESS ){
196 log << MSG::FATAL << "can not use MdcUtilitySvc" << endreq;
197 }*/
198
199
200 return StatusCode::SUCCESS;
201}
202
203// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
204StatusCode DotsConnection::execute() {
205
206 // Part 1: Get the messaging service, print where you are
207 MsgStream log(msgSvc(), name());
208 log << MSG::INFO << "DotsConnection execute()" << endreq;
209
210 // Part 2: Print out the different levels of messages
211 /*log << MSG::DEBUG << "A DEBUG message" << endreq;
212 log << MSG::INFO << "An INFO message" << endreq;
213 log << MSG::WARNING << "A WARNING message" << endreq;
214 log << MSG::ERROR << "An ERROR message" << endreq;
215 log << MSG::FATAL << "A FATAL error message" << endreq;
216 */
217
218 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
219 if (!eventHeader) {
220 log << MSG::WARNING << "Could not find Event Header" << endreq;
221 return StatusCode::FAILURE;
222 }
223 myRun = eventHeader->runNumber();
224 myEvt = eventHeader->eventNumber();
225 if(myDebug) {
226 cout<<endl<<"------------------------------ "<<endl;
227 //cout<<"run, evt = "<<myEvt<<", "<<myRun<<endl;
228 cout<<"run:"<<myRun<<" , event: "<<myEvt<<endl;
229 }
230
231 //if(myRun==763)
232 //testDotsHelixFitterAllHits();
233
234 //testDotsHelixFitterPartHits();
235
236
237 // --- ideal tracking
238 getMcFinalChargedStates();
239 bool bookTrkCol = registerRecMdcTrack();
240 if(!bookTrkCol)
241 {
242 cout<<"DotsConnection::execute(): failed to register RecMdcTrackCol!"<<endl;
243 return StatusCode::FAILURE;
244 }
245 associateDigisToMcParticles();
246
247
248
249 return StatusCode::SUCCESS;
250}
251
252// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
253
254StatusCode DotsConnection::finalize() {
255
256 // Part 1: Get the messaging service, print where you are
257 MsgStream log(msgSvc(), name());
258 log << MSG::INFO << "DotsConnection finalize()" << endreq;
259
260 return StatusCode::SUCCESS;
261}
262
263
264double DotsConnection::dPhi(double phi1, double phi2)
265{
266 double dphi=phi1-phi2;
267 while(dphi>M_PI) dphi-=M_PI*2.0;
268 while(dphi<-M_PI) dphi+=M_PI*2.0;
269 return dphi;
270}
271
272KalmanFit::Helix DotsConnection::getMCHelix()
273{
274
275 // particle information
276 std::string name;
277 double charge = 0;
278 double pos_x, pos_y, pos_z;
279 double px, py, pz;
280
281
282 // get PartPropSvc
283 IPartPropSvc* p_PartPropSvc;
284 static const bool CREATEIFNOTTHERE(true);
285 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
286 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
287 cout<< " Could not initialize Particle Properties Service" << endl;
288 }
289 HepPDT::ParticleDataTable* p_particleTable = p_PartPropSvc->PDT();
290
291 // get MC particle collection
292 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
293 if (mcParticleCol) {
294 McParticleCol::iterator iter_mc = mcParticleCol->begin();
295 for (;iter_mc != mcParticleCol->end(); iter_mc++ ) {
296 if(!(*iter_mc)->primaryParticle()) continue;
297 int pid = (*iter_mc)->particleProperty();
298 int pid_abs = abs(pid);
299 if( p_particleTable->particle(pid) ){
300 name = p_particleTable->particle(pid)->name();
301 charge= p_particleTable->particle(pid)->charge();
302 }else if( p_particleTable->particle(-pid) ){
303 name = "anti " + p_particleTable->particle(-pid)->name();
304 charge = (-1.)*p_particleTable->particle(-pid)->charge();
305 }
306 pos_x = (*iter_mc)->initialPosition().x();
307 pos_y = (*iter_mc)->initialPosition().y();
308 pos_z = (*iter_mc)->initialPosition().z();
309 px = (*iter_mc)->initialFourMomentum().px();
310 py = (*iter_mc)->initialFourMomentum().py();
311 pz = (*iter_mc)->initialFourMomentum().pz();
312 if(pid_abs==11||pid_abs==13||pid_abs==211||pid_abs==321||pid_abs==2212) break;
313 }
314 }
315
316 // --- get Helix ---
317 HepPoint3D posTruth(pos_x, pos_y, pos_z);
318 Hep3Vector p3Truth(px, py, pz);
319 KalmanFit::Helix helix_truth(posTruth,p3Truth,charge);
320 HepPoint3D origin(0, 0, 0);
321 helix_truth.pivot(origin);
322 if(myDebug) {
323 cout<<"DotsConnection::getMCHelix() finds a valid MC particle "<<name<<" at "<<posTruth<<" with p3="<<p3Truth<<endl;
324 cout<<"DotsConnection::getMCHelix() Helix = "<<helix_truth.a()<<endl;
325 }
326
327 return helix_truth;
328}
329
330void DotsConnection::getMcFinalChargedStates()
331{
332 resetFCSVec();
333
334 // particle information
335 std::string name;
336 double charge = 0;
337 double pos_x, pos_y, pos_z;
338 double px, py, pz;
339 HepPoint3D origin(0, 0, 0);
340
341
342 // get PartPropSvc
343 IPartPropSvc* p_PartPropSvc;
344 static const bool CREATEIFNOTTHERE(true);
345 StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
346 if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
347 cout<< " Could not initialize Particle Properties Service" << endl;
348 }
349 HepPDT::ParticleDataTable* p_particleTable = p_PartPropSvc->PDT();
350
351 // get MC particle collection
352 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
353 if (mcParticleCol) {
354 McParticleCol::iterator iter_mc = mcParticleCol->begin();
355 if(myDebug) cout<<"MC charged particles: "<<endl;
356 for (;iter_mc != mcParticleCol->end(); iter_mc++ )
357 {
358 int pdgid = (*iter_mc)->particleProperty();
359 if(myDebug) cout<<"idx, pdg, from "<<(*iter_mc)->trackIndex()<<", "<<pdgid<<", "<<((*iter_mc)->mother()).trackIndex()<<endl;
360 if( !( (*iter_mc)->primaryParticle() || (*iter_mc)->decayFromGenerator() || (*iter_mc)->decayInFlight() ) ) continue;
361 int pid = (*iter_mc)->particleProperty();
362 int pid_abs = abs(pid);
363 // --- e, mu, pi, K, p
364 if(pid_abs==11||pid_abs==13||pid_abs==211||pid_abs==321||pid_abs==2212)
365 {
366 if( p_particleTable->particle(pid) )
367 {
368 name = p_particleTable->particle(pid)->name();
369 charge= p_particleTable->particle(pid)->charge();
370 }
371 else if( p_particleTable->particle(-pid) )
372 {
373 name = "anti " + p_particleTable->particle(-pid)->name();
374 charge = (-1.)*p_particleTable->particle(-pid)->charge();
375 }
376 pos_x = (*iter_mc)->initialPosition().x();
377 pos_y = (*iter_mc)->initialPosition().y();
378 pos_z = (*iter_mc)->initialPosition().z();
379 px = (*iter_mc)->initialFourMomentum().px();
380 py = (*iter_mc)->initialFourMomentum().py();
381 pz = (*iter_mc)->initialFourMomentum().pz();
382 // --- get Helix ---
383 HepPoint3D posTruth(pos_x, pos_y, pos_z);
384 Hep3Vector p3Truth(px, py, pz);
385 KalmanFit::Helix helix_truth(posTruth,p3Truth,charge);
386 helix_truth.pivot(origin);
387 // --- get the flight length of the first half in MDC --
388 double dphi=helix_truth.IntersectCylinder(myRLayer[42]/10.);// mm -> cm
389 if(dphi==0) dphi=M_PI*charge;
390 HepPoint3D posOuter = helix_truth.x(dphi);
391 double lenOuter = helix_truth.flightLength(posOuter);
392 double lenInner = helix_truth.flightLength(posTruth);
393 double lenInMdc = lenOuter-lenInner;
394 myVecTrkLenFirstHalf.push_back(lenInMdc);
395 // --- fill charged final states
396 int trkIdx = (*iter_mc)->trackIndex();
397 myVecMCTrkId.push_back(trkIdx);
398 myVecPDG.push_back(pid);
399 vector<double> helix;
400 helix.push_back(helix_truth.dr());
401 helix.push_back(helix_truth.phi0());
402 helix.push_back(helix_truth.kappa());
403 helix.push_back(helix_truth.dz());
404 helix.push_back(helix_truth.tanl());
405 myVecHelix.push_back(helix);
406 if(myDebug)
407 cout<<" trk idx "<<trkIdx<<", pdg code "<<pid
408 <<", p3=("<<px<<", "<<py<<", "<<pz<<")"
409 <<", pos=("<<pos_x<<","<<pos_y<<","<<pos_z<<")"
410 <<", dphi="<<dphi<<", lenOuter="<<lenOuter<<", lenInner="<<lenInner
411 <<", lenInMdc = "<<lenInMdc
412 //<<", lenInMdc = "<<lenOuter
413 <<endl;
414 }
415 }
416 }
417}
418
419
420void DotsConnection::resetFCSVec()
421{
422 myVecMCTrkId.clear();
423 myVecPDG.clear();
424 myVecTrkLenFirstHalf.clear();
425 myVecHelix.clear();
426 myVecCgemMcHit.clear();
427 myVecCgemXcluster.clear();
428 myVecCgemVcluster.clear();
429}
430
431void DotsConnection::associateDigisToMcParticles()
432{
433 // --- get event T0 ---
434 double T0 = getEventStartTime();// in ns
435
436 // --- get CGEM MC hits
437 SmartDataPtr<Event::CgemMcHitCol> cgemMcHitCol(eventSvc(),"/Event/MC/CgemMcHitCol");
438 if(!cgemMcHitCol) cout<<"DotsConnection::associateDigisToMcParticles() does not find CgemMcHitCol!"<<endl;
439 //cout<<"CGEM MC hits obtained"<<endl;
440
441 // --- get CGEM cluster
442 SmartDataPtr<RecCgemClusterCol> aCgemClusterCol(eventSvc(),"/Event/Recon/RecCgemClusterCol");
443 if(!aCgemClusterCol) cout<<"DotsConnection::associateDigisToMcParticles() does not find RecCgemClusterCol!"<<endl;
444 RecCgemClusterCol::iterator iter_cluster_begin=aCgemClusterCol->begin();
445 RecCgemClusterCol::iterator iter_cluster;
446 //cout<<"CGEM clusters obtained"<<endl;
447
448 // --- get MDC MC hits
449 SmartDataPtr<Event::MdcMcHitCol> mdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
450 if(!mdcMcHitCol) cout<<"DotsConnection::associateDigisToMcParticles() does not find MdcMcHitCol!"<<endl;
451 if(myDebug)
452 cout<<"MDC MC hits obtained, n="<<mdcMcHitCol->size()<<endl;
453 Event::MdcMcHitCol::iterator iter_mdcMcHit = mdcMcHitCol->begin();
454 for(; iter_mdcMcHit!=mdcMcHitCol->end(); iter_mdcMcHit++ )
455 {
456 string creatorProcess = (*iter_mdcMcHit)->getCreatorProcess();
457 int trkId = (*iter_mdcMcHit)->getTrackIndex();
458 int isSec = (*iter_mdcMcHit)->getIsSecondary();
459 double trkLen = (*iter_mdcMcHit)->getFlightLength();
460 double px=(*iter_mdcMcHit)->getMomentumX();
461 double py=(*iter_mdcMcHit)->getMomentumY();
462 double pz=(*iter_mdcMcHit)->getMomentumZ();
463 double posx=(*iter_mdcMcHit)->getPositionX();
464 double posy=(*iter_mdcMcHit)->getPositionY();
465 double posz=(*iter_mdcMcHit)->getPositionZ();
466 int pdgCode = (*iter_mdcMcHit)->getCurrentTrackPID();
467 int digiIdx = (*iter_mdcMcHit)->getDigiIdx();
468 Identifier id = (*iter_mdcMcHit)->identify();
469 int layer = MdcID::layer(id);
470 int wire = MdcID::wire(id);
471 Hep3Vector pos(posx,posy,0);
472 Hep3Vector p3(px,py,0);
473 double dotProd = pos*p3;
474
475 if(myDebug)
476 cout<<" MDC MC hit from "<<creatorProcess
477 <<", trk "<<trkId
478 <<", isSec="<<isSec
479 <<", len="<<trkLen/10 // mm -> cm
480 <<", p3=("<<px<<", "<<py<<", "<<pz<<")"
481 <<", pos=("<<posx<<", "<<posy<<", "<<posz<<")"
482 <<", layer "<<layer<<" wire "<<wire
483 <<", dotProd="<<dotProd
484 <<", pdg code = "<<pdgCode
485 <<", digi index = "<<digiIdx
486 <<endl;
487 }
488
489 // --- get MDC digi vector
490 uint32_t getDigiFlag = 0;
491 MdcDigiVec mdcDigiVec = myRawDataProviderSvc->getMdcDigiVec(getDigiFlag);
492 if(myDebug)
493 cout<<"DotsConnection::associateDigisToMcParticles() get "<<mdcDigiVec.size()<<" MDC digis from RawDataProviderSvc"<<endl;
494 //vector<const MdcDigi*> aMdcDigiVec;
495 clearMdcDigiPointer();
496 vector<MdcDigi*>::iterator iter_mdcDigi = mdcDigiVec.begin();
497 for(; iter_mdcDigi!=mdcDigiVec.end(); iter_mdcDigi++)
498 {
499 // --- get id, layer, wire ---
500 Identifier id = (*iter_mdcDigi)->identify();
501 int layer = MdcID::layer(id);
502 //if( (layer>=8&&layer<=19) || (layer>=36&&layer<=42) ) // --- only axial hits kept
503 int wire = MdcID::wire(id);
504
505 double rawTime = RawDataUtil::MdcTime((*iter_mdcDigi)->getTimeChannel());
506 double tprop = myMdcCalibFunSvc->getTprop(layer, 0);
507 double charge = (*iter_mdcDigi)->getChargeChannel();
508 double T0Walk = myMdcCalibFunSvc->getT0(layer,wire) + myMdcCalibFunSvc->getTimeWalk(layer, charge);
509 double TOF = 0.;
510 double driftT = rawTime - T0 - TOF - T0Walk - tprop;
511 if(driftT>myDriftTimeUpLimit) continue;
512
513 myMdcDigiPointer[layer][wire]=*iter_mdcDigi;
514 }
515
516 // --- loop MC tracks to get the proper digits (CGEM clusters and MDC hits)
517 HepPoint3D origin(0, 0, 0);
518 int nMcTrk = myVecMCTrkId.size();
519 for(int i=0; i<nMcTrk; i++)
520 {
521 // --- get track index in mc particle collection
522 int trkIdx = myVecMCTrkId[i];
523 double trkLenInMdc = myVecTrkLenFirstHalf[i];
524 if(myDebug)
525 cout<<"MC particle "<<trkIdx<<", pdg code "<<myVecPDG[i]<<endl;
526
527 // --- associate CGEM clusters through CGEM MC hits
528 int CgemCluster[3][2]={0,0,0,0,0,0};// [layer][x/v]
529 myVecCgemXcluster.clear();
530 myVecCgemVcluster.clear();
531 myVecCgem1DCluster.clear();
532 for(int l=0; l<3; l++)
533 {
534 for(int m=0; m<2; m++)
535 {
536 myVecCgemXCluIdx[l][m].clear();
537 myVecCgemVCluIdx[l][m].clear();
538 }
539 }
540 Event::CgemMcHitCol::iterator iter_CgemMcHit = cgemMcHitCol->begin();
541 for(; iter_CgemMcHit!=cgemMcHitCol->end(); iter_CgemMcHit++ )
542 {
543 string creatorProcess = (*iter_CgemMcHit)->GetCreatorProcess();
544 if(myDebug)
545 cout<<" CGEM MC hit from process "<<creatorProcess;
546 if(creatorProcess=="Generator"||creatorProcess=="Decay")// only from generator or decay
547 {
548 //cout<<", trk id "<<(*iter_CgemMcHit)->GetTrackID()<<endl;
549 if((*iter_CgemMcHit)->GetTrackID()!=trkIdx) continue; // keep only matched MC hits
550 //cout<<" is secondary "<<(*iter_CgemMcHit)->GetIsSecondary()<<endl;
551 if((*iter_CgemMcHit)->GetIsSecondary()) continue;
552 const vector<int> & vec_xCluster = (*iter_CgemMcHit)->GetXclusterIdxVec();// find the X-clusters
553 int nXCluster=vec_xCluster.size();
554 for(int j=0; j<nXCluster; j++)
555 {
556 iter_cluster=iter_cluster_begin+vec_xCluster[j];
557 int clusterId = (*iter_cluster)->getclusterid();
558 int layer=(*iter_cluster)->getlayerid();
559 int sheet=(*iter_cluster)->getsheetid();
560 if(myDebug)
561 cout<<" find one X-cluster on layer "<<layer;
562 if(CgemCluster[layer][0]==0)
563 {
564 myVecCgemXcluster.push_back(*iter_cluster);
565 myVecCgem1DCluster.push_back(*iter_cluster);
566 myVecCgemXCluIdx[layer][sheet].push_back(clusterId);
567 CgemCluster[layer][0]++;
568 if(myDebug)
569 cout<<", associated. "<<endl;
570 }
571 }
572 //cout<<" X-clusters done "<<endl;
573 const vector<int> & vec_vCluster = (*iter_CgemMcHit)->GetVclusterIdxVec();// find the V-clusters
574 int nVCluster=vec_vCluster.size();
575 for(int j=0; j<nVCluster; j++)
576 {
577 iter_cluster=iter_cluster_begin+vec_vCluster[j];
578 int clusterId = (*iter_cluster)->getclusterid();
579 int layer=(*iter_cluster)->getlayerid();
580 int sheet=(*iter_cluster)->getsheetid();
581 if(myDebug)
582 cout<<" find one V-cluster on layer "<<layer;
583 if(CgemCluster[layer][1]==0)
584 {
585 myVecCgemVcluster.push_back(*iter_cluster);
586 myVecCgem1DCluster.push_back(*iter_cluster);
587 myVecCgemVCluIdx[layer][sheet].push_back(clusterId);
588 CgemCluster[layer][1]++;
589 if(myDebug)
590 cout<<", associated. "<<endl;
591 }
592 }
593 }
594 if(myDebug) cout<<endl;
595 }// end of looping CGEM MC hits
596
597 myVecMdcDigi.clear();
598 int nMdcXHits(0), nMdcVHits(0);
599
600 // --- associate MDC hits through MC track index
601 /*
602 vector<MdcDigi*>::iterator iter_mdcDigi = mdcDigiVec.begin();
603 for(; iter_mdcDigi!=mdcDigiVec.end(); iter_mdcDigi++)
604 {
605 int mcTrkIdx = (*iter_mdcDigi)->getTrackIndex();
606 if(mcTrkIdx==trkIdx)
607 {
608 //myVecMdcDigi.push_back((*iter_mdcDigi));
609 Identifier id = (*iter_mdcDigi)->identify();
610 int layer = MdcID::layer(id);
611 int wire = MdcID::wire(id);
612 //cout<<" MDC digi on layer "<<layer<<" wire "<<wire<<" associated"<<endl;
613 //cout<<" MDC digi on layer "<<layer<<" wire "<<wire<<" matched"<<endl;
614 }
615 }*/
616
617 // --- associate MDC hits through MdcMcHits
618 Event::MdcMcHitCol::iterator
619 iter_mdcMcHit = mdcMcHitCol->begin();
620 //int nXhits(0), nVhits(0);
621 for(; iter_mdcMcHit!=mdcMcHitCol->end(); iter_mdcMcHit++ )
622 {
623 int trkId = (*iter_mdcMcHit)->getTrackIndex();
624 if(trkId!=trkIdx) continue;
625 if((*iter_mdcMcHit)->getDigiIdx()==-9999) continue;
626 Identifier id = (*iter_mdcMcHit)->identify();
627 int layer = MdcID::layer(id);
628 int wire = MdcID::wire(id);
629 int isSec = (*iter_mdcMcHit)->getIsSecondary();
630 if(isSec!=0) {
631 /*cout<<" MDC digi on layer "<<layer<<" wire "<<wire<<" isSec"
632 <<", pid = "<<(*iter_mdcMcHit)->getCurrentTrackPID()
633 <<endl;*/
634 continue;
635 }
636 double trkLenMcHit = ((*iter_mdcMcHit)->getFlightLength())/10.;
637 if(trkLenMcHit>1.5*trkLenInMdc) {
638 //cout<<" MDC digi on layer "<<layer<<" wire "<<wire<<" trkLenMcHit>1.5*trkLenInMdc"<<endl;
639 continue;
640 }
641 double px=(*iter_mdcMcHit)->getMomentumX();
642 double py=(*iter_mdcMcHit)->getMomentumY();
643 double pz=(*iter_mdcMcHit)->getMomentumZ();
644 double posx=(*iter_mdcMcHit)->getPositionX();
645 double posy=(*iter_mdcMcHit)->getPositionY();
646 double posz=(*iter_mdcMcHit)->getPositionZ();
647 Hep3Vector pos(posx,posy,0);
648 Hep3Vector p3(px,py,0);
649 double dotProd = pos*p3;
650 if(dotProd<0) {
651 //cout<<" pos="<<pos<<", p3="<<p3<<endl;
652 //cout<<" MDC digi on layer "<<layer<<" wire "<<wire<<" coming back"<<endl;
653 continue;// reject hits from track coming back
654 }
655 //Identifier id = (*iter_mdcMcHit)->identify();
656 //int layer = MdcID::layer(id);
657 //int wire = MdcID::wire(id);
658 const MdcDigi* aMdcDigiPt = myMdcDigiPointer[layer][wire];
659 if(aMdcDigiPt!=NULL)
660 {
661 myVecMdcDigi.push_back(aMdcDigiPt);
662 myMdcDigiPointer[layer][wire]=NULL;
663 if(myDebug)
664 cout<<" MDC digi on layer "<<layer<<" wire "<<wire<<" associated"<<endl;
665 if(layer<8||(layer>19&&layer<36)) nMdcVHits++;
666 else nMdcXHits++;
667 }
668 }
669
670 // --- fit to Helix
671 //cout<<" start to fit "<<endl;
672 int fitFlag=99;
673 if((myVecCgemVcluster.size()+nMdcVHits)>=2&&(myVecCgemXcluster.size()+nMdcXHits+myVecCgemVcluster.size()+nMdcVHits)>5)
674 {
675 HepVector a_helix(5,0);
676 a_helix(1)=myVecHelix[i][0];
677 a_helix(2)=myVecHelix[i][1];
678 a_helix(3)=myVecHelix[i][2];
679 a_helix(4)=myVecHelix[i][3];
680 a_helix(5)=myVecHelix[i][4];
681 KalmanFit::Helix ini_helix(origin,a_helix);
682 myDotsHelixFitter.setInitialHelix(ini_helix);
683 myDotsHelixFitter.setCgemClusters(myVecCgem1DCluster);
684 myDotsHelixFitter.setDChits(myVecMdcDigi,T0);
685
686 myDotsHelixFitter.fitCircleOnly();
687 myDotsHelixFitter.useAxialHitsOnly();
688 int nIter=0;
689 while(1)
690 {
691 fitFlag = myDotsHelixFitter.calculateNewHelix();
692 nIter++;
693 //if(fitFlag!=0) break;
694 if(fitFlag==0)
695 {
696 if(myNtProd&1 && nIter==1) {
697 myNHitsCircle=0;
698 // --- fill CGEM clusters
699 vector<double> vecChiCgem = myDotsHelixFitter.getVecChiCgemCluster();
700 vector<const RecCgemCluster*>::iterator it_cluster=myVecCgem1DCluster.begin();
701 int i_cluster=0;
702 for(; it_cluster!=myVecCgem1DCluster.end(); it_cluster++, i_cluster++)
703 {
704 int flag=(*it_cluster)->getflag();
705 if(flag!=0) continue;
706 int layer=(*it_cluster)->getlayerid();
707 //if(flag==1) layer=-1*(layer+1);// for V-cluster
708 if(myNHitsCircle<100)
709 {
710 myLayerHitsCircle[myNHitsCircle]=layer;
711 myChiHitsCircle[myNHitsCircle] =vecChiCgem[i_cluster];
712 myNHitsCircle++;
713 }
714 }
715
716 // --- fill MDC hits
717 vector<const MdcDigi*>::iterator it_digi=myVecMdcDigi.begin();
718 vector<double> vecChiMdc = myDotsHelixFitter.getVecChiMdcDigi();
719 vector<int> vecIsActMdc = myDotsHelixFitter.getVecMdcDigiIsAct();
720 int i_digi=0;
721 for(; it_digi!=myVecMdcDigi.end(); it_digi++, i_digi++)
722 {
723 //if(vecIsActMdc[i_digi]==0) continue;
724 Identifier id = (*it_digi)->identify();
725 int layer = MdcID::layer(id);
726 if(myWireFlag[layer]!=0) continue;
727 if(myNHitsCircle<100)
728 {
729 myLayerHitsCircle[myNHitsCircle]=layer;
730 myChiHitsCircle[myNHitsCircle] =vecChiMdc[i_digi];
731 myNHitsCircle++;
732 }
733 }
734 }
735 }
736 else break;
737 if(myDotsHelixFitter.deactiveHits(myChiCut_circle,myNmaxDeact_circle)==0) break;
738 if(nIter>100) break;
739 }
740 if(myDebug) {
741 cout<<"initial helix "<<ini_helix.a()<<endl;
742 cout<<"fitFlag="<<fitFlag<<endl;
743 cout<<"nIter="<<nIter<<endl;
744 cout<<"fitted circle "<<myDotsHelixFitter.getHelix()<<endl;
745 cout<<"circle chi2="<<myDotsHelixFitter.getChi2()<<endl<<endl;
746 }
747 if(fitFlag==0) {
748 myDotsHelixFitter.fitModeHelix();
749 nIter=0;
750 while(1)
751 {
752 fitFlag = myDotsHelixFitter.calculateNewHelix();
753 if(fitFlag!=0) break;
754 if(myDotsHelixFitter.deactiveHits(myChiCut_helix,myNmaxDeact_helix)==0) break;
755 }
756 while(1)
757 {
758 if(fitFlag!=0) break;
759 if(myDotsHelixFitter.activeHits(myChiCut_helix)==0) break;
760 fitFlag = myDotsHelixFitter.calculateNewHelix();
761 }
762 if(myDebug) {
763 cout<<"fitFlag="<<fitFlag<<endl;
764 cout<<"fitted helix "<<myDotsHelixFitter.getHelix()<<endl;
765 cout<<"helix chi2="<<myDotsHelixFitter.getChi2()<<endl<<endl;
766 }
767 }
768 }
769
770 // --- save good track to RecMdcTrackCol and fill ntuple
771 //cout<<"before fill myArrayHelixMC"<<endl;
772 if(myNtProd&1) {
773 myNHits=0;
774 myRUN=myRun;
775 myEVT=myEvt;
776 //cout<<"fill 1"<<endl;
777 myPID=myVecPDG[i];
778 //cout<<"fill 2"<<endl;
779 myNPar=5;
780 myNXHits=nMdcXHits;
781 myNVHits=nMdcVHits;
782 myNXCluster=myVecCgemXcluster.size();
783 myNVCluster=myVecCgemVcluster.size();
784 //cout<<"fill 3"<<endl;
785 myArrayHelixMC[0]=myVecHelix[i][0];
786 //cout<<"fill 4"<<endl;
787 myArrayHelixMC[1]=myVecHelix[i][1];
788 myArrayHelixMC[2]=myVecHelix[i][2];
789 myArrayHelixMC[3]=myVecHelix[i][3];
790 myArrayHelixMC[4]=myVecHelix[i][4];
791 myTrkIdBest=-1;
792 myNHitsBestTrk=0;
793 myNSameHitsBestTrk=0;
794 }
795 //cout<<"after fill myArrayHelixMC"<<endl;
796 if(fitFlag==0)
797 {
798 saveARecMdcTrack(i);
799
800 if(myNtProd&1) {
801 // --- fill CGEM clusters
802 vector<double> vecChiCgem = myDotsHelixFitter.getVecChiCgemCluster();
803 vector<const RecCgemCluster*>::iterator it_cluster=myVecCgem1DCluster.begin();
804 int i_cluster=0;
805 for(; it_cluster!=myVecCgem1DCluster.end(); it_cluster++, i_cluster++)
806 {
807 int flag=(*it_cluster)->getflag();
808 int layer=(*it_cluster)->getlayerid();
809 if(flag==1) layer=-1*(layer+1);// for V-cluster
810 if(myNHits<100)
811 {
812 myLayerHits[myNHits]=layer;
813 myChiHits[myNHits] =vecChiCgem[i_cluster];
814 myNHits++;
815 }
816 }
817
818 // --- fill MDC hits
819 vector<const MdcDigi*>::iterator it_digi=myVecMdcDigi.begin();
820 vector<double> vecChiMdc = myDotsHelixFitter.getVecChiMdcDigi();
821 vector<int> vecIsActMdc = myDotsHelixFitter.getVecMdcDigiIsAct();
822 int i_digi=0;
823 for(; it_digi!=myVecMdcDigi.end(); it_digi++, i_digi++)
824 {
825 if(vecIsActMdc[i_digi]==0) continue;
826 Identifier id = (*it_digi)->identify();
827 int layer = MdcID::layer(id);
828 if(myNHits<100)
829 {
830 myLayerHits[myNHits]=layer;
831 myChiHits[myNHits] =vecChiMdc[i_digi];
832 myNHits++;
833 }
834
835 }
836 }// if myNtProd&1
837 }
838
839 // --- write out ntuple
840 if(myNtProd&1) myNtHelixFitter->write();
841
842 }// loop MC tracks
843
844 return;
845}
846
847vector<const MdcDigi*> DotsConnection::getMdcDigiVec()
848{
849 uint32_t getDigiFlag = 0;
850 /*
851 getDigiFlag += m_maxMdcDigi;
852 if(m_dropHot) getDigiFlag |= MdcRawDataProvider::b_dropHot;
853 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
854 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
855 */
856 MdcDigiVec mdcDigiVec = myRawDataProviderSvc->getMdcDigiVec(getDigiFlag);
857 if(myDebug)
858 cout<<"DotsConnection::getMdcDigiVec() get "<<mdcDigiVec.size()<<" MDC digis from RawDataProviderSvc"<<endl;
859 ///*
860 // MdcDigiVec aMdcDigiVec;
861 vector<const MdcDigi*> aMdcDigiVec;
862 vector<MdcDigi*>::iterator iter_mdcDigi = mdcDigiVec.begin();
863 for(; iter_mdcDigi!=mdcDigiVec.end(); iter_mdcDigi++) {
864 // --- get id, layer, wire ---
865 Identifier id = (*iter_mdcDigi)->identify();
866 int layer = MdcID::layer(id);
867
868 // --- only axial hits kept
869 //if( (layer>=8&&layer<=19) || (layer>=36&&layer<=42) )
870 aMdcDigiVec.push_back(*iter_mdcDigi);
871
872 }
873 return aMdcDigiVec;
874 //*/
875 //return mdcDigiVec;
876}
877
878double DotsConnection::getEventStartTime()
879{
880 double T0=0;
881
882 // --- get event start time ---
883 SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
884 if(evTimeCol){
885 RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
886 if(iter_evt != evTimeCol->end()){
887 //T0 = (*iter_evt)->getTest()*1.e-9;//m_bunchT0-s, getTest-ns
888 T0 = (*iter_evt)->getTest();// getTest-ns
889 }
890 }
891 else cout<<"error: DotsConnection::getEventStartTime() failed to access event start time"<<endl;
892
893 return T0;
894}
895
896void DotsConnection::testDotsHelixFitterAllHits()
897{
898 // --- set MC helix
899 KalmanFit::Helix ini_helix = getMCHelix();
900 myDotsHelixFitter.setInitialHelix(ini_helix);
901
902 // --- get all MDC digi
903 vector<const MdcDigi*> vecDigi = getMdcDigiVec();
904
905 // --- get event start time
906 double T0 = getEventStartTime();// in ns
907 //cout<<"DotsConnection::testDotsHelixFitterAllHits() T0 = "<<T0<<endl;
908 //myDotsHelixFitter.setT0(T0);
909
910 // --- set DC digi to DotsHelixFitter
911 myDotsHelixFitter.setDChits(vecDigi,T0);
912
913 // --- get DC digi ordered with the fligth length
914 vector<const MdcDigi*>::iterator iter_mdcDigi = vecDigi.begin();
915 map<double, const MdcDigi*> aMapDcDigi;
916 for(; iter_mdcDigi!=vecDigi.end(); iter_mdcDigi++)
917 {
918 aMapDcDigi[myDotsHelixFitter.getFlightLength(*iter_mdcDigi)] = *iter_mdcDigi;
919 }
920
921 vector<const RecCgemCluster*> vecCgemCluster = getCgemClusterVec();
922 myDotsHelixFitter.setCgemClusters(vecCgemCluster);
923
924 myDotsHelixFitter.fitCircleOnly();
925 myDotsHelixFitter.calculateNewHelix();
926
927 // --- fill a NTuple ---
928 myNPar=5;
929 HepVector a = ini_helix.a();
930 HepVector a_new = myDotsHelixFitter.getHelix();
931 for(int i=0; i<5; i++)
932 {
933 myArrayHelixMC[i]=a[i];
934 myArrayHelixFitted[i]=a_new[i];
935 }
936 myNtHelixFitter->write();
937}
938
939vector<const RecCgemCluster*> DotsConnection::getCgemClusterVec(int view)
940{
941 vector<const RecCgemCluster*> aVecCgemCluster;
942 SmartDataPtr<RecCgemClusterCol> aCgemClusterCol(eventSvc(),"/Event/Recon/RecCgemClusterCol");
943 if(aCgemClusterCol)
944 {
945 RecCgemClusterCol::iterator iter_cluster=aCgemClusterCol->begin();
946 int nCluster = aCgemClusterCol->size();
947 //cout<<"DotsConnection::getCgemClusterVec() finds a RecCgemClusterCol with "<<nCluster<<" clusters!"<<endl;
948 // cout<<"~~~~~~~~~~~~~~~~~~~~~~~~ check RecCgemClusterCol:"<<endl;
949 // cout <<setw(10)<<"idx"
950 // <<setw(10)<<"layer"
951 // <<setw(10)<<"sheet"
952 // <<setw(10)<<"XVFlag"
953 // <<setw(10)<<"id1 ~"
954 // <<setw(10)<<"id2"
955 // <<setw(15)<<"phi"
956 // <<setw(15)<<"V"
957 // <<setw(15)<<"z"
958 // <<endl;
959 for(; iter_cluster!=aCgemClusterCol->end(); iter_cluster++)
960 {
961 // cout <<setw(10)<<(*iter_cluster)->getclusterid()
962 // <<setw(10)<<(*iter_cluster)->getlayerid()
963 // <<setw(10)<<(*iter_cluster)->getsheetid()
964 // <<setw(10)<<(*iter_cluster)->getflag()
965 // <<setw(10)<<(*iter_cluster)->getclusterflagb()
966 // <<setw(10)<<(*iter_cluster)->getclusterflage()
967 // <<setw(15)<<setprecision(10)<<(*iter_cluster)->getrecphi()
968 // <<setw(15)<<setprecision(10)<<(*iter_cluster)->getrecv()
969 // <<setw(15)<<setprecision(10)<<(*iter_cluster)->getRecZ()
970 // <<endl;
971 int flag = (*iter_cluster)->getflag();
972 if(view==0||view==1)
973 {
974 if(flag==view) aVecCgemCluster.push_back(*iter_cluster);
975 }else if(view=2)
976 {
977 if(flag==0||flag==1) aVecCgemCluster.push_back(*iter_cluster);
978 }
979 }
980 // cout<<"~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
981 }
982 else cout<<"DotsConnection::getCgemClusterVec() does not find RecCgemClusterCol!"<<endl;
983 return aVecCgemCluster;
984}
985
986void DotsConnection::testDotsHelixFitterPartHits()
987{
988 // --- set MC helix
989 KalmanFit::Helix ini_helix = getMCHelix();
990 myDotsHelixFitter.setInitialHelix(ini_helix);
991
992 // --- get all MDC digi
993 vector<const MdcDigi*> vecDigi = getMdcDigiVec();
994
995 // --- get event start time
996 double T0 = getEventStartTime();// in ns
997 //cout<<"DotsConnection::testDotsHelixFitterPartHits() T0 = "<<T0<<endl;
998 myDotsHelixFitter.setT0(T0);
999
1000 // --- set DC digi to DotsHelixFitter
1001 //myDotsHelixFitter.setDChits(vecDigi,T0);
1002
1003 // --- get DC digi ordered with the fligth length
1004 vector<const MdcDigi*>::iterator iter_mdcDigi = vecDigi.begin();
1005 map<double, const MdcDigi*> aMapDcDigi;
1006 for(; iter_mdcDigi!=vecDigi.end(); iter_mdcDigi++)
1007 {
1008 aMapDcDigi[myDotsHelixFitter.getFlightLength(*iter_mdcDigi)] = *iter_mdcDigi;
1009 }
1010
1011 int layerUdStudy=8;
1012 int layerUsed=6;
1013 vector<const MdcDigi*> smallVecDigi;
1014 map<double, const MdcDigi*>::iterator iter_digi = aMapDcDigi.begin();
1015 int nLayerCount=-1;
1016 const MdcDigi* digiUdStudy=NULL;
1017 for(; iter_digi!=aMapDcDigi.end(); iter_digi++)
1018 {
1019 const MdcDigi* aDigi = iter_digi->second;
1020 int layer = myDotsHelixFitter.getLayer(aDigi);
1021 if(layer==layerUdStudy)
1022 {
1023 digiUdStudy=aDigi;
1024 nLayerCount=0;
1025 }
1026 if(layer>layerUdStudy&&nLayerCount>=0&&nLayerCount<layerUsed)
1027 {
1028 smallVecDigi.push_back(aDigi);
1029 nLayerCount++;
1030 }
1031 if(nLayerCount==layerUsed) break;
1032 }
1033 if( digiUdStudy!=NULL && nLayerCount==layerUsed )
1034 {
1035 myDotsHelixFitter.setDChits(smallVecDigi,T0);
1036 myDotsHelixFitter.fitCircleOnly();
1037 myDotsHelixFitter.calculateNewHelix();
1038 double doca = myDotsHelixFitter.getDocaFromTrk(digiUdStudy);
1039 //cout<<"layer "<<layerUdStudy<<" doca prediction from outer "<<layerUsed<<" digis : "<<doca<<endl;
1040 }
1041
1042
1043 // --- fill a NTuple ---
1044 /*
1045 myNPar=5;
1046 HepVector a = ini_helix.a();
1047 HepVector a_new = myDotsHelixFitter.getHelix();
1048 for(int i=0; i<5; i++)
1049 {
1050 myArrayHelixMC[i]=a[i];
1051 myArrayHelixFitted[i]=a_new[i];
1052 }
1053 myNtHelixFitter->write();
1054 */
1055}
1056
1057
1058void DotsConnection::clearMdcDigiPointer()
1059{
1060 for(int i=0; i<43; i++)
1061 for(int j=0; j<288; j++)
1062 myMdcDigiPointer[i][j]=NULL;
1063}
1064
1065bool DotsConnection::registerRecMdcTrack()
1066{
1067 MsgStream log(msgSvc(), name());
1068 StatusCode sc;
1069 IDataProviderSvc* eventSvc = NULL;
1070 service("EventDataSvc", eventSvc);
1071 if (eventSvc) {
1072 log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
1073 } else {
1074 log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
1075 return false;
1076 }
1077 IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*>(eventSvc);;
1078
1079 // --- register RecMdcTrackCol
1080 myRecMdcTrackCol = NULL; // new RecMdcTrackCol;
1081 DataObject *aRecMdcTrackCol;
1082 eventSvc->findObject("/Event/Recon/RecMdcTrackCol",aRecMdcTrackCol);
1083 if(aRecMdcTrackCol != NULL)
1084 {
1085 myRecMdcTrackCol = dynamic_cast<RecMdcTrackCol*> (aRecMdcTrackCol);
1086 //dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
1087 //eventSvc->unregisterObject("/Event/Recon/RecMdcTrackCol");
1088 }
1089 else {
1090 myRecMdcTrackCol = new RecMdcTrackCol;
1091 sc = eventSvc->registerObject("/Event/Recon/RecMdcTrackCol", myRecMdcTrackCol);
1092 if(sc.isFailure()) {
1093 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
1094 return false;
1095 }
1096 log << MSG::INFO << "RecMdcTrackCol registered successfully!" <<endreq;
1097 }
1098
1099 DataObject *aRecMdcHitCol;
1100 eventSvc->findObject("/Event/Recon/RecMdcHitCol",aRecMdcHitCol);
1101 if(aRecMdcHitCol != NULL) {
1102 //dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
1103 //eventSvc->unregisterObject("/Event/Recon/RecMdcHitCol");
1104 myRecMdcHitCol=dynamic_cast<RecMdcHitCol*> (aRecMdcHitCol);
1105 }
1106 else {
1107 myRecMdcHitCol= new RecMdcHitCol;
1108
1109 sc = eventSvc->registerObject("/Event/Recon/RecMdcHitCol", myRecMdcHitCol);
1110 if(sc.isFailure()) {
1111 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
1112 return false;
1113 }
1114 log << MSG::INFO << "RecMdcHitCol registered successfully!" <<endreq;
1115 }
1116
1117
1118 return true;
1119}
1120
1121bool sortCluster(const RecCgemCluster* clusterA , const RecCgemCluster* clusterB){
1122 return clusterA->getlayerid()<clusterB->getlayerid();
1123}
1124
1125bool DotsConnection::saveARecMdcTrack(int iTrk)
1126{
1127 int tkStat =4;
1128
1129 RecMdcTrack* recMdcTrack = new RecMdcTrack();
1130
1131 int trackId = myRecMdcTrackCol->size();
1132 recMdcTrack->setTrackId(trackId);
1133
1134 double helixPar[5];
1135 HepVector aHelixVec = myDotsHelixFitter.getHelix();
1136 helixPar[0]=aHelixVec[0];
1137 helixPar[1]=aHelixVec[1];
1138 helixPar[2]=aHelixVec[2];
1139 helixPar[3]=aHelixVec[3];
1140 helixPar[4]=aHelixVec[4];//helixPar[4]+=0.1;
1141 recMdcTrack->setHelix(helixPar);
1142
1143 // --- for ntuple
1144 if(myNtProd&1) {
1145 myArrayHelixFitted[0]=aHelixVec[0];
1146 myArrayHelixFitted[1]=aHelixVec[1];
1147 myArrayHelixFitted[2]=aHelixVec[2];
1148 myArrayHelixFitted[3]=aHelixVec[3];
1149 myArrayHelixFitted[4]=aHelixVec[4];
1150 }
1151
1152 int q = helixPar[2]>0? 1:-1;
1153 KalmanFit::Helix aHelix = myDotsHelixFitter.getClassHelix();
1154 double pxy = aHelix.pt();
1155 double px = aHelix.momentum(0).x();
1156 double py = aHelix.momentum(0).y();
1157 double pz = aHelix.momentum(0).z();
1158 double p = aHelix.momentum(0).mag();
1159 double theta = aHelix.direction(0).theta();
1160 double phi = aHelix.direction(0).phi();
1161 HepPoint3D poca = aHelix.x(0);
1162 HepPoint3D pivot = aHelix.pivot();
1163 double r = poca.perp();
1164 HepSymMatrix Ea = aHelix.Ea();
1165 //cout<<"Ea="<<Ea<<endl;
1166 double errorMat[15];
1167 int k = 0;
1168 for (int ie = 0 ; ie < 5 ; ie ++){
1169 for (int je = ie ; je < 5 ; je ++){
1170 errorMat[k] = Ea[ie][je];
1171 k++;
1172 }
1173 }
1174 double chisq = myDotsHelixFitter.getChi2();
1175 recMdcTrack->setCharge(q);
1176 recMdcTrack->setPxy(pxy);
1177 recMdcTrack->setPx(px);
1178 recMdcTrack->setPy(py);
1179 recMdcTrack->setPz(pz);
1180 recMdcTrack->setP(p);
1181 recMdcTrack->setTheta(theta);
1182 recMdcTrack->setPhi(phi);
1183 recMdcTrack->setPoca(poca);
1184 recMdcTrack->setX(poca.x());//poca
1185 recMdcTrack->setY(poca.y());
1186 recMdcTrack->setZ(poca.z());
1187 recMdcTrack->setR(sqrt(poca.x()*poca.x() + poca.y()*poca.y()));
1188 HepPoint3D apivot(0.,0.,0.);
1189 recMdcTrack->setPivot(apivot);
1190 recMdcTrack->setVX0(0.);//pivot
1191 recMdcTrack->setVY0(0.);
1192 recMdcTrack->setVZ0(0.);
1193 recMdcTrack->setError(errorMat);
1194 recMdcTrack->setError(Ea);
1195 recMdcTrack->setChi2(chisq);
1196 //recMdcTrack->setStat(tkStat);
1197 recMdcTrack->setStat(myVecPDG[iTrk]);
1198
1199 int maxLayerId = -1;
1200 int minLayerId = 43;
1201 double fiTerm = 0.;
1202 double fltLen = -0.00001;
1203 int layerMaxFltLen=-1;
1204
1205 // --- CGEM clusters
1206 ClusterRefVec clusterRefVec;
1207 map<int,int> clusterFitStat;
1208 SmartDataPtr<RecCgemClusterCol> aCgemClusterCol(eventSvc(),"/Event/Recon/RecCgemClusterCol");
1209 if(!aCgemClusterCol) cout<<"DotsConnection::saveARecMdcTrack() does not find RecCgemClusterCol!"<<endl;
1210 RecCgemClusterCol::iterator iter_cluster_begin=aCgemClusterCol->begin();
1211 RecCgemClusterCol::iterator iter_cluster=iter_cluster_begin;
1212 for(; iter_cluster!=aCgemClusterCol->end(); iter_cluster++)
1213 {
1214 int flag = (*iter_cluster)->getflag();
1215 if(flag!=2) continue;// skip 1D clusters
1216 int layer=(*iter_cluster)->getlayerid();
1217 int sheet=(*iter_cluster)->getsheetid();
1218 int idXClu = (*iter_cluster)->getclusterflagb();
1219 bool matchX=false;
1220 vector<int>::iterator iter=myVecCgemXCluIdx[layer][sheet].begin();
1221 for(; iter!=myVecCgemXCluIdx[layer][sheet].end(); iter++)
1222 if((*iter)==idXClu) matchX=true;
1223 if(!matchX) continue;
1224 int idVClu = (*iter_cluster)->getclusterflage();
1225 bool matchV=false;
1226 iter=myVecCgemVCluIdx[layer][sheet].begin();
1227 for(; iter!=myVecCgemVCluIdx[layer][sheet].end(); iter++)
1228 if((*iter)==idVClu) matchV=true;
1229 if(matchV)
1230 {
1231 const RecCgemCluster* recCgemCluster = (*iter_cluster);
1232 int clusterid = recCgemCluster->getclusterid();
1233 clusterRefVec.push_back(recCgemCluster);
1234 clusterFitStat[clusterid] = 1;
1235 if(maxLayerId<layer)
1236 {
1237 maxLayerId=layer;
1238 }
1239 }
1240 }
1241
1242
1243 // --- MDC hits
1244 set<int> setMdcIdt;
1245 int hitId = 0;
1246 HitRefVec hitRefVec;
1247 vector<RecMdcHit> aRecMdcHitVec=myDotsHelixFitter.makeRecMdcHitVec(1);
1248 int nMdcHits=aRecMdcHitVec.size();
1249 int nMdcHitsKept=0;
1250 vector<RecMdcHit>::iterator iter_recMdcHit = aRecMdcHitVec.begin();
1251 for(; iter_recMdcHit!=aRecMdcHitVec.end(); iter_recMdcHit++)
1252 {
1253 if(iter_recMdcHit->getChisqAdd()>myMdcHitChi2Cut) // skip hit with chi2>myMdcHitChi2Cut
1254 continue;
1255
1256 RecMdcHit* recMdcHit = new RecMdcHit(*iter_recMdcHit);
1257 recMdcHit->setId(hitId);
1258 recMdcHit->setTrkId(trackId);
1259 recMdcHit->setStat(1);
1260 myRecMdcHitCol->push_back(recMdcHit);
1261 SmartRef<RecMdcHit> refHit(recMdcHit);
1262 hitRefVec.push_back(refHit);
1263 nMdcHitsKept++;
1264
1265 Identifier mdcid = recMdcHit->getMdcId();
1266 int layer = MdcID::layer(mdcid);
1267 int wire = MdcID::wire(mdcid);
1268 setMdcIdt.insert(layer*1000+wire);
1269 if(layer>maxLayerId)
1270 {
1271 maxLayerId = layer;
1272 }
1273 if(layer<minLayerId)
1274 {
1275 minLayerId = layer;
1276 }
1277 if(fltLen<recMdcHit->getFltLen()) {
1278 fltLen=recMdcHit->getFltLen();
1279 layerMaxFltLen=layer;
1280 }
1281 hitId++;
1282 }
1283 if(myDebug)
1284 cout<<"track "<<trackId<<", "<<nMdcHitsKept<<"/"<<nMdcHits<<" hits kept"<<endl;
1285
1286 // --- phi term (phi for the outmost hit/cluster)
1287 if(maxLayerId>=0&&maxLayerId<3) {
1288 double rmax=myDotsHelixFitter.getRmidGapCgem(maxLayerId);
1289 fltLen=aHelix.flightLength(rmax);
1290 }
1291 if(fltLen>0) fiTerm=-fltLen*sin(theta)/aHelix.radius();
1292 else cout<<"fltLen<0!"<<endl;
1293 recMdcTrack->setFiTerm(fiTerm);
1294
1295 // --- check fiTerm
1296 HepPoint3D posOut = aHelix.x(fiTerm);
1297 if(myDebug)
1298 cout<<"fiTerm, layer, fltLen, pos= "<<fiTerm<<", "<<layerMaxFltLen<<", "<<fltLen<<", "<<posOut<<endl;
1299
1300 // --- setN* functions called in setVec* functions
1301 recMdcTrack->setVecHits(hitRefVec);
1302 std::sort(clusterRefVec.begin(),clusterRefVec.end(),sortCluster);
1303 recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
1304
1305 // --- check if any track share the same hits
1306 int n_CgemCluster(0), n_MdcHits(0), n_hits_tot(0);
1307 int n_sameCgemCluster(0), n_sameMdcHits(0), n_sameHit_tot(0);
1308 int bestTrkId(-1);
1309 RecMdcTrackCol::iterator iter_bestTrk;
1310 RecMdcTrackCol::iterator iter_trk = myRecMdcTrackCol->begin();
1311 for(; iter_trk!=myRecMdcTrackCol->end(); iter_trk++)
1312 {
1313 int stat =(*iter_trk)->stat();
1314 if(stat!=4) break;
1315 n_sameCgemCluster=0;
1316 ClusterRefVec aClusterRefVec=(*iter_trk)->getVecClusters();
1317 n_CgemCluster=aClusterRefVec.size();
1318 ClusterRefVec::iterator itCluster=aClusterRefVec.begin();
1319 for(; itCluster!=aClusterRefVec.end(); itCluster++) {
1320 int clusterid = (*itCluster)->getclusterid();
1321 if(clusterFitStat.count(clusterid)) n_sameCgemCluster++;
1322 }
1323
1324 n_sameMdcHits=0;
1325 HitRefVec aHitRefVec = (*iter_trk)->getVecHits();
1326 n_MdcHits=aHitRefVec.size();
1327 HitRefVec::iterator it_hit = aHitRefVec.begin();
1328 for(; it_hit!=aHitRefVec.end(); it_hit++) {
1329 Identifier mdcid = (*it_hit)->getMdcId();
1330 int layid = MdcID::layer(mdcid);
1331 int localwid = MdcID::wire(mdcid);
1332 if(setMdcIdt.find(layid*1000+localwid)!=setMdcIdt.end()) n_sameMdcHits++;
1333 }
1334 if((n_sameMdcHits+2*n_sameCgemCluster)>n_sameHit_tot) {
1335 n_sameHit_tot=n_sameMdcHits+2*n_sameCgemCluster;
1336 n_hits_tot=n_MdcHits+2*n_CgemCluster;
1337 iter_bestTrk=iter_trk;
1338 bestTrkId=(*iter_trk)->trackId();
1339 }
1340 }
1341 myTrkIdBest=bestTrkId;
1342 myNHitsBestTrk=n_hits_tot;
1343 myNSameHitsBestTrk=n_sameHit_tot;
1344 if(bestTrkId!=-1) (*iter_bestTrk)->setStat(400000*myVecPDG[iTrk]/abs(myVecPDG[iTrk])+myVecPDG[iTrk]);
1345
1346 // --- put a new track into RecMdcTrackCol
1347 myRecMdcTrackCol->push_back(recMdcTrack);
1348
1349 return true;
1350}
1351
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
sortCluster
bool sortCluster(const RecCgemCluster *clusterA, const RecCgemCluster *clusterB)
Definition: CgemMdcFitAlg.cxx:768
phi2
Double_t phi2
Definition: DataBase/tau_mode.c:7
phi1
Double_t phi1
Definition: DataBase/tau_mode.c:7
origin
const DifPoint origin
sortCluster
bool sortCluster(const RecCgemCluster *clusterA, const RecCgemCluster *clusterB)
Definition: DotsConnection.cxx:1121
abs
double abs(const EvtComplex &c)
Definition: EvtComplex.hh:212
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition: EvtStreamInputIterator.hh:101
q
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition: KKsem.h:33
MdcDigiVec
std::vector< MdcDigi * > MdcDigiVec
Definition: MdcRawDataProvider.h:11
RecMdcHitCol
ObjectVector< RecMdcHit > RecMdcHitCol
Definition: RecMdcHit.h:99
ClusterRefVec
SmartRefVector< RecCgemCluster > ClusterRefVec
Definition: RecMdcTrack.h:28
RecMdcTrackCol
ObjectVector< RecMdcTrack > RecMdcTrackCol
Definition: RecMdcTrack.h:102
HitRefVec
SmartRefVector< RecMdcHit > HitRefVec
Definition: RecMdcTrack.h:26
ntupleSvc
INTupleSvc * ntupleSvc()
Definition: ServiceAccessor.h:55
msgSvc
IMessageSvc * msgSvc()
Definition: ServiceAccessor.h:13
M_PI
#define M_PI
Definition: TConstant.h:4
DotsConnection::DotsConnection
DotsConnection(const std::string &name, ISvcLocator *pSvcLocator)
Definition: DotsConnection.cxx:23
DotsConnection::execute
StatusCode execute()
Definition: DotsConnection.cxx:204
DotsConnection::finalize
StatusCode finalize()
Definition: DotsConnection.cxx:254
DotsConnection::initialize
StatusCode initialize()
Definition: DotsConnection.cxx:44
DotsHelixFitter::setInitialHelix
void setInitialHelix(KalmanFit::Helix aHelix)
Definition: DotsHelixFitter.cxx:702
DotsHelixFitter::setChi2Diverge
void setChi2Diverge(double cut=1000000)
Definition: DotsHelixFitter.h:30
DotsHelixFitter::getFlightLength
double getFlightLength()
Definition: DotsHelixFitter.h:61
DotsHelixFitter::deactiveHits
int deactiveHits(double chi_cut=10, int nMax=1)
Definition: DotsHelixFitter.cxx:600
DotsHelixFitter::getLayer
int getLayer(const MdcDigi *aDcDigi)
Definition: DotsHelixFitter.h:63
DotsHelixFitter::setDChits
void setDChits(vector< const MdcDigi * > aVecMdcDigi, double T0)
Definition: DotsHelixFitter.cxx:139
DotsHelixFitter::setCgemClusters
void setCgemClusters(vector< const RecCgemCluster * > aVecCgemCluster)
Definition: DotsHelixFitter.cxx:159
DotsHelixFitter::useAxialHitsOnly
void useAxialHitsOnly(bool x=true)
Definition: DotsHelixFitter.h:40
DotsHelixFitter::getRmidGapCgem
double getRmidGapCgem(int i)
Definition: DotsHelixFitter.h:82
DotsHelixFitter::getVecChiCgemCluster
vector< double > getVecChiCgemCluster()
Definition: DotsHelixFitter.h:77
DotsHelixFitter::getClassHelix
KalmanFit::Helix getClassHelix()
Definition: DotsHelixFitter.h:42
DotsHelixFitter::getHelix
HepVector getHelix()
Definition: DotsHelixFitter.h:43
DotsHelixFitter::fitCircleOnly
void fitCircleOnly(bool x=true)
Definition: DotsHelixFitter.h:37
DotsHelixFitter::makeRecMdcHitVec
vector< RecMdcHit > makeRecMdcHitVec(int sel=1)
Definition: DotsHelixFitter.cxx:885
DotsHelixFitter::getVecChiMdcDigi
vector< double > getVecChiMdcDigi()
Definition: DotsHelixFitter.h:72
DotsHelixFitter::setT0
void setT0(double T0)
Definition: DotsHelixFitter.h:33
DotsHelixFitter::getVecMdcDigiIsAct
vector< int > getVecMdcDigiIsAct()
Definition: DotsHelixFitter.h:73
DotsHelixFitter::getDocaFromTrk
double getDocaFromTrk()
Definition: DotsHelixFitter.h:64
DotsHelixFitter::activeHits
int activeHits(double chi_cut=10)
Definition: DotsHelixFitter.cxx:650
DstMdcTrack::setPxy
void setPxy(const double pxy)
Definition: DstMdcTrack.h:86
DstMdcTrack::setTrackId
void setTrackId(const int trackId)
Definition: DstMdcTrack.h:84
DstMdcTrack::setPy
void setPy(const double py)
Definition: DstMdcTrack.h:88
DstMdcTrack::setZ
void setZ(const double z)
Definition: DstMdcTrack.h:95
DstMdcTrack::setX
void setX(const double x)
Definition: DstMdcTrack.h:93
DstMdcTrack::setError
void setError(double err[15])
Definition: DstMdcTrack.cxx:118
DstMdcTrack::setTheta
void setTheta(const double theta)
Definition: DstMdcTrack.h:91
DstMdcTrack::setStat
void setStat(const int stat)
Definition: DstMdcTrack.h:97
DstMdcTrack::setP
void setP(const double p)
Definition: DstMdcTrack.h:90
DstMdcTrack::setHelix
void setHelix(double helix[5])
Definition: DstMdcTrack.cxx:104
DstMdcTrack::setPoca
void setPoca(double poca[3])
Definition: DstMdcTrack.cxx:111
DstMdcTrack::setR
void setR(const double r)
Definition: DstMdcTrack.h:96
DstMdcTrack::setCharge
void setCharge(const int charge)
Definition: DstMdcTrack.h:85
DstMdcTrack::setY
void setY(const double y)
Definition: DstMdcTrack.h:94
DstMdcTrack::setChi2
void setChi2(const double chi)
Definition: DstMdcTrack.h:98
DstMdcTrack::setPhi
void setPhi(const double phi)
Definition: DstMdcTrack.h:92
DstMdcTrack::setPz
void setPz(const double pz)
Definition: DstMdcTrack.h:89
DstMdcTrack::setPx
void setPx(const double px)
Definition: DstMdcTrack.h:87
KalmanFit::Helix::momentum
Hep3Vector momentum(double dPhi=0.) const
returns momentum vector after rotating angle dPhi in phi direction.
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-15-14/src/helix/Helix.cxx:206
KalmanFit::Helix::Ea
const HepSymMatrix & Ea(void) const
returns error matrix.
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-15-14/KalFitAlg/helix/Helix.h:271
KalmanFit::Helix::x
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-15-14/src/helix/Helix.cxx:151
KalmanFit::Helix::a
const HepVector & a(void) const
returns helix parameters.
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-15-14/KalFitAlg/helix/Helix.h:265
KalmanFit::Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-15-14/KalFitAlg/helix/Helix.h:211
KalmanFit::Helix::direction
Hep3Vector direction(double dPhi=0.) const
returns direction vector after rotating angle dPhi in phi direction.
Definition: Reconstruction/KalFitAlg/KalFitAlg-00-15-14/KalFitAlg/helix/Helix.h:223
MdcCalibFunSvc::getT0
double getT0(int layid, int cellid) const
Definition: MdcCalibFunSvc.cxx:507
MdcCalibFunSvc::getTprop
double getTprop(int lay, double z) const
Definition: MdcCalibFunSvc.cxx:143
MdcCalibFunSvc::getTimeWalk
double getTimeWalk(int layid, double Q) const
Definition: MdcCalibFunSvc.cxx:514
MdcGeoLayer::Radius
double Radius(void) const
Definition: MdcGeoLayer.h:160
MdcGeoLayer::nomShift
double nomShift(void) const
Definition: MdcGeoLayer.h:169
MdcGeoLayer::NCell
int NCell(void) const
Definition: MdcGeoLayer.h:165
MdcGeoWire::Forward
HepPoint3D Forward(void) const
Definition: MdcGeoWire.h:129
MdcGeomSvc::Wire
const MdcGeoWire *const Wire(unsigned id)
Definition: MdcGeomSvc.cxx:770
MdcGeomSvc::getSuperLayerSize
const int getSuperLayerSize()
Definition: MdcGeomSvc.cxx:681
MdcGeomSvc::Layer
const MdcGeoLayer *const Layer(unsigned id)
Definition: MdcGeomSvc.cxx:786
MdcGeomSvc::getSegmentNo
const int getSegmentNo()
Definition: MdcGeomSvc.cxx:691
MdcGeomSvc::getWireSize
const int getWireSize()
Definition: MdcGeomSvc.cxx:671
MdcGeomSvc::getLayerSize
const int getLayerSize()
Definition: MdcGeomSvc.cxx:676
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
RawDataProviderSvc::getMdcDigiVec
MdcDigiVec & getMdcDigiVec(uint32_t control=0)
Definition: RawDataProviderSvc.h:29
RawDataUtil::MdcTime
static double MdcTime(int timeChannel)
Definition: RawDataUtil.h:8
RecCgemCluster::getclusterid
int getclusterid(void) const
Definition: RecCgemCluster.h:53
RecCgemCluster::getlayerid
int getlayerid(void) const
Definition: RecCgemCluster.h:55
RecMdcHit::getMdcId
const Identifier getMdcId(void) const
Definition: RecMdcHit.h:49
RecMdcHit::setStat
void setStat(int stat)
Definition: RecMdcHit.h:66
RecMdcHit::getFltLen
const double getFltLen(void) const
Definition: RecMdcHit.h:56
RecMdcHit::setTrkId
void setTrkId(int trkid)
Definition: RecMdcHit.h:59
RecMdcHit::setId
void setId(int id)
Definition: RecMdcHit.h:58
RecMdcTrack::setVecClusters
void setVecClusters(ClusterRefVec vecclusters)
Definition: RecMdcTrack.cxx:86
RecMdcTrack::setPivot
void setPivot(const HepPoint3D &pivot)
Definition: RecMdcTrack.h:79
RecMdcTrack::setVZ0
void setVZ0(double z0)
Definition: RecMdcTrack.h:76
RecMdcTrack::setVY0
void setVY0(double y0)
Definition: RecMdcTrack.h:75
RecMdcTrack::setVecHits
void setVecHits(HitRefVec vechits)
Definition: RecMdcTrack.cxx:80
RecMdcTrack::setFiTerm
void setFiTerm(double fiterm)
Definition: RecMdcTrack.h:77
RecMdcTrack::setVX0
void setVX0(double x0)
Definition: RecMdcTrack.h:74
Event
Definition: Event.h:21