source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/tools/src/EsafSpectrum.cc @ 114

Last change on this file since 114 was 114, checked in by moretto, 11 years ago

actual version of ESAF at CCin2p3
File size: 11.9 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: EsafSpectrum.cc 2927 2011-06-16 18:02:41Z mabl $
3// R.Pesce created Oct, 31 2003
4
5#include "EsafSpectrum.hh"
6#include "EsafRandom.hh"
7#include <stdexcept>
8#include <TMath.h>
9using namespace TMath;
10using namespace sou;
11
12
13ClassImp(EsafSpectrum)
14
15//______________________________________________________________________________
16EsafSpectrum::EsafSpectrum() : fGraph(0), fName("Empty") {
17    //
18    // Empty spectrum constructor
19    //
20    SetArraysSize(0);
21}
22
23//______________________________________________________________________________
24EsafSpectrum::EsafSpectrum(Float_t wl) : fGraph(0), fName("") {
25    //
26    // constructor
27    // call with monochromatic wavelength
28    //
29
30    BuildMono(wl);
31}
32
33//______________________________________________________________________________
34EsafSpectrum::EsafSpectrum(Float_t wlmin, Float_t wlmax) : fGraph(0), fName("") {
35    //
36    // constructor
37    // call with monochromatic wavelength
38    //
39
40    BuildFlat(wlmin,wlmax);
41}
42
43//______________________________________________________________________________
44EsafSpectrum::EsafSpectrum(const char* name) : fGraph(0), fName(""){
45    // call with string
46   
47    fName = "";
48
49    if (strncmp(name,"test",4)==0) BuildTest();
50    else if (strncmp(name,"nitro10km",9)==0) {
51        TFormula nitro("nitrogen - 10km","gaus(0)+gaus(3)+gaus(6)");
52        nitro.SetParameter(0, 49);
53        nitro.SetParameter(1, 337*nm);
54        nitro.SetParameter(2, 0.05*nm);
55        nitro.SetParameter(3, 36);
56        nitro.SetParameter(4, 357*nm);
57        nitro.SetParameter(5, 0.05*nm);
58        nitro.SetParameter(6, 15);
59        nitro.SetParameter(7, 391*nm);
60        nitro.SetParameter(8, 0.05*nm);
61        BuildFormula(&nitro, 100, 300*nm, 400*nm);
62    }
63    else
64        MsgForm(EsafMsg::Panic,"EsafSpectrum: Cannot build myself ");
65}
66
67//______________________________________________________________________________
68EsafSpectrum::EsafSpectrum(const TH1F* histo) : fGraph(0), fName(""){
69    // call with histogram
70   
71    BuildHisto(histo);
72}
73
74//______________________________________________________________________________
75EsafSpectrum::EsafSpectrum(const TF1* func) : fGraph(0), fName(""){
76    //
77    // call with function
78    //
79   
80    BuildFunction(func);
81}
82
83
84//______________________________________________________________________________
85EsafSpectrum::EsafSpectrum(TFormula* formula, Int_t nbins, Float_t xmin, Float_t xmax) : fGraph(0), fName(""){
86    //
87    // Call with formula
88    //
89   
90    BuildFormula(formula, nbins, xmin, xmax);
91}
92
93//______________________________________________________________________________
94EsafSpectrum::EsafSpectrum(const EsafSpectrum& o) : fGraph(0), fName(""){
95    //
96    // copy ctor
97    //
98 
99    o.Copy(*this);
100
101}
102
103//______________________________________________________________________________
104void EsafSpectrum::Copy(EsafSpectrum& o) const{
105    //
106    // copy ctor
107    //
108    fSpectrum.Copy(o.fSpectrum);
109    fIntegral.Copy(o.fIntegral);
110    fAxis.Copy(o.fAxis);
111   
112    SafeDelete(o.fGraph);
113    if (fGraph) o.fGraph = (TGraph*)fGraph->Clone();
114
115    o.fName = fName;
116}
117
118
119//______________________________________________________________________________
120void EsafSpectrum::Reset(Float_t wl) {
121    //
122    // reset monochromatic
123    //
124
125    BuildMono(wl);
126}
127
128//______________________________________________________________________________
129void EsafSpectrum::Reset(Float_t wlmin, Float_t wlmax) {
130    //
131    // reset flat
132    //
133
134    BuildFlat(wlmin, wlmax);
135}
136
137//______________________________________________________________________________
138void EsafSpectrum::Reset(const TH1F* histo) {
139    //
140    // reset histogram
141    //
142   
143    BuildHisto(histo);
144}
145
146//______________________________________________________________________________
147void EsafSpectrum::Reset(const TF1* func) {
148    //
149    // reset function
150    //
151   
152    BuildFunction(func);
153}
154
155//______________________________________________________________________________
156void EsafSpectrum::Reset(TFormula* formula, Int_t nbins, Float_t xmin, Float_t xmax) {
157    //
158    // reset formula
159    //
160   
161    BuildFormula(formula, nbins, xmin, xmax);
162}
163
164//______________________________________________________________________________
165void EsafSpectrum::BuildMono( Float_t wl) {
166    //
167    // monochromatic
168    //
169   
170    fAxis.Set(1, wl - 0.1*nm, wl + 0.1*nm);
171
172    SetArraysSize(1);
173
174    fSpectrum[0] = 1.0 / fAxis.GetBinWidth(1);
175    fIntegral[0] = 0.0;
176    fIntegral[1] = 1.0;
177    fIntegral[2] = 1.0;
178
179    fName = Form("Monochromatic %.1f nm",wl/nm);
180}
181
182
183//______________________________________________________________________________
184void EsafSpectrum::BuildFlat( Float_t wlmin, Float_t wlmax ) {
185    //
186    // monochromatic
187    //
188   
189    if ( wlmax < wlmin) {
190        Float_t x = wlmin;
191        wlmin = wlmax;
192        wlmax = x;
193        MsgForm(EsafMsg::Debug,"Wlmin > Wlmax, values switched."); 
194    }
195   
196    Int_t nbins = 2; // must be > 1
197    fAxis.Set(nbins, wlmin, wlmax);
198
199    SetArraysSize(nbins);
200
201    fIntegral[0] = 0. ;
202   
203    for ( Int_t i(0); i<nbins; i++ ) {
204        fSpectrum[i] = 1.0 / (fAxis.GetBinWidth(i+1)*nbins);
205        fIntegral[i+1] = (i+1) / (Float_t)nbins ;
206    }
207   
208    fIntegral[nbins+1] = 1. ;
209
210    fName = Form("Flat %.1f-%.1f nm",wlmin/nm,wlmax/nm);
211}
212
213//______________________________________________________________________________
214void EsafSpectrum::BuildHisto(const TH1F* histo) {
215    //
216    // histogram
217    //
218
219    (*histo->GetXaxis()).Copy(fAxis);
220    Int_t nbins = fAxis.GetNbins();
221    Float_t integral = histo->Integral(1, nbins, "width");
222    SetArraysSize(nbins);
223
224    fIntegral[0] = 0. ;
225   
226    for (Int_t i(0); i<nbins; i++) {
227        fSpectrum[i] = histo->GetBinContent(i+1)/integral;
228        fIntegral[i+1] = histo->Integral(1,i+1,"width")/integral;
229    }
230   
231    fIntegral[nbins+1] = 1. ;
232
233    fName = "Histogram ";
234    fName += histo->GetName();
235}
236
237//______________________________________________________________________________
238void EsafSpectrum::BuildFunction(const TF1* func) {
239    //
240    // function
241    //
242
243    TH1* HistoFunc = func->GetHistogram();
244    (*HistoFunc->GetXaxis()).Copy(fAxis);
245    Int_t nbins = fAxis.GetNbins();
246    Float_t integral = HistoFunc->Integral(1, nbins, "width");
247    SetArraysSize(nbins);
248
249    fIntegral[0] = 0. ;
250
251    for (Int_t i(0); i<nbins; i++) {
252        fSpectrum[i] = HistoFunc->GetBinContent(i+1)/integral;
253        fIntegral[i+1] = HistoFunc->Integral(1,i+1,"width")/integral;
254    } 
255
256    fIntegral[nbins+1] = 1. ;
257
258    fName = "Function ";
259    fName += func->GetName();
260}
261
262//______________________________________________________________________________
263void EsafSpectrum::BuildFormula(TFormula* formula, Int_t nbins, Float_t xmin, Float_t xmax) {
264    //
265    // formula
266    //
267   
268    Float_t integral=0;
269    fAxis.Set(nbins, xmin, xmax);
270    SetArraysSize(nbins);
271
272    for (Int_t i=0; i<nbins; i++) integral+=formula->Eval(fAxis.GetBinCenter(i+1))*fAxis.GetBinWidth(i+1);
273
274    fIntegral[0] = 0. ;
275   
276    for (Int_t i(0); i<nbins; i++) {
277        fSpectrum[i] = formula->Eval(fAxis.GetBinCenter(i+1))/integral;
278        fIntegral[i+1] =  fIntegral[i] + fSpectrum[i]*fAxis.GetBinWidth(i+1);
279    }
280
281    fIntegral[nbins+1] = 1. ;
282
283    fName = "Formula ";
284    fName += formula->GetName();
285}
286
287
288
289//______________________________________________________________________________
290EsafSpectrum::~EsafSpectrum() {
291    //
292    // destructor
293    //
294}
295
296//______________________________________________________________________________
297Float_t EsafSpectrum::GetLambda() const {
298    //
299    // generate a random number from fIntegral
300    //
301    Int_t nbins = fAxis.GetNbins();
302
303    if (nbins == 1) return fAxis.GetBinCenter(1);
304    else {
305        TRandom *rndm = EsafRandom::Get();
306        Float_t r1 = rndm->Rndm();
307
308        // +1 added because fIntegral's size is nbins, not nbins+2 like in TH1
309        Int_t ibin = BinarySearch(nbins,fIntegral.GetArray(),r1);
310        Float_t a = fAxis.GetBinLowEdge(ibin+1);
311        Stat_t bw = fAxis.GetBinWidth(ibin+1);
312        Stat_t i0 = fIntegral[ibin];
313        Stat_t i1 = fIntegral[ibin+1];
314
315        Stat_t d = 0;
316
317        if ( ibin == -1 ) { 
318            // FIXME should not happen
319            // this may happen when fIntegral[0] > 0
320            // it is somehow equivalent to state that fIntegral[-1]=0
321            for(Int_t i(0); i<fIntegral.GetSize(); i++)
322              d = bw*r1/i0;
323        } else {
324            // standard case
325            if ( (i1-i0) < kTolerance ) 
326                d = 0;
327            else 
328                d = bw*(r1-i0)/(i1-i0);
329//            d = fAxis.GetBinWidth(ibin+1)*(r1-fIntegral[ibin])/(fIntegral[ibin+1]-fIntegral[ibin]);
330        }
331        if ( IsNaN(a) || IsNaN(d) ) {
332            FatalError("EsafSpectrum::GetLambda() does not converge?!");
333        }
334       
335        return a+d;
336    }
337}
338
339//______________________________________________________________________________
340Float_t EsafSpectrum::GetLambda(Int_t i) const {
341    //
342    // return a random wavelength within the considered bin (i = number of the bin - 1 !!... tables convention)
343    //
344    if(i < 0 || i > fAxis.GetNbins()-1)
345        FatalError("Index out of range in EsafSpectrum::GetLambda(Int_t)");
346    TRandom* rndm = EsafRandom::Get();
347    return GetLambdaEntry(i) + (2*rndm->Rndm() - 1)/2*fAxis.GetBinWidth(i+1);
348}
349
350
351//______________________________________________________________________________
352void EsafSpectrum::BuildTest() {
353    //
354    // Build a test spectrum
355    //
356    Float_t xmin=300.0;
357    Float_t xmax=400.0;
358    Int_t nbins=100;
359    Float_t integral=0;
360
361    fAxis.Set(nbins, xmin, xmax);
362
363    SetArraysSize(nbins);
364
365    for (Int_t i=0; i<nbins; i++)
366        integral+=TMath::Abs(TMath::Sin((GetLambdaEntry(i)-xmin)*(TMath::TwoPi()/(xmax-xmin))))*fAxis.GetBinWidth(i);
367
368    fIntegral[0] = 0. ;
369
370    for (Int_t i=0; i<nbins; i++) {
371        fSpectrum[i] = TMath::Abs(TMath::Sin((GetLambdaEntry(i)-xmin)*(TMath::TwoPi()/(xmax-xmin))))/integral;
372        fIntegral[i+1] =  fIntegral[i] + fSpectrum[i+1]*fAxis.GetBinWidth(i);
373    }
374
375    fIntegral[nbins+1] = 1. ;
376
377    fName = "Test Spectrum";
378}
379
380//______________________________________________________________________________
381void EsafSpectrum::Clear(Float_t val) {
382    //
383    // fill all the bins with the given value
384    //
385    for(Int_t i=0; i<fAxis.GetNbins(); i++) {
386        fSpectrum[i] = val / fAxis.GetBinWidth(i+1);
387//        if (i == 0)
388//            fIntegral[i+1] = fSpectrum[i]*fAxis.GetBinWidth(i+1);
389//        else
390            fIntegral[i+1] =  fIntegral[i] + fSpectrum[i]*fAxis.GetBinWidth(i+1);
391    }
392}
393
394//______________________________________________________________________________
395void EsafSpectrum::Draw() {
396    //
397    // draw fSpectrum
398    //
399    Int_t nbins = fAxis.GetNbins();
400
401    // empty spectrum, quit
402    if (!nbins) return;
403
404    SafeDelete(fGraph);
405    fGraph = new TGraph;
406   
407    fGraph->Set(nbins);
408
409    for(Int_t i=0; i<nbins; i++) 
410        fGraph->SetPoint(i, GetLambdaEntry(i), fSpectrum[i]);
411
412    fGraph->Draw("ALP");
413}
414
415
416//______________________________________________________________________________
417Float_t EsafSpectrum::ConvoluteAndNormalize(const Double_t* coeff) {
418    //
419    // Convolute the spectrum with an array of numbers. Returns the normalization coeff
420    //
421    Int_t i;
422    Int_t nb = fAxis.GetNbins();
423    Float_t integral = 0;
424 //   Float_t prev = 0;
425    for(i=0; i<nb; i++) integral += fSpectrum[i]*Float_t(coeff[i])*fAxis.GetBinWidth(i+1);
426    for(i=0; i<nb; i++) {
427        if(integral == 0) FatalError("EsafSpectrum::Convolute() impossible if weights are equal to zero");
428        fSpectrum[i] *= Float_t(coeff[i])/integral;
429//        if(i) prev = fIntegral[i-1];
430//        prev = fIntegral[i-1];
431        fIntegral[i+1] = fIntegral[i] + fSpectrum[i]*fAxis.GetBinWidth(i+1);
432    }
433    return integral;
434}
435
436//______________________________________________________________________________
437void EsafSpectrum::Convolute(const Double_t* coeff) {
438    //
439    // convolute without renormalization
440    //
441    for(Int_t i=0; i<fAxis.GetNbins(); i++)
442        fSpectrum[i] *= Float_t(coeff[i]);
443}
Note: See TracBrowser for help on using the repository browser.