FdEnergyFinder.cc

Go to the documentation of this file.

#include "FdEnergyFinder.h"

#include <evt/Event.h>
#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/EyeRecData.h>
#include <evt/ShowerFRecData.h>
#include <evt/GaisserHillas4Parameter.h>
#include <fwk/CentralConfig.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/AugerUnits.h>
#include <utl/PhysicalConstants.h>
#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
//BRD added 28/2/05
#include <utl/PhysicalFunctions.h>

#include <iostream>
#include <TMinuit.h>



using namespace FdEnergyFinderOG;
using namespace utl;
using namespace evt;
using namespace std;

const utl::TabulatedFunctionErrors * FdEnergyFinder::fLongProfile =0;
double FdEnergyFinder::fChisqGH =0;
double FdEnergyFinder::fNdof=0;

fwk::VModule::ResultFlag FdEnergyFinder::Init(){

  fwk::CentralConfig* cc = fwk::CentralConfig::GetInstance();

  Branch topB =
     cc->GetTopBranch("FdEnergyFinder");

  topB.GetChild("fixX0").GetData(fFixX0);
  topB.GetChild("composition").GetData(fComposition);
  topB.GetChild("unseenEnergyModel").GetData(fUnseenEnergyModel);

  // SYB: Checks are now performed by XML schema.
//  if ( fUnseenEnergyModel > 2 || fUnseenEnergyModel < 1  ) {
//    ERROR("Error - unknown unseen energy model (1=QGSJET, 2=SIBYLL");
//    return eFailure;
//  }
//  //DV for unsigned this is always false
//  if (/*fComposition < 0 ||*/ fComposition > 2) {
//    ERROR("Error - unknown composition (0 = proton-iron mixed; 1 = protons; 2 = iron)");
//    return eFailure;
//  }

  fNmaxGH = 0;
  fXmaxGH = 0;
  fX0GH = 0;
  fEemGH = 0;
  fNmaxGH_error = 0;
  fXmaxGH_error = 0;
  fX0GH_error = 0;
  fChisqGH = 0;
  fNdof =0;
  fLongProfile =0;

  return eSuccess;
}

fwk::VModule::ResultFlag FdEnergyFinder::Run(evt::Event& event){

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

  fevt::FEvent::EyeIterator eyeiter;

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


    fNmaxGH = 0;
    fXmaxGH = 0;
    fX0GH = 0;
    fEemGH = 0;
    fNmaxGH_error = 0;
    fXmaxGH_error = 0;
    fX0GH_error = 0;
    fChisqGH = 0;
    fNdof =0;
    fLongProfile =0;
    fMinuitFailed=false;

    if (! eyeiter->HasRecData() ) continue;
    if (! eyeiter->GetRecData().HasFRecShower() ) continue;
    ShowerFRecData& frecshower = eyeiter->GetRecData().GetFRecShower();

    // Getting energy cutoff
    fEnergyCutoff = frecshower.GetEnergyCutoff();
    cout << "Energy cutoff = " << fEnergyCutoff/MeV << " [MeV]" << endl;

    if (!frecshower.HasLongitudinalProfile()) continue;

    const TabulatedFunctionErrors& long_profile =
      frecshower.GetLongitudinalProfile();

    fLongProfile = &long_profile;


    if (!frecshower.HasGHParameters()) {
        evt::GaisserHillas4Parameter gh_tmp;
        frecshower.MakeGHParameters (gh_tmp);
    }

    fGHParameters = dynamic_cast <evt::GaisserHillas4Parameter*>
        (&(frecshower.GetGHParameters()));

    if (fGHParameters==0) {
        WARNING ("ShowerFRecData of this eye already contains GH parameters.");
        continue;
    }

    fGHParameters->SetShapeParameter(gh::eLambda, kLambdaGH, 0.);

    this->FitProfile();

    fGHParameters->SetXMax(fXmaxGH,fXmaxGH_error );
    fGHParameters->SetShapeParameter(gh::eX0, fX0GH, fX0GH_error);
    fGHParameters->SetChiSquare( fChisqGH, (unsigned int)fNdof);
    fGHParameters->SetNMax(fNmaxGH,fNmaxGH_error );

    if (fMinuitFailed){
      fGHParameters->SetXMax (-999,0);
      fGHParameters->SetShapeParameter(gh::eX0, -999, 0);
      fGHParameters->SetChiSquare( -999, 0);
      fGHParameters->SetNMax(-999, 0 );
    }

    this->FindEmEnergy();

    this->FindEmEnergyError();

    this->FindTotalEnergy();

    frecshower.SetEmEnergy(fEemGH,fEemGH_error);

    frecshower.SetTotalEnergy(fEtotalGH,fEemGH_error );


    cout << "Em energy [GeV]" << fEemGH/GeV << " +- "
         << fEemGH_error/GeV << endl;
    cout << "Total Energy [GeV] " << fEtotalGH/GeV << endl;

  } // loop on eyes

  return eSuccess;
}

void FdEnergyFinder::FitProfile() {

  const int npar = 3; // number of parameters

  TMinuit theMinuit(npar); // 3 parameter fit

  theMinuit.SetPrintLevel(-1); // minuit verbosity
  theMinuit.SetFCN(FdEnergyFinder::MinuitFitFunction);

  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

  double log_n_max = 20.7;
  double x_max = 700 * (g/cm2);
  double x0 = 0.0 * (g/cm2);


  // something better has to be found for starting points
  vstart[0] = log_n_max;
  vstart[1] = x_max;
  vstart[2] = x0;

  step[0] = 0.05;
  // step[1]=0.5;
  // step[2]=0.1;
  step[1] = 10 * (g/cm2);
  step[2] = 10 * (g/cm2);

  if (fFixX0)  step[2] = 0.0;

  cout << " x0 step " << step[2] << endl;

  theMinuit.mnparm(0,"log_n_max",vstart[0], step[0],0,0,ierflag);
  theMinuit.mnparm(1,"x_max",  vstart[1], step[1],0,0,ierflag);
  theMinuit.mnparm(2,"x0",  vstart[2], step[2],0,0,ierflag);

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

  if (ierflag) {
    ERROR ("Minuit MIGRAD failed");
    fMinuitFailed=true;
  }

  double p[6];
  theMinuit.GetParameter(0, p[0], p[3]);
  theMinuit.GetParameter(1, p[1], p[4]);
  theMinuit.GetParameter(2, p[2], p[5]);

  fNmaxGH = exp(p[0]);
  fXmaxGH = p[1];
  fX0GH = p[2];

  fNmaxGH_error = p[3]* fNmaxGH;
  fXmaxGH_error = p[4];
  fX0GH_error =   p[5];

  cout << endl;
  cout << endl << "Gaisser-Hillas Fit result:" << endl;
  cout << "n_max " << fNmaxGH << " +/- " << fNmaxGH_error << "\n"
       <<  " X_max [g/cm2] " << fXmaxGH / (g/cm2)
       << " +/- " << fXmaxGH_error  / (g/cm2) << "\n"
       << " X0 " << fX0GH / (g/cm2)
    << " +/- " << fX0GH_error / (g/cm2) << endl;

  cout << "Energy Estimate from Nmax[Gev] " << fNmaxGH / 1.3 << endl;

}


fwk::VModule::ResultFlag
FdEnergyFinder::Finish()
{
  return eSuccess;
}


void
FdEnergyFinder::MinuitFitFunction(int& /*npar*/, double* /*gin*/, double& f,
                                  double* par, int /*iflag*/)
{
  // Eventually need maximum likelihood.  But for now,
  // calculate chisq bin by bin, and ignore empty bins.  Statistics
  // may not be quite correct for small pe numbers, but near enough.

  double chisq = 0;

  double Nmax= exp(par[0]);;
  double Xmax= par[1]; ;
  double X0  = par[2];;

  int    ndof_chisq =0;
  double size = 0.;

  const unsigned int npoints =  fLongProfile->GetNPoints();
  for (unsigned int i=4; i<npoints; i++){

    double X       =  fLongProfile->GetX(i);
    double n_e     =  fLongProfile->GetY(i);
    double n_e_err =  fLongProfile->GetYErr(i);

    if (X>X0 && Xmax>0) {
      size = Nmax*pow((X-X0)/(Xmax-X0),(Xmax-X0)/kLambdaGH)*exp((Xmax-X)/kLambdaGH);
      chisq += pow((n_e - size)/ n_e_err,2);
      ndof_chisq++;
    } // if X

  }// for points

  if (ndof_chisq < 1) chisq = 9999;

  f = chisq;

  fChisqGH = f;
  fNdof = (ndof_chisq - 3);
}

void FdEnergyFinder::FindEmEnergy() {

    if (fMinuitFailed) {
        fEemGH = 0;
        return;
    }

    /*
  double alpha = (fXmaxGH - fX0GH)/kLambdaGH;

  if(fNmaxGH <= 0. || fXmaxGH <= 0. || alpha <= 0.) {
    fEemGH = 0.;
    return;
  }

  if(alpha >= 40.)
    alpha = 40.; // avoids overflows in gh_gamma

  double gamma = fGHParameters->GammaIntegral(alpha);
    */

  //BRD New Eem calculation using mean dE/dX
  double size = 0;
  double sum_size=0;
  double sum_dEdX=0;
  for (unsigned int i=0; i<50; i++){
    double X = i*40.*(g/cm2);
    if (X>fX0GH && fXmaxGH>0) {
      size = GaisserHillas(X, fX0GH, fXmaxGH, fNmaxGH, kLambdaGH);
      sum_size += size;
      //BRD NOTE!
      //double fEnergyCutoff = 1.0*MeV;
      sum_dEdX += size*EnergyDeposit(X, fXmaxGH, fEnergyCutoff);
    }
  }
  double mean_dEdX=0.;
  if (sum_size)
    mean_dEdX = sum_dEdX/sum_size;
  cout << "Mean dEdX " << mean_dEdX/(MeV/(g/cm2)) << endl;
  //fEemGH = gamma * fNmaxGH * kLambdaGH * mean_dEdX;
  //

  fEemGH = fGHParameters->GetIntegral() * mean_dEdX;

}


void FdEnergyFinder::FindEmEnergyError() {

    if (fMinuitFailed) {
        fEemGH_error = 9.9e99;
        return;
    }

    fEemGH_error = fEemGH * fGHParameters->GetIntegralError();

    /*
  double alpha = (fXmaxGH - fX0GH)/kLambdaGH;

  if (fXmaxGH_error == 0. && fX0GH_error == 0. ){
    fEemGH_error = 9.9e99;
    return;
  }

  double alpha_error = sqrt(fXmaxGH_error*fXmaxGH_error
                            + fX0GH_error*fX0GH_error) / kLambdaGH;


  double alpha1 = alpha + alpha_error;
  double alpha2 = alpha - alpha_error;

  if(fNmaxGH <= 0. || fXmaxGH <= 0. || alpha <= 0. || alpha2 <= 0. )
    {
      fEemGH_error = 9.9e99;
      return;
    }

  if(alpha >= 40.) alpha = 40.; //avoids overflows in gh_gamma

  if(alpha1 >= 40.) alpha1 = 40.; //avoids overflows in gh_gamma

  double gamma = fGHParameters->GammaIntegral(alpha);

  double gamma1 = fGHParameters->GammaIntegral(alpha1);

  double gamma2 = fGHParameters->GammaIntegral(alpha2);


  double f1 = (gamma1-gamma2)/gamma/2;
  double f2 = fNmaxGH_error/fNmaxGH;
  fEemGH_error = sqrt(f1*f1 + f2*f2);
  fEemGH_error *= fEemGH;
    */

}


void
FdEnergyFinder::FindTotalEnergy()
{
  using namespace utl::InvisibleEnergy;
  const EInteractionModel model =
    (fUnseenEnergyModel == 1) ? eQGSJET01 : eSIBYLL21;

  ECompositionModel compo;
  if (fComposition == 1)
    compo = eProton;
  else if (fComposition == 2)
    compo = eIron;
  else
    compo = eMixedComposition;

  double dummycostheta=0.7071067; //this module does not use the data based invisible energy, so theta is not required
  fEtotalGH = fEemGH * Factor(fEemGH, model, compo,dummycostheta);
}


// /*!

//   RETURN total energy BY USING Eem/ETotal BASED ON CORSIKA SHOWER MONTE CARLO
//   (QGSJET MODEL).   CHS 6/18/98
//   NEW PARAMERIZATION HAVE BEEN DONE. THE THRESHOLD ENERGY OF ELECTRON
//   AND PHOTON BOTH IS 100 KeV(FROM 3 MeV ).   CHS 11/09/98
//   Ref: "Primary Energy Reconstruction through Fluorescence Light Detection"
//   tech note by Chihwa Song.
//   1 Dec 1998:  Mod by BRD
//   Reduce "missing energy" by 5% for all energies/primaries.  Chihwa
//   has shown that about 10% of primary energy ends up as sub-Mev gamma
//   -rays.  These will not participate in cascading, but will produce
//   energetic e's thru Compton scattering and photoelectric effect.
//   We have not evaluated how much of this ~10% will make it into
//   visible energy - guess 5%, then we can only be a few % wrong.

// @note takes energy in GeV, which is not auger-convention correct
// */
// void FdEnergyFinder::FindTotalEnergy(double Eem){

//   const int nx = 6;
//   //std::vector<double> xa;
//   //std::vector<double> ya;
//   double xa[nx+1];
//   double ya[nx+1];
//   double minlog10energy = 0.;
//   double maxlog10energy = 0.;

//   // Paramerization
// // Mean between proton and iron
// const double mixed_parameters[3][7]={{ 0,0,0,0,0,0,0},
//                           { 0, 7.342, 7.877, 8.364, 8.895, 9.379,  9.907},
//                           { 0, 0.733, 0.754, 0.770, 0.786, 0.797,  0.807}};
// // For proton
// const double proton_parameters[3][7]={{0,0,0,0,0,0,0},
//                            {0, 7.366,  7.895,  8.379,  8.908,  9.389,  9.915},
//                            {0, 0.775,  0.786,  0.797,  0.810,  0.816,  0.823}};
// // For iron
// const double iron_parameters[3][7]={{0,0,0,0,0,0,0},
//                          {0, 7.317,  7.858,  8.348,  8.881,  9.368,  9.898},
//                          {0, 0.691,  0.721,  0.742,  0.761,  0.777,  0.790}};


//   // Get energy conversion factor for given primary particle type
//   if (fComposition == 0)
//     {
//       for (int i=1; i<=nx; i++)
//      {
//        xa[i] = mixed_parameters[1][i];
//        ya[i] = mixed_parameters[2][i];
//      }
//       minlog10energy = 7.342;
//       maxlog10energy = 9.907;
//     }
//   else if (fComposition == 1)
//     {
//       for (int i=1; i<=nx; i++)
//      {
//        xa[i] = proton_parameters[1][i];
//        ya[i] = proton_parameters[2][i];
//      }
//       minlog10energy = 7.366;
//       maxlog10energy = 9.915;
//     }
//   else if (fComposition == 2)
//     {
//       for (int i=1; i<=nx; i++)
//      {
//        xa[i] = iron_parameters[1][i];
//        ya[i] = iron_parameters[2][i];
//      }
//       minlog10energy = 7.317;
//       maxlog10energy = 9.898;
//     }
//   else
//     {
//       cerr << "Error by FdEnergyFinder: primary type must be 0, 1 or 2"
//         << endl;
//       return;
//     }


//   //  If Eem is out of range, use the value at end point

//   double log10E = log10(Eem);

//   if (log10E < minlog10energy)
//     {
//       log10E = minlog10energy;
//       WARNING("Eem is out of range");
//     }
//   else if (log10E > maxlog10energy)
//     {
//       log10E = maxlog10energy;
//       WARNING("Eem is out of range");
//     }

//   double Eem_to_Etot = PolinomialInterpolation(xa, ya, nx, log10E);

//   //  Add extra 5% to visible energy (see note at top.  BRD 1 Dec 98)

//   Eem_to_Etot *= 1.05;

//   fEtotalGH = Eem/Eem_to_Etot * GeV;


// }


// Polynomial function interpolation and extrapolation from
// Numerical Recipes
// double FdEnergyFinder::PolinomialInterpolation(double xa[], double ya[],
//                                                  int n, double exp){

//   const int nmax = 10;
//   int  ns = 0;
//   double diff_i = 0.;
//   double c[nmax];
//   double d[nmax];
//   double Pol_Int = 0.;
//   double dPol_Int = 0.;
//   double den, ho, hp, w;

//   ns = 1;
//   double diff_1 = abs(exp-xa[1]);
//   for (int i=1; i<=n; i++)
//     {
//       diff_i = abs(exp-xa[i]);
//       if (diff_i < diff_1)
//      {
//        ns = i;
//        diff_1 = diff_i;
//      }
//       c[i] = ya[i];
//       d[i] = ya[i];
//     }

//   Pol_Int  = ya[ns];
//   ns = ns - 1;

//   for (int m=1; m<=n-1; m++)
//     {
//       for (int i=1; i<=n-m; i++)
//      {
//        ho  = xa[i] - exp;
//        hp  = xa[i+m] - exp;
//        w   = c[i+1] - d[i];
//        den = ho - hp;
//        if (den == 0)
//          {
//            cerr << "Fatal error by PolinomialInterpolation: denominator = 0" << endl;
//            exit(1);
//          }
//        den  = w/den;
//        d[i] = hp*den;
//        c[i] = ho*den;
//      }
//       if (2*ns < n-m)
//      dPol_Int = c[ns+1];
//       else
//      {
//        dPol_Int = d[ns];
//        ns = ns - 1;
//      }

//       Pol_Int = Pol_Int + dPol_Int;
//     }

//     return Pol_Int;

// }

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