ChopEvaluation.cxx

Go to the documentation of this file.

#include "ChopEvaluation.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "Calibrator/Calibrator.h"

//ClassImp(ChopEvaluation);
//ClassImp(NeutrinoInfo)
//ClassImp(ChopNeutrinoLeaf)
//ClassImp(ChopTreeLeaf)

ChopEvaluation::ChopEvaluation() :
  fNuMatch(50,0),
  fNuMatch_cal(50,0),
  fNuMatch_fid(50,0),
  fNuMatch_calt(50,0),
  fNuMatch_caltfid(50,0)
{
  Reset();
}


void ChopEvaluation::Reset()
{
  fTdcBeg = -1;
  fTdcEnd = -1;; 
  fSigcor = 0;
  fSigcor_cal = 0;
  fSigcor_fid = 0; 
  fMeanTime = 0; 
  fBigTdc = 0; 
  fBigTdcE = 0; 
  fBigPlane = 0; 
  fBigPlaneE = 0; 
  fPlanes = 0; 
  fStrips = 0;
  fCalStrips = 0;
  fBestNeutrino1 = -1;
  fBestNeutrino2 = -1;

  for(int i=0;i<kNProf;i++) {
    fProfile[i].i=i;
    fProfile[i].p=0;
  }
  /*
  fTdcProfile.clear();
  fPlaneProfile.clear();
  fNuMatch.clear();
  fNuMatch_cal.clear();
  fNuMatch_fid.clear();
  fNuMatch_calt.clear();
  fNuMatch_caltfid.clear();
  fStrips.clear();
  */
}

void ChopEvaluation::Evaluate(
                          const std::vector<CandDigitHandle>& inList, // digits to evaluate
                          const Truthifier& truth,                    // truth.
                          const std::vector<int>&    inNeutrinos,     // indecies of neutrinos in truth info
                          const std::vector<double>& inNeutrinoTZero, // time of first interaction of each nu
                          double inTCut
                          )  
{ 
  //
  //
  //
  if(inList.size()==0) return; // Null list.
  
  UgliGeomHandle ugli(*(inList[0].GetVldContext()));

  fTdcBeg = truth.GetRawDigit(inList[0])->GetTDC();
  fTdcEnd = truth.GetRawDigit(inList[0])->GetTDC();
  fPlaneBeg = inList[0].GetPlexSEIdAltL().GetPlane();
  fPlaneEnd = inList[0].GetPlexSEIdAltL().GetPlane();
  fDigits = inList.size();
  std::map<PlexStripEndId,int> hitStrips;    
  std::map<PlexStripEndId,int>::iterator hitStripsItr;    
    

  // First pass. Get start and end times.

  for(UInt_t idig=0;idig<inList.size();idig++) {
    CandDigitHandle cdh = inList[idig];
    int tdc = truth.GetRawDigit(inList[idig])->GetTDC();
    int plane = cdh.GetPlexSEIdAltL().GetPlane();
    PlexStripEndId seid = cdh.GetPlexSEIdAltL().GetBestSEId();

    if(tdc<fTdcBeg) fTdcBeg = tdc;
    if(tdc>fTdcEnd) fTdcEnd = tdc;
    if(plane<fPlaneBeg) fPlaneBeg = plane;
    if(plane>fPlaneEnd) fPlaneEnd = plane;
    hitStrips[seid]+=1;
  }
  fStrips = hitStrips.size();
  for(hitStripsItr=hitStrips.begin();  hitStripsItr!=hitStrips.end();  hitStripsItr++)
    if(hitStripsItr->first.GetPlane()<121) fCalStrips +=1;

  // resize arrays.
  fTdcProfile.resize(fTdcEnd-fTdcBeg+1,0);
  fPlaneProfile.resize(300,0);

  fNuMatch.        resize(inNeutrinos.size(),(double)0.);
  fNuMatch_cal.    resize(inNeutrinos.size(),(double)0.);
  fNuMatch_fid.    resize(inNeutrinos.size(),(double)0.);
  fNuMatch_calt.   resize(inNeutrinos.size(),(double)0.);
  fNuMatch_caltfid.resize(inNeutrinos.size(),(double)0.);
  fNuMatch_strips.   resize(inNeutrinos.size(),(double)0.);
  fNuMatch_calstrips.resize(inNeutrinos.size(),(double)0.);
  
  std::vector<std::map<PlexStripEndId,int> >  nuStrips(inNeutrinos.size());     
  std::map<PlexStripEndId,int>::iterator stripItr;


  for(UInt_t idig=0;idig<inList.size();idig++) {
    CandDigitHandle cdh = inList[idig];
    int tdc = truth.GetRawDigit(cdh)->GetTDC();
    float sig = cdh.GetCharge(CalDigitType::kSigCorr);
    int plane = cdh.GetPlexSEIdAltL().GetPlane();
    PlexStripEndId seid = cdh.GetPlexSEIdAltL().GetBestSEId();
    PlaneView::PlaneView_t view = cdh.GetPlexSEIdAltL().GetBestSEId().GetPlaneView();
    float tpos = ugli.GetStripHandle(cdh.GetPlexSEIdAltL().GetBestSEId()).GetTPos();
    
    bool cal = (plane<121);
    bool fid = false;
    if(plane>20&&plane<90) {
      if(  ( (view==PlaneView::kU) && (tpos > -0.3) && (tpos <  2.4) )  ||
           ( (view==PlaneView::kV) && (tpos > -2.4) && (tpos <  0.3) ) ) {
        fid = true;
      }
    }
    
    
    fSigcor += sig;
    if(cal) fSigcor_cal += sig;
    if(fid) fSigcor_fid += sig;

    fMeanTime = fMeanTime + sig*(double)tdc;
    
    fTdcProfile[tdc-fTdcBeg] += sig;
    fPlaneProfile[plane] += sig;

    for(UInt_t inu=0;inu<inNeutrinos.size();inu++) {
      int neu = inNeutrinos[inu];
      float frac  = truth.IsDigitFromNeutrino(neu,cdh);//,neuTrueStartTime[inu]+fParticleTimeout);        
      float fract = truth.IsDigitFromNeutrino(neu,cdh,
                                              inNeutrinoTZero[inu]+inTCut);//,neuTrueStartTime[inu]+fParticleTimeout);    
      
      fNuMatch[inu]+=frac*sig;
      if(cal) fNuMatch_cal[inu]   += frac*sig;
      if(fid) fNuMatch_fid[inu]   += frac*sig;
      if(cal) fNuMatch_calt[inu]    += fract*sig;
      if(fid) fNuMatch_caltfid[inu] += fract*sig;

      if(fract>0) nuStrips[inu][seid]++;
    }

    
  }

  //
  // Now finish evaluation..
  //

  fMeanTime = fMeanTime/fSigcor;

  // Loop on time buckets:
  for(UInt_t itdc=0; itdc<fTdcProfile.size(); itdc++) {
    int tdc = itdc+fTdcBeg;

    if(fTdcProfile[itdc]>fBigTdcE) {
      fBigTdcE = fTdcProfile[itdc];
      fBigTdc  = tdc;
    }
  }

  // Create profile.
  for(UInt_t itdc=0; itdc<fTdcProfile.size(); itdc++) {
    int tdc = itdc+fTdcBeg;
    int i = tdc - fBigTdc;
    int elem = i+5;
    if( (elem>=0) &&(elem<kNProf) ) {      
      fProfile[elem].p = fTdcProfile[itdc];
      fProfile[elem].i = elem;
    }
  }
  
  // Loop on planes:
  for(UInt_t p=0;p<fPlaneProfile.size();p++) {
    if(fPlaneProfile[p]>fBigPlaneE) {
      fBigPlaneE = fPlaneProfile[p];
      fBigPlane  = p;
    }
    if(fPlaneProfile[p]>0) fPlanes++;
  }

  // Loop on neutrinos:
  for(UInt_t inu=0;inu<inNeutrinos.size();inu++) {
    for(stripItr=nuStrips[inu].begin(); stripItr!=nuStrips[inu].end(); stripItr++) {
      fNuMatch_strips[inu] += 1;
      if(stripItr->first.GetPlane()<121) 
        fNuMatch_calstrips[inu] += 1;
    }
  }

  
  fBestNeutrino1 = -1;
  fBestNeutrino2 = -1;
  double best_size = 0;
  double best_size2 = 0;
  for(UInt_t inu=0;inu<inNeutrinos.size();inu++) {
    if(fNuMatch_calt[inu]>best_size) {
      fBestNeutrino1 = inu;
      best_size = fNuMatch_calt[inu];
    }
  }
   
  for(UInt_t inu=0;inu<inNeutrinos.size();inu++) {
    if(
       (fNuMatch_calt[inu]>best_size2)
       &&((int)inu!=fBestNeutrino1)
       ) {
      fBestNeutrino2 = inu;
      best_size2 = fNuMatch_calt[inu];
    }
  }
}

---