loon: AlgShowerSRList Class Reference

void AlgShowerSRList::RunAlg	(	AlgConfig &	ac,
		CandHandle &	ch,
		CandContext &	cx
	)			`[virtual]`

We begin with a nested loop over all CandSlices, and over all CandClusters within each CandSlice, generating list of all clusters in the U and V planes. Next, we execute a nested loop over U-view and V-view CandClusters, checking for matching 2D clusters. The following criteria are used: Beginning planes for the two 2D 'long' CandClusters (defined by PlaneScale) must be within DiffViewPlaneMatch., or DiffViewPlaneMatchShort if one or more clusters is short. Beginning times for the two 2D CandClusters must be within DiffViewTimeMatch. The fractional overlap of U and V clusters must exceed DiffViewPlaneCoverage, if both clusters are long. The total pulse heights of the U and V clusters must agree to within DiffViewPulseHeightCut if at least one of the clusters is long. (Note: This should probably be changed to a charge ratio) If a U/V cluster set is found which satisfies these requirements, a second loop over U clusters is performed, searching for alternative U/V pairs involving the original V cluster satisfying the requirements above, but higher in total energy. If no better alternative is found a CandShower is constructed. In this way, a given CandCluster will only be used once, in the shower with the highest total energy.

Implements AlgBase.

Reimplemented in AlgShowerSSList.

Definition at line 70 of file AlgShowerSRList.cxx.

References CandHandle::AddDaughterLink(), Registry::Get(), CandContext::GetCandRecord(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), Registry::GetDouble(), AlgFactory::GetInstance(), Registry::GetInt(), CandContext::GetMom(), RecMinos::GetVldContext(), Msg::kError, CandClusterHandle::KeyFromSlice(), PlaneView::kU, PlaneView::kV, CandShower::MakeCandidate(), max, min, MSG, and CandRecoHandle::SetCandSlice().

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 }


Public Member Functions
	AlgShowerSRList ()
virtual	~AlgShowerSRList ()
virtual void	RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx)
virtual void	Trace (const char *c) const

AlgShowerSRList Class Reference

Public Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation