AlgDeMuxGolden.cxx

Go to the documentation of this file.

// $Id: AlgDeMuxGolden.cxx,v 1.15 2003/05/14 19:28:02 brebel Exp $
//
// AlgDeMuxGolden
//
// An Algorithm class that DeMultiplexes a cosmic ray event and adjusts the
//weights in the PlexSEIdAltLists of the CandDeMuxDigits from that event
//
// Author:  B. Rebel 11/2000

#include <cassert>

#include "Candidate/CandContext.h"
#include "Candidate/CandHandle.h"
#include "CandDigit/CandDeMuxDigitHandle.h"
#include "CandDigit/CandDeMuxDigitListHandle.h"
#include "MessageService/MsgService.h"
#include "Algorithm/AlgConfig.h"
#include "DeMux/DmxStatus.h"
#include "DeMux/DmxUtilities.h"
#include "DeMux/DmxHypothesis.h"
#include "DeMux/DmxPlane.h"
#include "DeMux/DmxShowerPlane.h"
#include "DeMux/DmxMuonPlane.h"
#include "Algorithm/AlgFactory.h"
#include "Algorithm/AlgHandle.h"
#include "DeMux/AlgDeMuxGolden.h"
#include "Navigation/NavKey.h"
#include "Navigation/NavSet.h"
#include "TObjArray.h"
#include "TH1.h"
#include "TFolder.h"
#include "TROOT.h"


//define the NavKey to sort the planes
static NavKey KeyPlane( const DmxPlane *currentPlane){
  return currentPlane->GetPlaneNumber();
}

//define the NavKeys to select U and V planes
NavKey KeyU( const DmxPlane *currentPlane){
  return currentPlane->GetPlaneView() == PlaneView::kU;
}

NavKey KeyV( const DmxPlane *currentPlane){
  return currentPlane->GetPlaneView() == PlaneView::kV;
}

//define the NavKeys to select only golden U and V planes
NavKey KeyOnGoldenU( const DmxPlane *currentPlane){
  return currentPlane->IsGolden() && currentPlane->GetPlaneView() == PlaneView::kU;
}

NavKey KeyOnGoldenV( const DmxPlane *currentPlane){
  return currentPlane->IsGolden() && currentPlane->GetPlaneView() == PlaneView::kV;
}

//define the NavKeys to select only lead U and V planes
NavKey KeyOnLeadU( const DmxPlane *currentPlane){
  return !currentPlane->IsGolden() && currentPlane->GetPlaneView() == PlaneView::kU;
}

NavKey KeyOnLeadV( const DmxPlane *currentPlane){
  return !currentPlane->IsGolden() && currentPlane->GetPlaneView() == PlaneView::kV;
}

ClassImp(AlgDeMuxGolden)
  

//-----------------------------------------------------------------------

CVSID("$Id: AlgDeMuxGolden.cxx,v 1.15 2003/05/14 19:28:02 brebel Exp $");

//-----------------------------------------------------------------------

AlgDeMuxGolden::AlgDeMuxGolden() :
  fStrayPlanes(0),
  fPlanesInSet(6),
  fStrayCut(6)
{
  //default constructor
}

//-----------------------------------------------------------------------

AlgDeMuxGolden::~AlgDeMuxGolden()
{

}

//-----------------------------------------------------------------------
//run the demuxing on an event
void AlgDeMuxGolden::RunAlg(AlgConfig &acd, CandHandle &ch, CandContext & /* cx */)
{

  assert( ch.InheritsFrom("CandDeMuxDigitListHandle") );
  CandDeMuxDigitListHandle &cdlh = dynamic_cast<CandDeMuxDigitListHandle&>(ch);
  //get the AlgConfigDeMux object
  fPlanesInSet = acd.GetInt("PlanesInSet");
  //MSG("DMX", Msg::kInfo) << "window size = " << fPlanesInSet << "\tcut = " << fStrayCut << endl;

  //get the DmxStatus object
  // Find the DmxStatus object or create one for needed scratch space
  DmxStatus *status = dynamic_cast<DmxStatus *>(gROOT->GetRootFolder()->FindObject("Loon/DeMux/DmxStatus"));
  bool tempStatus = false;
  if (status == 0) {
    MSG("DMXX", Msg::kDebug) << "//root/Loon/DeMux/DmxStatus not found."
                              << " Create a temporary DmxStatus." << endl;
    status = new DmxStatus;      // Make a temporary DmxStatus if needed
    tempStatus = true;
  }
  else status->ResetStatus(); 
      
  //instantiate the DmxUtilities object and then fill DmxStatus with the event information      
  DmxUtilities util;
  
  util.FillPlaneArray(status, cdlh, acd);
  const TObjArray *planeArray = status->GetPlaneArray();
  //DmxPlane *plane = dynamic_cast<DmxPlane *>(planeArray->First());
  //MSG("DMX", Msg::kInfo) <<"Event = " <<  status->GetEventNumber() << endl;
  
  //create the DmxPlaneItr over planes and program it to sort the planes
  DmxPlaneItr planeItr(planeArray);
  DmxPlaneKeyFunc *pnKF = planeItr.CreateKeyFunc();
  pnKF->SetFun(KeyPlane);
  planeItr.GetSet()->AdoptSortKeyFunc(pnKF);
  pnKF = 0;
  
  planeItr.ResetFirst();
  status->SetNumberOfPlanes(planeItr.SizeSelect());
  status->SetVertexPlaneNumber(util.FindVertexPlane(planeItr));
  
  //demux the event if there is an identified vertex
  if( status->GetVertexPlaneNumber() != -1){
    
    planeItr.ResetFirst();
    planeItr.GetSet()->Slice(status->GetVertexPlaneNumber(), 500);
    //set the end plane number since there was a vertex plane
    status->SetEndPlaneNumber(util.FindEndPlane(planeItr));
    planeItr.GetSet()->Slice();
    planeItr.ResetFirst();

    //set the z position of the vertex plane
    planeItr.GetSet()->Slice(status->GetVertexPlaneNumber());
    planeItr.ResetFirst();
    status->SetVertexPlaneZPosition(dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetZPosition());

    planeItr.GetSet()->Slice();
    planeItr.ResetFirst();

    //DmxPlane *plane = dynamic_cast<DmxPlane *>(planeArray->First());
    //MSG("DMX", Msg::kInfo) <<"Event = " <<  status->GetEventNumber() << endl;
    
    //create two DmxPlaneItrs
    DmxPlaneItr leadPlaneItr(planeArray);
    DmxPlaneItr vertexPlaneItr(planeArray);
    DmxPlaneItr goldenPlaneItr(planeArray);
    
    //create a KeyFunc to sort planes by number
    DmxPlaneKeyFunc *pnKF = leadPlaneItr.CreateKeyFunc();
    DmxPlaneKeyFunc *vertexPlaneKF = vertexPlaneItr.CreateKeyFunc();
    DmxPlaneKeyFunc *goldenPlaneKF = goldenPlaneItr.CreateKeyFunc();
    
    //program the KeyFunc with the sort function
    pnKF->SetFun(KeyPlane);
    vertexPlaneKF->SetFun(KeyPlane);
    goldenPlaneKF->SetFun(KeyPlane);
    
    //get the NavSet from the iterator and pass the KeyFunc to it
    leadPlaneItr.GetSet()->AdoptSortKeyFunc(pnKF);
    vertexPlaneItr.GetSet()->AdoptSortKeyFunc(vertexPlaneKF);
    goldenPlaneItr.GetSet()->AdoptSortKeyFunc(goldenPlaneKF);
    
    //clear the KF pointer because we no longer own the KeyFunc
    
    pnKF = 0;
    vertexPlaneKF = 0;
    goldenPlaneKF = 0;
    leadPlaneItr.ResetFirst();
    goldenPlaneItr.ResetFirst();
    vertexPlaneItr.ResetFirst();

    //now program a key function to select on plane orientation - U first
    DmxPlaneKeyFunc *orientGoldenUKF = goldenPlaneItr.CreateKeyFunc();
    DmxPlaneKeyFunc *orientLeadUKF = leadPlaneItr.CreateKeyFunc();
    DmxPlaneKeyFunc *orientGoldenVKF = goldenPlaneItr.CreateKeyFunc();
    DmxPlaneKeyFunc *orientLeadVKF = leadPlaneItr.CreateKeyFunc();
    
    Int_t viewCtr = 0;
    while( viewCtr < 2 ){
     
      if( viewCtr == 0){
        //program this function with the select function
        orientLeadUKF->SetFun(KeyOnLeadU);
        orientGoldenUKF->SetFun(KeyOnGoldenU);
                
        //adopt it as a selection function
        leadPlaneItr.GetSet()->AdoptSelectKeyFunc(orientLeadUKF);
        orientLeadUKF = 0;
        leadPlaneItr.Reset();
        goldenPlaneItr.GetSet()->AdoptSelectKeyFunc(orientGoldenUKF);
        orientGoldenUKF = 0;
        goldenPlaneItr.Reset();
      }
      else if(viewCtr == 1){
        //program this function with the select function
        orientLeadVKF->SetFun(KeyOnLeadV);
        orientGoldenVKF->SetFun(KeyOnGoldenV);
                
        //adopt it as a selection function
        leadPlaneItr.GetSet()->AdoptSelectKeyFunc(orientLeadVKF);
        orientLeadVKF = 0;
        leadPlaneItr.Reset();
        goldenPlaneItr.GetSet()->AdoptSelectKeyFunc(orientGoldenVKF);
        orientGoldenVKF = 0;
        goldenPlaneItr.Reset();
      }

      //slice to those planes from the vertex on
      goldenPlaneItr.GetSet()->Slice(status->GetVertexPlaneNumber(), status->GetEndPlaneNumber());

      if(goldenPlaneItr.SizeSelect() >=2 ){

        //MSG("DMX", Msg::kInfo) << "event " << status->GetEventNumber() << endl;
        //loop over the golden planes and set them
        while(goldenPlaneItr.IsValid() ){

          //MSG("DMX", Msg::kInfo) << dynamic_cast<DmxPlane*>(goldenPlaneItr.Ptr())->GetPlaneNumber() << endl;
          if(dynamic_cast<DmxPlane*>(goldenPlaneItr.Ptr())->GetPlaneType() == DmxPlaneTypes::kShower){
            dynamic_cast<DmxShowerPlane *>(goldenPlaneItr.Ptr())->SetStrips("best");
          }
          goldenPlaneItr.Next();
        }//end setting golden planes

        goldenPlaneItr.Reset();

        UseGoldenFit(leadPlaneItr, goldenPlaneItr, status->GetVertexPlaneNumber(), 
                     status->GetEndPlaneNumber());

        //check the fit
        fStrayPlanes -= util.CheckFit(leadPlaneItr);

        //what about multiple muons?  look for another vertex after the end plane
        leadPlaneItr.GetSet()->Slice();
        goldenPlaneItr.GetSet()->Slice();
        vertexPlaneItr.GetSet()->Slice(status->GetEndPlaneNumber()+1, 500);
        //MSG("DMX1", Msg::kInfo)<< "event =  " << status->GetEventNumber() 
        //                   << "\tend plane = " << status->GetEndPlaneNumber() << endl;
        
        Int_t nextVertex = -1;
        if(vertexPlaneItr.SizeSelect() > 3) nextVertex = util.FindVertexPlane(vertexPlaneItr);
        
        while( nextVertex != -1){
          
          status->SetMultipleMuon(true); 
          vertexPlaneItr.GetSet()->Slice();
          vertexPlaneItr.GetSet()->Slice(nextVertex, 500);
          
          //get the next end plane
          Int_t endPlane = util.FindEndPlane(vertexPlaneItr);

          //MSG("DMX1", Msg::kInfo)<< "\tin double muon loop, vertex at plane " << nextVertex 
          //             << "\tend plane = " << endPlane << "\tplanes in slice = " 
          //             << vertexPlaneItr.SizeSelect() << endl;
          
          goldenPlaneItr.GetSet()->Slice(nextVertex, endPlane);
          
          if(goldenPlaneItr.SizeSelect() >= 2 ){
            //loop over the golden planes and set them
            while(goldenPlaneItr.IsValid() ){
              
              if(dynamic_cast<DmxPlane*>(goldenPlaneItr.Ptr())->GetPlaneType() == DmxPlaneTypes::kShower){
                dynamic_cast<DmxShowerPlane *>(goldenPlaneItr.Ptr())->SetStrips("best");
              }
              goldenPlaneItr.Next();
            }
            
            goldenPlaneItr.Reset();
            
            UseGoldenFit(leadPlaneItr, goldenPlaneItr, status->GetVertexPlaneNumber(), 
                         status->GetEndPlaneNumber());
            
            //check the fit
            fStrayPlanes -= util.CheckFit(leadPlaneItr);
            
            //CheckFitWithTiming(ch, planeItr);
            
          }//end if enough golden planes to try
          
          //look for another vertex after the end plane
          leadPlaneItr.GetSet()->Slice();
          goldenPlaneItr.GetSet()->Slice();
          vertexPlaneItr.GetSet()->Slice();
          vertexPlaneItr.GetSet()->Slice(endPlane+1, 500);
          //MSG("DMX1", Msg::kInfo)<< "\tlook for new vertex, slice from " << endPlane+1 
          //             << "\tto plane 500, planes in slice = " 
          //             << vertexPlaneItr.SizeSelect() << endl;
          
          vertexPlaneItr.Reset();
          
          if(vertexPlaneItr.SizeSelect() > 0){
            nextVertex = util.FindVertexPlane(vertexPlaneItr);
          }
          else{ nextVertex = -1; }
        }//end loop over new verticies
      }//end if at least 2 golden planes
      else{ 
        MSG("DMX", Msg::kWarning)<< "Event " << status->GetEventNumber() 
                                 << " not demuxed; not enough golden planes in view" << endl;
        status->SetValidPlanesFailure(true); 
      } //not enough planes to demux

      status->SetFigureOfMerit(goldenPlaneItr.SizeSelect(), fStrayPlanes);
      if(viewCtr == 0 && status->GetEventDeMuxed()){ 
        status->SetUStrayPlanes(fStrayPlanes);
        status->SetUValidPlanes(goldenPlaneItr.SizeSelect());
        if( status->GetFigureOfMeritFailure() )
        //see if the event failed the figure of Merit test, if so, see if it is an 
        //overlapping muon
        status->SetUOverlappingMultiple(util.IsOverlappingMultiple(goldenPlaneItr, status->GetVertexZPosition(), 
                                                              acd.GetDouble("HoughInterceptRMSLimit"), 
                                                              acd.GetDouble("HoughSlopeRMSLimit"),
                                                              acd.GetDouble("HoughPeakLimit")));
      }
      if(viewCtr == 1 && status->GetEventDeMuxed()){ 
        status->SetVStrayPlanes(fStrayPlanes);
        status->SetVValidPlanes(goldenPlaneItr.SizeSelect());
        //see if the event failed the figure of Merit test, if so, see if it is an 
        //overlapping muon
        if( status->GetFigureOfMeritFailure() )
          status->SetVOverlappingMultiple(util.IsOverlappingMultiple(goldenPlaneItr, status->GetVertexZPosition(),
                                                                acd.GetDouble("HoughInterceptRMSLimit"), 
                                                                acd.GetDouble("HoughSlopeRMSLimit"),
                                                                acd.GetDouble("HoughPeakLimit")));
      }

      //done using the first view, so clear the selection function
      //MSG("Dmx", Msg::kInfo) <<"Done with View " << viewCtr << endl;
      vertexPlaneItr.GetSet()->Slice();
      leadPlaneItr.GetSet()->Slice();
      leadPlaneItr.GetSet()->AdoptSelectKeyFunc(0);
      leadPlaneItr.Reset();
      goldenPlaneItr.GetSet()->Slice();
      goldenPlaneItr.GetSet()->AdoptSelectKeyFunc(0);
      goldenPlaneItr.Reset();

      ++viewCtr;
    } //end loop over views
    
  }//end if event has a vertex
  return;
}

//______________________________________________________________________
void AlgDeMuxGolden::Trace(const char * /* c */) const
{
}

//-----------------------------------------------------------------------------------
//function for implementing the golden plane demux algorithm
void AlgDeMuxGolden::UseGoldenFit(DmxPlaneItr &leadPlaneItr, DmxPlaneItr &goldenPlaneItr,
                                  Int_t vertexPlane, Int_t endPlane)
{
  
  //get the number of loops you have to iterate over.  you have (# golden planes - 1)
  //loops to go through them each
  Int_t numLoops = goldenPlaneItr.SizeSelect() - 1;
  Int_t loopCtr = 0;
  DmxUtilities util = DmxUtilities();

  goldenPlaneItr.Reset();
  
  while(loopCtr < numLoops){
    
    
    Int_t goldCtr = 0;
    Float_t goldCoG1 = -1.;
    Float_t goldCoG2 = -1.;
    Float_t goldZ1 = -1.;
    Float_t goldZ2 = -1.;
    Int_t firstPlane = -1;
    Int_t lastPlane = -1;
    //get the cog's for the next two golden planes
    while(goldCtr < 2 && goldenPlaneItr.IsValid()){

      //advance the iterator only if on the first plane of the set - that way the first plane
      //in the next set is the last in this one
      if(goldCtr == 0){
        firstPlane = dynamic_cast<DmxPlane*>(goldenPlaneItr.Ptr())->GetPlaneNumber();
        goldCoG1 = dynamic_cast<DmxPlane*>(goldenPlaneItr.Ptr())->GetCoG();
        goldZ1 = dynamic_cast<DmxPlane*>(goldenPlaneItr.Ptr())->GetZPosition();
        goldenPlaneItr.Next();
      }
      else if(goldCtr == 1){
        lastPlane = dynamic_cast<DmxPlane*>(goldenPlaneItr.Ptr())->GetPlaneNumber();
        goldCoG2 = dynamic_cast<DmxPlane*>(goldenPlaneItr.Ptr())->GetCoG();
        goldZ2 = dynamic_cast<DmxPlane*>(goldenPlaneItr.Ptr())->GetZPosition();
      }
      ++goldCtr;
    }//end loop over golden planes

    //find the slope and intercept between the two golden plane centers of gravity
    Float_t slope = (goldCoG2 - goldCoG1)/(goldZ2 - goldZ1);
    Float_t intercept = goldCoG2 - goldZ2*slope;
    
    //MSG("DMX", Msg::kInfo) << "first gold plane " << goldCoG1 << " " << goldZ1 
    //             << " second " << goldCoG2 << " " << goldZ2 << " " << slope << " " << intercept << endl; 

    //slice the lead plane iterator to be between the golden planes.
    //if first loop, slice from vertex to second golden plane
    //if last loop, slice from first golden plane to end plane
    
    if(loopCtr == 0) firstPlane = vertexPlane;
    else if(loopCtr == numLoops-1) lastPlane = endPlane;

    leadPlaneItr.GetSet()->Slice(firstPlane, lastPlane);

    SetPlanesToDeterminedFit(leadPlaneItr, intercept, slope, firstPlane, lastPlane);

    fStrayPlanes += util.FindPlanesOffFit(leadPlaneItr, fStrayCut);

    //clear the lead plane slice
    leadPlaneItr.GetSet()->Slice();

    ++loopCtr;
  }//end loop over sets of goldenPlanes

  leadPlaneItr.Reset();
  goldenPlaneItr.Reset();
  return;

}

//--------------------------------------------------------------------------
//function to set the planes in the event
void AlgDeMuxGolden::SetPlanesToDeterminedFit(DmxPlaneItr &planeItr, Double_t a, 
                                              Double_t b,
                                              Int_t firstPlane, Int_t lastPlane)
{
  planeItr.ResetFirst();
 
  while( planeItr.IsValid() ){

    DmxPlane *plane = dynamic_cast<DmxPlane *>(planeItr.Ptr());
    //MSG("DMX", Msg::kInfo)<<"\tin SetPlanesToDeterminedFit " << plane->GetPlaneNumber() << endl;
      
    if(plane->GetPlaneNumber()>=firstPlane && plane->GetPlaneNumber()<=lastPlane){
      Double_t fitCoG = a + (plane->GetZPosition() * b);
      
      plane->SetStrips(fitCoG); 
    }//end if plane is in the window bounds
  
    planeItr.Next();
  }

  planeItr.Reset();
  
  return;
}

---