loon: AlgShowerSRList.cxx Source File

00001 
00002 // $Id: AlgShowerSRList.cxx,v 1.24 2006/12/01 18:50:24 rhatcher Exp $
00003 //
00004 // AlgShowerSRList.cxx
00005 //
00006 // Begin_Html<img src="../../pedestrians.gif" align=center>
00007 // <a href="../source_warning.html">Warning for beginners</a>.<br> 
00008 //
00009 // AlgShowerSRList is a concrete ShowerSRList Algorithm class.
00010 //
00011 // Author:  R. Lee 2001.02.21
00012 //
00013 // Also see <a href="../../root_crib/index.html">The ROOT Crib</a> and 
00014 // <a href="../CandShowerSR.html"> CandShowerSR 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 "CandShowerSR/AlgShowerSRList.h"
00026 #include "RecoBase/CandShower.h"
00027 #include "RecoBase/CandShower.h"
00028 #include "RecoBase/CandShowerHandle.h"
00029 #include "RecoBase/CandShowerList.h"
00030 #include "RecoBase/CandShowerListHandle.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/CandClusterListHandle.h"
00038 #include "RecoBase/CandSliceHandle.h"
00039 #include "RecoBase/CandSliceListHandle.h"
00040 #include "UgliGeometry/UgliGeomHandle.h"
00041 #include "Validity/VldContext.h"
00042 #include <string.h>
00043 
00044 ClassImp(AlgShowerSRList)
00045 
00046 
00047 //______________________________________________________________________
00048 CVSID("$Id: AlgShowerSRList.cxx,v 1.24 2006/12/01 18:50:24 rhatcher Exp $");
00049 
00050 //______________________________________________________________________
00051 AlgShowerSRList::AlgShowerSRList()
00052 {
00053 }
00054 
00055 //______________________________________________________________________
00056 AlgShowerSRList::~AlgShowerSRList()
00057 {
00058 }
00059 
00060 //______________________________________________________________________
00070 void AlgShowerSRList::RunAlg(AlgConfig &ac, CandHandle &ch,CandContext &cx)
00071 {
00072   assert(cx.GetDataIn());
00073 
00074   if (!(cx.GetDataIn()->InheritsFrom("TObjArray"))) {
00075     return;
00076   }
00077 
00078   Int_t diffViewPlaneMatch = ac.GetInt("DiffViewPlaneMatch");
00079   Double_t diffViewTimeMatch = ac.GetDouble("DiffViewTimeMatch");
00080   Int_t planeScale = ac.GetInt("PlaneScale");
00081   Int_t diffViewPlaneMatchShort = ac.GetInt("DiffViewPlaneMatchShort");
00082   Double_t diffViewPlaneCoverage = ac.GetDouble("DiffViewPlaneCoverage");
00083   Double_t diffViewPulseHeightCut = ac.GetDouble("DiffViewPulseHeightCut");
00084 
00085   const CandSliceListHandle *slicelist = 0;
00086   const CandClusterListHandle *clusterlist = 0;
00087   const TObjArray *cxin = dynamic_cast<const TObjArray *>(cx.GetDataIn());
00088   for (Int_t i=0; i<=cxin->GetLast(); i++) {
00089     TObject *tobj = cxin->At(i);
00090     if (tobj->InheritsFrom("CandSliceListHandle")) {
00091       slicelist = dynamic_cast<CandSliceListHandle*>(tobj);
00092     }
00093     if (tobj->InheritsFrom("CandClusterListHandle")) {
00094       clusterlist = dynamic_cast<CandClusterListHandle*>(tobj);
00095     }
00096   }
00097   if (!slicelist || !clusterlist) {
00098     MSG("error",Msg::kError) <<
00099       "CandSliceListHandle or CandClusterListHandle missing\n";
00100   }
00101 
00102   CandContext cxx(this,cx.GetMom());
00103   AlgFactory &af = AlgFactory::GetInstance();
00104 
00105   const char *showerAlgConfig = 0;
00106   ac.Get("ShowerAlgConfig",showerAlgConfig);
00107   AlgHandle ah = af.GetAlgHandle("AlgShowerSR",showerAlgConfig);
00108    
00109   const CandRecord *candrec = cx.GetCandRecord();
00110   assert(candrec);
00111   const VldContext *vldcptr = candrec->GetVldContext();
00112   assert(vldcptr);
00113   VldContext vldc = *vldcptr;
00114 
00115   UgliGeomHandle ugh(vldc);
00116 
00117   // Iterate over slices
00118   CandSliceHandleItr sliceItr(slicelist->GetDaughterIterator());
00119   while (CandSliceHandle *slice = sliceItr()) {
00120 
00121     // Select Clusters within this slice
00122     CandClusterHandleItr clusterItr(clusterlist->GetDaughterIterator());
00123     CandClusterHandleKeyFunc *clusterKf = clusterItr.CreateKeyFunc();
00124     clusterKf->SetFun(CandClusterHandle::KeyFromSlice);
00125     clusterItr.GetSet()->AdoptSortKeyFunc(clusterKf);
00126     clusterKf = 0;
00127     clusterItr.GetSet()->Slice(static_cast<Int_t>(slice->GetUidInt()));
00128 
00129     // construct arrays for U, V clusters and showers
00130     TObjArray showerlist;
00131     TObjArray uclusterlist;
00132     TObjArray vclusterlist;
00133 
00134     // iterate over clusters, and place in appropriate array based on view
00135     while (CandClusterHandle *cluster = clusterItr()) {
00136       if (*cluster->GetCandSlice()==*slice) {
00137         switch (cluster->GetPlaneView()) {
00138         case PlaneView::kU:
00139           uclusterlist.Add(cluster);
00140           break;
00141         case PlaneView::kV:
00142           vclusterlist.Add(cluster);
00143           break;
00144         default:
00145           break;
00146         }
00147       }
00148     }
00149 
00150     // now do a nested loop over U and V view clusters, looking for matches
00151     for (Int_t iu=0; iu<=uclusterlist.GetLast(); iu++) {
00152       CandClusterHandle *ucluster =
00153         dynamic_cast<CandClusterHandle*>(uclusterlist.At(iu));
00154       Float_t highesttotE=0;
00155       Int_t bestiv=0;
00156       Bool_t found=false;
00157       for (Int_t iv=0; iv<=vclusterlist.GetLast(); iv++) {
00158         CandClusterHandle *vcluster =
00159           dynamic_cast<CandClusterHandle*>(vclusterlist.At(iv));
00160        
00161         // determine fractional plane overlap, defined to be the ratio of the number of
00162         // overlapping planes to the full exent in U and V
00163         Int_t ubeg = ucluster->GetBegPlane();
00164         Int_t uend = ucluster->GetEndPlane();
00165         Int_t vbeg = vcluster->GetBegPlane();
00166         Int_t vend = vcluster->GetEndPlane();
00167         Int_t idir=1;
00168         if (ubeg>uend) idir=-1;
00169         Int_t dplanecover = idir*(min(uend*idir,vend*idir)-max(ubeg*idir,vbeg*idir));
00170         Int_t mincoverage = min(abs(1+uend-ubeg),abs(1+vend-vbeg));
00171         Float_t dplanecoverfrac = (Float_t)(dplanecover)/(Float_t)(mincoverage); 
00172 
00173         // change to rough energy unit - 330 PEs/GeV (WHY?)
00174         Float_t totE=max(ucluster->GetCharge()+vcluster->GetCharge(),1.)/330.;
00175         Float_t dE = fabs(ucluster->GetCharge()-vcluster->GetCharge())/330.;   
00176         Float_t dpulseheight = dE/sqrt(totE);
00177 
00178         // check of rshower clusters, as defined by parameter planeScale
00179         Int_t ishort=0; // 0 means both are short, 1 means one is short, 2 means none are short
00180         if (abs(vbeg-vend)>planeScale && abs(ubeg-uend)>planeScale) {
00181           ishort = 2;
00182         } else
00183         if (abs(vbeg-vend)>planeScale || abs(ubeg-uend)>planeScale) {
00184           ishort = 1;
00185         } else {
00186           ishort = 0;
00187         }
00188 
00189         // begin plane match criteria is:
00190         // if neither is short and difference between beginning plane 
00191         //in each view is less than diffViewPlaneMatch
00192         // or one or more is shower and difference in start 
00193         //plane is less than diffViewPlaneMatchShort
00194         Bool_t passplanematch(0);
00195         if ((ishort==2 && abs(ubeg-vbeg)<=diffViewPlaneMatch) || (ishort<2 && abs(ubeg-vbeg)<=diffViewPlaneMatchShort)) passplanematch=1;
00196 
00197         // time match criteria that begin times be within diffViewTimeMatch
00198         Bool_t passtimematch(0);
00199         if (fabs(ucluster->GetBegTime()-vcluster->GetBegTime())<diffViewTimeMatch) passtimematch=1;
00200 
00201         // plane coverage match criteria is that fractional plane coverage 
00202         // is greater than diffViewPlaneCoverage
00203         Bool_t passplanecoverage(1);
00204         if (ishort==2) {
00205           if (dplanecover<0 || dplanecoverfrac<diffViewPlaneCoverage) passplanecoverage=0;
00206         }
00207 
00208         // pulse heights must match to within diffViewPulseHeightCut
00209         Bool_t passpulseheight(1);
00210         if (ishort>=1) {
00211           if (dpulseheight>diffViewPulseHeightCut) passpulseheight=0;
00212         }
00213 
00214         // if all criteria are met, then....
00215         if (passplanematch && passtimematch && passplanecoverage && passpulseheight) {
00216 
00217           // see if this is the best cluster u match for this v cluster
00218           // recalculate all quantities for this v cluster, and all other u clusters
00219           // if this is not the best u cluster for the v cluster, don't generate shower
00220           Bool_t foundbetteriu=false;
00221           for (Int_t iu2=0; iu2<=uclusterlist.GetLast(); iu2++) {
00222             CandClusterHandle *ucluster2 =
00223               dynamic_cast<CandClusterHandle*>(uclusterlist.At(iu2));
00224             Int_t ubeg2 = ucluster2->GetBegPlane();
00225             Int_t uend2 = ucluster2->GetEndPlane();
00226             Int_t idir2=1;
00227             if (ubeg2>uend2) idir2=-1;
00228             Int_t dplanecover2 = idir2*(min(uend2*idir,vend*idir)-max(ubeg2*idir,vbeg*idir));
00229             Int_t mincoverage2 = min(abs(1+uend2-ubeg2),abs(1+vend-vbeg));
00230             Float_t dplanecoverfrac2 = (Float_t)(dplanecover2)/(Float_t)(mincoverage2); 
00231 
00232             // change to rough energy unit - 330 PEs/GeV            
00233             Float_t totE2=max(ucluster2->GetCharge()+vcluster->GetCharge(),1.)/330.;
00234             Float_t dE2 = fabs(ucluster2->GetCharge()-vcluster->GetCharge())/330.;   
00235             Float_t dpulseheight2 = dE2/sqrt(totE2);
00236             
00237             Int_t ishort2=0; // 0 means both are short, 1 means one is short, 2 means none are short
00238             if (abs(vbeg-vend)>planeScale && abs(ubeg2-uend2)>planeScale) {
00239               ishort2 = 2;
00240             } else
00241               if (abs(vbeg-vend)>planeScale || abs(ubeg2-uend2)>planeScale) {
00242                 ishort2 = 1;
00243               } else {
00244                 ishort2 = 0;
00245               }
00246 
00247             Bool_t passplanematch2(0);
00248             if ((ishort2==2 && abs(ubeg2-vbeg)<=diffViewPlaneMatch) || 
00249                 (ishort2<2 && abs(ubeg2-vbeg)<=diffViewPlaneMatchShort)) 
00250               passplanematch2=1;
00251 
00252             Bool_t passtimematch2(0);
00253             if (fabs(ucluster2->GetBegTime()-vcluster->GetBegTime())<
00254                 diffViewTimeMatch) passtimematch2=1;
00255 
00256             Bool_t passplanecoverage2(1);
00257             if (ishort2==2) {
00258               if (dplanecover2<0 || 
00259                   dplanecoverfrac2<diffViewPlaneCoverage) passplanecoverage2=0;
00260             }
00261 
00262             Bool_t passpulseheight2(1);
00263             if (ishort2>=1) {
00264               if (dpulseheight2>diffViewPulseHeightCut) passpulseheight2=0;
00265             }
00266 
00267             // check to see whether this UV pair meets requirements, and is higher in total
00268             // energy match found so r.
00269             if (passplanematch2 && 
00270                 passtimematch2 && 
00271                 passplanecoverage2 && 
00272                 passpulseheight2 && 
00273                 totE2>totE) {
00274               foundbetteriu=true;
00275             }
00276           }
00277           // if the U cluster is indeed the best match to this V cluster
00278           // this is a possible best match.  Store energy and continue loop over V cluster
00279           if(!foundbetteriu){
00280             found=true;
00281             if(totE>highesttotE){
00282               highesttotE=totE;
00283               bestiv=iv;
00284             }
00285           }
00286         }
00287       }
00288       // if match found, make shower from the two clusters.  
00289       // clusters should be used once only, with the best match partner
00290       if(found){
00291         CandClusterHandle *vcluster =  dynamic_cast<CandClusterHandle*>(vclusterlist.At(bestiv));
00292         TObjArray newshower;
00293         newshower.Add(ucluster);
00294         newshower.Add(vcluster);
00295         cxx.SetDataIn(&newshower);
00296         CandShowerHandle showerhandle = CandShower::MakeCandidate(ah,cxx);
00297         showerhandle.SetCandSlice(slice);
00298         ch.AddDaughterLink(showerhandle);
00299       }
00300     }
00301   }
00302 }
00303 
00304 
00305 
00306 
00307 //______________________________________________________________________
00308 
00309 void AlgShowerSRList::Trace(const char * /* c */) const
00310 {
00311 }
00312 
00313