AlgFarDetEventList.cxx

Go to the documentation of this file.

#include <cstdlib>

#include "AlgFarDetEventList.h"
#include "FarDetEventHandle.h"
#include "FarDetEventListHandle.h"

#include "Algorithm/AlgFactory.h"
#include "Algorithm/AlgHandle.h"
#include "Algorithm/AlgConfig.h"

#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"

#include "RecoBase/CandFitTrackListHandle.h"
#include "RecoBase/CandFitTrackHandle.h"
#include "RecoBase/CandTrackListHandle.h"
#include "RecoBase/CandTrackHandle.h"
#include "RecoBase/CandShowerListHandle.h"
#include "RecoBase/CandShowerHandle.h"
#include "RecoBase/CandStripListHandle.h"
#include "RecoBase/CandStripHandle.h"

#include "Candidate/CandContext.h"

ClassImp(AlgFarDetEventList)    

CVSID("$Id: AlgFarDetEventList.cxx,v 1.4 2009/02/28 21:46:10 gmieg Exp $");

AlgFarDetEventList::AlgFarDetEventList() 
{

}   

AlgFarDetEventList::~AlgFarDetEventList()
{

}
   
void AlgFarDetEventList::RunAlg(AlgConfig& /*ac*/, CandHandle &ch, CandContext &cx)
{
  MSG("FarDetEvent",Msg::kDebug) << " AlgFarDetEventList::RunAlg(...) " << endl;

  FarDetEventListHandle& fdeventlist = dynamic_cast<FarDetEventListHandle&>(ch);

  const TObjArray* arr = dynamic_cast<const TObjArray*>(cx.GetDataIn());

  AlgFactory &af = AlgFactory::GetInstance();
  AlgHandle ah = af.GetAlgHandle("AlgFarDetEvent","default");
  CandRecord* candrec = (CandRecord*)(cx.GetCandRecord());
  CandContext cxnew(this, cx.GetMom());
  cxnew.SetCandRecord(candrec);

  // Create TObjArray to store the event information
  TObjArray* eventinfo = new TObjArray();
  
  // Start to find a new event...
  // N.B: Assume there is only one event for now!
  eventinfo->Clear();
    
  // Get the track information
  MSG("FarDetEvent",Msg::kVerbose) << " GET TRACK INFORMATION " << endl;
  CandFitTrackListHandle* fitlist = (CandFitTrackListHandle*)(arr->At(0));
  CandFitTrackHandle* bestFitTrack=0;

  CandTrackListHandle* trklist = (CandTrackListHandle*)(arr->At(1));
  CandTrackHandle* bestTrack=0;

  // Look in the fitted tracks first
  if( fitlist )
    {
      Double_t bestFitTracklength=-1.0;
      TIter fititr(fitlist->GetDaughterIterator());
      while(CandFitTrackHandle* fit = (CandFitTrackHandle*)(fititr()))
        {
          if(fit) 
            {
              if(fit->GetRange()>bestFitTracklength) 
                {
                  bestFitTrack = fit;
                  bestFitTracklength = bestFitTrack->GetRange();
                }
            }
        }
      
    }

  // Otherwise fall back to the found tracks.
  else if( trklist )
    {
      Double_t bestTracklength=-1.0;
      TIter trkitr(trklist->GetDaughterIterator());
      while(CandTrackHandle* trk = (CandTrackHandle*)(trkitr()))
        {
          if(trk) 
            {
              if(trk->GetRange()>bestTracklength) 
                {
                  bestTrack = trk;
                  bestTracklength = bestTrack->GetRange();
                }
            }
        }
    }

  eventinfo->Add(bestFitTrack); 
  eventinfo->Add(bestTrack); // Add bestTrack even if it is null.

  CandTrackHandle* primaryTrack = 0;
  if( bestFitTrack ) primaryTrack=(CandTrackHandle*)bestFitTrack;
  else if( bestTrack ) primaryTrack=(CandTrackHandle*)bestTrack;

  if( primaryTrack ){
    MSG("FarDetEvent",Msg::kVerbose) << "   ...found a primary track " << endl;
  }
  else{
    MSG("FarDetEvent",Msg::kVerbose) << "   ...there is no primary track " << endl;
  }

  // Now we've stored the track we want to use in our event, 
  // we should try and find a shower to go with it.
  MSG("FarDetEvent",Msg::kVerbose) << " GET SHOWER INFORMATION " << endl;
  CandShowerListHandle* shwlist = (CandShowerListHandle*)(arr->At(2));

  CandShowerHandle* vtxshower = 0;
  TObjArray* HadronicShowers = new TObjArray();
  TObjArray* BremsAndDeltas = new TObjArray();
  
  if( shwlist )
    { // Loop over the shower list and look for a vertex shower.
      
      Double_t separation(999.),temp(999.),shwenergy(0.0);
      if( primaryTrack )
        { // We have a track so find the shower closest to it
          this->CompareShowersWithTrack(primaryTrack, 
                                        shwlist, 
                                        vtxshower, 
                                        HadronicShowers, 
                                        BremsAndDeltas, 
                                        separation);
        }
      
      else
        { // We have no track so find the largest shower
          TIter shwitr(shwlist->GetDaughterIterator());
          while(CandShowerHandle* shw = (CandShowerHandle*)(shwitr()))
            { 
              if( shw )
                {
                  temp=shw->GetEnergy();
                  if(temp>shwenergy)
                    {
                      vtxshower=shw;
                      shwenergy=temp;
                    }
                }       
            }  
        }
    }
  
  if( vtxshower ){
    MSG("FarDetEvent",Msg::kVerbose) << "   ...found a primary shower " << endl;
  }
  else{
    MSG("FarDetEvent",Msg::kVerbose) << "   ...there is no primary shower " << endl;
  }

  eventinfo->Add(vtxshower);
  eventinfo->Add(HadronicShowers);
  eventinfo->Add(BremsAndDeltas);

  // Add the strip list
  MSG("FarDetEvent",Msg::kVerbose) << " GET STRIP INFORMATION " << endl;
  CandStripListHandle* strplist = (CandStripListHandle*)(arr->At(3));

  TObjArray* AssociatedStrips = new TObjArray();
  TObjArray* AllStrips = new TObjArray();
  
  if( strplist )
    {
      TIter strpitr(strplist->GetDaughterIterator());
      while(CandStripHandle* strp = (CandStripHandle*)(strpitr()))
        { 
          if( strp )
            {

              // record all strips
              AllStrips->Add(strp);

              // some association could go here...
              AssociatedStrips->Add(strp);

            }
        }
    }

  eventinfo->Add(AssociatedStrips);
  eventinfo->Add(AllStrips);
  
  // Add event info to CandContext
  cxnew.SetDataIn(eventinfo);
    
  // Make new FarDetEventHandle
  MSG("FarDetEvent",Msg::kVerbose) << " MAKING FARDETEVENT " << endl;
  FarDetEventHandle fdevent = FarDetEvent::MakeCandidate(ah,cxnew);
  fdevent.SetName("FarDetEvent");
  fdevent.SetTitle(TString("Created by AlgFarDetEventList"));
  fdeventlist.AddDaughterLink(fdevent);

  // Delete temporary objects
  HadronicShowers->Clear();
  delete HadronicShowers;

  BremsAndDeltas->Clear();
  delete BremsAndDeltas;

  AssociatedStrips->Clear();
  delete AssociatedStrips;

  AllStrips->Clear();
  delete AllStrips;

  //
  // Could turn this into a loop and make another event now...
  // 

  eventinfo->Clear();
  delete eventinfo;
}
   
void AlgFarDetEventList::Trace(const char* /*c*/ ) const
{

}

void AlgFarDetEventList::CompareShowersWithTrack(CandTrackHandle* bestTrack, CandShowerListHandle* shwlist, CandShowerHandle*& vtxshower, TObjArray* HadronicShowers, TObjArray* BremsAndDeltas, Double_t& separation)
{
  MSG("FarDetEvent",Msg::kVerbose) << " CompareShowersWithTrack(...) " << endl;

  if(!bestTrack)
    {
      MSG("FarDetEvent",Msg::kError) << " AlgFarDetEventList::CompareShowersWithTrack * Not Given A Track To Compare With Showers!" << endl;
      return;
    }
  if(!shwlist)
    {
      MSG("FarDetEvent",Msg::kError) << " AlgFarDetEventList::CompareShowersWithTrack * Not Given A List of Showers To Compare With Track!" << endl;
      return;
    }
  if(shwlist->GetNDaughters()<=0)
    {
      MSG("FarDetEvent",Msg::kWarning) << " AlgFarDetEventList::CompareShowersWithTrack * Given An Empty Showerlist To Compare With Track!" << endl;
      return;
    }

  HadronicShowers->Clear(); BremsAndDeltas->Clear();

  int vtxpln = bestTrack->GetVtxPlane();
  double trku = bestTrack->GetVtxU();
  double trkv = bestTrack->GetVtxV();
  double trkz = bestTrack->GetVtxZ();

  TIter shwitr(shwlist->GetDaughterIterator());
  while(CandShowerHandle* shw = (CandShowerHandle*)(shwitr()))
    { 
      if( shw )
        {
          // Look for possible hadronic showers
          // (offset by less than 10 planes from track vertex)
          if( abs(shw->GetVtxPlane()-vtxpln)<20 )
            {
              HadronicShowers->Add(shw);
              double temp = (shw->GetVtxU()-trku)*(shw->GetVtxU()-trku) 
                          + (shw->GetVtxV()-trkv)*(shw->GetVtxV()-trkv) 
                          + (shw->GetVtxZ()-trkz)*(shw->GetVtxZ()-trkz);
              if(temp<separation)
                {
                  vtxshower=shw;
                  separation=temp;
                }
            }
          // Look for possible brems and deltas
          // ( loop over the planes in the shower and see if 
          //   it is associated with the track in that plane)
          else
            {

            }
        }
    }
  return;
}

//
// void AlgFarDetEventList:::RemoveTracksinShowers(AlgConfig &ac, TObjArray * recolist)
// {
//
//   Double_t pMinTrackIsolation = ac.GetDouble("MinTrackIsolation");
//   Double_t pMinTrackIsolationDist = ac.GetDouble("MinTrackIsolationDist");
//
//   for (Int_t i = 0; i <= recolist->GetLast(); i++) {
//     CandTrackHandle *track=0;
//     Int_t nisolated = 0;
//     Int_t nisolatedU = 0;
//     Int_t nisolatedV = 0;
//     Bool_t isolated = true;
//     Int_t longestcontig = 0;
//     Int_t longestcontigU = 0;
//     Int_t longestcontigV = 0;
//
//     if (recolist->At(i)->InheritsFrom("CandTrackHandle")) {
//       track = dynamic_cast<CandTrackHandle*>(recolist->At(i));
//    
//       // count number of planes in each view in which
//       // track is outside shower
//       Int_t firstplane = TMath::Min(track->GetBegPlane(),
//                                  track->GetEndPlane());
//       Int_t lastplane = TMath::Max(track->GetBegPlane(),
//                                 track->GetEndPlane());
//       Int_t firstuplane = TMath::Min(track->GetBegPlane(PlaneView::kU),
//                                   track->GetEndPlane(PlaneView::kU));
//       Int_t firstvplane = TMath::Min(track->GetBegPlane(PlaneView::kV),
//                                   track->GetEndPlane(PlaneView::kV));
//       Int_t lastuplane = TMath::Max(track->GetBegPlane(PlaneView::kU),
//                                track->GetEndPlane(PlaneView::kU) );
//       Int_t lastvplane = TMath::Max(track->GetBegPlane(PlaneView::kV),
//                                track->GetEndPlane(PlaneView::kV) );      
//
//       MSG("FarDetEvent",Msg::kDebug) << " track extent " << firstplane << " " << lastplane <<  " " << track->GetNDaughters() << " " << track->GetBegPlane() << " " << track->GetEndPlane() << endl;
//
//       for(Int_t iplane = firstplane;iplane <= lastplane;iplane++){
//      PlexPlaneId plnid(track->GetVldContext()->GetDetector(),iplane,false);
//      // require that the plane be within the track 
//      // extent for this planeview
//      isolated = true;
//      if((plnid.GetPlaneView() == PlaneView::kU  && 
//          iplane >= firstuplane && iplane <= lastuplane) ||
//         (plnid.GetPlaneView() == PlaneView::kV  && 
//          iplane >= firstvplane && 
//          iplane<=lastvplane)){
//      
//        for (Int_t j = 0; j<=recolist->GetLast(); j++) {
//          if(isolated){
//            if (recolist->At(j)->InheritsFrom("CandShowerHandle")) {
//              CandShowerHandle * shower = dynamic_cast<CandShowerHandle*>(recolist->At(j));
//              if(iplane >= shower->GetBegPlane() &&
//                 iplane <= shower->GetEndPlane()){
//                // if plane within shower extent, compare track position
//                // with shower width. If track outside shower by
//                // a distance pMinTrackIsolationDist, or shower width
//                // is less than this distance, inc isolated.
//                Double_t minPE=2.0;
//                Float_t minUshw = shower->GetMinU(iplane,minPE)+0.02;
//                Float_t minVshw = shower->GetMinV(iplane,minPE)+0.02;
//                Float_t maxUshw = shower->GetMaxU(iplane,minPE)-0.02;
//                Float_t maxVshw = shower->GetMaxV(iplane,minPE)-0.02;
//
//                MSG("FarDetEvent",Msg::kDebug) << " shower extent Plane: " << iplane << " U " << minUshw << "/" <<  track->GetU(iplane) << "/" << maxUshw << " V  " << minVshw << "/" << track->GetV(iplane) << "/" << maxVshw << endl;
//                
//                if(plnid.GetPlaneView() == PlaneView::kU && 
//                   (track->GetU(iplane) >= (minUshw-pMinTrackIsolationDist) && 
//                    track->GetU(iplane) <= (maxUshw+pMinTrackIsolationDist) &&
//                    maxUshw-minUshw > (pMinTrackIsolationDist))){
//                  isolated = false;
//                }
//                if(plnid.GetPlaneView() == PlaneView::kV &&
//                   (track->GetV(iplane) >= (minVshw-pMinTrackIsolationDist)&& 
//                    track->GetV(iplane) <= (maxVshw+pMinTrackIsolationDist)&&
//                    maxVshw-minVshw > (pMinTrackIsolationDist))){
//                  isolated = false;
//                }
//                if(!track->IsTPosValid(iplane))isolated = false;
//              }
//            }      
//          }
//        }
//        if(track->IsTPosValid(iplane)){
//          nisolated++;
//          if(plnid.GetPlaneView() == PlaneView::kU){
//            nisolatedU++;
//          }
//          else{
//            nisolatedV++;
//          }
//        }
//        if(!isolated) {
//          nisolated=0;  // reset contig. isolation counter
//          if(plnid.GetPlaneView() == PlaneView::kU){
//            nisolatedU=0;
//          }
//          else{
//            nisolatedV=0;
//          }
//        }
//        if(nisolated > longestcontig)longestcontig = nisolated;
//        if(nisolatedU > longestcontigU)longestcontigU = nisolatedU;
//        if(nisolatedV > longestcontigV)longestcontigV = nisolatedV;
//         
//        MSG("FarDetEvent",Msg::kDebug) << " nisolated, longest " << nisolated << " " << longestcontig << endl;
//
//        // if the number of contiguous isolated track planes is less than threshold
//        // we remove this track from the event         
//      }
//       } 
//     }
//     MSG("FarDetEvent",Msg::kDebug) << " longest contig:" << longestcontig  << " U " << longestcontigU << " V " << longestcontigV<< "/" << pMinTrackIsolation << endl;
//     if(track && longestcontig < pMinTrackIsolation &&
//        longestcontigU < pMinTrackIsolation &&
//        longestcontigV < pMinTrackIsolation){
//       MSG("FarDetEvent",Msg::kDebug) << " REMOVING TRACK" << endl;
//       recolist->RemoveAt(i);     
//       recolist->Compress();
//       //If you compress here then decrement i so that you 
//       //don't miss any tracks in recolist
//       i-=1;
//     }
//   }
// }
//
//

---