loon: CandTrackHandle.cxx Source File

00001 
00002 // $Id: CandTrackHandle.cxx,v 1.51 2007/03/23 16:00:41 musser Exp $
00003 //
00004 // CandTrackHandle.cxx
00005 //
00006 // Begin_Html<img src="../../pedestrians.gif" align=center>
00007 // <a href="../source_warning.html">Warning for beginners</a>.<br> 
00008 //
00009 // CandTrackHandle is the specialized access handle to CandTrack.
00010 //
00011 // Each concrete CandHandle must define a DupHandle function.
00012 //
00013 // Author:  R. Lee
00014 //
00015 // Also see <a href="../../root_crib/index.html">The ROOT Crib</a> and 
00016 // <a href="../CandTrack.html"> CandTrack Classes</a> (part of
00017 // <a href="../index.html">The MINOS Class User Guide</a>)End_Html
00019 
00020 #include "TMath.h"
00021 
00022 #include <cassert>
00023 #include <cmath>
00024 #include <iostream>
00025 
00026 #include "MessageService/MsgService.h"
00027 #include "RecoBase/CandSliceHandle.h"
00028 #include "RecoBase/CandStripHandle.h"
00029 #include "RecoBase/CandStrip.h"
00030 #include "RecoBase/CandTrackHandle.h"
00031 #include "RecoBase/CandShowerHandle.h"
00032 #include "Validity/VldContext.h"
00033 #include "Algorithm/AlgConfig.h"
00034 #include "UgliGeometry/UgliGeomHandle.h"
00035 #include "UgliGeometry/UgliScintPlnHandle.h"
00036 #include "RecoBase/CandTrack.h"
00037 
00038 ClassImp(CandTrackHandle)
00039 
00040 //______________________________________________________________________
00041 CVSID("$Id: CandTrackHandle.cxx,v 1.51 2007/03/23 16:00:41 musser Exp $");
00042 
00043 //______________________________________________________________________
00044 CandTrackHandle::CandTrackHandle()
00045 {
00046 }
00047 
00048 //______________________________________________________________________
00049 CandTrackHandle::CandTrackHandle(const CandTrackHandle &cdh) :
00050   CandRecoHandle(cdh)
00051 {
00052 }
00053 
00054 //______________________________________________________________________
00055 CandTrackHandle::CandTrackHandle(CandTrack *cd) :
00056   CandRecoHandle(cd)
00057 {
00058 }
00059 
00060 //______________________________________________________________________
00061 CandTrackHandle::~CandTrackHandle()
00062 {
00063 }
00064 
00065 //______________________________________________________________________
00066 CandTrackHandle *CandTrackHandle::DupHandle() const
00067 {
00068    return (new CandTrackHandle(*this));
00069 }
00070 
00071 
00072 //______________________________________________________________________
00073 void CandTrackHandle::Trace(const char *c) const
00074 {
00075   MSG("Cand", Msg::kDebug)
00076     << "**********Begin CandTrackHandle::Trace(\"" << c << "\")" << endl
00077            << "Information from CandTrackHandle's CandHandle: " << endl;
00078   CandHandle::Trace(c);
00079   MSG("Cand", Msg::kDebug)
00080      << "**********End CandTrackHandle::Trace(\"" << c << "\")" << endl;
00081 }
00082 
00083 
00084 NavKey CandTrackHandle::KeyFromSlice(const CandTrackHandle *reco)
00085 {
00086   if (reco->GetCandSlice()) {
00087     return static_cast<Int_t>(reco->GetCandSlice()->GetUidInt());
00088   }
00089   return 0;
00090 
00091 }
00092 
00093 void CandTrackHandle::SetTrackPointError(Int_t plane, Float_t tpos)
00094 {
00095   CandTrack * track=dynamic_cast<CandTrack*>(GetOwnedCandBase());
00096   track->fTPosError[plane] = tpos;
00097 }
00098 
00099 void CandTrackHandle::SetU(Int_t plane, Float_t tpos) 
00100 {
00101   CandTrack *track = dynamic_cast<CandTrack*>(GetOwnedCandBase());
00102   track->fUPos[plane] = tpos;
00103 }
00104 
00105 void CandTrackHandle::SetV(Int_t plane, Float_t tpos) 
00106 {
00107   CandTrack *track = dynamic_cast<CandTrack*>(GetOwnedCandBase());
00108   track->fVPos[plane] = tpos;
00109 }
00110 
00111 void CandTrackHandle::SetdS(Int_t plane, Float_t ds)
00112 {
00113   CandTrack *track = dynamic_cast<CandTrack*>(GetOwnedCandBase());
00114   track->fdS[plane] = ds;
00115 }
00116 
00117 void CandTrackHandle::SetRange(Int_t plane, Float_t range)
00118 {
00119   CandTrack *track = dynamic_cast<CandTrack*>(GetOwnedCandBase());
00120   track->fRange[plane] = range;
00121 }
00122 
00123 void CandTrackHandle::SetT(Int_t plane, StripEnd::StripEnd_t stripend_t, Double_t time)
00124 {
00125   CandTrack *track = dynamic_cast<CandTrack*>(GetOwnedCandBase());
00126   switch (stripend_t) {
00127   case StripEnd::kNegative:
00128     track->fTime[0][plane] = time;
00129     break;
00130   case StripEnd::kPositive:
00131     track->fTime[1][plane] = time;
00132     break;
00133   default:
00134     track->fTime[0][plane] = time;
00135     break;
00136   }
00137 }
00138 
00139 void CandTrackHandle::ClearMaps()
00140 {
00141   CandTrack *track = dynamic_cast<CandTrack*>(GetOwnedCandBase());
00142   track->fUPos.clear();
00143   track->fVPos.clear();
00144   track->fTime[0].clear();
00145   track->fTime[1].clear();
00146   track->fdS.clear();
00147   track->fRange.clear();
00148 }
00149 
00150 void CandTrackHandle::ClearUVT()
00151 {
00152   ClearMaps();
00153 }
00154 
00155 Bool_t CandTrackHandle::IsTPosValid(Int_t plane) const
00156 {
00157   const CandTrack *track = dynamic_cast<const CandTrack*>(GetCandBase());
00158   if ((track->fUPos).count(plane) && (track->fVPos).count(plane)) {
00159     return kTRUE;
00160   }
00161   return kFALSE;
00162 }
00163 
00164 Float_t CandTrackHandle::GetTrackPointError(Int_t plane) const
00165 {
00166   const CandTrack *track = dynamic_cast<const CandTrack*>(GetCandBase());
00167   if((track->fTPosError).count(plane)){
00168     return track->fTPosError[plane];
00169   }
00170   return -99999.;
00171 
00172 }
00173 
00174 Float_t CandTrackHandle::GetU(Int_t plane) const
00175 {
00176   const CandTrack *track = dynamic_cast<const CandTrack*>(GetCandBase());
00177   if ((track->fUPos).count(plane)) {
00178     return track->fUPos[plane];
00179   }
00180   return -99999.;
00181 }
00182 
00183 Float_t CandTrackHandle::GetV(Int_t plane) const
00184 {
00185   const CandTrack *track = dynamic_cast<const CandTrack*>(GetCandBase());
00186   if ((track->fVPos).count(plane)) {
00187     return track->fVPos[plane];
00188   }
00189   return -99999.;
00190 }
00191 
00192 Float_t CandTrackHandle::GetZ(Int_t plane) const
00193 {
00194   TIter stripItr(GetDaughterIterator());
00195   while (CandStripHandle *strip = dynamic_cast<CandStripHandle*>(stripItr())) {
00196     if (strip->GetPlane()==plane) {
00197       return strip->GetZPos();
00198     }
00199   }
00200   return -1.;
00201 }
00202 
00203 Double_t CandTrackHandle::GetT(Int_t plane,StripEnd::StripEnd_t stripend_t) const
00204 {
00205   const CandTrack *track = dynamic_cast<const CandTrack*>(GetCandBase());
00206   if (stripend_t==StripEnd::kNegative) {
00207     if ((track->fTime[0]).count(plane)) {
00208       return track->fTime[0][plane];
00209     }
00210     return -99999.;
00211   }
00212   if (stripend_t==StripEnd::kPositive) {
00213     if ((track->fTime[1]).count(plane)) {
00214       return track->fTime[1][plane];
00215     }
00216     return -99999.;
00217   }
00218   if ((track->fTime[0]).count(plane) && (track->fTime[1]).count(plane)) {
00219     return min(track->fTime[0][plane],track->fTime[1][plane]);
00220   }
00221   else if ((track->fTime[0]).count(plane)) {
00222     return track->fTime[0][plane];
00223   }
00224   else if ((track->fTime[1]).count(plane)) {
00225     return track->fTime[1][plane];
00226   }
00227   return -99999.;
00228 }
00229 
00230 Double_t CandTrackHandle::GetT(StripEnd::StripEnd_t stripend_t,Int_t plane) const
00231 {
00232   return GetT(plane,stripend_t);
00233 }
00234 
00235 Double_t CandTrackHandle::GetT(Int_t plane) const
00236 {
00237   return GetT(plane,StripEnd::kWhole);
00238 }
00239 
00240 
00241 Float_t CandTrackHandle::GetdS(Int_t plane) const
00242 {
00243   if (plane<0) {
00244     plane = GetVtxPlane();
00245   }
00246   const CandTrack *track = dynamic_cast<const CandTrack*>(GetCandBase());
00247   if ((track->fdS).count(plane)) {
00248     return track->fdS[plane];
00249   }
00250   return -1.;
00251 }
00252 
00253 Float_t CandTrackHandle::GetRange(Int_t plane) const
00254 {
00255   if (plane<0) {
00256     plane = GetVtxPlane();
00257   }
00258   const CandTrack *track = dynamic_cast<const CandTrack*>(GetCandBase());
00259   if ((track->fRange).count(plane)) {
00260     return track->fRange[plane];
00261   }
00262   return -1.;
00263 }
00264 
00265 
00266 Double_t CandTrackHandle::GetScore() const
00267 {
00268   Double_t score = 0.;
00269   score = (Double_t)(GetNDaughters());
00270   Int_t dbegpln = GetBegPlane(PlaneView::kU)-GetBegPlane(PlaneView::kV);
00271   Int_t dendpln = GetEndPlane(PlaneView::kU)-GetEndPlane(PlaneView::kV);
00272   Double_t dbegpln2 = (Double_t)(dbegpln*dbegpln);
00273   Double_t dendpln2 = (Double_t)(dendpln*dendpln);
00274   score -= (dbegpln2+dendpln2);
00275   return score;
00276 }
00277 
00278 
00279 Int_t CandTrackHandle::GetNTrackPlane(PlaneView::PlaneView_t planeview_t) const
00280 {
00281   CandStripHandleItr stripItr(GetDaughterIterator());
00282   CandStripHandleKeyFunc *stripKf = stripItr.CreateKeyFunc();
00283   stripKf->SetFun(CandStripHandle::KeyFromPlane);
00284   stripItr.GetSet()->AdoptSortKeyFunc(stripKf);
00285   stripKf = 0;
00286 
00287   Int_t plane=0;
00288   Int_t oldplane=0;
00289   CandStripHandle *strip;
00290   while ((strip = dynamic_cast<CandStripHandle*>(stripItr()))) {
00291     if (IsInShower(strip)<=1) {
00292       PlaneView::PlaneView_t planeview = strip->GetPlaneView();
00293       if (planeview!=PlaneView::kA && planeview!=PlaneView::kB &&
00294           (planeview_t==PlaneView::kUnknown || planeview==planeview_t)) {
00295         if (!plane || strip->GetPlane()!=oldplane) {
00296           plane++;
00297         }
00298         oldplane = strip->GetPlane();
00299       }
00300     }
00301   }
00302   return plane;
00303 }
00304 
00305 
00306 Double_t CandTrackHandle::GetMomentum() const
00307 {
00308   return dynamic_cast<const CandTrack*>(GetCandBase())->fMomentum;
00309 }
00310 
00311 void CandTrackHandle::SetMomentum(Double_t momentum)
00312 {
00313   dynamic_cast<CandTrack*>(GetOwnedCandBase())->fMomentum = momentum;
00314 }
00315 
00316 Int_t CandTrackHandle::IsInShower(CandStripHandle *striphandle) const
00317 {
00318 
00319   const CandTrack *track = dynamic_cast<const CandTrack*>(GetCandBase());
00320   const CandStrip *strip = dynamic_cast<const CandStrip*>(striphandle->GetCandBase());
00321   if (track->fInShower.count(strip)>0) {
00322     return track->fInShower[strip];
00323   }
00324   return 0;
00325 }
00326 
00327 void CandTrackHandle::SetInShower(CandStripHandle *striphandle, Int_t ival)
00328 {
00329   CandTrack *track = dynamic_cast<CandTrack*>(GetOwnedCandBase());
00330   const CandStrip *strip = dynamic_cast<const CandStrip*>(striphandle->GetCandBase());
00331   track->fInShower[strip] = ival;
00332 }
00333 
00334 
00335 
00336 void CandTrackHandle::SetVtxTrace(Double_t dvar)
00337 {
00338   dynamic_cast<CandTrack *>(GetOwnedCandBase())->fVtxTrace = dvar;
00339 }
00340 
00341 void CandTrackHandle::SetVtxTraceZ(Double_t dvar)
00342 {
00343   dynamic_cast<CandTrack *>(GetOwnedCandBase())->fVtxTraceZ = dvar;
00344 }
00345 
00346 void CandTrackHandle::SetVtxnActiveUpstream(int ivar)
00347 {
00348   dynamic_cast<CandTrack *>(GetOwnedCandBase())->fVtxnActiveUpstream = ivar;
00349 }
00350 
00351 
00352 void CandTrackHandle::SetEndTrace(Double_t dvar)
00353 {
00354   dynamic_cast<CandTrack *>(GetOwnedCandBase())->fEndTrace = dvar;
00355 }
00356 
00357 void CandTrackHandle::SetEndTraceZ(Double_t dvar)
00358 {
00359   dynamic_cast<CandTrack *>(GetOwnedCandBase())->fEndTraceZ = dvar;
00360 }
00361 
00362 void CandTrackHandle::SetEndnActiveDownstream(int ivar)
00363 {
00364   dynamic_cast<CandTrack *>(GetOwnedCandBase())->fEndnActiveDownstream = ivar;
00365 }
00366 
00367 void CandTrackHandle::SetVtxDistToEdge(Double_t dvar)
00368 {
00369   dynamic_cast<CandTrack *>(GetOwnedCandBase())->fVtxDistToEdge = dvar;
00370 }
00371 
00372 void CandTrackHandle::SetEndDistToEdge(Double_t dvar)
00373 {
00374   dynamic_cast<CandTrack *>(GetOwnedCandBase())->fEndDistToEdge = dvar;
00375 }
00376 
00377 
00378 Double_t CandTrackHandle::GetVtxTrace() const
00379 {
00380   return dynamic_cast<const CandTrack *>(GetCandBase())->fVtxTrace;
00381 }
00382 
00383 Double_t CandTrackHandle::GetVtxTraceZ() const
00384 {
00385   return dynamic_cast<const CandTrack *>(GetCandBase())->fVtxTraceZ;
00386 }
00387 
00388 Int_t CandTrackHandle::GetVtxnActiveUpstream() const
00389 {
00390   return dynamic_cast<const CandTrack *>(GetCandBase())->fVtxnActiveUpstream;
00391 }
00392 
00393 Double_t CandTrackHandle::GetEndTrace() const
00394 {
00395   return dynamic_cast<const CandTrack *>(GetCandBase())->fEndTrace;
00396 }
00397 
00398 Double_t CandTrackHandle::GetEndTraceZ() const
00399 {
00400   return dynamic_cast<const CandTrack *>(GetCandBase())->fEndTraceZ;
00401 }
00402 
00403 Int_t CandTrackHandle::GetEndnActiveDownstream() const
00404 {
00405   return dynamic_cast<const CandTrack *>(GetCandBase())->fEndnActiveDownstream;
00406 }
00407 
00408 Double_t CandTrackHandle::GetVtxDistToEdge() const
00409 {
00410   return dynamic_cast<const CandTrack *>(GetCandBase())->fVtxDistToEdge;
00411 }
00412 
00413 Double_t CandTrackHandle::GetEndDistToEdge() const
00414 {
00415   return dynamic_cast<const CandTrack *>(GetCandBase())->fEndDistToEdge;
00416 }
00417 
00418 //_____________________________________________________________________
00419 Bool_t CandTrackHandle::BelongsWithTrack(CandTrackHandle * trk, 
00420                                           AlgConfig & ac, 
00421                                           const VldContext * vldcptr, 
00422                                           Double_t tolTPos2, Double_t tolZPos, Double_t tolTime){
00423 
00424   if(!trk)return false;
00425 
00426   Int_t pSMPlaneLast = ac.GetInt("SMPlaneLast");
00427   Int_t pSMPlaneFirst = ac.GetInt("SMPlaneFirst");
00428   VldContext vldc = *vldcptr;
00429   Double_t zGapSM=0;
00430   UgliGeomHandle ugh(vldc); 
00431   if(vldc.GetDetector() == Detector::kFar){  
00432     // calculate Z gap between SM for later use
00433     PlexPlaneId scintlastid(vldc.GetDetector(),pSMPlaneLast,kFALSE);
00434     PlexPlaneId scintfirstid(vldc.GetDetector(),pSMPlaneFirst,kFALSE);
00435     UgliScintPlnHandle scintlast = ugh.GetScintPlnHandle(scintlastid);
00436     UgliScintPlnHandle scintfirst = ugh.GetScintPlnHandle(scintfirstid);
00437     
00438     if (scintlast.IsValid() && scintfirst.IsValid()) {  
00439       zGapSM=scintfirst.GetZ0()-scintlast.GetZ0()-0.0594;
00440     } 
00441   } 
00442   Double_t dz =this->GetVtxZ()-trk->GetVtxZ();
00443   Double_t du= this->GetVtxU()-trk->GetVtxU();
00444   Double_t dv= this->GetVtxV()-trk->GetVtxV();  
00445   Double_t dt= this->GetVtxT()-trk->GetVtxT();
00446 
00447  //compensate dz for SM gap if necessary
00448   if (vldcptr->GetDetector()==Detector::kFar &&
00449       this->GetVtxPlane()<=pSMPlaneLast &&
00450       trk->GetVtxPlane()>=pSMPlaneFirst) dz+=zGapSM;
00451   else if (vldcptr->GetDetector()==Detector::kFar &&
00452            this->GetVtxPlane()>=pSMPlaneLast &&
00453            trk->GetVtxPlane()<=pSMPlaneFirst) dz-=zGapSM;
00454 
00455   if (du*du+dv*dv<tolTPos2 &&
00456       fabs(dz)<tolZPos && 
00457       fabs(dt)<tolTime) return true;
00458 
00459   return false;
00460 }
00461 //_____________________________________________________________________
00462 Bool_t CandTrackHandle::BelongsWithShower(CandShowerHandle * shw, 
00463                                           AlgConfig & ac, 
00464                                           const VldContext * vldcptr, 
00465                                           Double_t tolTPos2, Double_t tolZPos, Double_t tolTime){
00466   if(!shw)return false;
00467   
00468   Int_t pSMPlaneLast = ac.GetInt("SMPlaneLast");
00469   Int_t pSMPlaneFirst = ac.GetInt("SMPlaneFirst");
00470   Double_t zGapSM=0;
00471   VldContext vldc = *vldcptr;
00472   UgliGeomHandle ugh(vldc); 
00473   // calculate Z gap between SM for later use
00474   if(vldc.GetDetector() == Detector::kFar){     
00475     PlexPlaneId scintlastid(vldc.GetDetector(),pSMPlaneLast,kFALSE);
00476     PlexPlaneId scintfirstid(vldc.GetDetector(),pSMPlaneFirst,kFALSE);
00477     UgliScintPlnHandle scintlast = ugh.GetScintPlnHandle(scintlastid);  
00478     UgliScintPlnHandle scintfirst = ugh.GetScintPlnHandle(scintfirstid);
00479     if (scintlast.IsValid() && scintfirst.IsValid()) {  
00480       zGapSM=scintfirst.GetZ0()-scintlast.GetZ0()-0.0594;
00481     } 
00482   }
00483 
00484   Float_t tolTPos=TMath::Sqrt(tolTPos2);
00485   Int_t fwdTrkPlane=min(this->GetVtxPlane(),this->GetEndPlane());
00486   Int_t bckTrkPlane=max(this->GetVtxPlane(),this->GetEndPlane());
00487   Int_t fwdShwPlane=min(shw->GetVtxPlane(),shw->GetEndPlane());
00488   Int_t bckShwPlane=max(shw->GetVtxPlane(),shw->GetEndPlane());
00489   Double_t dt= this->GetVtxT()-shw->GetVtxT();
00490   if(this->IsTPosValid(shw->GetVtxPlane())){
00491     dt = shw->GetVtxT()-this->GetT(shw->GetVtxPlane());
00492   }
00493   if( bckShwPlane>=fwdTrkPlane && fwdShwPlane<=bckTrkPlane && fabs(dt)<tolTime){   // trk overlaps in Z and Time
00494     MSG("RecoBase",Msg::kDebug) << " trk and shower overlap in Z and time " << endl;
00495     Bool_t matchu=false;
00496     Bool_t matchv=false;
00497     for(Int_t iloop=0;iloop<2;iloop++){
00498 
00499       // if on second loop through and have at least one matching view, loosen tolerance
00500       if(iloop >0 && (matchu || matchv))tolTPos=tolTPos*2;
00501       // if on second loop and no matches so far, give up
00502       if(iloop >0 && !matchu && !matchv) break;
00503 
00504       for (Int_t iplane=fwdShwPlane; iplane<=bckShwPlane;iplane++){  // look for trk passing through shower
00505         Double_t trkU = 0, trkV = 0;
00506         //if(matchu || matchv)tolTPos = tolTPos*2;
00507         PlexPlaneId plnid(GetVldContext()->GetDetector(),iplane,false);
00508         if(plnid.GetPlaneView()==PlaneView::kU){
00509           trkU=this->GetU(iplane);
00510           if(!this->IsTPosValid(iplane)){
00511             if(abs(iplane-this->GetVtxPlane())<abs(iplane-this->GetEndPlane())){
00512               UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(plnid);
00513               Double_t dz=scintpln.GetZ0()-this->GetVtxZ();
00514               trkU = (this->GetVtxU() + dz*this->GetVtxDirCosU()/this->GetVtxDirCosZ());
00515             }
00516             else{
00517               UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(plnid);
00518               Double_t dz=scintpln.GetZ0()-this->GetEndZ();
00519               trkU = (this->GetEndU() + dz*this->GetEndDirCosU()/this->GetEndDirCosZ());
00520             }
00521           }
00522           if(shw->GetNStrips(iplane)>0 && 
00523              trkU>=shw->GetMinU(iplane)-tolTPos && 
00524              trkU<=shw->GetMaxU(iplane)+tolTPos) {
00525             MSG("RecoBase", Msg::kDebug)<< " u match! (details follow)" << endl;
00526             matchu=true;  
00527           }
00528         }
00529 
00530         else if(plnid.GetPlaneView()==PlaneView::kV){
00531           trkV=this->GetV(iplane);
00532           if(!this->IsTPosValid(iplane)){
00533             if(abs(iplane-this->GetVtxPlane())<abs(iplane-this->GetEndPlane())){
00534               UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(plnid);
00535               Double_t dz=scintpln.GetZ0()-this->GetVtxZ();
00536               trkV = (this->GetVtxV() + dz*this->GetDirCosV()/this->GetDirCosZ());
00537             }
00538             else{
00539               UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(plnid);
00540               Double_t dz=scintpln.GetZ0()-this->GetEndZ();
00541               trkV= (this->GetEndV() + dz*this->GetEndDirCosV()/this->GetEndDirCosZ());
00542             }
00543           }
00544 
00545           if(shw->GetNStrips(iplane)>0 && 
00546              trkV>=shw->GetMinV(iplane)-tolTPos && 
00547              trkV<=shw->GetMaxV(iplane)+tolTPos) {
00548             MSG("RecoBase", Msg::kDebug)<< " v match! ((details follow)" << endl;
00549             matchv=true;
00550           }
00551         }
00552 
00553         MSG("RecoBase",Msg::kDebug) << " plane  " << iplane 
00554                                     << " trk u " << trkU 
00555                                     << " min/max U " 
00556                                     << shw->GetMinU(iplane) << "/" 
00557                                     << shw->GetMaxU(iplane)
00558                                     << " trk v " << trkV 
00559                                     << " min/max V " 
00560                                     << shw->GetMinV(iplane) << "/" 
00561                                     << shw->GetMaxV(iplane) << endl; 
00562 
00563         if(matchu && matchv) {
00564           MSG("RecoBase", Msg::kDebug)<< " 3D  match! " << endl;
00565           return true;
00566         }
00567       }
00568     }
00569   }
00570 
00571   if( shw->GetVtxZ()>min(this->GetVtxZ(),this->GetEndZ())-tolZPos  &&
00572       shw->GetVtxZ()<max(this->GetVtxZ(),this->GetEndZ())+tolZPos) {
00573     Double_t du= this->GetVtxU()-shw->GetVtxU();
00574     Double_t dv= this->GetVtxV()-shw->GetVtxV();        
00575     dt= this->GetVtxT()-shw->GetVtxT();
00576     if(this->IsTPosValid(shw->GetVtxPlane())){
00577       du = shw->GetVtxU()-this->GetU(shw->GetVtxPlane());
00578       dv = shw->GetVtxV()-this->GetV(shw->GetVtxPlane());
00579       dt = shw->GetVtxT()-this->GetT(shw->GetVtxPlane());
00580     }
00581     MSG("RecoBase",Msg::kDebug)
00582       << "    at plane " << shw->GetVtxPlane() << " dt2  " << du*du+dv*dv << "/" << tolTPos2 << " dt " << dt*1.e9 << "/" << tolTime*1e9 <<  "\n";
00583     
00584     if(dt<tolTime && 
00585        du*du+dv*dv<tolTPos2) return true;   // check whether vertex within minimum distance from trk 
00586   }
00587 
00588   // finally, check distance between vertices 
00589   Double_t trkZ=this->GetVtxZ();
00590   Double_t trkU=this->GetVtxU();
00591   Double_t trkV=this->GetVtxV();
00592   Double_t trkT=this->GetVtxT();
00593   Int_t trkPlane=this->GetVtxPlane();
00594   Double_t trkDirU=this->GetVtxDirCosU();
00595   Double_t trkDirV=this->GetVtxDirCosV();
00596   Double_t trkDirZ=this->GetVtxDirCosZ();
00597 
00598   if(fabs(shw->GetVtxZ()-this->GetEndZ())<fabs(shw->GetVtxZ()-this->GetVtxZ())){
00599     trkZ=this->GetEndZ();
00600     trkU=this->GetEndU();
00601     trkV=this->GetEndV();
00602     trkT=this->GetEndT();
00603     trkDirU=this->GetEndDirCosU();
00604     trkDirV=this->GetEndDirCosV();
00605     trkDirZ=this->GetEndDirCosZ();
00606     trkPlane=this->GetEndPlane();
00607   }
00608 
00609   Double_t dz=trkZ-shw->GetVtxZ();
00610   dt= trkT-shw->GetVtxT(); 
00611   //compensate dz for SM gap if necessary
00612   if (vldcptr->GetDetector()==Detector::kFar &&
00613       trkPlane<=pSMPlaneLast &&
00614       shw->GetVtxPlane()>=pSMPlaneFirst) dz+=zGapSM;
00615   else if (vldcptr->GetDetector()==Detector::kFar &&
00616            trkPlane>=pSMPlaneLast &&
00617            shw->GetVtxPlane()<=pSMPlaneFirst) dz-=zGapSM;
00618   
00619   Double_t du = (trkU - dz*trkDirU/trkDirZ) - shw->GetVtxU();
00620   Double_t dv = (trkV - dz*trkDirV/trkDirZ) - shw->GetVtxV();
00621   MSG("RecoBase",Msg::kDebug)
00622     << "    dvertex shower/track " << du
00623     << " " << dv << " " << dz <<" "  << dt*1.e9 << "\n";
00624   if (du*du+dv*dv<tolTPos2 &&
00625       fabs(dz)<tolZPos && 
00626       fabs(dt)<tolTime) return true;
00627 
00628   // now try shower End
00629   dz=trkZ-shw->GetEndZ();
00630   dt= trkT-shw->GetVtxT(); 
00631   //compensate dz for SM gap if necessary
00632   if (vldcptr->GetDetector()==Detector::kFar &&
00633       trkPlane<=pSMPlaneLast &&
00634       shw->GetVtxPlane()>=pSMPlaneFirst) dz+=zGapSM;
00635   else if (vldcptr->GetDetector()==Detector::kFar &&
00636            trkPlane>=pSMPlaneLast &&
00637            shw->GetEndPlane()<=pSMPlaneFirst) dz-=zGapSM;
00638   
00639   du = (trkU - dz*trkDirU/trkDirZ) - shw->GetEndU();
00640   dv = (trkV - dz*trkDirV/trkDirZ) - shw->GetEndV();
00641 
00642   MSG("RecoBase",Msg::kDebug)
00643     << "    dvertex shower/track " << du
00644     << " " << dv << " " << dz <<" "  << dt*1.e9 << "\n";
00645   if (du*du+dv*dv<tolTPos2 &&
00646       fabs(dz)<tolZPos && 
00647       fabs(dt)<tolTime) return true;
00648 
00649   return false;
00650 }
00651 
00652 //_____________________________________________________________________
00653 
00654 Bool_t CandTrackHandle::IsUnphysical(Float_t trkFrac,Float_t asymCut,
00655                                      Float_t xtalkFrac,Float_t xtalkCut)
00656 {
00657   // loop through track planes and calculate #gaps/#planes 
00658   // also consider balance of U/V hit planes
00659   // (for ND consider only area covered by partial plane)
00660 
00661   Float_t totTrkPlanes = 0;   // total number of hit planes in track
00662   Float_t nTrkPlanesU  = 0;   // total number of valid planes 
00663   Float_t nTrkPlanesV  = 0;   // between beg/end in U and V
00664   Float_t nHitPlanesU  = 0;   // total number of valid hit planes 
00665   Float_t nHitPlanesV  = 0;   // between beg/end in U and V
00666   Float_t nXTalkPlanes = 0;   // total number of xtalk-like track planes
00667 
00668   UgliGeomHandle ugh(*this->GetVldContext());
00669   for(int ipln=this->GetBegPlane();ipln<=this->GetEndPlane();ipln++){
00670     MSG("RecoBase",Msg::kDebug) << " plane = " << ipln << endl;
00671 
00672     //add up total number of track planes containing hits:
00673     if(this->IsTPosValid(ipln)) totTrkPlanes += 1;
00674 
00675     PlexPlaneId scintid(this->GetVldContext()->GetDetector(),ipln,kFALSE);
00676     UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(scintid);
00677     if(scintpln.IsValid()){
00678       MSG("RecoBase",Msg::kDebug) << "Valid plane (has scintillator)" << endl;
00679       //add up total number of valid track planes:
00680       if(this->GetVldContext()->GetDetector()==Detector::kFar) {        
00681         if(scintid.GetPlaneView()==PlaneView::kU) nTrkPlanesU += 1;
00682         else if(scintid.GetPlaneView()==PlaneView::kV) nTrkPlanesV += 1;
00683       }
00684       //for neardet, only consider the partial plane region:
00685       else if(this->GetVldContext()->GetDetector()==Detector::kNear){
00686         Float_t z = scintpln.GetZ0();
00687         Float_t u = this->GetVtxU();
00688         Float_t v = this->GetVtxV();
00689         if(this->IsTPosValid(ipln)) {
00690           u = this->GetU(ipln); 
00691           v = this->GetV(ipln);
00692         }
00693         else{
00694           Int_t count = 1;
00695           while(!this->IsTPosValid(ipln+count) && 
00696                 ipln+count<this->GetEndPlane()) count +=1;
00697           PlexPlaneId scintid_upp(this->GetVldContext()->GetDetector(),
00698                                   ipln+count,kFALSE);
00699           UgliScintPlnHandle scintpln_upp = ugh.GetScintPlnHandle(scintid_upp);
00700           Float_t upp_u = this->GetU(ipln+count);
00701           Float_t upp_v = this->GetV(ipln+count);
00702           Float_t upp_z = scintpln_upp.GetZ0();
00703           count = 1;
00704           while(!this->IsTPosValid(ipln-count) && 
00705                 ipln-count>this->GetBegPlane()) count +=1;
00706           PlexPlaneId scintid_low(this->GetVldContext()->GetDetector(),
00707                                   ipln-count,kFALSE);
00708           UgliScintPlnHandle scintpln_low = ugh.GetScintPlnHandle(scintid_low);
00709           Float_t low_u = this->GetU(ipln-count);
00710           Float_t low_v = this->GetV(ipln-count);
00711           Float_t low_z = scintpln_low.GetZ0();
00712           u = (z - low_z)*(upp_u - low_u)/(upp_z - low_z);
00713           v = (z - low_z)*(upp_v - low_v)/(upp_z - low_z);
00714         }
00715         
00716         //define a region for the ND to do this test in:
00717         const Float_t part_u_min = -0.2 * Munits::m;
00718         const Float_t part_v_min = -2.4 * Munits::m;
00719         const Float_t part_u_max = 2.4 * Munits::m;
00720         const Float_t part_v_max = 0.2 * Munits::m;
00721         if( u<part_u_min || u>part_u_max || //partial planes
00722             v<part_v_min || v>part_v_max ) continue;
00723         
00724         Float_t r2 = u*u+v*v;
00725         const Float_t coil_r     =  0.5 * Munits::m;
00726         const Float_t detect_r   =  2.8 * Munits::m;
00727         const Float_t coil_r2    = coil_r * coil_r;
00728         const Float_t detect_r2  = detect_r * detect_r; 
00729         if( r2<=coil_r2 || r2>=detect_r2 ) continue;  //coil hole + detector edge   
00730 
00731         if(scintid.GetPlaneView()==PlaneView::kU) nTrkPlanesU += 1;
00732         else if(scintid.GetPlaneView()==PlaneView::kV) nTrkPlanesV += 1;
00733         
00734       }
00735       if(this->IsTPosValid(ipln)) {
00736         //add upp number of valid planes that have hits:
00737         if(scintid.GetPlaneView()==PlaneView::kU) nHitPlanesU += 1;
00738         else if(scintid.GetPlaneView()==PlaneView::kV) nHitPlanesV += 1;
00739         //also count number of planes that have only xtalk-like hits
00740         if(this->GetPlaneCharge(ipln,CalStripType::kPE)<xtalkCut) 
00741           nXTalkPlanes+=1;
00742       }
00743     }
00744   }
00745 
00746   MSG("RecoBase",Msg::kDebug) 
00747     << endl << "Total Number of track planes:" 
00748     << ( this->GetEndPlane() - this->GetBegPlane() + 1) << endl
00749     << "Total Number of hit track planes:" << totTrkPlanes << endl
00750     << "Number of (valid) track planes U: " << nTrkPlanesU 
00751     << " V: " << nTrkPlanesV << endl
00752     << "Number of (valid) hit track planes U: " << nHitPlanesU 
00753     << " V: " << nHitPlanesV << endl
00754     << "Number of xtalk-like (valid) hit track planes:  " 
00755     << nXTalkPlanes << endl;
00756   
00757   //if track hits are mainly in spectrometer in ND do not perform tests
00758   if(this->GetVldContext()->GetDetector()==Detector::kNear){
00759     if(nHitPlanesU+nHitPlanesV<totTrkPlanes-nHitPlanesU-nHitPlanesV) 
00760       return false;
00761   }
00762 
00763   //check for nonsensical values:
00764   if(nTrkPlanesU+nTrkPlanesV<=0) return true; 
00765   if(nHitPlanesU+nHitPlanesV<=0) return true;
00766 
00767   //now check for gap planes, xtalk and asymmetry:
00768   if((nHitPlanesU+nHitPlanesV)/(nTrkPlanesU+nTrkPlanesV)<trkFrac) {
00769     MSG("RecoBase",Msg::kDebug) 
00770       << "IsUnphysical because trkFrac = "
00771       << (nHitPlanesU+nHitPlanesV)/(nTrkPlanesU+nTrkPlanesV)
00772       << " and cut-off is = " << trkFrac << endl;
00773     return true;
00774   }
00775   if(nXTalkPlanes/(nHitPlanesU+nHitPlanesV)>xtalkFrac) {
00776     MSG("RecoBase",Msg::kDebug) 
00777       << "IsUnphysical because xtalkFrac = "
00778       << nXTalkPlanes/(nHitPlanesU+nHitPlanesV)
00779       << " and cut-off is = " << xtalkFrac << endl;
00780     return true;
00781   }
00782   if(TMath::Abs(nHitPlanesU-nHitPlanesV)/
00783      TMath::Max(nHitPlanesU,nHitPlanesV)>asymCut) {
00784     MSG("RecoBase",Msg::kDebug) 
00785       << "IsUnphysical because asym = "
00786       << ( TMath::Abs(nHitPlanesU-nHitPlanesV) / 
00787            TMath::Max(nHitPlanesU,nHitPlanesV) )
00788       << " and cut-off is = " << asymCut << endl;
00789     return true;
00790   }
00791 
00792   return false;
00793 }
00794 
00795 void CandTrackHandle::SetNTrackStrip(Int_t n)
00796 {
00797   dynamic_cast<CandTrack *>(GetOwnedCandBase())->fNTrackStrip = n;
00798 }
00799 
00800 void CandTrackHandle::SetNTrackDigit(Int_t n)
00801 {
00802   dynamic_cast<CandTrack *>(GetOwnedCandBase())->fNTrackDigit = n;
00803 }
00804 
00805 void CandTrackHandle::SetNTimeFitDigit(Int_t n)
00806 {
00807   dynamic_cast<CandTrack *>(GetOwnedCandBase())->fNTimeFitDigit = n;
00808 }
00809 
00810 void CandTrackHandle::SetTimeFitChi2(Double_t x)
00811 {
00812   dynamic_cast<CandTrack *>(GetOwnedCandBase())->fTimeFitChi2 = x;
00813 }
00814 
00815 void CandTrackHandle::SetTimeForwardFitRMS(Double_t x)
00816 {
00817   dynamic_cast<CandTrack *>(GetOwnedCandBase())->fTimeForwardFitRMS = x;
00818 }
00819 void CandTrackHandle::SetTimeForwardFitNDOF(Int_t x)
00820 {
00821   dynamic_cast<CandTrack *>(GetOwnedCandBase())->fTimeForwardFitNDOF = x;
00822 }
00823 
00824 void CandTrackHandle::SetTimeBackwardFitRMS(Double_t x)
00825 {
00826   dynamic_cast<CandTrack *>(GetOwnedCandBase())->fTimeBackwardFitRMS = x;
00827 }
00828 
00829 void CandTrackHandle::SetTimeBackwardFitNDOF(Int_t x)
00830 {
00831   dynamic_cast<CandTrack *>(GetOwnedCandBase())->fTimeBackwardFitNDOF = x;
00832 }
00833 
00834 
00835 Int_t CandTrackHandle::GetNTrackStrip() const
00836 {
00837   return dynamic_cast<const CandTrack *>(GetCandBase())->fNTrackStrip;
00838 }
00839 
00840 Int_t CandTrackHandle::GetNTrackDigit() const
00841 {
00842   return dynamic_cast<const CandTrack *>(GetCandBase())->fNTrackDigit;
00843 }
00844 
00845 Int_t CandTrackHandle::GetNTimeFitDigit() const
00846 {
00847   return dynamic_cast<const CandTrack *>(GetCandBase())->fNTimeFitDigit;
00848 }
00849 
00850 Double_t CandTrackHandle::GetTimeFitChi2() const
00851 {
00852   return dynamic_cast<const CandTrack *>(GetCandBase())->fTimeFitChi2;
00853 }
00854 
00855 
00856 Double_t CandTrackHandle::GetTimeForwardFitRMS() const
00857 {
00858   return dynamic_cast<const CandTrack *>(GetCandBase())->fTimeForwardFitRMS;
00859 }
00860 
00861 
00862 Int_t CandTrackHandle::GetTimeForwardFitNDOF() const
00863 {
00864   return dynamic_cast<const CandTrack *>(GetCandBase())->fTimeForwardFitNDOF;
00865 }
00866 
00867 
00868 Double_t CandTrackHandle::GetTimeBackwardFitRMS() const
00869 {
00870   return dynamic_cast<const CandTrack *>(GetCandBase())->fTimeBackwardFitRMS;
00871 }
00872 
00873 
00874 Int_t CandTrackHandle::GetTimeBackwardFitNDOF() const
00875 {
00876   return dynamic_cast<const CandTrack *>(GetCandBase())->fTimeBackwardFitNDOF;
00877 }
00878 
00879 bool CandTrackHandle::IsContained() 
00880 {
00881 
00882   bool contained = (GetVtxTrace()>0.1 && GetEndTrace()>0.1 &&
00883                     GetVtxnActiveUpstream()>2 && GetEndnActiveDownstream()>2 &&
00884                     GetVtxDistToEdge()>0.1 && GetEndDistToEdge()>0.1 );
00885   return contained;
00886 }
00887 
00888 XXXITRIMP(CandTrackHandle)