loon: AlgClusterSRList.cxx Source File

00001 
00002 // $Id: AlgClusterSRList.cxx,v 1.21 2007/11/11 08:26:57 rhatcher Exp $
00003 //
00004 // AlgClusterSRList.cxx
00005 //
00006 // Begin_Html<img src="../../pedestrians.gif" align=center>
00007 // <a href="../source_warning.html">Warning for beginners</a>.<br> 
00008 //
00009 // AlgClusterSRList is a concrete ClusterList Algorithm class.
00010 //
00011 // Author:  R. Lee 2001.02.13
00012 //
00013 // Also see <a href="../../root_crib/index.html">The ROOT Crib</a> and 
00014 // <a href="../CandClusterSR.html"> CandClusterSR Classes</a> (part of
00015 // <a href="../index.html">The MINOS Class User Guide</a>)End_Html
00017 
00018 #include <cassert>
00019 #include <cmath>
00020 
00021 #include "Algorithm/AlgFactory.h"
00022 #include "Algorithm/AlgHandle.h"
00023 #include "Algorithm/AlgConfig.h"
00024 #include "Candidate/CandContext.h"
00025 #include "CandClusterSR/AlgClusterSRList.h"
00026 #include "RecoBase/CandCluster.h"
00027 #include "RecoBase/CandCluster.h"
00028 #include "RecoBase/CandClusterHandle.h"
00029 #include "RecoBase/CandClusterList.h"
00030 #include "RecoBase/CandClusterListHandle.h"
00031 #include "Conventions/PlaneView.h"
00032 #include "MessageService/MsgService.h"
00033 #include "MinosObjectMap/MomNavigator.h"
00034 #include "Navigation/NavKey.h"
00035 #include "Navigation/NavSet.h"
00036 #include "RecoBase/CandClusterHandle.h"
00037 #include "RecoBase/CandSliceHandle.h"
00038 #include "RecoBase/CandSliceListHandle.h"
00039 #include "Validity/VldContext.h"
00040 
00041 ClassImp(AlgClusterSRList)
00042 
00043 //______________________________________________________________________
00044 // CVSID("$Id: AlgClusterSRList.cxx,v 1.21 2007/11/11 08:26:57 rhatcher Exp $");
00045 //______________________________________________________________________
00046 
00047 AlgClusterSRList::AlgClusterSRList()
00048 {
00049 }
00050 
00051 //______________________________________________________________________
00052 AlgClusterSRList::~AlgClusterSRList()
00053 {
00054 }
00055 
00056 //______________________________________________________________________
00057 
00061 void AlgClusterSRList::RunAlg(AlgConfig &ac, CandHandle &ch,CandContext &cx)
00062 {
00063 
00064   Double_t maxTimeDiff = ac.GetDouble("StripNeighborTimeDiff");
00065   Int_t maxStripDiff = ac.GetInt("StripNeighborStripDiff");
00066   Double_t maxTposDiffshower = ac.GetDouble("StripNeighborTPosDiffShower");
00067   Double_t maxTposDifftrack = ac.GetDouble("StripNeighborTPosDiffTrack");
00068   Int_t maxPlaneDiff = ac.GetInt("StripNeighborPlaneDiff");
00069   Double_t minpulseheight = ac.GetDouble("MinPulseHeight");
00070   Double_t minshwstripPE = ac.GetDouble("MinShwStripPE");
00071   Int_t minplanecoverageshower = ac.GetInt("MinPlaneCoverageShower");
00072   Int_t minplanecoveragetrack = ac.GetInt("MinPlaneCoverageTrack");
00073   Int_t minplaneneighbor = ac.GetInt("MinPlaneNeighbor");
00074   Int_t sm1LastPlane = ac.GetInt("SMPlaneLast");
00075   Int_t sm2FirstPlane = ac.GetInt("SMPlaneFirst");
00076   assert(cx.GetDataIn());
00077 
00078 // Check for CandSliceListHandle input
00079   if (!(cx.GetDataIn()->InheritsFrom("CandSliceListHandle"))) {
00080     return;
00081   }
00082 
00083   const CandSliceListHandle *cslh =
00084         dynamic_cast<const CandSliceListHandle*>(cx.GetDataIn());
00085 
00086   CandContext cxx(this,cx.GetMom());
00087   AlgFactory &af = AlgFactory::GetInstance();
00088   AlgHandle ah = af.GetAlgHandle("AlgClusterSR","default");
00089 
00090   CandSliceHandleItr sliceItr(cslh->GetDaughterIterator());
00091   while (CandSliceHandle *slice = sliceItr()) {
00092     TObjArray * sliceClusters = new TObjArray(1,0);    
00093     Int_t planeExtent = slice->GetEndPlane()-slice->GetBegPlane()+1;
00094     
00095     /* Generate the neighbor map, which for each CandStrip in the CandSlice holds the number of neighbor CandStrips, defined to be the number of CandStrips on the same plane having times within StripNeighborTimeDiff, and which lie within StripNeighborStripDiff strips.*/    
00096     GenNeighborMap(ac,*slice);
00097     
00098     // loop over views
00099     
00100     for (Int_t view=PlaneView::kU; view<=PlaneView::kV; view++) {
00101       /* For each non-spectrometer CandStrip, we iterate over previously constructed CandClusters, determining whether the strip should be added to an existing cluster. To be added, the strip must have neighbors, and the cluster must be showerlike. . */
00102       CandStripHandleItr stripItr(slice->GetDaughterIterator());
00103       while (CandStripHandle *strip = stripItr()) {
00104         // don't cluster in muon spectrometer
00105         if (strip->GetPlaneView()==view &&
00106             (cslh->GetVldContext()->GetDetector()!=Detector::kNear ||
00107              strip->GetPlane()<=121)) {
00108           Int_t striphasneighbors(0);
00109           if (fNNeighbors[strip]>=minplaneneighbor) {
00110             striphasneighbors = 1;
00111           }
00112           // iterate over clusters
00113           Int_t nfound = 0;
00114           CandClusterHandle *foundcluster = 0;
00115           TIter clusterItr(sliceClusters);
00116           while (CandClusterHandle *cluster = dynamic_cast<CandClusterHandle *>(clusterItr()) )  {
00117             if (fClusterView[cluster]==view &&
00118                 ((striphasneighbors && 
00119                   cluster->IsShowerLike()) ||
00120                  planeExtent<minplanecoveragetrack)) {
00121               
00122               /*  For each CandCluster, we select daughter CandStrips which lie within +/- StripNeighborPlaneDiff of the outer loop CandStrip.  In the far detector, this interval is extended by SMPlaneFirst-SMPlaneLast-1 (3) planes on the low plane side if the CandStrip plane is StripNeighborPlaneDiff or more lower than the first plane in SM2 (SMPlaneFirst) or by the same number of planes on the high plane side if the CandStrip plane is StripNeighborPlaneDiff or more higher than the last plane in SM1 (SMPlaneLast). */
00123               
00124               CandStripHandleItr clsstripItr(cluster->GetDaughterIterator());
00125               /* If a Cluster CandStrip lies in the same plane as the primary CandStrip, we compare the timing and strip difference  */
00126               Bool_t lfound(0);
00127               for (CandStripHandle *clsstrip = clsstripItr();
00128                    !lfound && clsstrip; clsstrip = clsstripItr()) {
00129                 Double_t dtime = clsstrip->GetBegTime()-strip->GetBegTime();
00130                 if (clsstrip->GetPlane() == strip->GetPlane()) {
00131                   Int_t dstrip = clsstrip->GetStrip()-strip->GetStrip();
00132                   if (fabs(dtime)<=maxTimeDiff && 
00133                       abs(dstrip)<=maxStripDiff) {
00134                     lfound = 1;
00135                   }
00136                 }
00137                 
00138                 /*  If a Cluster CandStrip in a different plane from the primary CandStrip, we use these criteria, and in addition require that the CandStrip lie within a transverse distance of StripNeighborTPosDiffShower if the primary CandStrip was found to have neighbors in the initial step, or within StripNeighborTPosDiffTrack if the CandStrip was found to be isolated, and to be within planediff planes */
00139                 
00140                 else {
00141                   Int_t dplane = clsstrip->GetPlane()-strip->GetPlane();
00142                   Double_t dtpos = strip->GetTPos()-clsstrip->GetTPos();
00143                   Int_t maxPlaneDiffcorr = maxPlaneDiff;
00144                   // if cluster crossed SM boundary, increase allowable plane separation
00145                   if (cslh->GetVldContext()->GetDetector()==Detector::kFar &&
00146                       min(clsstrip->GetPlane(),strip->GetPlane())<=
00147                       sm1LastPlane &&
00148                       max(clsstrip->GetPlane(),strip->GetPlane())>=
00149                       sm2FirstPlane) {
00150                     maxPlaneDiffcorr += sm2FirstPlane-sm1LastPlane-1;
00151                   }
00152                   
00153                   if (fabs(dtime)<=maxTimeDiff && 
00154                       abs(dplane)<=maxPlaneDiffcorr &&
00155                       ((striphasneighbors && fabs(dtpos)<=maxTposDiffshower) ||
00156                        (!striphasneighbors && fabs(dtpos)<=maxTposDifftrack))) {
00157                     lfound = 1;
00158                   }
00159                 }
00160               }
00161               if (lfound) {
00162                 if (!nfound) {            
00163                   // this is the first cluster match - add strip to daughter list
00164                   cluster->AddDaughterLink(*strip);
00165                   foundcluster = cluster;
00166                 }
00167                 else {
00168                   // multiple matchs were found - merge this one with the first
00169                   TObjArray trashbin;
00170                   CandStripHandleItr
00171                     newstripItr(cluster->GetDaughterIterator());
00172                   while (CandStripHandle *newstrip = newstripItr()) {
00173                     foundcluster->AddDaughterLink(*newstrip);
00174                     trashbin.Add(newstrip);
00175                   }
00176                   for (Int_t i=0; i<=trashbin.GetLast(); i++) {
00177                     CandStripHandle *trashstrip =
00178                       dynamic_cast<CandStripHandle*>(trashbin.At(i));
00179                     cluster->RemoveDaughter(trashstrip);
00180                   }
00181                 }
00182                 nfound++;
00183               }
00184             }
00185           }
00186           // remove merged clusters here for efficiency
00187           TObjArray trashbin;
00188           CandClusterHandleItr cluster2Itr(ch.GetDaughterIterator());
00189           while (CandClusterHandle *cluster = cluster2Itr()) {
00190             if (cluster->GetNDaughters()==0) {
00191               trashbin.Add(cluster);
00192               sliceClusters->Remove(cluster);
00193             }
00194           }
00195           sliceClusters->Compress();
00196           for (Int_t i=0; i<=trashbin.GetLast(); i++) {
00197             CandClusterHandle *trashcluster =
00198               dynamic_cast<CandClusterHandle*>(trashbin.At(i));
00199             ch.RemoveDaughter(trashcluster);
00200           }
00201           // no match was found - start a new cluster
00202           //If the primary CandStrip has strip neighbors, the new CandCluster is labeled 'showerlike'.  If not, it is labeled 'tracklike'.
00203           if (!nfound) {
00204             TObjArray striparray;
00205             striparray.Add(strip);
00206             cxx.SetDataIn(&striparray);
00207             CandClusterHandle clusterhandle = 
00208               CandCluster::MakeCandidate(ah,cxx);
00209             clusterhandle.SetCandSlice(slice);
00210             if (striphasneighbors || 
00211                 (planeExtent<minplanecoveragetrack && strip->GetCharge()>minshwstripPE)) {
00212               clusterhandle.IsShowerLike(1);
00213               clusterhandle.IsTrackLike(0);
00214             }
00215             else {
00216               clusterhandle.IsShowerLike(0);
00217               clusterhandle.IsTrackLike(1);
00218             }
00219             ch.AddDaughterLink(clusterhandle);
00220             CandClusterHandle * daucluster = dynamic_cast<CandClusterHandle *>(ch.FindDaughter(&clusterhandle)); 
00221             if(daucluster){
00222               fClusterView[daucluster]=clusterhandle.GetPlaneView();
00223               sliceClusters->Add(daucluster);
00224             }
00225             else{
00226               cout << " new cluster not found on daughter list " << endl;
00227             }
00228           }
00229         }
00230       }
00231     }
00232     delete sliceClusters;
00233   }
00234 
00235   // remove clusters which are below threshold
00236   TObjArray trashbin;
00237   CandClusterHandleItr clusterItr(ch.GetDaughterIterator());
00238   while (CandClusterHandle *cluster = clusterItr()) {
00239     Int_t dplane = cluster->GetEndPlane()-cluster->GetBegPlane()+1;
00240     if (cluster->GetNDaughters()==0 || cluster->GetCharge()<minpulseheight ||
00241         (cluster->IsShowerLike() && dplane<minplanecoverageshower) ||
00242         (cluster->IsTrackLike() && dplane<minplanecoveragetrack)) {
00243       trashbin.Add(cluster);
00244     }
00245   }
00246   for (Int_t i=0; i<=trashbin.GetLast(); i++) {
00247     CandClusterHandle *trashcluster =
00248       dynamic_cast<CandClusterHandle*>(trashbin.At(i));
00249     ch.RemoveDaughter(trashcluster);
00250   }
00251 }
00252 // ----------------------------------------------------------
00253 
00257 void AlgClusterSRList::GenNeighborMap(AlgConfig &ac, CandSliceHandle &csh){
00258   
00259   CandStripHandleItr strip1Itr(csh.GetDaughterIterator());
00260   CandStripHandleKeyFunc *stripKF = strip1Itr.CreateKeyFunc();
00261   stripKF->SetFun(AlgClusterSRList::StripKeyFromPlane);
00262   strip1Itr.GetSet()->AdoptSortKeyFunc(stripKF);
00263   stripKF=0;
00264   
00265   TObjArray planeStripList(1000);
00266   Int_t iplane = 0;
00267   Bool_t first = true;
00268   CandStripHandle * strip;
00269   TObjArray *planeStrips;
00270   while( (strip = strip1Itr())) {
00271     if(strip->GetPlane()<1000){
00272       if(strip->GetPlane()!= iplane || first){
00273         iplane=strip->GetPlane();
00274         planeStrips = new TObjArray(1,0);
00275         planeStripList.AddAt(planeStrips,iplane);
00276         planeStrips->Add(strip);
00277         first = false;
00278       }
00279       else{
00280         planeStrips = (TObjArray *)planeStripList.At(iplane);
00281         if(strip) planeStrips->Add(strip);
00282       }
00283     }
00284   }    
00285   
00286   Double_t maxTimeDiff = ac.GetDouble("StripNeighborTimeDiff");
00287   Int_t maxStripDiff = ac.GetInt("StripNeighborStripDiff");
00288 
00289   strip1Itr.Reset();
00290   while (CandStripHandle *strip1 = strip1Itr()) {
00291     Int_t neighbor = 0;
00292     Int_t dstrip;
00293     Double_t dtime;
00294     TObjArray * planeStrips=(TObjArray *)planeStripList.At(strip1->GetPlane());
00295     TIter strip2Itr(planeStrips);
00296     while (CandStripHandle *strip2 = (CandStripHandle *)(strip2Itr.Next()) )  {
00297       if ((strip1 != strip2) && (strip1->GetPlane()==strip2->GetPlane()) ) {
00298         dtime = strip2->GetBegTime()-strip1->GetBegTime();
00299         dstrip = strip2->GetStrip()-strip1->GetStrip();
00300         if (fabs(dtime)<=maxTimeDiff && abs(dstrip)<=maxStripDiff) {
00301           neighbor++;
00302         }
00303       }
00304     }
00305     fNNeighbors[strip1] = neighbor;
00306   }
00307   for (Int_t iplane = 0;iplane < 1000;iplane++){
00308     planeStrips = (TObjArray *)planeStripList.At(iplane);
00309     if(planeStrips){
00310       delete planeStrips;
00311     }
00312   }
00313 }
00314 
00315 NavKey AlgClusterSRList::StripKeyFromPlane(const CandStripHandle *strip)
00316 {
00317   return const_cast<CandStripHandle *>(strip)->GetPlane();
00318 }
00319 
00320 
00321 
00322 //______________________________________________________________________
00323 void AlgClusterSRList::Trace(const char * /* c */) const
00324 {
00325 }
00326