FdSDPFinder.cc

Go to the documentation of this file.

#include <iostream>

#include <evt/Event.h>

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

#include <fwk/CentralConfig.h>

#include <utl/AxialVector.h>
#include <utl/ErrorLogger.h>
#include <utl/MathConstants.h>
#include <utl/Reader.h>

#include <det/Detector.h>

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

#include <TMinuit.h>

#include "FdSDPFinder.h"

using namespace std;
using namespace FdSDPFinderOG;
using namespace evt;
using namespace fevt;
using namespace fwk;
using namespace utl;
using namespace det;
using namespace fdet;


fevt::Eye * FdSDPFinder::fCurEye =0;
double FdSDPFinder::fChi2=0;
unsigned int FdSDPFinder::fNDof=0;

VModule::ResultFlag FdSDPFinder::Init(){

  Branch topB =
     CentralConfig::GetInstance()->GetTopBranch("FdSDPFinder");

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

  topB.GetChild("minPixels").GetData(fMinPixels);
  topB.GetChild("maxDistance").GetData(fMaxDistance);
  topB.GetChild("firstGuessFromCharge").GetData(fFirstGuessFromCharge);

  // info output
  ostringstream info;
  info << "\n"
          " Version: " << GetVersionInfo(VModule::eRevisionNumber) << "\n"
          " Parameters:\n"
          "       min no. pixels: " << fMinPixels << "\n"
          "         max distance: " << fMaxDistance << "\n"
          "       first guess from track bed "
       << (fFirstGuessFromCharge?"charge":"pixels") << '\n';

  INFO(info);

  return eSuccess;
}


VModule::ResultFlag FdSDPFinder::Run(evt::Event& event){

  if (!event.HasFEvent()) return eSuccess;

  FEvent&  fevent = event.GetFEvent();

  FEvent::EyeIterator eyeiter;

  unsigned int neyes =0; // number of eyes with enough pixels

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

    fevt::Eye& eye = *eyeiter;

    // skip virtual eyes
    if (det::Detector::GetInstance().GetFDetector().GetEye(eye).IsVirtual())
      continue;

    if (!eye.HasRecData())
      eye.MakeRecData();
    EyeRecData& eyerecdata = eye.GetRecData();

    ostringstream info;
    info << eyerecdata.GetNPulsedPixels() << " pixels had signal in Eye " << eye.GetId();
    INFO(info);

    if (eyerecdata.GetNPulsedPixels() < fMinPixels) {
      ostringstream warn;
      warn << " **** too few pixels with signal (minNPixels=" << fMinPixels << "). Skipping eye. **** ";
      INFO(warn);
      continue;
    }



    fIsolationLimit = 0;

    fGoodEvent = true;
    fSDPFirstGuess = FindSDPFirstGuess(eye);

    if (!fGoodEvent ) continue;

    // Finds a trial SDP.
    // Only pixels that
    // are not isolated in time or space
    // are considered.

    ReadmitPixel(eye);
    info.str("");
    info << eyerecdata.GetNSDPPixels()
         << " pixels remaining after removing noise";
    INFO(info);
    // Estimating Chi_i for each Pixel using the  SDP first guess

    if (eyerecdata.GetNSDPPixels() < 3 ) continue;

    neyes++;

    EstimateChi_i(eye);

    RemoveOutliers(eye);  // Reject pixels according a linear
                                // time fit.
    info.str("");
    info << eyerecdata.GetNSDPPixels()
         << " pixels remaining after removing outliers";
    INFO(info);

    FindSDP(eye);  // Fits a SDP



    if ( eyerecdata.GetNSDPPixels() < eyerecdata.GetNPulsedPixels() ) {

      unsigned int PixelsReadmited = true;
      while ( PixelsReadmited ) {
        PixelsReadmited = false;
        unsigned int nopixels = eyerecdata.GetNSDPPixels();
        // Now we use the fitted SDP
        // to check whether we need to readmit
        // (or remove) any pixel that was rejected (or admitted) when
        // using the TRIAL SDP
        fSDPFirstGuess = fSDP;
        ReadmitPixel(eye);
        int nReadmitted = eyerecdata.GetNSDPPixels() -  nopixels;

        info.str("");
        info << nReadmitted << " pixels readmitted.";
        INFO(info);

        unsigned int nopixels2 = eyerecdata.GetNSDPPixels();
        RemoveNoise(eye);           // Reject pixels far from SDP
        int nRemoved =  nopixels2 - eyerecdata.GetNSDPPixels();

        info.str("");
        info << nRemoved << " pixels removed.";
        INFO(info);

        if ( nReadmitted != 0 || nRemoved !=  0) {
          FindSDP(eye);  // Fits a new SDP
          PixelsReadmited = true;
        }
      }
    }
    eyerecdata.SetSDP(fSDP);
    eyerecdata.SetSDPThetaError(fSDPThetaError);
    eyerecdata.SetSDPPhiError(fSDPPhiError);
    eyerecdata.SetSDPCorrThetaPhi (fSDPCorrThetaPhi);
    eyerecdata.SetSDPFitChiSquare(fChi2,fNDof);

    //--------------------------------------------
    // Also fill the TelescopeRecData with the same values
    // since this module isn't able to deal with compound eyes
    // that consist of telescopes at multiple different positions
    for (fevt::Eye::TelescopeIterator telIt = eye.TelescopesBegin(ComponentSelector::eHasData);
         telIt != eye.TelescopesEnd(ComponentSelector::eHasData); ++telIt)
    {
      if (!telIt->HasRecData())
        telIt->MakeRecData();
      fevt::TelescopeRecData& telRecData = telIt->GetRecData();
      telRecData.SetSDP(fSDP);
      telRecData.SetSDPThetaError(fSDPThetaError);
      telRecData.SetSDPPhiError(fSDPPhiError);
      telRecData.SetSDPCorrThetaPhi(fSDPCorrThetaPhi);
      telRecData.SetSDPFitChiSquare(fChi2, fNDof);
    }
  }// for eyeiter

  if (!neyes) return eContinueLoop; //no eyes worth reconstructing, next event
  return eSuccess;
}

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

/*
 * \brief Find an inital guess of the SDP.
 *
 * Find the cross products of the pixels with a pulse
 * Pick the one that accomodates most trackbed pixels
 *
 */

utl::AxialVector FdSDPFinder::FindSDPFirstGuess(fevt::Eye& eye){

  const double trackbed = 2.0 * degree;
  const double trackbedCut = cos(trackbed)*cos(trackbed);

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

  CoordinateSystemPtr CS = deteye.GetEyeCoordinateSystem();

  // get the candidate normal vectors (cross products of all pixels in pairs(
  // and put them in a vector
  // order the cross product earlier cross later

  const EyeRecData& recdata = eye.GetRecData();

  const FDetector& fdetector = Detector::GetInstance().GetFDetector();

  // transform all directions to eyeCS to avoid useless coordinate transformations
  CoordinateSystemPtr eyeCS = fdetector.GetEye(eye).GetEyeCoordinateSystem();

  std::vector<double> candidatePixelTimes;
  std::vector<Vector > candidatePixelDirections;
  std::vector<std::vector<double> > candidatePixDirComponents;
  std::vector<double> candidatePixelCharge;

  // fill vectors of pixels, which are not isolated
  EyeRecData::ConstPixelIterator pixeliter;
  for (pixeliter = recdata.PulsedPixelsBegin();
       pixeliter != recdata.PulsedPixelsEnd();
       ++pixeliter) {

    if ( ! IsIsolated(*pixeliter,eye) ) {

      const fdet::Pixel& detpixel     = fdetector.GetPixel(*pixeliter);
      const fdet::Channel& detChannel = fdetector.GetChannel(*pixeliter);
      const int telid = pixeliter->GetTelescopeId();
      const unsigned int timeoffset = eye.GetTelescope(telid).GetTimeOffset();

      const double pixelTime = pixeliter->GetRecData().GetCentroid() +
                               double(timeoffset/detChannel.GetFADCBinSize());
      const double charge    = pixeliter->GetRecData().GetTotalCharge();

      const Vector& pixelDirection=detpixel.GetDirection();
      pixelDirection.TransformTo(eyeCS);

      vector<double> components;
      components.push_back(pixelDirection.GetX(eyeCS));
      components.push_back(pixelDirection.GetY(eyeCS));
      components.push_back(pixelDirection.GetZ(eyeCS));

      candidatePixDirComponents.push_back(components);
      candidatePixelDirections.push_back(pixelDirection);
      candidatePixelCharge.push_back(charge);
      candidatePixelTimes.push_back(pixelTime);

    }

  }


  if( candidatePixelTimes.size() < 3) {
    fGoodEvent = false;
    ostringstream warn;
    warn << " *************  This event has less than 3 good pixels. Skip!  ****************";
    WARNING(warn);
    return utl::AxialVector(0.,0.,1.,eyeCS);
  }

  // loop on pixels to fill candidate SDPs

  std::vector <utl::AxialVector> crossProducts;

  for (unsigned int i=0;i<candidatePixelTimes.size(); i++ ) {

    const double t_1   = candidatePixelTimes[i];
    const Vector& dir1 = candidatePixelDirections[i];

    for (unsigned int j = i; j < candidatePixelTimes.size(); ++j) {

      if (j != i) {

        const double t_2   = candidatePixelTimes[j];
        const Vector& dir2 = candidatePixelDirections[j];

        // order the cross product : earlier cross later
        utl::AxialVector crossprod = (t_1 < t_2) ?
          cross(dir1,dir2) :
          cross(dir2,dir1);
        crossprod.Normalize();

        if ( crossprod.GetMag() != 0. )
          crossProducts.push_back(crossprod);

      }
    }
  }

  // look for best trial between candidate vectors:
  // it will be the one that accomodates most trackbed pixels


  ostringstream info;
  info << " track bed maximization (N=" << crossProducts.size() << "*"
       << candidatePixelTimes.size() << ")="
       << crossProducts.size()*candidatePixelTimes.size();
  INFO(info);

  double maxTrackBed=-1.;
  utl::AxialVector trialSDP;

  for(unsigned int k=0;k<crossProducts.size(); k++) {

    double thisTrackBed=0.;

    const utl::Vector& tmpSDP = crossProducts[k];

    const double x1 = tmpSDP.GetX(eyeCS);
    const double y1 = tmpSDP.GetY(eyeCS);
    const double z1 = tmpSDP.GetZ(eyeCS);

    for (unsigned int i=0;i<candidatePixelTimes.size(); i++ ) {

      const double * pixDir = &(candidatePixDirComponents[i].front());
      const double cosTheta = x1*(*(pixDir))
                             +y1*(*(pixDir+1))
                             +z1*(*(pixDir+2));
      const double sin2Theta = 1. - cosTheta*cosTheta;

      // --> original cut:
      //  const double d = abs (kPi/2 - acos(cosTheta));
      //  if ( d <trackbed)

      // --> same thing without acos():
      if ( sin2Theta > trackbedCut ) {
        if ( fFirstGuessFromCharge )
          thisTrackBed+=candidatePixelCharge[i];
        else
          thisTrackBed+=1.;
      }

    }

    if (thisTrackBed > maxTrackBed) {
      maxTrackBed = thisTrackBed;
      trialSDP   = tmpSDP;
    }

  }

  if( maxTrackBed <=0. ) {
    fGoodEvent = false;
    ostringstream warn;
    warn << " *************  no signal/pixels in track bed!!! Skip!  ****************";
    WARNING(warn);
    return utl::AxialVector(0.,0.,1.,eyeCS);
  }

  return trialSDP;

}

/* \brief call Minuit using MinuitFunction as function to be minimized
 */

void  FdSDPFinder::FindSDP(fevt::Eye& eye){

  fCurEye = &eye;

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

  //   use the eye coordinate system to do calculations
  CoordinateSystemPtr CS = deteye.GetEyeCoordinateSystem();


  const int npar = 2; // number of parameters

  TMinuit theMinuit(npar); // 2 parameter fit

  theMinuit.SetPrintLevel(-1); // minuit verbosity
  theMinuit.SetFCN(FdSDPFinder::MinuitFitFunc);

  double arglist[2];  // arguments to pass to the Minuit interpreter
  int ierflag =0;

  arglist[0]=1; // [up]
  theMinuit.mnexcm("SET ERR", arglist,1,ierflag);

  if (ierflag) ERROR ("Minuit SET ERR failed");

  double vstart[npar]; // inital point
  double step[npar];   // step size


  vstart[0] = fSDPFirstGuess.GetTheta(CS);
  vstart[1] = fSDPFirstGuess.GetPhi(CS);

  ostringstream info;
  info << "SDPn Initial guess (theta, phi): ("
       <<  vstart[0]/degree << ", " << vstart[1]/degree <<")";
  INFO(info);

  step[0]=0.01;
  step[1]=0.01;

  theMinuit.mnparm(0,"Theta",vstart[0], step[0],0,0,ierflag);
  theMinuit.mnparm(1,"Phi",  vstart[1], step[1],0,0,ierflag);

  arglist[0]=500;  // [maxcalls]
  arglist[1]=1;    // [tolerance]
  theMinuit.mnexcm("MIGRAD", arglist ,2,ierflag);

  if (ierflag) ERROR ("Minuit MIGRAD failed");

  theMinuit.mnexcm("HESSE", arglist, 2, ierflag);

  // Check the status of the covariance matrix
  // if bad use MINOS
  double fmin, fedm, errdef; // unused dummies
  int npari, nparx; //unused dummies
  int istat = 0;
  theMinuit.mnstat(fmin, fedm, errdef, npari, nparx, istat);
  if (istat < 3) { // no accurate covariance matrix
    theMinuit.mnexcm("MINOS", arglist, 2, ierflag);
    theMinuit.mnstat(fmin, fedm, errdef, npari, nparx, istat);
    if (istat < 3)
      WARNING("Neither HESSE nor MINOS found accurate covariance matrix");
  }

  double sdpTheta, sdpThetaErr;
  double sdpPhi,sdpPhiErr;

  theMinuit.GetParameter(0,sdpTheta,sdpThetaErr);
  theMinuit.GetParameter(1,sdpPhi,sdpPhiErr);

  info.str("");
  info << "Fitted SDPn (theta,phi): (" << sdpTheta/degree
       << "," << sdpPhi/degree << ")";
  INFO(info);

  fSDP = Vector(1,sdpTheta,sdpPhi,CS, Vector::kSpherical);
  fSDPThetaError = sdpThetaErr;
  fSDPPhiError = sdpPhiErr;

  //Retrieve the covariance matrix
  double covMatrix[npar][npar];
  theMinuit.mnemat( &covMatrix[0][0], npar );
  fSDPCorrThetaPhi = covMatrix[0][1]/sdpThetaErr/sdpPhiErr;

   fChi2 = theMinuit.fAmin;

}


void
FdSDPFinder::RemoveNoise(fevt::Eye& eye)
{
  EyeRecData& recdata = eye.GetRecData();

  // iterate over pixels with a pulse
  for (EyeRecData::PixelIterator pixeliter = recdata.SDPPixelsBegin();
       pixeliter != recdata.SDPPixelsEnd();) {

     fevt::Pixel& evtpixel = *pixeliter;
     const fdet::Pixel& detpixel =
       Detector::GetInstance().GetFDetector().GetPixel(evtpixel);
     const Vector& dir = detpixel.GetDirection();

     const double d = abs(kPi/2 - Angle(dir, fSDPFirstGuess));
     if (d > fMaxDistance )
       pixeliter =  recdata.RemoveSDPPixel(pixeliter);
     else ++pixeliter;
  }
}


void
FdSDPFinder::ReadmitPixel(fevt::Eye& eye)
{
  EyeRecData& recdata = eye.GetRecData();


  // iterate over pixels with a pulse
  for (EyeRecData::PixelIterator pixeliter = recdata.PulsedPixelsBegin();
       pixeliter != recdata.PulsedPixelsEnd(); ++pixeliter) {

    fevt::Pixel& evtpixel = *pixeliter;
    const fdet::Pixel& detpixel =
      Detector::GetInstance().GetFDetector().GetPixel(evtpixel);
    const unsigned int pixel_id = evtpixel.GetId();
    //-----------------------------
    // Checking whether this pixel was already
    // selected for SDP fitting.
    bool SDPfit_pixel = false;
    for (EyeRecData::ConstPixelIterator iter = recdata.SDPPixelsBegin();
         iter != recdata.SDPPixelsEnd(); ++iter)
      if (pixel_id == iter->GetId()) {
        SDPfit_pixel = true;
        break;
      }
    //-------------------------------------
    if (SDPfit_pixel)
      continue;
    Vector dir = detpixel.GetDirection();
    const double d = abs(kPi/2 - Angle(dir, fSDPFirstGuess));
    if (d < fMaxDistance && !IsIsolated(evtpixel, eye))
      recdata.AddSDPPixel(evtpixel);

  } // for pixel
}


bool
FdSDPFinder::IsIsolated(const fevt::Pixel& pixel, const fevt::Eye& eye)
{
  const EyeRecData& recdata = eye.GetRecData();
  const fdet::FDetector& detFD = det::Detector::GetInstance().GetFDetector();

  fevt::EyeRecData::ConstPixelIterator pixeliter;
  //----------------------------------------------------------
  //checking for wheter the pixel is isolated in time
  //or in space.

  //- For each triggered pixel the time distance to each of the rest
  //  of the triggered pixels is computed.
  //- If the time distance to all triggered pixels is greater
  //  than fIsolationLimit, the pixel is considered isolated
  //  in time and excluded.
  //- fIsolationLimit is estimated for every event. Depending
  //  on the event geometry fIsolationLimit can be as small as 1000ns.
  //  Far away showers need a bigger fIsolationLimit.
  //- If none of the neighbour pixels triggered, then the corresponding
  //  pixel is considered isolated in space. This pixel will not be used
  // to estimate the SDP first guess. However, it could be admitted
  // for the final SDP fitting.

  //---------------------------------------------------
  // Searching for the best  fIsolationLimi for this
  // particular event.
  //
  // The time differences among all pixels are distributed in the
  // array timediff_bin_population[400]. The array size is 400 bins
  // and the bin size is 500 ns.

  const int nbins = 400;
  int timediff_bin_population[nbins]={0};
  const double timediff_bin_size = 500.0;

  if ( fIsolationLimit == 0 ) {

    std::vector<double> candidatePixelTimes;
    for (pixeliter = recdata.PulsedPixelsBegin();
         pixeliter != recdata.PulsedPixelsEnd();
         ++pixeliter){

      const fdet::Channel& detChannel = detFD.GetChannel(*pixeliter);
      const int telid = pixeliter->GetTelescopeId();
      const unsigned int timeoffset = eye.GetTelescope(telid).GetTimeOffset();

      const double pixelTime = pixeliter->GetRecData().GetCentroid()*detChannel.GetFADCBinSize()
                               + timeoffset;

      candidatePixelTimes.push_back(pixelTime);
    }


    for (unsigned int i=0;i<candidatePixelTimes.size(); i++ ) {

      const double t1 = candidatePixelTimes[i];

      for (unsigned int j=0;j<candidatePixelTimes.size(); j++ ) {

        if ( i==j )
          continue;

        const double t2 = candidatePixelTimes[j];

        double timediff = std::abs(t1-t2);

        int timediff_bin = int(timediff / timediff_bin_size );
        if (timediff_bin >= nbins )
          continue;
        timediff_bin_population[timediff_bin]++ ;
      }
    }

    // Selecting the best fIsolationLimit

    int max_population = 0;
    int bin_with_maximum = 0;
    for ( int i = 0 ;  i < nbins  ; i++ ) {
      if ( timediff_bin_population[i] > max_population) {
        bin_with_maximum = i;
        max_population = timediff_bin_population[i];
      }
    }
    fIsolationLimit  = (bin_with_maximum+6) * timediff_bin_size ;

  }


  const fdet::Pixel& detpixel1     = detFD.GetPixel(pixel);
  const fdet::Channel& detChannel1 = detFD.GetChannel(pixel);
  const int telid1                 = pixel.GetTelescopeId();
  const unsigned int timeoffset1   = eye.GetTelescope(telid1).GetTimeOffset();

  const double pixelTime1 =  pixel.GetRecData().GetCentroid()*detChannel1.GetFADCBinSize()
                           + timeoffset1;

  const Vector& PixelDir1 = detpixel1.GetDirection();

  for (pixeliter = recdata.PulsedPixelsBegin();
       pixeliter != recdata.PulsedPixelsEnd();
       ++pixeliter){

    if (*pixeliter==pixel)
      continue;

    const fdet::Pixel& detpixel2     = detFD.GetPixel(*pixeliter);
    const fdet::Channel& detChannel2 = detFD.GetChannel(*pixeliter);
    const int telid2                 = pixeliter->GetTelescopeId();
    const unsigned int timeoffset2   = eye.GetTelescope(telid2).GetTimeOffset();

    const double pixelTime2 =  pixeliter->GetRecData().GetCentroid()*detChannel2.GetFADCBinSize()
                             + timeoffset2;

    const double timediff =
      std::abs(pixelTime1 - pixelTime2);

    const Vector& PixelDir2 = detpixel2.GetDirection();
    const double angdiff =  Angle (PixelDir1, PixelDir2);

    // Checking whether the trial pixel is
    // isolated in time or in space
    if (timediff < fIsolationLimit && angdiff < 5* degree)
      return false;

  } // for

  return true;
}



void FdSDPFinder::RemoveOutliers(fevt::Eye& eye) {


  if (!eye.HasRecData()) return;
  EyeRecData& recdata = eye.GetRecData();


  double TrackBed= 1.5 *degree;
  //use all pairs to generate candidate lines
  EyeRecData::PixelIterator iter;
  EyeRecData::PixelIterator iter2;

  vector<double> slope;
  vector<double> intcpt;

  for (iter=recdata.SDPPixelsBegin();
       iter!=recdata.SDPPixelsEnd(); ++iter){

    for (iter2= iter + 1;
         iter2!=recdata.SDPPixelsEnd();
         ++iter2){

      const fevt::Pixel& pix1 = *iter;
      const fevt::Pixel& pix2 = *iter2;

      double t1 = pix1.GetRecData().GetT_i().GetNanoSecond();
      double t2 = pix2.GetRecData().GetT_i().GetNanoSecond();
      double chi1 = pix1.GetRecData().GetChi_i();
      double chi2 = pix2.GetRecData().GetChi_i();

      double tmpSlope  =  (t1-t2) / (chi1-chi2);
      double tmpIntcpt  = t1 - tmpSlope*chi1;
      slope.push_back(tmpSlope);
      intcpt.push_back(tmpIntcpt);
    }

  }// for

  if (! slope.size()) return;

  // pick the line that accomodates more trackbed pixels

  unsigned int maxNPixels=0;
  unsigned int bestIndex=0;

  for (unsigned int i=0; i<slope.size(); i++){

    unsigned int tmpNPixels=0;
    for (iter= recdata.SDPPixelsBegin();
         iter!=recdata.SDPPixelsEnd();
         ++iter){

      const fevt::Pixel& pix = *iter;

      const double chi_i =  pix.GetRecData().GetChi_i();
      const double t_i   =  pix.GetRecData().GetT_i().GetNanoSecond();



      double dchi = chi_i - (t_i - intcpt[i])/slope[i];

      if (std::abs(dchi) < TrackBed * 2) tmpNPixels++;

    } // for pixel


    if (tmpNPixels> maxNPixels) {bestIndex=i;  maxNPixels= tmpNPixels;}


  }// for slope

  const double theSlope  =  slope[bestIndex];
  const double theIntcpt =  intcpt[bestIndex];

  if (  theSlope > 0 )  fSDPFirstGuess *= -1;

  // Remove outliers
  //  bool pixel_removed = false; // unused
  for (iter=recdata.SDPPixelsBegin();
       iter!=recdata.SDPPixelsEnd(); ){

    const fevt::Pixel& pix = *iter;
    const double chi_i =  pix.GetRecData().GetChi_i();
    const double t_i   =  pix.GetRecData().GetT_i().GetNanoSecond();

    double dchi = chi_i - (t_i - theIntcpt)/theSlope;

    if (std::abs(dchi) > TrackBed * 2  && recdata.GetNSDPPixels() > 3) {
      iter = recdata.RemoveSDPPixel(iter);
      //      pixel_removed = true; // unused

    } else  ++iter;



  }// for pixels
}



void
FdSDPFinder::MinuitFitFunc(int& /*npar*/, double* /*gin*/, double& f,
                           double* par, int /*iflag*/)
{
  // fCurEye is an escamotage : we need something static for TMinuit,
  // which needs a static function ! Ahhh ! (sa,23122003)
  const EyeRecData& recdata = fCurEye->GetRecData();

  EyeRecData::ConstPixelIterator pixeliter;

  double chisq =0;

  const fdet::Eye& deteye =
    Detector::GetInstance().GetFDetector().GetEye(*fCurEye);

  // use the eye coordinate system (must be the same we used before)
  CoordinateSystemPtr CS= deteye.GetEyeCoordinateSystem();

  double& sdpTheta = par[0];
  double& sdpPhi   = par[1];

  // the trial SDP
  Vector sdp (1,sdpTheta ,sdpPhi ,CS, Vector::kSpherical);

  // iterate over pixels with a pulse
  unsigned int ndof =0;
  double total_charge=0;

  for (pixeliter = recdata.SDPPixelsBegin();
       pixeliter != recdata.SDPPixelsEnd();
       ++pixeliter){

     const fevt::Pixel& evtpixel = *pixeliter;
     const fdet::Pixel& detpixel =
       Detector::GetInstance().GetFDetector().GetPixel(evtpixel);

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

     // distance from trial plane
     double tmp = kPi/2. - Angle( dir , sdp) ;

     // weight to assign
     double charge = evtpixel.GetRecData().GetTotalCharge();

     // 0.35 degree from Jose's studies on laser shots
     chisq += tmp*tmp  * charge /(0.35 * degree * 0.35 * degree);
     ndof++;
     total_charge += charge;
  }// for pixeliter

  // normalize chisq
  chisq /= total_charge;
  chisq *= ndof;
  //
  f= chisq;

  fNDof = ndof -2 ;
}

void
FdSDPFinder::EstimateChi_i(fevt::Eye& eye)
{
  fevt::EyeRecData& recdata = eye.GetRecData();
  const fdet::Eye& deteye =
    det::Detector::GetInstance().GetFDetector().GetEye(eye);

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

  // iterate over pulsed pixels
  for (fevt::EyeRecData::PixelIterator pixeliter = recdata.PulsedPixelsBegin();
       pixeliter != recdata.PulsedPixelsEnd(); ++pixeliter) {

    fevt::Pixel& evtpixel = *pixeliter;
    const fdet::Pixel& detpixel =
      det::Detector::GetInstance().GetFDetector().GetPixel(evtpixel);
    const fdet::Channel& detChannel =
      det::Detector::GetInstance().GetFDetector().GetChannel(evtpixel);
    const double timebinsize = detChannel.GetFADCBinSize();


    // time correction for different eyes
    const int telid = evtpixel.GetTelescopeId();
    const unsigned int timeoffset = eye.GetTelescope(telid).GetTimeOffset();
    const double t_i =
      timebinsize * evtpixel.GetRecData().GetCentroid() + timeoffset;

    //----------------------------------------------
    //estimating Chi_i

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

    Vector vertical(0,0,1, cs, Vector::kCartesian);
    Vector horizontalWithinSDP = cross(fSDPFirstGuess, vertical);
    horizontalWithinSDP /= horizontalWithinSDP.GetR(cs);


    // subtract component parallel to sdp
    Vector chi_i_vector = pixeldir - (fSDPFirstGuess * pixeldir) *  fSDPFirstGuess;
    chi_i_vector.TransformTo(cs);
    //     chi_i_vector.Normalize();
    chi_i_vector /= chi_i_vector.GetR(cs);
    const double chi_i = Angle(chi_i_vector, horizontalWithinSDP);

    evtpixel.GetRecData().SetChi_i(chi_i);
    evtpixel.GetRecData().SetT_i(t_i, 100);

  }
}


// Configure (x)emacs for this file ...
// Local Variables:
// mode:c++
// compile-command: "make -C .. FdSDPFinderOG/FdSDPFinder.o  -k"
// End: