source: trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc@ 7825

Last change on this file since 7825 was 7784, checked in by tbretz, 18 years ago
*** empty log message ***
File size: 15.1 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz 5/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
19! Author(s): Marcos Lopez 10/2003 <mailto:marcos@gae.ucm.es>
20!
21! Copyright: MAGIC Software Development, 2000-2006
22!
23!
24\* ======================================================================== */
25
26//////////////////////////////////////////////////////////////////////////////
27//
28// MMcSpectrumWeight
29//
30// Change the spectrum of the MC showers simulated with Corsika (a power law)
31// to a new one, which can be either, again a power law but with a different
32// spectral index, or a generalizeed spectrum. The new spectrum can be
33// pass to this class in different ways:
34//
35// 1. Is the new spectrum will be a power law, just introduce the slope
36// of this power law.
37// 2. Is the new spectrum will have a general shape:
38// The new spectrum is passed as a char* (SetFormula())
39//
40// Method:
41// -------
42//
43// - Corsika spectrun: dN/dE = A * E^(a)
44// with a = fOldSlope, and A = N/integral{E*de} from ELowLim to EUppLim
45//
46// - New spectrum: dN/dE = B * g(E)
47// where B = N/integral{g*dE} from ELowLim to EUppLim, and N=NumEvents
48//
49// For converting the spectrum simulated with Corsika to the new one, we
50// apply a weight to each event, given by:
51//
52// W(E) = B/A * g(E)/E^(a)
53//
54// In the case the new spectrum is simply a power law: dN/dE = B * E^(b), we
55// have:
56//
57// W(E) = B/A * E^(b-a)
58//
59// (The factor B/A is used in order both the original and new spectrum have
60// the same area (i.e. in order they represent the same number of showers))
61//
62//
63// If using SetFormula you can specify formulas accepted by TF1, eg:
64// pow(X, -2.6)
65// (Rem: all capital (!) 'X' are replaced by the corresponding %s.fEnergy
66// automatically)
67//
68//
69// Input Containers:
70// MMcEvt
71// MMcCorsikaRunHeader
72//
73// Output Container:
74// MWeight [MParameterD]
75//
76//////////////////////////////////////////////////////////////////////////////
77#include "MMcSpectrumWeight.h"
78
79#include <TF1.h>
80#include <TH1.h>
81#include <TSpline.h>
82
83#include "MLog.h"
84#include "MLogManip.h"
85
86#include "MParList.h"
87#include "MParameters.h"
88
89#include "MHillas.h"
90#include "MPointingPos.h"
91
92#include "MMcEvt.hxx"
93#include "MMcCorsikaRunHeader.h"
94
95ClassImp(MMcSpectrumWeight);
96
97using namespace std;
98
99void MMcSpectrumWeight::Init(const char *name, const char *title)
100{
101 fName = name ? name : "MMcSpectrumWeight";
102 fTitle = title ? title : "Task to calculate weights to change the energy spectrum";
103
104 AddToBranchList("MMcEvt.fEnergy");
105
106 fNameWeight = "MWeight";
107 fNameMcEvt = "MMcEvt";
108
109 fNewSlope = -1;
110 fOldSlope = -1;
111
112 fEnergyMin = -1;
113 fEnergyMax = -2;
114
115 fNorm = 1;
116 fNormEnergy = 1;
117
118 fAllowChange = kFALSE;
119
120 fFunc = NULL;
121 fMcEvt = NULL;
122 fHillas = NULL;
123 fWeight = NULL;
124 fWeightsZd = NULL;
125 fWeightsSize = NULL;
126 fPointing = NULL;
127}
128
129// ---------------------------------------------------------------------------
130//
131// Default Constructor.
132//
133MMcSpectrumWeight::MMcSpectrumWeight(const char *name, const char *title)
134{
135 Init(name,title);
136}
137
138// ---------------------------------------------------------------------------
139//
140// Destructor. If necessary delete fFunc
141//
142MMcSpectrumWeight::~MMcSpectrumWeight()
143{
144 if (fFunc)
145 delete fFunc;
146// if (fWeightsSize)
147// delete fWeightsSize;
148}
149
150// ---------------------------------------------------------------------------
151//
152// Search for
153// - fNameMcEvt [MMcEvtBasic]
154//
155// Find/Create:
156// - fNameWeight [MWeight]
157//
158Int_t MMcSpectrumWeight::PreProcess(MParList *pList)
159{
160 fMcEvt = (MMcEvt*)pList->FindObject(fNameMcEvt, "MMcEvtBasic");
161 if (!fMcEvt)
162 {
163 *fLog << err << fNameMcEvt << " [MMcEvtBasic] not found... abort." << endl;
164 return kFALSE;
165 }
166
167 fWeight = (MParameterD*)pList->FindCreateObj("MParameterD", fNameWeight);
168 if (!fWeight)
169 return kFALSE;
170
171 if (fWeightsZd)
172 {
173 fPointing = (MPointingPos*)pList->FindObject("MPointingPos");
174 if (!fPointing)
175 {
176 *fLog << err << "MPointingPos not found... abort." << endl;
177 return kFALSE;
178 }
179 }
180
181 if (fWeightsSize)
182 {
183 fHillas = (MHillas*)pList->FindObject("MHillas");
184 if (!fHillas)
185 {
186 *fLog << err << "MHillas not found... abort." << endl;
187 return kFALSE;
188 }
189 }
190
191 return kTRUE;
192}
193
194// ---------------------------------------------------------------------------
195//
196// Replace {fNameMcEvt}.fEnergy by "(x)" and return the result.
197//
198TString MMcSpectrumWeight::ReplaceX(TString str) const
199{
200 return str.ReplaceAll(Form("%s.fEnergy", fNameMcEvt.Data()), "(x)");
201}
202
203// ---------------------------------------------------------------------------
204//
205// Return the function corresponding to the mc spectrum with
206// slope fOldSlope: pow({fNameMcEvt}.fEnergy, fOldSlope)
207//
208// The slope is returned as %.3f
209//
210TString MMcSpectrumWeight::GetFormulaSpecOld() const
211{
212 return Form("pow(%s.fEnergy, %.3f)", fNameMcEvt.Data(), fOldSlope);
213}
214
215// ---------------------------------------------------------------------------
216//
217// Return the function corresponding to the new spectrum with
218// slope fNewSlope: pow({fNameMcEvt}.fEnergy, fNewSlope)
219//
220// The slope is returned as %.3f
221//
222// If a different formula is set (SetFormula()) this formula is returned
223// unchanged.
224//
225TString MMcSpectrumWeight::GetFormulaSpecNew() const
226{
227 TString str = fFormula.IsNull() ? Form("pow(%s.fEnergy, %.3f)", fNameMcEvt.Data(), fNewSlope) : fFormula.Data();
228 if (!fFormula.IsNull())
229 str.ReplaceAll("X", Form("(%s.fEnergy)", fNameMcEvt.Data()));
230
231 return str;
232}
233
234// ---------------------------------------------------------------------------
235//
236// Return the formula to calculate weights.
237// Is is compiled by
238// o1 = integral(fEnergyMin, fEnergyMax, GetFormulaSpecOldX());
239// n1 = integral(fEnergyMin, fEnergyMax, GetFormulaSpecNewX());
240// o2 = CalcSpecOld(fNormEnergy);
241// n2 = CalcSpecNew(fNormEnergy);
242//
243// result (fNormEnergy<0):
244// fNorm*o1/n1*GetFormulaNewSpec()/GetFormulaOldSpec()
245//
246// result (fNormEnergy>0):
247// fNorm*o2/n2*GetFormulaNewSpec()/GetFormulaOldSpec()
248//
249// fNorm is 1 by default but can be overwritten using SetNorm()
250//
251// If the formulas GetFormulaSpecOldX() and GetFormulaSpecNewX()
252// are equal only fNorm is returned.
253//
254// The normalization constant is returned as %.16f
255//
256// Example: 0.3712780019*(pow(MMcEvt.fEnergy,-2.270))/(pow(MMcEvt.fEnergy,-2.600))
257//
258TString MMcSpectrumWeight::GetFormulaWeights() const
259{
260 if (GetFormulaSpecOld()==GetFormulaSpecNew())
261 return Form("%.16f", fNorm);
262
263 const Double_t iold = fNormEnergy<0 ? GetSpecOldIntegral() : CalcSpecOld(fNormEnergy);
264 const Double_t inew = fNormEnergy<0 ? GetSpecNewIntegral() : CalcSpecNew(fNormEnergy);
265
266 const Double_t norm = fNorm*iold/inew;
267
268 return Form("%.16f*(%s)/(%s)", norm, GetFormulaSpecNew().Data(), GetFormulaSpecOld().Data());
269}
270
271// ---------------------------------------------------------------------------
272//
273// Returns the integral between fEnergyMin and fEnergyMax of
274// GetFormulaSpecNewX() describing the destination spectrum
275//
276Double_t MMcSpectrumWeight::GetSpecNewIntegral() const
277{
278 TF1 funcnew("Dummy", GetFormulaSpecNewX());
279 return funcnew.Integral(fEnergyMin, fEnergyMax);
280}
281
282// ---------------------------------------------------------------------------
283//
284// Returns the integral between fEnergyMin and fEnergyMax of
285// GetFormulaSpecOldX() describing the simulated spectrum
286//
287Double_t MMcSpectrumWeight::GetSpecOldIntegral() const
288{
289 TF1 funcold("Dummy", GetFormulaSpecOldX());
290 return funcold.Integral(fEnergyMin, fEnergyMax);
291}
292
293// ---------------------------------------------------------------------------
294//
295// Returns the value of GetFormulaSpecNewX() at the energy e describing
296// the destination spectrum
297//
298Double_t MMcSpectrumWeight::CalcSpecNew(Double_t e) const
299{
300 TF1 funcnew("Dummy", GetFormulaSpecNewX());
301 return funcnew.Eval(e);
302}
303
304// ---------------------------------------------------------------------------
305//
306// Returns the value of GetFormulaSpecOldX() at the energy e describing
307// the simulated spectrum
308//
309Double_t MMcSpectrumWeight::CalcSpecOld(Double_t e) const
310{
311 TF1 funcnew("Dummy", GetFormulaSpecOldX());
312 return funcnew.Eval(e);
313}
314
315void MMcSpectrumWeight::SetWeightsSize(TH1D *h)
316{
317 fWeightsSize=h;
318 /*
319 if (h==0)
320 {
321 fWeightsSize=0;
322 return;
323 }
324
325 if (fWeightsSize)
326 delete fWeightsSize;
327
328 const Double_t xmin = TMath::Log10(h->GetXaxis()->GetXmin());
329 const Double_t xmax = TMath::Log10(h->GetXaxis()->GetXmax());
330 const Double_t xnum = h->GetNbinsX()+1;
331
332 fWeightsSize = new TSpline3("WeightsSize", xmin, xmax,
333 h->GetArray()+1, xnum);*/
334}
335
336// ---------------------------------------------------------------------------
337//
338// Initialize fEnergyMin, fEnergymax and fOldSlope from MMcCorsikaRunHeader
339// by GetELowLim(), GetEUppLim() and GetSlopeSpec().
340//
341// If fEnergyMax>fEnergyMin (means: the values have already been
342// initialized) and !fAllowChange the consistency of the new values
343// with the present values is checked with a numerical precision of 1e-10.
344// If one doesn't match kFALSE is returned.
345//
346// If the mc slope is -1 kFALSE is returned.
347//
348// If the new slope for the spectrum is -1 it is set to the original MC
349// slope.
350//
351// fFunc is set to the formula returned by GetFormulaWeightsX()
352//
353Bool_t MMcSpectrumWeight::Set(const MMcCorsikaRunHeader &rh)
354{
355 if (fEnergyMax>fEnergyMin && !fAllowChange)
356 {
357 if (TMath::Abs(fOldSlope-rh.GetSlopeSpec())>1e-10)
358 {
359 *fLog << err;
360 *fLog << "ERROR - The slope of the Monte Carlo is not allowed to change ";
361 *fLog << "(" << fOldSlope << " --> " << rh.GetSlopeSpec() << ")... abort." << endl;
362 return kFALSE;
363 }
364 if (TMath::Abs(fEnergyMin-rh.GetELowLim())>1e-10)
365 {
366 *fLog << err;
367 *fLog << "ERROR - The minimum simulated Monte Carlo energy is not allowed to change ";
368 *fLog << "(" << fEnergyMin << " --> " << rh.GetELowLim() << ")... abort." << endl;
369 return kFALSE;
370 }
371 if (TMath::Abs(fEnergyMax-rh.GetEUppLim())>1e-10)
372 {
373 *fLog << err;
374 *fLog << "ERROR - The maximum simulated Monte Carlo energy is not allowed to change (";
375 *fLog << "(" << fEnergyMax << " --> " << rh.GetEUppLim() << ")... abort." << endl;
376 return kFALSE;
377 }
378 return kTRUE;
379 }
380
381 //
382 // Sanity checks to be sure that we won't divide by zero later on
383 //
384 if (rh.GetSlopeSpec()==-1)
385 {
386 *fLog << err << "The MC Slope of the power law must be different of -1... exit" << endl;
387 return kFALSE;
388 }
389
390 fOldSlope = rh.GetSlopeSpec();
391 fEnergyMin = rh.GetELowLim();
392 fEnergyMax = rh.GetEUppLim();
393
394 if (fNewSlope==-1)
395 {
396 *fLog << inf << "The new slope of the power law is -1... no weighting applied." << endl;
397 fNewSlope = fOldSlope;
398 }
399
400 TString form(GetFormulaWeightsX());
401
402 if (fFunc)
403 delete fFunc;
404
405 fFunc = new TF1("", form);
406 gROOT->GetListOfFunctions()->Remove(fFunc);
407 fFunc->SetName("SpectralWeighs");
408
409 return kTRUE;
410}
411
412// ---------------------------------------------------------------------------
413//
414// The current contants are printed
415//
416void MMcSpectrumWeight::Print(Option_t *o) const
417{
418 *fLog << all << GetDescriptor() << endl;
419 *fLog << " Simulated energy range: " << fEnergyMin << "GeV - " << fEnergyMax << "GeV" << endl;
420 *fLog << " Simulated spectral slope: " << fOldSlope << endl;
421 *fLog << " New spectral slope: " << fNewSlope << endl;
422 *fLog << " User normalization: " << fNorm << endl;
423 *fLog << " Spectra are normalized: " << (fNormEnergy<0?"by integral":Form("at %.1fGeV", fNormEnergy)) << endl;
424 *fLog << " Old Spectrum: " << GetFormulaSpecOldX() << " (I=" << GetSpecOldIntegral() << ")" << endl;
425 *fLog << " New Spectrum: " << GetFormulaSpecNewX() << " (I=" << GetSpecNewIntegral() << ")" << endl;
426 if (fFunc)
427 *fLog << " Weight function: " << fFunc->GetTitle() << endl;
428}
429
430// ----------------------------------------------------------------------------
431//
432// Executed each time a new root file is loaded
433// We will need fOldSlope and fE{Upp,Low}Lim to calculate the weights
434//
435Bool_t MMcSpectrumWeight::ReInit(MParList *plist)
436{
437 MMcCorsikaRunHeader *rh = (MMcCorsikaRunHeader*)plist->FindObject("MMcCorsikaRunHeader");
438 if (!rh)
439 {
440 *fLog << err << "MMcCorsikaRunHeader not found... abort." << endl;
441 return kFALSE;
442 }
443
444 return Set(*rh);
445}
446
447// ----------------------------------------------------------------------------
448//
449// Fill the result of the evaluation of fFunc at fEvEvt->GetEnergy
450// into the weights container.
451//
452Int_t MMcSpectrumWeight::Process()
453{
454 const Double_t e = fMcEvt->GetEnergy();
455
456 Double_t w = 1;
457
458 if (fWeightsZd)
459 {
460 const Int_t i = fWeightsZd->GetXaxis()->FindFixBin(fPointing->GetZd());
461 w = fWeightsZd->GetBinContent(i);
462 }
463 if (fWeightsSize)
464 {
465 const Int_t i = fWeightsSize->GetXaxis()->FindFixBin(fHillas->GetSize());
466 w *= fWeightsSize->GetBinContent(i);
467 // w *= fWeightsSize->Eval(TMath::Log10(fHillas->GetSize()));
468 }
469
470 fWeight->SetVal(fFunc->Eval(e)*w);
471
472 return kTRUE;
473}
474
475// --------------------------------------------------------------------------
476//
477// Read the setup from a TEnv, eg:
478// MMcSpectrumWeight.NewSlope: -2.6
479// MMcSpectrumWeight.Norm: 1.0
480// MMcSpectrumWeight.NormEnergy: 200
481// MMcSpectrumWeight.Formula: pow(X, -2.6)
482//
483Int_t MMcSpectrumWeight::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
484{
485 Bool_t rc = kFALSE;
486 if (IsEnvDefined(env, prefix, "NewSlope", print))
487 {
488 rc = kTRUE;
489 SetNewSlope(GetEnvValue(env, prefix, "NewSlope", fNewSlope));
490 }
491 if (IsEnvDefined(env, prefix, "Norm", print))
492 {
493 rc = kTRUE;
494 SetNorm(GetEnvValue(env, prefix, "Norm", fNorm));
495 }
496 if (IsEnvDefined(env, prefix, "NormEnergy", print))
497 {
498 rc = kTRUE;
499 SetNormEnergy(GetEnvValue(env, prefix, "NormEnergy", fNormEnergy));
500 }
501 if (IsEnvDefined(env, prefix, "Formula", print))
502 {
503 rc = kTRUE;
504 SetFormula(GetEnvValue(env, prefix, "Formula", fFormula));
505 }
506
507 return rc;
508}
Note: See TracBrowser for help on using the repository browser.