VertexFinder.cxx

Go to the documentation of this file.

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

#include "VertexFinder.h"
#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 <cmath>

ClassImp(VertexFinder)
CVSID("$Id: VertexFinder.cxx,v 1.13 2008/09/02 22:28:49 boehm Exp $");

//Only Constructor
VertexFinder::VertexFinder(CandEventHandle * ceh, const VldContext * vldptr)
{
   kUVPlanesBefore = 4;
   kUVPlanesAfter = 5;
   kNoiseCut = 2.0;

   fTotalEnergy = 0;

   fCeh=ceh;
   fVld = vldptr;
   fX = fY = fU = fV = fPlane = 0;
   kNPlanes = 0;


   if(vldptr->GetDetector() == Detector::kFar)
       kNPlanes = 486;

   if(vldptr->GetDetector() == Detector::kNear)
       kNPlanes = 285;

   if(vldptr->GetDetector() == Detector::kCalDet)
      kNPlanes = 64;

   if(kNPlanes == 0){
      MSG("VertexFinder",Msg::kError)
         <<"Detector Type is  "<<vldptr->GetDetector()
         <<":  VertexFinder only designed for work in Near or Far Det"<<endl;
   } 
}

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


Int_t VertexFinder::FindVertex()
{
  if(!fCeh || kNPlanes == 0) return -1;
 
  Int_t retval = FindVtxPlane();
  if(retval < 0) 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){
    UgliGeomHandle ugh(*fVld);
    PlexPlaneId plnid(fVld->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 VertexFinder::FindVtxPlane()
{
  if(!fCeh) return -1;
  Int_t retval = 1;
                                                                             
  fTotalEnergy = 0;
  Float_t trimEnergy = 0;
                                                                             
  vector<Float_t> longEnergy(kNPlanes, 0);
  TIter stripItr(fCeh->GetDaughterIterator());
  while (CandStripHandle *strip = 
                 dynamic_cast<CandStripHandle*> (stripItr())) {
      Double_t charge = strip->GetCharge(CalDigitType::kPE);
      if(charge > kNoiseCut){
           longEnergy[strip->GetPlane()] += charge;
           trimEnergy += charge;
     }
     fTotalEnergy += charge;
  }
    
  if(trimEnergy == 0){
    MSG("VertexFinder",Msg::kError)
           <<"Zero Energy event, I declare that bad"<<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 VertexFinder::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;

  TIter stripItr(fCeh->GetDaughterIterator());    
                                                            
  stripItr.Reset();
  while (CandStripHandle *strip = dynamic_cast<CandStripHandle*>
         (stripItr())) {
    Double_t charge = strip->GetCharge(CalDigitType::kPE);
    MSG("VertexFinder",Msg::kDebug)<<"charge deposit:"
                 <<charge<<"    "<<strip->GetPlane()<<"    "<<endl;

    if(charge > kNoiseCut){
      delta = strip->GetPlane()-fPlane;
      inRange = false;
      if((delta <= 0 && fabs(delta) < kUVPlanesBefore) ||
           (delta >= 0 && fabs(delta) < kUVPlanesAfter)){
              inRange = true;
      }
                                                                       
      if(strip->GetPlaneView()==PlaneView::kU){
        tenestuAll->Fill(strip->GetTPos(),charge);
        if(inRange)   tenestu->Fill(strip->GetTPos(),charge);
      }
      else if(strip->GetPlaneView()==PlaneView::kV){
        tenestvAll->Fill(strip->GetTPos(),charge);
        if(inRange)  tenestv->Fill(strip->GetTPos(),charge);
      }
    }
                                                                    
    if(strip->GetPlaneView()==PlaneView::kU){
       tenestuBackup->Fill(strip->GetTPos(),charge);
    }
    else if(strip->GetPlaneView()==PlaneView::kV){
       tenestvBackup->Fill(strip->GetTPos(),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 VertexFinder::CheckTrackInfo()
{
  CandTrackHandle * track = fCeh->GetPrimaryTrack();
  CandFitTrackHandle *fittrack = 0;
                                                                                
  Int_t pass = 0;
  Int_t plane = 0;
  if(track){
     if (track->InheritsFrom("CandFitTrackHandle")) {
        fittrack = dynamic_cast<CandFitTrackHandle*>(track);
        pass = fittrack->GetPass();
     }

     Float_t nTrackPlane = static_cast<Float_t>(track->GetNPlane());
     Float_t nTrackLike = static_cast<Float_t>(track->GetNTrackPlane());
     plane = TMath::Min(track->GetBegPlane(), track->GetEndPlane());
       
     if(plane < 5 && nTrackPlane > 10) return 1;

     if(pass==1){
       if(track->GetNPlane()>45) return 1;
                                                                       
       Float_t ratio=0;
       if(nTrackPlane>0)  ratio = nTrackLike/nTrackPlane;

       if(track->GetNPlane() > 10  && ratio > 0.92)
          return 1;

     }
  }

  return 0;
}

Int_t VertexFinder::AcceptTrackInfo()
{
   CandTrackHandle * track = fCeh->GetPrimaryTrack();
   if(track == 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(track->GetBegPlane(), track->GetEndPlane());
   fU = track->GetVtxU();
   fV = track->GetVtxV();
   fPlane = plane;

   return 1;
}

---