Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

KakimotoFluoCalculator - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: KakimotoFluoCalculator.cc,v 1.20 2005/11/15 11:57:48 moreggia Exp $
// Anne Stutz created Mar, 23 2004

#include <stdexcept>
#include "KakimotoFluoCalculator.hh"
#include "EsafSpectrum.hh"
#include "Atmosphere.hh"
#include <TROOT.h>
#include <TProfile.h>
#include "TGraph.h"

using namespace sou;

ClassImp(KakimotoFluoCalculator)

//________________________________________________________________________________________________________________________
 KakimotoFluoCalculator::KakimotoFluoCalculator() : FluoCalculator() {
    //
    // ctor
    //
    
    fName = "Kakimoto";
    Msg(EsafMsg::Info) << "KakimotoluoCalculator built" << MsgDispatch;
    //    PlotYield();
}

//________________________________________________________________________________________________________________________
 KakimotoFluoCalculator::~KakimotoFluoCalculator() {
    //
    // dtor
    //
}

//________________________________________________________________________________________________________________________
 Double_t KakimotoFluoCalculator::GetFluoYield(const Double_t KF_alt, const Double_t KF_energy, EsafSpectrum* FluoSpectrum) const {
    //
    // get the wavelenght spectrum of the fluorescence emission
    //

    // Get Atmosphere
    const Atmosphere* atmo = Atmosphere::Get();

    // Kakimoto et al. measurements of yield at 760 mmHg, 15°C for 1.4 MeV electrons
    Int_t nbin = 3;  // number of fluorescence bins
   
    Double_t Wavelenght[] = {337*nm,357*nm,391*nm};  
    Double_t RefYield[]= {1.109/m,1.19/m,0.495/m};
    Double_t total2P = RefYield[0]+RefYield[1];
    Double_t total1N = RefYield[2];
 
    Double_t density = atmo->Air_Density(KF_alt);
    Double_t temp = atmo->Temperature(KF_alt);
    if(temp==0) 
       Msg(EsafMsg::Panic) << "Invalid Temperature in KakimotoFluoCalculator"<<MsgDispatch;
    
    // energy dependance 
    Double_t eneref = 1.4*MeV;
    Double_t energy = KF_energy;
    for( Int_t i=0; i<nbin; i++) RefYield[i] = GetdEdX (energy )/ GetdEdX( eneref ) * RefYield[i];

    // density and temperature dependance from Kakimoto et al. NIMA372(1996)527-533
    Double_t A1 = 89.0*m2/kg;   // +/- 1.7 m(2)/kg
    Double_t A2 = 55.0*m2/kg;   // 
    Double_t B1 = 1.85*m3/kg;   // +/- 0.04 m(3)/kg.K(-1/2)
    Double_t B2 = 6.50*m3/kg;   //

    Double_t Yield[nbin];
    Double_t fTotalYield = 0.0;

    string KF_string = "gaus(0)";
    char gaussian[10];
    for (Int_t i=1; i<nbin; i++) {
      sprintf(gaussian,"+gaus(%d)",i*3);
      KF_string += gaussian;
    }
    const char* KF_char = KF_string.c_str();
    TFormula* kakimoto = new TFormula("fluo - kakimoto",KF_char);

    for (Int_t i=0; i<nbin; i++) {
        if ( Wavelenght[i] != 391*nm ) 
	  Yield[i] =RefYield[i]/total2P
                     *density*A1/(1+density*B1*sqrt(temp));
        else 
	  Yield[i] = RefYield[i]/total1N
                     *density*A2/(1+density*B2*sqrt(temp));

        fTotalYield = fTotalYield + Yield[i];

        kakimoto->SetParameter(i*3,Yield[i]);
        kakimoto->SetParameter(i*3+1,Wavelenght[i]);
        kakimoto->SetParameter(i*3+2,0.05*nm);
    }
    if ( FluoSpectrum !=0 ) FluoSpectrum->Reset(kakimoto, 100, 300*nm, 400*nm);
    SafeDelete(kakimoto);
    
    return fTotalYield;
}

//________________________________________________________________________________________________________________________
 Double_t  KakimotoFluoCalculator::GetFluoYield(const Double_t KF_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 = (Emax-Emin)/1000.;
    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);

    // integral over the energy spectrum
    while(true) {
        fTotalYield += RefYield *  GetdEdX (E*MeV)/ GetdEdX( eneref ) * EnergyDistribution->Eval(E) * dE;
        E += dE;
        if(( Emax-Emin )<( E-Emin )) break;
    }
    fTotalYield /= EnergyDistribution->Integral(Emin,Emax);   // energy spectrum normalization

    return fTotalYield;
}

//________________________________________________________________________________________________________________________
 Double_t  KakimotoFluoCalculator::GetFluoYieldHisto(const Double_t KF_alt, const TH1F* ehisto, EsafSpectrum* FluoSpectrum) const {
    Double_t fTotalYield =0.;
    Double_t Int=ehisto->Integral();
    if (Int==0) return 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);

    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)/ GetdEdX( eneref ) * ehisto->GetBinContent(i) * dE;
    }
    fTotalYield /= Int;   // energy spectrum normalization

    return fTotalYield;
}

//________________________________________________________________________________________________________________________
 void KakimotoFluoCalculator::PlotYield() {
    //
    // plot the fluorescence yield along Nadir
    //
/*
    TProfile* TotalYield = (TProfile*)gROOT->FindObject("TotalYield");
    if ( !TotalYield ) TotalYield = new TProfile("TotalYield","Total Yield along Nadir",100,0,50,0,6);
    TProfile* Yield_337 = (TProfile*)gROOT->FindObject("Yield_337");
    if ( !Yield_337 ) Yield_337 = new TProfile("Yield_337","337nm Yield along Nadir",100,0,50,0,6);
    TProfile* Yield_357 = (TProfile*)gROOT->FindObject("Yield_357");
    if ( !Yield_357 ) Yield_357 = new TProfile("Yield_357","357nm Yield along Nadir",100,0,50,0,6);
    TProfile* Yield_391 = (TProfile*)gROOT->FindObject("Yield_391");
    if ( !Yield_391 ) Yield_391 = new TProfile("Yield_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 KakimotoFluoCalculator::GetdEdX(const Double_t KF_energy) const {
    //
    // get the dE/dX for electrons of energy KF_energy
    //
    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<=KF_energy && KF_energy<= energy[i+1]*MeV ) 
	   dEdXE = dEdX[i] + (dEdX[i+1]-dEdX[i])*(KF_energy-energy[i]*MeV)/(energy[i+1]*MeV-energy[i]*MeV); 
    }
    return dEdXE;
}
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