// ESAF : Euso Simulation and Analysis Framework
// $Id: HmaxForProtonModule.cc,v 1.8 2005/03/25 18:32:43 naumov Exp $
// 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 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
tmpHMaxProtonPointer = this;
}
//_____________________________________________________________________________
HmaxForProtonModule::~HmaxForProtonModule() {
// dtor
}
//_____________________________________________________________________________
Bool_t HmaxForProtonModule::Init() {
Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;
return kTRUE;
}
//_____________________________________________________________________________
Bool_t HmaxForProtonModule::PreProcess() {
return kTRUE;
}
//_____________________________________________________________________________
Bool_t HmaxForProtonModule::Process(RecoEvent *ev) {
const RecoModuleData *td = ev->GetModuleData("TrackDirection");
if ( td == NULL ) {
Msg(EsafMsg::Panic) << "TrackDirectionModule is not found." << MsgDispatch;
throw runtime_error( "HmaxForProtonModule failed.");
}
fTheta = td->GetDouble("Theta");
fPhi = td->GetDouble("Phi");
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() {
}