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