Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

EsafSpectrum - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: EsafSpectrum.cc,v 1.14 2005/10/11 13:24:01 thea Exp $
// R.Pesce created Oct, 31 2003

#include "EsafSpectrum.hh"
#include "EsafRandom.hh"
#include <stdexcept>

using namespace TMath;
using namespace sou;


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);
        Stat_t bw = fAxis.GetBinWidth(ibin+1);
        Stat_t i0 = fIntegral[ibin];
        Stat_t i1 = fIntegral[ibin+1];

        Stat_t d;

        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++)
            d = bw*r1/i0;
        } else {
            // standard case
            if ( (i1-i0) < kTolerance ) 
                d = 0;
            else 
                d = bw*(r1-i0)/(i1-i0);
//            d = fAxis.GetBinWidth(ibin+1)*(r1-fIntegral[ibin])/(fIntegral[ibin+1]-fIntegral[ibin]);
        }
        if ( IsNaN(a) || IsNaN(d) ) {
            cout << "mavafffancuuuulo" << endl;
        }
        
        return a+d;
    }
}

//______________________________________________________________________________
 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]);
}
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 Wed Nov 16 16:10:10 2005