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

// ESAF : Euso Simulation and Analysis Framework
// $Id: HmaxForProtonModule.cc 2633 2006-03-28 09:13:43Z thea $
// Naumov Dmitry created Jul,  4 2004

#include "HmaxForProtonModule.hh"
#include "RecoRootEvent.hh"
#include "RecoEvent.hh"
#include "RecoPixelData.hh"
#include "utils.hh"
#include "Atmosphere.hh"

using namespace sou;
using namespace TMath;

ClassImp(HmaxForProtonModule)
HmaxForProtonModule* tmpHMaxProtonPointer = NULL;

//_____________________________________________________________________________
float HmaxForProton_depthToAltitude(float h) {
    // Root of this function is the Altitude along the direction Theta which correspond to GIL estimation of the
    // Shower Maximum
    Double_t Theta = tmpHMaxProtonPointer->GetTheta();
    return Atmosphere::Get()->Depth(h*km,Theta) - tmpHMaxProtonPointer->XMaxEstimate();
}
    
//_____________________________________________________________________________
Double_t HmaxForProtonModule::XMaxEstimate() {
    Double_t t0 = 37*g/cm2;
    return t0*(1.7+0.76*TMath::Log(fEnergy/0.81e8));        
}

//_____________________________________________________________________________
HmaxForProtonModule::HmaxForProtonModule() : RecoModule("HmaxForProton") {
    // ctor
    fUseTrueAngles = kFALSE;
    tmpHMaxProtonPointer = this;
}

//_____________________________________________________________________________
HmaxForProtonModule::~HmaxForProtonModule() {
    // dtor
}

//_____________________________________________________________________________
Bool_t HmaxForProtonModule::Init() {
    Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;
    fUseTrueAngles     = (Int_t)Conf()->GetNum("HmaxForProtonModule.fUseTrueAngles");
    return kTRUE;
}

//_____________________________________________________________________________
Bool_t HmaxForProtonModule::PreProcess() {
    return kTRUE;
}

//_____________________________________________________________________________
Bool_t HmaxForProtonModule::Process(RecoEvent *ev) {
    const RecoModuleData *td = ev->GetModuleData("TrackDirection");
    if ( td == NULL ) td = ev->GetModuleData("TrackDirection2"); // try to get TrackDirection2
    if (fUseTrueAngles) {
        fTheta = ev->GetHeader().GetTrueTheta();
        fPhi   = ev->GetHeader().GetTruePhi();
    }
    else {
        if ( td == NULL ) Msg(EsafMsg::Panic) << "None of TrackDirectionModule and TrackDirectionModule2 is not found." << MsgDispatch;
        fTheta = td->GetDouble("Theta");
        fPhi   = td->GetDouble("Phi");
        if(IsNaN(fTheta) || IsNaN(fPhi)) return kFALSE;
    }

    const RecoModuleData *en = ev->GetModuleData("Energy");
    fEnergy = 1.e20; // default guess
    if ( en != NULL ) 
        fEnergy = en->GetDouble("Energy");
    Msg(EsafMsg::Debug) << "Processing Theta " << fTheta*RadToDeg() << " deg, Phi " << fPhi*RadToDeg() << " deg" << MsgDispatch;    
    
    if (fTheta*RadToDeg() < 0 || fTheta*RadToDeg() > 92) return kFALSE;
    
    Float_t t1(0),t2(100);
    Float_t (*f)(float)= &HmaxForProton_depthToAltitude;
    zbrac(f,&t1,&t2);
    fAltitude =  rtbis(f,t1,t2,0.01)*km;
    MyData()->Add("Altitude",fAltitude);
    Msg(EsafMsg::Debug) << "Hmax = " << fAltitude/km << " km " << MsgDispatch;    
    
    return kTRUE;
}

//_____________________________________________________________________________
Bool_t HmaxForProtonModule::PostProcess() {
    return kTRUE;
}

//_____________________________________________________________________________
Bool_t HmaxForProtonModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
    fRecoRootEvent->GetRecoHmaxForProton().SetHmax(fAltitude);
    fRecoRootEvent->GetRecoHmaxForProton().SetQuality(1);
    Double_t HmaxError  = fRecoRootEvent->GetTruth().GetHmax() - fAltitude;
    fRecoRootEvent->GetRecoHmaxForProton().SetErrorHmax(HmaxError);
    return kTRUE;
}
//_____________________________________________________________________________
Bool_t HmaxForProtonModule::Done() {
    Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
    return kTRUE;
}

//______________________________________________________________________________
void HmaxForProtonModule::UserMemoryClean()  {
    
}

//______________________________________________________________________________
void HmaxForProtonModule::TmpMemoryClean()  {
    
}