Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

AlgEventSSList Class Reference

#include <AlgEventSSList.h>

Inheritance diagram for AlgEventSSList:

AlgBase List of all members.

Public Member Functions

 AlgEventSSList ()
virtual ~AlgEventSSList ()
virtual void RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx)
void BuildEventFromUnassoc (AlgConfig &ac, CandHandle &ch, CandContext &cx, CandSliceHandle *slice)
virtual void Trace (const char *c) const
void SetPrimaryShowers (CandHandle &ch, AlgConfig &ac)
void ReConstructShowers (AlgConfig &ac, CandHandle &ch, CandContext &cx)
void FillRecoList (CandSliceHandle *, CandShowerHandleItr *, CandTrackHandleItr *, TObjArray &)
void AddObjectToEvent (CandRecoHandle *reco, TObjArray *, TObjArray *, Bool_t)
void AddStripToEvent (CandEventHandle *ev, CandShowerHandle *shw, CandSubShowerSRListHandle *subshwlist, CandStripHandle *closeststrip)
void CreatePrimaryShower (AlgConfig &ac, CandContext &cxx, CandEventHandle *closestevent, CandShowerListHandle *showerlist, CandSubShowerSRListHandle *subshowerlist, CandStripHandle *closeststrip)
void FindUnassociated (CandSliceHandleItr &sliceItr, CandEventHandleItr *, TObjArray &unassociated)
void ReFillDist2Map (AlgConfig &ac, CandStripHandle *closeststrip, CandEventHandle *closestevent, TObjArray &unassociated, std::vector< std::map< const CandEventHandle *, Double_t > > &dist2map)
void FillDist2Map (AlgConfig &ac, TObjArray &unassociated, CandHandle &ch, std::vector< std::map< const CandEventHandle *, Double_t > > &dist2map)

Constructor & Destructor Documentation

AlgEventSSList::AlgEventSSList  ) 
 

Definition at line 62 of file AlgEventSSList.cxx.

00063 {
00064 }

AlgEventSSList::~AlgEventSSList  )  [virtual]
 

Definition at line 67 of file AlgEventSSList.cxx.

00068 {
00069 }


Member Function Documentation

void AlgEventSSList::AddObjectToEvent CandRecoHandle reco,
TObjArray *  ,
TObjArray *  ,
Bool_t 
 

Definition at line 1821 of file AlgEventSSList.cxx.

References MSG.

01822                                                                            {
01823   if (!merge) {
01824     // if this is the first match we have found for this object, 
01825     //add object to list
01826     // and point prevlist to this object list
01827     objectlist->Add(reco);
01828     MSG("EventSS",Msg::kVerbose) << "    adding to list\n";
01829   }
01830   else {                  
01831     MSG("EventSS",Msg::kVerbose) << "    combining lists\n";
01832     // this object matches more than one event - merge their member
01833     // objects to list prevlist, and add the new object
01834     for (Int_t ireco2=0; ireco2<=objectlist->GetLast();ireco2++)
01835       prevlist->Add(objectlist->At(ireco2));
01836     prevlist->Add(reco);
01837     objectlist->Clear();
01838   }
01839 }

void AlgEventSSList::AddStripToEvent CandEventHandle closestevent,
CandShowerHandle shower,
CandSubShowerSRListHandle subshowerlist,
CandStripHandle closeststrip
 

Add a strip to the event daughter list, the primary shower, and place the modified cluster on the cluster list

Definition at line 966 of file AlgEventSSList.cxx.

References CandHandle::AddDaughterLink(), CandShowerSRHandle::AddSubShower(), CandSubShowerSRHandle::DupHandle(), CandSubShowerSRHandle::GetAvgDev(), CandRecoHandle::GetCandSlice(), CandRecoHandle::GetCharge(), CandHandle::GetDaughterIterator(), CandSubShowerSRHandle::GetEnergy(), CandShowerSRHandle::GetLastSubShower(), Calibrator::GetMIP(), CandStripHandle::GetPlane(), CandStripHandle::GetPlaneView(), CandSubShowerSRHandle::GetPlaneView(), CandSubShowerSRHandle::GetSlope(), CandShowerSRHandle::GetSubShower(), CandStripHandle::GetTPos(), CandHandle::GetUidInt(), CandRecoHandle::GetVtxU(), CandRecoHandle::GetVtxV(), CandRecoHandle::GetVtxZ(), CandStripHandle::GetZPos(), Calibrator::Instance(), CandHandle::IsEquivalent(), MSG, CandHandle::RemoveDaughter(), CandShowerSRHandle::RemoveSubShower(), CandRecoHandle::SetCandSlice(), and CandSubShowerSRHandle::SetEnergy().

Referenced by RunAlg().

00970 {
00971   MSG("EventSS",Msg::kVerbose) << "In AddStripToEvent:" << endl;
00972   closestevent->AddDaughterLink(*closeststrip);
00973   shower->AddDaughterLink(*closeststrip);
00974   if(subshowerlist && shower->InheritsFrom("CandShowerSRHandle")) {
00975     Calibrator& calibrator=Calibrator::Instance();
00976     Int_t foundsubshower=0;
00977     Int_t closestSubShower=-1;
00978     Double_t BestStripDistance = 99999;
00979     Double_t StripDistance = 99999;
00980     CandShowerSRHandle *showerSR = 
00981       dynamic_cast<CandShowerSRHandle*> (shower);
00982     MSG("EventSS",Msg::kVerbose) << "Shower UID = " 
00983                                  << shower->GetUidInt() << endl;
00984 
00985     for (Int_t isubshower=0;
00986          isubshower<showerSR->GetLastSubShower()+1;
00987          isubshower++) {
00988       CandSubShowerSRHandle *ssl = const_cast<CandSubShowerSRHandle*>
00989         (showerSR->GetSubShower(isubshower));
00990       //add strip to cluster of the same view
00991       if (ssl->GetPlaneView()==closeststrip->GetPlaneView()) {
00992         //decide which is most appropriate subshower to add strip to
00993         //based on distance to/spread of subshower:
00994         Double_t tposVtx = 0;
00995         if(ssl->GetPlaneView()==PlaneView::kU || 
00996            ssl->GetPlaneView()==PlaneView::kX) 
00997           tposVtx = ssl->GetVtxU();
00998         else if(ssl->GetPlaneView()==PlaneView::kV || 
00999                 ssl->GetPlaneView()==PlaneView::kY)
01000           tposVtx = ssl->GetVtxV();
01001         StripDistance = 99999;
01002         if(ssl->GetAvgDev()>0) {
01003           StripDistance = 
01004             (TMath::Abs(closeststrip->GetTPos() - 
01005                         (tposVtx + ssl->GetSlope() * 
01006                          (closeststrip->GetZPos() - 
01007                           ssl->GetVtxZ()))) * 
01008              TMath::Cos(TMath::ATan(ssl->GetSlope()))) * 
01009             calibrator.GetMIP(closeststrip->
01010                               GetCharge(CalDigitType::kSigCorr)) / 
01011             ssl->GetAvgDev();
01012         }
01013         if(StripDistance>9999){ //bad slope
01014           StripDistance = 9999 + 
01015             TMath::Sqrt(TMath::Power(TMath::Abs(closeststrip->GetZPos() - 
01016                                                 ssl->GetVtxZ()),2) +
01017                         TMath::Power(TMath::Abs(closeststrip->GetTPos() - 
01018                                                 tposVtx),2));
01019         }
01020         if(StripDistance<BestStripDistance || foundsubshower==0){
01021           foundsubshower=1;
01022           BestStripDistance = StripDistance;
01023           closestSubShower = isubshower;
01024         }
01025       }
01026     }
01027     if(foundsubshower) {
01028       MSG("EventSS",Msg::kVerbose) << "adding strip to subshower "
01029                                    << closeststrip->GetPlane() 
01030                                    << endl;
01031       CandSubShowerSRHandle *ssl = const_cast<CandSubShowerSRHandle*>
01032         (showerSR->GetSubShower(closestSubShower));
01033       //extract mip -> gev conversion:
01034       Double_t GeVperMIP = 0;
01035       if(ssl->GetEnergy()!=0) GeVperMIP = ssl->GetEnergy()/
01036         Double_t(calibrator.GetMIP(ssl->GetCharge(CalStripType::kSigCorr)));
01037       TIter csshItr(subshowerlist->GetDaughterIterator());
01038       while (CandSubShowerSRHandle *cssh = 
01039              dynamic_cast<CandSubShowerSRHandle*>(csshItr())) { 
01040         if(cssh->IsEquivalent(ssl) && ssl->GetUidInt()==cssh->GetUidInt()) {
01041           MSG("EventSS",Msg::kVerbose) << "SubShower (old,new) UID " 
01042                                        << ssl->GetUidInt()
01043                                        << " " << cssh->GetUidInt() << endl; 
01044           showerSR->RemoveSubShower(ssl);
01045           CandSubShowerSRHandle *newss = cssh->DupHandle();
01046           newss->SetCandSlice(cssh->GetCandSlice());
01047           newss->AddDaughterLink(*closeststrip);
01048           newss->SetEnergy(GeVperMIP*Double_t(calibrator.GetMIP(closeststrip->
01049                  GetCharge(CalDigitType::kSigCorr))) + newss->GetEnergy());
01050           subshowerlist->RemoveDaughter(cssh);
01051           subshowerlist->AddDaughterLink(*newss);
01052           showerSR->AddSubShower(newss);
01053           delete newss;
01054           break;
01055         }
01056       }
01057     }
01058     else {
01059       MSG("EventSS",Msg::kWarning) << "Strip added to shower but not to subshower!"
01060                                    << " StripDistance = " << StripDistance
01061                                    << endl;
01062       MSG("EventSS",Msg::kWarning) << "Number of subshowers in shower = " 
01063                                    << showerSR->GetLastSubShower()+1
01064                                    << endl;
01065       for(int i=0;i<showerSR->GetLastSubShower()+1;i++){
01066         CandSubShowerSRHandle *ssl = const_cast<CandSubShowerSRHandle*>
01067           (showerSR->GetSubShower(i));
01068         MSG("EventSS",Msg::kWarning) << "SS " << i << " UID = " 
01069                                      << ssl->GetUidInt() << endl;
01070       }
01071     }
01072   }
01073 }

void AlgEventSSList::BuildEventFromUnassoc AlgConfig ac,
CandHandle ch,
CandContext cx,
CandSliceHandle slice
 

this method is called for slices in which no shower or track was previously found. We see whether there is a 3D cluster in this slice which could be an event

Definition at line 1389 of file AlgEventSSList.cxx.

References abs(), CandHandle::AddDaughterLink(), CandHandle::FindDaughter(), Registry::Get(), AlgFactory::GetAlgHandle(), CandStripHandle::GetBegTime(), CandHandle::GetCandRecord(), CandContext::GetCandRecord(), CandStripHandle::GetCharge(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), Registry::GetDouble(), AlgFactory::GetInstance(), Registry::GetInt(), CandContext::GetMom(), CandStripHandle::GetPlane(), CandStripHandle::GetPlaneView(), CandStripHandle::GetTPos(), RecMinos::GetVldContext(), CandHandle::GetVldContext(), CandRecoHandle::GetVtxT(), CandRecoHandle::GetVtxU(), CandRecoHandle::GetVtxV(), CandRecoHandle::GetVtxZ(), CandShowerSRHandle::IsUnphysical(), CandEvent::MakeCandidate(), CandShowerSR::MakeCandidate(), CandSubShowerSR::MakeCandidate(), MSG, CandHandle::SetCandRecord(), CandEventHandle::SetCandSlice(), CandRecoHandle::SetCandSlice(), CandSubShowerSRHandle::SetClusterID(), and CandEventHandle::SetPrimaryShower().

Referenced by RunAlg().

01393 {
01394   assert(cx.GetDataIn());
01395   MSG("EventSS",Msg::kDebug) <<" Building Event - no showers or tracks " 
01396                              << endl;
01397   if (!(cx.GetDataIn()->InheritsFrom("TObjArray"))) {
01398     MSG("EventSS",Msg::kDebug)<<" dataIn inherits from TObjArray " << endl;
01399     return;
01400   }
01401 
01402   Double_t pShwXTalkFracCut = ac.GetDouble("ShwXTalkFracCut");
01403   Double_t pShwXTalkDef = ac.GetDouble("ShwXTalkDef");
01404 
01405   const CandSliceListHandle *slicelist = 0;
01406   CandShowerListHandle *showerlist = 0;
01407   CandSubShowerSRListHandle *subshowerlist = 0;
01408   const TObjArray *cxin = dynamic_cast<const TObjArray *>(cx.GetDataIn());
01409   for (Int_t i=0; i<=cxin->GetLast(); i++) {
01410     TObject *tobj = cxin->At(i);
01411     if (tobj->InheritsFrom("CandSliceListHandle")) {
01412       slicelist = dynamic_cast<CandSliceListHandle*>(tobj);
01413     }
01414     // Retrieve persistant modifiable CandShowerList
01415     if (tobj->InheritsFrom("CandShowerListHandle")) {
01416       showerlist = dynamic_cast<CandShowerListHandle*>(tobj);
01417     }
01418     if (tobj->InheritsFrom("CandSubShowerSRListHandle")) {
01419       subshowerlist = dynamic_cast<CandSubShowerSRListHandle*>(tobj);
01420     }
01421   }
01422   
01423   if (!showerlist || !slicelist || !subshowerlist) {
01424     MSG("EventSS",Msg::kDebug)
01425       << "Missing either slicelist, showerlist or subshowerlist" << endl;
01426     return;
01427   }
01428   Int_t detectorType = slicelist->GetVldContext()->GetDetector();
01429   
01430   CandContext cxx(this,cx.GetMom());
01431   AlgFactory &af = AlgFactory::GetInstance();
01432   
01433   const char *pEventAlgConfig = 0;
01434   ac.Get("EventAlgConfig",pEventAlgConfig);
01435   AlgHandle ah = af.GetAlgHandle("AlgEventSR",pEventAlgConfig);
01436   
01437   Double_t pHitAssocDPlane = ac.GetInt("HitAssocDPlane");
01438   Double_t pHitAssocDTime0 = ac.GetDouble("HitAssocDTime0");
01439   Double_t pHitAssocDTime1 = ac.GetDouble("HitAssocDTime1");
01440   Double_t pHitAssocMaxDist2 = ac.GetDouble("HitAssocMaxDist2");
01441   Double_t pHitAssocZParm = ac.GetDouble("HitAssocZParm");
01442   Double_t pHitAssocPParm = ac.GetDouble("HitAssocPParm");
01443   Double_t pHitAssocTParm = ac.GetDouble("HitAssocTParm");
01444   Double_t pShwShwDtpos2 = ac.GetDouble("ShwShwDtpos2");
01445   Double_t pShwShwDz = ac.GetDouble("ShwShwDz");
01446   Double_t pShwShwDt = ac.GetDouble("ShwShwDt");
01447   Double_t pMinShwEFract = ac.GetDouble("MinShwEFract");
01448 
01449 
01450   const CandRecord *candrec = cx.GetCandRecord();
01451   assert(candrec);
01452   const VldContext *vldcptr = candrec->GetVldContext();
01453   assert(vldcptr);
01454   VldContext vldc = *vldcptr;
01455   UgliGeomHandle ugh(vldc);
01456 
01457   // Iterate over strips , finding sets which are 'in time', and within
01458   // a given distance of the nearest neighbor in the set, in both views.
01459   // Each of these sets is then used to construct a shower, which becomes
01460   // the primary shower of a new event.  We continue until an iteration
01461   // in which no new event is created.
01462   
01463   // fill array of unique unassociated strips
01464   TObjArray unassociated;
01465   TIter stripItr(slice->GetDaughterIterator());
01466   while (CandStripHandle *strip =
01467          dynamic_cast<CandStripHandle*>(stripItr())) {
01468     Bool_t found=false;
01469     for (Int_t i=0; i<=unassociated.GetLast() && !found; i++) {
01470       CandStripHandle *strip0 =
01471         dynamic_cast<CandStripHandle*>(unassociated.At(i));
01472       if (*strip == *strip0) {
01473         found = true;
01474       }
01475     }
01476     if (!found) {
01477       if(!(detectorType==Detector::kNear && strip->GetPlane()>120)){
01478         unassociated.Add(strip);
01479       }
01480     }
01481   }
01482   
01483   // loop until no new matches are found in the loop over unassoc hits.
01484   Bool_t found=false;
01485   while (unassociated.GetLast()>0 && !found) {
01486     MSG("EventSS",Msg::kDebug)<< " # unassoc " 
01487                               << unassociated.GetLast()+1 << endl; 
01488     Bool_t firstu=true;
01489     Bool_t firstv=true;
01490     TObjArray neweventu;
01491     TObjArray neweventv;
01492     Bool_t addedstrip = true;
01493     Float_t totcharge=0;
01494     while (addedstrip) {
01495       addedstrip=false;      
01496       // loop over unassoc strips.
01497       for (Int_t i=0; i<=unassociated.GetLast(); i++) {
01498         CandStripHandle *strip =
01499           dynamic_cast<CandStripHandle*>(unassociated.At(i));
01500         switch (strip->GetPlaneView()) {
01501         case PlaneView::kU:
01502           // strip is in U plane
01503           // if this is the first strip, add to newevent strip list
01504           if (firstu) {
01505             // if we don't have v strips yet, just add this strip to u list
01506             // seed the event
01507             if (firstv) {
01508               neweventu.Add(strip);
01509               MSG("EventSS",Msg::kDebug)<< " add  first  " << endl;
01510               totcharge +=strip->GetCharge();
01511               addedstrip = true;
01512               firstu=false;
01513             }
01514             // if we already have one or more v strips check that they are
01515             // in time, and  agree in Z position.
01516             else {
01517               Bool_t found2=false;
01518               for (Int_t j=0; j<=neweventv.GetLast() && !found2;
01519                    j++) {
01520                 CandStripHandle *estrip =
01521                   dynamic_cast<CandStripHandle*>(neweventv.At(j));
01522                 Double_t dtime =
01523                   strip->GetBegTime()-estrip->GetBegTime(); 
01524                 MSG("EventSS",Msg::kVerbose) << " dtime " << dtime*1e9 
01525                                              << " " <<pHitAssocDTime0*1e9 
01526                                              << "/" << pHitAssocDTime1*1e9 
01527                                              << " plane diff " 
01528                                              <<  abs(strip->GetPlane() - 
01529                                                      estrip->GetPlane()) 
01530                                              << "/" << pHitAssocDPlane << endl;
01531                 if (dtime>pHitAssocDTime0 &&
01532                     dtime<pHitAssocDTime1 &&
01533                     abs(strip->GetPlane()-estrip->GetPlane()) <=
01534                     pHitAssocDPlane) {
01535                   MSG("EventSS",Msg::kVerbose)<< " add  " << endl;
01536                   // match is found, add this strip to array
01537                   neweventu.Add(strip);
01538                   totcharge +=strip->GetCharge();
01539                   addedstrip = true;
01540                   firstu=0;
01541                   found2=true;
01542                 } // if match found
01543               }  // for loop over V strips in new event
01544             } // end not first V strip
01545           } // end first U strip
01546           else {  
01547             // if we already have u strips, make sure that the one we
01548             // add is adjaccent, using distance criteria based
01549             // on plane, time, and transverse postion separation
01550             Bool_t found2=false;
01551             for (Int_t j=0; j<=neweventu.GetLast() && !found2; j++) {
01552               CandStripHandle *estrip =
01553                 dynamic_cast<CandStripHandle*>(neweventu.At(j));
01554               Double_t dtime =
01555                 strip->GetBegTime()-estrip->GetBegTime(); 
01556               MSG("EventSS",Msg::kVerbose) << " dtime " << dtime*1e9 
01557                                            << " " <<pHitAssocDTime0*1e9 
01558                                            << "/" << pHitAssocDTime1*1e9 
01559                                            << " plane diff " 
01560                                            << abs(strip->GetPlane() - 
01561                                                   estrip->GetPlane())
01562                                            << " dtpos " << (strip->GetTPos() - 
01563                                                             estrip->GetTPos())
01564                                            << endl;
01565               if (dtime>pHitAssocDTime0 &&
01566                   dtime<pHitAssocDTime1) {
01567                 Double_t dplane =
01568                   (Double_t)(strip->GetPlane()-estrip->GetPlane());
01569                 Double_t dtpos = strip->GetTPos()-estrip->GetTPos();
01570                 Double_t dist2 = (pHitAssocZParm*dplane*dplane) +
01571                   (pHitAssocPParm*dtpos*dtpos) +
01572                   (pHitAssocTParm*dtime*dtime);
01573                 // if dist2 2 parameter meets match criteria
01574                 // add strip to array
01575                 if (dist2<pHitAssocMaxDist2*4) {
01576                   MSG("EventSS",Msg::kVerbose)<< " add " << endl;
01577                   neweventu.Add(strip);
01578                   totcharge +=strip->GetCharge();
01579                   addedstrip = true;
01580                   firstu=0;
01581                   found2=true;
01582                 } // if positions match
01583               } // if in time
01584             } // end loop over event u strips
01585           } // end if already have u strip
01586           break;
01587           
01588         case PlaneView::kV:
01589           
01590           // if this is the first strip, add to newevent strip list
01591           if (firstv) {
01592             // strip is in V plane
01593             // if this is the first strip, add to newevent strip list
01594             if (firstu) {
01595               neweventv.Add(strip);
01596               totcharge +=strip->GetCharge();
01597               MSG("EventSS",Msg::kDebug)<< " add  first" << endl;
01598               addedstrip = true;
01599               firstv = 0;
01600             }
01601             else {
01602               // if we already have one or more u strips check that they are
01603               // in time, and require agreement in Z position.
01604               Bool_t found2=false;
01605               for (Int_t j=0; j<=neweventu.GetLast() && !found2;
01606                    j++) {
01607                 CandStripHandle *estrip =
01608                   dynamic_cast<CandStripHandle*>(neweventu.At(j));
01609                 Double_t dtime =
01610                   strip->GetBegTime()-estrip->GetBegTime();
01611                 MSG("EventSS",Msg::kVerbose) << " dtime " << dtime*1e9 
01612                                              << " " << pHitAssocDTime0*1e9 
01613                                              << "/" << pHitAssocDTime1*1e9 
01614                                              << " plane diff " 
01615                                              <<  abs(strip->GetPlane() - 
01616                                                      estrip->GetPlane()) 
01617                                              << "/" << pHitAssocDPlane << endl;
01618                 if (dtime>pHitAssocDTime0 &&
01619                     dtime<pHitAssocDTime1 &&
01620                     abs(strip->GetPlane()-estrip->GetPlane()) <=
01621                     pHitAssocDPlane) {
01622                   // if dist2 2 parameter meets match criteria, 
01623                   //add strip to array
01624                   MSG("EventSS",Msg::kVerbose) << " add  " << endl;
01625                   neweventv.Add(strip);
01626                   totcharge +=strip->GetCharge();
01627                   addedstrip = true;
01628                   firstv = false;
01629                   found2 = true;
01630                 }
01631               } // end loop over u strips in event
01632             } // end if already have u strips
01633           } // end if this is the first v strip
01634           else {
01635             // if we already have v strips, make sure that the one we
01636             // add is adjaccent, using same distance criteria base
01637             // on plane, time, and transverse postion separation
01638             Bool_t found2=false;
01639             for (Int_t j=0; j<=neweventv.GetLast() && !found2; j++) {
01640               CandStripHandle *estrip =
01641                 dynamic_cast<CandStripHandle*>(neweventv.At(j));
01642               Double_t dtime = strip->GetBegTime() -
01643                 estrip->GetBegTime();              
01644               MSG("EventSS",Msg::kVerbose)<< " dtime " << dtime*1e9 << " " 
01645                                           <<  pHitAssocDTime0*1e9 << "/" 
01646                                           << pHitAssocDTime1*1e9 
01647                                           << " plane diff " 
01648                                           <<  abs(strip->GetPlane() - 
01649                                                   estrip->GetPlane())
01650                                           << " dtpos " << (strip->GetTPos() - 
01651                                                            estrip->GetTPos())
01652                                           << endl;
01653               if (dtime>pHitAssocDTime0 &&
01654                   dtime<pHitAssocDTime1) {
01655                 Double_t dplane =
01656                   (Double_t)(strip->GetPlane()-estrip->GetPlane());
01657                 Double_t dtpos = strip->GetTPos()-estrip->GetTPos();
01658                 Double_t dist2 = (pHitAssocZParm*dplane*dplane) +
01659                   (pHitAssocPParm*dtpos*dtpos) +
01660                   (pHitAssocTParm*dtime*dtime);
01661                 // if dist2 2 parameter meets match criteria, 
01662                 // add strip to array
01663                 if (dist2<pHitAssocMaxDist2*4) {
01664                   MSG("EventSS",Msg::kVerbose) << " add  " << endl;
01665                   neweventv.Add(strip);
01666                   totcharge +=strip->GetCharge();
01667                   addedstrip = true;
01668                   firstv = false;
01669                   found2 = true;
01670                 } // end if add strip
01671               } //end if in time
01672             }  // end loop over v strips
01673           }
01674           break;
01675           
01676         default:
01677           break;
01678         }
01679       }
01680       
01681       // end of loop over hits - remove new 'associations' and compress
01682       for (Int_t i=0; i<=neweventu.GetLast(); i++) {
01683         CandStripHandle *estrip =
01684           dynamic_cast<CandStripHandle*>(neweventu.At(i));
01685         unassociated.Remove(estrip);
01686       }
01687       for (Int_t i=0; i<=neweventv.GetLast(); i++) {
01688         CandStripHandle *estrip =
01689           dynamic_cast<CandStripHandle*>(neweventv.At(i));
01690         unassociated.Remove(estrip);
01691       }
01692       unassociated.Compress();
01693     }
01694     
01695     // if we have hits in both views, build an event from them
01696     // if pulse height>10 pes
01697     MSG("EventSS",Msg::kDebug)<< " considering event formation  u/v: " 
01698                               << !firstu  << "/" << !firstv << " charge " 
01699                               << totcharge << endl;
01700     if (!firstu && !firstv && totcharge>20.0) {      
01701       found=true;
01702 
01703       //construct CandSubShowerSRs
01704       cxx.SetDataIn(&neweventu);
01705       AlgHandle subshowerah = af.GetAlgHandle("AlgSubShowerSR","default");
01706       CandSubShowerSRHandle usubshowerhandle =
01707         CandSubShowerSR::MakeCandidate(subshowerah,cxx);
01708       usubshowerhandle.SetCandSlice(slice);
01709       usubshowerhandle.SetClusterID(ClusterType::kUnknown);
01710       subshowerlist->AddDaughterLink(usubshowerhandle);
01711       cxx.SetDataIn(&neweventv);
01712       CandSubShowerSRHandle vsubshowerhandle =
01713         CandSubShowerSR::MakeCandidate(subshowerah,cxx);
01714       vsubshowerhandle.SetCandSlice(slice);
01715       vsubshowerhandle.SetClusterID(ClusterType::kUnknown);
01716       subshowerlist->AddDaughterLink(vsubshowerhandle);
01717       CandHandle *usubshower=subshowerlist->FindDaughter(&usubshowerhandle);
01718       CandHandle *vsubshower=subshowerlist->FindDaughter(&vsubshowerhandle);
01719       
01720       //build shower
01721       TObjArray newshower;
01722       newshower.Add(usubshower);
01723       newshower.Add(vsubshower);
01724       cxx.SetDataIn(&newshower);
01725       AlgHandle showerah = af.GetAlgHandle("AlgShowerSS","default");
01726       CandShowerSRHandle showerhandle =
01727         CandShowerSR::MakeCandidate(showerah,cxx);
01728       showerhandle.SetCandSlice(slice);
01729       showerhandle.SetCandRecord(slicelist->GetCandRecord());
01730 
01731       Bool_t newShowerOverlapsOld=false;
01732       CandShowerHandleItr showerItr(showerlist->GetDaughterIterator());      
01733       while (CandShowerHandle *shower =
01734              dynamic_cast<CandShowerHandle*>(showerItr())){
01735 
01736         Double_t dz = showerhandle.GetVtxZ()-shower->GetVtxZ();
01737         Double_t du = showerhandle.GetVtxU()-shower->GetVtxU();
01738         Double_t dv = showerhandle.GetVtxV()-shower->GetVtxV();
01739         Double_t dt = showerhandle.GetVtxT()-shower->GetVtxT();
01740 
01741         if( ( (du*du+dv*dv<pShwShwDtpos2 &&fabs(dz)<pShwShwDz) ||
01742               (du*du<pShwShwDtpos2/4. && dv*dv<pShwShwDtpos2*4. && 
01743                fabs(dz)<pShwShwDz/2.) ||
01744               (du*du<pShwShwDtpos2*4 && dv*dv<pShwShwDtpos2/4. && 
01745                fabs(dz)<pShwShwDz/2.) ) &&  
01746             fabs(dt)<pShwShwDt ) {
01747           newShowerOverlapsOld = true;
01748           MSG("EventSS",Msg::kDebug) 
01749             << " new shower overlaps with previous - don't make new event " 
01750             << endl;
01751           break;
01752         }
01753       }
01754         
01755       if(!newShowerOverlapsOld){
01756         if(!showerhandle.IsUnphysical(pShwXTalkFracCut,pShwXTalkDef)) {
01757           MSG("EventSS",Msg::kDebug)<< " Building Event! " << endl;
01758           showerlist->AddDaughterLink(showerhandle);
01759           TObjArray recolist;
01760           recolist.Add(&showerhandle);
01761           cxx.SetDataIn(&recolist);
01762           AlgHandle eventah = af.GetAlgHandle("AlgEventSR",pEventAlgConfig);
01763           CandEventHandle eventhandle =
01764             CandEvent::MakeCandidate(eventah,cxx);
01765           eventhandle.SetCandSlice(slice);
01766           eventhandle.SetPrimaryShower(pMinShwEFract,pShwShwDz);
01767           ch.AddDaughterLink(eventhandle);      
01768         }
01769       }
01770     }
01771   }
01772 }

void AlgEventSSList::CreatePrimaryShower AlgConfig ac,
CandContext cxx,
CandEventHandle closestevent,
CandShowerListHandle showerlist,
CandSubShowerSRListHandle subshowerlist,
CandStripHandle closeststrip
 

Definition at line 1077 of file AlgEventSSList.cxx.

References abs(), CandHandle::AddDaughterLink(), CandEventHandle::AddShower(), CandHandle::FindDaughter(), Registry::Get(), AlgFactory::GetAlgHandle(), CandRecoHandle::GetBegPlane(), CandEventHandle::GetCandSlice(), CandStripHandle::GetCharge(), CandHandle::GetDaughterIterator(), Registry::GetDouble(), AlgFactory::GetInstance(), CandStripHandle::GetPlane(), CandSubShowerSRHandle::GetPlaneView(), CandStripHandle::GetPlaneView(), CandEventHandle::GetPrimaryTrack(), CandShowerSR::MakeCandidate(), CandSubShowerSR::MakeCandidate(), MSG, CandRecoHandle::SetCandSlice(), CandSubShowerSRHandle::SetClusterID(), CandContext::SetDataIn(), and CandEventHandle::SetPrimaryShower().

Referenced by RunAlg().

01082 {
01083   // need to create a shower -- include strips from first plane in each
01084   // view plus this previously unassociated hit to generate a 3D object
01085   MSG("EventSS",Msg::kVerbose) << "In CreatePrimaryShower" << endl;
01086   // Get singleton instance of AlgFactory
01087   const char *pEventAlgConfig = 0;
01088   ac.Get("EventAlgConfig",pEventAlgConfig);
01089 
01090   //Double_t shwshwdz = ac.GetDouble("ShwShwDz");
01091   //Double_t minShwEFract = ac.GetDouble("MinShwEFract");
01092   Double_t pMinShwStripPE = ac.GetDouble("MinShwStripPE");
01093   Double_t pMaxNewShwLen = ac.GetDouble("MaxNewShwLen");
01094 
01095   AlgFactory &af = AlgFactory::GetInstance();
01096   
01097   CandTrackHandle *track = closestevent->GetPrimaryTrack();
01098   if (!showerlist) {
01099     MSG("EventSS",Msg::kError) <<"CandShowerListHandle missing\n" << endl;
01100     return;
01101   }
01102   
01103   closestevent->AddDaughterLink(*closeststrip);
01104   
01105   // create arrays of strips for each planeview to create
01106   // shower from
01107   
01108   MSG("EventSS",Msg::kVerbose) <<  "creating shower \n";
01109   TObjArray ustriparray;
01110   TObjArray vstriparray;
01111   Int_t uplane=track->GetBegPlane(PlaneView::kU);
01112   Int_t vplane=track->GetBegPlane(PlaneView::kV);
01113   // seed this shower with strips in most upstream track hits 
01114   // in each view, along with hits above a charge threshold and
01115   // within pMaxNewShwLen of  the track start
01116   TIter stripItr(track->GetDaughterIterator());
01117   while (CandStripHandle *tstrip =
01118          dynamic_cast<CandStripHandle*>(stripItr())) {
01119     if ((tstrip->GetPlane() == uplane) || 
01120         (tstrip->GetPlaneView()==PlaneView::kU &&
01121             tstrip->GetCharge()>pMinShwStripPE && 
01122             abs(tstrip->GetPlane()-uplane)<pMaxNewShwLen)) {
01123       ustriparray.Add(tstrip);
01124     }
01125     if ((tstrip->GetPlane() == vplane) || 
01126         (tstrip->GetPlaneView()==PlaneView::kV &&
01127             tstrip->GetCharge()>pMinShwStripPE && 
01128             abs(tstrip->GetPlane()-vplane)<pMaxNewShwLen)) {
01129       vstriparray.Add(tstrip);
01130     }
01131   }
01132   // add unassoc strip to appropriate array
01133   switch (closeststrip->GetPlaneView()) {
01134   case PlaneView::kU:
01135     ustriparray.Add(closeststrip);
01136     break;
01137   case PlaneView::kV:
01138     vstriparray.Add(closeststrip);
01139     break;
01140   default:
01141     break;
01142   }
01143   // build shower from U/V clusters or subshowers
01144   if(ustriparray.GetLast()>=0 && vstriparray.GetLast()>=0){
01145     if(showerlist->InheritsFrom("CandShowerSRListHandle") && subshowerlist){
01146       MSG("EventSS",Msg::kVerbose) << "U striparray has " 
01147                                    << ustriparray.GetLast()+1 << " entries" << endl;
01148       MSG("EventSS",Msg::kVerbose) << "V striparray has " 
01149                                    << vstriparray.GetLast()+1 << " entries" << endl;  
01150 
01151       // now construct CandSubShowers
01152       AlgHandle subshowerah = af.GetAlgHandle("AlgSubShowerSR","default");      
01153       cxx.SetDataIn(&ustriparray);
01154       CandSubShowerSRHandle usubshowerhandle =
01155         CandSubShowerSR::MakeCandidate(subshowerah,cxx);
01156       usubshowerhandle.SetCandSlice(closestevent->GetCandSlice());
01157       usubshowerhandle.SetClusterID(ClusterType::kUnknown);
01158       subshowerlist->AddDaughterLink(usubshowerhandle);
01159       MSG("EventSS",Msg::kVerbose) << "U subshower view: " 
01160                                    << usubshowerhandle.GetPlaneView() <<endl;
01161       
01162       cxx.SetDataIn(&vstriparray);
01163       CandSubShowerSRHandle vsubshowerhandle =
01164         CandSubShowerSR::MakeCandidate(subshowerah,cxx);
01165       vsubshowerhandle.SetCandSlice(closestevent->GetCandSlice());
01166       vsubshowerhandle.SetClusterID(ClusterType::kUnknown);
01167       subshowerlist->AddDaughterLink(vsubshowerhandle);
01168       MSG("EventSS",Msg::kVerbose) << "V subshower view: " 
01169                                    << vsubshowerhandle.GetPlaneView() <<endl;
01170       
01171       //add subshowers to list
01172       CandShowerSRListHandle* showerlistSR = 
01173         dynamic_cast<CandShowerSRListHandle*>(showerlist);
01174       TObjArray newshower;
01175       CandSubShowerSRHandle * usubshower = dynamic_cast<CandSubShowerSRHandle*>
01176         (subshowerlist->FindDaughter(&usubshowerhandle));
01177       
01178       CandSubShowerSRHandle * vsubshower = dynamic_cast<CandSubShowerSRHandle*>
01179         (subshowerlist->FindDaughter(&vsubshowerhandle));           
01180       newshower.Add(usubshower);
01181       newshower.Add(vsubshower);
01182       cxx.SetDataIn(&newshower);
01183       
01184       //get ShowerSS alghandle and make shower from subshowers
01185       AlgHandle showerahSS = af.GetAlgHandle("AlgShowerSS","default");
01186       CandShowerSRHandle showerhandleSR = CandShowerSR::MakeCandidate(showerahSS,cxx);
01187       showerhandleSR.SetCandSlice(closestevent->GetCandSlice());
01188       showerlistSR->AddDaughterLink(showerhandleSR);
01189       CandHandle *shw = showerlistSR->FindDaughter(&showerhandleSR);
01190       CandShowerHandle *shwh = dynamic_cast<CandShowerHandle*>(shw);
01191       closestevent->AddShower(shwh);
01192       // Change: cbs 9Aug06
01193       // NEED TO FORCE CREATED SHOWER TO BE PRIMARY
01194       // due to a change in assignment of primary shower, can now create
01195       // a shower here which is *not* called primary by SetPrimaryShower()
01196       // Causes multiple CreatePrimaryShower calls in event builder and 
01197       // large numbers of reco'd showers per events. 
01198       // SetPrimaryShowers is called at the end of RunAlg so a correct
01199       // assignment will be made later
01200       //closestevent->SetPrimaryShower(minShwEFract,shwshwdz);
01201       closestevent->SetPrimaryShower(shwh);
01202     }
01203   }
01204 }

void AlgEventSSList::FillDist2Map AlgConfig ac,
TObjArray &  unassociated,
CandHandle ch,
std::vector< std::map< const CandEventHandle *, Double_t > > &  dist2map
 

Definition at line 1890 of file AlgEventSSList.cxx.

References CandRecoHandle::GetBegPlane(), CandStripHandle::GetBegTime(), CandHandle::GetDaughterIterator(), Registry::GetDouble(), CandStripHandle::GetPlane(), CandStripHandle::GetPlaneView(), CandEventHandle::GetPrimaryShower(), CandEventHandle::GetPrimaryTrack(), CandStripHandle::GetStrip(), CandStripHandle::GetTPos(), CandRecoHandle::GetVtxT(), CandRecoHandle::GetVtxZ(), and MSG.

Referenced by RunAlg().

01890                                                                                                                                                    {
01891   // grab all needed algorithm parameters  
01892   
01893   Double_t pShwTrkDz = ac.GetDouble("ShwTrkDz");
01894   Double_t pHitAssocDTime0 = ac.GetDouble("HitAssocDTime0");
01895   Double_t pHitAssocDTime1 = ac.GetDouble("HitAssocDTime1");
01896   Double_t pHitAssocZParm = ac.GetDouble("HitAssocZParm");
01897   Double_t pHitAssocPParm = ac.GetDouble("HitAssocPParm");
01898   Double_t pHitAssocTParm = ac.GetDouble("HitAssocTParm");
01899 
01900   // construct the map, first looping over strips
01901   for (Int_t i=0; i<=unassociated.GetLast(); i++) {
01902     CandStripHandle *strip =
01903       dynamic_cast<CandStripHandle*>(unassociated.At(i));
01904     TIter eventItr(ch.GetDaughterIterator());
01905     // inner loop is over events
01906     while (CandEventHandle *event =
01907            dynamic_cast<CandEventHandle*>(eventItr())) {
01908       
01909       // calculate separation in space/time
01910       Double_t dtime = (strip->GetBegTime()-event->GetVtxT());
01911       Double_t bestdist2=-1.;
01912       Double_t bestdplane=-1;
01913       Double_t bestdtpos=-1;
01914       Double_t bestdtime=-1;
01915       dist2map[i][event] = 999999999.;
01916       MSG("EventSS",Msg::kVerbose) << i << " " << strip->GetPlane() << " " 
01917                                    << strip->GetStrip() << " stp beg time "
01918                                    << strip->GetBegTime() << " evt vtx time "
01919                                    << event->GetVtxT()
01920                                    << " dtime " 
01921                                    << dtime*1e9 << "/" << pHitAssocDTime0*1e9 
01922                                    << "/" << pHitAssocDTime1*1e9 << endl;  
01923       // if there is a rough time match proceed
01924       if (dtime>pHitAssocDTime0 && dtime<pHitAssocDTime1) {
01925 
01926         Float_t vertexsep = 0;
01927         // If  event has both shower
01928         // and track calc separation between the two.  If this 
01929         // separation is greater
01930         // than pSHwTrkDz, we won't bother adding the strip to the
01931         // primary shower, and will instead consider whether
01932         // to create a new shower at track vertex.
01933         if (event->GetPrimaryShower() && event->GetPrimaryTrack())
01934           vertexsep = fabs(event->GetPrimaryShower()->GetVtxZ() -
01935                           event->GetPrimaryTrack()->GetVtxZ());
01936 
01937         // the event has a shower, and either no track, or track and
01938         // shower vertices are in agreement.  In this case, we
01939         // consider adding this strip to this shower, so calc separation.
01940         if (event->GetPrimaryShower() && 
01941                (!event->GetPrimaryTrack() || vertexsep<pShwTrkDz)) {
01942 
01943           // create iterator of shower daughter strips in same view
01944           // as unassoc strip
01945           CandStripHandleItr shwstripItr(event->GetPrimaryShower()->
01946                                          GetDaughterIterator());
01947           CandStripHandleKeyFunc *shwstripKf = shwstripItr.CreateKeyFunc();
01948           shwstripKf->SetFun(CandStripHandle::KeyFromView);
01949           shwstripItr.GetSet()->AdoptSortKeyFunc(shwstripKf);
01950           shwstripKf = 0;
01951           shwstripItr.GetSet()->Slice();
01952           shwstripItr.GetSet()->Slice(strip->GetPlaneView(),
01953                                       strip->GetPlaneView());
01954           
01955           // iterate over the daughter list, and find minimum distane
01956           // between daughter and unassoc strip
01957           while (CandStripHandle *shwstrip =
01958                         dynamic_cast<CandStripHandle*>(shwstripItr())) {
01959             Double_t dplane = (Double_t)(strip->GetPlane()-shwstrip->GetPlane());
01960             Double_t dtpos = strip->GetTPos()-shwstrip->GetTPos();
01961             Double_t dist2 = (pHitAssocZParm*dplane*dplane) +
01962                              (pHitAssocPParm*dtpos*dtpos) +
01963                              (pHitAssocTParm*dtime*dtime);
01964             if (bestdist2<0. || dist2<bestdist2) {
01965               bestdist2 = dist2;
01966               bestdplane=dplane;
01967               bestdtpos=dtpos;
01968               bestdtime=dtime;
01969             }
01970           }
01971           MSG("EventSS",Msg::kVerbose)
01972             << "primary shower." << " dplane:" << bestdplane <<" dtpos:"
01973             << bestdtpos <<" dtime:" << bestdtime <<" dist2:" << bestdist2<< "\n";
01974           dist2map[i][event] = bestdist2;
01975         } // end if event has shower, etc.
01976         // else if criteria above was not satisfied, consider whether
01977         // strip is near track vertex, in which case we start new shower
01978         else if (event->GetPrimaryTrack()) {
01979 
01980           // calc distance to vertex
01981           Double_t dplane = (Double_t)(strip->GetPlane() - 
01982                                        event->GetPrimaryTrack()->
01983                                        GetBegPlane());
01984           Double_t dtpos=0;
01985           switch (strip->GetPlaneView()) {
01986           case PlaneView::kU:
01987             dtpos=strip->GetTPos()-event->GetPrimaryTrack()->
01988               GetU(event->GetPrimaryTrack()->GetBegPlane());
01989             break;
01990           case PlaneView::kV:
01991             dtpos=strip->GetTPos()-event->GetPrimaryTrack()->
01992               GetV(event->GetPrimaryTrack()->GetBegPlane());
01993             break;
01994           default:
01995             break;
01996           }
01997           Double_t dist2 = ( (pHitAssocZParm*dplane*dplane) +
01998                              (pHitAssocPParm*dtpos*dtpos) +
01999                              (pHitAssocTParm*dtime*dtime) );
02000           dist2map[i][event] = dist2;
02001           MSG("EventSS",Msg::kVerbose)
02002             << "primary track. begin:" << " dplane:" << dplane
02003             << " dtpos:" << dtpos <<" dtime:" << dtime <<" dist2:"
02004             << dist2 << "\n";
02005         } // end if primary track
02006       }  // end if time match OK
02007     }  // end loop over events
02008   }  // end loop over unassoc strips
02009 }

void AlgEventSSList::FillRecoList CandSliceHandle ,
CandShowerHandleItr *  ,
CandTrackHandleItr *  ,
TObjArray & 
 

Definition at line 1776 of file AlgEventSSList.cxx.

References CandRecoHandle::GetCandSlice(), CandHandle::GetNDaughters(), CandHandle::IsCloneOf(), and MSG.

Referenced by RunAlg().

01777                                                                                      {
01778 
01779   if (showerItr) {
01780     showerItr->Reset();
01781     while (CandShowerHandle *shower =
01782            dynamic_cast<CandShowerHandle*>((*showerItr)())) {
01783       if (shower->GetCandSlice()) {
01784         if (slice->IsCloneOf(*(shower->GetCandSlice()))) {
01785           recolist.Add(shower);
01786         }
01787       }
01788       else {
01789         MSG("EventSS", Msg::kError)
01790           << "Shower doesn't contain a slice.  Not added to recolist."
01791           << endl;
01792       }
01793     }
01794   }
01795   
01796   if (trackItr) {
01797     trackItr->Reset();
01798     while (CandTrackHandle *track =
01799            dynamic_cast<CandTrackHandle*>((*trackItr)())) {
01800       if (track->GetCandSlice()) {
01801         if (slice->IsCloneOf(*track->GetCandSlice())) {
01802           if(track->GetNDaughters()==0 && 
01803              track->InheritsFrom("CandFitTrackHandle") &&  
01804              dynamic_cast<CandFitTrackHandle*>(track)->GetFinderTrack()){
01805             track = dynamic_cast<CandFitTrackHandle*>(track)->GetFinderTrack();    
01806             MSG("EventSS", Msg::kDebug)<< "Finder Track being used " << endl;
01807           }
01808           recolist.Add(track);
01809         }
01810       }
01811       else {
01812         MSG("EventSS", Msg::kError)
01813           << "Track doesn't contain a slice.  Not added to recolist."
01814           << endl;
01815       }
01816     }
01817   }
01818 }

void AlgEventSSList::FindUnassociated CandSliceHandleItr &  sliceItr,
CandEventHandleItr *  ,
TObjArray &  unassociated
 

Definition at line 1843 of file AlgEventSSList.cxx.

References CandHandle::FindDaughter(), CandEventHandle::GetCandSlice(), CandHandle::GetDaughterIterator(), and CandHandle::IsCloneOf().

Referenced by RunAlg().

01845                                                                {
01846 
01847   while (CandSliceHandle *slice =
01848          dynamic_cast<CandSliceHandle*>(sliceItr())) {
01849     TIter stripItr(slice->GetDaughterIterator());
01850     while (CandStripHandle *strip =
01851            dynamic_cast<CandStripHandle*>(stripItr())) {
01852       Bool_t found=false;
01853       // check for duplicates      
01854       for (Int_t i=0; i<=unassociated.GetLast() && !found; i++) {
01855         CandStripHandle *strip0 =
01856           dynamic_cast<CandStripHandle*>(unassociated.At(i));
01857         if (*strip == *strip0) {
01858           found = true;
01859           break;
01860         }
01861       }
01862       // if this strip not already in unnassociated list, then
01863       // loop over objects in this slice, and check against all daughter strips
01864       if (!found) {
01865         if (eventItr) {
01866           eventItr->Reset();
01867           while (CandEventHandle *event =
01868                  dynamic_cast<CandEventHandle*>((*eventItr)())) {
01869             if (event->GetCandSlice()) {
01870               if (slice->IsCloneOf(*(event->GetCandSlice()))) { 
01871                 if (event->FindDaughter(strip)) {
01872                   found = true;
01873                   break;
01874                 }
01875               }
01876             }
01877           }
01878         }
01879       }
01880       // this strip is not in any event, add to unassociated list
01881       if (!found) {
01882         unassociated.Add(strip);
01883       }
01884     }
01885   }
01886 }

void AlgEventSSList::ReConstructShowers AlgConfig ac,
CandHandle ch,
CandContext cxx
 

This method fills a cluster list for every shower in the snarl, and then reruns the shower algorithm. Following this, it tracks the primary track in each event through the primary shower, and subtracts 1 MIP for each strip passed through.

Definition at line 1235 of file AlgEventSSList.cxx.

References CandHandle::AddDaughterLink(), CandShowerHandle::CalibrateEnergy(), CandSubShowerSRHandle::DupHandle(), Registry::Get(), AlgHandle::GetAlgConfig(), AlgFactory::GetAlgHandle(), CandEventHandle::GetCandSlice(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), AlgFactory::GetInstance(), CandEventHandle::GetLastShower(), CandShowerSRHandle::GetLastSubShower(), CandEventHandle::GetPrimaryTrack(), CandEventHandle::GetShowerWritable(), CandShowerSRHandle::GetSubShower(), CandHandle::GetUidInt(), CandHandle::IsCloneOf(), CandHandle::IsEqual(), MSG, CandHandle::RemoveDaughter(), AlgHandle::RunAlg(), CandRecoHandle::SetCandSlice(), and CandContext::SetDataIn().

Referenced by RunAlg().

01238 {
01239   MSG("EventSS",Msg::kVerbose) << "In ReConstructShowers: " << endl;
01240   // Retrieve persistant modifiable CandShowerList
01241   CandShowerListHandle *showerlist = const_cast<CandShowerListHandle*>
01242     (dynamic_cast<const CandShowerListHandle*>(cxx.GetDataIn()));
01243 
01244   if(!showerlist->InheritsFrom("CandShowerSRListHandle")) {
01245     MSG("EventSS",Msg::kWarning)
01246       << "ListHandle is not a CandShowerSRListHandle."
01247       << "  ShowerList not processed in ReConstructShowers." << endl;
01248     return;
01249   }
01250 
01251   // Get singleton instance of AlgFactory
01252   AlgFactory &af = AlgFactory::GetInstance();
01253   
01254   // Get an AlgHandle to AlgShowerSR with default AlgConfig
01255   const char *pEventAlgConfig = 0;
01256   ac.Get("EventAlgConfig",pEventAlgConfig);
01257 
01258   AlgHandle ahss = af.GetAlgHandle("AlgShowerSS","default");
01259   AlgConfig acss = ahss.GetAlgConfig();
01260   Int_t eventCounter = 0;
01261   TIter evItr(ch.GetDaughterIterator());
01262   while (CandEventHandle *ev =
01263          dynamic_cast<CandEventHandle*>(evItr())) {  //event loop
01264     MSG("EventSS",Msg::kVerbose) << "event = " << eventCounter << endl;
01265     for (Int_t ishower=0; ishower<ev->GetLastShower()+1; ishower++) { //shower loop
01266       MSG("EventSS",Msg::kVerbose) << "shower = " << ishower << endl;
01267       CandShowerHandle *shower =
01268         dynamic_cast<CandShowerHandle*>(ev->GetShowerWritable(ishower));
01269       CandShowerSRHandle *showerSR = 0;
01270       if(shower->InheritsFrom("CandShowerSRHandle")) {
01271         showerSR = dynamic_cast<CandShowerSRHandle*>(shower);
01272       }
01273       MSG("EventSS",Msg::kVerbose) << "Shower UID = " << showerSR->GetUidInt() << endl;
01274       //assume subshower chain
01275       if(showerSR) {  
01276         TObjArray newshower;
01277         // fill array of subshowers for this shower
01278         for (Int_t isubshower=0;
01279              isubshower<showerSR->GetLastSubShower()+1;
01280              isubshower++) {
01281           MSG("EventSS",Msg::kVerbose) << "subshower = " << isubshower << endl;
01282           const CandSubShowerSRHandle *ssh = 
01283             showerSR->GetSubShower(isubshower);
01284           if (ssh) {
01285             // New handle on heap
01286             CandSubShowerSRHandle *sshd = ssh->DupHandle();
01287             // Add new heap-based CandHandle to TObjArray
01288             newshower.Add(sshd);
01289             MSG("EventSS",Msg::kDebug) << "Uid's (old,new) " 
01290                                        << ssh->GetUidInt() << " " 
01291                                        << sshd->GetUidInt() << endl;
01292           }
01293         }
01294         // re-run algorithm for this shower
01295         if (newshower.GetEntriesFast() > 0) {
01296           cxx.SetDataIn(&newshower);
01297           ahss.RunAlg(*showerSR,cxx);
01298         }
01299         else {
01300           MSG("EventSS", Msg::kError)
01301             << "Attempt to pass empty TObjArray newshower to AlgShowerSS"
01302             << " in CandContext.  AlgShowerSS::RunAlg() call ignored."
01303             << endl;
01304         }
01305         // Delete heap-based CandHandles in TObjArray
01306         newshower.Delete();
01307       }
01308       else {
01309         MSG("EventSS",Msg::kWarning)
01310           << "Handle is not a CandShowerSRHandle."
01311           << "  Shower not processed in ReConstructShowers."
01312           << endl;
01313         continue;
01314       }
01315       MSG("EventSS",Msg::kVerbose) << "Modified shower UID = " 
01316                                    << showerSR->GetUidInt() << endl; 
01317       // now loop through strips in this shower, and determine whether primary
01318       // track in this event intercepts this strip.  If so, subtract 1 mip
01319       // from shower energy.
01320       
01321       // if event has primary track, extrapolate through shower 
01322       // and subtract energy loss to recalibrate 
01323       CandTrackHandle *primaryTrack = ev->GetPrimaryTrack();
01324       if (primaryTrack) {
01325         shower->CalibrateEnergy(primaryTrack,acss);
01326       }
01327       
01328       // Update modified showers in persistent modifiable CandShowerList.
01329       // CandHandle could potentially supply a ReplaceDaughter() function.
01330       // Interface of ReplaceDaughter() would be debatable.
01331       if (showerlist) {
01332         TIter shiter(showerlist->GetDaughterIterator());
01333         CandShowerHandle *target;
01334         Bool_t found = kFALSE;
01335         while (!found &&
01336                (target = dynamic_cast<CandShowerHandle*>(shiter()))) {
01337           if (shower->IsCloneOf(*target)) {  // Tests clone or Cand ==
01338             found = kTRUE;
01339             
01340             // Replace target only if shower is a modified version of target
01341             if (!(shower->IsEqual(target))) {   // CandBase inequality
01342               shower->SetCandSlice(ev->GetCandSlice());
01343               if (!showerlist->RemoveDaughter(target)) {    // Failure
01344                 MSG("EventSS",Msg::kError) 
01345                   << endl
01346                   << "Failure of ShowerList daughter removal " << endl
01347                   << "attempted during replacement by modified Shower."
01348                   << endl << "Will result in double counted Shower."
01349                   << endl;
01350               }
01351               showerlist->AddDaughterLink(*shower);
01352             }
01353           }
01354         }
01355         
01356         // Add new (unfound) shower to persistent modifiable CandShowerList.
01357         if (!found){
01358           shower->SetCandSlice(ev->GetCandSlice());
01359           showerlist->AddDaughterLink(*shower);
01360         }
01361       }
01362     }
01363     eventCounter++;
01364   }
01365 }

void AlgEventSSList::ReFillDist2Map AlgConfig ac,
CandStripHandle closeststrip,
CandEventHandle closestevent,
TObjArray &  unassociated,
std::vector< std::map< const CandEventHandle *, Double_t > > &  dist2map
 

Definition at line 1207 of file AlgEventSSList.cxx.

References CandStripHandle::GetBegTime(), Registry::GetDouble(), CandStripHandle::GetPlane(), CandStripHandle::GetPlaneView(), CandStripHandle::GetTPos(), and CandRecoHandle::GetVtxT().

Referenced by RunAlg().

01208 {
01209   Double_t pHitAssocZParm = ac.GetDouble("HitAssocZParm");
01210   Double_t pHitAssocPParm = ac.GetDouble("HitAssocPParm");
01211   Double_t pHitAssocTParm = ac.GetDouble("HitAssocTParm");
01212 
01213   for (Int_t i=0; i<=unassociated.GetLast(); i++) {
01214     CandStripHandle *strip =
01215       dynamic_cast<CandStripHandle*>(unassociated.At(i));
01216     if(strip->GetPlaneView()==closeststrip->GetPlaneView()){
01217       Double_t dtime = strip->GetBegTime()-closestevent->GetVtxT();
01218       Double_t dplane =
01219         (Double_t)(strip->GetPlane()-closeststrip->GetPlane());
01220       Double_t dtpos = strip->GetTPos()-closeststrip->GetTPos();
01221       Double_t dist2 = (pHitAssocZParm*dplane*dplane) +
01222         (pHitAssocPParm*dtpos*dtpos) +
01223         (pHitAssocTParm*dtime*dtime);
01224       if (dist2<dist2map[i][closestevent]) {
01225         dist2map[i][closestevent] = dist2;
01226       }
01227     }
01228   } // loop over unassoc hits
01229 }

void AlgEventSSList::RunAlg AlgConfig ac,
CandHandle ch,
CandContext cx
[virtual]
 

Implements AlgBase.

Definition at line 72 of file AlgEventSSList.cxx.

References abs(), CandHandle::AddDaughterLink(), AddStripToEvent(), CandShowerSRHandle::BelongsWithShower(), CandTrackHandle::BelongsWithTrack(), CandShowerSRHandle::BelongsWithTrack(), BuildEventFromUnassoc(), CandShowerSRHandle::BuriedTrack(), CreatePrimaryShower(), FillDist2Map(), FillRecoList(), FindUnassociated(), Registry::Get(), AlgFactory::GetAlgHandle(), CandRecoHandle::GetBegPlane(), CandStripHandle::GetBegTime(), CandContext::GetCandRecord(), CandRecoHandle::GetCharge(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), CandRecoHandle::GetDirCosU(), CandRecoHandle::GetDirCosV(), Registry::GetDouble(), CandRecoHandle::GetEndPlane(), AlgFactory::GetInstance(), Registry::GetInt(), CandEventHandle::GetLastShower(), CandContext::GetMom(), CandHandle::GetNDaughters(), CandRecoHandle::GetNStrip(), CandStripHandle::GetPlane(), CandEventHandle::GetPrimaryShower(), CandEventHandle::GetPrimaryTrack(), CandEventHandle::GetShower(), CandEventHandle::GetShowerWritable(), CandStripHandle::GetTPos(), RecMinos::GetVldContext(), CandRecoHandle::GetVtxPlane(), CandRecoHandle::GetVtxT(), CandHandle::IsCloneOf(), CandTrackHandle::IsUnphysical(), CandShowerSRHandle::IsUnphysical(), CandEvent::MakeCandidate(), max, min, MSG, ChopHelp::nchop, PITCHINMETRES, ChopHelp::Print(), ReConstructShowers(), ReFillDist2Map(), CandEventHandle::SetCandSlice(), and SetPrimaryShowers().

00073 {
00074 
00075   assert(cx.GetDataIn());
00076   if (!(cx.GetDataIn()->InheritsFrom("TObjArray"))) return;
00077 
00078   CandContext cxx(this,cx.GetMom());
00079   AlgFactory &af = AlgFactory::GetInstance();
00080   const char *pEventAlgConfig = 0;
00081   ac.Get("EventAlgConfig",pEventAlgConfig);
00082   AlgHandle ah = af.GetAlgHandle("AlgEventSR",pEventAlgConfig);
00083   // grab all needed algorithm parameters  
00084 
00085   //generic parameters:
00086   Double_t pTrkTrkDtpos2 = ac.GetDouble("TrkTrkDtpos2");
00087   Double_t pTrkTrkDz = ac.GetDouble("TrkTrkDz");
00088   Double_t pTrkTrkDt = ac.GetDouble("TrkTrkDt");
00089   Double_t pShwTrkDtpos2 = ac.GetDouble("ShwTrkDtpos2");
00090   Double_t pShwTrkDz = ac.GetDouble("ShwTrkDz");
00091   Double_t pShwTrkDt = ac.GetDouble("ShwTrkDt");
00092   Double_t pShwShwDtpos2 = ac.GetDouble("ShwShwDtpos2");
00093   Double_t pShwShwDz = ac.GetDouble("ShwShwDz");
00094   Double_t pShwShwDt = ac.GetDouble("ShwShwDt");
00095   Double_t pMaxNewShwLen = ac.GetDouble("MaxNewShwLen");
00096   Double_t pHitAssocDTime0 = ac.GetDouble("HitAssocDTime0");
00097   Double_t pHitAssocDTime1 = ac.GetDouble("HitAssocDTime1");
00098   Double_t pHitAssocMaxDist2 = ac.GetDouble("HitAssocMaxDist2");
00099   Double_t pMinTrackIsolation = ac.GetDouble("MinTrackIsolation");
00100 
00101   //SS specific parameters
00102   Double_t pHardBuriedFrac = ac.GetDouble("HardBuriedFrac");
00103   Double_t pSoftBuriedFrac = ac.GetDouble("SoftBuriedFrac");
00104   Double_t pDCosUVCut = ac.GetDouble("DCosUVCut");
00105   Double_t pTrkGapFracCut = ac.GetDouble("TrkGapFracCut");
00106   Double_t pTrkAsymCut = ac.GetDouble("TrkAsymCut");
00107   Double_t pTrkXTalkFracCut = ac.GetDouble("TrkXTalkFracCut");
00108   Double_t pTrkXTalkDef = ac.GetDouble("TrkXTalkDef");
00109   Double_t pShwXTalkFracCut = ac.GetDouble("ShwXTalkFracCut");
00110   Double_t pShwXTalkDef = ac.GetDouble("ShwXTalkDef");
00111   Int_t useChopper = ac.GetInt("UseChopper");
00112 
00113   MSG("EventSS",Msg::kDebug) << pHardBuriedFrac << " " << pSoftBuriedFrac << " "
00114                             << pDCosUVCut << " " << pTrkGapFracCut << " "
00115                             << pTrkAsymCut << " " << pTrkXTalkFracCut << " "
00116                             << pTrkXTalkDef << " " << pShwXTalkFracCut << " "
00117                             << pShwXTalkDef << endl;
00118 
00119   // grab all needed input lists
00120   const CandRecord *candrec = cx.GetCandRecord();
00121   assert(candrec);
00122   const VldContext *vldcptr = candrec->GetVldContext();
00123   assert(vldcptr);
00124   VldContext vldc = *vldcptr;
00125   const CandSliceListHandle *slicelist = 0;
00126   CandShowerListHandle *showerlist = 0;
00127   const CandTrackListHandle *tracklist = 0;
00128   CandSubShowerSRListHandle *subshowerlist = 0;
00129   const TObjArray *cxin = dynamic_cast<const TObjArray *>(cx.GetDataIn());
00130   for (Int_t i=0; i<=cxin->GetLast(); i++) {
00131     TObject *tobj = cxin->At(i);
00132     if (tobj->InheritsFrom("CandSliceListHandle")) {
00133       slicelist = dynamic_cast<CandSliceListHandle*>(tobj);
00134     }
00135     if (tobj->InheritsFrom("CandShowerListHandle")) {
00136       showerlist = dynamic_cast<CandShowerListHandle*>(tobj);
00137     }
00138     if (tobj->InheritsFrom("CandTrackListHandle")) {
00139       tracklist = dynamic_cast<CandTrackListHandle*>(tobj);
00140     }
00141     if (tobj->InheritsFrom("CandSubShowerSRListHandle")) {
00142       subshowerlist = dynamic_cast<CandSubShowerSRListHandle*>(tobj);
00143     }
00144   }
00145   if (!slicelist) {
00146     MSG("EventSS",Msg::kError) << "CandSliceListHandle missing\n";
00147     return;
00148   }
00149 
00150   CandShowerHandleItr * showerItr=0;
00151   // create iterators for track, shower, and cluster lists, keying on slice
00152   if (showerlist && showerlist->GetNDaughters()>0) {
00153     showerItr = new CandShowerHandleItr(showerlist->GetDaughterIterator());
00154     CandShowerHandleKeyFunc *showerKf = showerItr->CreateKeyFunc();
00155     showerKf->SetFun(CandShowerHandle::KeyFromSlice);
00156     showerItr->GetSet()->AdoptSortKeyFunc(showerKf);
00157     showerKf = 0;
00158   }
00159 
00160   CandTrackHandleItr * trackItr = 0;
00161   if (tracklist && tracklist->GetNDaughters()>0) {
00162     trackItr = new CandTrackHandleItr(tracklist->GetDaughterIterator());
00163     CandTrackHandleKeyFunc *trackKf = trackItr->CreateKeyFunc();
00164     trackKf->SetFun(CandTrackHandle::KeyFromSlice);
00165     trackItr->GetSet()->AdoptSortKeyFunc(trackKf);
00166     trackKf = 0;
00167   }
00168 
00169   CandSubShowerSRHandleItr * subshowerItr = 0;
00170   if (subshowerlist && subshowerlist->GetNDaughters()>0) {
00171     subshowerItr = new CandSubShowerSRHandleItr(subshowerlist->
00172                                                 GetDaughterIterator());
00173     CandSubShowerSRHandleKeyFunc *subshowerKf = subshowerItr->CreateKeyFunc();
00174     subshowerKf->SetFun(CandSubShowerSRHandle::KeyFromSlice);
00175     subshowerItr->GetSet()->AdoptSortKeyFunc(subshowerKf);
00176     subshowerKf = 0;
00177   }
00178 
00179   Int_t nslice=0;
00180 
00181   // loop over slices
00182   CandSliceHandleItr sliceItr(slicelist->GetDaughterIterator());
00183   while (CandSliceHandle *slice =
00184          dynamic_cast<CandSliceHandle*>(sliceItr())) { //loop over slices
00185   
00186     ++nslice;
00187     MSG("EventSS", Msg::kDebug)
00188       << " ****** Slice " << nslice << " *******"  
00189       << endl;   
00190 
00191     //TObjArray newevent;
00192     TObjArray recolist;
00193     TObjArray eventlist;
00194 
00195     //construct list of all showers and tracks in this slice
00196     FillRecoList(slice,showerItr,trackItr,recolist);
00197 
00198     //first try to remove tracks buried in showers    
00199     std::vector<Bool_t> removeMe(recolist.GetLast()+1);
00200     for (Int_t i=0;i<=recolist.GetLast();i++) {
00201       removeMe[i] = false;
00202       if(recolist.At(i)->InheritsFrom("CandTrackHandle")) {
00203         CandTrackHandle *track = 
00204           dynamic_cast<CandTrackHandle*>(recolist.At(i));
00205         Bool_t removeTrack = false;
00206         for (Int_t j=0;j<=recolist.GetLast();j++) {         
00207           if (recolist.At(j)->InheritsFrom("CandShowerSRHandle")) {
00208             CandShowerSRHandle * shower = 
00209               dynamic_cast<CandShowerSRHandle*>(recolist.At(j));
00210             if(shower->BuriedTrack(track,pMinTrackIsolation)) {
00211               removeTrack = true;
00212               break;
00213             }
00214           }
00215         }
00216         if(removeTrack){
00217           MSG("EventSS",Msg::kDebug) 
00218             << " REMOVING TRACK" << endl;
00219           removeMe[i] = true;
00220         }
00221       }
00222     }
00223     for (Int_t i=0;i<=recolist.GetLast();i++) {
00224       if(removeMe[i]) recolist.RemoveAt(i);
00225     }
00226     recolist.Compress();
00227 
00228     MSG("EventSS",Msg::kDebug) << "Number of reco objects in slice: " 
00229                                << recolist.GetLast()+1 << endl;
00230     //now make matrix to hold match results between all objects.
00231     //this also includes buried tracks that have not been removed
00232     //NOTE: matching requires both views match;
00233     //      burial just counts up matching planes in either view;
00234     //  One view being significantly buried may indicate a bad track, which
00235     //  should be associated with the shower, not used to make a new event
00236     std::vector<std::vector<Int_t> > objectMatch(recolist.GetLast()+1,0);
00237     //initialize:
00238     for(int i=0;i<=recolist.GetLast();i++){
00239       for(int j=0;j<=recolist.GetLast();j++) {
00240         objectMatch[i].push_back(0);
00241       }
00242     }
00243 
00244     //use this vector to keep track of objects already assigned
00245     std::vector<Bool_t> objectUsed(recolist.GetLast()+1,0);
00246     
00247     //loop over all objects
00248     for(Int_t i=0;i<=recolist.GetLast();i++) { //loop over i
00249       objectUsed[i] = false;  //initialization
00250       
00251       CandShowerSRHandle *shower=0;
00252       CandTrackHandle *track=0;
00253       if (recolist.At(i)->InheritsFrom("CandShowerSRHandle"))
00254         shower = dynamic_cast<CandShowerSRHandle*>(recolist.At(i));
00255       if (recolist.At(i)->InheritsFrom("CandTrackHandle"))
00256         track = dynamic_cast<CandTrackHandle*>(recolist.At(i));
00257       
00258       //loop over all later objects only (since matching matrix is symmetric)
00259       for(Int_t j=i+1;j<=recolist.GetLast();j++) { //loop over j
00260 
00261         MSG("EventSS",Msg::kDebug) << "Checking objectMatch between " 
00262                                    << i << " and " << j << endl;
00263 
00264         CandShowerSRHandle *shower2=0;
00265         CandTrackHandle *track2=0;
00266         if (recolist.At(j)->InheritsFrom("CandShowerSRHandle"))
00267           shower2 = dynamic_cast<CandShowerSRHandle*>(recolist.At(j));
00268         if (recolist.At(j)->InheritsFrom("CandTrackHandle"))
00269           track2 = dynamic_cast<CandTrackHandle*>(recolist.At(j));
00270         
00271         if(shower) {
00272           if(shower2) {
00273             if(!shower->IsUnphysical(pShwXTalkFracCut,pShwXTalkDef) && 
00274                !shower2->IsUnphysical(pShwXTalkFracCut,pShwXTalkDef)) {
00275               if(shower->BelongsWithShower(shower2,ac,vldcptr,pShwShwDtpos2,
00276                                            pShwShwDz,pShwShwDt)) objectMatch[i][j] = 1;
00277             }
00278             objectMatch[j][i] = objectMatch[i][j];
00279           }
00280           else if(track2) {
00281             if(shower->BelongsWithTrack(track2,ac,vldcptr,pShwTrkDtpos2,
00282                                         pShwTrkDz,pShwTrkDt)) objectMatch[i][j] = 1;
00283 
00284             if(objectMatch[i][j]!=1) {
00285               //check track isn't buried:
00286               Int_t *stats = new Int_t[8];
00287               shower->BuriedTrack(track2,1000,stats); //1000 forces all calculations
00288               if(stats[3]>0) { //at least one buried track planes
00289                 Int_t nplanes = track2->GetEndPlane() - track2->GetBegPlane() + 1;
00290                 MSG("EventSS",Msg::kDebug) 
00291                   << "BuriedTrack Stats: " << "\n"
00292                   << "Number of track planes = " << nplanes << "\n"
00293                   << "nSharedPlanes = " << stats[0] << "\n" 
00294                   << "nMissingIsoPlanes = " << stats[1] << "\n" 
00295                   << "nMissingSharedPlanes = " << stats[2] << "\n" 
00296                   << "nBuriedPlanes = " << stats[3] << "\n" 
00297                   << "nPhysBuriedPlanes = " << stats[4] << "\n" 
00298                   << "nTrkLike = " << stats[5] << "\n" 
00299                   << "nTrkLikeTopol = " << stats[6] << "\n" 
00300                   << "nXTalk = " << stats[7] << endl;
00301                 nplanes -= (stats[1] + stats[2]); //subtract missing track planes
00302                 //if more than half of the track planes are in the shower core: match
00303                 if(Float_t(stats[4])/Float_t(nplanes) > pHardBuriedFrac)
00304                   objectMatch[i][j] = -1;
00305                 //else if more than half the track planes are anywhere in 
00306                 //the shower or are x-talk like and are not "trklike": match
00307                 else if(Float_t(stats[3]-stats[5]-stats[6])/Float_t(nplanes) > 
00308                         pHardBuriedFrac) objectMatch[i][j] = -1;
00309                 //check whether track passes through shower at high angle in either view
00310                 //expect significant number of buried planes
00311                 else if(stats[3]>=Int_t(nplanes*pSoftBuriedFrac) && 
00312                         ( TMath::Abs(track2->GetDirCosU())>pDCosUVCut || 
00313                           TMath::Abs(track2->GetDirCosV())>pDCosUVCut) ) {
00314                   objectMatch[i][j] = -1;
00315                 }
00316               }
00317               delete [] stats;
00318             }
00319             objectMatch[j][i] = objectMatch[i][j];
00320           }
00321         }
00322         
00323         if(track) {
00324           if(shower2) {
00325             if(shower2->BelongsWithTrack(track,ac,vldcptr,pShwTrkDtpos2,
00326                                          pShwTrkDz,pShwTrkDt)) objectMatch[i][j] = 1;
00327             if(objectMatch[i][j]!=1) {
00328               //check track isn't buried:
00329               Int_t *stats = new Int_t[8];
00330               shower2->BuriedTrack(track,1000,stats); //1000 forces all calculations
00331               if(stats[3]>0) { //at least one buried track planes
00332                 Int_t nplanes = track->GetEndPlane() - track->GetBegPlane() + 1;
00333                 MSG("EventSS",Msg::kDebug) 
00334                   << "BuriedTrack Stats: " << "\n"
00335                   << "Number of track planes = " << nplanes << "\n"
00336                   << "nSharedPlanes = " << stats[0] << "\n" 
00337                   << "nMissingIsoPlanes = " << stats[1] << "\n" 
00338                   << "nMissingSharedPlanes = " << stats[2] << "\n" 
00339                   << "nBuriedPlanes = " << stats[3] << "\n" 
00340                   << "nPhysBuriedPlanes = " << stats[4] << "\n" 
00341                   << "nTrkLike = " << stats[5] << "\n" 
00342                   << "nTrkLikeTopol = " << stats[6] << "\n" 
00343                   << "nXTalk = " << stats[7] << endl;
00344                 nplanes -= (stats[1] + stats[2]); //subtract missing track planes
00345                 //if more than half of the track planes are in the shower core: match
00346                 if(Float_t(stats[4])/Float_t(nplanes) > pHardBuriedFrac)
00347                   objectMatch[i][j] = -1;
00348                 //else if more than half the track planes are anywhere in 
00349                 //the shower or are x-talk like and are not "trklike": match
00350                 else if(Float_t(stats[3]-stats[5]-stats[6])/Float_t(nplanes) > 
00351                         pHardBuriedFrac) objectMatch[i][j] = -1;
00352                 //check whether track passes through shower at high angle in either view
00353                 //expect significant number of buried planes
00354                 else if(stats[3]>=Int_t(nplanes*pSoftBuriedFrac) && 
00355                         (TMath::Abs(track2->GetDirCosU())>pDCosUVCut || 
00356                          TMath::Abs(track2->GetDirCosV())>pDCosUVCut) ) {
00357                   objectMatch[i][j] = -1;
00358                 }
00359               }
00360               delete [] stats;
00361             }
00362             objectMatch[j][i] = objectMatch[i][j];
00363           }
00364           else if(track2) {
00365             if(track->BelongsWithTrack(track2,ac,vldcptr,pTrkTrkDtpos2,
00366                                        pTrkTrkDz,pTrkTrkDt)) objectMatch[i][j] = 1;
00367             objectMatch[j][i] = objectMatch[i][j];
00368           }
00369         }
00370         MSG("EventSS",Msg::kDebug) << "Match for objects " << i << " and " 
00371                                    << j << " : " << objectMatch[i][j] << endl;
00372       }
00373     }
00374 
00375     MSG("EventSS",Msg::kDebug) << "Starting to form events" << endl;
00376 
00377     //start to make events:
00378     Int_t nRecoObjects = recolist.GetLast()+1;
00379     //get total charge from each object
00380     std::vector<Float_t> ph(nRecoObjects,0);
00381     for(int i=0;i<nRecoObjects;i++){
00382       if (recolist.At(i)->InheritsFrom("CandShowerSRHandle")) {
00383         CandShowerSRHandle *shower = dynamic_cast<CandShowerSRHandle*>(recolist.At(i));
00384         ph[i] = shower->GetCharge(CalStripType::kMIP);
00385       }
00386       if (recolist.At(i)->InheritsFrom("CandTrackHandle")) {
00387         CandTrackHandle *track = dynamic_cast<CandTrackHandle*>(recolist.At(i));
00388         //ph[i] = track->GetCharge(CalStripType::kMIP)+1e7; //force tracks to come first
00389         //force tracks to come first, order tracks by number of strips
00390         ph[i] = track->GetNStrip()+1e7; 
00391       }
00392     }
00393 
00394     Int_t nEvents = 0;
00395     Int_t nUnusedObjects = nRecoObjects;
00396     while(nUnusedObjects>0) { //loop until all reco objects are in events
00397       
00398       //find largest object:
00399       Float_t maxPH = 0;
00400       Int_t maxPHIndex = -1;
00401       for(int i=0;i<nRecoObjects;i++){
00402         if(ph[i]>maxPH) {
00403           maxPH = ph[i];
00404           maxPHIndex = i;
00405         }
00406       }
00407       
00408       //add largest objects first:
00409       TObjArray *event = new TObjArray;
00410       event->Add(recolist.At(maxPHIndex));     
00411       ph[maxPHIndex] = 0;
00412       objectUsed[maxPHIndex] = true;
00413       nUnusedObjects-=1;
00414 
00415       MSG("EventSS",Msg::kDebug) << "Added object " << maxPHIndex 
00416                                  << " to event " << nEvents << endl;
00417       
00418       if(recolist.At(maxPHIndex)->InheritsFrom("CandTrackHandle")) {
00419         MSG("EventSS",Msg::kDebug) << "Object " << maxPHIndex << " is a track" << endl;
00420 
00421         //maxPH object is a track
00422         CandTrackHandle *track = dynamic_cast<CandTrackHandle*>(recolist.At(maxPHIndex));
00423         
00424         //now add objects clearly associated with largest object:
00425         for(int i=0;i<nRecoObjects;i++){
00426           if(objectMatch[maxPHIndex][i]!=0 && !objectUsed[i]) {       
00427             if(recolist.At(i)->InheritsFrom("CandTrackHandle")) { 
00428               MSG("EventSS",Msg::kDebug) << "Track match. Adding object " << i 
00429                                          << " to event " << nEvents << endl;
00430               //if this match is a track, add it to event
00431               event->Add(recolist.At(i));
00432               ph[i] = 0;
00433               objectUsed[i] = true;
00434               nUnusedObjects-=1;
00435             } // end if for CandTrackHandle
00436             else if(recolist.At(i)->InheritsFrom("CandShowerSRHandle")) {           
00437               MSG("EventSS",Msg::kDebug) << "Shower match. Testing whether to add object " 
00438                                          << i << " to event " << nEvents << endl;
00439               CandShowerSRHandle *shower2 = dynamic_cast<CandShowerSRHandle*>(recolist.At(i));
00440               //if this match is a shower, check where shower is relative to track vertex
00441               //if(TMath::Abs(shower2->GetVtxZ()-track->GetVtxZ())<pShwTrkDz/2.){
00442               if(TMath::Abs(shower2->GetBegPlane(PlaneView::kU) - 
00443                             track->GetBegPlane(PlaneView::kU)) < 
00444                  pShwTrkDz/(PITCHINMETRES*2.) || 
00445                  TMath::Abs(shower2->GetBegPlane(PlaneView::kV) - 
00446                             track->GetBegPlane(PlaneView::kV)) < 
00447                  pShwTrkDz/(PITCHINMETRES*2.)){
00448                 MSG("EventSS",Msg::kDebug) << "Shower " << i << " is close to track vertex. "
00449                                            << "Adding to event " << nEvents << endl;
00450                 event->Add(recolist.At(i));
00451                 ph[i] = 0;
00452                 objectUsed[i] = true;
00453                 nUnusedObjects-=1;
00454                 //look for objects associated with this vertex shower
00455                 for(int j=0;j<nRecoObjects;j++){
00456                   if(j!=maxPHIndex && objectMatch[i][j]!=0 && !objectUsed[j]) { //got a match
00457                     if(objectMatch[maxPHIndex][j]==0) { // j does not match maxPHIndex
00458                       if(recolist.At(j)->InheritsFrom("CandTrackHandle")) {
00459                         MSG("EventSS",Msg::kDebug) << "Shower " << i 
00460                                                    << " has another track match. " 
00461                                                    << "Test how buried track " << j
00462                                                    << " is in shower" << endl; 
00463                         if(objectMatch[i][j]==-1){
00464                           event->Add(recolist.At(j));
00465                           ph[j] = 0;
00466                           objectUsed[j] = true;
00467                           nUnusedObjects-=1;
00468                         }
00469                         else {
00470                           CandTrackHandle *track2 = 
00471                             dynamic_cast<CandTrackHandle*>(recolist.At(j));
00472                           //check if track2 is buried
00473                           Int_t *stats = new Int_t[8];
00474                           shower2->BuriedTrack(track2,1000,stats);
00475                           if(stats[3]>0){
00476                             Int_t nplanes = track2->GetEndPlane() - track2->GetBegPlane() + 1;
00477                             MSG("EventSS",Msg::kDebug) 
00478                               << "BuriedTrack Stats: " << "\n"
00479                               << "Number of track planes = " << nplanes << "\n"
00480                               << "nSharedPlanes = " << stats[0] << "\n" 
00481                               << "nMissingIsoPlanes = " << stats[1] << "\n" 
00482                               << "nMissingSharedPlanes = " << stats[2] << "\n" 
00483                               << "nBuriedPlanes = " << stats[3] << "\n" 
00484                               << "nPhysBuriedPlanes = " << stats[4] << "\n" 
00485                               << "nTrkLike = " << stats[5] << "\n" 
00486                               << "nTrkLikeTopol = " << stats[6] << "\n" 
00487                               << "nXTalk = " << stats[7] << endl;
00488                             nplanes -= (stats[1] + stats[2]);
00489                             if( (Float_t(stats[4])/Float_t(nplanes) > pHardBuriedFrac) ||
00490                                 (Float_t(stats[3]-stats[5]-stats[6])/Float_t(nplanes) > 
00491                                  pHardBuriedFrac) ||
00492                                 (stats[3]>=Int_t(nplanes*pSoftBuriedFrac) && 
00493                                  (TMath::Abs(track2->GetDirCosU())>pDCosUVCut || 
00494                                   TMath::Abs(track2->GetDirCosV())>pDCosUVCut) ) ) {
00495                               MSG("EventSS",Msg::kDebug) << "Adding track " << j 
00496                                                          << " to event " << nEvents << endl;
00497                               event->Add(recolist.At(j));
00498                               ph[j] = 0;
00499                               objectUsed[j] = true;
00500                               nUnusedObjects-=1;
00501                             }
00502                           }
00503                           delete [] stats;
00504                         }
00505                       }
00506                       else if(recolist.At(j)->InheritsFrom("CandShowerSRHandle")) {
00507                         MSG("EventSS",Msg::kDebug) << "Shower " << i 
00508                                                    << " has another shower match. " 
00509                                                    << "Test whether shower " << j 
00510                                                    << " has a better association" << endl; 
00511                         Bool_t isAssocTrack = false; //is this shower assoc with a track
00512                         for(int k=0;k<nRecoObjects;k++) { 
00513                           //check for match
00514                           if(k!=maxPHIndex && k!=i && objectMatch[j][k]!=0){
00515                             //check for track
00516                             if(recolist.At(k)->InheritsFrom("CandTrackHandle"))
00517                               //check it does not match with shower2
00518                               if(objectMatch[i][k]==0) isAssocTrack = true;
00519                           }
00520                         }
00521                         if(!isAssocTrack) { //no track associated
00522                           MSG("EventSS",Msg::kDebug) << "Shower " << j 
00523                                                      << " has no track association. " 
00524                                                      << "Add to event " << nEvents << endl; 
00525                           event->Add(recolist.At(j));
00526                           ph[j] = 0;
00527                           objectUsed[j] = true;
00528                           nUnusedObjects-=1;
00529                         }
00530                       }
00531                     }
00532                     else { //j matches maxPHIndex, so add it to event
00533                       MSG("EventSS",Msg::kDebug) << "Shower " << j 
00534                                                  << " matches orginal track " << maxPHIndex 
00535                                                  << " Add to event " << nEvents << endl; 
00536                       event->Add(recolist.At(j));
00537                       ph[j] = 0;
00538                       objectUsed[j] = true;
00539                       nUnusedObjects-=1;
00540                     }
00541                   }
00542                 }
00543               } //end if for vertex shower
00544               else { //if not a vertex shower, check other associations
00545                 MSG("EventSS",Msg::kDebug) << "Testing whether shower " << i 
00546                                            << " has better associations. " << endl;
00547                 Bool_t otherTrackVertexMatch = false;
00548                 Bool_t otherTrackMatch = false;
00549                 Bool_t otherShowerMatch = false;
00550                 for(int j=0;j<nRecoObjects;j++){
00551                   if(j!=maxPHIndex && objectMatch[i][j]!=0 && 
00552                      !objectUsed[j] && objectMatch[maxPHIndex][j]==0) {
00553                     if(recolist.At(j)->InheritsFrom("CandTrackHandle")) {
00554                       CandTrackHandle *track2 = 
00555                         dynamic_cast<CandTrackHandle*>(recolist.At(j));
00556                       if(TMath::Abs(shower2->GetBegPlane(PlaneView::kU) - 
00557                                     track->GetBegPlane(PlaneView::kU)) < 
00558                          pShwTrkDz/(PITCHINMETRES*2.) || 
00559                          TMath::Abs(shower2->GetBegPlane(PlaneView::kV) -
00560                                     track->GetBegPlane(PlaneView::kV)) < 
00561                          pShwTrkDz/(PITCHINMETRES*2.) )
00562                         otherTrackVertexMatch = true;
00563                       else {
00564                         //first test whether track2 is buried
00565                         //if so then can still associate shower to track
00566                         if(objectMatch[i][j]!=-1) {
00567                           //object match is not flagged as buried
00568                           //perform check:
00569                           Bool_t track2Buried = false;
00570                           Int_t *stats2 = new Int_t[8];
00571                           shower2->BuriedTrack(track2,1000,stats2);                       
00572                           if(stats2[3]>0){
00573                             Int_t nplanes = track2->GetEndPlane() - track2->GetBegPlane() + 1;
00574                             MSG("EventSS",Msg::kDebug) 
00575                               << "BuriedTrack Stats: " << "\n"
00576                               << "Number of track planes = " << nplanes << "\n"
00577                               << "nSharedPlanes = " << stats2[0] << "\n" 
00578                               << "nMissingIsoPlanes = " << stats2[1] << "\n" 
00579                               << "nMissingSharedPlanes = " << stats2[2] << "\n" 
00580                               << "nBuriedPlanes = " << stats2[3] << "\n" 
00581                               << "nPhysBuriedPlanes = " << stats2[4] << "\n" 
00582                               << "nTrkLike = " << stats2[5] << "\n" 
00583                               << "nTrkLikeTopol = " << stats2[6] << "\n" 
00584                               << "nXTalk = " << stats2[7] << endl;
00585                             nplanes -= (stats2[1] + stats2[2]);
00586                             if( (Float_t(stats2[4])/Float_t(nplanes) > pHardBuriedFrac) ||
00587                                 (Float_t(stats2[3]-stats2[5]-stats2[6])/Float_t(nplanes) > 
00588                                  pHardBuriedFrac) ||
00589                                 (stats2[3]>=Int_t(nplanes*pSoftBuriedFrac) && 
00590                                  (TMath::Abs(track2->GetDirCosU())>pDCosUVCut || 
00591                                   TMath::Abs(track2->GetDirCosV())>pDCosUVCut) ) ) {
00592                                track2Buried = true;
00593                               objectMatch[i][j]=-1;
00594                               objectMatch[j][i]=-1;
00595                             }
00596                           }
00597                           if(!track2Buried){
00598                             //test whether track2 matches better than track
00599                             Int_t *stats = new Int_t[8];
00600                             shower2->BuriedTrack(track,1000,stats);
00601                             if( ( stats2[3] - stats2[5] - stats2[6] ) >=
00602                                 ( stats[3]  - stats[5]  - stats[6]  ) && 
00603                                 ( stats2[4] >= stats[4] ) ) {
00604                               otherTrackMatch = true;
00605                             }
00606                             delete [] stats;
00607                           }
00608                           delete [] stats2;
00609                         }
00610                       }
00611                     }
00612                     else if(recolist.At(j)->InheritsFrom("CandShowerSRHandle")) {
00613                       Int_t *stats = new Int_t[8];
00614                       shower2->BuriedTrack(track,1000,stats);
00615                       if(stats[4]==0) { //how convincing is first track/shower match?
00616                         otherShowerMatch = true;
00617                       }
00618                       delete [] stats;
00619                     }
00620                   }
00621                 }
00622                 
00623                 //if there are no other viable matches:
00624                 if(otherTrackVertexMatch || otherShowerMatch || 
00625                    otherTrackMatch) {
00626                   MSG("EventSS",Msg::kDebug) << "Shower " << i << " has better association: "
00627                                              << " otherTrackVertexMatch = " 
00628                                              << otherTrackVertexMatch 
00629                                              << " otherShowerMatch = " 
00630                                              << otherShowerMatch 
00631                                              << " otherTrackMatch = "
00632                                              << otherTrackMatch << endl;
00633                   continue;
00634                 }
00635                 //then add this shower to this track
00636                 MSG("EventSS",Msg::kDebug) << "Shower " << i << " has no better association. "
00637                                            << "Add to event " << nEvents << endl;
00638                 event->Add(recolist.At(i));
00639                 ph[i] = 0;
00640                 objectUsed[i] = true;
00641                 nUnusedObjects-=1;                
00642               } //end if for non-vertex shower
00643             } //end if for CandShowerSRHandle
00644           } //end if for valid association
00645         } //end for loop over all objects
00646       } //end if maxPH object is a track
00647       
00648       else if(recolist.At(maxPHIndex)->InheritsFrom("CandShowerSRHandle")) {
00649         MSG("EventSS",Msg::kDebug) << "Object " << maxPHIndex << " is a shower" << endl;
00650         MSG("EventSS",Msg::kDebug) << "Check whether shower " << maxPHIndex 
00651                                    << " has other associations" << endl;
00652         //maxPH object is a shower
00653         //CandShowerSRHandle *shower = dynamic_cast<CandShowerSRHandle*>(recolist.At(maxPHIndex));
00654         //now add objects clearly associated with largest object:
00655         for(int i=0;i<nRecoObjects;i++){
00656           if(objectMatch[maxPHIndex][i]!=0 && !objectUsed[i]) {   
00657             if(recolist.At(i)->InheritsFrom("CandTrackHandle")) {
00658               MSG("EventSS",Msg::kWarning) << "Arrgh - Logic flaw in eventSR. "
00659                                            << "Seeing tracks in event building, "
00660                                            << "when all tracks should already be "
00661                                            << "assigned!" << endl;
00662             } //end if CandTrackHandle
00663             else if(recolist.At(i)->InheritsFrom("CandShowerSRHandle")) {
00664               //check for two more levels of associated showers
00665               //CandShowerSRHandle *shower2 = dynamic_cast<CandShowerSRHandle*>(recolist.At(i));
00666               MSG("EventSS",Msg::kDebug) << "Adding shower " << i << " to event " 
00667                                          << nEvents << endl;
00668               MSG("EventSS",Msg::kDebug) << "Check whether shower " << i
00669                                          << " has other associations" << endl;
00670               event->Add(recolist.At(i));
00671               ph[i] = 0;
00672               objectUsed[i] = true;
00673               nUnusedObjects-=1;
00674               for(int j=0;j<nRecoObjects;j++){
00675                 if(j!=maxPHIndex && objectMatch[i][j]!=0 && 
00676                    !objectUsed[j] && objectMatch[maxPHIndex][j]==0) {
00677                   if(recolist.At(j)->InheritsFrom("CandTrackHandle")) {
00678                     MSG("EventSS",Msg::kWarning) << "Arrgh again - Logic flaw in eventSR. "
00679                                                  << "Seeing tracks in event building, "
00680                                                  << "when all tracks should already be "
00681                                                  << "assigned!" << endl;
00682                   }
00683                   else if(recolist.At(j)->InheritsFrom("CandShowerSRHandle")) {
00684                     MSG("EventSS",Msg::kDebug) << "Adding shower " << j << " to event " 
00685                                                << nEvents << endl;
00686                     event->Add(recolist.At(j));
00687                     ph[j] = 0;
00688                     objectUsed[j] = true;
00689                     nUnusedObjects-=1;
00690                   }
00691                 }
00692               }
00693             } //end if CandShowerSRHandle
00694           } //end if valid association
00695         } //end loop over all objects
00696       } //end if maxPH object is a shower
00697         //add event to eventlist
00698       eventlist.Add(event);
00699       nEvents+=1;
00700     }
00701     
00702     // remove all empty events from event list
00703 
00704     for (Int_t i=0; i<=eventlist.GetLast(); i++) {
00705       TObjArray *objectlist = dynamic_cast<TObjArray*>(eventlist.At(i));
00706       if (objectlist->GetLast()<0) eventlist.Remove(objectlist);
00707       else if(objectlist->GetLast()==0) { //consider single shower/track events
00708         if(objectlist->At(0)->InheritsFrom("CandShowerSRHandle")) {
00709           MSG("EventSS",Msg::kDebug) << "Have a single candshower event " << endl;
00710           CandShowerSRHandle *showerSR = 
00711             dynamic_cast<CandShowerSRHandle*> (objectlist->At(0));
00712           //if shower seems unphysical, remove it
00713           if(showerSR->IsUnphysical(pShwXTalkFracCut,pShwXTalkDef)) {
00714             eventlist.Remove(objectlist);
00715             MSG("EventSS",Msg::kDebug) << "REMOVING SHOWER-ONLY EVENT" << endl;
00716             continue;
00717           }
00718         }
00719         else if(objectlist->At(0)->InheritsFrom("CandTrackHandle")) {     
00720           CandTrackHandle *track = dynamic_cast <CandTrackHandle*> (objectlist->At(0));
00721           if(track->IsUnphysical(pTrkGapFracCut,pTrkAsymCut,pTrkXTalkFracCut,pTrkXTalkDef)) {
00722             //if track seems unphysical, remove it
00723            
00724             eventlist.Remove(objectlist);
00725             MSG("EventSS",Msg::kDebug) << "REMOVING TRACK-ONLY EVENT" << endl;
00726             continue;
00727           }
00728         }
00729       }
00730       else if(objectlist->GetLast()==1) { //multiple shower/track events
00731         for(int ind_loop = 0;ind_loop<2;ind_loop++){
00732           Int_t ind_loop1 = ind_loop;
00733           Int_t ind_loop2 = ind_loop+1; if(ind_loop2==2) ind_loop2 = 0;
00734           if(objectlist->At(ind_loop1)->InheritsFrom("CandShowerSRHandle")){
00735             CandShowerSRHandle *showerSR = 
00736               dynamic_cast<CandShowerSRHandle*> (objectlist->At(ind_loop1));
00737             if(showerSR->IsUnphysical(pShwXTalkFracCut,pShwXTalkDef)) {
00738               if(objectlist->At(ind_loop2)->InheritsFrom("CandTrackHandle")){
00739                 CandTrackHandle *track = dynamic_cast <CandTrackHandle*> (objectlist->At(ind_loop2));
00740                 if(track->IsUnphysical(pTrkGapFracCut,pTrkAsymCut,pTrkXTalkFracCut,pTrkXTalkDef)) {
00741                   eventlist.Remove(objectlist);
00742                   MSG("EventSS",Msg::kDebug) << "REMOVING UNPHYSICAL SHOWER+TRACK EVENT" << endl;
00743                   continue;
00744                 }
00745               }
00746             }
00747           }
00748         }
00749       }
00750       if(useChopper==1) {
00751         std::vector<CandHandle> evtVec;
00752         for(int i=0;i<=objectlist->GetLast();i++){
00753           CandHandle *event = dynamic_cast<CandHandle*>(objectlist->At(i));
00754           evtVec.push_back(*event);       
00755         }
00756         ChopHelper chopper(cx.GetMom());
00757         ChopHelp *chop = chopper.GetChopHelp(evtVec);           
00758         if(chop->nchop==1) continue;
00759         else {
00760           chop->Print();
00761           //do something with chop info
00762           //for(int i=0;i<ncandidates;i++){
00763           //}
00764         }
00765       }
00766     }
00767     eventlist.Compress();
00768 
00769 
00770     // now make CandEvents for each event we have found
00771     Bool_t buildEvent=false;
00772     for (Int_t i=0; i<=eventlist.GetLast(); i++) {
00773       TObjArray *objectlist = dynamic_cast<TObjArray*>(eventlist.At(i));
00774       MSG("EventSS",Msg::kDebug)  << "making event of " 
00775                                     << objectlist->GetLast()+1 
00776                                     << " objects \n";
00777       cxx.SetDataIn(objectlist);
00778       buildEvent=true;
00779       CandEventHandle eventhandle = CandEvent::MakeCandidate(ah,cxx);
00780       MSG("EventSS",Msg::kDebug) << " # of strips:  " 
00781                                    << eventhandle.GetNDaughters() << "\n";
00782       eventhandle.SetCandSlice(slice);
00783       ch.AddDaughterLink(eventhandle);
00784       delete eventlist.At(i);
00785     }
00786     if (!buildEvent) BuildEventFromUnassoc(ac,ch,cx,slice);
00787   } // end loop over slice list
00788 
00789   delete trackItr;
00790   delete showerItr;
00791   delete subshowerItr;
00792 
00793   MSG("EventSS",Msg::kDebug) << "starting unassociated hits assignment \n";
00794 
00795   // find unassociated hits, and place them in tobjarray
00796   sliceItr.Reset();
00797   TObjArray unassociated;
00798   CandEventListHandle &eventlist = dynamic_cast<CandEventListHandle &>(ch);
00799   CandEventHandleItr * eventItr2 = 0;
00800   if (eventlist.GetNDaughters()>0) {    
00801     eventItr2 = new CandEventHandleItr(eventlist.GetDaughterIterator());
00802   }
00803   FindUnassociated(sliceItr,eventItr2,unassociated);
00804   delete eventItr2;
00805   SetPrimaryShowers(ch,ac);
00806   
00807   // now loop over unnassociated hits, and check for adjacency with all
00808   // showers, and track vertices (in the case that no shower exists).  
00809   // If a match is found in the former case, add strip to shower. 
00810   // In the latter case, create a new shower.
00811  
00812   // this map contains the separation between each event/unassociated hit
00813   vector<map<const CandEventHandle*, Double_t> > dist2map(unassociated.GetLast()+1);
00814   FillDist2Map(ac,unassociated,ch,dist2map);
00815   
00816   // initialize a vector which holds a Bool_t
00817   // indicating whether a hit has already been used
00818   vector<Bool_t> alreadyUsed(unassociated.GetLast()+1);
00819   for (Int_t i=0; i<=unassociated.GetLast(); i++) alreadyUsed[i] = false;  
00820 
00821   // loop over unnassociated hits, adding them to events.
00822   // After each loop the affected distance map entries are recalculated, 
00823   // and the  loop is repeated.  This continues until a loop adds no strips.
00824   Bool_t found=true;
00825   Int_t n_found=0;
00826   while (found) {
00827     
00828     //found indicates whether a strip was added to an event in this loop
00829     found=false;
00830     
00831     CandStripHandle *closeststrip= 0;
00832     CandEventHandle *closestevent= 0;
00833     Double_t closestdist2 = 9999999.;
00834     Int_t closesti=-1;
00835     
00836     // loop over unassoc strips
00837     for (Int_t i=0; i<=unassociated.GetLast(); i++) {
00838       if(alreadyUsed[i]) continue;
00839 
00840       CandStripHandle *strip =
00841         dynamic_cast<CandStripHandle*>(unassociated.At(i));
00842       
00843       CandEventHandle *bestevent = 0;
00844       Double_t bestdist2 = 9999999.;
00845       // loop over events, and find event which is closest to strip
00846       // and in time
00847       TIter eventItr(ch.GetDaughterIterator());
00848       while (CandEventHandle *event =
00849              dynamic_cast<CandEventHandle*>(eventItr())) {
00850         Double_t dtime = strip->GetBegTime()-event->GetVtxT();
00851         if ( dtime>pHitAssocDTime0 && dtime<pHitAssocDTime1 &&
00852              (!bestevent || dist2map[i][event]<bestdist2) ) {
00853           bestevent = event;
00854           bestdist2 = dist2map[i][event];
00855         }
00856       }
00857       if (bestevent && (!closeststrip || bestdist2<closestdist2)) {
00858         closeststrip = strip;
00859         closesti=i;
00860         closestdist2 = bestdist2;
00861         closestevent = bestevent;
00862       }
00863     }
00864     MSG("EventSS",Msg::kVerbose)<< "closest strip " << closestdist2
00865                                 << "/"<<pHitAssocMaxDist2<< endl;
00866 
00867     // if distance to closest event less than maximum allowable add to event   
00868     if(closeststrip && closestevent && closestdist2<pHitAssocMaxDist2){
00869       MSG("EventSS",Msg::kVerbose) << "closest strip "
00870                                    << closeststrip->GetPlane() << " "
00871                                    << closeststrip->GetTPos()
00872                                    << " " << closestdist2 << "\n";       
00873       alreadyUsed[closesti] = true;
00874       found = true;
00875       ++n_found;
00876       
00877       // now either add this strip to existing primary shower, or
00878       // create a new shower at the track vertex if one does not exist
00879       // Get PrimaryShower CandHandle from fShowerList for modification.
00880       CandShowerHandle *shower = 0;
00881       if (CandShowerHandle *primsh = closestevent->GetPrimaryShower()) {
00882         for (Int_t ishower=0; 
00883              !shower && ishower<closestevent->GetLastShower()+1; 
00884              ishower++) {
00885           const CandShowerHandle *target = closestevent->GetShower(ishower);
00886           if (primsh->IsCloneOf(*target))     // Tests clone or Cand ==
00887             shower = closestevent->GetShowerWritable(ishower);
00888         }
00889       }
00890 
00891       Float_t minShwPlane=0;
00892       Float_t maxShwPlane=0;
00893       if(shower){
00894         minShwPlane=min(shower->GetVtxPlane(),shower->GetEndPlane());
00895         maxShwPlane=max(shower->GetVtxPlane(),shower->GetEndPlane());
00896       }
00897       Float_t vertexsep = 0;
00898       CandTrackHandle *track = closestevent->GetPrimaryTrack();
00899       if (shower && track) {
00900         vertexsep = abs(shower->GetBegPlane() - track->GetBegPlane());
00901         vertexsep *= PITCHINMETRES;
00902         MSG("EventSS",Msg::kVerbose) << "Found track and primary shower. \n" 
00903                                      << "Vertex separation (limit) = " 
00904                                      << vertexsep << "(" 
00905                                      << pShwTrkDz << ")" << " \n";
00906       }
00907       // if event has shower and either no track, or shower is close
00908       // to track vertex, add this hit to the shower which was
00909       // compared in position to this strip earler.
00910       if (shower && (!track || vertexsep<pShwTrkDz)) {
00911         // NOTE: CandShowerHandle *shower, used to modify CandShower, is owned
00912         // by fShowerList.
00913         // add new strip to this shower, and the appropriate cluster
00914         MSG("EventSS",Msg::kVerbose) << "Shower min/max plane (limit): "
00915                                      << minShwPlane << "/" << maxShwPlane
00916                                      << " (" << pMaxNewShwLen << ")" << "\n";
00917         if(closeststrip->GetPlane()>minShwPlane-pMaxNewShwLen &&
00918            closeststrip->GetPlane()<maxShwPlane+pMaxNewShwLen ){
00919           MSG("EventSS",Msg::kVerbose) << "adding strip to shower \n"; 
00920           AddStripToEvent(closestevent,shower,subshowerlist,closeststrip);
00921         }
00922       }
00923       // no shower consistent with track vertex, and this strip
00924       // is close to vertex, so start new shower.    
00925       else if (track) CreatePrimaryShower(ac,cxx,closestevent,showerlist,
00926                                           subshowerlist,closeststrip);      
00927       
00928       // need to recalculate distances bewteen remaining
00929       // unassoc hits and this modified event, if the added strip
00930       // results in a smaller separation, replace value in dist2map 
00931       ReFillDist2Map(ac,closeststrip,closestevent,unassociated,dist2map);      
00932     } // if adding strip to event in form of new shower
00933   } // if found (if an unassoc hit was matched to event in last loop over them
00934 
00935   // because we have added strips to previously constructed showers, their
00936   // runalg method should be run on them again.  In the future, an
00937   // optimization would be to do this only on showers which have been
00938   // touched.
00939   MSG("EventSS",Msg::kVerbose) << "added " << n_found << " of "
00940                                << unassociated.GetLast()+1 
00941                                << " unassociated hits " << endl;
00942   
00943   cxx.SetDataIn(showerlist);  // Pass along persistable CandShowerList
00944   ReConstructShowers(ac,ch,cxx);
00945 
00946   //Final step: use the chopper to combine events across slices:
00947   if(useChopper==1) {
00948     std::vector<CandHandle> evtVec;
00949     CandEventHandleItr evitr(ch.GetDaughterIterator());
00950     while(CandEventHandle *event = dynamic_cast<CandEventHandle*>(evitr())){
00951       evtVec.push_back(*event);
00952     }
00953     ChopHelper chopper(cx.GetMom());
00954     ChopHelp *chop = chopper.GetChopHelp(evtVec);
00955     chop->Print();
00956   }
00957 
00958   // For each event, set primary shower based on 
00959   // distance from vertex and energy
00960   SetPrimaryShowers(ch,ac);
00961 }

void AlgEventSSList::SetPrimaryShowers CandHandle ch,
AlgConfig ac
 

Redetermine primary shower for each event. If the largest shower is greater than shwshwdz from the track vertex, and the closest shower has energy greater than minshwEfract call the closest shower the primary, otherwise the largest shower is the primary

Definition at line 1370 of file AlgEventSSList.cxx.

References CandHandle::GetDaughterIterator(), Registry::GetDouble(), and CandEventHandle::SetPrimaryShower().

Referenced by RunAlg().

01371 {
01372   
01373   // This method loops over events in the eventlist, finding the largest
01374   // shower and setting that to the primary shower in the event.
01375   
01376   Double_t pShwShwDz = ac.GetDouble("ShwShwDz");
01377   Double_t pMinShwEFract = ac.GetDouble("MinShwEFract");
01378   
01379   TIter evItr(ch.GetDaughterIterator());
01380   while (CandEventHandle *ev =
01381                               dynamic_cast<CandEventHandle*>(evItr())) {
01382     ev->SetPrimaryShower(pMinShwEFract,pShwShwDz);
01383   }
01384 }

void AlgEventSSList::Trace const char *  c  )  const [virtual]
 

Reimplemented from AlgBase.

Definition at line 2012 of file AlgEventSSList.cxx.

02013 {
02014 }


The documentation for this class was generated from the following files:
Generated on Fri Mar 28 15:53:05 2008 for loon by  doxygen 1.3.9.1