CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
MdcxCosmicSewer Class Reference

#include <MdcxCosmicSewer.h>

+ Inheritance diagram for MdcxCosmicSewer:

Public Member Functions

 MdcxCosmicSewer (const std::string &name, ISvcLocator *pSvcLocator)
 
virtual ~MdcxCosmicSewer ()
 
StatusCode initialize ()
 
StatusCode execute ()
 
StatusCode finalize ()
 
StatusCode beginRun ()
 
void MdcxHitsToHots (HepVector &helix, TrkRecoTrk *trk, HitRefVec &hits, HitRefVec &skiped)
 
void store (TrkRecoTrk *tk, HitRefVec &skip)
 
void clearRecMdcTrackHit ()
 
void dumpTdsTrack (RecMdcTrackCol *trackList)
 

Detailed Description

Definition at line 43 of file MdcxCosmicSewer.h.

Constructor & Destructor Documentation

◆ MdcxCosmicSewer()

MdcxCosmicSewer::MdcxCosmicSewer ( const std::string &  name,
ISvcLocator *  pSvcLocator 
)

Definition at line 58 of file MdcxCosmicSewer.cxx.

58 :
59 Algorithm(name, pSvcLocator)
60{
61
62 declareProperty("debug", m_debug = false);
63 declareProperty("hist", m_hist = false);
64
65 declareProperty("doSag", m_doSag = false);
66
67 declareProperty("cosmicSewPar", m_cosmicSewPar);
68 declareProperty("cosmicSewSkip", m_cosmicSewSkip=false);
69
70 declareProperty("countPropTime", m_countPropTime = true);
71 declareProperty("lineFit", m_lineFit = false);
72}

◆ ~MdcxCosmicSewer()

MdcxCosmicSewer::~MdcxCosmicSewer ( )
virtual

Definition at line 77 of file MdcxCosmicSewer.cxx.

77 {
78 delete m_bfield;
79 delete m_context;
80}

Member Function Documentation

◆ beginRun()

StatusCode MdcxCosmicSewer::beginRun ( )

Definition at line 82 of file MdcxCosmicSewer.cxx.

82 {
83 //Initailize MdcDetector
84 m_gm = MdcDetector::instance(m_doSag);
85 if(NULL == m_gm) return StatusCode::FAILURE;
86
87 return StatusCode::SUCCESS;
88}
static MdcDetector * instance()
Definition: MdcDetector.cxx:21

◆ clearRecMdcTrackHit()

void MdcxCosmicSewer::clearRecMdcTrackHit ( )

Definition at line 639 of file MdcxCosmicSewer.cxx.

639 {
640 MsgStream log(msgSvc(), name());
641 StatusCode sc;
642 //-------------------------------------------------------
643 //clear TDS
644 //-------------------------------------------------------
645 //IDataManagerSvc *dataManSvc;
646 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
647 DataObject *aTrackCol;
648 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
649 if(aTrackCol != NULL) {
650 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
651 eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
652 }
653 DataObject *aRecHitCol;
654 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
655 if(aRecHitCol != NULL) {
656 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
657 eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
658 }
659 SmartDataPtr<RecMdcTrackCol> trackListTemp(eventSvc(),EventModel::Recon::RecMdcTrackCol);
660 if (!trackListTemp) {
661 RecMdcTrackCol *trackListTemp= new RecMdcTrackCol;
662 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackListTemp);
663 if(!sc.isSuccess()) {
664 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
665 }
666 }
667 SmartDataPtr<RecMdcHitCol> hitListTemp(eventSvc(),EventModel::Recon::RecMdcHitCol);
668 if (!hitListTemp) {
669 RecMdcHitCol *hitListTemp= new RecMdcHitCol;
670 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitListTemp);
671 if(!sc.isSuccess()) {
672 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
673 }
674 }
675}
ObjectVector< RecMdcHit > RecMdcHitCol
Definition: RecMdcHit.h:99
ObjectVector< RecMdcTrack > RecMdcTrackCol
Definition: RecMdcTrack.h:102
IMessageSvc * msgSvc()
_EXTERN_ std::string RecMdcTrackCol
Definition: EventModel.h:92
_EXTERN_ std::string RecMdcHitCol
Definition: EventModel.h:91

Referenced by execute().

◆ dumpTdsTrack()

void MdcxCosmicSewer::dumpTdsTrack ( RecMdcTrackCol trackList)

◆ execute()

StatusCode MdcxCosmicSewer::execute ( )

Definition at line 170 of file MdcxCosmicSewer.cxx.

170 {
171 MsgStream log(msgSvc(), name());
172 log << MSG::INFO << "in execute()" << endreq;
173
174
175 setFilterPassed(false);
176 // Get event No.
177 SmartDataPtr<Event::EventHeader> evtHead(eventSvc(),"/Event/EventHeader");
178 if (!evtHead) {
179 log << MSG::FATAL<< "Could not retrieve event header" << endreq;
180 return StatusCode::FAILURE;
181 }
182 long t_runNo = evtHead->runNumber();
183 long t_evtNo = evtHead->eventNumber();
184 if (m_debug) std::cout << "sew "<<++i_evt<<" run "<<t_runNo<<" evt " << t_evtNo << std::endl;
185
186
187 // Get bunch time t0 (ns)
188 m_bunchT0 = -999.;
189 SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
190 if (!aevtimeCol || aevtimeCol->size()==0) {
191 log << MSG::WARNING<< "Could not find RecEsTimeCol"<< endreq;
192 return StatusCode::SUCCESS;
193 }
194 RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
195 for(; iter_evt!=aevtimeCol->end(); iter_evt++){
196 m_bunchT0 = (*iter_evt)->getTest();
197 }
198
199 //read event data
200 SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(),EventModel::Recon::RecMdcTrackCol);
201 SmartDataPtr<RecMdcHitCol> recMdcHitCol(eventSvc(),EventModel::Recon::RecMdcHitCol);
202 if (!recMdcTrackCol || ! recMdcHitCol) return StatusCode::SUCCESS;
203 if (2!=recMdcTrackCol->size()){
205 return StatusCode::SUCCESS;
206 }
207 if(m_debug)std::cout<<name()<<" nTk==2 begin sewing"<< std::endl;
208
209
210 //get best track
211 double dParam[5]={0.,0.,0.,0.,0.};
212 float bprob =0.;
213 float chisum =0.;
214 RecMdcTrackCol::iterator besthel;
215
216 int besthelId = -999;
217 RecMdcTrackCol::iterator it = recMdcTrackCol->begin();
218 for (int iTk=0; it != recMdcTrackCol->end(); it++,iTk++) {
219 RecMdcTrack* tk = *it;
220 double Bz = m_bfield->bFieldNominal();
221 double d0 = -tk->helix(0); //cm
222 double phi0 = tk->helix(1)+ CLHEP::halfpi;
223 double omega = Bz*tk->helix(2)/-333.567;
224 double z0 = tk->helix(3); //cm
225 double tanl = tk->helix(4); //cm
226
227 if(phi0>Constants::pi)phi0-=Constants::twoPi;
228 if(phi0<-Constants::pi)phi0+=Constants::twoPi;
229
230 if(m_debug)std::cout<<__FILE__<<" "<<__LINE__<< "tk"<<iTk
231 <<"("<<d0<<","<<phi0<<"," <<","<<omega<<","<<z0<<","<<tanl<<")"<<std::endl;
232
233 if(iTk==0){//store param1
234 dParam[0]=d0;
235 dParam[1]=phi0;
236 dParam[2]=omega;
237 dParam[3]=z0;
238 dParam[4]=tanl;
239 }else{//calc. diff between param1 and param2
240 dParam[0] += d0;
241 dParam[1] -= (phi0+Constants::pi);
242 dParam[2] -= omega;
243 dParam[3] -= z0;
244 dParam[4] += tanl;
245 }
246
247 if(phi0<0) {
248 besthelId = iTk;
249 besthel = it;
250 }
251 float bchisq = tk->chi2();
252 int bndof = tk->ndof();
253 bprob=Mdcxprobab(bndof,bchisq);
254 chisum += bchisq;
255 //if (bprob > bestprob){
256 //bestprob = bprob;
257 //besthel = it;
258 //}
259 }//end for recMdcTrackCol
260
261 if(besthelId == -999) return StatusCode::SUCCESS;
262
263
264 TrkHelixMaker helixfactory;
265 TrkLineMaker linefactory;
266
267 if(m_debug){
268 std::cout<<__FILE__<<" param diff: " << "\n D0 " << dParam[0]
269 << "\n Phi0 " << dParam[1] << "\n Omega " << dParam[2]
270 << "\n Z0 " << dParam[3] << "\n Tanl " << dParam[4] << std::endl;
271 }
272
273 if(m_hist){
274 m_csmcD0= dParam[0];
275 m_csmcPhi0=dParam[1];
276 m_csmcOmega=dParam[2];
277 m_csmcZ0=dParam[3];
278 m_csmcTanl=dParam[4];
279 m_xtupleCsmcSew->write();
280 }
281
282 //get track 5 parameters
283 if( (fabs(dParam[0])>m_cosmicSewPar[0]) ||
284 (fabs(dParam[1])>m_cosmicSewPar[1]) ||
285 (fabs(dParam[2])>m_cosmicSewPar[2]) ||
286 (fabs(dParam[3])>m_cosmicSewPar[3]) ||
287 (fabs(dParam[4])>m_cosmicSewPar[4]) ){
288 if(m_debug)std::cout<<__FILE__<<" 2 trk param not satisfied skip!"<<std::endl;
289 if(m_debug){
290 if(fabs(dParam[0])>m_cosmicSewPar[0]) std::cout<<" cut by d0 "<<std::endl;
291 if(fabs(dParam[1])>m_cosmicSewPar[1]) std::cout<<" cut by phi0 "<<std::endl;
292 if(fabs(dParam[2])>m_cosmicSewPar[2]) std::cout<<" cut by omega "<<std::endl;
293 if(fabs(dParam[3])>m_cosmicSewPar[3]) std::cout<<" cut by z0"<<std::endl;
294 if(fabs(dParam[4])>m_cosmicSewPar[4]) std::cout<<" cut by tanl"<<std::endl;
295 }
297 return StatusCode::SUCCESS;
298 }
299
300 //got 2 sew-able trks, get helix param
301 RecMdcTrack* tk = *besthel;
302 double Bz = m_bfield->bFieldNominal();
303 double d0 = -tk->helix(0); //cm
304 double phi0 = tk->helix(1)+ CLHEP::halfpi;
305 double omega = Bz*tk->helix(2)/-333.567;
306 double z0 = tk->helix(3); //cm
307 double tanl = tk->helix(4); //cm
308
309 TrkExchangePar tt(d0,phi0,omega,z0,tanl);
310 if(m_debug)std::cout<<__FILE__<<" best helix: No."<<tk->trackId()<<" Pat param=("<<d0<<" "<<phi0
311 <<" "<<omega<<" "<<z0<<" "<<tanl<<")"<< std::endl;
312
313
314 //new track
315 TrkRecoTrk* newTrack;
316 if(m_lineFit){
317 newTrack = linefactory.makeTrack(tt, chisum, *m_context, m_bunchT0*1.e-9);
318 linefactory.setFlipAndDrop(*newTrack, true, true);
319 }else{
320 newTrack = helixfactory.makeTrack(tt, chisum, *m_context, m_bunchT0*1.e-9);
321 helixfactory.setFlipAndDrop(*newTrack, true, true);
322 }
323
324
325 // combine hit list
326 HitRefVec skippedHits;
327 it = recMdcTrackCol->begin();
328 for (int iTk=0; it != recMdcTrackCol->end(); it++,iTk++) {
329 RecMdcTrack* tk = *it;
330 HitRefVec hits = tk->getVecHits();
331 HepVector helixPatPar = m_mdcUtilitySvc->besPar2PatPar(tk->helix());
332 MdcxHitsToHots(helixPatPar,newTrack,hits,skippedHits);
333 }
334
335
336 //------------------------------------
337 //do fit
338 //------------------------------------
339 TrkErrCode err=newTrack->hits()->fit();
340 const TrkFit* newFit = newTrack->fitResult();
341
342
343 if (!newFit) {
344 // Fit bad
345 if(m_debug)std::cout<<__FILE__<<" sew fit failed!!!"<< std::endl;
346 HitRefVec::iterator it= skippedHits.begin();
347 for (;it!=skippedHits.end();++it){ delete it->data();}
348 delete newTrack;
350 } else {
351 //------------------------------------
352 // Fit good
353 //------------------------------------
354 if(m_lineFit){
355 linefactory.setFlipAndDrop(*newTrack, false, false);
356 }else{
357 helixfactory.setFlipAndDrop(*newTrack, false, false);
358 }
359 if(m_debug){
360 std::cout<<__FILE__<<" sew cosmic fit good "<< std::endl;
361 newTrack->print(std::cout);
362 }
363
364 SmartDataPtr<MdcHitCol> m_hitCol(eventSvc(), "/Event/Hit/MdcHitCol");
365 if (!m_hitCol){
366 m_hitCol= new MdcHitCol;
367 StatusCode sc = eventSvc()->registerObject("/Event/Hit/MdcHitCol",m_hitCol);
368 if(!sc.isSuccess()) {
369 std::cout<< " Could not register MdcHitCol" <<endreq;
370 return StatusCode::FAILURE;
371 }
372 }
373 uint32_t getDigiFlag = 0;
374 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
375 const MdcDigi* m_digiMap[43][288];
376 MdcDigiVec::iterator iter = mdcDigiVec.begin();
377 for (; iter != mdcDigiVec.end(); iter++ ) {
378 const MdcDigi* aDigi = *iter;
379 const Identifier id= aDigi->identify();
380 int layer = MdcID::layer(id);
381 int wire = MdcID::wire(id);
382 m_digiMap[layer][wire] = aDigi;
383 }
384
385 //calc. doca of skipped hits
386 HitRefVec::iterator itHitSkipped = skippedHits.begin();
387 for (;itHitSkipped!=skippedHits.end();++itHitSkipped){
388 RecMdcHit* h = *itHitSkipped;
389 Identifier id(h->getMdcId());
390 int layer = MdcID::layer(id);
391 int wire = MdcID::wire(id);
392 double fltLen = 0.;
393 HepVector helix = newFit->helix(fltLen).params();
394 HepSymMatrix err = newFit->helix(fltLen).covariance();
395 double docaFltLen = h->getFltLen();
396 //double doca = m_mdcUtilitySvc->docaPatPar(layer,wire,docaFltLen,helix,err);//yzhang 2012-07-19
397 double doca = m_mdcUtilitySvc->docaPatPar(layer,wire,helix,err);
398 double newDoca = fabs(doca);
399 int flagLR = h->getFlagLR();
400 if ( flagLR == 0 ){ newDoca = -fabs(doca); }//left driftDist <0
401 h->setDoca(newDoca);
402 if(m_debug>=3)std::cout<<"("<<layer<<","<<wire<<") new doca "<<newDoca<<" resid "<<fabs(fabs(newDoca)-fabs(h->getDriftDistLeft()))<<std::endl;
403 if(m_debug>=3&&fabs(fabs(newDoca)-fabs(h->getDriftDistLeft()))>0.02) std::cout<<" bad resid "<<fabs(fabs(newDoca)-fabs(h->getDriftDistLeft()))<<" doca "<<newDoca<<" dd "<<h->getDriftDistLeft()<<std::endl;
404
405 //calc drift dist of skipped hit
406 MdcHit *thehit = new MdcHit(m_digiMap[layer][wire], m_gm);
407 thehit->setCosmicFit(true);
408 thehit->setCalibSvc(m_mdcCalibFunSvc);
409 thehit->setCountPropTime(m_countPropTime);
410 m_hitCol->push_back(thehit);
411
412 HepPoint3D poca; Hep3Vector tempDir;
413 getInfo(helix,0,poca,tempDir);
414 if(m_debug>3)std::cout<<"track poca "<<poca<<std::endl;
415 HepPoint3D pos; Hep3Vector dir;
416 getInfo(helix,docaFltLen,pos,dir);
417 //std::cout<<" pos "<<pos<<" dir "<<dir<<" "<<std::endl;
418 if(pos.y()>poca.y()){
419 int wireAmb = patAmbig(flagLR);
420 //double tofold = (docaFltLen/Constants::c + m_bunchT0)*1.e-9;
421 double tof = docaFltLen/Constants::c + (m_bunchT0)*1.e-9;
422 //std::cout<<layer<<" "<<wire<<" tofold "<<tofold*1.e9-m_bunchT0 <<" tof "<<tof*1.e9-m_bunchT0<<" diff "<<(tof-tofold)*1.e9<<std::endl;
423 //tof = tofold;
424 double eAngle = EntranceAngle(dir.phi() - thehit->phi(pos.z()));
425 double dAngle = Constants::pi/2 - 1.*dir.theta();
426 double z = pos.z();
427 double dist = thehit->driftDist(tof, wireAmb, eAngle, dAngle, z);
428 double sigma = thehit->sigma(dist,wireAmb,eAngle,dAngle,z);
429
430 if(m_debug>3&&fabs(fabs(h->getDoca())-fabs(dist))>0.02)
431 std::cout<<"("<<layer<<","<<wire<<") old doca "<<h->getDoca()<<" old ddr "<<h->getDriftDistRight()
432 <<" new dd "<<dist <<" newAmbig "<<wireAmb
433 //<<" tof "<<tof <<" z "<<z <<" eAngle "<<eAngle <<" dTime "<<thehit->driftTime(tof,z)
434 <<" old sigma "<<h->getErrDriftDistLeft()
435 <<" new sigma "<<sigma<<" "<<std::endl;
436 h->setDriftDistLeft(-1.*abs(dist));
437 h->setDriftDistRight(abs(dist));
438 h->setErrDriftDistLeft(sigma);
439 h->setErrDriftDistRight(sigma);
440 }
441 }
442
443 //-------------------------------------------------------
444 //Store TDS
445 //-------------------------------------------------------
446 store(newTrack, skippedHits);//newTrack have been deleted
447
448 setFilterPassed(true);
449 m_nSewed++;
450 }// end if newFit
451
452 if(m_debug) { m_mdcPrintSvc->printRecMdcTrackCol(); }
453
454 return StatusCode::SUCCESS;
455}
double abs(const EvtComplex &c)
Definition: EvtComplex.hh:212
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
ObjectVector< MdcHit > MdcHitCol
Definition: MdcHit.h:129
std::vector< MdcDigi * > MdcDigiVec
float Mdcxprobab(int &ndof, float &chisq)
Definition: Mdcxprobab.cxx:4
SmartRefVector< RecMdcHit > HitRefVec
Definition: RecMdcTrack.h:26
static const double pi
Definition: Constants.h:38
static const double twoPi
Definition: Constants.h:39
static const double c
Definition: Constants.h:43
const double chi2() const
Definition: DstMdcTrack.h:66
const int trackId() const
Definition: DstMdcTrack.h:52
const int ndof() const
Definition: DstMdcTrack.h:67
const HepVector helix() const
......
Definition: MdcHit.h:44
void setCalibSvc(const MdcCalibFunSvc *calibSvc)
Definition: MdcHit.cxx:136
void setCosmicFit(const bool cosmicfit)
Definition: MdcHit.h:87
double phi() const
Definition: MdcHit.h:75
double sigma(double, int, double, double, double) const
Definition: MdcHit.cxx:184
void setCountPropTime(const bool count)
Definition: MdcHit.h:86
double driftDist(double, int, double, double, double) const
Definition: MdcHit.cxx:156
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
Definition: MdcID.cxx:49
static int wire(const Identifier &id)
Definition: MdcID.cxx:54
void printRecMdcTrackCol() const
HepVector besPar2PatPar(const HepVector &helixPar) const
double docaPatPar(int layer, int cell, const HepVector helixPat, const HepSymMatrix errMatPat, bool passCellRequired=true, bool doSag=true) const
void MdcxHitsToHots(HepVector &helix, TrkRecoTrk *trk, HitRefVec &hits, HitRefVec &skiped)
void store(TrkRecoTrk *tk, HitRefVec &skip)
MdcDigiVec & getMdcDigiVec(uint32_t control=0)
virtual Identifier identify() const
Definition: RawData.cxx:15
const int getFlagLR(void) const
Definition: RecMdcHit.h:47
void setErrDriftDistRight(double erddr)
Definition: RecMdcHit.h:63
void setErrDriftDistLeft(double erddl)
Definition: RecMdcHit.h:62
void setDriftDistLeft(double ddl)
Definition: RecMdcHit.h:60
const Identifier getMdcId(void) const
Definition: RecMdcHit.h:49
void setDoca(double doca)
Definition: RecMdcHit.h:71
const double getDriftDistRight(void) const
Definition: RecMdcHit.h:43
const double getDriftDistLeft(void) const
Definition: RecMdcHit.h:42
const double getFltLen(void) const
Definition: RecMdcHit.h:56
const double getDoca(void) const
Definition: RecMdcHit.h:53
void setDriftDistRight(double ddr)
Definition: RecMdcHit.h:61
const double getErrDriftDistLeft(void) const
Definition: RecMdcHit.h:44
const HitRefVec getVecHits(void) const
Definition: RecMdcTrack.h:67
const HepVector & params() const
const HepSymMatrix & covariance() const
Definition: TrkFit.h:23
virtual TrkExchangePar helix(double fltL) const =0
TrkErrCode fit()
Definition: TrkHitList.cxx:59
TrkHitList * hits()
Definition: TrkRecoTrk.h:107
virtual void print(std::ostream &) const
Definition: TrkRecoTrk.cxx:166
const TrkFit * fitResult() const
Definition: TrkRecoTrk.cxx:387
bool setFlipAndDrop(TrkRecoTrk &, bool allowFlips, bool allowDrops) const
TrkRecoTrk * makeTrack(const TrkExchangePar &helix, const double chi2, const TrkContext &, double trackT0) const

◆ finalize()

StatusCode MdcxCosmicSewer::finalize ( )

Definition at line 457 of file MdcxCosmicSewer.cxx.

457 {
458 MsgStream log(msgSvc(), name());
459 log << MSG::INFO << "in finalize()" << endreq;
460 std::cout<<name()<<" after sewed got "<<m_nSewed<<" cosmic tracks"<<std::endl;
461 return StatusCode::SUCCESS;
462}

◆ initialize()

StatusCode MdcxCosmicSewer::initialize ( )

Definition at line 93 of file MdcxCosmicSewer.cxx.

93 {
94 MsgStream log(msgSvc(), name());
95 StatusCode sc;
96 log << MSG::INFO << "in initialize()" << endreq;
97
98 i_evt=0;
99 m_nSewed = 0;
100
101 //Initailize magnetic filed
102 sc = service ("MagneticFieldSvc",m_pIMF);
103 if(sc!=StatusCode::SUCCESS) {
104 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
105 }
106 m_bfield = new BField(m_pIMF);
107 log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
108 m_context = new TrkContextEv(m_bfield);
109
110
111 // Get MdcCalibFunSvc
112 IMdcCalibFunSvc* imdcCalibSvc;
113 sc = service ("MdcCalibFunSvc", imdcCalibSvc);
114 m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
115 if ( sc.isFailure() ){
116 log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endreq;
117 return StatusCode::FAILURE;
118 }
119
120 // Get RawDataProviderSvc
121 IRawDataProviderSvc* iRawDataProvider;
122 sc = service ("RawDataProviderSvc", iRawDataProvider);
123 if ( sc.isFailure() ){
124 log << MSG::FATAL << "Could not load RawDataProviderSvc!" << endreq;
125 return StatusCode::FAILURE;
126 }
127 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*>(iRawDataProvider);
128
129 //Initailize MdcUtilitySvc
130 IMdcUtilitySvc * imdcUtilitySvc;
131 sc = service ("MdcUtilitySvc", imdcUtilitySvc);
132 m_mdcUtilitySvc = dynamic_cast<MdcUtilitySvc*> (imdcUtilitySvc);
133 if ( sc.isFailure() ){
134 log << MSG::FATAL << "Could not load MdcUtilitySvc!" << endreq;
135 return StatusCode::FAILURE;
136 }
137
138 //Initailize MdcPrintSvc
139 IMdcPrintSvc* imdcPrintSvc;
140 sc = service ("MdcPrintSvc", imdcPrintSvc);
141 m_mdcPrintSvc = dynamic_cast<MdcPrintSvc*> (imdcPrintSvc);
142 if ( sc.isFailure() ){
143 log << MSG::FATAL << "Could not load MdcPrintSvc!" << endreq;
144 return StatusCode::FAILURE;
145 }
146
147 if(m_hist){
148 //book ntuple
149 NTuplePtr ntCsmcSew(ntupleSvc(), "MdcxReco/csmc");
150 if ( ntCsmcSew ) { m_xtupleCsmcSew = ntCsmcSew;}
151 else {
152 m_xtupleCsmcSew = ntupleSvc()->book ("MdcxReco/csmc", CLID_ColumnWiseTuple, "MdcxReco reconsturction results");
153 if ( m_xtupleCsmcSew ) {
154 sc = m_xtupleCsmcSew->addItem ("dD0", m_csmcD0);
155 sc = m_xtupleCsmcSew->addItem ("dPhi0", m_csmcPhi0);
156 sc = m_xtupleCsmcSew->addItem ("dZ0", m_csmcZ0);
157 sc = m_xtupleCsmcSew->addItem ("dOmega", m_csmcOmega);
158 sc = m_xtupleCsmcSew->addItem ("dPt", m_csmcPt);
159 sc = m_xtupleCsmcSew->addItem ("dTanl", m_csmcTanl);
160 }else{
161 log << MSG::FATAL << "Could not book MdcxReco/csmc!" << endreq;
162 return StatusCode::FAILURE;
163 }
164 }//end of book
165 }
166
167 return StatusCode::SUCCESS;
168}
INTupleSvc * ntupleSvc()

◆ MdcxHitsToHots()

void MdcxCosmicSewer::MdcxHitsToHots ( HepVector &  helix,
TrkRecoTrk trk,
HitRefVec hits,
HitRefVec skiped 
)

Definition at line 464 of file MdcxCosmicSewer.cxx.

464 {
465
466 TrkHitList* trkHitList = newTrack->hits();
467 //store new MdcHit
468 SmartDataPtr<MdcHitCol> m_hitCol(eventSvc(), "/Event/Hit/MdcHitCol");
469 if (!m_hitCol){
470 m_hitCol= new MdcHitCol;
471 StatusCode sc = eventSvc()->registerObject("/Event/Hit/MdcHitCol",m_hitCol);
472 if(!sc.isSuccess()) {
473 std::cout<< " Could not register MdcHitCol" <<endreq;
474 return;
475 }
476 }
477
478 //get MdcDigi pointer
479 uint32_t getDigiFlag = 0;
480 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
481 if (0 == mdcDigiVec.size()){
482 std::cout << " No hits in MdcDigiVec" << std::endl;
483 return;
484 }
485 const MdcDigi* m_digiMap[43][288];
486 MdcDigiVec::iterator iter = mdcDigiVec.begin();
487 for (; iter != mdcDigiVec.end(); iter++ ) {
488 const MdcDigi* aDigi = *iter;
489 const Identifier id= aDigi->identify();
490 int layer = MdcID::layer(id);
491 int wire = MdcID::wire(id);
492 m_digiMap[layer][wire] = aDigi;
493 }
494
495
496 HepPoint3D poca; Hep3Vector tempDir;
497 getInfo(helix,0,poca,tempDir);
498 if(m_debug>3)std::cout<<"track poca "<<poca<<std::endl;
499 // new MdcRecoHitOnTracks
500 HitRefVec::iterator itHit = recMdcHits.begin();
501 for (;itHit!=recMdcHits.end();++itHit){
502 RecMdcHit* h = *itHit;
503 int layer = MdcID::layer(h->getMdcId());
504 int wire = MdcID::wire(h->getMdcId());
505 // new fltLen and ambig
506 double newFltLen = h->getFltLen();
507 HepPoint3D pos; Hep3Vector dir;
508 getInfo(helix,h->getFltLen(),pos,dir);
509 int newAmbig = patAmbig(h->getFlagLR());
510 if(pos.y()>poca.y()){//yzhang TEMP 2012-07-17
511 newFltLen *= -1.;
512 newAmbig *= -1;
513 if(m_debug>3) std::cout<<" Up track, change sign of fltLen "<<std::endl;
514 }
515 int newFlagLR = bes3FlagLR(newAmbig);
516 // calc. position of this point
517 if(m_cosmicSewSkip && layer<20){
518 RecMdcHit* newSkippedHit = new RecMdcHit();
519 newSkippedHit->setDriftDistLeft(h->getDriftDistLeft());
520 newSkippedHit->setDriftDistRight(h->getDriftDistRight());
521 newSkippedHit->setErrDriftDistLeft(h->getErrDriftDistLeft());
522 newSkippedHit->setErrDriftDistRight(h->getErrDriftDistRight());
523 newSkippedHit->setChisqAdd(h->getChisqAdd());
524 newSkippedHit->setFlagLR(newFlagLR);
525 newSkippedHit->setStat(h->getStat());
526 newSkippedHit->setMdcId(h->getMdcId());
527 newSkippedHit->setTdc(h->getTdc());
528 newSkippedHit->setAdc(h->getAdc());
529 newSkippedHit->setDriftT(h->getDriftT());
530 newSkippedHit->setDoca(999.);
531 //newSkippedHit->setDoca(h->getDoca());
532 newSkippedHit->setEntra(h->getEntra());
533 newSkippedHit->setZhit(h->getZhit());
534 newSkippedHit->setFltLen(newFltLen);
535 skippedHits.push_back(newSkippedHit);
536 }else{
537 MdcHit *thehit = new MdcHit(m_digiMap[layer][wire], m_gm);
538 thehit->setCosmicFit(true);
539 thehit->setCalibSvc(m_mdcCalibFunSvc);
540 thehit->setCountPropTime(m_countPropTime);
541 m_hitCol->push_back(thehit);
542
543 MdcRecoHitOnTrack temp(*thehit, newAmbig, m_bunchT0);//m_bunchT0 nano second here
544 MdcHitOnTrack* newhot = &temp;
545 newhot->setFltLen(newFltLen);
546 if(m_debug>3) std::cout<<name()<<" ("<<layer<<","<<wire<<") fltLen "<<h->getFltLen()<<" newFltLen "<<newFltLen<<" newAmbig "<<newAmbig<<" pos.y() "<<pos.y()<<std::endl;
547
548 trkHitList->appendHot(newhot);
549 }
550 }//end of loop hits
551}
const double getZhit(void) const
Definition: RecMdcHit.h:55
void setMdcId(Identifier mdcid)
Definition: RecMdcHit.h:67
void setFltLen(double fltLen)
Definition: RecMdcHit.h:74
const int getStat(void) const
Definition: RecMdcHit.h:48
const double getChisqAdd(void) const
Definition: RecMdcHit.h:46
const double getAdc(void) const
Definition: RecMdcHit.h:51
const double getTdc(void) const
Definition: RecMdcHit.h:50
const double getDriftT(void) const
Definition: RecMdcHit.h:52
const double getEntra(void) const
Definition: RecMdcHit.h:54
void setStat(int stat)
Definition: RecMdcHit.h:66
void setTdc(double tdc)
Definition: RecMdcHit.h:68
void setAdc(double adc)
Definition: RecMdcHit.h:69
void setFlagLR(int lr)
Definition: RecMdcHit.h:65
void setChisqAdd(double pChisq)
Definition: RecMdcHit.h:64
void setZhit(double zhit)
Definition: RecMdcHit.h:73
void setDriftT(double driftT)
Definition: RecMdcHit.h:70
const double getErrDriftDistRight(void) const
Definition: RecMdcHit.h:45
void setEntra(double entra)
Definition: RecMdcHit.h:72
TrkHitOnTrk * appendHot(const TrkHitOnTrk *theHot)
Definition: TrkHitList.cxx:91
void setFltLen(double f)
Definition: TrkHitOnTrk.h:147

Referenced by execute().

◆ store()

void MdcxCosmicSewer::store ( TrkRecoTrk tk,
HitRefVec skip 
)

Definition at line 553 of file MdcxCosmicSewer.cxx.

553 {
554 StatusCode sc;
555 MsgStream log(msgSvc(), name());
556
557 assert (aTrack != NULL);
558
559 //IDataManagerSvc *dataManSvc;
560 SmartIF<IDataManagerSvc> dataManSvc(eventSvc());
561 DataObject *aTrackCol;
562 eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
563 if(aTrackCol != NULL) {
564 dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
565 eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
566 }
567
568 DataObject *aRecHitCol;
569 eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
570 if(aRecHitCol != NULL) {
571 dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
572 eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
573 }
574
575 SmartDataPtr<RecMdcTrackCol> trackList(eventSvc(),EventModel::Recon::RecMdcTrackCol);
576 if (!trackList) {
577 RecMdcTrackCol *trackList = new RecMdcTrackCol;
578 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList);
579 if(!sc.isSuccess()) {
580 log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
581 }
582 }
583
584 SmartDataPtr<RecMdcHitCol> hitList(eventSvc(),EventModel::Recon::RecMdcHitCol);
585 if (!hitList) {
586 hitList = new RecMdcHitCol;
587 sc = eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList);
588 if(!sc.isSuccess()) {
589 log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
590 }
591 }
592 int trackId = trackList->size() - 1;
593
594 if(m_debug)std::cout<<__FILE__<<" "<<__LINE__<<" STORED"<< std::endl;
595 MdcTrack mdcTrack(aTrack);//aTrack have been deleted in ~MdcTrack()
596
597 //tkStat: 0,PatRec; 1,MdcxReco; 2,Tsf; 3,CurlFinder; -1,Combined cosmic
598 int tkStat = -1;
599 mdcTrack.storeTrack(trackId, trackList, hitList, tkStat);
600
601 //store hits on track
602 RecMdcTrackCol::iterator it = trackList->begin();
603 HitRefVec hl = (*it)->getVecHits();
604 HitRefVec hitRefVec;
605 HitRefVec::iterator itHit = hl.begin();
606 int ihl=0;
607 for (;itHit!=hl.end();++itHit){
608 ihl++;
609 RecMdcHit* h = *itHit;
610 SmartRef<RecMdcHit> refHit(h);
611 hitRefVec.push_back(refHit);
612 }
613
614 if(m_debug) std::cout<<__FILE__<<" "<<__LINE__<<"refHit stored "<< ihl<< std::endl;
615 //store skipped hits
616 HitRefVec::iterator itHitSkipped = skippedHits.begin();
617 int iskipped=0;
618 for (;itHitSkipped!=skippedHits.end();++itHitSkipped){
619 iskipped++;
620 (*itHitSkipped)->setTrkId(trackId);
621 hitList->push_back(*itHitSkipped);
622 SmartRef<RecMdcHit> refHitSkipped(*itHitSkipped);
623 hitRefVec.push_back(refHitSkipped);
624 }
625 if(m_debug) std::cout<<__FILE__<<" "<<__LINE__<<"skippedHits stored "<< iskipped<< std::endl;
626 (*it)->setVecHits(hitRefVec);
627
628 (*it)->setNhits((*it)->getVecHits().size());
629 (*it)->setNster(-1);
630 //omega==0, line fit :ndof = nhit-4
631 //omega>0, helix fit :ndof = nhit-5
632 if((*it)->helix(2)<0.00001)(*it)->setNdof((*it)->getNhits()-4);
633 else (*it)->setNdof((*it)->getNhits()-5);
634
635 if(m_debug) std::cout<<__FILE__<<" "<<__LINE__<<" stored "<< (*it)->getVecHits().size()<< std::endl;
636
637}

Referenced by execute().


The documentation for this class was generated from the following files: