BOSS 7.1.0
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 59 of file MdcxCosmicSewer.cxx.

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

◆ ~MdcxCosmicSewer()

MdcxCosmicSewer::~MdcxCosmicSewer ( )
virtual

Definition at line 78 of file MdcxCosmicSewer.cxx.

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

Member Function Documentation

◆ beginRun()

StatusCode MdcxCosmicSewer::beginRun ( )

Definition at line 83 of file MdcxCosmicSewer.cxx.

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

◆ clearRecMdcTrackHit()

void MdcxCosmicSewer::clearRecMdcTrackHit ( )

Definition at line 640 of file MdcxCosmicSewer.cxx.

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

Referenced by execute().

◆ dumpTdsTrack()

void MdcxCosmicSewer::dumpTdsTrack ( RecMdcTrackCol trackList)

◆ execute()

StatusCode MdcxCosmicSewer::execute ( )

Definition at line 171 of file MdcxCosmicSewer.cxx.

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

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

◆ initialize()

StatusCode MdcxCosmicSewer::initialize ( )

Definition at line 94 of file MdcxCosmicSewer.cxx.

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

◆ MdcxHitsToHots()

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

Definition at line 465 of file MdcxCosmicSewer.cxx.

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

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

Referenced by execute().


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