Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

ProfileModule - source file

// $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;
}



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