170 {
171 MsgStream log(
msgSvc(), name());
172 log << MSG::INFO << "in execute()" << endreq;
173
174
175 setFilterPassed(false);
176
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
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
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
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++) {
221 double d0 = -tk->
helix(0);
222 double phi0 = tk->
helix(1)+ CLHEP::halfpi;
223 double omega = Bz*tk->
helix(2)/-333.567;
224 double z0 = tk->
helix(3);
225 double tanl = tk->
helix(4);
226
229
230 if(m_debug)std::cout<<__FILE__<<" "<<__LINE__<< "tk"<<iTk
231 <<"("<<d0<<","<<phi0<<"," <<","<<omega<<","<<z0<<","<<tanl<<")"<<std::endl;
232
233 if(iTk==0){
234 dParam[0]=d0;
235 dParam[1]=phi0;
236 dParam[2]=omega;
237 dParam[3]=z0;
238 dParam[4]=tanl;
239 }else{
240 dParam[0] += d0;
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();
254 chisum += bchisq;
255
256
257
258
259 }
260
261 if(besthelId == -999) return StatusCode::SUCCESS;
262
263
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
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
303 double d0 = -tk->
helix(0);
304 double phi0 = tk->
helix(1)+ CLHEP::halfpi;
305 double omega = Bz*tk->
helix(2)/-333.567;
306 double z0 = tk->
helix(3);
307 double tanl = tk->
helix(4);
308
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
316 if(m_lineFit){
317 newTrack = linefactory.
makeTrack(tt, chisum, *m_context, m_bunchT0*1.e-9);
319 }else{
320 newTrack = helixfactory.
makeTrack(tt, chisum, *m_context, m_bunchT0*1.e-9);
322 }
323
324
325
327 it = recMdcTrackCol->begin();
328 for (int iTk=0; it != recMdcTrackCol->end(); it++,iTk++) {
333 }
334
335
336
337
338
341
342
343 if (!newFit) {
344
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
353
354 if(m_lineFit){
356 }else{
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){
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;
375 const MdcDigi* m_digiMap[43][288];
376 MdcDigiVec::iterator
iter = mdcDigiVec.begin();
377 for (;
iter != mdcDigiVec.end();
iter++ ) {
382 m_digiMap[layer][wire] = aDigi;
383 }
384
385
386 HitRefVec::iterator itHitSkipped = skippedHits.begin();
387 for (;itHitSkipped!=skippedHits.end();++itHitSkipped){
392 double fltLen = 0.;
396
397 double doca = m_mdcUtilitySvc->
docaPatPar(layer,wire,helix,err);
398 double newDoca = fabs(doca);
400 if ( flagLR == 0 ){ newDoca = -fabs(doca); }
402 if(m_debug>=3)std::cout<<
"("<<layer<<
","<<wire<<
") new doca "<<newDoca<<
" resid "<<fabs(fabs(newDoca)-fabs(h->
getDriftDistLeft()))<<std::endl;
404
405
406 MdcHit *thehit =
new MdcHit(m_digiMap[layer][wire], m_gm);
410 m_hitCol->push_back(thehit);
411
413 getInfo(helix,0,poca,tempDir);
414 if(m_debug>3)std::cout<<"track poca "<<poca<<std::endl;
416 getInfo(helix,docaFltLen,pos,dir);
417
418 if(pos.y()>poca.y()){
419 int wireAmb = patAmbig(flagLR);
420
421 double tof = docaFltLen/
Constants::c + (m_bunchT0)*1.e-9;
422
423
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)
432 <<" new dd "<<dist <<" newAmbig "<<wireAmb
433
435 <<
" new sigma "<<
sigma<<
" "<<std::endl;
440 }
441 }
442
443
444
445
446 store(newTrack, skippedHits);
447
448 setFilterPassed(true);
449 m_nSewed++;
450 }
451
453
454 return StatusCode::SUCCESS;
455}
ObjectVector< MdcHit > MdcHitCol
std::vector< MdcDigi * > MdcDigiVec
float Mdcxprobab(int &ndof, float &chisq)
SmartRefVector< RecMdcHit > HitRefVec
double bFieldNominal() const
static const double twoPi
const double chi2() const
const int trackId() const
const HepVector helix() const
......
void setCalibSvc(const MdcCalibFunSvc *calibSvc)
void setCosmicFit(const bool cosmicfit)
double sigma(double, int, double, double, double) const
void setCountPropTime(const bool count)
double driftDist(double, int, double, double, double) const
static int layer(const Identifier &id)
Values of different levels (failure returns 0)
static int wire(const Identifier &id)
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 clearRecMdcTrackHit()
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
const int getFlagLR(void) const
void setErrDriftDistRight(double erddr)
void setErrDriftDistLeft(double erddl)
void setDriftDistLeft(double ddl)
const Identifier getMdcId(void) const
void setDoca(double doca)
const double getDriftDistRight(void) const
const double getDriftDistLeft(void) const
const double getFltLen(void) const
const double getDoca(void) const
void setDriftDistRight(double ddr)
const double getErrDriftDistLeft(void) const
const HitRefVec getVecHits(void) const
const HepVector & params() const
const HepSymMatrix & covariance() const
virtual TrkExchangePar helix(double fltL) const =0
virtual void print(std::ostream &) const
const TrkFit * fitResult() const
bool setFlipAndDrop(TrkRecoTrk &, bool allowFlips, bool allowDrops) const
TrkRecoTrk * makeTrack(const TrkExchangePar &helix, const double chi2, const TrkContext &, double trackT0) const