Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

NaganoFluoCalculator - source file

// ESAF : Euso Simulation and Analysis Framework
// contact person : Naoto SAKAKI <n-sakaki@riken.jp>
// $Id: NaganoFluoCalculator.cc,v 1.20 2005/10/21 13:31:30 moreggia Exp $

#include <stdexcept>
#include "NaganoFluoCalculator.hh"
#include "EsafSpectrum.hh"
#include "Atmosphere.hh"
#include "EConst.hh"

#include <TMath.h>
#include <TProfile.h>
#include "TGraph.h"

using namespace TMath;
using namespace sou;
using namespace EConst;

ClassImp(NaganoFluoCalculator)

// number of fluorescence bins
const Int_t NaganoFluoCalculator::fNumWavelengths = 15;
const Double_t hPa = 100.*pascal;

// 391nm and 428nm are 1N bands and others are 2P bands
const Double_t NaganoFluoCalculator::fWavelength[]  = {
  316*nm, 329*nm, 337*nm, 354*nm, 358*nm,
  376*nm, 381*nm, 391*nm, 394*nm, 400*nm,
  406*nm, 414*nm, 420*nm, 427*nm, 428*nm};

const Double_t NaganoFluoCalculator::fPhi0[]  = {
  4.80e-4, 0.880e-4, 10.01e-4, 0.769e-4, 7.82e-4,
  1.20e-4, 2.46e-4,   9.60e-4, 0.42e-4,  0.847e-4,
  1.49e-4, 0.327e-4,  0.86e-4, 0.069e-4, 4.57e-4};

const Double_t NaganoFluoCalculator::fRefPressure[] = {
    23.0*hPa, 40.2*hPa, 19.2*hPa, 30.6*hPa, 18.1*hPa,
    34.1*hPa, 19.4*hPa,  5.02*hPa,24.2*hPa, 24.2*hPa,
    12.3*hPa, 19.3*hPa,  7.3*hPa, 72.*hPa, 3.86*hPa};


//______________________________________________________________________________
 NaganoFluoCalculator::NaganoFluoCalculator() {
    // 
    // Constructor
    //
    fName = "Nagano";
    Msg(EsafMsg::Info) << "NaganoFluoCalculator built" << MsgDispatch;

    TString formula = "gaus(0)";
    for (Int_t i(1); i<fNumWavelengths; i++)
        formula += Form("+gaus(%d)",i*3);
     
    fNaganoFluo = new TFormula("fluo - nagano", formula ); 
}

//______________________________________________________________________________
 NaganoFluoCalculator::~NaganoFluoCalculator() {
    // 
    // Destructor
    // 

    SafeDelete(fNaganoFluo);
}


//______________________________________________________________________________
 Double_t NaganoFluoCalculator::GetFluoYield(const Double_t NF_alt,
                   const Double_t NF_energy, EsafSpectrum* FluoSpectrum) const {
    //
    // Get the wavelength spectrum of the fluorescence emission as a function
    // of electron energy
    //
    // Input:
    //     Double_t NF_alt         the altitude where the electron is moving
    //     Double_t NF_energy      the energy of the electron
    // Output:
    //     fTotalYield             total fluorescence yield between 300 and 430 nm
    //                             per electron per unit length
    //     EsafSpectrum* FluoSpectrum fluorescence light spectrum per electron
    //                             per unit length
    //

    // [1] Nagano et al. Astroparticle Physics, 20 (2003) 293.(astro-ph/0303193)
    // [2] Nagano et al. astro-ph/0406474 (2004)
    // measurements of yield for 0.85MeV electrons

#ifdef DEBUG
    Msg(EsafMsg::Debug) << "NF "<<NF_alt/km << " " <<NF_energy/MeV << MsgDispatch;
#endif
    // Get Atmosphere
    const Atmosphere* atmo = Atmosphere::Get();

    Double_t density     = atmo->Air_Density(NF_alt);
    Double_t temperature = atmo->Temperature(NF_alt);
    if(temperature<=0)
        FatalError("Invalid Temperature in NaganoFluoCalculator");

    const Double_t dedx_ref = GetdEdX( 0.85*MeV );

    // density and temperature dependance from Nagano et al.
    const static Double_t A[] = {
        20.5*m2/kg,   3.91*m2/kg, 45.6*m2/kg, 3.68*m2/kg, 37.8*m2/kg,
        6.07*m2/kg, 12.7*m2/kg,  50.8*m2/kg, 2.25*m2/kg,  4.58*m2/kg,
        8.18*m2/kg,  1.83*m2/kg,  4.9*m2/kg, 0.40*m2/kg, 26.5*m2/kg
    };
    const static Double_t B[] = {
        2.14*m3/kg, 1.22*m3/kg, 2.56*m3/kg, 1.60*m3/kg,  2.72*m3/kg,
        1.44*m3/kg, 2.53*m3/kg, 9.80*m3/kg, 2.03*m3/kg,  2.03*m3/kg,
        3.99*m3/kg, 2.55*m3/kg, 6.8 *m3/kg, 0.68*m3/kg, 12.7*m3/kg
    };

    Double_t fTotalYield = 0.0;
    Double_t Yield[fNumWavelengths];

    for (Int_t i=0; i<fNumWavelengths; i++) {
        Yield[i] = GetdEdX( NF_energy ) / dedx_ref
            *density*A[i]/(1+density*B[i]*sqrt(temperature));
        fTotalYield += Yield[i];

        fNaganoFluo->SetParameter(i*3,Yield[i]);
        fNaganoFluo->SetParameter(i*3+1,fWavelength[i]);
        fNaganoFluo->SetParameter(i*3+2,0.05*nm);
    }
    if ( FluoSpectrum !=0 )
        FluoSpectrum->Reset(fNaganoFluo, 130, 300*nm, 430*nm);

#ifdef DEBUG
    Msg(EsafMsg::Debug) << "NF TotalYield/m "<<fTotalYield*m  << MsgDispatch;
#endif

    return fTotalYield;
}

//______________________________________________________________________________
 Double_t NaganoFluoCalculator::GetFluoYield_dE(const Double_t NF_alt, 
                       const Double_t NF_dE, EsafSpectrum* FluoSpectrum) const {
    // 
    // Get the wavelength spectrum of the fluorescence emission as a function
    // of energy deposit
    //
    // Input:
    //     Double_t NF_alt         the altitude where the electron(s) is(are) moving
    //     Double_t NF_dE         the energy deposit of particle(s)
    // Output:
    //     fTotalYield             total fluorescence yield between 300 and 430 nm
    //                             per electron per unit length
    //     EsafSpectrum* FluoSpectrum fluorescence light spectrum per electron
    //     per unit length
    //


    // Get Atmosphere
    const Atmosphere* atmo = Atmosphere::Get();
    Double_t temperature = atmo->Temperature(NF_alt);
    Double_t pressure    = atmo->Pressure(NF_alt);
    if(temperature<=0) FatalError("Invalid Temperature in NaganoFluoCalculator");
    
    // const Int_t nbin = 15;  // number of fluorescence bins
    Double_t Yield[fNumWavelengths];
    Double_t fTotalYield = 0.0;

    for (Int_t nwl=0; nwl<fNumWavelengths; nwl++) {
        Double_t RefPressure=fRefPressure[nwl]*sqrt(temperature/(293.0*kelvin));
        Double_t Phi=fPhi0[nwl]/(1+pressure/RefPressure);
        Yield[nwl]=NF_dE*Phi/(Hplanck()*Clight()/fWavelength[nwl]);
        fTotalYield+=Yield[nwl];

        fNaganoFluo->SetParameter(nwl*3,Yield[nwl]);
        fNaganoFluo->SetParameter(nwl*3+1,fWavelength[nwl]);
        fNaganoFluo->SetParameter(nwl*3+2,0.05*nm);
    }
    if ( FluoSpectrum !=0 )
        FluoSpectrum->Reset(fNaganoFluo, 130, 300*nm, 430*nm);

#ifdef DEBUG
    Msg(EsafMsg::Debug) << "NF TotalYield/m "<<fTotalYield*m  << MsgDispatch;
#endif

    return fTotalYield;
}

/* to be completed
//______________________________________________________________________________
Double_t NaganoFluoCalculator::GetdE(const Double_t NF_alt, 
                                     const EsafSpectrum* FluoSpectrum) const {
    // 
    // To be completed
    // 

    // Get Atmosphere
    const Atmosphere* atmo = Atmosphere::Get();
    Double_t temperature = atmo->Temperature(NF_alt);
    Double_t pressure    = atmo->Pressure(NF_alt);

    Double_t fEnergyDeposit = 0.0;
    Int_t SpectNbin = FluoSpectrum->GetAxis().GetNbins();
    for(Int_t i=0; i<SpectNbin; i++) {
        for(Int_t nwl=0; nwl<fNumWavelengths; nwl++) {
            if(FluoSpectrum->GetAxis().GetBinLowEdge(i) <= fWavelength[nwl] &&
                    fWavelength[nwl] < FluoSpectrum->GetAxis().GetBinLowEdge(i+1)  ) {
                Double_t RefPressure=fRefPressure[nwl]*sqrt(temperature/(293.0*kelvin));
                Double_t Phi=fPhi0[nwl]/(1+pressure/RefPressure);
                fEnergyDeposit += FluoSpectrum->GetWeight(i) * h_Planck * c_light / fWavelength[nwl] / Phi;
            }
        }
    }
    return fEnergyDeposit;
}
*/

//______________________________________________________________________________
 Double_t  NaganoFluoCalculator::GetFluoYield(const Double_t NF_alt, TF12* EnergyDistribution,
                                             EsafSpectrum* FluoSpectrum) const {
    //
    //
    //

    Double_t Emax = EnergyDistribution->GetXmax();  // unit to be fixed in shower (MeV)
    Double_t Emin = EnergyDistribution->GetXmin();
    Double_t E=Emin;
    Double_t dE = 0.5*MeV; 
    Double_t fTotalYield =0.;

    // get the wavelength spectrum and total yield for 1.4 MeV
    Double_t eneref = 1.4*MeV;
    Double_t RefYield = GetFluoYield(NF_alt,eneref,FluoSpectrum);
    Double_t Integral = EnergyDistribution->Integral(Emin,Emax);

    if (Integral == 0) return 0; 

    Float_t dEdX_reference = GetdEdX( eneref ); 
    
    // integral over the energy spectrum
    while(kTRUE) {
        fTotalYield += GetdEdX (E*MeV)*EnergyDistribution->Eval(E);
        E += dE;
        if(( Emax-Emin )<( E-Emin )) break;
    }
    fTotalYield *= RefYield/dEdX_reference*dE/Integral;
    return fTotalYield;
}
//_____________________________________________________________________________
 Double_t  NaganoFluoCalculator::GetFluoYieldHisto(const Double_t KF_alt, const TH1F* ehisto, EsafSpectrum* FluoSpectrum) const {
    Double_t fTotalYield =0.;

    // get the wavelenght spectrum and total yield for 1.4 MeV
    Double_t eneref = 1.4*MeV;
    Double_t RefYield = GetFluoYield(KF_alt,eneref,FluoSpectrum);
    Double_t Int = ehisto->Integral();

    if (Int == 0) return 0; 
    Float_t dEdX_reference = GetdEdX( eneref ); 
    for (Int_t i(0); i < ehisto->GetNbinsX(); i++) {
        Double_t E  = ehisto->GetBinCenter(i);
        Double_t dE = ehisto->GetBinWidth(i);  
        fTotalYield += RefYield *  GetdEdX (E*MeV)/dEdX_reference * (ehisto->GetBinContent(i))/Int * dE;
    }
    return fTotalYield;
}


//______________________________________________________________________________
 void NaganoFluoCalculator::PlotYield() {
    //
    // Plot the fluorescence yield along Nadir
    //
/*
    TProfile* TotalYield = (TProfile*)gROOT->FindObject("TotalYield_N");
    if ( !TotalYield ) TotalYield = new TProfile("TotalYield_N","Total Yield along Nadir",100,0,50,0,6);
    TProfile* Yield_337 = (TProfile*)gROOT->FindObject("YieldN_337");
    if ( !Yield_337 ) Yield_337 = new TProfile("YieldN_337","337nm Yield along Nadir",100,0,50,0,6);
    TProfile* Yield_357 = (TProfile*)gROOT->FindObject("YieldN_357");
    if ( !Yield_357 ) Yield_357 = new TProfile("YieldN_357","357nm Yield along Nadir",100,0,50,0,6);
    TProfile* Yield_391 = (TProfile*)gROOT->FindObject("YieldN_391");
    if ( !Yield_391 ) Yield_391 = new TProfile("YieldN_391","391nm Yield along Nadir",100,0,50,0,6);
   
    EsafSpectrum* spectrum = new EsafSpectrum(357*nm);    
   
    Double_t alt,total,y337,y357,y391;

    for (Double_t u=0; u<100;u++) {
        alt=u/2*km;
        Double_t energy = 80.*MeV;
        total=GetFluoYield(alt,energy,spectrum)*m;
        TotalYield->Fill(u/2,total);
        Int_t nb = spectrum->GetAxis().GetNbins();
        y337 = 0.;
        y357 = 0.;
        y391 = 0.;
        for ( Int_t i=0; i<nb; i++) {
          if ( 334*nm<=spectrum->GetAxis().GetBinLowEdge(i) && spectrum->GetAxis().GetBinLowEdge(i+1)<=338*nm )
          y337 = y337 + spectrum->GetWeight(i);     
          if ( 354*nm<=spectrum->GetAxis().GetBinLowEdge(i) && spectrum->GetAxis().GetBinLowEdge(i+1)<=358*nm ) 
          y357 = y357 + spectrum->GetWeight(i);     
          if ( 388*nm<=spectrum->GetAxis().GetBinLowEdge(i) && spectrum->GetAxis().GetBinLowEdge(i+1)<=392*nm ) 
          y391 = y391 + spectrum->GetWeight(i);     
        }  
        Yield_337->Fill(u/2 , y337);       
        Yield_357->Fill(u/2 , y357);
        Yield_391->Fill(u/2 , y391);
    }
*/
}

//______________________________________________________________________________
 Double_t NaganoFluoCalculator::GetdEdX(const Double_t NF_energy) const {
    //
    // Get the dE/dX for electrons of energy KF_energy
    // same function as in KakimotoFluoCalculator
    // 
   Double_t energy[33]={.01,.02,.03,.04,.06,.08,.1,.2,.3,.4,.5,.6,.7,.8,.9,1.,
                        1.25,1.5,2.,4.,6.,8.,10.,20.,40.,60.,80.,100.,200.,400.,
                        600.,800.,1000.};
   Double_t dEdX[33]={19.95,11.68,8.564,6.904,5.15,4.229,3.66,2.486,2.097,1.914,
                      1.813,1.753,1.716,1.693,1.679,1.67,1.665,1.67,1.693,1.799,
                      1.879,1.94,1.988,2.144,2.29,2.355,2.395,2.424,2.51,2.59,
                      2.633,2.66,2.681};

   Double_t dEdXE = 0;
   for (Int_t i=0; i<33; i++) {
       if ( energy[i]*MeV<=NF_energy && NF_energy<= energy[i+1]*MeV ) 
      dEdXE = dEdX[i] + (dEdX[i+1]-dEdX[i])*(NF_energy-energy[i]*MeV)/(energy[i+1]*MeV-energy[i]*MeV); 
   }
   return dEdXE*MeV/(g/cm2);
}
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