Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

N2dataFluoCalculator - source file

// $Id: N2dataFluoCalculator.cc,v 1.5 2005/10/13 14:00:17 berat Exp $
// Author: Berat   2005/01/10

/*****************************************************************************
 * ESAF: Euso Simulation and Analysis Framework                              *
 *                                                                           *
 *  Id: N2dataFluoCalculator                                                           *
 *  Package: <packagename>                                                   *
 *  Coordinator: <coordinator>                                               *
 *                                                                           *
 *****************************************************************************/

//_____________________________________________________________________________
//
// N2dataFluoCalculator
//
// To compute yield from data on molecular...
//
//   Config file parameters
//   ======================
//
//   <parameter name>: <parameter description>
//   -Valid options: <available options>
//

#include "N2dataFluoCalculator.hh"
#include "Atmosphere.hh"
#include "NumbersFileParser.hh"
#include "ConfigFileParser.hh"
#include <TProfile.h>
#include <TMath.h>
#include "TGraph.h"
#include "EsafSpectrum.hh"
#include "EConst.hh"

using namespace sou;
using namespace TMath;
using EConst::STP_Pressure; 

ClassImp(N2dataFluoCalculator)

//_____________________________________________________________________________
 N2dataFluoCalculator::N2dataFluoCalculator() {
    //
    // Constructor
    //
    fName = "N2data";
    Msg(EsafMsg::Info) << "N2dataFluoCalculator built" << MsgDispatch;
    SetFluoLines();


}

//_____________________________________________________________________________
 N2dataFluoCalculator::~N2dataFluoCalculator() {
    //
    // Destructor
    //
    SafeDelete(fN2dataFluo);
}
//____________________________________________________________________________
 void N2dataFluoCalculator::SetFluoLines() {
    //
    // Set tables using data file
    //

    fLines = new Double_t[100];
    Double_t xl,w;

    fWlim[0]=300.;fWlim[1]=430;
    

    NumbersFileParser n1n("config/LightSource/N2dataFluo/band1n",4);
    vector<Double_t> nu_1n   = n1n.GetCol(0);
    vector<Double_t> nup_1n  = n1n.GetCol(1);
    vector<Double_t> w_1n    = n1n.GetCol(2);
    vector<Double_t> xl_1n   = n1n.GetCol(3);

    NumbersFileParser n2p("config/LightSource/N2dataFluo/band2p",4);
    vector<Double_t> nu_2p   = n2p.GetCol(0);
    vector<Double_t> nup_2p  = n2p.GetCol(1);
    vector<Double_t> w_2p    = n2p.GetCol(2);
    vector<Double_t> xl_2p   = n2p.GetCol(3);
 
    Int_t nlines = 0, nu=0,nup=0;
    for(Int_t il=0; il<50; il++) {
             xl = xl_1n[il]/10.;
             if (xl<=0.) break;
             nu = (Int_t)nu_1n[il];
             nup= (Int_t)nup_1n[il];
             w = w_1n[il];
             fL1n[nu][nup]=xl;
             fW1n[nu][nup]=w_1n[il];
             if(xl>=fWlim[0] &&  xl<=fWlim[1]) {
	        fLines[nlines]=xl/1000000.;
                nlines++;}
             }

    for(Int_t il=0; il<50; il++) {
             xl = xl_2p[il]/10.;
             if (xl<=0.) break;
             nu = (Int_t)nu_2p[il];
             nup= (Int_t)nup_2p[il];
             w = w_2p[il];
             fL2p[nu][nup]=xl;
             fW2p[nu][nup]=w_2p[il];
             if(xl>=fWlim[0] &&  xl<=fWlim[1]) {
	        fLines[nlines]=xl/1000000.;
                nlines++;}
             }

    fN2Wavelengths=nlines;
        Msg(EsafMsg::Info)<<"Number of fluorescence lines ="<<nlines<< 
      " between " <<fWlim[0] <<" and " <<fWlim[1] <<" nm" << MsgDispatch;
	 
        TString formula = "gaus(0)";
        for (Int_t i(0); i<fN2Wavelengths; i++){
        formula += Form("+gaus(%d)",i*3);
	//    Msg(EsafMsg::Debug) << "  i " <<i<<"  fLines " << fLines[i]*1e06<<MsgDispatch; 
       }

        fN2dataFluo = new TFormula("fluo - n2data", formula ); 

	//	 Msg(EsafMsg::Debug)<<" formula :             "<<formula<<MsgDispatch;

}
 
//______________________________________________________________________________
 Double_t N2dataFluoCalculator::GetFluoYield(const Double_t N2_alt,
                   const Double_t N2_energy, EsafSpectrum* FluoSpectrum) const {
  //Msg(EsafMsg::Debug) << "N2 : alt= "<<N2_alt/km 
  //                      << " energy= " <<N2_energy/MeV << MsgDispatch;


    Double_t fTotalYield = 0.0;

    Double_t pcol2p[4],pcol1n[2];
    Double_t Tmes = 300*kelvin;
    TString coltype="NAG";
    TString atm_hum="USDRY";
    Double_t conversion= TwoPi() * EConst::HbarC()/(nm*eV);

 
    collision_quench_input(coltype, pcol2p, pcol1n);

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

    Double_t density=atmo->Air_Density(N2_alt);
    Double_t T=atmo->Temperature(N2_alt);
    Double_t P=atmo->Pressure(N2_alt);
    if(T<=0) FatalError("Invalid Temperature in N2dataFluoCalculator");
    //           Msg(EsafMsg::Debug) << " N2_alt "<<N2_alt<<" temp "<< T 
    //                 << " rho " << density <<"  p "<<P<<MsgDispatch;

    Double_t z = N2_alt/km; 
    Double_t dair=density/g*cm3;
    Double_t press=P/bar*1000.;      // conversion in bars
    Double_t p0 = STP_Pressure()/bar*1000.;// conversion in bars
    p0=p0/1000.;  // for compatibility P0 
    Double_t pr=press*1.e-03*760./p0;   //pression in mmHg

    Double_t quen1=0,quen2=0;
    quench2p(p0,z,T,atm_hum,quen2);
    quench1n(p0,z,T,atm_hum,quen1);
 

    Double_t depot,e_elec=32.3;             // To be confirmed
    depot=GetdEdX(e_elec)*1e6*dair*100.;  //deposit in eV over 1metre 
    if(depot>=e_elec*1e6)depot=e_elec*1e6;

    //        Msg(EsafMsg::Debug) << " altit "<<z<< "  P = " << press 
    //      << "  T = " << T << " rho = " <<dair<< " Depot = "<<depot
    //      <<"  quen2 = "<<quen2<< "  quen1 = "<<quen1<<MsgDispatch;   

 
    Double_t Yield[fN2Wavelengths];
    Int_t nl=0,nlvalid;         
    Double_t xphitot=0;
    Double_t xl,w,ephi,xphi,y2p,y1n;

	  // ***bande 2P***

	  Double_t snorm2p=5200*1.e-7;		//norme davidson 337
	  Double_t xphitot2p=0;
          TString flag="2P";
          for(Int_t nu=0;nu<5;nu++){
             for(Int_t nup=0;nup<10;nup++){
                   xl=fL2p[nu][nup];
                   if(xl>=fWlim[0] && xl<=fWlim[1]) {
                   w=fW2p[nu][nup];
	           ephi=conversion/xl;
	           xphi=w*depot*snorm2p/(quen2*ephi);
	           y2p=yield(pr,T,flag,nu,pcol2p, pcol1n, Tmes)
                         /yield(600.,293.,flag, nu, pcol2p, pcol1n, Tmes);
	           xphi=xphi*y2p;
	           xphitot2p=xphitot2p+xphi;
                   Yield[nl]=xphi/m;
                   nl++;}
	     }
	  }

	  // ***bande 1N***

	  Double_t snorm1n=102.e-7;
          Double_t xphitot1n=0;
	  flag="1N";	      
          for(Int_t nu=0;nu<5;nu++){
             for(Int_t nup=0;nup<5;nup++){
                   xl=fL1n[nu][nup];
                   if(xl>=fWlim[0] && xl<=fWlim[1]) {
                   w=fW1n[nu][nup];
	           ephi=conversion/xl;
	           xphi=w*depot*snorm1n/(quen1*ephi);
		   y1n=yield(pr,T,flag,nu,pcol2p, pcol1n, Tmes)
		     /yield(600.,293.,flag,nu,pcol2p, pcol1n, Tmes);
  	           xphi=xphi*y1n;
	           xphitot1n=xphitot1n+xphi;   
                   Yield[nl]=xphi/m;
                   nl++;}
	     }
	  }
          nlvalid=nl++;

	  xphitot=xphitot1n+xphitot2p;

	  //printf("%10.4e %10.4e %10.4e %10.4e %10.4e %10.4e %10.4e %10.4e %10.4e n",z,T ,press,pr,dair,quen1,quen2,
	  //              depot,xphitot);

	fTotalYield=xphitot/m;


        for(Int_t i=0; i<nlvalid; i++){
	  // 	  Msg(EsafMsg::Debug) << "  i " <<i<<"  fLines " << fLines[i]<< " yield   " <<Yield[i]<<MsgDispatch; 
	  fN2dataFluo->SetParameter(i*3,Yield[i]);
	  fN2dataFluo->SetParameter(i*3+1,fLines[i]);
          fN2dataFluo->SetParameter(i*3+2,0.05*nm);
        }

        if ( FluoSpectrum !=0 ) FluoSpectrum->Reset(fN2dataFluo, 130, 300*nm, 430*nm);

        //   Msg(EsafMsg::Debug) << "fTotalYield 1  "<< fTotalYield 
        //              <<"  alt "<< N2_alt/km <<"   energy "<<N2_energy  
        //              << "  N2 TotalYield/m "<<fTotalYield*m << MsgDispatch;

    return fTotalYield;
}

// __________________________________________________________________________
 Double_t N2dataFluoCalculator::GetFluoYield(const Double_t N2_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 wavelength spectrum and total yield for 1.4 MeV
    Double_t eneref = 1.4*MeV;
    Double_t RefYield = GetFluoYield(N2_alt,eneref,FluoSpectrum);

    Double_t Integral = EnergyDistribution->Integral(Emin,Emax);
   if (Integral == 0) return 0; 

    Double_t dEdX_reference = GetdEdX( eneref ); 
    

    // integral over the energy spectrum
    while(true) {
      //        fTotalYield += GetFluoYield(N2_alt,E,FluoSpectrum) * EnergyDistribution->Eval(E) * dE;
       fTotalYield += GetdEdX (E*MeV)*EnergyDistribution->Eval(E);
       E += dE;
       if(( Emax-Emin )<( E-Emin )) break;
    }
    fTotalYield *= RefYield/dEdX_reference*dE/Integral;
 
    //Msg(EsafMsg::Debug) << "fTotalYield  2 "<< fTotalYield  
    //                    << "N2 TotalYield/m "<<fTotalYield*m << MsgDispatch;

     return fTotalYield;
}
   

// __________________________________________________________________________
 Double_t  N2dataFluoCalculator::GetFluoYieldHisto(const Double_t N2_alt, const TH1F* ehisto , EsafSpectrum* FluoSpectrum) const{

    Double_t fTotalYield =0.;
    Double_t Int=ehisto->Integral();
    if (Int==0) return 0;


    for (Int_t i(0); i < ehisto->GetNbinsX(); i++) {
        Double_t E  = ehisto->GetBinCenter(i);
	Double_t dE = ehisto->GetBinWidth(i);  
        fTotalYield += GetFluoYield(N2_alt,E,FluoSpectrum) * (ehisto->GetBinContent(i))/Int * dE;
    }
       //Msg(EsafMsg::Debug) << "fTotalYield  (Histo) "<< fTotalYield 
       //                  << "N2 TotalYield/m "<<fTotalYield*m << MsgDispatch;

    return fTotalYield;
}


//___________________________________________________________________
// get the data necessary to compute collision quenching
 void N2dataFluoCalculator::collision_quench_input(TString typ, Double_t *pcol2p, Double_t *pcol1n) const{

	if(typ=="AIR"){
	  pcol2p[0]=15.;		// norm air 300K Bunner tab 2.2
	  pcol2p[1]=6.5;		 
	  pcol2p[2]=4.6;		 
	  pcol2p[3]=2.5;		 
	  pcol1n[0]=1.08;
	  pcol1n[1]=1.08;
	}
	if(typ=="N90"){		//norm N2 300K
	  pcol2p[0]=90.;		 
	  pcol2p[1]=24.5;		 
	  pcol2p[2]=10.9;	 
	  pcol2p[3]=5.4; 
	  pcol1n[0]=1.49;
	  pcol1n[1]=1.49;
	}
	if(typ=="NAG"){		//norm N2 300K
	  pcol2p[0]=115.85*760/1013.25;            //0x
	  pcol2p[1]=94.21*760./1013.25;            //1x		 
	  pcol2p[2]=66.8*760/1013.25;		 
	  pcol2p[3]=5.4*760/1013.25;
	  pcol1n[0]=4.42*760/1013.25;
	  pcol1n[1]=3.32*760/1013.25;
       }
	return;
}

//_________________________________________________________________________
 Double_t N2dataFluoCalculator::yield(Double_t p, Double_t T, TString flag, Int_t n, Double_t *pcol2p, Double_t *pcol1n, Double_t Tmes) const{
	
        Int_t i=n;
	Double_t pcol;
        pcol=pcol2p[i]*sqrt(Tmes/T);		
	if(flag=="1N") pcol=pcol1n[i]*sqrt(Tmes/T);
	return 1./(1.+p/pcol);
}
//__________________________________________________________________________
 void N2dataFluoCalculator::quench2p(Double_t p0,Double_t z,Double_t T,
				    TString atm, Double_t &quench) const{

        const Double_t t0=47.0e-3,c_o2=0.209,c_n2=0.781;
        const Double_t attach=2000.;
        Double_t pw,t_o2,t_h2o, qo2, qh2o,qo2h2o;

        pw=2.504e6*exp(-5417/T)/p0;
	t_o2=attach*c_o2/c_n2;		//microsec-1
	t_h2o=attach*4.6*pw*hum(z,atm);
	qo2=1+t_o2*t0;
	qh2o=1+t_h2o*t0;
	qo2h2o=1+t0*(t_o2+t_h2o);
	quench=qo2h2o;

	return;
}
//__________________________________________________________________________
   void N2dataFluoCalculator::quench1n(Double_t p0,Double_t z,Double_t T,
				      TString atm, Double_t &quench) const{

        const Double_t c_o2=0.209,c_n2=0.781;
        const Double_t wn2=35.8,wo2=32.2;
        Double_t wair,air_ion,azot_ion;

        wair=c_n2*wn2+c_o2*wo2;
        air_ion=1/wair;		//nombre d'ions crees dans l'air
        azot_ion=1/wn2;		//nombre d'ions crees dans l'azote
        quench=((c_o2+c_n2)/c_n2);

        return;
}
//________________________________________________________________________
 Double_t N2dataFluoCalculator::hum(Double_t z, TString atm) const{

        if(atm=="USDRY"||atm=="OSVAL"){
             return 0;}

	if(atm=="UST76"){
	    if(z<=10.) return 0.5;			//standard
	    if(z>=10 && z<=15) return 0.5-(z-10)*(0.5)/5.;
	    if(z>=15)return 0;
	}

	if(atm=="USTRO"){				//dry
	  if(z<=3) return 0.75;			//tropical
	  if(z>=3. && z<13) return 0.75-(z-3)*(0.75)/10.; 
	  if(z>=13 && z<16) return (z-13)*0.3/3;
	  if(z>=16 && z<20)return (20-z)*0.3/4;
	  if(z>=20) return 0.;
        }
	
	if(atm=="MIDLS"){				
	  if(z<5)return 0.75-z*(0.45)/5.; 
	  if(z>=5 && z<10)return 0.3;
	  if(z>=10 && z<15)return (15-z)*0.3/5.;
	  if(z>=15)return 0.;
	}

	if(atm=="MIDLW"){			
	  if(z<12) return 0.75-z*0.75/12 ;
	  if(z>=12) return 0.;
        }
        return 0.;
}

//__________________________________________________________________________
// get the dE/dX for electrons of energy N2_energy
 Double_t N2dataFluoCalculator::GetdEdX(const Double_t N2_energy) const {
   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;
   for (Int_t i=0; i<33; i++) {
       if ( energy[i]*MeV<=N2_energy && N2_energy<= energy[i+1]*MeV ) 
	  dEdXE = dEdX[i] + (dEdX[i+1]-dEdX[i])*(N2_energy-energy[i]*MeV)/(energy[i+1]*MeV-energy[i]*MeV); 
   }
   //   Msg(EsafMsg::Debug) << " N2_energy  "<< N2_energy/MeV <<" dEdXE  "<<dEdXE 
   //                  << " In Esaf units :" <<dEdXE*MeV/(g/cm2)<<MsgDispatch;

   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