HdAxisFinder.cc

Go to the documentation of this file.

#include "HdAxisFinder.h"

#include <det/Detector.h>
#include <fdet/FDetector.h>
#include <sdet/SDetector.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerFRecData.h>

#include <fdet/Channel.h>
#include <fdet/Eye.h>
#include <fdet/Pixel.h>

#include <fevt/Eye.h>
#include <fevt/EyeHeader.h>
#include <fevt/EyeRecData.h>
#include <fevt/FEvent.h>
#include <fevt/Pixel.h>
#include <fevt/PixelRecData.h>
#include <fevt/Telescope.h>

#include <fwk/CoordinateSystemRegistry.h>
#include <fwk/CentralConfig.h>

#include <sevt/SEvent.h>                // define type `struct sevt::SEvent'
#include <sevt/StationRecData.h>        // to get the station RecData
#include <sevt/Station.h>

#include <utl/ErrorLogger.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/Reader.h>
#include <utl/Transformation.h>

#include <TMinuit.h>

#include <iostream>

using namespace det;
using namespace std;
using namespace fwk;
using namespace utl;
using namespace fevt;
using namespace sevt;
using namespace evt;
using namespace HdAxisFinderUU;

fevt::Eye * HdAxisFinder::fCurEye = 0;
sevt::SEvent * HdAxisFinder::fCurSEvent = 0 ;
double HdAxisFinder::fChi2 = 0;
double HdAxisFinder::fMyChiSquare = 0;
unsigned int HdAxisFinder::fNDof = 0;

utl::Vector HdAxisFinder::fSDP = utl::Vector();

VModule::ResultFlag HdAxisFinder::Init(){
// Note: I need to check all eyes were recon.
// ie. may be an stereo event reconstructs only one of the eyes...
// Anyway, use the same cuts as the axis guesser:
  CentralConfig* cc = CentralConfig::GetInstance(); 
 
  Branch topB = 
    cc->GetTopBranch("HdAxisFinder");

  if (!topB) {
    ERROR("Could not find branch HdAxisFinder");
    return eFailure;
  }


  Branch minPixB =  topB.GetChild("MinPixels");
  if (!minPixB) {
    ERROR("Could not find branch MinPixels");
    return eFailure;
  }
  minPixB.GetData(fMinPixels);

  
  return eSuccess;
}

VModule::ResultFlag HdAxisFinder::Run(evt::Event& event){
//  INFO("HdAxisFinder::Run()");

  fMinuitFailed = false;

  fevt::FEvent& fevent = event.GetFEvent();
  sevt::SEvent& sevent= event.GetSEvent();

  fCurSEvent = &sevent; // acrobacy for TMinuit

  fevt::FEvent::EyeIterator eyeiter;

/* I'm NOT filling the "event" core, so forget about this:
  // Before anything, check if it's an stereo event
  int MaxNPixels = 0, numEyes = 0;
  int fHotEyeId = 0;
  for (eyeiter= fevent.EyesBegin(); 
       eyeiter != fevent.EyesEnd(); 
       ++eyeiter){

    fevt::Eye& eye = *eyeiter;
    fevt::EyeRecData& eyerecdata = eye.GetRecData();    
    int numPix = eyerecdata.GetNTimeFitPixels() ;

//    cout << "Eye #" << eye.GetId() << " has " << numPix << " pixels." << endl ;

    if (numPix > MaxNPixels) {
            MaxNPixels = numPix ;
            fHotEyeId = eye.GetId();
    }
    ++numEyes;
  }

  bool Stereo = false;
  if (numEyes > 1) Stereo = true ;
  // End of checking on stereo event
 */

  for (eyeiter= fevent.EyesBegin(); 
       eyeiter != fevent.EyesEnd(); 
       ++eyeiter){

    fevt::Eye& eye = *eyeiter;

    fevt::EyeRecData& eyerecdata = eye.GetRecData();    

    if(eyerecdata.GetNTimeFitPixels() < fMinPixels) continue ;

    fSDP = eyerecdata.GetSDP();

    fCurEye = &eye; // acrobacy for TMinuit

    this->FillPoints(eye);

    fChi0First = eyerecdata.GetChiZero() ;
    fRpFirst = eyerecdata.GetRp() ;
    fT0First = eyerecdata.GetTZero() ;
  
    this->FindAxis(eye);

// change this. If Minuit failed, skip all the following
/*
    if (fMinuitFailed && (!Stereo || (Stereo && (eye.GetId()==fHotEyeId)))) 
      return eContinueLoop;
//    if (fMinuitFailed) return eContinueLoop;
 */
    if (fMinuitFailed) continue;
// if it's an stereo event and both eyes failed
// the event will be cut next in the Final Cut (or it should...)
// right, but what if this happens in the iteration AFTER the final cut...

    eyerecdata.SetChiZero(fChi0,fChi0Error);
    eyerecdata.SetTZero(fT0,fT0Error);
    eyerecdata.SetRp(fRp,fRpError);
// Do this later!
//    eyerecdata.SetTimeFitChiSquare(fChi2,fNDof);

    // calculate axis parameters. 
    CoordinateSystemPtr eyeCS =
      det::Detector::GetInstance().GetFDetector().GetEye(eye).GetEyeCoordinateSystem();
    
    Vector  vertical(0,0,1,eyeCS);
    Vector horizontalInSDP = cross(fSDP,vertical); // horizontal
    horizontalInSDP.Normalize();

    // check the horizontal points outward
    if ( horizontalInSDP.GetPhi(eyeCS) < 0 || 
          horizontalInSDP.GetPhi(eyeCS) > kPi) {
      fSDP *= -1;
      horizontalInSDP *= -1;
    }

    Transformation rot (Transformation::Rotation(-fChi0,fSDP, eyeCS));

    Vector axis (rot* horizontalInSDP);             // apply rotation    
    if (axis.GetTheta(eyeCS) > (kPi/2)) axis *= -1; // reflect if necessary
    axis.Normalize();

    CoordinateSystemPtr siteCS =
      det::Detector::GetInstance().GetSiteCoordinateSystem();
    
    /* Minimize print outs...
    cout << "(Hybrid) Axis(Theta,Phi)(siteCS) : " << axis.GetTheta( siteCS) / degree
         << "," << axis.GetPhi( siteCS) /degree << " [deg]" <<endl; ;
    cout << "(Hybrid) Axis(Theta,Phi)(eyeCS) : " << axis.GetTheta( eyeCS) / degree
         << "," << axis.GetPhi( eyeCS) /degree << " [deg]" <<endl; ;
     */

    Vector core_eye_vec = fRp / sin( kPi - fChi0) * horizontalInSDP;
    
    Point eyeposition = 
      det::Detector::GetInstance().GetFDetector().GetEye(eye).GetPosition();

    // Calc. the core at the stations height! and NOT at the Fd height.
    // Note: still each eye keeps its right Rp, Chi0 and T0
    //
    // loop on stations
    const sdet::SDetector& theSDetector = det::Detector::GetInstance().GetSDetector();

    SEvent::StationIterator iter;
    double station_level = 0.0 ;
    int ntanks = 0 ;
    for (SEvent::ConstCandidateStationIterator iter = sevent.CandidateStationsBegin();
         iter != sevent.CandidateStationsEnd(); ++iter) {

      ++ntanks ; // Actual candidate tanks!
      const sevt::Station& station = *iter;
      const sdet::Station& currentSt = theSDetector.GetStation(station);
      const Point& SdPos = currentSt.GetPosition();
      // cout << " Station " << iter->GetId() << " has elevation: " << SdPos.GetZ(eyeCS) << " m" << endl;

      station_level += SdPos.GetZ(eyeCS) ;
    }

    if (ntanks > 0)
      station_level /= ntanks;
    else 
      station_level =  core_eye_vec.GetZ(eyeCS) ;
// This is wrong: (Remember I might leave tanks there... but they're just not "candidates")
//    station_level /= sevent.GetNumberOfStations() ;
//    cout << "Core level = " << station_level << " m." << endl ;

    if (station_level==0) station_level = core_eye_vec.GetZ(eyeCS) ;
    double coord_along_axis = (station_level - core_eye_vec.GetZ(eyeCS))/axis.GetZ(eyeCS) ;

    Point core = eyeposition + core_eye_vec + coord_along_axis * axis;

      /*
    cout << "(Hybrid) Core (x,y)(eyeCS) " << core.GetX(eyeCS)<<","
         <<  core.GetY(eyeCS)<< " [m]" << endl;

    cout << "(Hybrid) Core (x,y)(siteCS) " << core.GetX(siteCS)<<","
         <<  core.GetY(siteCS)<< " [m]" << endl;
       */
    
// Set all the EYE rec data:
    if (!eyerecdata.HasFRecShower()) eyerecdata.MakeFRecShower();
    ShowerFRecData& eyerecshower = eyerecdata.GetFRecShower();
    eyerecshower.SetAxis(axis);
    eyerecshower.SetCorePosition(core);
    eyerecshower.SetCoreTime(fCurEye->GetHeader().GetTimeStamp() - 
                             TimeInterval(1000*100*ns), fRp, fChi0, fT0);

/* I moved these lines to the Final Cut because they are being executed in the first pass to the Axis Finder
So, all tanks are there... I guess...
// Add these lines for the RecDataProducer:
    for(iter=sevent.StationsBegin(); iter != sevent.StationsEnd(); ++iter)
    {
      if (iter->IsCandidate()){
        eyerecshower.AddStationId( iter->GetId() );
      }
    }
 */


// Fill out the "real" chi square!
/*
  cout << "Compare Chi square: " << fChi2 << " (old) vs. " << fMyChiSquare << " (new) and ndof = " 
       << fNDof << " (old) vs. " << ntanks+eyerecdata.GetNTimeFitPixels()-3 << " (new) " << endl ;
 */
    eyerecdata.SetTimeFitChiSquare(fMyChiSquare,ntanks+eyerecdata.GetNTimeFitPixels()-3);

/* This crap doesn't really make too much sense.
 * Don't fill this no more!

    // Note: For the EVENT core it uses ONLY the "Hot" eye, ie the one with more pixels...
    // when it's an stereo event.
    if (!Stereo || (fHotEyeId == eye.GetId()) )
    {
      if (!event.HasRecShower()) event.MakeRecShower();
      ShowerRecData& recshower = event.GetRecShower();
      recshower.SetAxis(axis);
      recshower.SetCorePosition(core);
    }
 */

  }// for eyeiter

  return eSuccess;
}


void
HdAxisFinder::FindAxis(fevt::Eye& /*eye*/)
{
  const int npar = 3; // number of parameters
  
  TMinuit theMinuit(npar); // 3 parameter fit
  // this line was leaking memory! :O(
//  TMinuit * theMinuit = new TMinuit (npar); // 3 parameter fit
  
  theMinuit.SetPrintLevel(-1); // minuit verbosity
  theMinuit.SetFCN(HdAxisFinder::MinuitFitFunc);
   
  double arglist[npar];  // arguments to pass to the Minuit interpreter
  int ierflag =0;
  
  arglist[0]=1; // [up]
  arglist[1]=1; //  unused
  theMinuit.mnexcm("SET ERR", arglist,1,ierflag);
  
  if (ierflag) ERROR ("Minuit SET ERR failed");
  
  static double vstart[npar]; // inital point
  static double step[npar];   // step size

  vstart[0] = fChi0First;
  vstart[1] = fT0First;
  vstart[2] = fRpFirst;
  
  step[0]=0.001;
  step[1]=0.001;
  step[2]=0.001;
  
  theMinuit.mnparm(0,"Chi0",vstart[0], step[0],0,0,ierflag);    
  theMinuit.mnparm(1,"T0",  vstart[1], step[1],0,0,ierflag);
  theMinuit.mnparm(2,"Rp",  vstart[2], step[2],0,0,ierflag);
  
  arglist[0]=5000;  // [maxcalls]
  arglist[1]=1;    // [tolerance]
  theMinuit.mnexcm("MIGRAD", arglist ,2,ierflag);

//    theMinuit.mncomd("SHOW COV", ierflag);
//    theMinuit.mncomd("SHOW COR", ierflag);

  
  if (ierflag) {
    ERROR ("Minuit MIGRAD failed");
    fMinuitFailed=true;
  }
  theMinuit.GetParameter(0,fChi0,fChi0Error);
  theMinuit.GetParameter(1,fT0,fT0Error);
  theMinuit.GetParameter(2,fRp,fRpError);
  
  /*
  cout << "(Hybrid) TO   "  << fT0  << " +/- " << fT0Error << endl;
  cout << "(Hybrid) ChiO "  << fChi0/degree << " +/- " << fChi0Error/degree << endl;
  cout << "(Hybrid) Rp   "  << fRp   << " +/- " << fRpError << endl;
   */
  
}

void HdAxisFinder::FillPoints(fevt::Eye& eye){

  fevt::EyeRecData& recdata = eye.GetRecData();

  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();
  const fdet::Eye& deteye = detFD.GetEye(eye);

  // do all calculation in the eye reference system
  CoordinateSystemPtr CS= deteye.GetEyeCoordinateSystem();

  // old stuff:
//  const double timebinsize = fdet::Channel::GetFADCBinSize();
//  // now this stuf is moved into the pixel loop

//there's no need to do this for EVERY pixel!
  fSDP.TransformTo(CS);
  fSDP.Normalize();

  Vector vertical(0,0,1,CS,Vector::kCartesian);
  Vector horizontalWithinSDP = cross(fSDP, vertical);

  // check that it is the one looking outward
  if ( (horizontalWithinSDP.GetPhi(CS) < 0) || (horizontalWithinSDP.GetPhi(CS) > kPi) ) //is looking outward!
  {
    horizontalWithinSDP*=-1;
    fSDP *= -1;
  }
  horizontalWithinSDP.Normalize();

  // iterate over pixels used for time fit, fill chi_0 and t_0
  fevt::EyeRecData::PixelIterator pixeliter;
   
  for (pixeliter = recdata.TimeFitPixelsBegin();
       pixeliter != recdata.TimeFitPixelsEnd(); 
       ++pixeliter){

     fevt::Pixel& evtpixel = *pixeliter; 
     const fdet::Pixel& detpixel = detFD.GetPixel(evtpixel);
     const fdet::Channel& detChannel = detFD.GetChannel(evtpixel);
  //const double timebinsize = fdet::Channel::GetFADCBinSize();
  const double timebinsize = detChannel.GetFADCBinSize();

     PixelRecData& pixrec = evtpixel.GetRecData();

     // time correction for different eyes
     int telid = evtpixel.GetTelescopeId();
     unsigned int timeoffset = eye.GetTelescope(telid).GetTimeOffset();
     
     double t_i = pixrec.GetCentroid()  ;

     t_i *= timebinsize;  // convert to nanosecond
     t_i += timeoffset;   // add offset
     
/* New stuff. Replace centroid error (too small) with "my" function
 * that parametrizes the funny steps of GAP 2003-072
     double t_i_Err = evtpixel.GetRecData().GetCentroidError()   ; 
     t_i_Err *= timebinsize;
 */
     double charge   = pixrec.GetTotalCharge() ;
     double noise    = pixrec.GetRMS() ;
     double duration = (pixrec.GetPulseStop() - pixrec.GetPulseStart())* timebinsize ;
     double ratio    = (charge/noise)/(duration/timebinsize) ;
//     double t_i_Err  = (200*exp(-log(4.)*duration/500)+50) * (4*exp(-log(8.)*ratio/7)+0.5) ;
// Implement latest parameterization
     double t_i_Err  = (180*exp(-duration/280)+50) * (3.6*exp(-ratio/3)+0.4) ;
// end of new parametrization of errors

     Vector pixeldir = detpixel.GetDirection();
     pixeldir.TransformTo(CS);

     // subtract component parallel to sdp
     Vector chi_i_vector = pixeldir - (fSDP * pixeldir) *  fSDP;
     chi_i_vector.TransformTo(CS);
     chi_i_vector.Normalize();

     double chi_i = Angle(chi_i_vector,horizontalWithinSDP); 
//cout << "fit values: chi_i = " << chi_i/deg << " ti = " << t_i << endl ;
//cout << "Pixel #" << evtpixel.GetId() << " fit values: chi_i = " << chi_i/deg << " ti = " << t_i 
//     << "; t_i error = " << t_i_Err << " [ns]" << endl ;
     pixrec.SetChi_i(chi_i);
     pixrec.SetT_i(t_i,t_i_Err);
    
  } // for pixel
    
}


void
HdAxisFinder::MinuitFitFunc(int& /*npar*/, double* /*gin*/, double& f,
                            double* par, int /*iflag*/)
{  
  fevt::EyeRecData& recdata = fCurEye->GetRecData();
  
  fMyChiSquare = 0.0;

  double chisq =0;
  double sdchisq = 0;
  unsigned int ndof =0;

  double& chi_0 = par[0];
  double& t_0   = par[1] ;
  double& r_p   = par[2] ;

  EyeRecData::ConstPixelIterator iter;
  double totalcharge = 0. ;
  double max_error = 0.0 ;

  for (iter= recdata.TimeFitPixelsBegin(); 
       iter!=recdata.TimeFitPixelsEnd(); ++iter){  
  
    const Pixel& pixel = *iter;
    const PixelRecData& pixrec = pixel.GetRecData();
    
    double t_i       = pixrec.GetT_i().GetNanoSecond();
    double t_i_Err   = pixrec.GetT_iError().GetInterval();
//    cout << "(In Fit) Pixel #" << pixel.GetId() << " t_i = " << t_i << " error = " << t_i_Err << " [ns]" << endl ;
    if (t_i_Err > max_error) max_error = t_i_Err ; // keep it to weight the tank-core distance error

    double chi_i     = pixrec.GetChi_i();
//    double chi_i_Err = 0;

    double t_expected = 
      t_0 + r_p/kSpeedOfLight *  std::tan((chi_0 - chi_i) /2. );

    double tmp = t_i - t_expected  ; 

    double charge = pixrec.GetTotalCharge() ;
    totalcharge += charge ;
/*
 */
// Try just the time uncertainty for now...
//    chisq += tmp*tmp/t_i_Err /t_i_Err ;
    chisq += tmp*tmp/ t_i_Err /  t_i_Err * charge;
    fMyChiSquare += tmp*tmp/t_i_Err/t_i_Err ;
    ++ndof;

  }// for iter

  int npixels = recdata.GetNTimeFitPixels() ; 
  chisq *= npixels ; 
//  chisq /= recdata.GetNTimeFitPixels() ; 
  chisq /= totalcharge ;

//cout << "fd chi sq. = " << chisq << endl ;
  // Get all the Eye info: Rp vector, axis, etc.
  CoordinateSystemPtr CS = 
          Detector::GetInstance().GetFDetector().GetEye(*fCurEye).GetLocalCoordinateSystem();

  fSDP.TransformTo(CS) ;
  Vector vertical(0,0,1,CS,Vector::kCartesian) ;
  AxialVector horizontal = cross(fSDP, vertical) ;

  CoordinateSystemPtr EyeCS = 
          Detector::GetInstance().GetFDetector().GetEye(*fCurEye).GetEyeCoordinateSystem();

  if ( (horizontal.GetPhi(EyeCS) < 0) || (horizontal.GetPhi(EyeCS) > kPi) ) //is looking outward!
  {
     fSDP *= -1 ;     // I need this to get the rotation right => the axis
     horizontal *= -1 ;
  }
  horizontal.Normalize() ;

  Transformation rot (Transformation::Rotation(-chi_0, fSDP, CS));

  Vector axis (rot*horizontal);             // apply rotation
  if (axis.GetTheta(CS) > (kPi/2)) axis *= -1; // reflect if necessary
  axis.Normalize();

  double theta = Angle(axis,vertical); 

  Vector FdCore = r_p / sin(kPi - chi_0) * horizontal;
  FdCore.TransformTo(CS) ;

/*
    cout << "(Hybrid) Axis(Theta,Phi)(eyeCS) : " << axis.GetTheta(CS) / degree
         << "," << axis.GetPhi(CS) /degree << " [deg]" <<endl; ;

    cout << "(FD) Core (x,y)(eyeCS) " << FdCore.GetX(CS)<<","
         <<  FdCore.GetY(CS)<< "," << FdCore.GetZ(CS) << " [m]" << endl;
 */

  double EyeTriggerTime = fCurEye->GetHeader().GetTimeStamp().GetGPSNanoSecond() - 1e5 ;

  const sdet::SDetector& theSDetector = det::Detector::GetInstance().GetSDetector();

  // loop on stations now
  totalcharge = 0. ;
  //double max_signal = 0.0 ;
  //double min_distance = 0.0 ;

  double baryWeight, baryWeightSum = 0.0, baryX = 0.0, baryY = 0.0, baryZ = 0.0 ;

  for (SEvent::ConstCandidateStationIterator Stiter = fCurSEvent->CandidateStationsBegin();
       Stiter != fCurSEvent->CandidateStationsEnd(); ++Stiter) {

    const sevt::Station& station = *Stiter;
    const sdet::Station& currentSt = theSDetector.GetStation(station);
    const Point& SdPos = currentSt.GetPosition();
    //cout << "Tank position: " << SdPos.GetX(CS)<<"," <<  SdPos.GetY(CS)<< "," << SdPos.GetZ(CS) << " [m]" << endl;
    const Vector TankToEye(SdPos.GetX(CS), SdPos.GetY(CS), SdPos.GetZ(CS), CS);
    const Vector FdCoreToTank = TankToEye - FdCore;
    /*
    cout << "(fit) Vector fd core to tank X = " << FdCoreToTank.GetX(CS) 
         << ", Y = " << FdCoreToTank.GetY(CS) <<", Z = " << FdCoreToTank.GetZ(CS) << " [m]" <<endl;

    cout << "Distance core-to-tank: " << FdCoreToTank.GetMag() << " [m]" << endl;
    */

    const double station_level = SdPos.GetZ(CS); // for the distance to the core I need the "real" core.
    //cout << "Core level = " << station_level << " m." << endl ;

    const double coord_along_axis = (station_level - FdCore.GetZ(CS))/axis.GetZ(CS);
    // this "coordinate" is the delta(altitude)/cos(theta). (In the unit axis the Z coord. is equal to cos(theta)

    const Vector RealCore = FdCore + coord_along_axis * axis; //The plus sign depends only on the "direction" of the (altitude) subtraction.
    /*
    cout << "(real) Core (x,y)(eyeCS) " << RealCore.GetX(CS)<<","
         <<  RealCore.GetY(CS)<< "," << RealCore.GetZ(CS) << " [m]" << endl;
    */
    const Vector CoreToTank = TankToEye - RealCore;

    const double distance = CoreToTank.GetMag(); // distance between the tank and the real core
    //cout << "Distance to core " << distance << " [m]" << endl;

    // project tank into axis
    const double TankProj = FdCoreToTank * axis;
    //cout << "Distance (in axis) to core = " << TankProj << " [m] angle = " << acos(TankProj/FdCoreToTank.GetMag())/deg << endl;
    const Vector TankToAxis = FdCoreToTank - (TankProj * axis);
    const double d = TankToAxis.GetMag(); // distance from the tank to the shower axis

    // To get the time uncertainties copy
    // from GAP 2003-012:
    double sigma_t2 = (2800. + 0.0064*distance*distance)/3.5;

    const StationRecData& recdata = station.GetRecData();
    const double t_i = recdata.GetSignalStartTime().GetGPSNanoSecond();

    // time expected at the tank:
    double t_expected = t_0 - r_p/(kSpeedOfLight * tan(chi_0)); // time at the (Fd) core.
    t_expected += EyeTriggerTime;
    t_expected -= TankProj/kSpeedOfLight;
    //cout << " First time comp. = " << TankProj/kSpeedOfLight << " [ns]" << endl ;

    // Finally, the so expected curvature correction! MAM 01/10/05
    const double Rinit = 6500./cos(theta) ;// fixed for now
    // add asymmetry effect GAP 2003-108 => variable (adjustable) curvature radius
    // except for the sign of 2*R*r*cos(alpha) that I think they got wrong.
    const double alpha = Angle(axis,CoreToTank) ;
    const double Rcurv = sqrt(Rinit*Rinit + distance*distance - 2*Rinit*distance*cos(alpha));
    // -------------------------------------------------------^ this sign is + in their gap note...
    // Next step: fit the damn thing!
    const double gamma = 1/(2*Rcurv) ;
    t_expected += (d*d*gamma/kSpeedOfLight);
    //cout << " Second time comp. = " << d*d*gamma/kSpeedOfLight << " [ns]" << endl;

    // also add effect to variance
    const double crit_dist = 700./sqrt(cos(theta));
    sigma_t2 *= 1 + pow(d/crit_dist, 4);

    const double tmp = t_i - t_expected;
    /*
    cout << "\tTime at tank: " << int(t_i) << " - " << int(t_expected) << " (expected) = " << (tmp) << " [ns] +/- " 
         << int(sqrt(sigma_t2)) << " trigger time: " << int(EyeTriggerTime) << endl;
    */

    const double StTotSig = recdata.GetTotalSignal();
    /* not doing this anymore, now I use the barycenter for minimization
    if (StTotSig > max_signal) {
      min_distance = distance ;
      max_signal = StTotSig ;
    }
    */

    // Calc. barycenter (to cut on the hybrid core -> barycenter distance
    baryWeight = sqrt(StTotSig);
    baryWeightSum += baryWeight;
    baryX += baryWeight * SdPos.GetX(CS);
    baryY += baryWeight * SdPos.GetY(CS);
    baryZ += baryWeight * station_level;

    totalcharge += StTotSig;

    // also here. Try only time uncertainty for the chi^2
    //sdchisq += tmp*tmp/sigma_t2;
    sdchisq += tmp*tmp/sigma_t2 * StTotSig;
    fMyChiSquare += tmp*tmp/sigma_t2;

    ++ndof;

  }  // end of tanks loop

  const int ntanks = ndof - npixels;
  sdchisq *= ntanks;
  sdchisq /= totalcharge;
  //int ntanks = fCurSEvent->GetNumberOfStations();
  //sdchisq /= fCurSEvent->GetNumberOfStations();
  //sdchisq /= totalcharge;
  //cout << "sd chi sq. = " << sdchisq << endl;

  // Add one term for minimizing the core-tank distance for the one with the largest signal.
  // 10*m is too low, and they drive the fit always.
  // 100*m fixed is still not optimum (but much better than 10.)
  // Do this only on wimpy events! (MAM 03/25/05)
  //if (ntanks < 4) 
  // Now instead always minimize the distance between the barycenter and the hybrid core:
  //double dist_err = 100*m + max_error*kSpeedOfLight;
  baryX /= baryWeightSum;
  baryY /= baryWeightSum;
  baryZ /= baryWeightSum;
  const Vector BaryCenter(baryX, baryY, baryZ, CS);

  const double coord_along_axis = (baryZ - FdCore.GetZ(CS))/axis.GetZ(CS);
  const Vector RealCore = FdCore + coord_along_axis * axis; //The plus sign depends on the altitude subtraction.
  const Vector CoreToBary = BaryCenter - RealCore;
  const double min_distance = CoreToBary.GetMag(); // distance between the tank and the real core

  const double dist_err = 500*m;  // from Billoir's cut
  const double dist_chi = min_distance*min_distance / (dist_err*dist_err);
  ++ndof;

  sdchisq += dist_chi;
  //chisq *= (npixels + ntanks);
  // move on:

  chisq += sdchisq;

  //cout << "\t\t\tTotal chi^2 = " << chisq << " contribution from barycenter distance = " << dist_chi << endl;

  f = chisq;
  fChi2 = chisq;
  fNDof = ndof - 3;
}


VModule::ResultFlag
HdAxisFinder::Finish()
{
  return eSuccess;
}