source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/reconstruction/modules/shower/energy/src/HmaxByShapeMethod.cc @ 117

// ESAF : Euso Simulation and Analysis Framework
// $Id: HmaxByShapeMethod.cc 2946 2011-06-25 13:37:26Z mabl $
// Dmitry V. Naumov created Jun,  6 2004

#include "HmaxByShapeMethod.hh"
#include "Atmosphere.hh"
#include "utils.hh"
#include "RecoEvent.hh"
#include "EDetector.hh"
#include "EConst.hh"
#include "EnergyTypes.hh"
#include "EnergyViewer.hh"

#include <TH1F.h>
#include <TF1.h>
// #include <TMinuit.h>
#include <TVirtualFitter.h>

#include <iostream>

using namespace sou;
using namespace TMath;
using namespace EConst;
using namespace std;

ClassImp(HmaxByShapeMethod)
// fFluorescence function
//_____________________________________________________________________________
Double_t HSM_fFluorescence(Double_t *x, Double_t *par) {
    return par[0]*exp(-0.5*pow(x[0]-par[1],2)/pow(par[2],2));
}


// Sum of background and peak functions
//_____________________________________________________________________________
Double_t HSM_fitFunction(Double_t *x, Double_t *par) {
    return HSM_fFluorescence(x,&par[0]);
}

//_____________________________________________________________________________
float HSM_DensityToAltitude(float h) {
//     HmaxByShapeMethod* hm = (HmaxByShapeMethod*)gMinuit->GetObjectFit();
    HmaxByShapeMethod* hm = dynamic_cast<HmaxByShapeMethod*>( (TVirtualFitter::GetFitter())->GetUserFunc() );
    return Atmosphere::Get()->Air_Density(h) - hm->GetDensity();
}

//_____________________________________________________________________________
HmaxByShapeMethod::HmaxByShapeMethod(EnergyModule* em): EsafMsgSource() {
    //
    // ctor
    //
    if (!em) Msg(EsafMsg::Panic)<< "Pointer to EnergyModule is NULL. I do not want to work like this" <<MsgDispatch;

    fEnergyModule  = em;
    fEnergy        = 1.e14;
    fFitFcn        = NULL;
    fSigma         = -HUGE;
    fAltitude      = -HUGE;
    fDensity       = -HUGE;
}

//_____________________________________________________________________________
HmaxByShapeMethod::HmaxByShapeMethod(): EsafMsgSource() {
    //
    // ctor
    //
    fEnergyModule  = NULL;
    fEnergy        = 1.e14;
    fFitFcn        = NULL;
    fSigma         = -HUGE;
    fAltitude      = -HUGE;
    fDensity       = -HUGE;
}

//_____________________________________________________________________________
HmaxByShapeMethod::~HmaxByShapeMethod() {
    //
    // dtor
    //
    Clear();
}

//_____________________________________________________________________________
void HmaxByShapeMethod::Clear()  {
    //
    // clean 
    //
    SafeDelete(fFitFcn);
}

//_____________________________________________________________________________
Bool_t HmaxByShapeMethod::Fit(TH1F* data) {
    //
    // fit entry histogramm and estimate altitude of the shower maximum. The entry histo is assumed to be unfolded of inefficients  
    //
    Clear();
    RecoEvent *ev = fEnergyModule->GetRecoEvent();
    Int_t fFirstHit = ev->GetStartGtu();
    Int_t fLastHit  = ev->GetEndGtu();
    Float_t gtulength = fEnergyModule->GetRecoEvent()->GetHeader().GetGtuLength();
    FitHisto(data,fFirstHit,fLastHit,gtulength);
    Int_t debug = fEnergyModule->GetDebugLevel();
#ifdef USE_VISUALIZATION
    if (debug>=1) {
        EnergyViewer *fViewer = fEnergyModule->GetEnergyViewer();
        fViewer->DisplayAltitude(data,fFitFcn);
    }
#endif
    return kTRUE;
}

//_____________________________________________________________________________
void HmaxByShapeMethod::FitHisto(TH1F* data, Int_t first, Int_t last, Float_t gtulength) {
    //
    // fit entry histogramm and estimate altitude of the shower maximum. The entry histo is assumed to be unfolded of inefficients  
    //
    //gMinuit->SetObjectFit(this);

    Double_t par[3]={1,1,1};
    fFitFcn = new TF1("fFitFcn",HSM_fitFunction,first,last,3);
    fFitFcn->SetNpx(500);
    fFitFcn->SetLineWidth(2);
    fFitFcn->SetLineColor(kMagenta);

    fFitFcn->SetParLimits(0,0,data->GetMaximum());
    fFitFcn->SetParLimits(1,0,last);
    fFitFcn->SetParLimits(2,0,0.5*last); 
    fFitFcn->SetParameter(1,fFluorMaxTime);

    data->Fit("fFitFcn","RQ");
    fChi2 = data->GetFunction("fFitFcn")->GetChisquare()/data->GetFunction("fFitFcn")->GetNDF();
    data->GetFunction("fFitFcn")->GetParameters(&par[0]);
    fFitFcn->SetParameters(par);
    fFluorMaxTime = par[1];
    TVirtualFitter *min = TVirtualFitter::Fitter(0,4);
    min->SetUserFunc(this);

    // Compute altitude of the shower maximum
    Float_t x0 = 37.15*g/cm2;
    Float_t t0 = 1.7+0.76*Log(fEnergy/81.);                   // energy in MeV !!!
    fSigma = par[2]*gtulength;
    fDensity = Sqrt(2*t0)*x0*(1.+fShowerDir*fShowerMaxDir)/Clight()/fSigma;
    DensityToAltitude();
    Msg(EsafMsg::Debug) << "Reconstructed altitude of the shower maximum "<< fAltitude/km <<MsgDispatch;
    Msg(EsafMsg::Debug) << "Reconstructed density at the shower maximum " << fDensity/g*cm3 << " sigma " << fSigma/microsecond << MsgDispatch;
}

//_____________________________________________________________________________
void HmaxByShapeMethod::DensityToAltitude() {
    //
    // Finds the altitude which corresponds to the given density
    Float_t h1 = 0*km, h2 = 100*km; 
    Float_t (*f)(float)= &HSM_DensityToAltitude;
    zbrac(f,&h1,&h2);
    fAltitude = rtbis(f,h1,h2,0.01*km);
}


//_____________________________________________________________________________
void HmaxByShapeMethod::Set(Float_t energy,Float_t timemax, TVector3 showerdir, TVector3 showermaxdir){
    //
    // Set up the module
    // 
    SetEnergy(energy);
    SetFluoMaxTime(timemax);
    SetShowerDir(showerdir);
    SetShowerMaxDir(showermaxdir);
}