// ESAF : Euso Simulation and Analysis Framework
// $Id: XmaxModule.cc,v 1.6 2004/11/25 18:49:32 naumov Exp $
// Naumov Dmitry created Jul, 4 2004
#include "XmaxModule.hh"
#include "RecoRootEvent.hh"
#include "RecoEvent.hh"
#include "EarthVector.hh"
#include "Atmosphere.hh"
ClassImp(XmaxModule)
//_____________________________________________________________________________
XmaxModule::XmaxModule() : RecoModule("Xmax"){
// ctor
}
//_____________________________________________________________________________
XmaxModule::~XmaxModule() {
// dtor
}
//_____________________________________________________________________________
Bool_t XmaxModule::Init() {
Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;
return kTRUE;
}
//_____________________________________________________________________________
Bool_t XmaxModule::PreProcess() {
fAltitude = 100;
fAltitudeProton = 100;
return kTRUE;
}
//_____________________________________________________________________________
Bool_t XmaxModule::Process(RecoEvent *ev) {
const RecoModuleData *td = ev->GetModuleData("TrackDirection");
if ( td==NULL ) {
Msg(EsafMsg::Panic) << "TrackDirectionModule is not found." << MsgDispatch;
throw runtime_error( "XmaxModule failed.");
}
fTheta = td->GetDouble("Theta");
fPhi = td->GetDouble("Phi");
const RecoModuleData *gcm = ev->GetModuleData("GTUClustering");
if ( gcm==NULL ) {
Msg(EsafMsg::Panic) << "GTUClusteringModule is not found." << MsgDispatch;
throw runtime_error("TrackDirectionModule failed.");
}
else {
if ( gcm->GetObj("CluPixels") == NULL ) {
Msg(EsafMsg::Panic) << "No data in event" << MsgDispatch;
return kTRUE;
}
fPointsId = *(vector<Int_t>*)gcm->GetObj("CluPixels");
if (fPointsId.size() == 0) {
fNumPoints = 0;
Msg(EsafMsg::Panic) << "No pixels selected: terminated." << MsgDispatch;
return kTRUE;
} else {
fNumPoints = fPointsId.size();
Msg(EsafMsg::Debug) << "Processing " << fNumPoints << " points" << MsgDispatch;
}
}
const RecoModuleData *hsm = ev->GetModuleData("HmaxByShapeMethod");
const RecoModuleData *hpr = ev->GetModuleData("HmaxForProton");
if ( hsm == NULL && hpr == NULL) {
Msg(EsafMsg::Panic) << "Not found neither HmaxForProtonModule nor HmaxByShapeMethodModule." << MsgDispatch;
throw runtime_error("XmaxModule failed.");
}
if ( hsm != NULL)
fAltitude = hsm->GetDouble("Altitude");
if ( hpr != NULL)
fAltitudeProton = hpr->GetDouble("Altitude");
fXmax = Atmosphere::Get()->Depth(fAltitude*km,fTheta)*cm2/g;
fXmaxP = Atmosphere::Get()->Depth(fAltitudeProton*km,fTheta)*cm2/g;
MyData()->Add("Xmax",fXmax);
MyData()->Add("XmaxProton",fXmaxP);
Msg(EsafMsg::Debug) << "Xmax: " << fXmax << "fAltitude: " << fAltitude << "XmaxProton: " << fXmaxP << "fAltitudeProton:" << fAltitudeProton << MsgDispatch;
return kTRUE;
}
//_____________________________________________________________________________
Bool_t XmaxModule::PostProcess() {
return kTRUE;
}
//_____________________________________________________________________________
Bool_t XmaxModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
fRecoRootEvent->GetRecoXmax().SetXmax(fXmax);
fRecoRootEvent->GetRecoXmax().SetXmaxProton(fXmaxP);
Double_t XmaxError = fRecoRootEvent->GetTruth().GetXmax() - fXmax;
Double_t XmaxErrorP = fRecoRootEvent->GetTruth().GetXmax() - fXmaxP;
fRecoRootEvent->GetRecoXmax().SetErrorXmax(XmaxError);
fRecoRootEvent->GetRecoXmax().SetErrorXmaxProton(XmaxErrorP);
fRecoRootEvent->GetRecoXmax().SetQuality(1);
return kTRUE;
}
//_____________________________________________________________________________
Bool_t XmaxModule::Done() {
Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
return kTRUE;
}
void XmaxModule::UserMemoryClean() {
}