// $Id: ProfileModule.cc,v 1.5 2005/10/21 13:30:08 moreggia Exp $
// Author: Anne Stutz Jun, 6 2005
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: ProfileModule *
* Package: <packagename> *
* Coordinator: <coordinator> *
* *
*****************************************************************************/
//_____________________________________________________________________________
//
// ProfileModule
//
// <extensive class description>
//
// Config file parameters
// ======================
//
// <parameter name>: <parameter description>
// -Valid options: <available options>
//
#include "ProfileModule.hh"
#include "RecoEvent.hh"
#include "RecoPhotonOnPupilData.hh"
#include <TMath.h>
#include <TF1.h>
#include "EConst.hh"
#include "KakimotoFluoCalculator.hh"
#include "EsafSpectrum.hh"
#include "O1_ClearSkyPropagator.hh" //DELETE
#include "TestGround.hh"
#include "Atmosphere.hh"
ClassImp(ProfileModule)
using namespace sou;
using namespace EConst;
// parametrization definitions
//_________________________________________________________________________________________
Double_t GIL(Double_t* x, Double_t* par) {
//
// number of electrons at given depth according to GIL parameterization
// GIL stands for Greizen-Ilina-Linsley which is the QGSJET model parameterization.
//
// par[0] : Nmax
// par[1] : Xmax
//
Double_t Nmax = par[0];
Double_t Xmax = par[1];
Double_t t0 = 37.15;
Double_t t = x[0]/t0;
Double_t t_max = Xmax/t0;
Double_t Age = 2*t / (t+t_max);
Double_t Ne = -1;
if(Age <= 0) cout << "<WARNING> Negative depth/age in GIL calculations" << endl;
Ne = Nmax * exp(t - t_max - 2*t*log(Age));
return Ne;
}
//_________________________________________________________________________________________
Double_t GFA(Double_t* x, Double_t* par) {
//
// number of electrons at given depth according to Gaussian Function in Age parameterization.
// Formulae from C. Song, Astropart. Physics 22(2004)151
//
// par[0] : Nmax
// par[1] : Xmax
// par[2] : sigma
//
Double_t Nmax = par[0];
Double_t Xmax = par[1];
Double_t sigma = par[2];
Double_t Age = 3*x[0] / (x[0] + 2*Xmax);
Double_t Ne = -1;
if(Age <= 0) cout << "<WARNING> Negative depth/age in GFA calculations" << endl;
Ne = Nmax * exp((-1./(2.*sigma*sigma)) * pow(Age - 1.,2));
return Ne;
}
//_________________________________________________________________________________________
Double_t GHF(Double_t* x, Double_t* par) {
//
// number of electrons at given depth according to Gaisser Hillas parametrisation
// Formulae from C. Song, Astropart. Physics 22(2004)151
//
// par[0] : Nmax
// par[1] : Xmax
// par[2] : lambda
// par[3] : X0
//
Double_t Nmax = par[0];
Double_t Xmax = par[1];
Double_t lambda = par[2];
Double_t X0 = par[3];
Double_t X = x[0];
Double_t Ne = -1;
Ne = Nmax * pow((X-X0)/(Xmax-X0),(Xmax-X0)/lambda) * exp((Xmax-X)/lambda);
return Ne;
}
//_____________________________________________________________________________
ProfileModule::ProfileModule() : RecoModule("Profile") {
//
// ctor
//
}
//_____________________________________________________________________________
ProfileModule::~ProfileModule() {
//
// dtor
//
}
//_____________________________________________________________________________
Bool_t ProfileModule::Init() {
//
// initializations
//
Msg(EsafMsg::Info) << "Initializing " << MsgDispatch;
fPrevious = Conf()->GetStr( "ProfileModule.fPrevious" );
fParamFormula = Conf()->GetStr("ProfileModule.fParamFormula");
if ( fPrevious == "truth" )
Msg(EsafMsg::Info) << " Data needed for ProfileModule read from Truth " << MsgDispatch;
else {
Msg(EsafMsg::Warning) << "Wrong config value ProfileModule.fPrevious " << fPrevious << " Setted to truth" << MsgDispatch;
fPrevious = "truth";
}
return kTRUE;
}
//_____________________________________________________________________________
Bool_t ProfileModule::PreProcess() {
//
//
//
fEvent = NULL;
fTheta = 0.;
fPhi = 0.;
fNGtu = 0;
fTmin = 0;
fTmax = 0;
fTimeAtMax = 0;
fLhmax = 0;
fSh = 0;
fLh = 0;
fNmax = 0;
fXmax = 0;
fSigma = 0;
fX0 = 0;
fLambda = 0;
fChiSquare = 0;
fFluoPhP = NULL;
fFluoAtTrack = NULL;
fEusoRadius = 1250*mm;
fEUSO.SetXYZ(0,0,430*km);
Msg(EsafMsg::Warning) <<"SHOULD check if true EUSO position at (0,0,430km) and Radius = 1250mm"<<MsgDispatch;
return kTRUE;
}
//_____________________________________________________________________________
Bool_t ProfileModule::Process(RecoEvent *ev) {
//
//
//
// manage the profile reconstruction chain
fEvent = ev;
// load information from previous modules
// amd buil
LoadModules();
// fit fluo profile on pupil (gaussian fit) -> get timeAtmax
FitFluoOnPupil();
// get track direction in MES (fTheta is not in MES)
// //TOFIX : fTheta should be defined at TOA, not at first interaction point
/*
EarthVector v1 = fEvent->GetHeader().GetTrueInitPos().Unit();
EarthVector vrot;
fTrackDir.SetXYZ(1,0,0);
EarthVector Uz(0,0,1);
vrot = Uz.Cross(v1);
fTrackDir.Rotate(v1.Theta(),vrot);
fTrackDir.Rotate(fPhi,v1);
vrot = v1.Cross(fTrackDir);
fTrackDir.Rotate(TMath::Pi()/2+fTheta,vrot);
fTrackDir.SetTheta(fTrackDir.Theta()+TMath::Pi());
*/
//TOFIX
EarthVector initpos = fEvent->GetHeader().GetTrueInitPos();
fTrackDir.SetMagThetaPhi(1.,(fFluoMaxPos - initpos).Theta(),(fFluoMaxPos - initpos).Phi());
// Define some quantities used to translate time on pupil into position along track
// origin of track coordinates == point corresponding to up-edge of first GTU (on pupil)
// firstly : calculate ShowerMax coordinate along track
// then : the rest follows
fLhmax = (EUSO() - fFluoMaxPos).Mag();
Double_t scalar = fTrackDir.Dot(EUSO() - fFluoMaxPos);
Double_t deltaT_c = (fTmin - fTimeAtMax) * Clight();
fShmax = (deltaT_c*deltaT_c + 2*deltaT_c*fLhmax) / (2*(fLhmax + deltaT_c - scalar));
fShmax *= -1.;
fSh = fShmax + scalar;
fImpactParameter = fFluoMaxPos + (fSh-fShmax)*fTrackDir;
fLh = (EUSO() - fImpactParameter).Mag();
cout<<"true theta phi = "<<fTheta/deg<<" "<<fPhi/deg<<endl;
cout<<"theta phi mag = "<<fTrackDir.Theta()/deg<<" "<<fTrackDir.Phi()/deg<<" "<<fTrackDir.Mag() <<endl;
cout<<"Shmax = "<<fShmax/km<<endl;
cout<<"Sh1 = "<<GetTrackCoordinate(fTmin)/km<<endl;
cout<<"Shend = "<<GetTrackCoordinate(fTmax)/km<<endl;
cout<<"showermaxpos.Zv = "<<fFluoMaxPos.Zv()/km<<endl;
cout<<"showermaxpos = "<<fFluoMaxPos<<endl;
// translate fluo transmitted on pupil into fluo generated along track
MakeFluoAtTrack();
// Iteration in order to estimate fluorescence photons emitted at shower track
//fPhFluoOnPupil = fPhOnPupil;
LightToSize();
// fit shower profile with required parametrization type
FitShowerProfile(fParamFormula);
return kTRUE;
}
//_____________________________________________________________________________
void ProfileModule::LoadModules() {
//
//
//
// load information from previous modules or from truth
if ( fPrevious == "truth" ) {
fTheta = fEvent->GetHeader().GetTrueTheta();
fPhi = fEvent->GetHeader().GetTruePhi();
SetFluoMaxPos(fEvent->GetHeader().GetTrueFluoMaxPos());
Msg(EsafMsg::Debug)<<" Nb of true photons on pupil "<<fEvent->GetHeader().GetTruePhotonsOnPupil()<< MsgDispatch;
RecoPhotonOnPupilData *rph = 0;
fTmin = HUGE;
fTmax = 0;
Double_t GTUlength = fEvent->GetHeader().GetGtuLength();
for ( Int_t i=0;i< fEvent->GetHeader().GetTruePhotonsOnPupil();i++) {
rph = fEvent->GetRecoPhotonOnPupilData(i);
// fluo only, in a first time
if(rph->GetType()==0) {
if ( rph->GetTime()>fTmax ) fTmax = rph->GetTime();
if ( rph->GetTime()<fTmin ) fTmin = rph->GetTime();
}
}
fNGtu = Int_t((fTmax-fTmin)/GTUlength) + 1;
Msg(EsafMsg::Debug)<<" Nb of GTU "<<fNGtu <<" from "<< fTmin/microsecond <<" to "<<fTmax/microsecond<<MsgDispatch;
SetArraysSize(fNGtu+1);
//fGtuAxis.Set(fNGtu, fTmin, fTmin + fNGtu*GTUlength);
fFluoPhP = (TH1F*)gROOT->FindObject("fluoOnPupil");
if(!fFluoPhP) fFluoPhP = new TH1F("fluoOnPupil","fluoOnPupil",fNGtu,fTmin/microsecond,(fTmin + fNGtu*GTUlength)/microsecond);
for ( Int_t i=0;i< fEvent->GetHeader().GetTruePhotonsOnPupil();i++) {
rph = fEvent->GetRecoPhotonOnPupilData(i);
// fluo only, in a first time
if(rph->GetType()==0) fFluoPhP->Fill(rph->GetTime()/microsecond);
}
}
}
//_____________________________________________________________________________
void ProfileModule::MakeFluoAtTrack() {
//
// transform photons on pupil to photons emitted isotropically at shower track
//
// define Xaxis of fluo histo along track
for ( Int_t i=0; i<fNGtu;i++ ) {
fFluoAtTrack_array[i] = GetTrackCoordinate(fFluoPhP->GetBinLowEdge(i+1)*microsecond)/km;
}
fFluoAtTrack_array[fNGtu] = GetTrackCoordinate((fFluoPhP->GetBinLowEdge(fNGtu)+fFluoPhP->GetBinWidth(fNGtu))*microsecond)/km;
// build output histo
fFluoAtTrack = (TH1F*)gROOT->FindObject("fluoAtTrack");
if(!fFluoAtTrack) fFluoAtTrack = new TH1F("fluoAtTrack","fluoAtTrack",fNGtu,fFluoAtTrack_array.GetArray());
// fill output histo
EarthVector pos;
Double_t trans = 1.;
for ( Int_t i=0; i<fNGtu;i++ ) {
// Get mean position in space corresponding to the GTU bin
pos = fFluoMaxPos + (fFluoAtTrack->GetBinCenter(i+1)*km - fShmax)*fTrackDir;
trans = FluoMeanTransmission(pos);
fFluoAtTrack->SetBinContent(i+1,fFluoPhP->GetBinContent(i+1)/EusoOmega(pos)/trans);
}
fFluoAtTrack->Write();//DELETE
// for debugging
/*
TArrayD fluoalt_array;
fluoalt_array.Set(fNGtu+1);
EarthVector alt;
for ( Int_t i=0; i<fNGtu;i++ ) {
alt = fFluoMaxPos + (fFluoAtTrack_array[i]*km - fShmax)*fTrackDir;
fluoalt_array[fNGtu - i] = alt.Zv()/km;
cout<<"altitude["<<fNGtu - i<<"] = "<<fluoalt_array[fNGtu - i]<<endl;
}
alt = fFluoMaxPos + (fFluoAtTrack_array[fNGtu]*km - fShmax)*fTrackDir;
fluoalt_array[0] = alt.Zv()/km;
cout<<"altitude[0] = "<<fluoalt_array[0]<<endl;
TH1F* fFluoAltitude = 0;
fFluoAltitude = (TH1F*)gROOT->FindObject("FluoAltitude");
if(!fFluoAltitude) fFluoAltitude = new TH1F("FluoAltitude","FluoAltitude",fNGtu,fluoalt_array.GetArray());
for ( Int_t i=0; i<fNGtu;i++ ) {
fFluoAltitude->SetBinContent(fNGtu - i,fFluoAtTrack->GetBinContent(i+1));
}
fFluoAltitude->Write();
*/
}
//_____________________________________________________________________________
void ProfileModule::LightToSize() {
//
// calculate number of charged particles in each depth bin
//
// get atmosphere description
const Atmosphere* atmo = Atmosphere::Get();
// build output histo Xaxis
EarthVector pos;
for ( Int_t i=0; i<fNGtu+1;i++ ) {
pos = fFluoMaxPos + (fFluoAtTrack_array[i]*km - fShmax)*fTrackDir;
fShowerSize_array[i] = atmo->Grammage(pos,-fTrackDir,"dir")*cm2/g;
}
// build output histo
fShowerSize = (TH1F*)gROOT->FindObject("ShowerSize");
if(!fShowerSize) fShowerSize = new TH1F("ShowerSize","ShowerSize",fNGtu,fShowerSize_array.GetArray());
// fill output histo
for ( Int_t i=0; i<fNGtu;i++ ) {
fShowerSize->SetBinContent(i+1,fFluoAtTrack->GetBinContent(i+1)/(4 * fFluoAtTrack->GetBinWidth(i+1)/m));
}
fShowerSize->Write();
}
//_____________________________________________________________________________
Double_t ProfileModule::GetTrackCoordinate(Double_t t_on_p) {
//
// Get track coordinate corresponding to a given time on pupil
// Coordinate origin corresponds to the up-edge of first GTU
//
Double_t rtn = 0;
Double_t deltaT_c = 0;
// get position on track as function of time on pupil
deltaT_c = (t_on_p - fTimeAtMax) * Clight();
rtn = (fLhmax + fShmax + deltaT_c)*(fLhmax + fShmax + deltaT_c) - fLh*fLh - fSh*fSh;
rtn /= 2*(deltaT_c + fLhmax + fShmax - fSh);
return rtn;
}
//_____________________________________________________________________________
void ProfileModule::FitFluoOnPupil() {
//
// fit fluo profile on pupil with a gaussian
// get time of maximum, nb of photons at maximum, profile width (sigma)
//
Double_t* par = new Double_t[3];
TF1* mygaus = new TF1("mygaus","gaus"); //TOFIX : to delete
fFluoPhP->Fit("mygaus","Q");
mygaus->GetParameters(par);
fFluoPhP->Write(); //DELETE
fTimeAtMax = par[1]*microsecond;
delete[] par;
par = 0;
}
//_____________________________________________________________________________
void ProfileModule::FitShowerProfile(string type) {
//
// fit shower profile along track
//
Double_t* par = new Double_t[4];
TF1* func = 0;
Double_t uplim = fShowerSize->GetBinLowEdge(fShowerSize->GetNbinsX()) + fShowerSize->GetBinWidth(fShowerSize->GetNbinsX());
Double_t Nmax = fShowerSize->GetMaximum(); //init value
Double_t Xmax = fShowerSize->GetBinCenter(fShowerSize->GetMaximumBin()); //init value
if(fParamFormula == "GIL") func = new TF1("ShowerFit",GIL,0.,uplim,2);
//if(fParamFormula == "GHF") func = new TF1("ShowerFit",GHF,0.,uplim,3);
//if(fParamFormula == "GFA") func = new TF1("ShowerFit",GFA,0.,uplim,4);
cout <<"INITIAL Nmax Xmax = "<<Nmax<<" "<<Xmax<<endl;
func->SetParameters(Nmax,Xmax);
func->SetParLimits(0,1e10,1e20); //TOFIX
func->SetParLimits(1,400,2000); //TOFIX
fShowerSize->Fit("ShowerFit","");
func->GetParameters(par);
fNmax = par[0];
fXmax = par[1];
fChiSquare = func->GetChisquare();
fShowerSize->Write(); //DELETE
cout <<"Nmax Xmax Chisquare = "<<fNmax<<" "<<fXmax<<" "<<fChiSquare<<endl;
delete[] par;
par = 0;
}
//_____________________________________________________________________________
Double_t ProfileModule::FluoMeanTransmission(const EarthVector& pos) {
//
// calculate mean transmission from a given position on showertrack
//
// initializations
// use of simulation classes -> //TOFIX
KakimotoFluoCalculator fluocalcul;
EsafSpectrum fluospec;
TestGround gnd;
O1_ClearSkyPropagator* transcalcul = new O1_ClearSkyPropagator(&gnd);
DetectorGeometry dum;
dum.SetRadius(EusoRadius());
dum.SetPos(EUSO());
transcalcul->CopyDetectorGeometry(&dum,false,false);
// get fluo spectrum and yield at pos
// electrons energy taken constant at 80MeV
Double_t yield = fluocalcul.GetFluoYield(pos.Zv(),80*MeV,&fluospec);
Double_t meantrans = transcalcul->GoToDetector(fluospec,pos);
delete transcalcul;
transcalcul = 0;
return meantrans;
}
//_____________________________________________________________________________
Bool_t ProfileModule::PostProcess() {
//
//
//
return kTRUE;
}
//_____________________________________________________________________________
Bool_t ProfileModule::Done() {
//
//
//
cout << "ProfileModule done. " << endl;
return kTRUE;
}
//_____________________________________________________________________________
void ProfileModule::UserMemoryClean() {
//
// delete user objects
//
SafeDelete(fFluoPhP);
SafeDelete(fFluoAtTrack);
}
//_____________________________________________________________________________
Bool_t ProfileModule::SaveRootData(RecoRootEvent *fRecoRootEvent) {
//
//
//
return kTRUE;
}