Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

HmaxForProtonModule - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: HmaxForProtonModule.cc,v 1.10 2005/10/02 14:18:08 thea 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 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
    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 ) td = ev->GetModuleData("TrackDirection2"); // try to get TrackDirection2
    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()  {
    
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005