Module/NtpVtxFinder.cxx

Go to the documentation of this file.

//
//                          NtpVtxFinder.cxx  -  description
//                             -------------------
//    begin                : Wed Oct 19 2005
//    copyright            : (C) 2005 by Joshua Boehm
//    email                : boehm@physics.harvard.edu
//    last modified        : October 19, 2005
//    NtpVtxFinder is a reconstruction class used to determine neutrino 
//     interaction vertices.

#include <Validity/VldContext.h>
#include <UgliGeometry/UgliGeomHandle.h>
#include <RecoBase/CandStripHandle.h>
#include <CandTrackSR/CandTrackSRHandle.h>
#include <CandFitTrackSR/CandFitTrackSRHandle.h>
#include <Plex/PlexPlaneId.h>
#include "MessageService/MsgService.h"
#include "Conventions/Detector.h"

#include "AnalysisNtuples/Module/ANtpRecoNtpManipulator.h"
#include "StandardNtuple/NtpStRecord.h"
#include "VertexFinder/Module/NtpVtxFinder.h"

CVSID("$Id: NtpVtxFinder.cxx,v 1.3 2007/11/11 04:58:36 rhatcher Exp $");

//Only Constructor
NtpVtxFinder::NtpVtxFinder(int event, RecRecordImp<RecCandHeader> *srobj)
{
   kUVPlanesBefore = 4;
   kUVPlanesAfter = 5;
   kNoiseCut = 2.0;

   fTotalEnergy = 0;

   fX = fY = fU = fV = fPlane = 0;

   fEventNo = event;
   fSt = dynamic_cast<NtpStRecord *>(srobj);
   fSrObj = srobj;
   if(fSt)
   {
       int detType = 
             fSt->GetHeader().GetVldContext().GetDetector();

       if(detType == Detector::kFar)
         kNPlanes = 486;

       if(detType == Detector::kNear)
         kNPlanes = 285;
  }
}

//Nothing to specifically delete
NtpVtxFinder::~NtpVtxFinder(){
}


Int_t NtpVtxFinder::FindVertex()
{
  if(!fSt || !fSrObj) return -1;
 
  ANtpRecoNtpManipulator ntpManipulator(fSt);
  ntpManipulator.SetPrimaryTrackCriteria(0,1,0);
      // nplanes, length, total pulse height
      //get the primary track for the event - if no track is present it
  ANtpEventManipulator * ntpEventManipulator =
          ntpManipulator.GetEventManipulator();

  ntpEventManipulator->SetEventInSnarl(fEventNo);
  fEvent = ntpEventManipulator->GetEvent();
  fTrack = ntpEventManipulator->GetPrimaryTrack();
                                           
  Int_t retval = FindVtxPlane();
  if(retval < -1) return retval; 
  FindVtxUV();
  
  if(CheckTrackInfo()) {
      AcceptTrackInfo();
      MSG("VertexFinder",Msg::kDebug)<<"Taking trk vtx"<<endl;
  }else{
    MSG("VertexFinder",Msg::kDebug)<<"Rejected trk vtx"<<endl;
  }

  //One Plane subtraction  
  if(fPlane>0){
    fPlane--;
  }
    
  if(fPlane >=0){
    VldContext vc=fSrObj->GetHeader().GetVldContext(); 
    UgliGeomHandle ugh(vc);
    PlexPlaneId plnid(vc.GetDetector(),fPlane,true);
    UgliSteelPlnHandle usph = ugh.GetSteelPlnHandle(plnid);
    fZ = usph.GetZ0();
                                                                               
    fX = (sqrt(2.0)/2.0)*(fU - fV);
    fY = (sqrt(2.0)/2.0)*(fU + fV);
  }
  else return -1;
  
  return 1;
}


Int_t NtpVtxFinder::FindVtxPlane()
{
  if(!fSt) return -1;
  Int_t retval = 1;
                                                                             
  fTotalEnergy = 0;
  Float_t trimEnergy = 0;
                                                                             
  vector<Float_t> longEnergy(kNPlanes, 0);
                             
  for(int i=0;i<fEvent->nstrip;i++){
      Int_t index = SntpHelpers::GetStripIndex(i,fEvent);
      NtpSRStrip *strip = SntpHelpers::GetStrip(index,fSrObj);

      Double_t charge = strip->ph0.pe+strip->ph1.pe;

      if(charge > kNoiseCut){
           longEnergy[strip->plane] += charge;
           trimEnergy += charge;
     }
     fTotalEnergy += charge;
  }

  MSG("VertexFinder",Msg::kDebug)
          <<"Total Energy in the event: "<<fTotalEnergy<<"\n"
          <<"Over Threshold Energy in the event: "<<trimEnergy<<endl;
    
  if(trimEnergy < 40){
    MSG("VertexFinder",Msg::kDebug)
           <<"Returning due to insufficient energy"<<endl;
    return -5;
  }

  VtxClusterList* vcl = new VtxClusterList(longEnergy, kNPlanes, kNoiseCut);
  fPlane = vcl->Process();
//  MSG("VertexFinder",Msg::kDebug)<<"Reporting "<<vcl->Report();
  vcl->Clear();
  delete vcl;                                                            

  return retval;
}


Int_t NtpVtxFinder::FindVtxUV()
{
  Int_t retval = 1;

  Int_t nbins = 200;
  Float_t start = -5.0;
  Float_t end = 5.0;
  
  TH1F *tenestu = new TH1F("tenestu","TE (U view)", nbins, start, end);
  TH1F *tenestv = new TH1F("tenestv","TE (V view)", nbins, start, end);
  TH1F *tenestuAll = new TH1F("tenestu all","TE(U)", nbins, start, end);
  TH1F *tenestvAll = new TH1F("tenestv all","TE(V)", nbins, start, end);
  TH1F *tenestuBackup = new TH1F("tenestu backup","TE(U)", nbins, start, end);
  TH1F *tenestvBackup = new TH1F("tenestv backup","TE(V)", nbins, start, end);

  bool inRange = false;
  Float_t delta = -1;

  for(int i=0;i<fEvent->nstrip;i++){
        Int_t index = SntpHelpers::GetStripIndex(i,fEvent);
        NtpSRStrip *strip = SntpHelpers::GetStrip(index,fSrObj);
                                                                                
        Double_t charge = strip->ph0.pe+strip->ph1.pe;
                                                                                
        MSG("VertexFinder",Msg::kVerbose)<<"charge deposit:"
                 <<charge<<"    "<<strip->plane<<"    "<<endl;

       if(charge > kNoiseCut){
         delta = strip->plane-fPlane;
         inRange = false;
         if((delta <= 0 && fabs(delta) < kUVPlanesBefore) ||
             (delta >= 0 && fabs(delta) < kUVPlanesAfter)){
                inRange = true;
         }
                                                                       
       if(strip->planeview==PlaneView::kU){
          tenestuAll->Fill(strip->tpos,charge);
          if(inRange)   tenestu->Fill(strip->tpos,charge);
       }
      else if(strip->planeview==PlaneView::kV){
        tenestvAll->Fill(strip->tpos,charge);
        if(inRange)  tenestv->Fill(strip->tpos,charge);
      }
    }
                                                                    
    if(strip->planeview==PlaneView::kU){
       tenestuBackup->Fill(strip->tpos,charge);
    }
    else if(strip->planeview==PlaneView::kV){
       tenestvBackup->Fill(strip->tpos,charge);
    }          

  }//end of strip loop 

  //Fill parameters for U and V direction
  fU = tenestu->GetMean();
  fV = tenestv->GetMean(); 

  MSG("VertexFinder",Msg::kDebug)
          <<tenestu->GetMaximum()<<"    "<<tenestuAll->GetMaximum()<<endl;
  if(tenestu->GetMaximum() < 2) fU = tenestuAll->GetMean();
  if(tenestv->GetMaximum() < 2) fV = tenestvAll->GetMean();
  if(tenestuAll->GetMaximum() < 2) fU = tenestuBackup->GetMean();
  if(tenestvAll->GetMaximum() < 2) fV = tenestvBackup->GetMean();

  MSG("VertexFinder",Msg::kDebug)<<"plane guess "<<fPlane
          <<" U/V: "<<fU<<"   "<<fV << endl;

  //clean up after myself
  delete tenestu;
  delete tenestv;
  delete tenestuAll;
  delete tenestvAll;
  delete tenestuBackup;
  delete tenestvBackup;

  return retval;
}

Int_t NtpVtxFinder::CheckTrackInfo()
{
  Int_t pass = 0;
  Int_t plane = 0;
  if(fTrack){
     pass = fTrack->fit.pass;
     Float_t nTrackPlane = static_cast<Float_t>(fTrack->plane.n);
     Float_t nTrackLike = static_cast<Float_t>(fTrack->plane.ntrklike);
     plane = TMath::Min(fTrack->plane.beg, fTrack->plane.end);
       
    if(plane < 5 && nTrackPlane > 10) return 1;   //3 and 0.92

    if(pass==1){//  removed for pass 1
       if(nTrackPlane>45) return 1;
                                                                       
       Float_t ratio=0;
       if(nTrackPlane>0)  ratio = nTrackLike/nTrackPlane;

       if(nTrackPlane > 10 && ratio > 0.9)
          return 1;

    }
  }

  return 0;
}

Int_t NtpVtxFinder::AcceptTrackInfo()
{
   if(fTrack == 0) {
      MSG("VertexFinder",Msg::kError)<<"No track present, but trying to "
             <<"accept track info anyway - should NOT be possible"<<endl;
   } 

   Int_t plane = TMath::Min(fTrack->plane.beg, fTrack->plane.end);
   fU = fTrack->vtx.u;
   fV = fTrack->vtx.v;
   fPlane = plane;

   return 1;
}

---