loon: AlgTrack.cxx Source File

00001 
00002 // $$
00003 //
00004 // AlgTrack.cxx
00005 //
00006 // Begin_Html<img src="../../pedestrians.gif" align=center>
00007 // <a href="../source_warning.html">Warning for beginners</a>.<br> 
00008 //
00009 // This is an abstract Algorithm reconstruction object.
00010 // This class should and cannot be instantiated.
00011 //
00012 //
00013 // Author:  R. Lee
00014 //
00015 // Also see <a href="../../root_crib/index.html">The ROOT Crib</a> and 
00016 // <a href="../AlgTrack.html"> AlgTrack Classes</a> (part of
00017 // <a href="../index.html">The MINOS Class User Guide</a>)End_Html
00019 
00020 #include "Algorithm/AlgHandle.h"
00021 #include "Conventions/Munits.h"
00022 #include "Conventions/MinosMaterial.h"
00023 #include "Conventions/PlaneView.h"
00024 #include "Conventions/SimFlag.h"
00025 #include "MessageService/MsgService.h"
00026 #include "Navigation/NavKey.h"
00027 #include "Navigation/NavSet.h"
00028 #include "Navigation/XxxItr.h"
00029 #include "RecoBase/AlgTrack.h"
00030 #include "RecoBase/ArrivalTime.h"
00031 #include "RecoBase/CandRecoHandle.h"
00032 #include "RecoBase/CandStripHandle.h"
00033 #include "RecoBase/CandTrackHandle.h"
00034 #include "RecoBase/PropagationVelocity.h"
00035 #include "MessageService/MsgService.h"
00036 #include "Plex/PlexPlaneId.h"
00037 #include "UgliGeometry/UgliGeomHandle.h"
00038 #include "UgliGeometry/UgliScintPlnHandle.h"
00039 #include "UgliGeometry/UgliSteelPlnHandle.h"
00040 #include "DataUtil/PlaneOutline.h"
00041 #include "TMath.h"
00042 #include <cmath>
00043 
00044 ClassImp(AlgTrack)
00045 
00046 //______________________________________________________________________
00047 CVSID("$Id: AlgTrack.cxx,v 1.49 2007/11/11 05:26:09 rhatcher Exp $");
00048 
00049 
00050 void AlgTrack::SetUVZT(CandTrackHandle *candtrack)
00051 {
00052   SetUVZ(candtrack);
00053   SetT(candtrack);
00054   SetdS(candtrack);
00055 }
00056 
00057 void AlgTrack::SetUVZ(CandTrackHandle *candtrack)
00058 {
00059 
00060   const CandStripHandle *firststrip = 
00061     dynamic_cast<const CandStripHandle *>(candtrack->GetDaughter(0));
00062   if (!firststrip) return;  // if no strips, do nothing
00063 
00064   const VldContext *vldc = firststrip->GetVldContext();
00065   Detector::Detector_t det = vldc->GetDetector();
00066 
00067   UgliGeomHandle ugh = UgliGeomHandle(*vldc);
00068 
00069   CandStripHandleItr orderedstripItr(candtrack->GetDaughterIterator());
00070   CandStripHandleKeyFunc *stripKf = orderedstripItr.CreateKeyFunc();
00071   stripKf->SetFun(CandStripHandle::KeyFromPlane);
00072   orderedstripItr.GetSet()->AdoptSortKeyFunc(stripKf);
00073   stripKf = 0;
00074 
00075   Bool_t first(1);
00076   Int_t oldplane = -1;
00077   PlaneView::PlaneView_t oldview = PlaneView::kUnknown;
00078   Float_t tpos = 0.;
00079   Float_t ph = 0.;
00080   Int_t idir=1;
00081   Int_t vtxplane = candtrack->GetBegPlane();
00082   Int_t endplane = candtrack->GetEndPlane();
00083   if (vtxplane>endplane) idir=-1;
00084   Int_t begplane = min(vtxplane,endplane);
00085   endplane = max(vtxplane,endplane);
00086   begplane = -1;
00087   endplane = -1;
00088 
00089 // find begin and end planes
00090   while (CandStripHandle *strip = dynamic_cast<CandStripHandle*>(orderedstripItr())) {
00091     if (strip->GetPlaneView()==PlaneView::kU ||
00092         strip->GetPlaneView()==PlaneView::kV) {
00093       if (begplane<0 || strip->GetPlane()<begplane) {
00094         begplane = strip->GetPlane();
00095       }
00096       if (endplane<0 || strip->GetPlane()>endplane) {
00097         endplane = strip->GetPlane();
00098       }
00099     }
00100   }
00101   Int_t *planelist = new Int_t[endplane+1];
00102   Int_t *stripendlist = new Int_t[endplane+1];
00103   for (int i=0; i<=endplane; i++) {
00104     planelist[i] = 0;
00105     stripendlist[i] = 0;
00106   }
00107 
00108 // transverse positions
00109 
00110   orderedstripItr.Reset();
00111 
00112   while (CandStripHandle *strip = dynamic_cast<CandStripHandle*>(orderedstripItr())) {
00113     if (strip->GetPlaneView()==PlaneView::kU ||
00114         strip->GetPlaneView()==PlaneView::kV) {
00115       planelist[strip->GetPlane()] = strip->GetPlaneView();
00116       if (strip->GetNDigit(StripEnd::kNegative)) {
00117         stripendlist[strip->GetPlane()] += (1<<StripEnd::kNegative);
00118       }
00119       if (strip->GetNDigit(StripEnd::kPositive)) {
00120         stripendlist[strip->GetPlane()] += (1<<StripEnd::kPositive);
00121       }
00122       if (first || strip->GetPlane()==oldplane) {
00123         if (first) {
00124           first = 0;
00125           oldplane = strip->GetPlane();
00126           oldview = strip->GetPlaneView();
00127         }
00128         Float_t stripph = strip->GetCharge();
00129         tpos += strip->GetTPos()*stripph;
00130         if(stripph>0){
00131           ph += stripph;
00132         }
00133       }
00134       else {
00135         if(ph>0){
00136           tpos /= ph;
00137         }
00138             
00139         if (oldview==PlaneView::kU) {
00140           candtrack->SetU(oldplane,tpos);
00141         } else if (oldview==PlaneView::kV) {
00142           candtrack->SetV(oldplane,tpos);
00143         }
00144         oldplane = strip->GetPlane();
00145         oldview = strip->GetPlaneView();
00146         Float_t stripph = strip->GetCharge();
00147         tpos = strip->GetTPos()*stripph;
00148         if(ph>0){
00149           ph = stripph;
00150         }
00151       }
00152     }
00153   }
00154   if (ph>0.) {
00155     tpos /= ph;
00156     if (oldview==PlaneView::kU) {
00157       candtrack->SetU(oldplane,tpos);
00158     } else if (oldview==PlaneView::kV) {
00159       candtrack->SetV(oldplane,tpos);
00160     }
00161   }
00162 
00163   // create a starting and stopping PlexPlaneId
00164   Int_t idirAdjoin = (begplane>endplane) ? -1 : 1; // direction
00165   PlexPlaneId planeid(det,begplane,false);
00166   PlexPlaneId planeidend(det,endplane,false);
00167   PlexPlaneId planeid_toofar = planeidend.GetAdjoinScint(idirAdjoin);
00168   
00169   // after each step, if not reached the end, move to adjoining *scint* plane
00170   for ( ; planeid != planeid_toofar ; 
00171           planeid  = planeid.GetAdjoinScint(idirAdjoin) ) {
00172 
00173     if ( ! planeid.IsValid() ) break; // quit if hit something meaningless
00174     UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(planeid);
00175 
00176     if ( ! scintpln.IsValid() ) continue;
00177 
00178     Int_t iplane = planeid.GetPlane();
00179     Float_t zpos = scintpln.GetZ0();
00180     for (Int_t planeview=PlaneView::kU; planeview<=PlaneView::kV; planeview++) {
00181       if (planelist[iplane]!=planeview) {
00182         Int_t found = 0;
00183         Int_t nfound = 0;
00184         Int_t jplane[2];
00185         // search upstream
00186         for (Int_t iplane1=iplane-1; iplane1>=begplane && !found; iplane1--){
00187           if (planelist[iplane1]==planeview) {
00188             found++;
00189             jplane[nfound] = iplane1;
00190             nfound++;
00191           }
00192         }
00193         Int_t ndown = 1;
00194         if (!found) {
00195           ndown++;
00196         }
00197         found = 0;
00198         // search downstream
00199         for (Int_t iplane1=iplane+1; iplane1<=endplane && found<ndown; iplane1++) {
00200           if (planelist[iplane1]==planeview) {
00201             found++;
00202             jplane[nfound] = iplane1;
00203             nfound++;
00204           }
00205         }
00206         if (!found) {
00207           for (Int_t iplane1=jplane[0]-1; iplane1>=begplane && !found; iplane1--) {
00208             if (planelist[iplane1]==planeview) {
00209               found++;
00210               jplane[nfound] = iplane1;
00211               nfound++;
00212             }
00213           }
00214         }
00215         if (nfound==1) {
00216           if (planeview==PlaneView::kU) {
00217             candtrack->SetU(iplane,candtrack->GetU(jplane[0]));
00218           } else {
00219             candtrack->SetV(iplane,candtrack->GetV(jplane[0]));
00220           }
00221         }
00222         else {
00223           Float_t z[2];
00224           z[0] = candtrack->GetZ(jplane[0]);
00225           z[1] = candtrack->GetZ(jplane[1]);
00226           Float_t x[2];
00227           if (planeview==PlaneView::kU) {
00228             x[0] = candtrack->GetU(jplane[0]);
00229             x[1] = candtrack->GetU(jplane[1]);
00230             candtrack->SetU(iplane,x[0]+(x[1]-x[0])/(z[1]-z[0])*(zpos-z[0]));
00231           } else {
00232             x[0] = candtrack->GetV(jplane[0]);
00233             x[1] = candtrack->GetV(jplane[1]);
00234             candtrack->SetV(iplane,x[0]+(x[1]-x[0])/(z[1]-z[0])*(zpos-z[0]));
00235           }
00236         }
00237       }
00238     }
00239  
00240   }
00241 
00242   delete [] planelist;
00243   delete [] stripendlist;
00244 
00245 }
00246 
00247 void AlgTrack::SetT(CandTrackHandle *candtrack)
00248 {
00249 
00250 // we take a weighted average of hits in the same
00251 // plane.  the proper way to do this would be to consider hits coming
00252 // from the same PMT and use the earliest time.
00253 
00254   MSG("RecoBase",Msg::kDebug) << "starting SetT" << endl;
00255 
00256   TIter stripItr(candtrack->GetDaughterIterator());
00257   CandStripHandle *firststrip = 0;
00258   firststrip = dynamic_cast<CandStripHandle*>(stripItr());
00259   UgliGeomHandle *ugh = 0;
00260   SimFlag::SimFlag_t simFlag = SimFlag::kUnknown;
00261   if (firststrip) {
00262     ugh = new UgliGeomHandle(*firststrip->GetVldContext());
00263     simFlag = firststrip->GetVldContext()->GetSimFlag();
00264   }
00265 
00266   Double_t timewgt[2] = {0.,0.};
00267   Double_t time[2] = {0.,0.};
00268   StripEnd::StripEnd_t stripend[2] = {StripEnd::kNegative,StripEnd::kPositive};
00269 
00270   Double_t apos;
00271 
00272   ArrivalTime arrtime;
00273 
00274   CalTimeType::CalTimeType_t caltimetype = CalTimeType::kNone;
00275 
00276 // these constants should be taken from the database
00277 
00278   Double_t propagation_velocity = PropagationVelocity::Velocity(simFlag)*Munits::ns;
00279   //  Double_t c1 = -14.45;
00280   // Double_t c2 =  7.498;
00281   // Double_t c3 = -1.566;
00282 
00283   CandStripHandleItr orderedstripItr(candtrack->GetDaughterIterator());
00284   CandStripHandleKeyFunc *stripKf = orderedstripItr.CreateKeyFunc();
00285   stripKf->SetFun(CandStripHandle::KeyFromPlane);
00286   orderedstripItr.GetSet()->AdoptSortKeyFunc(stripKf);
00287   stripKf = 0;
00288 
00289   Int_t oldplane = -1;
00290   PlaneView::PlaneView_t oldview = PlaneView::kUnknown;
00291 
00292   Bool_t first(1);
00293   while (CandStripHandle *strip = dynamic_cast<CandStripHandle*>(orderedstripItr())) {
00294     if (strip->GetPlaneView()==PlaneView::kU ||
00295         strip->GetPlaneView()==PlaneView::kV) {
00296       if (strip->GetPlaneView()==PlaneView::kU) {
00297         apos = candtrack->GetV(strip->GetPlane());
00298       } else {
00299         apos = candtrack->GetU(strip->GetPlane());
00300       }
00301       apos=min(apos,4.0); 
00302       apos=max(apos,-4.0);  // sanity check on apos bounds
00303       UgliStripHandle striphandle = ugh->GetStripHandle(strip->GetStripEndId());
00304       caltimetype = strip->GetCalTimeType();
00305       if (first || strip->GetPlane()==oldplane) {
00306         if (first) {
00307           first = 0;
00308           oldplane = strip->GetPlane();
00309           oldview = strip->GetPlaneView();
00310         }
00311         Double_t thisph[2] = {
00312           strip->GetCharge(CalDigitType::kPE,stripend[0]),
00313           strip->GetCharge(CalDigitType::kPE,stripend[1])};
00314         Double_t thiswgt[2] = {
00315           min(0.4,arrtime.Weight(max(1.,thisph[0]))),
00316           min(0.4,arrtime.Weight(max(1.,thisph[1])))};
00317         for (int i=0; i<2; i++) {
00318           if (thisph[i]>0.) {
00319             Double_t plusminus = 1.;
00320             if ((strip->GetPlaneView()==PlaneView::kV && i==1) ||
00321                 (strip->GetPlaneView()==PlaneView::kU && i==0)) {
00322               plusminus = -1.;
00323             }
00324             Double_t thistime = strip->GetBegTime(stripend[i])-1.e-9*(striphandle.ClearFiber(stripend[i])+striphandle.WlsPigtail(stripend[i])+(striphandle.GetHalfLength()+plusminus*apos))/propagation_velocity;
00325             //        Double_t logph = log(thisph[i])/2.3;
00326             // thistime -= (c1*logph+c2*logph*logph+c3*logph*logph*logph)*1.e-9;
00327             MSG("RecoBase",Msg::kDebug) << "plane " << strip->GetPlane() << " end " << i << " raw time " << strip->GetBegTime(stripend[i])*1e9 << " pulse height " << thisph[i] << " weight " << thiswgt[i] << " clearfiber " << striphandle.ClearFiber(stripend[i]) << " wlspigtail " << striphandle.WlsPigtail(stripend[i]) << " halflength " << striphandle.GetHalfLength() << " apos " << apos << " 3D time " << thistime*1e9 << endl;
00328             time[i] += thiswgt[i]*thistime;
00329             timewgt[i] += thiswgt[i];
00330             MSG("RecoBase",Msg::kDebug) << "summed 3D time " << time[i]/timewgt[i] << endl;
00331           }
00332         }
00333       }
00334       else {
00335         for (int i=0; i<2; i++) {
00336           if (timewgt[i]>0.) {
00337             time[i] /= timewgt[i];
00338             candtrack->SetT(oldplane,stripend[i],time[i]);
00339             MSG("RecoBase",Msg::kDebug) << "setting time " << oldplane << " " << stripend[i] << " " << time[i]*1e9 << endl;
00340           }
00341         }
00342         oldplane = strip->GetPlane();
00343         oldview = strip->GetPlaneView();
00344         Double_t thisph[2] = {
00345           strip->GetCharge(CalDigitType::kPE,stripend[0]),
00346           strip->GetCharge(CalDigitType::kPE,stripend[1])};
00347         Double_t thiswgt[2] = {
00348           min(0.4,arrtime.Weight(max(1.,thisph[0]))),
00349           min(0.4,arrtime.Weight(max(1.,thisph[1])))};
00350         for (int i=0; i<2; i++) {
00351           time[i] = 0.;
00352           timewgt[i] = 0.;
00353           if (thisph[i]>0.) {
00354             Double_t plusminus = 1.;
00355             if ((strip->GetPlaneView()==PlaneView::kV && i==1) ||
00356                 (strip->GetPlaneView()==PlaneView::kU && i==0)) {
00357               plusminus = -1.;
00358             }
00359             Double_t thistime = strip->GetBegTime(stripend[i])-1.e-9*(striphandle.ClearFiber(stripend[i])+striphandle.WlsPigtail(stripend[i])+(striphandle.GetHalfLength()+plusminus*apos))/propagation_velocity;
00360             //            Double_t logph = log(thisph[i])/2.3;
00361             //  thistime -= (c1*logph+c2*logph*logph+c3*logph*logph*logph)*1.e-9;
00362             MSG("RecoBase",Msg::kDebug) << "plane " << strip->GetPlane() << " end " << i << " raw time " << strip->GetBegTime(stripend[i])*1e9 << " pulse height " << thisph[i] << " weight " << thiswgt[i] << " clearfiber " << striphandle.ClearFiber(stripend[i]) << " wlspigtail " << striphandle.WlsPigtail(stripend[i]) << " halflength " << striphandle.GetHalfLength() << " apos " << apos << " 3D time " << thistime*1e9 << endl;
00363             time[i] = thiswgt[i]*thistime;
00364             timewgt[i] = thiswgt[i];
00365             MSG("RecoBase",Msg::kDebug) << "summed 3D time " << time[i]/timewgt[i] << endl;
00366           }
00367         }
00368       }
00369     }
00370   }
00371   for (int i=0; i<2; i++) {
00372     if (timewgt[i]>0.) {
00373       time[i] /= timewgt[i];
00374       candtrack->SetT(oldplane,stripend[i],time[i]);
00375       MSG("RecoBaseAlgTrack",Msg::kDebug) << "setting time " << oldplane << " " << stripend[i] << " " << time[i]*1e9 << endl;
00376     }
00377   }
00378 
00379   if (ugh) delete ugh;
00380 
00381 }
00382 
00383 
00384 void AlgTrack::SetdS(CandTrackHandle *candtrack)
00385 {
00386   // Calculate path lengths from the vertex plane to each of the planes
00387   // from the vertex to the end.  The "dS" value is the sum of the 
00388   // line segments in meters, the "Range" is the distance in g/cm**2.
00389 
00390   TIter stripItr(candtrack->GetDaughterIterator());
00391   CandStripHandle *firststrip = 0;
00392   firststrip = dynamic_cast<CandStripHandle*>(stripItr());
00393   if (!firststrip) {
00394     MSG("RecoBaseAlgTrack",Msg::kDebug) << "no strips, returning" << endl;
00395   }
00396   if (!firststrip) return; // if no strips do nothing
00397   //RWH: is this really necesary? can't we get the VldContext
00398   //     from the CandTrackHandle itself?
00399 
00400   const VldContext *vldc = firststrip->GetVldContext();
00401   Detector::Detector_t det = vldc->GetDetector();
00402   UgliGeomHandle ugh = UgliGeomHandle(*vldc);
00403 
00404   Int_t vtxplane = candtrack->GetVtxPlane();
00405 
00406   Int_t endplane = candtrack->GetEndPlane();
00407   Int_t idir = (vtxplane>endplane) ? 1 : -1; 
00408 
00409   Bool_t  first = true;
00410   Float_t ds = 0.;
00411   Float_t range = 0.;
00412 
00413   Float_t uvzold[3] = { 0, 0, 0}, uvznew[3] = { 0, 0, 0};
00414 
00415   // correct range by direction cosine (start w/ initial value)
00416   // futher segments will adjust this
00417   Float_t dzds = fabs(candtrack->GetDirCosZ());
00418   if (dzds==0.) dzds=1.;
00419 
00420  MinosMaterial::StdMaterial_t fStdMatSc=MinosMaterial::ePolystyreneMinos;
00421  MinosMaterial::StdMaterial_t fStdMatSteel=MinosMaterial::eIronMinos;
00422 
00423  const double scint_density =  MinosMaterial::Density(fStdMatSc);// used to be 1.032  g/cm**2
00424  const double steel_density = MinosMaterial::Density(fStdMatSteel); // used to be  7.87 g/cm**2
00425 
00426  MSG("RecoBaseAlgTrack",Msg::kDebug) << " densitities " << scint_density << " " << steel_density << endl;
00427 
00428  
00429   // create a starting and stopping PlexPlaneId
00430   PlexPlaneId plnid(det,endplane,false);    // start w/ vtx scint plane
00431   PlexPlaneId plnidend(det,vtxplane,false); // last scint plane
00432 
00433   MSG("RecoBaseAlgTrack",Msg::kDebug) << "starting with plane " << endplane 
00434                                       << " and ending with plane " << vtxplane << endl;
00435 
00436   // build an id one too far for use in the iteration
00437   PlexPlaneId plnid_toofar = plnidend.GetAdjoin(idir);
00438 
00439   // advance one more at vertex and at end for steel plane
00440   PlexPlaneId plnid0 = plnid.GetAdjoin(-idir);
00441   if (plnid0.IsValid()) {
00442     plnid = plnid.GetAdjoin(-idir);
00443   }
00444   PlexPlaneId plnid_toofar0 = plnid_toofar.GetAdjoin(idir);
00445   if (plnid_toofar0.IsValid()) {
00446     plnid_toofar = plnid_toofar.GetAdjoin(idir);
00447   }
00448   Bool_t isfirst = 1;
00449 
00450   // after each step, if not reached the end, move to next plane (steel/scint)
00451   for ( ; plnid != plnid_toofar ; plnid  = plnid.GetAdjoin(idir) ) {
00452     Int_t iplane = plnid.GetPlane();
00453     if ( ! plnid.IsValid() ){ 
00454       MSG("RecoBaseAlgTrack",Msg::kDebug) << "plane " << plnid.GetPlane() 
00455                                           << " not valid, ending loop" << endl;
00456       break;
00457     }
00458 
00459     MSG("RecoBaseAlgTrack",Msg::kDebug) 
00460       << "plane " << iplane*idir << "/" << plnid_toofar.GetPlane()*idir << " dzds = " << dzds 
00461       << " uvzold = " << uvzold[0] << " " << uvzold[1] << " " << uvzold[2] 
00462       << " ds = " << ds << " range = " << range << endl;
00463 
00464     if ( ! plnid.IsSteel() ) {
00465       UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(plnid);
00466       if ( ! scintpln.IsValid() ) continue;
00467 
00468       Float_t zpos = scintpln.GetZ0();
00469       if (candtrack->IsTPosValid(iplane) ) {
00470         if (first) {
00471           first = 0;
00472           uvznew[0] = candtrack->GetU(iplane);
00473           uvznew[1] = candtrack->GetV(iplane);
00474           uvznew[2] = zpos;
00475         }
00476         else {
00477           uvzold[0] = uvznew[0]; uvznew[0] = candtrack->GetU(iplane);
00478           uvzold[1] = uvznew[1]; uvznew[1] = candtrack->GetV(iplane);
00479           uvzold[2] = uvznew[2]; uvznew[2] = zpos;
00480 
00481           if( TMath::IsNaN(uvzold[0]) || TMath::IsNaN(uvzold[1]) 
00482               || TMath::IsNaN(uvzold[2]) )
00483             MSG("AlgTrack", Msg::kWarning) << "IsNaN found during calculation of prev position "
00484                                            << "u = " << uvzold[0]
00485                                            << "v = " << uvzold[1]
00486                                            << "z = " << uvzold[2] << endl;
00487           if( TMath::IsNaN(uvznew[0]) || TMath::IsNaN(uvznew[1]) 
00488               || TMath::IsNaN(uvznew[2]) )
00489             MSG("AlgTrack", Msg::kWarning) << "IsNaN found during calculation of new position "
00490                                            << "u = " << uvznew[0]
00491                                            << "v = " << uvznew[1]
00492                                            << "z = " << uvznew[2] << endl;
00493 
00494           Float_t duvz[3] = {uvznew[0]-uvzold[0],
00495                              uvznew[1]-uvzold[1],
00496                              uvznew[2]-uvzold[2]};
00497           Float_t adds=1;
00498           Double_t insq=duvz[0]*duvz[0]+duvz[1]*duvz[1]+duvz[2]*duvz[2];
00499           if(!(TMath::IsNaN(insq)) && insq>0) adds=TMath::Sqrt(insq);
00500           ds += adds;
00501           dzds = fabs(duvz[2]) / adds;
00502           if(dzds<0.06)dzds=1;  // if dzds less that 1 plane/m, not physical
00503         }
00504         MSG("RecoBaseAlgTrack",Msg::kDebug) << "Setting dS " << iplane << " " << ds << endl;
00505         candtrack->SetdS(iplane,ds);
00506       } else {
00507         MSG("RecoBaseAlgTrack",Msg::kDebug) << "is not valid tpos plane" << endl;
00508       }
00509 
00510       // The dE/dx in scintillator (in units of GeV cm**2/g) is larger
00511       // than it would be in Iron.  So, the range in scintillator
00512       // should really be scaled up to reflect this.
00513       // This means that the range returned by this routine is
00514       // an equivalent range in Fe.
00515       // Unfortunately the dE/dx ratio for scint to iron is not constant,
00516       // and depnds on the energy (or range so far) of the muon.
00517       // The dE/dx ratio is well described by two polynomials of Log(range),
00518       // one for muons with a range below 4 g/cm**2, and one for muons 
00519       // above 4.  See NuMI-Note-Comp-1102 for further discussion.
00520 
00521 
00522       double dedx_ratio = 1; // ratio of scintillator dedx to iron
00523                              // dedx for a given muon energy
00524       
00525       if (range > 0.0) {
00526         
00527           double logrange = TMath::Log(range);
00528           if (range<=4) {
00529               dedx_ratio
00530                   = 1.419-0.0223*logrange
00531                   + 0.0055*logrange*logrange;
00532           }
00533           else {
00534               dedx_ratio
00535                   = 1.382+0.0124*logrange
00536                   - 0.0052*logrange*logrange
00537                   + 0.00019*logrange*logrange*logrange;
00538           }
00539        
00540       }
00541       
00542       // range is in g/cm**2
00543        
00544       range += (dedx_ratio * scint_density * 
00545                 2.0*scintpln.GetHalfThickness()/Munits::cm)/dzds;
00546  
00547      double drange = dedx_ratio * scint_density * 
00548         2.0*scintpln.GetHalfThickness()/Munits::cm/dzds;
00549    
00550       MSG("RecoBaseAlgTrack",Msg::kDebug)  << " scint plane " << iplane << " drange " << drange  << " dzds = " << dzds << " total range = " << range << endl;
00551 
00552    
00553       candtrack->SetRange(iplane,range);
00554       
00555     }
00556     else {
00557 
00558     UgliSteelPlnHandle steelpln = ugh.GetSteelPlnHandle(plnid);
00559       if ( ! steelpln.IsValid() ) continue;
00560       // range is in g/cm**2
00561 
00562       
00563  // if near detector and endplane is in spectrometer, add two planes to track
00564       float nsteelhalf = 2;
00565       if( isfirst){
00566         if ( det == Detector::kNear && 
00567              (unsigned int)endplane > PlexPlaneId::LastPlaneNearCalor() )nsteelhalf=5;
00568         else if(det == Detector::kFar)nsteelhalf = 1;
00569         
00570       } 
00571       range += (steel_density * 
00572                 nsteelhalf*steelpln.GetHalfThickness()/Munits::cm)/dzds;
00573 
00574       double drange = steel_density * 
00575         nsteelhalf*steelpln.GetHalfThickness()/Munits::cm/dzds;
00576       MSG("RecoBaseAlgTrack",Msg::kDebug) << " steel plane " << iplane << " drange " << drange  << " dzds = " << dzds << " total range = " << range << endl;
00577       
00578       
00579       candtrack->SetRange(iplane,range);
00580     }
00581     isfirst = 0;
00582   }
00583   // subtract 1/2 steel plane at vertex
00584   range -= (steel_density * 0.0127/Munits::cm)/dzds;
00585         
00586   candtrack->SetRange(vtxplane,range);
00587 
00588 
00589   MSG("RecoBaseAlgTrack",Msg::kDebug) << "Total range = " << range << " gm/cm**3" << endl;
00590 
00591 }
00592 
00594 void AlgTrack::CalculateTrace(CandTrackHandle &cth){
00595   
00596   double m[2]; double c[2]; 
00597   TIter stripItr(cth.GetDaughterIterator());
00598   CandStripHandle *firststrip = 0;
00599   firststrip = dynamic_cast<CandStripHandle*>(stripItr());
00600   if (!firststrip) {
00601     MSG("RecoBaseAlgTrack",Msg::kDebug) << "no strips, returning" << endl;
00602   }
00603   if (!firststrip) return; // if no strips do nothing
00604   //RWH: is this really necesary? can't we get the VldContext
00605   //     from the CandTrackHandle itself?
00606 
00607   const VldContext *vldc = firststrip->GetVldContext();
00608  
00609   UgliGeomHandle ugh(*vldc); 
00610   PlaneOutline pl;
00611   float detzmin=0;
00612   float detzmax=999;
00613   ugh.GetZExtent(detzmin,detzmax,-1);
00614 
00615   int dir=1;
00616   if(cth.GetVtxZ()>cth.GetEndZ())dir=-1;
00617 
00618   // For vertex - loop over planes upstream from vertex until exit detector
00619   // trace = linear extrapolation to track entry point
00620   bool IsOutside = false; bool IsInside = true;
00621   double trace=0;
00622   double tracez=0;
00623   double z = cth.GetVtxZ();
00624   double u=0; double v=0; double x=0; double y=0; double r=0;
00625   float xedge,yedge,dist;
00626   int nActiveUpstream = 0;
00627   Detector::Detector_t detector = vldc->GetDetector();
00628   if(cth.GetVtxDirCosZ()!=0.) {
00629     m[0]=cth.GetVtxDirCosU()/cth.GetVtxDirCosZ();
00630     m[1]=cth.GetVtxDirCosV()/cth.GetVtxDirCosZ();
00631     c[0]=cth.GetVtxU()-(m[0]*cth.GetVtxZ());
00632     c[1]=cth.GetVtxV()-(m[1]*cth.GetVtxZ());
00633     PlexPlaneId plnid(detector,cth.GetVtxPlane(),false);
00634     while(!IsOutside){
00635       plnid = plnid.GetAdjoinScint(-dir);
00636 
00637       if(plnid.IsValid()){
00638         UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(plnid);
00639         z = scintpln.GetZ0();
00640         u = c[0] + m[0]*z;
00641         v = c[1] + m[1]*z;
00642         x = 0.707*(u-v);
00643         y = 0.707*(u+v);
00644         r = pow(x*x+y*y,0.5);
00645         pl.DistanceToEdge(x, y,
00646                           plnid.GetPlaneView(),
00647                           plnid.GetPlaneCoverage(),
00648                           dist, xedge, yedge);
00649         
00650        IsInside = pl.IsInside( x, y,
00651                                plnid.GetPlaneView(),
00652                                plnid.GetPlaneCoverage(),true);
00653        IsInside |= r<0.5;
00654        if(IsInside){
00655          nActiveUpstream++;
00656          trace = pow(pow(u-cth.GetVtxU(),2)+pow(v-cth.GetVtxV(),2)+pow(z-cth.GetVtxZ(),2),0.5);
00657          tracez = fabs(z - cth.GetVtxZ());
00658        }
00659       }
00660       IsOutside = !IsInside | !plnid.IsValid();
00661     }
00662     if(plnid.IsValid()) trace += dist;
00663   }
00664   cth.SetVtxTrace(trace);
00665   cth.SetVtxTraceZ(tracez);
00666   cth.SetVtxnActiveUpstream(nActiveUpstream);
00667 
00668   // set distance from track vertext to nearest edge on that plane
00669   z = cth.GetVtxZ();
00670   u = cth.GetVtxU();
00671   v = cth.GetVtxV();
00672   x = 0.707*(u-v);
00673   y = 0.707*(u+v);
00674   r = pow(x*x+y*y,0.5);
00675   PlexPlaneId plnid(detector,cth.GetVtxPlane(),false);
00676   dist = 0.;
00677   IsInside=true;
00678   if(plnid.IsValid()){
00679     pl.DistanceToEdge(x, y,
00680                       plnid.GetPlaneView(),
00681                       plnid.GetPlaneCoverage(),
00682                       dist, xedge, yedge);
00683     IsInside = pl.IsInside( x, y,
00684                             plnid.GetPlaneView(),
00685                             plnid.GetPlaneCoverage(),true);
00686     
00687   }
00688   if(!IsInside && r<0.5) dist=0;
00689   cth.SetVtxDistToEdge(dist);
00690 
00691   // For vertex - loop over planes upstream from vertex until exit detector
00692   // trace = linear extrapolation to track entry point
00693   IsOutside = false; IsInside=true;
00694   trace=0; tracez=0;
00695   z = cth.GetEndZ();
00696   u=0; v=0; x=0; y=0; r=0;
00697   int nActiveDownstream = 0;
00698   if(cth.GetEndDirCosZ()!=0.) {
00699     m[0]=cth.GetEndDirCosU()/cth.GetEndDirCosZ();
00700     m[1]=cth.GetEndDirCosV()/cth.GetEndDirCosZ();
00701     c[0]=cth.GetEndU()-(m[0]*cth.GetEndZ());
00702     c[1]=cth.GetEndV()-(m[1]*cth.GetEndZ());
00703    
00704     PlexPlaneId plnid(detector,cth.GetEndPlane(),false);
00705 
00706     while(!IsOutside){
00707       plnid = plnid.GetAdjoinScint(dir);
00708       if(plnid.IsValid()){
00709         UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(plnid);
00710         z = scintpln.GetZ0();
00711         u = c[0] + m[0]*z;
00712         v = c[1] + m[1]*z;
00713         x = 0.707*(u-v);
00714         y = 0.707*(u+v);
00715         r = pow(x*x+y*y,0.5);
00716         pl.DistanceToEdge(x, y,
00717                           plnid.GetPlaneView(),
00718                           plnid.GetPlaneCoverage(),
00719                           dist, xedge, yedge);
00720         
00721         IsInside = pl.IsInside( x, y,
00722                                 plnid.GetPlaneView(),
00723                                 plnid.GetPlaneCoverage(),true);
00724         IsInside |= r<0.5;
00725         if(IsInside){
00726           nActiveDownstream++;
00727           trace = pow(pow(u-cth.GetEndU(),2)+pow(v-cth.GetEndV(),2)+pow(z-cth.GetEndZ(),2),0.5);
00728           tracez = fabs(z - cth.GetEndZ());
00729         }
00730       }
00731       IsOutside = !IsInside | !plnid.IsValid();
00732     }
00733     if(plnid.IsValid()) trace +=dist;
00734   }
00735   cth.SetEndTrace(trace);
00736   cth.SetEndTraceZ(tracez);
00737   cth.SetEndnActiveDownstream(nActiveDownstream);
00738 
00739   // set distance from track end to nearest edge on that plane
00740   z = cth.GetEndZ();
00741   u = cth.GetEndU();
00742   v = cth.GetEndV();
00743   x = 0.707*(u-v);
00744   y = 0.707*(u+v);
00745   r = pow(x*x+y*y,0.5);
00746        
00747   PlexPlaneId plnidend(detector,cth.GetEndPlane(),false);
00748   dist = 0.;
00749   IsInside=true;
00750   if(plnidend.IsValid()){
00751     pl.DistanceToEdge(x, y,
00752                       plnidend.GetPlaneView(),
00753                       plnidend.GetPlaneCoverage(),
00754                       dist, xedge, yedge);
00755     IsInside = pl.IsInside( x, y,
00756                             plnidend.GetPlaneView(),
00757                             plnidend.GetPlaneCoverage(),true);
00758 
00759   }
00760   // set distance to edge =0 if inside coil.  This forces tracks that end in
00761   // coil to be !contained. 
00762   if(!IsInside && r<0.5) dist=0; 
00763   
00764   cth.SetEndDistToEdge(dist);
00765 
00766 
00767 }
00768 
00769