DmxHypothesis.cxx

Go to the documentation of this file.

//$Id: DmxHypothesis.cxx,v 1.62 2007/02/04 05:55:06 rhatcher Exp $
//
//DmxHypothesis.cxx
//
//DmxHypothesis owns the description of signal to strip 
//correllations for the hypothesis.  it is passed a pointer to the
//stat to use to evaluate the goodness of the reconstruction
//
//Author:  B. Rebel 6/2000

//#include "DeMux/DmxHypothesis.h"
#include "DeMux/DmxHypothesis.h"
#include "MessageService/MsgService.h"
#include "MessageService/MsgFormat.h"
#include "DeMux/DmxChiSqrStat.h"
#include "DeMux/DmxRMSStat.h"
#include "DeMux/DmxStatTypes.h"
#include "Navigation/NavSet.h"
#include "UgliGeometry/UgliStripHandle.h"
#include "UgliGeometry/UgliScintPlnHandle.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "Algorithm/AlgConfig.h"
#include "CandDigit/CandDeMuxDigit.h"
#include "CandDigit/CandDeMuxDigitHandle.h"

#include "TMath.h"

ClassImp(DmxHypothesis)
XXXITRIMP(DmxHypothesis)

//______________________________________________________________________
CVSID("$Id: DmxHypothesis.cxx,v 1.62 2007/02/04 05:55:06 rhatcher Exp $");

//----------------------------------------------------------------------
DmxHypothesis::DmxHypothesis() :
  fStat(0)
{
}

//----------------------------------------------------------------------
DmxHypothesis::DmxHypothesis(AlgConfig &acd, CandDeMuxDigitHandleItr &cdhitr,
                             Int_t lowerStrip, Int_t upperStrip) :
  fAllUsed(true),
  fCdhit(cdhitr),
  fCompare(0),
  fCoG(0.),
  fLowerBound(lowerStrip),
  fMaxSeparation(23),
  fMatedSignalRatio(0.0),
  fRequiredMatedSignalRatio(acd.GetDouble("RatioMatedSignalForValid")),
  fStat(0),
  fStripsUsed(0),
  //fTimeOffset(-10.),
  fUpperBound(upperStrip)
{
  
  fCdhit.Reset();

  //loop over the digits to set the values

  //variables to calculate the center of gravity for the hypothesis
  Float_t numerator = 0.;       
  Float_t denominator = 0.;

  //initialize the signal and strip arrays
  Float_t signalE[192];
  Float_t signalW[192];

  //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++){
   }

  for(Int_t vv = 0; vv < 192; vv++){
    signalE[vv] = 0.;
    signalW[vv] = 0.;
    //westFiber[vv] = -1.;
    //eastFiber[vv] = -1.;
    //eastTime[vv] = -1.;
    //westTime[vv] = -1.;
 }

  while( fCdhit.IsValid() ){
    
    //get the next CandDeMuxDigitHanlde object
    CandDeMuxDigitHandle *currentDigit = fCdhit.Ptr();
    
    //set the alternatives list to the first one
    currentDigit->GetPlexSEIdAltL().SetFirst();
      
    //loop through the Strip End Alternative list and see which one is
    //in the current hypothesis
    
    bool digitUsed = false;

    while( currentDigit->GetPlexSEIdAltL().IsValid() ){
      Int_t strip = currentDigit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip();
      if(strip >= fLowerBound && strip <= fUpperBound){

        //if the flag isnt set, then use the digit - flags are set only if
        //the signal fraction is below the limit or the signal is too low
        if(currentDigit->GetDeMuxDigitFlagWord() == 0){
          
          digitUsed = true;
          
          if( currentDigit->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast ){
            signalE[strip] = currentDigit->GetCharge();
            //MSG("DmxHyp", Msg::kDebug)<<"\teast strip " << strip << " charge " << signalE[strip] << endl;
          }
          else if( currentDigit->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest ){
            signalW[strip] = currentDigit->GetCharge();
            //MSG("DmxHyp", Msg::kDebug)<<"\twest strip " << strip << " charge " << signalW[strip] << endl;
          }
          
          numerator += currentDigit->GetCharge() * strip;
          denominator += currentDigit->GetCharge();
                
        }//end if this is a good digit to use

        //you found the strip for this digit, so break the loop
        break;
      }//end if stripend is in hypothesis limit

      //if that Strip End Id wasnt in the hypothesis go to the next one
      currentDigit->GetPlexSEIdAltL().Next();

    }//end loop over strip ends
    
    if( !digitUsed && currentDigit->GetDeMuxDigitFlagWord() == 0){ fAllUsed = false; }
    //advance to the next digit
    fCdhit.Next(); 
  }//end loop over digits

  fCdhit.ResetFirst();

  //find the mated signal and get the timing offset for this hypothesis
  //also see how many strips are used in this hypothesis
  Float_t matches = 0.0;

  for(Int_t j = 0; j < 192; j++){
          if( signalW[j] > 0. && signalE[j] > 0.)
                  matches += signalW[j] + signalE[j]; 
          if(signalW[j]>0. || signalE[j]>0.) ++fStripsUsed;
  }
  
  if(denominator > 0.){
    fCoG = numerator / denominator;
    fMatedSignalRatio = matches / denominator;
  }

  //get the correct DmxStatistic object from the AlgDeMuxConfig object
      
  //if(acd.GetStatParameter() == DmxStatTypes::kChi2){ 
    fStat = new DmxChiSqrStat(signalW, signalE, fCoG); 
    //}
    //else if(acd.GetStatParameter() == DmxStatTypes::kRms){ 
    //fStat = new DmxRMSStat(signalW, signalE, fCoG); 
    //}
  
  return;
}

//-----------------------------------------------------------------------
DmxHypothesis::~DmxHypothesis()
{
  //MSG("Dmx", Msg::kVerbose) << "deleting DmxHypothesis with lowerbound = " 
  //              << fLowerBound << endl;
  delete fStat;
}

//-----------------------------------------------------------------------
Int_t DmxHypothesis::GetNumberOfStripsUsed() const
{
        return fStripsUsed;
}

//-----------------------------------------------------------------------
Float_t DmxHypothesis::GetCoG() const
{
  return fCoG;
}

//-----------------------------------------------------------------------
Int_t DmxHypothesis::GetCompareFlag() const 
{
  return fCompare;
}

//-----------------------------------------------------------------------
Int_t DmxHypothesis::GetLowerBound() const 
{
  return fLowerBound;
}

//-----------------------------------------------------------------------
Float_t DmxHypothesis::GetMatedSignalRatio() const  
{
  return fMatedSignalRatio;
}

//-----------------------------------------------------------------------
Float_t DmxHypothesis::GetStat() const  
{
  return fStat->GetGoodness();
}

//-----------------------------------------------------------------------
DmxStatistic *DmxHypothesis::GetStatObject() const  
{
  return fStat;
}

//-----------------------------------------------------------------------
Float_t DmxHypothesis::GetTieBreakerStat() const 
{
  return fStat->GetTieBreaker();
}

//-----------------------------------------------------------------------
Int_t DmxHypothesis::GetUpperBound() const 
{
  return fUpperBound;
}

//-----------------------------------------------------------------------
Bool_t DmxHypothesis::IsValid() const
{
  //MSG("Dmx", Msg::kDebug) << "used flag = " << fAllUsed 
  //              << " for lowerbound = " << fLowerBound << endl;
  //return fMaxSeparation<24 && fMaxSeparation>0 && fAllUsed>0 && fMatedSignalRatio>0.3;
  //return fAllUsed>0 && fMatedSignalRatio>fRequiredMatedSignalRatio;
  return fMatedSignalRatio>fRequiredMatedSignalRatio && fAllUsed;
}

//-----------------------------------------------------------------------
void DmxHypothesis::PrintRecon()
{
  MSG("DMX", Msg::kInfo) << "\tLower Bound = " << fLowerBound
                         << "\tGoodness of Fit = " << fStat->GetGoodness()
                         << "\tTieBreaker = " << fStat->GetTieBreaker() << endl;
  while( fCdhit.IsValid() ){
    while( fCdhit.Ptr()->GetPlexSEIdAltL().IsValid() ){
        if( fCdhit.Ptr()->GetPlexSEIdAltL().GetCurrentSEId().GetStrip() >= fLowerBound && 
            fCdhit.Ptr()->GetPlexSEIdAltL().GetCurrentSEId().GetStrip() <= fUpperBound){

          if( fCdhit.Ptr()->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest ){
            MSG("DMX", Msg::kInfo)<< "\t\tWest\t" 
                                  << fCdhit.Ptr()->GetPlexSEIdAltL().GetCurrentSEId().GetStrip() 
                                  << "\t" << fCdhit.Ptr()->GetCharge() << endl;
          }
          if( fCdhit.Ptr()->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast ){
            MSG("DMX", Msg::kInfo)<< "\t\tEast\t" 
                                  << fCdhit.Ptr()->GetPlexSEIdAltL().GetCurrentSEId().GetStrip() 
                                  << "\t" << fCdhit.Ptr()->GetCharge() << endl;
          }
        }
    }
    fCdhit.Next();
  }

  return;
}

//----------------------------------------------------------------------
void DmxHypothesis::SetStrips()
{
  fCdhit.Reset();
  while( fCdhit.IsValid() ){
    
    //get the next CandDigitHanlde 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 the veto flag to nonzero, if the digit has a SEId in this hypothesis
    //set it back 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() >= fLowerBound && 
            currentDigit->GetPlexSEIdAltL().GetCurrentSEId().GetStrip() <= fUpperBound){

          //set the Strip weight to 1 if it is in the hypothesis and veto flag to 0
          currentDigit->GetPlexSEIdAltLWritable().SetCurrentWeight(1.0);
          currentDigit->GetPlexSEIdAltLWritable().SetDemuxVetoFlag(0);

          //print the association
          //MSG("Dmx", Msg::kDebug) << "strip = " 
          //  << currentDigit->GetPlexSEIdAltLWritable().GetCurrentSEId().GetStrip() 
          //      << "\tside = " 
          //      << (Int_t)currentDigit->GetPlexSEIdAltLWritable().GetCurrentSEId().GetEnd() 
          //      << "\tsignal = " << currentDigit->GetCharge() 
          //  << endl;

          //you found the Strip End Id in the current hypothesis so quit looking
          break;
        }
        //if that Strip End Id wasnt in the hypothesis go to the next one
        currentDigit->GetPlexSEIdAltL().Next();
    }
    
    fCdhit.Next();
  }
  fCdhit.Reset();
}

//-------------------------------------------------------------------------------
// Float_t DmxHypothesis::GetTimingOffset()
// {
//   MSG("DmxHyp", Msg::kDebug) << "timing offset = " << fTimeOffset << endl;
//   return fTimeOffset;
// }

//--------------------------------------------------------------------------------
Float_t DmxHypothesis::GetTimingOffset(UgliGeomHandle *ugh)
{
  
  //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.;
  }

  //MSG("DmxHyp", Msg::kDebug) << "starting the timing offset in hypothesis" << endl;
    
  fCdhit.ResetFirst();
  //MSG("DmxHyp", Msg::kDebug) << "reset digit iterator in hypothesis" << endl;

  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 = ugh->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::kDebug) << "shower plane offset = " << offset << endl;
      avTimeNum += offset;
      avTimeDenom += 1.;
    }
  }
        
  if(avTimeDenom>0) offset = avTimeNum/avTimeDenom;
  //MSG("DmxTim", Msg::kDebug) << "shower plane offset = " << offset << endl;
       
  //clear the digit slice
  fCdhit.Reset();

  return offset;
} 

---