171 {
172 MsgStream log(
msgSvc(), name());
173 log << MSG::INFO << "in execute()" << endreq;
174
175
176 setFilterPassed(false);
177
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
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
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
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++) {
222 double d0 = -tk->
helix(0);
223 double phi0 = tk->
helix(1)+ CLHEP::halfpi;
224 double omega = Bz*tk->
helix(2)/-333.567;
225 double z0 = tk->
helix(3);
226 double tanl = tk->
helix(4);
227
230
231 if(m_debug)std::cout<<__FILE__<<" "<<__LINE__<< "tk"<<iTk
232 <<"("<<d0<<","<<phi0<<"," <<","<<omega<<","<<z0<<","<<tanl<<")"<<std::endl;
233
234 if(iTk==0){
235 dParam[0]=d0;
236 dParam[1]=phi0;
237 dParam[2]=omega;
238 dParam[3]=z0;
239 dParam[4]=tanl;
240 }else{
241 dParam[0] += d0;
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();
255 chisum += bchisq;
256
257
258
259
260 }
261
262 if(besthelId == -999) return StatusCode::SUCCESS;
263
264
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
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
304 double d0 = -tk->
helix(0);
305 double phi0 = tk->
helix(1)+ CLHEP::halfpi;
306 double omega = Bz*tk->
helix(2)/-333.567;
307 double z0 = tk->
helix(3);
308 double tanl = tk->
helix(4);
309
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
317 if(m_lineFit){
318 newTrack = linefactory.
makeTrack(tt, chisum, *m_context, m_bunchT0*1.e-9);
320 }else{
321 newTrack = helixfactory.
makeTrack(tt, chisum, *m_context, m_bunchT0*1.e-9);
323 }
324
325
326
328 it = recMdcTrackCol->begin();
329 for (int iTk=0; it != recMdcTrackCol->end(); it++,iTk++) {
334 }
335
336
337
338
339
342
343
344 if (!newFit) {
345
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
354
355 if(m_lineFit){
357 }else{
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){
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;
376 const MdcDigi* m_digiMap[43][288];
377 MdcDigiVec::iterator
iter = mdcDigiVec.begin();
378 for (;
iter != mdcDigiVec.end();
iter++ ) {
383 m_digiMap[layer][wire] = aDigi;
384 }
385
386
387 HitRefVec::iterator itHitSkipped = skippedHits.begin();
388 for (;itHitSkipped!=skippedHits.end();++itHitSkipped){
393 double fltLen = 0.;
397
398 double doca = m_mdcUtilitySvc->
docaPatPar(layer,wire,helix,err);
399 double newDoca = fabs(doca);
401 if ( flagLR == 0 ){ newDoca = -fabs(doca); }
403 if(m_debug>=3)std::cout<<
"("<<layer<<
","<<wire<<
") new doca "<<newDoca<<
" resid "<<fabs(fabs(newDoca)-fabs(h->
getDriftDistLeft()))<<std::endl;
405
406
407 MdcHit *thehit =
new MdcHit(m_digiMap[layer][wire], m_gm);
411 m_hitCol->push_back(thehit);
412
414 getInfo(helix,0,poca,tempDir);
415 if(m_debug>3)std::cout<<"track poca "<<poca<<std::endl;
417 getInfo(helix,docaFltLen,pos,dir);
418
419 if(pos.y()>poca.y()){
420 int wireAmb = patAmbig(flagLR);
421
422 double tof = docaFltLen/
Constants::c + (m_bunchT0)*1.e-9;
423
424
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)
433 <<" new dd "<<dist <<" newAmbig "<<wireAmb
434
436 <<
" new sigma "<<
sigma<<
" "<<std::endl;
441 }
442 }
443
444
445
446
447 store(newTrack, skippedHits);
448
449 setFilterPassed(true);
450 m_nSewed++;
451 }
452
454
455 return StatusCode::SUCCESS;
456}
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