| 1 | #ifndef MARS_MMcSpectrumWeight
|
|---|
| 2 | #define MARS_MMcSpectrumWeight
|
|---|
| 3 |
|
|---|
| 4 | #ifndef MARS_MTask
|
|---|
| 5 | #include "MTask.h"
|
|---|
| 6 | #endif
|
|---|
| 7 |
|
|---|
| 8 | class TF1;
|
|---|
| 9 | class TH1;
|
|---|
| 10 | class TH1D;
|
|---|
| 11 | class TSpline3;
|
|---|
| 12 | class MParList;
|
|---|
| 13 | class MMcEvt;
|
|---|
| 14 | class MHillas;
|
|---|
| 15 | class MParameterD;
|
|---|
| 16 | class MPointingPos;
|
|---|
| 17 | class MMcCorsikaRunHeader;
|
|---|
| 18 |
|
|---|
| 19 | class MMcSpectrumWeight : public MTask
|
|---|
| 20 | {
|
|---|
| 21 | private:
|
|---|
| 22 | const MMcEvt *fMcEvt; // Pointer to the container with the MC energy
|
|---|
| 23 | const MHillas *fHillas;
|
|---|
| 24 | MParameterD *fWeight; // Pointer to the output MWeight container
|
|---|
| 25 | MPointingPos *fPointing;
|
|---|
| 26 |
|
|---|
| 27 | TString fNameWeight; // Name of the MWeight container
|
|---|
| 28 | TString fNameMcEvt; // Name of the MMcEvt container
|
|---|
| 29 |
|
|---|
| 30 | TF1 *fFunc; // Function calculating the weights
|
|---|
| 31 | TH1 *fWeightsZd; // Set additional ZA weights
|
|---|
| 32 | TH1 *fWeightsSize; // Set additional ZA weights
|
|---|
| 33 | // TSpline3 *fWeightsSize;
|
|---|
| 34 |
|
|---|
| 35 | Double_t fOldSlope; // Slope of energy spectrum generated with Corsika
|
|---|
| 36 | Double_t fNewSlope; // Slope of the new spectrum (if it is a power law)
|
|---|
| 37 |
|
|---|
| 38 | Double_t fEnergyMin; // Minimum energy simulated
|
|---|
| 39 | Double_t fEnergyMax; // Maximum energy simulated
|
|---|
| 40 |
|
|---|
| 41 | Double_t fNorm; // Normalization constant (additional normalization constant)
|
|---|
| 42 | Double_t fNormEnergy; // Energy at which the spectra are normalized (default -1 means the integral is used)
|
|---|
| 43 |
|
|---|
| 44 | TString fFormula; // Text Formula for new spectrum: eg. "pow(MMcEvt.fEnergy, -2.0)"
|
|---|
| 45 |
|
|---|
| 46 | Bool_t fAllowChange;
|
|---|
| 47 |
|
|---|
| 48 | // MMcSpectrumWeight
|
|---|
| 49 | void Init(const char *name, const char *title);
|
|---|
| 50 | TString ReplaceX(TString) const;
|
|---|
| 51 |
|
|---|
| 52 | // MTask
|
|---|
| 53 | Bool_t ReInit(MParList *plist);
|
|---|
| 54 | Int_t PreProcess(MParList *pList);
|
|---|
| 55 | Int_t Process();
|
|---|
| 56 |
|
|---|
| 57 | // MParContainer
|
|---|
| 58 | Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
|
|---|
| 59 |
|
|---|
| 60 | public:
|
|---|
| 61 | MMcSpectrumWeight(const char *name=NULL, const char *title=NULL);
|
|---|
| 62 | ~MMcSpectrumWeight();
|
|---|
| 63 |
|
|---|
| 64 | // Setter
|
|---|
| 65 | void SetNameWeight(const char *n="MWeight") { fNameWeight = n; }
|
|---|
| 66 | void SetNameMcEvt(const char *n="MMcEvt") { fNameMcEvt = n; }
|
|---|
| 67 | void SetNewSlope(Double_t s=-1) { fNewSlope = s; }
|
|---|
| 68 | void SetNorm(Double_t s=1) { fNorm = s; }
|
|---|
| 69 | void SetNormEnergy(Double_t s=1) { fNormEnergy = s; }
|
|---|
| 70 | void SetFormula(const char *f="") { fFormula = f; }
|
|---|
| 71 | void SetEnergyRange(Double_t min=-2, Double_t max=-1) { fEnergyMin=min; fEnergyMax=max; }
|
|---|
| 72 | void SetOldSlope(Double_t s=-2.6) { fOldSlope=s; }
|
|---|
| 73 | void SetWeightsZd(TH1 *h=0) { fWeightsZd = h; }
|
|---|
| 74 | void SetWeightsSize(TH1D *h=0);
|
|---|
| 75 |
|
|---|
| 76 | Bool_t Set(const MMcCorsikaRunHeader &h);
|
|---|
| 77 |
|
|---|
| 78 | // Getter
|
|---|
| 79 | TString GetFormulaSpecOld(const char *name) const;
|
|---|
| 80 | TString GetFormulaSpecNew(const char *name) const;
|
|---|
| 81 | TString GetFormulaWeights(const char *name) const;
|
|---|
| 82 | TString GetFormulaSpecOld() const { return GetFormulaSpecOld(fNameMcEvt); }
|
|---|
| 83 | TString GetFormulaSpecNew() const { return GetFormulaSpecNew(fNameMcEvt); }
|
|---|
| 84 | TString GetFormulaWeights() const { return GetFormulaWeights(fNameMcEvt); }
|
|---|
| 85 |
|
|---|
| 86 | TString GetFormulaSpecOldX() const { return ReplaceX(GetFormulaSpecOld()); }
|
|---|
| 87 | TString GetFormulaSpecNewX() const { return ReplaceX(GetFormulaSpecNew()); }
|
|---|
| 88 | TString GetFormulaWeightsX() const { return ReplaceX(GetFormulaWeights()); }
|
|---|
| 89 |
|
|---|
| 90 | Double_t GetSpecNewIntegral(Double_t emin, Double_t emax) const;
|
|---|
| 91 | Double_t GetSpecOldIntegral(Double_t emin, Double_t emax) const;
|
|---|
| 92 |
|
|---|
| 93 | Double_t GetSpecNewIntegral() const { return GetSpecNewIntegral(fEnergyMin, fEnergyMax); }
|
|---|
| 94 | Double_t GetSpecOldIntegral() const { return GetSpecOldIntegral(fEnergyMin, fEnergyMax); }
|
|---|
| 95 |
|
|---|
| 96 | Double_t CalcSpecNew(Double_t e) const;
|
|---|
| 97 | Double_t CalcSpecOld(Double_t e) const;
|
|---|
| 98 |
|
|---|
| 99 | Double_t GetEnergyMin() const { return fEnergyMin; }
|
|---|
| 100 | Double_t GetEnergyMax() const { return fEnergyMax; }
|
|---|
| 101 |
|
|---|
| 102 | // Functions
|
|---|
| 103 | void CompleteEnergySpectrum(TH1 &h, Double_t emin) const;
|
|---|
| 104 |
|
|---|
| 105 | // TObject
|
|---|
| 106 | void Print(Option_t *o="") const;
|
|---|
| 107 |
|
|---|
| 108 | ClassDef(MMcSpectrumWeight, 0) // Task to calculate weights to change the energy spectrum
|
|---|
| 109 | };
|
|---|
| 110 |
|
|---|
| 111 | #endif
|
|---|
