loon: AlgSliceSRList.cxx Source File

00001 
00002 // $Id: AlgSliceSRList.cxx,v 1.34 2009/10/22 18:21:39 tagg Exp $
00003 //
00004 // AlgSliceSRList.cxx
00005 //
00006 // Begin_Html<img src="../../pedestrians.gif" align=center>
00007 // <a href="../source_warning.html">Warning for beginners</a>.<br> 
00008 //
00009 // AlgSliceSRList is a concrete SliceSRList Algorithm class.
00010 //
00011 // Author:  R. Lee 2001.01.29
00012 // Revised: T. Osiecki 2004.07.09 - osiecki@mail.hep.utexas.edu
00013 // Revised: T. Osiecki 2004.12.08 - "                           
00014 // 
00015 // Also see <a href="../../root_crib/index.html">The ROOT Crib</a> and 
00016 // <a href="../CandSliceSR.html"> CandSliceSR Classes</a> (part of
00017 // <a href="../index.html">The MINOS Class User Guide</a>)End_Html
00019 
00020 #include <cassert>
00021 #include <cmath>
00022 
00023 #include "Algorithm/AlgFactory.h"
00024 #include "Algorithm/AlgHandle.h"
00025 #include "Algorithm/AlgConfig.h"
00026 #include "Candidate/CandContext.h"
00027 #include "CandDigit/CandDigitHandle.h"
00028 #include "CandSliceSR/AlgSliceSRList.h"
00029 #include "RecoBase/CandSlice.h"
00030 #include "RecoBase/CandSlice.h"
00031 #include "RecoBase/CandSliceHandle.h"
00032 #include "RecoBase/CandSliceList.h"
00033 #include "RecoBase/CandSliceListHandle.h"
00034 #include "Conventions/Detector.h"
00035 #include "MessageService/MsgService.h"
00036 #include "MinosObjectMap/MomNavigator.h"
00037 #include "Navigation/NavKey.h"
00038 #include "Navigation/NavSet.h"
00039 #include "RawData/RawDigit.h"
00040 #include "RawData/RawHeader.h"
00041 #include "RawData/RawRecord.h"
00042 #include "RawData/RawChannelId.h"
00043 #include "RawData/RawDigitDataBlock.h"
00044 #include "RecoBase/CandSliceHandle.h"
00045 #include "RecoBase/CandStripHandle.h"
00046 #include "RecoBase/CandStripListHandle.h"
00047 #include "UgliGeometry/UgliGeomHandle.h"
00048 #include "Validity/VldContext.h"
00049 #include "DataUtil/GetRawBlock.h"
00050 #include <vector>
00051 ClassImp(AlgSliceSRList)
00052 
00053 static NavKey SliceSRKeyFromForwardTime(const CandSliceHandle *);
00054 static NavKey StripSRKeyFromCorrTime(const CandStripHandle *);
00055 static NavKey StripSRKeyFromBegTime(const CandStripHandle *);
00056 //______________________________________________________________________
00057 CVSID("$Id: AlgSliceSRList.cxx,v 1.34 2009/10/22 18:21:39 tagg Exp $");
00058 
00059 //______________________________________________________________________
00060 AlgSliceSRList::AlgSliceSRList()
00061 {
00062 }
00063 
00064 //______________________________________________________________________
00065 AlgSliceSRList::~AlgSliceSRList()
00066 {
00067 }
00068 
00069 //______________________________________________________________________
00073 void AlgSliceSRList::RunAlg(AlgConfig &ac, CandHandle &ch,CandContext &cx)
00074 {
00075    assert(cx.GetDataIn());
00076 
00077    Int_t passall  = ac.GetInt("PassAll");
00078    Int_t passasap = ac.GetInt("PassASAP");
00079    Int_t passmst  = ac.GetInt("PassMST");
00080    
00081    // Check for CandStripListHandle input
00082    if (cx.GetDataIn()->InheritsFrom("CandStripListHandle")) {
00083      const MomNavigator *mom = cx.GetMom();
00084 
00085      const RawDigitDataBlock *rddb = 
00086         DataUtil::GetRawBlock<RawDigitDataBlock>(mom);
00087      if (!rddb) {
00088        MSG("SliceSR", Msg::kWarning) << "No RawDigitDataBlock in RawRecord." << endl;
00089        return;
00090      }
00091 
00092      if(passall){
00093        PassAll(ac,ch,cx);
00094      }
00095      if(passasap){
00096        SlicetheSnarl_ASAP(ac,ch,cx);
00097      }
00098      if(passmst){
00099        SlicetheSnarl_MST(ac,ch,cx);
00100      }
00101      else if(!passall && !passasap && !passmst){
00102        SlicetheSnarl(ac,ch,cx);
00103      }
00104    }
00105 }
00106 //______________________________________________________________________________________
00113 void AlgSliceSRList::SlicetheSnarl(AlgConfig &ac, CandHandle &ch,CandContext &cx)
00114 {
00115 
00116   const CandStripListHandle *cslh =
00117         dynamic_cast<const CandStripListHandle*>(cx.GetDataIn());
00118 
00119   MSG("AlgSliceSR",Msg::kDebug) << "Starting SR Slicer!" << endl;
00120 
00121 // Sort by Plane and Strip
00122   CandStripHandleItr cshItr(cslh->GetDaughterIterator());
00123   CandStripHandleKeyFunc *cshKf = cshItr.CreateKeyFunc();
00124   cshKf->SetFun(StripSRKeyFromCorrTime);
00125   cshItr.GetSet()->AdoptSortKeyFunc(cshKf);
00126   cshKf = 0;
00127 
00128   AlgFactory &af = AlgFactory::GetInstance();
00129   AlgHandle ah = af.GetAlgHandle("AlgSliceSR","default");
00130   CandContext cxx(this,cx.GetMom());
00131 
00132   Double_t timewindow           = ac.GetDouble("TimeWindow");
00133   Int_t    minstrip             = ac.GetInt("MinStrip");
00134   Double_t mincharge            = ac.GetDouble("MinCharge");
00135   Double_t maxtimegap           = ac.GetDouble("MaxTimeGap");
00136   Double_t earlytimediff        = ac.GetDouble("EarlyTimeDiff");
00137   Double_t timediffspect        = ac.GetDouble("TimeDiffSpect");
00138   Double_t minpulseheight       = ac.GetDouble("MinPulseHeight");
00139   Double_t zsplitdistance       = ac.GetDouble("ZSplitDistance");
00140   Double_t trackiness           = ac.GetDouble("Trackiness");
00141   Double_t mincombiningdistance = ac.GetDouble("MinCombiningDistance");
00142   Int_t    combiningonoff       = ac.GetInt("CombiningOnOff");
00143   earlytimediff                 = fabs(earlytimediff); // make sure this is a positive quantity
00144   
00145   Int_t specplane = 999;
00146   if (cslh->GetVldContext()->GetDetector() == Detector::kNear) specplane = 121;
00147   
00148   // Iterate over CandStripList daughters
00149   TObjArray striparray;
00150   striparray.Clear();
00151   Int_t ifirst = 0;
00152   Int_t nstrip = 0;
00153   Double_t dtime = 0.;
00154   while (CandStripHandle *curr = cshItr()) {
00155     
00156     MSG("AlgSliceSR",Msg::kDebug) << " strip plane " << curr->GetPlane() << " strip " << curr->GetStrip() << " charge " << curr->GetCharge() << endl;
00157  
00158     if (curr->GetPlane()<specplane && curr->GetCharge()>=mincharge) {
00159       // Only consider upstream of muon spectrometer region and strips of a minimum pulse height
00160       if (!striparray.First()) {
00161         MSG("AlgSliceSR",Msg::kDebug) << " first strip " << endl; 
00162         striparray.Add(curr);
00163         nstrip++;
00164       } else if (curr->GetCorrBegTime()-dynamic_cast<CandStripHandle*> (striparray.Last())->GetCorrBegTime()>maxtimegap) {
00165         // If difference between this TOF corrected time and last TOF corrected
00166         // time in striparray is greater than maxtimegap, form a new candslice
00167         if (striparray.GetLast()+1>=minstrip) {
00168           MSG("AlgSliceSR",Msg::kDebug) << " starting new slice " << endl; 
00169           cxx.SetDataIn(&striparray);
00170           CandSliceHandle slicehandle =
00171             CandSlice::MakeCandidate(ah,cxx);
00172           ch.AddDaughterLink(slicehandle);
00173         }
00174         striparray.Clear();
00175         striparray.Add(curr);
00176         ifirst = 0;
00177         nstrip = 1;
00178       } else {
00179         if (curr->GetCorrBegTime()<=dynamic_cast<CandStripHandle*>
00180             (striparray.At(ifirst))->GetCorrBegTime()+timewindow) {
00181           // If difference between this TOF corrected time and first TOF corrected
00182           // time in striparray is less than timewindow, add to striparray
00183           MSG("AlgSliceSR",Msg::kDebug) << " adding strip " << endl; 
00184           striparray.Add(curr);
00185           nstrip++;
00186           dtime = curr->GetCorrBegTime()-dynamic_cast<CandStripHandle*>
00187             (striparray.At(0))->GetCorrBegTime();
00188           // dtime is the time extend of striparray
00189         } else {
00190           /*  Gap is less than MaxTimeGap, but the total list time duration is greater than TimeWindow.  In this case we assume that we have missed the true start of this CandSlice. We now check to see if, by removing the first entry on the list, and adding the new entry, we end up with a shorter time duration in the list.
00191            */ 
00192           if (striparray.At(ifirst+1) && 
00193               curr->GetCorrBegTime()-dynamic_cast<CandStripHandle*>
00194               (striparray.At(ifirst+1))->GetCorrBegTime()<dtime &&
00195               striparray.GetLast()+1-ifirst>=nstrip) {
00196             MSG("AlgSliceSR",Msg::kDebug) << " removing first/adding " << endl; 
00197             // Removing first strip and adding last results in a shorter slice, so assume this is the thing to do    
00198             striparray.RemoveAt(0);
00199             striparray.Compress();
00200             striparray.Add(curr);
00201             nstrip = striparray.GetLast()-ifirst;
00202             dtime = curr->GetCorrBegTime()-dynamic_cast<CandStripHandle*> (striparray.At(0))->GetCorrBegTime();
00203           } else {
00204             // Otherwise, make a new slice          
00205             if (striparray.GetLast()+1>=minstrip) {           
00206               MSG("AlgSliceSR",Msg::kDebug) << " make new slice  " << endl; 
00207               cxx.SetDataIn(&striparray);
00208               CandSliceHandle slicehandle =
00209                 CandSlice::MakeCandidate(ah,cxx);
00210               ch.AddDaughterLink(slicehandle);
00211             }
00212             // Prepare for next slice
00213             striparray.Clear();
00214             striparray.Add(curr);
00215             ifirst = 0;
00216             nstrip = 1;
00217           }
00218         }
00219       }
00220     }
00221   }
00222   
00223   // Cleanup - make slice out of leftovers
00224   if (striparray.GetLast()+1>=minstrip) {
00225   MSG("AlgSliceSR",Msg::kDebug) << " final slice make " << endl; 
00226     cxx.SetDataIn(&striparray);
00227     CandSliceHandle slicehandle = CandSlice::MakeCandidate(ah,cxx);
00228     ch.AddDaughterLink(slicehandle);
00229   }
00230   
00231   /* Fix eaten slices, doing it with Z first (only NEAR DET). Start off with 1st strip hit and find everything within a distance of it, for all hits, put left overs in a separate array, now have 2 slices.  In the slice that is further in Z (presumably a track), one checks how many planes have > 2 hits, if it truly is a track then it should have about 1 hit / plane.  Counting the number of planes that have at least a hit and those that have > 2 hits, I take the ratio, if it is less than 0.05 split, otherwise do not.
00232   */
00233 
00234   if (cslh->GetVldContext()->GetDetector() == Detector::kNear){
00235     MSG("AlgSliceSR",Msg::kDebug) << "Z Splitting!" << endl;
00236     CandSliceHandleItr sliceItr(ch.GetDaughterIterator());
00237     CandSliceHandleKeyFunc *sliceKf = sliceItr.CreateKeyFunc();
00238     sliceKf->SetFun(SliceSRKeyFromForwardTime);
00239     sliceItr.GetSet()->AdoptSortKeyFunc(sliceKf);
00240     sliceKf = 0;
00241     
00242     sliceItr.ResetLast();
00243     TObjArray z_split;
00244     z_split.Clear();
00245     TObjArray z_split2;
00246     z_split2.Clear();
00247     while (CandSliceHandle *curr = sliceItr()) {
00248       z_split.Clear(); z_split2.Clear();
00249       CandStripHandleItr currstripItr(curr->GetDaughterIterator());
00250       while(CandStripHandle *strip3 = currstripItr()) {
00251         z_split.Add(strip3);
00252       }
00253       CandStripHandle *temp = dynamic_cast<CandStripHandle*>(z_split.At(0));
00254       z_split2.Add(temp); z_split.RemoveAt(0); z_split.Compress();
00255       while(z_split.GetEntries() > 0) {
00256         Int_t sa_size = z_split.GetEntries();
00257         for(int i=0; i<=z_split2.GetLast(); i++) {
00258           CandStripHandle *sa2 = dynamic_cast<CandStripHandle*>(z_split2.At(i));
00259           for(int j=0; j<=z_split.GetLast(); j++) {
00260             CandStripHandle *sa = dynamic_cast<CandStripHandle*>(z_split.At(j));
00261             if(fabs(sa->GetZPos()-sa2->GetZPos()) < zsplitdistance) {
00262               z_split2.Add(sa);
00263               z_split.RemoveAt(j);
00264               z_split.Compress();
00265             }
00266           }
00267         }
00268         Int_t sa_size2 = z_split.GetEntries();
00269         if(sa_size2 == sa_size) {
00270           break;
00271         }
00272       }
00273       
00274       //Now need to check if I'm not splitting a shower and a track that went through the coil!
00275       int planer[300] = {0}; int planer2[300] = {0};
00276       int bigplane = 0, bigplane2 = 0;
00277       for (Int_t i=0; i<=z_split.GetLast(); i++) {
00278         CandStripHandle *csh = dynamic_cast<CandStripHandle*>(z_split.At(i));
00279         planer[csh->GetPlane()]++;
00280         if(csh->GetPlane() > bigplane) bigplane = csh->GetPlane();
00281       }
00282       for (Int_t i=0; i<=z_split2.GetLast(); i++) {
00283         CandStripHandle *csh = dynamic_cast<CandStripHandle*>(z_split2.At(i));
00284         planer2[csh->GetPlane()]++;
00285         if(csh->GetPlane() > bigplane2) bigplane2 = csh->GetPlane();
00286       }
00287       
00288       int isbad = 0, nump  = 0;
00289       if(bigplane > bigplane2) {
00290         for(int i=0; i<300; i++) {
00291           if(planer[i]>0) nump++;
00292           if(planer[i]>2) isbad++;
00293         }
00294         if(nump == 0) {
00295           nump = 100; isbad = 1;
00296         }
00297       } else {
00298         for(int i=0; i<300; i++) {
00299           if(planer2[i]>0) nump++;
00300           if(planer2[i]>2) isbad++;
00301         }
00302         if(nump == 0) {
00303           nump = 100; isbad = 1;
00304         }
00305       }
00306             
00307       if((float) isbad / nump > trackiness) { 
00308         if(z_split2.GetEntries() >= minstrip) {
00309           striparray.Clear();
00310           for(int i=0; i<=z_split2.GetLast(); i++) {
00311             CandStripHandle *temp = dynamic_cast<CandStripHandle*>(z_split2.At(i));
00312             striparray.Add(temp);
00313           } 
00314           cxx.SetDataIn(&striparray);
00315           CandSliceHandle slicehandle = CandSlice::MakeCandidate(ah,cxx);
00316           ch.AddDaughterLink(slicehandle); 
00317         }
00318         if(z_split.GetEntries() >= minstrip) {
00319           striparray.Clear();
00320           for(int i=0; i<=z_split.GetLast(); i++) {
00321             CandStripHandle *temp = dynamic_cast<CandStripHandle*>(z_split.At(i));
00322             striparray.Add(temp);
00323           } 
00324           cxx.SetDataIn(&striparray);
00325           CandSliceHandle slicehandle = CandSlice::MakeCandidate(ah,cxx);
00326           ch.AddDaughterLink(slicehandle);
00327         }
00328         ch.RemoveDaughter(curr);
00329       }
00330       z_split.Clear(); z_split2.Clear(); striparray.Clear();    
00331     }
00332   }
00333   
00334   // Prepare slice iterator for next phase
00335   CandSliceHandleItr sliceItr2(ch.GetDaughterIterator());
00336   CandSliceHandleKeyFunc *sliceKf2 = sliceItr2.CreateKeyFunc();
00337   sliceKf2->SetFun(SliceSRKeyFromForwardTime);
00338   sliceItr2.GetSet()->AdoptSortKeyFunc(sliceKf2);
00339   sliceKf2 = 0;
00340 
00341   /* We now pick up CandStrips in the spectrometer. To do this, we iterate over CandSlices constructed in the previous stage, and determine which remaining CandStrips should be added. For each CandSlice, a time interval bounded by the CandSlice start time minus EarlyTimeDiff and the CandSlice start time plus TimeDiffSpect is defined, and CandStrips in the spectrometer within this time range are added to the slice */
00342 
00343   if (cslh->GetVldContext()->GetDetector() == Detector::kNear){
00344     MSG("AlgSliceSR",Msg::kDebug) << "Adding Spectro Hits!" << endl;
00345     sliceItr2.Reset();
00346     while (CandSliceHandle *slice = sliceItr2()) {
00347       cshItr.GetSet()->Slice();
00348       cshItr.GetSet()->
00349         Slice(slice->GetCorrTime()-earlytimediff,
00350               slice->GetCorrTime()+timediffspect);
00351       while (CandStripHandle *strip1 = cshItr()) {
00352         if (strip1->GetPlane()>= specplane){
00353           slice->AddDaughterLink(*strip1);
00354         }   
00355       }
00356     }
00357   }
00358 
00359   /* we next pick up low pulse height strips. To do this, we iterate over CandSlices constructed in the previous stage, and determine which remaining CandStrips should be added. For each CandSlice, a time interval bounded by the CandSlice start time and the start time of the NEXT CandSlice is defined, and CandStrips with q<mincharge within this time range are added to the slice */
00360 
00361   MSG("AlgSliceSR",Msg::kDebug) << "Adding Small Hits" << endl;
00362   Double_t SliceTime[100] = {0.0};
00363   Int_t counter = 0;
00364   sliceItr2.Reset();
00365   while (CandSliceHandle *slice = sliceItr2()) {
00366     SliceTime[counter] = slice->GetCorrBegTime()*1e9;
00367     counter++;
00368   }
00369   SliceTime[counter] = 30000.0;
00370   counter = 0;
00371   sliceItr2.Reset();
00372   while (CandSliceHandle *slice = sliceItr2()) {
00373     MSG("AlgSliceSR",Msg::kDebug) << " ## slice ## " << counter << " time " << SliceTime[counter] << " next time " << SliceTime[counter+1] <<  endl;
00374     Float_t ticker = fabs(SliceTime[counter+1]/1e9-SliceTime[counter]/1e9);
00375     cshItr.GetSet()->Slice();
00376     cshItr.GetSet()->Slice(slice->GetCorrBegTime()-(earlytimediff/2.0), slice->GetCorrBegTime()+ticker-(earlytimediff/2.0));
00377     while (CandStripHandle *strip = cshItr()) {
00378       MSG("AlgSliceSR",Msg::kDebug) << " strip plane " << strip->GetPlane() << " strip " << strip->GetStrip() << " charge " << strip->GetCharge() << " time " << strip->GetCorrBegTime()*1e9 <<  endl;
00379       if (strip->GetPlane()<specplane && strip->GetCharge() < mincharge) {
00380         slice->AddDaughterLink(*strip);
00381         MSG("AlgSliceSR",Msg::kDebug) << " adding " << endl;    
00382       } 
00383     }
00384     counter++;
00385   }
00386   
00387   /* Here we get rid of slices with  < 2000 ADC counts (mostly junk anyway) */
00388 
00389   MSG("AlgSliceSR",Msg::kDebug) << "Killing Small Slices" << endl;
00390   sliceItr2.Reset();
00391   while( CandSliceHandle *slice = sliceItr2()) {
00392     if(slice->GetCharge(CalDigitType::kNone) < minpulseheight) {
00393       ch.RemoveDaughter(slice);
00394     }
00395   }
00396 
00397   // Prepare Slice iterator for next phase.
00398   CandSliceHandleItr sliceItr3(ch.GetDaughterIterator());
00399   CandSliceHandleKeyFunc *sliceKf3 = sliceItr3.CreateKeyFunc();
00400   sliceKf3->SetFun(SliceSRKeyFromForwardTime);
00401   sliceItr3.GetSet()->AdoptSortKeyFunc(sliceKf3);
00402   sliceKf3 = 0;
00403 
00404   /* Combining Step.  Calculate the Center of Pulse height for all slices < 10,000 ADC counts and loop over all slices that occur before said slice in time and have > 10,000 ADC counts.  Find the min distance between the small slice and the bigger slice.  If the dist in either view is less than a parameter (mincombiningdistance) then combine the 2 slices.  This significantly reduces the number of low-completeness slices, while not affecting the high-completeness regime negatively.  One can see the improvement on slices, but CandEventSR is at a stage where this doesn't help 'yet'.  It can be turned on or off.  12-15-04 */
00405 
00406   if(combiningonoff) {
00407     MSG("AlgSliceSR",Msg::kDebug) << "Starting Combining!" << endl;
00408     Int_t badindex = 0;
00409     TObjArray stripu; stripu.Clear();
00410     TObjArray stripv; stripv.Clear();
00411     while(CandSliceHandle *bad = sliceItr3()) {      
00412       stripu.Clear(); stripv.Clear();
00413       if(bad->GetCharge(CalDigitType::kNone) < 10000) {
00414         CandStripHandleItr stripitr(bad->GetDaughterIterator());
00415         Float_t BegTimeBad = 30000.0;
00416         while(CandStripHandle *curr = stripitr()) {
00417           if(curr->GetTime()*1e9 < BegTimeBad) BegTimeBad = curr->GetTime()*1e9;
00418           if(curr->GetPlane()<121 && curr->GetCharge()>2.0) {
00419             if(curr->GetPlaneView() == 2) {
00420               stripu.Add(curr);
00421             } else {
00422               stripv.Add(curr);
00423             }
00424           }
00425         }
00426         //Calculating Center of Pulse Height Coordinates
00427         Float_t phtotuz = 0.0, phtotvz = 0.0, phu = 0.0, phv = 0.0, phzu = 0.0, phzv = 0.0; 
00428         for(int i=0; i<stripu.GetEntriesFast(); i++) {
00429           CandStripHandle *strip = dynamic_cast<CandStripHandle*>(stripu.At(i));
00430           phtotuz += strip->GetCharge(CalDigitType::kNone);
00431           phu     += strip->GetCharge(CalDigitType::kNone)*strip->GetTPos();
00432           phzu    += strip->GetCharge(CalDigitType::kNone)*strip->GetZPos();
00433         }
00434         for(int i=0; i<stripv.GetEntriesFast(); i++) {
00435           CandStripHandle *strip = dynamic_cast<CandStripHandle*>(stripv.At(i));
00436           phtotvz += strip->GetCharge(CalDigitType::kNone);
00437           phv     += strip->GetCharge(CalDigitType::kNone)*strip->GetTPos();
00438           phzv    += strip->GetCharge(CalDigitType::kNone)*strip->GetZPos();
00439         }
00440         
00441         Float_t Upos = -10000.0, ZUpos = -10000.0, Vpos = -10000.0, ZVpos = -10000.0;
00442         if(phtotuz > 0) {
00443           Upos  = phu / phtotuz; 
00444           ZUpos = phzu / phtotuz;
00445         }
00446         if(phtotvz > 0) {
00447           Vpos  = phv / phtotvz; 
00448           ZVpos = phzv / phtotvz;
00449         }
00450         
00451         Float_t minu = 100000.0, minv = 100000.0; Int_t indexu = 0, indexv = 0;
00452         Float_t BegTimeParent = 30000.0;
00453         Int_t goodindex = 0;
00454         Double_t phdadu = -10.0, phdadv =  -10.0;
00455         
00456         CandSliceHandleItr sliceItr4(ch.GetDaughterIterator());
00457         CandSliceHandleKeyFunc *sliceKf4 = sliceItr4.CreateKeyFunc();
00458         sliceKf4->SetFun(SliceSRKeyFromForwardTime);
00459         sliceItr4.GetSet()->AdoptSortKeyFunc(sliceKf4);
00460         sliceKf4 = 0; 
00461         
00462         sliceItr4.Reset();
00463         while(CandSliceHandle *parent = sliceItr4()) {
00464           if(goodindex != badindex) { 
00465             if(parent->GetCharge(CalDigitType::kNone) >= 10000) {
00466               CandStripHandleItr parentitr(parent->GetDaughterIterator());
00467               while(CandStripHandle *parentstrip = parentitr()) {
00468                 if(parentstrip->GetTime()*1e9 < BegTimeParent) 
00469                   BegTimeParent = parentstrip->GetTime()*1e9;
00470               }
00471               if(BegTimeParent < BegTimeBad) {
00472                 stripu.Clear(); stripv.Clear();
00473                 parentitr.Reset();
00474                 while(CandStripHandle *curr = parentitr()) {
00475                   if(curr->GetPlane()<121 && curr->GetCharge()>2.0) {
00476                     if(curr->GetPlaneView() == 2) {
00477                       stripu.Add(curr);
00478                     } else {
00479                       stripv.Add(curr);
00480                     }
00481                   }
00482                 }
00483                 
00484                 Float_t tempu = 300000.0, tempv = 300000.0;
00485                                 
00486                 for(int i=0; i<stripu.GetEntriesFast(); i++) {
00487                   CandStripHandle *strip = dynamic_cast<CandStripHandle*>(stripu.At(i));
00488                   if(sqrt(pow(Upos-strip->GetTPos(),2)+pow(ZUpos-strip->GetZPos(),2)) < tempu) 
00489                     tempu = sqrt(pow(Upos-strip->GetTPos(),2)+pow(ZUpos-strip->GetZPos(),2));
00490                 }
00491                 for(int i=0; i<stripv.GetEntriesFast(); i++) {
00492                   CandStripHandle *strip = dynamic_cast<CandStripHandle*>(stripv.At(i));
00493                   if(sqrt(pow(Vpos-strip->GetTPos(),2)+pow(ZVpos-strip->GetZPos(),2)) < tempv) 
00494                     tempv = sqrt(pow(Vpos-strip->GetTPos(),2)+pow(ZVpos-strip->GetZPos(),2));
00495                 }
00496                 stripu.Clear(); stripv.Clear();
00497                 if(tempu < minu) {
00498                   minu = tempu;
00499                   indexu = goodindex;
00500                   phdadu = parent->GetCharge(CalDigitType::kNone);
00501                 }
00502                 if(tempv < minv) {
00503                   minv = tempv;
00504                   indexv = goodindex;
00505                   phdadv = parent->GetCharge(CalDigitType::kNone);
00506                 }
00507                 
00508               } 
00509             }  //Loop over slices with ph > 10k, trying to find parent
00510           }
00511           goodindex++;
00512         } //Loop over all slices
00513         
00514         Float_t minofuandv = 0.0;  Int_t minofindexuandv = 0;
00515         if(minu < minv) {
00516           minofuandv = minu;
00517           minofindexuandv = indexu;
00518         } else {
00519           minofuandv = minv;
00520           minofindexuandv = indexv;
00521         }
00522         
00523         if(minofuandv <= mincombiningdistance) {
00524           CandSliceHandleItr sliceItr5(ch.GetDaughterIterator());
00525           CandSliceHandleKeyFunc *sliceKf5 = sliceItr5.CreateKeyFunc();
00526           sliceKf5->SetFun(SliceSRKeyFromForwardTime);
00527           sliceItr5.GetSet()->AdoptSortKeyFunc(sliceKf5);
00528           sliceKf5 = 0; 
00529           
00530           goodindex = 0;
00531           sliceItr5.Reset();
00532           while(CandSliceHandle *slice = sliceItr5()) {
00533             if(goodindex == minofindexuandv) {
00534               CandStripHandleItr badstripitr(slice->GetDaughterIterator());
00535               while(CandStripHandle *badstrip = badstripitr()) {
00536                 slice->AddDaughterLink(*badstrip);
00537               }
00538               ch.RemoveDaughter(bad);
00539               break;
00540             }
00541             goodindex++;
00542           }
00543         } //Combining Loop
00544       } //Loop over slices with ph < 10k, trying to find low-C slices
00545       badindex++;
00546     } //end loop over all slices
00547   } // Combining On if statement
00548   CandSliceHandleItr sliceItr10(ch.GetDaughterIterator());
00549   Int_t nslice=0;
00550    MSG("AlgSliceSR",Msg::kDebug) << "Final Slice Contents " << endl;
00551    while(CandSliceHandle *slice = sliceItr10()) {
00552      MSG("AlgSliceSR",Msg::kDebug) << "  ** Slice **  " << nslice << endl;
00553      CandStripHandleItr liststripitr(slice->GetDaughterIterator());
00554       while(CandStripHandle *strip = liststripitr()) {
00555         MSG("AlgSliceSR",Msg::kDebug) << "strip plane  " << strip->GetPlane() << " strip " << strip->GetStrip() << " charge " << strip->GetCharge() <<" time " << strip->GetCorrBegTime()*1e9 <<  endl;
00556       }
00557       nslice++;
00558    }
00559 
00560 
00561 } //End Slice the Snarl
00562 
00563 //______________________________________________________________________ 
00570 void AlgSliceSRList::SlicetheSnarl_ASAP(AlgConfig &ac, CandHandle &ch,CandContext &cx) // NEED TO DEFINE THAT
00571 {
00572 
00573   cout << " IN SLICE ASAP " <<endl;
00574   const CandStripListHandle *cslh =
00575         dynamic_cast<const CandStripListHandle*>(cx.GetDataIn());
00576 
00577 // Sort by Plane and Strip
00578   CandStripHandleItr cshItr(cslh->GetDaughterIterator());
00579   CandStripHandleKeyFunc *cshKf = cshItr.CreateKeyFunc();
00580   cshKf->SetFun(StripSRKeyFromBegTime);
00581   cshItr.GetSet()->AdoptSortKeyFunc(cshKf);
00582   cshKf = 0;
00583 
00584   AlgFactory &af = AlgFactory::GetInstance();
00585   AlgHandle ah = af.GetAlgHandle("AlgSliceSR","default");
00586   CandContext cxx(this,cx.GetMom());
00587 
00588   Double_t nbuck = ac.GetDouble("NBuck"); 
00589   Int_t minstrip = ac.GetInt("MinStrip");
00590   Double_t mincharge = ac.GetDouble("MinCharge");
00591 
00592   TObjArray striparraysl;
00593   striparraysl.Clear();
00594   
00595   Double_t prevtime = 0;
00596   Bool_t first=true;
00597   Double_t timediff=0;
00598   Double_t timeu;
00599   Int_t ncount=-1;
00600   Int_t nentries=-1;
00601   Double_t cuttime=nbuck*19.*Munits::ns;
00602  // NBuck is be number of buckets 
00603  while(CandStripHandle *currstrip = cshItr()){ 
00604     if(currstrip->GetCharge()>=mincharge) nentries=nentries+1;
00605   }
00606   
00607   cshItr.Reset();
00608    
00609  
00610    
00611     while(CandStripHandle *currstrip = cshItr()){   
00612       if(currstrip->GetCharge()>=mincharge){ 
00613         ncount=ncount+1;
00614         if(first) {
00615          first=false;
00616          prevtime=currstrip->GetBegTime();      
00617         }
00618         timeu=currstrip->GetBegTime();          
00619         timediff  = timeu-prevtime;     
00620           if(timediff<cuttime){ 
00621            striparraysl.Add(currstrip);    
00622           }
00623           if(timediff>=cuttime || nentries==ncount){
00624             if(striparraysl.GetLast()+1>=minstrip){
00625             cxx.SetDataIn(&striparraysl); 
00626             CandSliceHandle slicehandle = CandSlice::MakeCandidate(ah,cxx);
00627             ch.AddDaughterLink(slicehandle);  
00628             }
00629            striparraysl.Clear();
00630            striparraysl.Add(currstrip);
00631            }         
00632          prevtime  = timeu;
00633       }
00634     }
00635 
00636 }
00637 
00638 //_________________________________________________________________________
00646 void AlgSliceSRList::SlicetheSnarl_MST(AlgConfig &ac, CandHandle &ch,CandContext &cx)
00647 {
00648   cout << " IN SLICE MST " << endl;
00649   const CandStripListHandle *cslh =
00650         dynamic_cast<const CandStripListHandle*>(cx.GetDataIn());
00651 // Sort by Plane and Strip
00652   CandStripHandleItr cshItr(cslh->GetDaughterIterator());
00653   CandStripHandleKeyFunc *cshKf = cshItr.CreateKeyFunc();
00654   cshKf->SetFun(StripSRKeyFromCorrTime);
00655   cshItr.GetSet()->AdoptSortKeyFunc(cshKf);
00656   cshKf = 0;
00657 
00658   AlgFactory &af = AlgFactory::GetInstance();
00659   AlgHandle ah = af.GetAlgHandle("AlgSliceSR","default");
00660   CandContext cxx(this,cx.GetMom());
00661  
00662   Int_t minstrip     = ac.GetInt("MinStrip");
00663   Double_t mincharge = ac.GetDouble("MinCharge");
00664   Double_t zfact     = ac.GetDouble("Zfact");
00665   Double_t maxlen    = ac.GetDouble("MaxLen");
00666 
00667   cout << " zfact " << zfact << " maxlen " << maxlen << endl;
00668   TObjArray striparrayn;
00669   striparrayn.Clear();
00670   Int_t n_strip=-1;
00671   while(CandStripHandle *currstrip = cshItr()){
00672     if(currstrip->GetCharge()>=mincharge){
00673       striparrayn.Add(currstrip);
00674       n_strip++;
00675     }
00676   }
00677   
00678   // Create lower triangular matrix of MST distances.
00679   
00680   Int_t n=-1;
00681   Double_t cc=3.*1.e8;     
00682   Int_t kMax_d=((n_strip+1)*n_strip)/2;
00683   Double_t *d_mst = new Double_t[kMax_d];
00684   for(Int_t i=0;i<=(n_strip-1);i++){ // 
00685     for(Int_t j=0;j<=n_strip;j++) {
00686       if((i+1)>j) {
00687         n++;
00688         Double_t zpos1=dynamic_cast<CandStripHandle*>(striparrayn.At(i+1))->GetZPos();
00689         Double_t zpos2=dynamic_cast<CandStripHandle*>(striparrayn.At(j))->GetZPos();
00690         Double_t t1=dynamic_cast<CandStripHandle*>(striparrayn.At(i+1))->GetCorrBegTime(); 
00691         Double_t t2=dynamic_cast<CandStripHandle*>(striparrayn.At(j))->GetCorrBegTime();
00692         d_mst[n]=sqrt((t1*cc-t2*cc)*(t1*cc-t2*cc))+((zpos1-zpos2)*(zpos1-zpos2)*zfact*zfact);   
00693       }
00694     }
00695   }
00696   //    
00697   // Create the MST
00698   Bool_t    *a_mst      = new Bool_t   [n_strip+1];
00699   Int_t     *b_mst      = new Int_t [n_strip+1];
00700   Double_t  *c_mst      = new Double_t [n_strip+1];
00701   Int_t     *hist_mst   = new Int_t [n_strip+1];
00702   Int_t     *route_mst  = new Int_t [n_strip+1];
00703   Int_t     *bb_mst     = new Int_t [n_strip+1];
00704   Int_t     *ii_mst     = new Int_t [n_strip+1];
00705   Double_t  *cc_mst     = new Double_t [n_strip+1];
00706   
00707   Double_t dlarge=99999999.;
00708   Double_t am,dist;
00709   Int_t next = 0;
00710   // Initialization
00711   for(Int_t i=1;i<=n_strip;i++){
00712     a_mst[i]=true;
00713     b_mst[i]=0;
00714     c_mst[i]=dlarge;
00715   } 
00716   // Make the tree
00717   Int_t j=0;
00718   Int_t l;
00719   for(Int_t i=1;i<=n_strip;i++){ 
00720     am=dlarge; 
00721     for(Int_t k=1;k<=n_strip;k++){
00722       if(a_mst[k]) {
00723         if(j>k) l=(j-0)*(j-1)/2+k;
00724         if(j<=k) l=(k-0)*(k-1)/2+j;
00725         dist=d_mst[l];
00726         if(dist<c_mst[k]) {
00727           c_mst[k]=dist;
00728           b_mst[k]=j;
00729         }
00730         if(am>c_mst[k]){ 
00731           am=c_mst[k];
00732           next=k;
00733         }
00734       }  
00735     }
00736     j=next;
00737     a_mst[j]=false;
00738   }
00739   
00740   delete [] d_mst;
00741   delete [] a_mst;
00742   
00743   // DO MST REORDERING
00744   
00745   for (Int_t i=0;i<=n_strip;i++){
00746     hist_mst[i]=0; 
00747   }
00748   for (Int_t i=1;i<=n_strip;i++){
00749     j=b_mst[i];
00750     hist_mst[j]=hist_mst[j]+1; 
00751   } 
00752   route_mst[0]=0;
00753   j=0;
00754   Int_t k=0;
00755   Int_t kkk=0;
00756   for(Int_t i=1;i<=n_strip;i++){
00757     while(hist_mst[k]==0){
00758       j=j-1;
00759       k=route_mst[j]; 
00760     }
00761     hist_mst[k]=hist_mst[k]-1;   
00762     
00763     Int_t m=1;
00764     Bool_t foundk=false;
00765     while(m<=n_strip && !foundk){
00766       if(k==b_mst[m]) {
00767         foundk=true;
00768       }
00769       else m=m+1; 
00770     }
00771     kkk=kkk+1;
00772     bb_mst[kkk]=k;
00773     ii_mst[kkk]=m;
00774     cc_mst[kkk]=c_mst[m];
00775     
00776     j=j+1;
00777     route_mst[j]=m;
00778     k=m;
00779     b_mst[m]=-b_mst[m];
00780     if(b_mst[m]==0) b_mst[m]=9999999;
00781   }
00782   delete [] b_mst;
00783   delete [] c_mst;
00784   delete [] hist_mst;
00785   delete [] route_mst;
00786   
00787   // NOW CUT BRANCHES AND MAKE THE CLUSTERS
00788   
00789   TObjArray striparraysl;
00790   striparraysl.Clear();
00791   
00792   for(Int_t i=1;i<=n_strip;i++){
00793     if(cc_mst[i]<maxlen) {
00794       if(i==1) {
00795         striparraysl.Add(dynamic_cast<CandStripHandle*>(striparrayn.At(ii_mst[i]))); 
00796         striparraysl.Add(dynamic_cast<CandStripHandle*>(striparrayn.At(bb_mst[i]))); 
00797       }
00798       else  striparraysl.Add(dynamic_cast<CandStripHandle*>(striparrayn.At(ii_mst[i])));  
00799     }
00800     else if(cc_mst[i]>=maxlen || i==n_strip) {
00801       if(i==1) {
00802         striparraysl.Add(dynamic_cast<CandStripHandle*>(striparrayn.At(ii_mst[i])));
00803       }
00804       else {
00805         // new slice
00806         if(striparraysl.GetLast()+1>=minstrip){    
00807           cxx.SetDataIn(&striparraysl);
00808           CandSliceHandle slicehandle = CandSlice::MakeCandidate(ah,cxx);
00809           ch.AddDaughterLink(slicehandle);  
00810         }
00811         striparraysl.Clear();  
00812         striparraysl.Add(dynamic_cast<CandStripHandle*>(striparrayn.At(ii_mst[i]))); 
00813       } 
00814     }
00815   }
00816   
00817   delete [] bb_mst;
00818   delete [] cc_mst;
00819   delete [] ii_mst;
00820 }
00821 
00822 //______________________________________________________________________
00824 void AlgSliceSRList::PassAll(AlgConfig & /* ac */, CandHandle &ch,CandContext &cx)
00825 {
00826 
00827    const CandStripListHandle *cslh =
00828          dynamic_cast<const CandStripListHandle*>(cx.GetDataIn());
00829 
00830    AlgFactory &af = AlgFactory::GetInstance();
00831    AlgHandle ah = af.GetAlgHandle("AlgSliceSR","default");
00832    CandContext cxx(this,cx.GetMom());
00833 
00834    CandStripHandleItr cshItr(cslh->GetDaughterIterator());
00835 
00836 // Iterate over CandStripList daughters
00837    TObjArray striparray;
00838    striparray.Clear();
00839    while (CandStripHandle *curr = cshItr()) {
00840       striparray.Add(curr);
00841    }
00842    cxx.SetDataIn(&striparray);
00843    CandSliceHandle slicehandle = CandSlice::MakeCandidate(ah,cxx);
00844    ch.AddDaughterLink(slicehandle);
00845 }
00846 
00847 
00848 //______________________________________________________________________
00849 void AlgSliceSRList::Trace(const char * /* c */) const
00850 {
00851 }
00852 
00853 NavKey SliceSRKeyFromForwardTime(const CandSliceHandle *csh)
00854 {
00855 //
00856 //  Purpose:  Generate a sort key based on the beginning time of CandSlice.
00857 //
00858 //  Arguments:
00859 //    csh        in    Handle to CandSlice to be sorted.
00860 //
00861 //  Return:   Begin time
00862 //
00863   Double_t time = (Double_t)csh->GetCorrBegTime();
00864   return time;
00865 }
00866 
00867 NavKey StripSRKeyFromCorrTime(const CandStripHandle *csh)
00868 {
00869 //
00870 //  Purpose:  Generate a sort key based on the beginning time of CandSlice.
00871 //
00872 //  Arguments:
00873 //    csh        in    Handle to CandSlice to be sorted.
00874 //
00875 //  Return:   Begin time
00876 //
00877 
00878    Double_t time = (Double_t)csh->GetCorrBegTime();
00879    return time;
00880 }
00881 
00882 NavKey StripSRKeyFromBegTime(const CandStripHandle *csh)
00883 {
00884 //
00885 //  Purpose:  Generate a sort key based on the beginning time of CandSlice.
00886 //
00887 //  Arguments:
00888 //    csh        in    Handle to CandSlice to be sorted.
00889 //
00890 //  Return:   Begin time
00891 //
00892 
00893    Double_t time = (Double_t)csh->GetBegTime(); 
00894    return time;
00895 }
00896 
00897 NavKey StripSRKeyFromPlane(const CandStripHandle *csh)
00898 {
00899 //
00900 //  Purpose:  Generate a sort key based on the beginning time of CandSlice.
00901 //
00902 //  Arguments:
00903 //    csh        in    Handle to CandSlice to be sorted.
00904 //
00905 //  Return:   Begin time
00906 //
00907 
00908 
00909    Int_t iplane = const_cast<CandStripHandle *>(csh)->GetPlane();
00910    return iplane;
00911 
00912 }