Index: trunk/MagicSoft/Mars/mbase/MTaskList.cc =================================================================== --- trunk/MagicSoft/Mars/mbase/MTaskList.cc (revision 7687) +++ trunk/MagicSoft/Mars/mbase/MTaskList.cc (revision 7688) @@ -338,4 +338,54 @@ // -------------------------------------------------------------------------- // +// Add all objects in list to the tasklist after task where. If some of +// them do not inherit from MTask return kFALSE, also if AddToListAfter +// returns an error for one of the tasks +// +Bool_t MTaskList::AddToListAfter(const TList &list, const MTask *where, const char *tType) +{ + TIter Next(&list); + TObject *obj=0; + while ((obj=Next())) + { + if (!obj->InheritsFrom(MTask::Class())) + { + *fLog << err << "ERROR - Object " << obj->GetName() << " doesn't inherit from MTask..." << endl; + return kFALSE; + } + + if (!AddToListAfter(static_cast(obj), where, tType)) + return kFALSE; + + where = static_cast(obj); + } + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Add all objects in list to the tasklist before task where. If some of +// them do not inherit from MTask return kFALSE, also if AddToListBefore +// returns an error for one of the tasks +// +Bool_t MTaskList::AddToListBefore(const TList &list, const MTask *where, const char *tType) +{ + TIter Next(&list); + TObject *obj=0; + while ((obj=Next())) + { + if (!obj->InheritsFrom(MTask::Class())) + { + *fLog << err << "ERROR - Object " << obj->GetName() << " doesn't inherit from MTask..." << endl; + return kFALSE; + } + + if (!AddToListBefore(static_cast(obj), where, tType)) + return kFALSE; + } + return kTRUE; +} + +// -------------------------------------------------------------------------- +// // Find an object in the list. // 'name' is the name of the object you are searching for. @@ -913,4 +963,30 @@ // -------------------------------------------------------------------------- // +// Removes all task of the TList from the tasklist. Returns kFALSE if any +// of the objects was not an MTask or not found in the list. +// +Bool_t MTaskList::RemoveFromList(const TList &list) +{ + Bool_t rc = kTRUE; + + TIter Next(&list); + TObject *obj=0; + while ((obj=Next())) + { + if (!obj->InheritsFrom(MTask::Class())) + { + *fLog << err << "ERROR - Object " << obj->GetName() << " doesn't inherit from MTask..." << endl; + rc = kFALSE; + continue; + } + + if (!RemoveFromList(static_cast(obj))) + rc = kFALSE; + } + return rc; +} + +// -------------------------------------------------------------------------- +// // Find an object with the same name in the list and replace it with // the new one. If the kIsOwner flag is set and the object was not Index: trunk/MagicSoft/Mars/mbase/MTaskList.h =================================================================== --- trunk/MagicSoft/Mars/mbase/MTaskList.h (revision 7687) +++ trunk/MagicSoft/Mars/mbase/MTaskList.h (revision 7688) @@ -48,4 +48,6 @@ Bool_t AddToList(MTask *task, const char *tType="All"); Bool_t AddToList(const TList &list, const char *tType="All"); + Bool_t AddToListAfter(const TList &list, const MTask *where, const char *tType="All"); + Bool_t AddToListBefore(const TList &list, const MTask *where, const char *tType="All"); void SetSerialNumber(Byte_t num); @@ -53,4 +55,5 @@ Bool_t Replace(MTask *obj); Bool_t RemoveFromList(MTask *task); + Bool_t RemoveFromList(const TList &list); TObject *FindObject(const char *name) const; Index: trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc (revision 7687) +++ trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc (revision 7688) @@ -19,5 +19,5 @@ ! Author(s): Marcos Lopez 10/2003 ! -! Copyright: MAGIC Software Development, 2000-2005 +! Copyright: MAGIC Software Development, 2000-2006 ! ! @@ -26,5 +26,5 @@ ////////////////////////////////////////////////////////////////////////////// // -// MMcWeightEnergySlopeCalc +// MMcSpectrumWeight // // Change the spectrum of the MC showers simulated with Corsika (a power law) @@ -104,14 +104,15 @@ AddToBranchList("MMcEvt.fEnergy"); - fNameWeight = "MWeight"; - fNameMcEvt = "MMcEvt"; - - fNewSlope = -1; - fOldSlope = -1; - - fEnergyMin = -1; - fEnergyMax = -2; - - fNorm = 1; + fNameWeight = "MWeight"; + fNameMcEvt = "MMcEvt"; + + fNewSlope = -1; + fOldSlope = -1; + + fEnergyMin = -1; + fEnergyMax = -2; + + fNorm = 1; + fNormEnergy = 1; fAllowChange = kFALSE; @@ -235,8 +236,14 @@ // Return the formula to calculate weights. // Is is compiled by -// o = integral(fEnergyMin, fEnergyMax, GetFormulaSpecOldX()); -// n = integral(fEnergyMin, fEnergyMax, GetFormulaSpecNewX()); -// -// result: fNorm*o/n*GetFormulaNewSpec()/GetFormulaOldSpec() +// o1 = integral(fEnergyMin, fEnergyMax, GetFormulaSpecOldX()); +// n1 = integral(fEnergyMin, fEnergyMax, GetFormulaSpecNewX()); +// o2 = CalcSpecOld(fNormEnergy); +// n2 = CalcSpecNew(fNormEnergy); +// +// result (fNormEnergy<0): +// fNorm*o1/n1*GetFormulaNewSpec()/GetFormulaOldSpec() +// +// result (fNormEnergy>0): +// fNorm*o2/n2*GetFormulaNewSpec()/GetFormulaOldSpec() // // fNorm is 1 by default but can be overwritten using SetNorm() @@ -254,6 +261,6 @@ return Form("%.16f", fNorm); - const Double_t iold = GetSpecOldIntegral(); - const Double_t inew = GetSpecNewIntegral(); + const Double_t iold = fNormEnergy<0 ? GetSpecOldIntegral() : CalcSpecOld(fNormEnergy); + const Double_t inew = fNormEnergy<0 ? GetSpecNewIntegral() : CalcSpecNew(fNormEnergy); const Double_t norm = fNorm*iold/inew; @@ -282,4 +289,26 @@ TF1 funcold("Dummy", GetFormulaSpecOldX()); return funcold.Integral(fEnergyMin, fEnergyMax); +} + +// --------------------------------------------------------------------------- +// +// Returns the value of GetFormulaSpecNewX() at the energy e describing +// the destination spectrum +// +Double_t MMcSpectrumWeight::CalcSpecNew(Double_t e) const +{ + TF1 funcnew("Dummy", GetFormulaSpecNewX()); + return funcnew.Eval(e); +} + +// --------------------------------------------------------------------------- +// +// Returns the value of GetFormulaSpecOldX() at the energy e describing +// the simulated spectrum +// +Double_t MMcSpectrumWeight::CalcSpecOld(Double_t e) const +{ + TF1 funcnew("Dummy", GetFormulaSpecOldX()); + return funcnew.Eval(e); } @@ -392,4 +421,5 @@ *fLog << " New spectral slope: " << fNewSlope << endl; *fLog << " User normalization: " << fNorm << endl; + *fLog << " Spectra are normalized: " << (fNormEnergy<0?"by integral":Form("at %.1fGeV", fNormEnergy)) << endl; *fLog << " Old Spectrum: " << GetFormulaSpecOldX() << " (I=" << GetSpecOldIntegral() << ")" << endl; *fLog << " New Spectrum: " << GetFormulaSpecNewX() << " (I=" << GetSpecNewIntegral() << ")" << endl; @@ -446,7 +476,8 @@ // // Read the setup from a TEnv, eg: -// MMcSpectrumWeight.NewSlope: -2.6 -// MMcSpectrumWeight.Norm: 1.0 -// MMcSpectrumWeight.Formula: pow(X, -2.6) +// MMcSpectrumWeight.NewSlope: -2.6 +// MMcSpectrumWeight.Norm: 1.0 +// MMcSpectrumWeight.NormEnergy: 200 +// MMcSpectrumWeight.Formula: pow(X, -2.6) // Int_t MMcSpectrumWeight::ReadEnv(const TEnv &env, TString prefix, Bool_t print) @@ -463,4 +494,9 @@ SetNorm(GetEnvValue(env, prefix, "Norm", fNorm)); } + if (IsEnvDefined(env, prefix, "NormEnergy", print)) + { + rc = kTRUE; + SetNormEnergy(GetEnvValue(env, prefix, "NormEnergy", fNormEnergy)); + } if (IsEnvDefined(env, prefix, "Formula", print)) { Index: trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.h (revision 7687) +++ trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.h (revision 7688) @@ -40,4 +40,5 @@ Double_t fNorm; // Normalization constant (additional normalization constant) + Double_t fNormEnergy; // Energy at which the spectra are normalized (default -1 means the integral is used) TString fFormula; // Text Formula for new spectrum: eg. "pow(MMcEvt.fEnergy, -2.0)" @@ -66,4 +67,5 @@ void SetNewSlope(Double_t s=-1) { fNewSlope = s; } void SetNorm(Double_t s=1) { fNorm = s; } + void SetNormEnergy(Double_t s=1) { fNormEnergy = s; } void SetFormula(const char *f="") { fFormula = f; } void SetEnergyRange(Double_t min=-2, Double_t max=-1) { fEnergyMin=min; fEnergyMax=max; } @@ -71,4 +73,5 @@ void SetWeightsZd(TH1 *h=0) { fWeightsZd = h; } void SetWeightsSize(TH1D *h=0); + Bool_t Set(const MMcCorsikaRunHeader &h); @@ -85,4 +88,7 @@ Double_t GetSpecOldIntegral() const; + Double_t CalcSpecNew(Double_t e) const; + Double_t CalcSpecOld(Double_t e) const; + Double_t GetEnergyMin() const { return fEnergyMin; } Double_t GetEnergyMax() const { return fEnergyMax; }