// ESAF : Euso Simulation and Analysis Framework
// $Id: EsafSpectrum.cc,v 1.12 2005/04/28 08:23:40 thea Exp $
// R.Pesce created Oct, 31 2003
#include "EsafSpectrum.hh"
#include "EsafRandom.hh"
#include <stdexcept>
using namespace TMath;
ClassImp(EsafSpectrum)
//______________________________________________________________________________
EsafSpectrum::EsafSpectrum() : fGraph(0), fName("Empty") {
//
// Empty spectrum constructor
//
SetArraysSize(0);
}
//______________________________________________________________________________
EsafSpectrum::EsafSpectrum(Float_t wl) : fGraph(0), fName("") {
//
// constructor
// call with monochromatic wavelength
//
BuildMono(wl);
}
//______________________________________________________________________________
EsafSpectrum::EsafSpectrum(Float_t wlmin, Float_t wlmax) : fGraph(0), fName("") {
//
// constructor
// call with monochromatic wavelength
//
BuildFlat(wlmin,wlmax);
}
//______________________________________________________________________________
EsafSpectrum::EsafSpectrum(const char* name) : fGraph(0), fName(""){
// call with string
fName = "";
if (strncmp(name,"test",4)==0) BuildTest();
else if (strncmp(name,"nitro10km",9)==0) {
TFormula nitro("nitrogen - 10km","gaus(0)+gaus(3)+gaus(6)");
nitro.SetParameter(0, 49);
nitro.SetParameter(1, 337*nm);
nitro.SetParameter(2, 0.05*nm);
nitro.SetParameter(3, 36);
nitro.SetParameter(4, 357*nm);
nitro.SetParameter(5, 0.05*nm);
nitro.SetParameter(6, 15);
nitro.SetParameter(7, 391*nm);
nitro.SetParameter(8, 0.05*nm);
BuildFormula(&nitro, 100, 300*nm, 400*nm);
}
else
MsgForm(EsafMsg::Panic,"EsafSpectrum: Cannot build myself ");
}
//______________________________________________________________________________
EsafSpectrum::EsafSpectrum(const TH1F* histo) : fGraph(0), fName(""){
// call with histogram
BuildHisto(histo);
}
//______________________________________________________________________________
EsafSpectrum::EsafSpectrum(const TF1* func) : fGraph(0), fName(""){
//
// call with function
//
BuildFunction(func);
}
//______________________________________________________________________________
EsafSpectrum::EsafSpectrum(TFormula* formula, Int_t nbins, Float_t xmin, Float_t xmax) : fGraph(0), fName(""){
//
// Call with formula
//
BuildFormula(formula, nbins, xmin, xmax);
}
//______________________________________________________________________________
EsafSpectrum::EsafSpectrum(const EsafSpectrum& o) : fGraph(0), fName(""){
//
// copy ctor
//
o.Copy(*this);
}
//______________________________________________________________________________
void EsafSpectrum::Copy(EsafSpectrum& o) const{
//
// copy ctor
//
fSpectrum.Copy(o.fSpectrum);
fIntegral.Copy(o.fIntegral);
fAxis.Copy(o.fAxis);
SafeDelete(o.fGraph);
if (fGraph) o.fGraph = (TGraph*)fGraph->Clone();
o.fName = fName;
}
//______________________________________________________________________________
void EsafSpectrum::Reset(Float_t wl) {
//
// reset monochromatic
//
if (wl <= 400.0*nm && wl >= 300.0*nm) BuildMono(wl);
else
FatalError("EsafSpectrum: Wavelenght specified is not in the allowed range");
}
//______________________________________________________________________________
void EsafSpectrum::Reset(Float_t wlmin, Float_t wlmax) {
//
// reset flat
//
BuildFlat(wlmin, wlmax);
}
//______________________________________________________________________________
void EsafSpectrum::Reset(const TH1F* histo) {
//
// reset histogram
//
BuildHisto(histo);
}
//______________________________________________________________________________
void EsafSpectrum::Reset(const TF1* func) {
//
// reset function
//
BuildFunction(func);
}
//______________________________________________________________________________
void EsafSpectrum::Reset(TFormula* formula, Int_t nbins, Float_t xmin, Float_t xmax) {
//
// reset formula
//
BuildFormula(formula, nbins, xmin, xmax);
}
//______________________________________________________________________________
void EsafSpectrum::BuildMono( Float_t wl) {
//
// monochromatic
//
fAxis.Set(1, wl - 0.1*nm, wl + 0.1*nm);
SetArraysSize(1);
fSpectrum[0] = 1.0 / fAxis.GetBinWidth(1);
fIntegral[0] = 0.0;
fIntegral[1] = 1.0;
fIntegral[2] = 1.0;
fName = Form("Monochromatic %.1f nm",wl/nm);
}
//______________________________________________________________________________
void EsafSpectrum::BuildFlat( Float_t wlmin, Float_t wlmax ) {
//
// monochromatic
//
if ( wlmax < wlmin) {
Float_t x = wlmin;
wlmin = wlmax;
wlmax = x;
MsgForm(EsafMsg::Debug,"Wlmin > Wlmax, values switched.");
}
Int_t nbins = 2; // must be > 1
fAxis.Set(nbins, wlmin, wlmax);
SetArraysSize(nbins);
fIntegral[0] = 0. ;
for ( Int_t i(0); i<nbins; i++ ) {
fSpectrum[i] = 1.0 / (fAxis.GetBinWidth(i+1)*nbins);
fIntegral[i+1] = (i+1) / (Float_t)nbins ;
}
fIntegral[nbins+1] = 1. ;
fName = Form("Flat %.1f-%.1f nm",wlmin/nm,wlmax/nm);
}
//______________________________________________________________________________
void EsafSpectrum::BuildHisto(const TH1F* histo) {
//
// histogram
//
fAxis = *histo->GetXaxis();
Int_t nbins = fAxis.GetNbins();
Float_t integral = histo->Integral(1, nbins, "width");
SetArraysSize(nbins);
fIntegral[0] = 0. ;
for (Int_t i(0); i<nbins; i++) {
fSpectrum[i] = histo->GetBinContent(i+1)/integral;
fIntegral[i+1] = histo->Integral(1,i+1,"width")/integral;
}
fIntegral[nbins+1] = 1. ;
fName = "Histogram ";
fName += histo->GetName();
}
//______________________________________________________________________________
void EsafSpectrum::BuildFunction(const TF1* func) {
//
// function
//
TH1* HistoFunc = func->GetHistogram();
fAxis = *HistoFunc->GetXaxis();
Int_t nbins = fAxis.GetNbins();
Float_t integral = HistoFunc->Integral(1, nbins, "width");
SetArraysSize(nbins);
fIntegral[0] = 0. ;
for (Int_t i(0); i<nbins; i++) {
fSpectrum[i] = HistoFunc->GetBinContent(i+1)/integral;
fIntegral[i+1] = HistoFunc->Integral(1,i+1,"width")/integral;
}
fIntegral[nbins+1] = 1. ;
fName = "Function ";
fName += func->GetName();
}
//______________________________________________________________________________
void EsafSpectrum::BuildFormula(TFormula* formula, Int_t nbins, Float_t xmin, Float_t xmax) {
//
// formula
//
Float_t integral=0;
fAxis.Set(nbins, xmin, xmax);
SetArraysSize(nbins);
for (Int_t i=0; i<nbins; i++) integral+=formula->Eval(fAxis.GetBinCenter(i+1))*fAxis.GetBinWidth(i+1);
fIntegral[0] = 0. ;
for (Int_t i(0); i<nbins; i++) {
fSpectrum[i] = formula->Eval(fAxis.GetBinCenter(i+1))/integral;
fIntegral[i+1] = fIntegral[i] + fSpectrum[i]*fAxis.GetBinWidth(i+1);
}
fIntegral[nbins+1] = 1. ;
fName = "Formula ";
fName += formula->GetName();
}
//______________________________________________________________________________
EsafSpectrum::~EsafSpectrum() {
//
// destructor
//
}
//______________________________________________________________________________
Float_t EsafSpectrum::GetLambda() const {
//
// generate a random number from fIntegral
//
Int_t nbins = fAxis.GetNbins();
if (nbins == 1) return fAxis.GetBinCenter(1);
else {
TRandom *rndm = EsafRandom::Get();
Float_t r1 = rndm->Rndm();
// +1 added because fIntegral's size is nbins, not nbins+2 like in TH1
Int_t ibin = BinarySearch(nbins,fIntegral.GetArray(),r1);
Float_t a = fAxis.GetBinLowEdge(ibin+1);
if ( ibin == -1 ) {
// FIXME should not happen
// this may happen when fIntegral[0] > 0
// it is somehow equivalent to state that fIntegral[-1]=0
for(Int_t i(0); i<fIntegral.GetSize(); i++)
a += fAxis.GetBinWidth(ibin+1)*r1/(fIntegral[0]);
} else
// standard case
a += fAxis.GetBinWidth(ibin+1)*(r1-fIntegral[ibin])/(fIntegral[ibin+1]-fIntegral[ibin]);
return a;
}
}
//______________________________________________________________________________
Float_t EsafSpectrum::GetLambda(Int_t i) const {
//
// return a random wavelength within the considered bin (i = number of the bin - 1 !!... tables convention)
//
if(i < 0 || i > fAxis.GetNbins()-1)
FatalError("Index out of range in EsafSpectrum::GetLambda(Int_t)");
TRandom* rndm = EsafRandom::Get();
return GetLambdaEntry(i) + (2*rndm->Rndm() - 1)/2*fAxis.GetBinWidth(i+1);
}
//______________________________________________________________________________
void EsafSpectrum::BuildTest() {
//
// Build a test spectrum
//
Float_t xmin=300.0;
Float_t xmax=400.0;
Int_t nbins=100;
Float_t integral=0;
fAxis.Set(nbins, xmin, xmax);
SetArraysSize(nbins);
for (Int_t i=0; i<nbins; i++)
integral+=TMath::Abs(TMath::Sin((GetLambdaEntry(i)-xmin)*(TMath::TwoPi()/(xmax-xmin))))*fAxis.GetBinWidth(i);
fIntegral[0] = 0. ;
for (Int_t i=0; i<nbins; i++) {
fSpectrum[i] = TMath::Abs(TMath::Sin((GetLambdaEntry(i)-xmin)*(TMath::TwoPi()/(xmax-xmin))))/integral;
fIntegral[i+1] = fIntegral[i] + fSpectrum[i+1]*fAxis.GetBinWidth(i);
}
fIntegral[nbins+1] = 1. ;
fName = "Test Spectrum";
}
//______________________________________________________________________________
void EsafSpectrum::Clear(Float_t val) {
//
// fill all the bins with the given value
//
for(Int_t i=0; i<fAxis.GetNbins(); i++) {
fSpectrum[i] = val / fAxis.GetBinWidth(i+1);
// if (i == 0)
// fIntegral[i+1] = fSpectrum[i]*fAxis.GetBinWidth(i+1);
// else
fIntegral[i+1] = fIntegral[i] + fSpectrum[i]*fAxis.GetBinWidth(i+1);
}
}
//______________________________________________________________________________
void EsafSpectrum::Draw() {
//
// draw fSpectrum
//
Int_t nbins = fAxis.GetNbins();
// empty spectrum, quit
if (!nbins) return;
SafeDelete(fGraph);
fGraph = new TGraph;
fGraph->Set(nbins);
for(Int_t i=0; i<nbins; i++)
fGraph->SetPoint(i, GetLambdaEntry(i), fSpectrum[i]);
fGraph->Draw("ALP");
}
//______________________________________________________________________________
Float_t EsafSpectrum::ConvoluteAndNormalize(const Double_t* coeff) {
//
// Convolute the spectrum with an array of numbers. Returns the normalization coeff
//
Int_t i;
Int_t nb = fAxis.GetNbins();
Float_t integral = 0;
// Float_t prev = 0;
for(i=0; i<nb; i++) integral += fSpectrum[i]*Float_t(coeff[i])*fAxis.GetBinWidth(i+1);
for(i=0; i<nb; i++) {
if(integral == 0) FatalError("EsafSpectrum::Convolute() impossible if weights are equal to zero");
fSpectrum[i] *= Float_t(coeff[i])/integral;
// if(i) prev = fIntegral[i-1];
// prev = fIntegral[i-1];
fIntegral[i+1] = fIntegral[i] + fSpectrum[i]*fAxis.GetBinWidth(i+1);
}
return integral;
}
//______________________________________________________________________________
void EsafSpectrum::Convolute(const Double_t* coeff) {
//
// convolute without renormalization
//
for(Int_t i=0; i<fAxis.GetNbins(); i++)
fSpectrum[i] *= Float_t(coeff[i]);
}