DmxMuonPlane.cxx

Go to the documentation of this file.

//$Id: DmxMuonPlane.cxx,v 1.52 2007/11/11 07:36:22 rhatcher Exp $
//
//DmxMuonPlane.cxx
//
//DmxMuonPlane 
//
//Author:  B. Rebel 7/2000

//#include "DeMux/DmxMuonPlane.h"
#include "DeMux/DmxMuonPlane.h"
#include "MessageService/MsgService.h"
#include "Navigation/NavSet.h"
#include "Conventions/PlaneView.h"
#include "Conventions/StripEnd.h"
#include "UgliGeometry/UgliStripHandle.h"
#include "UgliGeometry/UgliScintPlnHandle.h"
#include "Plex/PlexStripEndId.h"
#include "Algorithm/AlgConfig.h"
#include "CandDigit/CandDeMuxDigit.h"
#include "CandDigit/CandDeMuxDigitHandle.h"
#include "TMath.h"

ClassImp(DmxMuonPlane)

//______________________________________________________________________
CVSID("$Id: DmxMuonPlane.cxx,v 1.52 2007/11/11 07:36:22 rhatcher Exp $");

//----------------------------------------------------------------------
DmxMuonPlane::DmxMuonPlane()
{
}

//----------------------------------------------------------------------
DmxMuonPlane::DmxMuonPlane(AlgConfig & /* acd */, 
                           CandDeMuxDigitHandleItr &cdhitr, 
                           Int_t planeNumber) :
  fCdhit(cdhitr),
  fCoG(0.),
  fInitialCoG(-1.),
  fNumberOfStrips(1), //its a muon plane - can only have one strip hit
  fPlaneCharge(0.),
  fPlaneNumber(planeNumber),
  fPlaneType(DmxPlaneTypes::kMuon),
  fPlaneView(PlaneView::kUnknown),
  fStripsSet(0),
  fIsValid(false),
  fIsStray(false),
  fIsGolden(false),
  fUgh(0)
{
  fCdhit.Reset();
  
  //get the plane orientation
  fPlaneView = fCdhit.Ptr()->GetPlexSEIdAltL().GetPlaneView();

  //get the ugli geometry handle from a cand digit
  fUgh = new UgliGeomHandle(*fCdhit.Ptr()->GetVldContext());

  Int_t numDigits = 0;
  while( fCdhit.IsValid() ){
      if( fCdhit.Ptr()->GetDeMuxDigitFlagWord() == 0) ++numDigits;
      fCdhit.Next();
  }
  fCdhit.Reset();
  
  //MSG("Dmx", Msg::kDebug) << "Number of digits for plane " 
  //             << fPlaneNumber << " = " << numDigits << endl;
  
  if( numDigits == 2){
    //check to see that you have one from each side
    Int_t west = 0;
    Int_t east = 0;
    CandDeMuxDigitHandle *westDigit = fCdhit.Ptr();
    CandDeMuxDigitHandle *eastDigit = fCdhit.Ptr();
    
    while( fCdhit.IsValid() ){
      fCdhit.Ptr()->GetPlexSEIdAltLWritable().ClearWeights();
      fCdhit.Ptr()->GetPlexSEIdAltL().SetFirst();
      
      if( fCdhit.Ptr()->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast 
          && fCdhit.Ptr()->GetDeMuxDigitFlagWord() == 0){
        east = 1; 
        eastDigit = fCdhit.Ptr();
      }
      else if( fCdhit.Ptr()->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest 
               && fCdhit.Ptr()->GetDeMuxDigitFlagWord() == 0){
        west = 1; 
        westDigit = fCdhit.Ptr();
      }
      
      fCdhit.Next();
    }
    if(west + east == 2){
      //loop over the strip possibilities to see if they match up on each side
      Int_t westStrip = -1;
      Int_t eastStrip = -2;
      
      //reset the SEId alt lists
      westDigit->GetPlexSEIdAltL().SetFirst();
      
      while( westDigit->GetPlexSEIdAltL().IsValid() ){
        westStrip = westDigit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();
        //MSG("Dmx", Msg::kDebug) << fPlaneNumber << " west strip = " 
        //             << westStrip << endl;
        
        eastDigit->GetPlexSEIdAltL().SetFirst();
      
        while( eastDigit->GetPlexSEIdAltL().IsValid() ){
          eastStrip = eastDigit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();
          //MSG("Dmx", Msg::kDebug) << "\teast strip = " << eastStrip << endl;
          
          if(eastStrip == westStrip){
            westDigit->GetPlexSEIdAltLWritable().SetCurrentWeight(1.0);
            eastDigit->GetPlexSEIdAltLWritable().SetCurrentWeight(1.0);

            //set the veto flags to 0
            westDigit->GetPlexSEIdAltLWritable().SetDemuxVetoFlag(0);
            eastDigit->GetPlexSEIdAltLWritable().SetDemuxVetoFlag(0);
            
            //print the association
            // MSG("Dmx", Msg::kDebug) << "west strip = " 
//                                  << westDigit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip() 
//                                  << "\tsignal = " << westDigit->GetCharge() << endl;
//          MSG("Dmx", Msg::kDebug) << "east strip = " 
//                                 << eastDigit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip() 
//                                 << "\tsignal = " << eastDigit->GetCharge() << endl;
            
            //set the center of gravity
            fCoG = westDigit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();
            fInitialCoG = fCoG;
            
            //set the true muon plane flag but be careful - there could be 
            //multiple solutions for this plane so if the is valid has been
            //set to true before, then you have multiple solutions and should 
            //set it to false
            if(!fIsValid){
              fIsValid = true;
              fIsGolden = true;
            }
            else{
              fIsValid = false;
              fIsGolden = false;
            }
            //set the strips set flag
            fStripsSet = 1;
          }
          eastDigit->GetPlexSEIdAltL().Next();
        }

        westDigit->GetPlexSEIdAltL().Next();
      }
    }
  }
 
  MSG("DmxMuonPlane", Msg::kDebug) << "plane = " << fPlaneNumber 
                                                                  << " initial muon cog = " << fInitialCoG << "/" << fCoG << endl;
  
  //if this is a golden plane, set any digits that might not be set because
  //of cross talk or low signal
  if(fIsGolden) SetStrips();

  //get the charge for the plane
  fCdhit.ResetFirst();
  while(fCdhit.IsValid()){
    //if(fCdhit.GetSet()->GetMasks().GetMask(fCdhit.Ptr())){fPlaneCharge += fCdhit.Ptr()->GetCharge();}
    if( fCdhit.Ptr()->GetDeMuxDigitFlagWord() == 0) 
      fPlaneCharge += fCdhit.Ptr()->GetCharge();
    fCdhit.Next();
  }
  fCdhit.Reset();
  return;
}

//----------------------------------------------------------------------
DmxMuonPlane::~DmxMuonPlane()
{
  delete fUgh;
  //  MSG("Dmx", Msg::kDebug) << "deleting Muon Plane for plane " 
  //              << fPlaneNumber << endl;  
}

//----------------------------------------------------------------------
Int_t DmxMuonPlane::GetNumberOfStrips() const
{
        return fNumberOfStrips;
}
//----------------------------------------------------------------------
Float_t DmxMuonPlane::GetCoG() const
{
  //get the strip end id for the center of gravity for this hypothesis
  PlexStripEndId seid(Detector::kFar,fPlaneNumber,TMath::Nint(fCoG));
  UgliStripHandle ush = fUgh->GetStripHandle(seid);

  return ush.GetTPos();
  //return fCoG;
}

//----------------------------------------------------------------------
Float_t DmxMuonPlane::GetInitialCoG() const
{
        
        MSG("DmxMuonPlane", Msg::kDebug) << "plane = " << fPlaneNumber 
                                                                        << " initial muon cog = " << fInitialCoG << "/" << fCoG << endl;
                
        //get the strip end id for the center of gravity for this hypothesis
        PlexStripEndId seid(Detector::kFar,fPlaneNumber,TMath::Nint(fInitialCoG));
        UgliStripHandle ush = fUgh->GetStripHandle(seid);

        return ush.GetTPos();
        //return fCoG;
}

//----------------------------------------------------------------------
Int_t DmxMuonPlane::GetHypothesisLowerBound(Float_t tPos) const
{
  //change the transverse position for the fit into a strip number
  PlexStripEndId seid(Detector::kFar, fPlaneNumber, 1);
  UgliScintPlnHandle usph = fUgh->GetScintPlnHandle(seid);
  UgliStripHandle ush = usph.GetClosestStrip(tPos);
  Float_t fitCoG = 1.*ush.GetSEId().GetStrip();
  
  Int_t hypLowerBound = 0;
  if(fitCoG >= 11.5 && fitCoG <= 179.5){hypLowerBound = TMath::Nint(fitCoG-11.5);}
  else if(fitCoG < 11.5){hypLowerBound = 0;}
  else if(fitCoG > 179.5){hypLowerBound = 168;}
  
  return hypLowerBound;
}

//----------------------------------------------------------------------
Float_t DmxMuonPlane::GetStripCoG() const
{
  return fCoG;
}

//----------------------------------------------------------------------
Float_t DmxMuonPlane::GetInitialStripCoG() const
{
  return fInitialCoG;
}

//----------------------------------------------------------------------
Float_t DmxMuonPlane::GetPlaneCharge() const
{
  return fPlaneCharge;
}

//----------------------------------------------------------------------
Int_t DmxMuonPlane::GetPlaneNumber() const  
{
  return fPlaneNumber;
}

//----------------------------------------------------------------------
Float_t DmxMuonPlane::GetZPosition() const  
{
  PlexStripEndId seid(Detector::kFar, fPlaneNumber, 1);
  UgliScintPlnHandle usph = fUgh->GetScintPlnHandle(seid);
  return usph.GetZ0();
  //return 1.*fPlaneNumber;
}

//----------------------------------------------------------------------
DmxPlaneTypes::PlaneType_t DmxMuonPlane::GetPlaneType() const
{
  return fPlaneType;
}

//----------------------------------------------------------------------
PlaneView::PlaneView_t DmxMuonPlane::GetPlaneView() const 
{
  return fPlaneView;
}

//----------------------------------------------------------------------
Bool_t DmxMuonPlane::GetStripsSetFlag() const
{
  return fStripsSet>0;
}

//----------------------------------------------------------------------
Bool_t DmxMuonPlane::IsValid() const
{
  return fIsValid;
}

//----------------------------------------------------------------------
Bool_t DmxMuonPlane::IsStray() const
{
  return fIsStray;
}

//----------------------------------------------------------------------
Bool_t DmxMuonPlane::IsGolden() const
{
  return fIsGolden;
}

//-----------------------------------------------------------------------
void DmxMuonPlane::PrintRecon()
{
  fCdhit.Reset();
  MSG("Dmx", Msg::kInfo) << "Plane = " << fPlaneNumber 
                         << "\tIsValid = " << (Int_t)fIsValid << endl;
  while( fCdhit.IsValid() ){
    CandDeMuxDigitHandle *currentDigit = fCdhit.Ptr();
    
    if(currentDigit->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest){
      MSG("Dmx", Msg::kInfo) << "\tWest\t" ;
    }
    else if(currentDigit->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast){
      MSG("Dmx", Msg::kInfo) << "\tEast\t" ;
    }
    MSG("Dmx", Msg::kInfo) << currentDigit->GetPlexSEIdAltL().GetBestSEId().GetStrip() << "\t";
    MSG("Dmx", Msg::kInfo) << currentDigit->GetCharge() << endl;
    
    fCdhit.Next();
  }
  return;
}
//----------------------------------------------------------------------
void DmxMuonPlane::SetStrips()
{
  Float_t fitCoG = fInitialCoG;
  Int_t hypLower = 0;
  if(fitCoG >= 11.5 && fitCoG <= 179.5){hypLower = TMath::Nint(fitCoG-11.5);}
  else if(fitCoG < 11.5){hypLower = 0;}
  else if(fitCoG > 179.5){hypLower = 168;}
  
  fStripsSet = 1;
  fCdhit.Reset();
  
  Double_t cogNum = 0.;
  Double_t cogDenom = 0.;
  
  while( fCdhit.IsValid() ){
      
    //get the next CandDigitHandle object
    CandDeMuxDigitHandle *currentDigit = fCdhit.Ptr();
    
    //clear all weights for the digit
    currentDigit->GetPlexSEIdAltLWritable().ClearWeights();
    
    //set the alternatives list to the first one
    currentDigit->GetPlexSEIdAltL().SetFirst();

    //set veto flag to nonzero, if there is a SEId in this hypothesis set it to
    //0
    currentDigit->GetPlexSEIdAltLWritable().SetDemuxVetoFlag(1);
    
    //loop through the Strip End Alternative list and see which one is
    //in the current hypothesis
    while( currentDigit->GetPlexSEIdAltL().IsValid() ){
      
      if( currentDigit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip() 
          >= hypLower && currentDigit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip() 
          <= (hypLower+23)){
        
        //set the Strip weight to 1 if it is the first choice above the 
        //lowerbound and set veto flag to 0
        currentDigit->GetPlexSEIdAltLWritable().SetCurrentWeight(1.0);
        currentDigit->GetPlexSEIdAltLWritable().SetDemuxVetoFlag(0);

        //print the association
        //MSG("Dmx", Msg::kVerbose) << "strip = " 
        //                << currentDigit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip() 
        //                << "\tside = " 
        //                << (Int_t)currentDigit->GetPlexSEIdAltL().GetCurrentSEId().GetEnd() 
        //                << "\tsignal = " << currentDigit->GetCharge() << endl;
        
        //set the cog to this strip, if it isnt a cross talk digit
        if(currentDigit->GetDeMuxDigitFlagWord() == 0){
          //if(fCdhit.GetSet()->GetMasks().GetMask(fCdhit.Ptr())){
          cogNum += currentDigit->GetCharge()*currentDigit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();
          cogDenom += currentDigit->GetCharge();
        }
                
        //you found the Strip End Id so quit looking
        break;
      }
      //if that Strip End Id wasnt in the hypothesis go to the next one
      currentDigit->GetPlexSEIdAltL().Next();
    }

    fCdhit.Next();
  }
   
  if(cogDenom >0.) fCoG = cogNum/cogDenom;

  return;
}

//----------------------------------------------------------------------
void DmxMuonPlane::SetStrips(Float_t tPos)
{
        //change the transverse position for the fit into a strip number
        PlexStripEndId seid(Detector::kFar, fPlaneNumber, 1);
        UgliScintPlnHandle usph = fUgh->GetScintPlnHandle(seid);
        UgliStripHandle ush = usph.GetClosestStrip(tPos);
        Float_t fitCoG = 1.*ush.GetSEId().GetStrip();
        //Float_t fitCoG = tPos;
  
        Int_t hypLower = 0;
        if(fitCoG >= 11.5 && fitCoG <= 179.5) hypLower = TMath::Nint(fitCoG-11.5);
        else if(fitCoG < 11.5) hypLower = 0;
        else if(fitCoG > 179.5) hypLower = 168;
  
        MSG("DmxMuonPlane", Msg::kDebug) << "\tfit in strip vs plane space " << fPlaneNumber 
                                                                         << "\t" << fitCoG << " " << hypLower << endl; 
        MSG("DmxMuonPlane", Msg::kDebug) << "\tfit in tpos vs z space " << GetZPosition() 
                                                                         << "\t" << tPos << endl; 
  
        fStripsSet = 1;
        fCdhit.Reset();
  
        Double_t cogNum = 0.;
        Double_t cogDenom = 0.;
  
        while( fCdhit.IsValid() ){
      
                //get the next CandDigitHandle object
                CandDeMuxDigitHandle *currentDigit = fCdhit.Ptr();
    
                //clear all weights for the digit
                currentDigit->GetPlexSEIdAltLWritable().ClearWeights();
    
                //set the alternatives list to the first one
                currentDigit->GetPlexSEIdAltLWritable().SetFirst();

                //set veto flag to nonzero, if there is a SEId in this hypothesis set it to
                //0
                currentDigit->GetPlexSEIdAltLWritable().SetDemuxVetoFlag(1);
    
                //loop through the Strip End Alternative list and see which one is
                //in the current hypothesis
                while( currentDigit->GetPlexSEIdAltLWritable().IsValid() ){
                
                        MSG("DmxMuonPlane", Msg::kDebug) << "plane = " 
                                                                                         << currentDigit->GetPlexSEIdAltLWritable().GetPlane()
                                                                                         << " strip alt = " 
                                                                                         << currentDigit->GetPlexSEIdAltLWritable().GetCurrentSEId().GetStrip() 
                                                                                         << endl;
      
                        if( currentDigit->GetPlexSEIdAltLWritable().GetCurrentSEId().GetStrip() 
                                >= hypLower && currentDigit->GetPlexSEIdAltLWritable().GetCurrentSEId().GetStrip() 
                                <= (hypLower+23)){
        
                                //set the Strip weight to 1 if it is the first choice above the 
                                //lowerbound and set veto flag to 0
                                currentDigit->GetPlexSEIdAltLWritable().SetCurrentWeight(1.0);
                                currentDigit->GetPlexSEIdAltLWritable().SetDemuxVetoFlag(0);
                                
                                MSG("DmxMuonPlane", Msg::kDebug) << "strip = " 
                                                                                                 << currentDigit->GetPlexSEIdAltLWritable().GetCurrentSEId().GetStrip() 
                                                                                                 << " " << currentDigit->GetPlexSEIdAltLWritable().GetCurrentWeight() 
                                                                                                 << endl;
                                
                                //print the association
                                //MSG("Dmx", Msg::kVerbose) << "strip = " 
                                //                << currentDigit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip() 
                                //                << "\tside = " 
                                //                << (Int_t)currentDigit->GetPlexSEIdAltL().GetCurrentSEId().GetEnd() 
                                //                << "\tsignal = " << currentDigit->GetCharge() << endl;
        
                                //set the cog to this strip, if it isnt a cross talk digit
                                if(currentDigit->GetDeMuxDigitFlagWord() == 0){
                                        cogNum += currentDigit->GetCharge()*currentDigit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();
                                        cogDenom += currentDigit->GetCharge();
                                }
                
                                //you found the Strip End Id so quit looking
                                break;
                        }
                        //if that Strip End Id wasnt in the hypothesis go to the next one
                        currentDigit->GetPlexSEIdAltL().Next();
                }

                fCdhit.Next();
        }
   
        if(cogDenom >0.) fCoG = cogNum/cogDenom;

        return;
}

//----------------------------------------------------------------------
void DmxMuonPlane::SetStray(bool stray)
{
  fIsStray = stray;
  return;
}

//----------------------------------------------------------------------
void DmxMuonPlane::SetGolden(bool golden)
{
  fIsGolden = golden;
  return;
}

//------------------------------------------------------------------------
Float_t DmxMuonPlane::GetTimingOffset()
{
//make arrays to hold the clear fiber + wls pigtail lengths for each digit
  Float_t westFiber[192];
  Float_t eastFiber[192];
  Float_t eastTime[192];
  Float_t westTime[192];

  //reset the arrays
  for(Int_t i = 0; i<192; i++){
    westFiber[i] = -1.;
    eastFiber[i] = -1.;
    eastTime[i] = -1.;
    westTime[i] = -1.;
  }
    
  fCdhit.Reset();

  while(fCdhit.IsValid()){
    
    CandDeMuxDigitHandle *digit = fCdhit.Ptr();
    
    if(digit->GetPlexSEIdAltL().GetDemuxVetoFlag() == 0 && digit->GetDeMuxDigitFlagWord() == 0){
      
      Int_t strip = digit->GetPlexSEIdAltL().GetBestSEId().GetStrip();
      
      //get the side for this digit and fill the array
      UgliStripHandle ush = fUgh->GetStripHandle(digit->GetPlexSEIdAltL().GetBestSEId());
      
      if( digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast ){
        eastFiber[strip] = ush.ClearFiber(StripEnd::kEast) + ush.WlsPigtail(StripEnd::kEast);
        eastTime[strip] = 1e9*digit->GetTime(CalTimeType::kT0) - 2093./(102.+TMath::Power(1.*digit->GetCharge(),1.2));
      }
      else if( digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest ){
        westFiber[strip] = ush.ClearFiber(StripEnd::kWest) + ush.WlsPigtail(StripEnd::kWest);
        westTime[strip] = 1e9*digit->GetTime(CalTimeType::kT0) - 2093./(102.+TMath::Power(1.*digit->GetCharge(),1.2));
      }
      
    }
    fCdhit.Next();
  }
  
  //find the average offset from center for this plane based on timing and fill the array
  Float_t avTimeNum = 0.;
  Float_t avTimeDenom = 0.;
  Float_t offset = -10.;
  for(Int_t k=0; k<192; k++){
    
    //if the strip has doubled sided readout, then use it
    if(westTime[k] != -1.&& eastTime[k] != -1.){
      offset = (0.0825*(eastTime[k] - westTime[k]) + (westFiber[k]-eastFiber[k])/2.);
      //MSG("DmxTim", Msg::kInfo) << fPlaneNumber << " offset = " << offset << endl;
      avTimeNum += offset;
      avTimeDenom += 1.;
    }
  }
        
  if(avTimeDenom>0) offset = avTimeNum/avTimeDenom;

  //MSG("DmxTim", Msg::kDebug) << fPlaneNumber << " av offset = " << offset << endl;
  
  //clear the digit slice
  fCdhit.Reset();

  return offset;
} 

---