Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 7093) +++ trunk/MagicSoft/Mars/Changelog (revision 7094) @@ -29,4 +29,56 @@ - implement usage of extractor resolution (only if UseExtractorRes is set). + + + 2005/05/27 Thomas Bretz + + * Makefile: + - removed mmontecarlo directory + + * mmontecarlo/MMcEnergyEst.[h,cc], + mmontecarlo/MMcTimeGenerate.[h,cc], + mmontecarlo/MMcWeightEnergySpecCalc.[h,cc]: + - removed + + * sponde.rc: + - added new line for weighted spectral index + + * mbadpixels/MBadPixelsCalc.[h,cc]: + - added an option to perform the checks also in PostProcess + + * mhbase/MFillH.h: + - added default argument to SetWeight + + * mhbase/MH3.h: + - added Sumw2() member function + + * mhflux/MHCollectionArea.[h,cc]: + - added TLatex output to plots + - added some Getter + + * mjobs/MJSpectrum.[h,cc]: + - implemented the possibility to weight the monte carlo spectrum + to a new index or function. More details can be found + in MMcSpectrumWeight + - slightly changed the plot comparing the size distributions + - scale the comparsison plots by the resulting spectrum + + * mjobs/MJob.[h,cc]: + - added a member function to check ReadEnv of a single + container + + * mjobs/Makefile: + - added -I../mmc + + * mmc/MMcEvt.[hxx,cxx], mmc/McEvtBasic.[hxx,cxx]: + - changed the inheritance: MMcEvt now derives from MMcEvtBasic + so that both classes are interchangable + - increased both class versions + - chaged the default partictle in MMcEvtBasic from + kGAMMA to kUNDEFINED + - added new particle type: kUNDEFINED + + * mhflux/MMcSpectrumWeight.[h,cc]: + - added Index: trunk/MagicSoft/Mars/Makefile =================================================================== --- trunk/MagicSoft/Mars/Makefile (revision 7093) +++ trunk/MagicSoft/Mars/Makefile (revision 7094) @@ -56,5 +56,4 @@ msql \ mimage \ - mmontecarlo \ mhft \ mmc \ Index: trunk/MagicSoft/Mars/NEWS =================================================================== --- trunk/MagicSoft/Mars/NEWS (revision 7093) +++ trunk/MagicSoft/Mars/NEWS (revision 7094) @@ -2,4 +2,17 @@ *** Version + + - general: MMcEvt now derived from MMcEvtBasic which should + have no influence on compatibility with older camera files + + - sponde: The input MC spectrum can now be weighted to fake a + different spectrum. This is done via MMcSpectrumWeight. For + more details see the class description and sponde.rc + + - sponde: The paremetr comparsion plots are not scaled by + their entries anymore. Instead the MC plot is scaled by using + the result spectrum of the analysis. If the input MC spectrum + and the result spectrum has different slopes the absolut + normalization is normally wrong. Index: trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.cc (revision 7093) +++ trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.cc (revision 7094) @@ -86,5 +86,6 @@ // MBadPixelsCalc::MBadPixelsCalc(const char *name, const char *title) - : fPedestalLevel(3), fPedestalLevelVariance(5), fNamePedPhotCam("MPedPhotCam") + : fPedestalLevel(3), fPedestalLevelVariance(5), fNamePedPhotCam("MPedPhotCam"), + fCheckInProcess(kTRUE), fCheckInPostProcess(kFALSE) { fName = name ? name : gsDefName.Data(); @@ -129,42 +130,13 @@ // -------------------------------------------------------------------------- // -// Check the pedestal RMS of every pixel with respect to the mean pedestal RMS -// of the camera (Sigmabar). -// -// The pixels can be set as blind if the pedestalRMS is too big or 0. -// -// If you don't want to use this option set the PedestalLevel<=0; -// -// MBadPixelsCalc calc; -// calc.SetPedestalLevel(-1); -/* -void MBadPixelsCalc::CheckPedestalRMS() const -{ - const Int_t entries = fPedPhotCam->GetSize(); - - const Float_t meanPedRMS = fSigmabar->GetSigmabar(); - - for (Int_t i=0; iGetPixRatio(i); - const Double_t pixPedRms = (*fPedPhotCam)[i].GetRms(); - - if (pixPedRms*nratio > fPedestalLevel * meanPedRMS || pixPedRms == 0) - (*fBadPixels)[i].SetUnsuitable(MBadPixelsPix::kUnsuitableEvt); - } -} -*/ - -// -------------------------------------------------------------------------- -// // Check the pedestal Rms of the pixels: compute with 2 iterations the mean // for inner and outer pixels. Set as blind the pixels with too small or // too high pedestal Rms with respect to the mean. // -Bool_t MBadPixelsCalc::CheckPedestalRms() const -{ +Bool_t MBadPixelsCalc::CheckPedestalRms(MBadPixelsPix::UnsuitableType_t type) const +{ + if (fPedestalLevel<=0 && fPedestalLevelVariance<=0) + return kTRUE; + const Int_t entries = fPedPhotCam->GetSize(); @@ -266,5 +238,8 @@ continue; - (*fBadPixels)[i].SetUnsuitable(MBadPixelsPix::kUnsuitableEvt); + (*fBadPixels)[i].SetUnsuitable(type); + if (type==MBadPixelsPix::kUnsuitableRun) + (*fBadPixels)[i].SetUncalibrated(MBadPixelsPix::kDeadPedestalRms); + bads++; } @@ -281,5 +256,4 @@ } - // -------------------------------------------------------------------------- // @@ -287,8 +261,13 @@ Int_t MBadPixelsCalc::Process() { - if (fPedestalLevel>0 || fPedestalLevelVariance>0) - CheckPedestalRms(); - - return kTRUE; + return fCheckInProcess ? CheckPedestalRms(MBadPixelsPix::kUnsuitableEvt) : kTRUE; +} + +// -------------------------------------------------------------------------- +// +// +Int_t MBadPixelsCalc::PostProcess() +{ + return fCheckInPostProcess ? CheckPedestalRms(MBadPixelsPix::kUnsuitableRun) : kTRUE; } Index: trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.h =================================================================== --- trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.h (revision 7093) +++ trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.h (revision 7094) @@ -4,4 +4,7 @@ #ifndef MARS_MTask #include "MTask.h" +#endif +#ifndef MARS_MBadPixelsPix +#include "MBadPixelsPix.h" #endif @@ -21,10 +24,17 @@ TString fNamePedPhotCam; // name of the 'MPedPhotCam' container - - // void CheckPedestalRMS() const; - Bool_t CheckPedestalRms() const; + Bool_t fCheckInProcess; + Bool_t fCheckInPostProcess; + + // MBadPixelsCalc + Bool_t CheckPedestalRms(MBadPixelsPix::UnsuitableType_t t) const; + + // MTask Int_t PreProcess(MParList *pList); Int_t Process(); + Int_t PostProcess(); + + // MParContainer Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); @@ -32,7 +42,11 @@ MBadPixelsCalc(const char *name=NULL, const char *title=NULL); + // Setter void SetPedestalLevel(Float_t f) { fPedestalLevel=f; } void SetPedestalLevelVariance(Float_t f) { fPedestalLevelVariance=f; } void SetNamePedPhotCam(const char *name) { fNamePedPhotCam = name; } + + void EnableCheckInProcess(Bool_t b=kTRUE) { fCheckInProcess = b; } + void EnableCheckInPostProcess(Bool_t b=kTRUE) { fCheckInPostProcess = b; } ClassDef(MBadPixelsCalc, 1) // Task to find bad pixels (star, broken pixels, etc) Index: trunk/MagicSoft/Mars/mhbase/MFillH.h =================================================================== --- trunk/MagicSoft/Mars/mhbase/MFillH.h (revision 7093) +++ trunk/MagicSoft/Mars/mhbase/MFillH.h (revision 7094) @@ -63,5 +63,5 @@ void SetWeight(MParameterD *w) { fWeight = w; } - void SetWeight(const char *name) { fWeightName = name; } + void SetWeight(const char *name="MWeight") { fWeightName = name; } void SetDrawOption(Option_t *option=""); Index: trunk/MagicSoft/Mars/mhbase/MH3.h =================================================================== --- trunk/MagicSoft/Mars/mhbase/MH3.h (revision 7093) +++ trunk/MagicSoft/Mars/mhbase/MH3.h (revision 7094) @@ -26,7 +26,4 @@ Double_t fScale[3]; // Scale for the three axis (eg unit) - // TString fNameProfX; //! This should make sure, that gROOT doen't confuse the profile with something else - // TString fNameProfY; //! This should make sure, that gROOT doen't confuse the profile with something else - void StreamPrimitive(ofstream &out) const; @@ -36,7 +33,4 @@ kIsLogz = BIT(19) }; - -// Int_t DistancetoPrimitive(Int_t px, Int_t py); -// void ExecuteEvent(Int_t event, Int_t px, Int_t py); public: @@ -48,4 +42,5 @@ ~MH3(); + // Setter void SetScaleX(Double_t scale) { fScale[0] = scale; } void SetScaleY(Double_t scale) { fScale[1] = scale; } @@ -56,19 +51,20 @@ void SetLogz(Bool_t b=kTRUE) { b ? fHist->SetBit(kIsLogz) : fHist->ResetBit(kIsLogz); } + void Sumw2() const { if (fHist) fHist->Sumw2(); } + + // Getter Int_t GetDimension() const { return fDimension; } Int_t GetNbins() const; Int_t FindFixBin(Double_t x, Double_t y=0, Double_t z=0) const; - void SetName(const char *name); - void SetTitle(const char *title); - - Bool_t SetupFill(const MParList *pList); - Bool_t Fill(const MParContainer *par, const Stat_t w=1); + TH1 &GetHist() { return *fHist; } + const TH1 &GetHist() const { return *fHist; } TString GetDataMember() const; TString GetRule(const Char_t axis='x') const; - TH1 &GetHist() { return *fHist; } - const TH1 &GetHist() const { return *fHist; } + // MH + Bool_t SetupFill(const MParList *pList); + Bool_t Fill(const MParContainer *par, const Stat_t w=1); TH1 *GetHistByName(const TString name="") const { return fHist; } @@ -79,9 +75,14 @@ } + // MParContainer + MParContainer *New() const; + + // TObject + void SetName(const char *name); + void SetTitle(const char *title); + void SetColors() const; void Draw(Option_t *opt=NULL); void Paint(Option_t *opt=""); - - MParContainer *New() const; ClassDef(MH3, 1) // Generalized 1/2/3D-histogram for Mars variables Index: trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc (revision 7093) +++ trunk/MagicSoft/Mars/mhflux/MHCollectionArea.cc (revision 7094) @@ -31,4 +31,5 @@ #include "MHCollectionArea.h" +#include #include #include @@ -81,12 +82,12 @@ fHistAll.SetTitle("Number of events produced"); fHistAll.SetXTitle("\\Theta [deg]"); - fHistAll.SetYTitle("E [GeV]"); + fHistAll.SetYTitle("E_{mc} [GeV]"); fHistAll.SetDirectory(NULL); fHistAll.UseCurrentStyle(); fHEnergy.SetName("CollEnergy"); - fHEnergy.SetTitle("Collection Area vs Energy"); + fHEnergy.SetTitle("Collection Area"); fHEnergy.SetXTitle("E [GeV]"); - fHEnergy.SetYTitle("A [m^{2}]"); + fHEnergy.SetYTitle("A_{eff} [m^{2}]"); fHEnergy.SetDirectory(NULL); fHEnergy.UseCurrentStyle(); @@ -119,9 +120,5 @@ // read from run header, but it is not yet in!! // - const Float_t r1 = 0; - const Float_t r2 = fMcAreaRadius; - - //*fLog << warn << endl << dbginf << "WARNING! I will assume a maximum impact parameter of " << fMCAreaRadius << " meters for the MC events. Check that this is the true one!" <(gPad->FindObject("Efficiency")); + if (h) + { + const TString txt = Form("N/N_{0}=%.2f", + GetCollectionAreaEff(), + GetCollectionAreaAbs(), fMcAreaRadius); + + TLatex text(0.31, 0.95, txt); + text.SetBit(TLatex::kTextNDC); + text.SetTextSize(0.04); + text.Paint();*/ + return; + } + if (TString(option)=="paint4") + { + const TString txt = Form("A_{eff}=%.0fm^{2} A_{abs}=%.0fm^{2} r=%.0fm", + GetCollectionAreaEff(), + GetCollectionAreaAbs(), fMcAreaRadius); + + TLatex text(0.31, 0.95, txt); + text.SetBit(TLatex::kTextNDC); + text.SetTextSize(0.04); + text.Paint(); + return; + } + TVirtualPad *pad; @@ -405,4 +430,5 @@ h->SetNameTitle("Efficiency", "Combined cut and trigger efficiency"); h->SetDirectory(NULL); + AppendPad("paint3"); } else @@ -416,4 +442,5 @@ gPad->SetGridy(); fHEnergy.Draw(); + AppendPad("paint4"); } else @@ -428,5 +455,5 @@ //*fLog << energy << " " << theta << endl; - fHistSel.Fill(theta, energy); + fHistSel.Fill(theta, energy, weight); return kTRUE; Index: trunk/MagicSoft/Mars/mhflux/MHCollectionArea.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MHCollectionArea.h (revision 7093) +++ trunk/MagicSoft/Mars/mhflux/MHCollectionArea.h (revision 7094) @@ -56,4 +56,8 @@ const TH1D &GetHEnergy() const { return fHEnergy; } + Double_t GetEntries() const { return fHistAll.Integral(); } + Double_t GetCollectionAreaEff() const { return fHEnergy.Integral(); } + Double_t GetCollectionAreaAbs() const { return TMath::Pi()*fMcAreaRadius*fMcAreaRadius; } + /* void DrawAll(Option_t *option=""); Index: trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc =================================================================== --- trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc (revision 7094) +++ trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.cc (revision 7094) @@ -0,0 +1,396 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz 5/2005 +! Author(s): Marcos Lopez 10/2003 +! +! Copyright: MAGIC Software Development, 2000-2005 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MMcWeightEnergySlopeCalc +// +// Change the spectrum of the MC showers simulated with Corsika (a power law) +// to a new one, which can be either, again a power law but with a different +// spectral index, or a generalizeed spectrum. The new spectrum can be +// pass to this class in different ways: +// +// 1. Is the new spectrum will be a power law, just introduce the slope +// of this power law. +// 2. Is the new spectrum will have a general shape: +// The new spectrum is passed as a char* (SetFormula()) +// +// Method: +// ------- +// +// - Corsika spectrun: dN/dE = A * E^(a) +// with a = fOldSlope, and A = N/integral{E*de} from ELowLim to EUppLim +// +// - New spectrum: dN/dE = B * g(E) +// where B = N/integral{g*dE} from ELowLim to EUppLim, and N=NumEvents +// +// For converting the spectrum simulated with Corsika to the new one, we +// apply a weight to each event, given by: +// +// W(E) = B/A * g(E)/E^(a) +// +// In the case the new spectrum is simply a power law: dN/dE = B * E^(b), we +// have: +// +// W(E) = B/A * E^(b-a) +// +// (The factor B/A is used in order both the original and new spectrum have +// the same area (i.e. in order they represent the same number of showers)) +// +// +// Input Containers: +// MMcEvt +// MMcCorsikaRunHeader +// +// Output Container: +// MWeight [MParameterD] +// +////////////////////////////////////////////////////////////////////////////// +#include "MMcSpectrumWeight.h" + +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" +#include "MParameters.h" + +#include "MMcEvt.hxx" +#include "MMcCorsikaRunHeader.h" + +ClassImp(MMcSpectrumWeight); + +using namespace std; + +void MMcSpectrumWeight::Init(const char *name, const char *title) +{ + fName = name ? name : "MMcSpectrumWeight"; + fTitle = title ? title : "Task to calculate weights to change the energy spectrum"; + + AddToBranchList("MMcEvt.fEnergy"); + + fNameWeight = "MWeight"; + fNameMcEvt = "MMcEvt"; + + fNewSlope = -1; + fOldSlope = -1; + + fEnergyMin = -1; + fEnergyMax = -2; + + fNorm = 1; + + fAllowChange = kFALSE; + + fFunc = NULL; + fMcEvt = NULL; + fWeight = NULL; +} + +// --------------------------------------------------------------------------- +// +// Default Constructor. +// +MMcSpectrumWeight::MMcSpectrumWeight(const char *name, const char *title) +{ + Init(name,title); +} + +// --------------------------------------------------------------------------- +// +// Destructor. If necessary delete fFunc +// +MMcSpectrumWeight::~MMcSpectrumWeight() +{ + if (fFunc) + delete fFunc; +} + +// --------------------------------------------------------------------------- +// +// Search for +// - fNameMcEvt [MMcEvtBasic] +// +// Find/Create: +// - fNameWeight [MWeight] +// +Int_t MMcSpectrumWeight::PreProcess(MParList *pList) +{ + fMcEvt = (MMcEvt*)pList->FindObject(fNameMcEvt, "MMcEvtBasic"); + if (!fMcEvt) + { + *fLog << err << fNameMcEvt << " [MMcEvtBasic] not found... abort." << endl; + return kFALSE; + } + + fWeight = (MParameterD*)pList->FindCreateObj("MParameterD", fNameWeight); + if (!fWeight) + return kFALSE; + + return kTRUE; +} + +// --------------------------------------------------------------------------- +// +// Replace {fNameMcEvt}.fEnergy by "(x)" and return the result. +// +TString MMcSpectrumWeight::ReplaceX(TString str) const +{ + return str.ReplaceAll(Form("%s.fEnergy", fNameMcEvt.Data()), "(x)"); +} + +// --------------------------------------------------------------------------- +// +// Return the function corresponding to the mc spectrum with +// slope fOldSlope: pow({fNameMcEvt}.fEnergy, fOldSlope) +// +// The slope is returned as %.3f +// +TString MMcSpectrumWeight::GetFormulaSpecOld() const +{ + return Form("pow(%s.fEnergy, %.3f)", fNameMcEvt.Data(), fOldSlope); +} + +// --------------------------------------------------------------------------- +// +// Return the function corresponding to the new spectrum with +// slope fNewSlope: pow({fNameMcEvt}.fEnergy, fNewSlope) +// +// The slope is returned as %.3f +// +// If a different formula is set (SetFormula()) this formula is returned +// unchanged. +// +TString MMcSpectrumWeight::GetFormulaSpecNew() const +{ + return fFormula.IsNull() ? Form("pow(%s.fEnergy, %.3f)", fNameMcEvt.Data(), fNewSlope) : fFormula; +} + +// --------------------------------------------------------------------------- +// +// 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() +// +// fNorm is 1 by default but can be overwritten using SetNorm() +// +// If the formulas GetFormulaSpecOldX() and GetFormulaSpecNewX() +// are equal only fNorm is returned. +// +// The normalization constant is returned as %.16f +// +// Example: 0.3712780019*(pow(MMcEvt.fEnergy,-2.270))/(pow(MMcEvt.fEnergy,-2.600)) +// +TString MMcSpectrumWeight::GetFormulaWeights() const +{ + if (GetFormulaSpecOld()==GetFormulaSpecNew()) + return Form("%.16f", fNorm); + + const Double_t iold = GetSpecOldIntegral(); + const Double_t inew = GetSpecNewIntegral(); + + const Double_t norm = fNorm*iold/inew; + + return Form("%.16f*(%s)/(%s)", norm, GetFormulaSpecNew().Data(), GetFormulaSpecOld().Data()); +} + +// --------------------------------------------------------------------------- +// +// Returns the integral between fEnergyMin and fEnergyMax of +// GetFormulaSpecNewX() describing the destination spectrum +// +Double_t MMcSpectrumWeight::GetSpecNewIntegral() const +{ + TF1 funcnew("Dummy", GetFormulaSpecNewX()); + return funcnew.Integral(fEnergyMin, fEnergyMax); +} + +// --------------------------------------------------------------------------- +// +// Returns the integral between fEnergyMin and fEnergyMax of +// GetFormulaSpecOldX() describing the simulated spectrum +// +Double_t MMcSpectrumWeight::GetSpecOldIntegral() const +{ + TF1 funcold("Dummy", GetFormulaSpecOldX()); + return funcold.Integral(fEnergyMin, fEnergyMax); +} + +// --------------------------------------------------------------------------- +// +// Initialize fEnergyMin, fEnergymax and fOldSlope from MMcCorsikaRunHeader +// by GetELowLim(), GetEUppLim() and GetSlopeSpec(). +// +// If fEnergyMax>fEnergyMin (means: the values have already been +// initialized) and !fAllowChange the consistency of the new values +// with the present values is checked with a numerical precision of 1e-10. +// If one doesn't match kFALSE is returned. +// +// If the mc slope is -1 kFALSE is returned. +// +// If the new slope for the spectrum is -1 it is set to the original MC +// slope. +// +// fFunc is set to the formula returned by GetFormulaWeightsX() +// +Bool_t MMcSpectrumWeight::Set(const MMcCorsikaRunHeader &rh) +{ + if (fEnergyMax>fEnergyMin && !fAllowChange) + { + if (TMath::Abs(fOldSlope-rh.GetSlopeSpec())>1e-10) + { + *fLog << err; + *fLog << "ERROR - The slope of the Monte Carlo is not allowed to change "; + *fLog << "(" << fOldSlope << " --> " << rh.GetSlopeSpec() << ")... abort." << endl; + return kFALSE; + } + if (TMath::Abs(fEnergyMin-rh.GetELowLim())>1e-10) + { + *fLog << err; + *fLog << "ERROR - The minimum simulated Monte Carlo energy is not allowed to change "; + *fLog << "(" << fEnergyMin << " --> " << rh.GetELowLim() << ")... abort." << endl; + return kFALSE; + } + if (TMath::Abs(fEnergyMax-rh.GetEUppLim())>1e-10) + { + *fLog << err; + *fLog << "ERROR - The maximum simulated Monte Carlo energy is not allowed to change ("; + *fLog << "(" << fEnergyMax << " --> " << rh.GetEUppLim() << ")... abort." << endl; + return kFALSE; + } + return kTRUE; + } + + // + // Sanity checks to be sure that we won't divide by zero later on + // + if (rh.GetSlopeSpec()==-1) + { + *fLog << err << "The MC Slope of the power law must be different of -1... exit" << endl; + return kFALSE; + } + + fOldSlope = rh.GetSlopeSpec(); + fEnergyMin = rh.GetELowLim(); + fEnergyMax = rh.GetEUppLim(); + + if (fNewSlope==-1) + { + *fLog << inf << "The new slope of the power law is -1... no weighting applied." << endl; + fNewSlope = fOldSlope; + } + + TString form(GetFormulaWeightsX()); + + if (fFunc) + delete fFunc; + + fFunc = new TF1("", form); + gROOT->GetListOfFunctions()->Remove(fFunc); + fFunc->SetName("SpectralWeighs"); + + return kTRUE; +} + +// --------------------------------------------------------------------------- +// +// The current contants are printed +// +void MMcSpectrumWeight::Print(Option_t *o) const +{ + *fLog << all << GetDescriptor() << endl; + *fLog << " Simulated energy range: " << fEnergyMin << "GeV - " << fEnergyMax << "GeV" << endl; + *fLog << " Simulated spectral slope: " << fOldSlope << endl; + *fLog << " New slope: " << fNewSlope << endl; + *fLog << " User normalization: " << fNorm << endl; + *fLog << " Old Spectrum: " << GetFormulaSpecOldX() << " (I=" << GetSpecOldIntegral() << ")" << endl; + *fLog << " New Spectrum: " << GetFormulaSpecNewX() << " (I=" << GetSpecNewIntegral() << ")" << endl; + if (fFunc) + *fLog << " Weight function: " << fFunc->GetTitle() << endl; +} + +// ---------------------------------------------------------------------------- +// +// Executed each time a new root file is loaded +// We will need fOldSlope and fE{Upp,Low}Lim to calculate the weights +// +Bool_t MMcSpectrumWeight::ReInit(MParList *plist) +{ + MMcCorsikaRunHeader *rh = (MMcCorsikaRunHeader*)plist->FindObject("MMcCorsikaRunHeader"); + if (!rh) + { + *fLog << err << "MMcCorsikaRunHeader not found... abort." << endl; + return kFALSE; + } + + return Set(*rh); +} + +// ---------------------------------------------------------------------------- +// +// Fill the result of the evaluation of fFunc at fEvEvt->GetEnergy +// into the weights container. +// +Int_t MMcSpectrumWeight::Process() +{ + const Double_t e = fMcEvt->GetEnergy(); + + fWeight->SetVal(fFunc->Eval(e)); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Read the setup from a TEnv, eg: +// MMcSpectrumWeight.NewSlope: -2.6 +// MMcSpectrumWeight.Norm: 1.0 +// MMcSpectrumWeight.Formula: pow(MMcEvt.fEnergy, -2.6) +// +Int_t MMcSpectrumWeight::ReadEnv(const TEnv &env, TString prefix, Bool_t print) +{ + Bool_t rc = kFALSE; + if (IsEnvDefined(env, prefix, "NewSlope", print)) + { + rc = kTRUE; + SetNewSlope(GetEnvValue(env, prefix, "NewSlope", fNewSlope)); + } + if (IsEnvDefined(env, prefix, "Norm", print)) + { + rc = kTRUE; + SetNorm(GetEnvValue(env, prefix, "Norm", fNorm)); + } + if (IsEnvDefined(env, prefix, "Formula", print)) + { + rc = kTRUE; + SetFormula(GetEnvValue(env, prefix, "Formula", fFormula)); + } + + return rc; +} Index: trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.h =================================================================== --- trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.h (revision 7094) +++ trunk/MagicSoft/Mars/mhflux/MMcSpectrumWeight.h (revision 7094) @@ -0,0 +1,84 @@ +#ifndef MARS_MMcSpectrumWeight +#define MARS_MMcSpectrumWeight + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class TF1; +class MParList; +class MMcEvt; +class MParameterD; +class MMcCorsikaRunHeader; + +class MMcSpectrumWeight : public MTask +{ +private: + const MMcEvt *fMcEvt; // Pointer to the container with the MC energy + MParameterD *fWeight; // Pointer to the output MWeight container + + TString fNameWeight; // Name of the MWeight container + TString fNameMcEvt; // Name of the MMcEvt container + + TF1 *fFunc; // Function calculating the weights + + Double_t fOldSlope; // Slope of energy spectrum generated with Corsika + Double_t fNewSlope; // Slope of the new spectrum (if it is a power law) + + Double_t fEnergyMin; // Minimum energy simulated + Double_t fEnergyMax; // Maximum energy simulated + + Double_t fNorm; // Normalization constant (additional normalization constant) + + TString fFormula; // Text Formula for new spectrum: eg. "pow(MMcEvt.fEnergy, -2.0)" + + Bool_t fAllowChange; + + // MMcSpectrumWeight + void Init(const char *name, const char *title); + TString ReplaceX(TString) const; + + // MTask + Bool_t ReInit(MParList *plist); + Int_t PreProcess(MParList *pList); + Int_t Process(); + + // MParContainer + Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); + +public: + MMcSpectrumWeight(const char *name=NULL, const char *title=NULL); + ~MMcSpectrumWeight(); + + // Setter + void SetNameWeight(const char *n="MWeight") { fNameWeight = n; } + void SetNameMcEvt(const char *n="MMcEvt") { fNameMcEvt = n; } + void SetNewSlope(Double_t s=-1) { fNewSlope = s; } + void SetNorm(Double_t s=1) { fNorm = s; } + void SetFormula(const char *f="") { fFormula = f; } + void SetEnergyRange(Double_t min=-2, Double_t max=-1) { fEnergyMin=min; fEnergyMax=max; } + void SetOldSlope(Double_t s=-2.6) { fOldSlope=s; } + Bool_t Set(const MMcCorsikaRunHeader &h); + + // Getter + TString GetFormulaSpecOld() const; + TString GetFormulaSpecNew() const; + TString GetFormulaWeights() const; + + TString GetFormulaSpecOldX() const { return ReplaceX(GetFormulaSpecOld()); } + TString GetFormulaSpecNewX() const { return ReplaceX(GetFormulaSpecNew()); } + TString GetFormulaWeightsX() const { return ReplaceX(GetFormulaWeights()); } + + Double_t GetSpecNewIntegral() const; + Double_t GetSpecOldIntegral() const; + + Double_t GetEnergyMin() const { return fEnergyMin; } + Double_t GetEnergyMax() const { return fEnergyMax; } + + // TObject + void Print(Option_t *o="") const; + + ClassDef(MMcSpectrumWeight, 0) // Task to calculate weights to change the energy spectrum +}; + +#endif Index: trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (revision 7093) +++ trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (revision 7094) @@ -52,4 +52,5 @@ #include "MBinning.h" #include "MDataSet.h" +#include "MMcCorsikaRunHeader.h" // Spectrum @@ -58,4 +59,5 @@ #include "../mhflux/MHCollectionArea.h" #include "../mhflux/MHEnergyEst.h" +#include "../mhflux/MMcSpectrumWeight.h" // Eventloop @@ -155,4 +157,48 @@ } +// -------------------------------------------------------------------------- +// +// Read the first MMcCorsikaRunHeader from the RunHeaders tree in +// the dataset. +// The simulated energy range and spectral slope is initialized from +// there. +// In the following eventloops the forced check in MMcSpectrumWeight +// ensures, that the spectral slope and energy range doesn't change. +// +Bool_t MJSpectrum::InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const +{ + fLog->Separator("Initialize energy weighting"); + + MParList l; + l.AddToList(&w); + if (!CheckEnv(l)) + { + *fLog << err << "ERROR - Reading resources for MMcSpectrumWeight failed." << endl; + return kFALSE; + } + + TChain chain("RunHeaders"); + set.AddFilesOn(chain); + + MMcCorsikaRunHeader *h=0; + chain.SetBranchAddress("MMcCorsikaRunHeader.", &h); + chain.GetEntry(1); + + if (!h) + { + *fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from DataSet." << endl; + return kFALSE; + } + + if (!w.Set(*h)) + { + *fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl; + return kFALSE; + } + + w.Print(); + return kTRUE; +} + Float_t MJSpectrum::ReadInput(MParList &plist, TH1D &h1, TH1D &h2) { @@ -211,10 +257,15 @@ } -Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h) const +Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const { // Some debug output fLog->Separator("Compiling original MC distribution"); - *fLog << inf << "Please stand by, this may take a while..." << flush; + weight.SetNameMcEvt("MMcEvtBasic"); + const TString w(weight.GetFormulaWeights()); + weight.SetNameMcEvt(); + + *fLog << inf << "Using weights: " << w << endl; + *fLog << "Please stand by, this may take a while..." << flush; if (fDisplay) @@ -236,5 +287,5 @@ h.SetYTitle("E [GeV]"); h.SetZTitle("Counts"); - chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", "", "goff"); + chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", w, "goff"); } else @@ -243,5 +294,5 @@ h.SetXTitle("\\Theta [\\circ]"); h.SetYTitle("Counts"); - chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", "", "goff"); + chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", w, "goff"); } h.SetDirectory(0); @@ -413,5 +464,5 @@ } -Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set) const +Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &weight) const { MTaskList tlist1; @@ -439,4 +490,5 @@ } tlist1.AddToList(&readmc); + tlist1.AddToList(&weight); temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg"); @@ -477,8 +529,9 @@ bins3->SetName("BinningTheta"); } + return kTRUE; } -void MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const +TArrayD MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const { TH1D collarea(area.GetHEnergy()); @@ -546,5 +599,5 @@ f.SetParameter(1, 1.9e-6); f.SetLineColor(kGreen); - spectrum.Fit(&f, "NI", "", 55, 2e4); + spectrum.Fit(&f, "NIM", "", 55, 2e4); f.DrawCopy("same"); @@ -558,9 +611,13 @@ tex.DrawLatex(2e2, 7e-5, str); */ -} - -Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot) const -{ - cout << name << endl; + + TArrayD res(2); + res[0] = f.GetParameter(0); + res[1] = f.GetParameter(1); + return res; +} + +Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const +{ TString same(name); same += "Same"; @@ -587,13 +644,13 @@ const MAlphaFitter *fit = (MAlphaFitter*)plist.FindObject("MAlphaFitter"); - const Double_t scale = fit ? fit->GetScaleFactor() : 1; + const Double_t ascale = fit ? fit->GetScaleFactor() : 1; gPad->SetBorderMode(0); h2->SetLineColor(kBlack); h3->SetLineColor(kBlue); - h2->Add(h1, -scale); - - h2->Scale(1./h2->Integral()); - h3->Scale(1./h3->Integral()); + h2->Add(h1, -ascale); + + //h2->Scale(1./ontime); //h2->Integral()); + h3->Scale(scale); //h3->Integral()); h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum())); @@ -609,5 +666,5 @@ } -Bool_t MJSpectrum::DisplaySize(MParList &plist) const +Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) const { *fLog << inf << "Reading from file: " << fPathIn << endl; @@ -644,21 +701,33 @@ excess->SetTitle("Number of excess events vs Size (data, mc/blue)"); - excess->Scale(1./excess->Integral()); - excess = excess->DrawCopy(); + excess = excess->DrawCopy("E2"); // Don't do this on the original object! excess->SetStats(kFALSE); + excess->SetMarkerStyle(kFullDotMedium); + excess->SetFillColor(kBlack); + excess->SetFillStyle(0); + excess->SetName("Excess "); + excess->SetDirectory(0); TObject *o=0; - if ((o=plist.FindObject("ExcessSize"))) + if ((o=plist.FindObject("ExcessMC"))) { TH1 *histsel = (TH1F*)o->FindObject(""); if (histsel) { - histsel->Scale(1./histsel->Integral()); + if (scale<0) + scale = excess->Integral()/histsel->Integral(); + + histsel->Scale(scale); histsel->SetLineColor(kBlue); histsel->SetBit(kCanDelete); - histsel = histsel->DrawCopy("same"); + histsel = histsel->DrawCopy("E1 same"); // Don't do this on the original object! histsel->SetStats(kFALSE); + + fLog->Separator("Kolmogorov Test"); + histsel->KolmogorovTest(excess, "DX"); + fLog->Separator("Chi^2 Test"); + histsel->Chi2Test(excess, "P"); } } @@ -666,17 +735,17 @@ // -------------- Comparison of Image Parameters -------------- c.cd(2); - PlotSame(arr, plist, "Dist", "HilSrc", "MHHilSrcMCPost"); + PlotSame(arr, plist, "Dist", "HilSrc", "MHHilSrcMCPost", scale); c.cd(3); - PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost"); + PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost", scale); c.cd(4); - PlotSame(arr, plist, "M3l", "HilExt", "MHHilExtMCPost"); + PlotSame(arr, plist, "M3l", "HilExt", "MHHilExtMCPost", scale); c.cd(5); - PlotSame(arr, plist, "Conc1", "NewPar", "MHNewParMCPost"); + PlotSame(arr, plist, "Conc1", "NewPar", "MHNewParMCPost", scale); c.cd(6); - PlotSame(arr, plist, "Width", "PostCut", "MHHillasMCPost"); + PlotSame(arr, plist, "Width", "PostCut", "MHHillasMCPost", scale); return kTRUE; @@ -734,9 +803,13 @@ } + MMcSpectrumWeight weight; + if (!InitWeighting(set, weight)) + return kFALSE; + PrintSetup(fit); bins3.SetEdges(temp1, 'x'); TH1D temp2(temp1); - if (!ReadOrigMCDistribution(set, temp2)) + if (!ReadOrigMCDistribution(set, temp2, weight)) return kFALSE; @@ -754,5 +827,5 @@ hist.UseCurrentStyle(); MH::SetBinning(&hist, &bins3/*temp1.GetXaxis()*/, &bins2/*excess.GetXaxis()*/); - if (!ReadOrigMCDistribution(set, hist)) + if (!ReadOrigMCDistribution(set, hist, weight)) return kFALSE; @@ -762,9 +835,14 @@ for (int x=0; xSaveAsRoot(fPathOut); Index: trunk/MagicSoft/Mars/mjobs/MJSpectrum.h =================================================================== --- trunk/MagicSoft/Mars/mjobs/MJSpectrum.h (revision 7093) +++ trunk/MagicSoft/Mars/mjobs/MJSpectrum.h (revision 7094) @@ -18,4 +18,5 @@ class MStatusArray; class MHCollectionArea; +class MMcSpectrumWeight; class MJSpectrum : public MJob @@ -35,15 +36,16 @@ Bool_t ReadTask(MTask* &task, const char *name) const; Float_t ReadInput(MParList &plist, TH1D &h1, TH1D &size); - Bool_t ReadOrigMCDistribution(const MDataSet &set, TH1 &h) const; + Bool_t ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &w) const; Bool_t GetThetaDistribution(TH1D &temp1, TH1D &temp2) const; Bool_t Refill(MParList &plist, TH1D &h) const; + Bool_t InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const; // Display Output void PrintSetup(const MAlphaFitter &fit) const; void DisplayResult(const TH2D &mh1) const; - Bool_t IntermediateLoop(MParList &plist, MH3 &h1, TH1D &temp1, const MDataSet &set) const; - void DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const; - Bool_t DisplaySize(MParList &plist) const; - Bool_t PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot) const; + Bool_t IntermediateLoop(MParList &plist, MH3 &h1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &w) const; + TArrayD DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const; + Bool_t DisplaySize(MParList &plist, Double_t scale) const; + Bool_t PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const; public: Index: trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.cc =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.cc (revision 7093) +++ (revision ) @@ -1,350 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * This file is part of MARS, the MAGIC Analysis and Reconstruction -! * Software. It is distributed to you in the hope that it can be a useful -! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. -! * It is distributed WITHOUT ANY WARRANTY. -! * -! * Permission to use, copy, modify and distribute this software and its -! * documentation for any purpose is hereby granted without fee, -! * provided that the above copyright notice appear in all copies and -! * that both that copyright notice and this permission notice appear -! * in supporting documentation. It is provided "as is" without express -! * or implied warranty. -! * -! -! -! Author(s): Thomas Bretz 1/2002 -! Author(s): Wolfgang Wittek 1/2002 -! Author(s): Abelardo Moralejo 4/2003 -! Author(s): Marcos Lopez 5/2004 -! -! Copyright: MAGIC Software Development, 2000-2003 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// -// MMcEnergyEst -// -// Class for otimizing the parameters of the energy estimator -// -// FIXME: the class must be made more flexible, allowing for different -// parametrizations to be used. Also a new class is needed, a container -// for the parameters of the energy estimator. -// FIXME: the starting values of the parameters are now fixed. -// -///////////////////////////////////////////////////////////////////////////// -#include "MMcEnergyEst.h" - -#include // fabs on Alpha - -#include -#include -#include - -#include "MParList.h" -#include "MTaskList.h" -#include "MGeomCamCT1.h" -#include "MFEventSelector.h" -#include "MReadTree.h" -#include "MFCT1SelFinal.h" -#include "MHMatrix.h" -#include "MEnergyEstParam.h" -#include "MMatrixLoop.h" -#include "MChisqEval.h" -#include "MEvtLoop.h" -#include "MDataElement.h" -#include "MDataMember.h" -#include "MLog.h" -#include "MLogManip.h" -#include "MParameters.h" - -ClassImp(MMcEnergyEst); - -using namespace std; - -//------------------------------------------------------------------------ -// -// fcn calculates the function to be minimized (using TMinuit::Migrad) -// -static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) -{ - MEvtLoop *evtloop = (MEvtLoop*)gMinuit->GetObjectFit(); - - MParList *plist = evtloop->GetParList(); - MTaskList *tlist = (MTaskList*)plist->FindObject("MTaskList"); - - - // Pass current minuit parameters to the energy estimation class - MEnergyEstParam *eest = (MEnergyEstParam*)tlist->FindObject("MEnergyEstParam"); - eest->SetCoeff(TArrayD(eest->GetNumCoeff(), par)); - - - evtloop->Eventloop(); - - - // Get result of the ChiSquare - MParameterD *eval = (MParameterD*)plist->FindObject("MFitResult", "MParameterD"); - - f = eval->GetVal(); -} - - -// -------------------------------------------------------------------------- -// -// Default constructor. -// -MMcEnergyEst::MMcEnergyEst(const char *name, const char *title) -{ - fName = name ? name : "MMcEnergyEst"; - fTitle = title ? title : "Optimizer of the energy estimator"; - - SetHillasName("MHillas"); - SetHillasSrcName("MHillasSrc"); - - // - // Set initial values of the parameters (close enough to the final - // ones, taken from previous runs of the procedure). Parameter - // fA[4] is not used in the default energy estimation model (from - // D. Kranich). - // - fA.Set(5); - fB.Set(7); - - fA[0] = 21006.2; - fA[1] = -43.2648; - fA[2] = -690.403; - fA[3] = -0.0428544; - fA[4] = 1.; - fB[0] = -3391.05; - fB[1] = 136.58; - fB[2] = 0.253807; - fB[3] = 254.363; - fB[4] = 61.0386; - fB[5] = -0.0190984; - fB[6] = -421695; -} - -Bool_t MMcEnergyEst::SetCoeff(TArrayD &coeff) -{ - if (coeff.GetSize() != fA.GetSize()+fB.GetSize()) - { - *fLog << err << dbginf << "Wrong number of parameters!" << endl; - return(kFALSE); - } - - for (Int_t i = 0; i < fA.GetSize(); i++) - fA[i] = coeff[i]; - for (Int_t i = 0; i < fB.GetSize(); i++) - fB[i] = coeff[i+fA.GetSize()]; - - return(kTRUE); - -} - -//------------------------------------------------------------------------ -// -// Optimization (via migrad minimization) of parameters of energy estimation. -// -void MMcEnergyEst::FindParams() -{ - MParList parlist; - - MFEventSelector selector; - selector.SetNumSelectEvts(fNevents); - *fLog << inf << "Read events from file '" << fInFile << "'" << endl; - - MReadTree read("Events"); - read.AddFile(fInFile); - read.DisableAutoScheme(); - read.SetSelector(&selector); - - *fLog << inf << "Define columns of matrix" << endl; - MHMatrix matrix; - Int_t colenergy = matrix.AddColumn("MMcEvt.fEnergy"); - Int_t colimpact = matrix.AddColumn("MMcEvt.fImpact"); - - if (colenergy < 0 || colimpact < 0) - { - *fLog << err << dbginf << "colenergy, colimpact = " << colenergy << ", " - << colimpact << endl; - return; - } - - MEnergyEstParam eest(fHillasName); - eest.Add(fHillasSrcName); - eest.InitMapping(&matrix); - - *fLog << inf << "--------------------------------------" << endl; - *fLog << inf << "Fill events into the matrix" << endl; - if ( !matrix.Fill(&parlist, &read, fEventFilter) ) - return; - *fLog << inf << "Matrix was filled with " << matrix.GetNumRows() - << inf << " events" << endl; - - //----------------------------------------------------------------------- - // - // Setup the eventloop which will be executed in the fcn of MINUIT - // - *fLog << inf << "--------------------------------------" << endl; - *fLog << inf << "Setup event loop to be used in 'fcn'" << endl; - - MTaskList tasklist; - - MMatrixLoop loop(&matrix); - - MChisqEval eval; - eval.SetY1(new MDataElement(&matrix, colenergy)); - eval.SetY2(new MDataMember("MEnergyEst.fEnergy")); - eval.SetOwner(); - - // - // Entries in MParList - - parlist.AddToList(&tasklist); - - // - // Entries in MTaskList - - tasklist.AddToList(&loop); - tasklist.AddToList(&eest); - tasklist.AddToList(&eval); - - - *fLog << inf << "Event loop was setup" << endl; - MEvtLoop evtloop; - evtloop.SetParList(&parlist); - - // - //---------- Start of minimization part -------------------- - // - // Do the minimization with MINUIT - // - // Be careful: This is not thread safe - // - *fLog << inf << "--------------------------------------" << endl; - *fLog << inf << "Start minimization for the energy estimator" << endl; - - gMinuit = new TMinuit(12); - gMinuit->SetPrintLevel(-1); - - gMinuit->SetFCN(fcn); - gMinuit->SetObjectFit(&evtloop); - - // Ready for: gMinuit->mnexcm("SET ERR", arglist, 1, ierflg) - - if (gMinuit->SetErrorDef(1)) - { - *fLog << err << dbginf << "SetErrorDef failed." << endl; - return; - } - - // - // Set starting values and step sizes for parameters - // - for (Int_t i=0; iDefineParameter(i, name, vinit, step, limlo, limup); - if (!rc) - continue; - - *fLog << err << dbginf << "Error in defining parameter #" << i << endl; - return; - } - - for (Int_t i=0; iDefineParameter(i+fA.GetSize(), name, vinit, step, limlo, limup); - if (!rc) - continue; - - *fLog << err << dbginf << "Error in defining parameter #" << i+fA.GetSize() << endl; - return; - } - - TStopwatch clock; - clock.Start(); - - // Now ready for minimization step: - - gLog.SetNullOutput(kTRUE); - Bool_t rc = gMinuit->Migrad(); - gLog.SetNullOutput(kFALSE); - - if (rc) - { - *fLog << err << dbginf << "Migrad failed." << endl; - return; - } - - *fLog << inf << endl; - clock.Stop(); - clock.Print(); - *fLog << inf << endl; - - *fLog << inf << "Resulting Chisq: " << gMinuit->fAmin << endl; - - // - // Update values of fA, fB: - // - for (Int_t i = 0; i < fA.GetSize(); i++) - { - Double_t x1, x2; - gMinuit->GetParameter(i,x1,x2); - fA[i] = x1; - } - for (Int_t i = fA.GetSize(); i < fA.GetSize()+fB.GetSize(); i++) - { - Double_t x1, x2; - gMinuit->GetParameter(i,x1,x2); - fB[i-fA.GetSize()] = x1; - } - - // eest.Print(); - eest.StopMapping(); - *fLog << inf << "Minimization for the energy estimator finished" << endl; - -} - -//------------------------------------------------------------------------ -// -// Print current values of parameters -// -void MMcEnergyEst::Print(Option_t *o) const -{ - for (Int_t i = 0; i < fA.GetSize(); i++) - *fLog << inf << "Parameter " << i << ": " << const_cast(fA)[i] << endl; - - for (Int_t i = fA.GetSize(); i < fA.GetSize()+fB.GetSize(); i++) - *fLog << inf << "Parameter " << i << ": " << const_cast(fB)[i-fA.GetSize()] << endl; - - /* - // Print results - Double_t amin, edm, errdef; - Int_t nvpar, nparx, icstat; - gMinuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat); - gMinuit->mnprin(3, amin); - */ - -} Index: trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.h =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.h (revision 7093) +++ (revision ) @@ -1,58 +1,0 @@ -#ifndef MARS_MMcEnergyEst -#define MARS_MMcEnergyEst - -#ifndef MARS_MParContainer -#include "MParContainer.h" -#endif - -#ifndef ROOT_TArrayD -#include -#endif - -#include "MFilter.h" - -class MMcEnergyEst : public MParContainer -{ -private: - - TString fInFile, fOutFile; - TString fHillasName; - TString fHillasSrcName; - Int_t fNevents; - - MFilter *fEventFilter; //! - - TArrayD fA; - TArrayD fB; - -public: - MMcEnergyEst(const char *name=NULL, const char *title=NULL); - - void SetInFile(const TString &name) {fInFile = name;} - void SetOutFile(const TString &name) {fOutFile = name;} - void SetHillasName(const TString &name) {fHillasName = name;} - void SetHillasSrcName(const TString &name) {fHillasSrcName = name;} - void SetEventFilter(MFilter *filter) {fEventFilter = filter;} - void SetNevents(Int_t n) {fNevents = n;} - - TString GetInFile() const {return fInFile;} - TString GetOutFile() const {return fOutFile;} - TString GetHillasName() const {return fHillasName;} - TString GetHillasSrcName() const {return fHillasSrcName;} - Int_t GetNevents() const {return fNevents;} - - Int_t GetNumCoeffA() const {return fA.GetSize(); } - Int_t GetNumCoeffB() const {return fB.GetSize(); } - - Double_t GetCoeff(Int_t i) { return i -! Author(s): Harald Kornmayer 1/2001 -! -! Copyright: MAGIC Software Development, 2000-2001 -! -! -\* ======================================================================== */ - -#include "MMcTimeGenerate.h" - -#include "MLog.h" -#include "MLogManip.h" - -#include "MParList.h" -#include "MTime.h" - -ClassImp(MMcTimeGenerate); - -// -------------------------------------------------------------------------- -// -MMcTimeGenerate::MMcTimeGenerate(const char *name, const char *title) -{ - fName = name ? name : "MMcTimeGenerate"; - fTitle = title ? title : "Task to generate a random event time"; - - const Double_t lambda = 100; // [Hz] - - fFunc = new TF1("Poisson", "[0] * exp(-[0]*x)", 0, 1); - fFunc->SetParameter(0, lambda); - - fDeadTime = 0.1/lambda; -} - -MMcTimeGenerate::~MMcTimeGenerate() -{ - delete fFunc; -} - - -// -------------------------------------------------------------------------- -// -// The PreProcess connects the raw data with this task. It checks if the -// input containers exist, if not a kFalse flag is returned. It also checks -// if the output contaniers exist, if not they are created. -// This task can read either Montecarlo files with multiple trigger -// options, either Montecarlo files with a single trigger option. -// -Int_t MMcTimeGenerate::PreProcess (MParList *pList) -{ - // connect the raw data with this task - - fTime = (MTime*)pList->FindCreateObj("MTime"); - if (!fTime) - return kFALSE; - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// -Int_t MMcTimeGenerate::Process() -{ - Double_t dt; - - do dt = fFunc->GetRandom(); - while (dt < fDeadTime); - -#warning SetTime not valid anymore! -// fTime->SetTime(*fTime+dt/1e4); // [0.1ms] - - return kTRUE; -} Index: trunk/MagicSoft/Mars/mmontecarlo/MMcTimeGenerate.h =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcTimeGenerate.h (revision 7093) +++ (revision ) @@ -1,33 +1,0 @@ -#ifndef MARS_MMcTimeGenerate -#define MARS_MMcTimeGenerate - -#ifndef MARS_MTask -#include "MTask.h" -#endif -#ifndef ROOT_TF1 -#include "TF1.h" -#endif - -class MParList; -class MTime; - -class MMcTimeGenerate : public MTask -{ -private: - MTime *fTime; //! - TF1 *fFunc; //! - - Double_t fDeadTime; - - Int_t PreProcess(MParList *pList); - Int_t Process(); - -public: - MMcTimeGenerate(const char *name=NULL, const char *title=NULL); - - ~MMcTimeGenerate(); - - ClassDef(MMcTimeGenerate, 0) // To generate a random time -}; - -#endif Index: trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.cc (revision 7093) +++ (revision ) @@ -1,358 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * This file is part of MARS, the MAGIC Analysis and Reconstruction -! * Software. It is distributed to you in the hope that it can be a useful -! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. -! * It is distributed WITHOUT ANY WARRANTY. -! * -! * Permission to use, copy, modify and distribute this software and its -! * documentation for any purpose is hereby granted without fee, -! * provided that the above copyright notice appear in all copies and -! * that both that copyright notice and this permission notice appear -! * in supporting documentation. It is provided "as is" without express -! * or implied warranty. -! * -! -! -! Author(s): Marcos Lopez 10/2003 -! -! Copyright: MAGIC Software Development, 2000-2003 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// -// MMcWeightEnergySlopeCalc -// -// Change the spectrum of the MC showers simulated with Corsika (a power law) -// to a new one, which can be either, again a power law but with a different -// spectral index, or a generalizeed spectrum. The new spectrum can be -// pass to this class in different ways: -// 1. Is the new spectrum will be a power law, just introduce the slope -// of this power law. -// 2. Is the new spectrum will have a general shape, different options are -// available: -// a) The new spectrum is pass as a TF1 function -// b) The new spectrum is pass as a char* -// c) The new spectrum is pass as a "interpreted function", i.e., a -// function defined inside a ROOT macro, which will be invoked by the -// ROOT Cint itself. This is the case when we use ROOT macros. -// d) The new spectrum is pass as a "real function", i.e., a -// function defined inside normal c++ file. -// -// Method: -// ------ -// -// -Corsika spectrun: dN/dE = A * E^(a) -// with a = fCorsikaSlope, and A = N/integral{E*de} from ELowLim to EUppLim -// -// -New spectrum: dN/dE = B * g(E) -// where B = N/integral{g*dE} from ELowLim to EUppLim, and N=NumEvents -// -// For converting the spectrum simulated with Corsika to the new one, we apply -// a weight to each event, given by: -// -// W(E) = B/A * g(E)/E^(a) -// -// In the case the new spectrum is simply a power law: dN/dE = B * E^(b), we -// have: -// -// W(E) = B/A * E^(b-a) -// -// (The factor B/A is used in order both the original and new spectrum have -// the same area (i.e. in order they represent the same number of showers)) -// -// Note: -// ------ -// -If the the new spectrum is just a power law (i.e. the user only specify -// the slope), the needed calculations (such as the integral of the -// spectrum) are done analytically. But if the new spectrum is given as a -// TF1 object, the calculations is done numerically. -// -// ToDo: -// ----- -// -Give to the user also the possibility to specify the integral of the -// spectrum as another TF1 object (or C++ function) -// -// -// Input Containers: -// MMcEvt, MMcRunHeader, MMcCorsikaRunHeader -// -// Output Container: -// MWeight [MParameterD] -// -////////////////////////////////////////////////////////////////////////////// - -#include "MMcWeightEnergySpecCalc.h" - -#include "MParList.h" -#include "MLog.h" -#include "MLogManip.h" -#include "MMcEvt.hxx" -#include "MMcRunHeader.hxx" -#include "MMcCorsikaRunHeader.h" -#include "MParameters.h" - -#include "TF1.h" -#include "TGraph.h" - -ClassImp(MMcWeightEnergySpecCalc); - -using namespace std; - - - -void MMcWeightEnergySpecCalc::Init(const char *name, const char *title) -{ - - fName = name ? name : "MMcWeightEnergySpecCalc"; - fTitle = title ? title : "Task to calculate weights to change the energy spectrum"; - - AddToBranchList("MMcEvt.fEnergy"); - - fAllEvtsTriggered = kFALSE; - fTotalNumSimulatedShowers = 0; -} - - -// --------------------------------------------------------------------------- -// -// Constructor. The new spectrum will be just a power law. -// -MMcWeightEnergySpecCalc::MMcWeightEnergySpecCalc(Float_t slope, - const char *name, const char *title) -{ - fNewSpecIsPowLaw = kTRUE; - fNewSlope = slope; - - fNewSpectrum = NULL; - - Init(name,title); -} - -// --------------------------------------------------------------------------- -// -// Constructor. The new spectrum will have a general shape, given by the user -// as a TF1 function. -// -MMcWeightEnergySpecCalc::MMcWeightEnergySpecCalc(const TF1& spectrum, - const char *name, const char *title) -{ - fNewSpecIsPowLaw = kFALSE; - fNewSpectrum = (TF1*)spectrum.Clone(); - - Init(name,title); -} - -// --------------------------------------------------------------------------- -// -// As before, but the function which represent the new spectrum is given as -// a char* . Starting from it, we build a TF1 function -// -MMcWeightEnergySpecCalc::MMcWeightEnergySpecCalc(const char* spectrum, - const char *name, const char *title) -{ - fNewSpecIsPowLaw = kFALSE; - fNewSpectrum = new TF1("NewSpectrum",spectrum); - - Init(name,title); -} - -// --------------------------------------------------------------------------- -// -// As before, but the new spectrum is given as a intrepreted C++ function. -// Starting from it we build a TF1 function. -// This constructor is called for interpreted functions by CINT, i.e., when -// the functions are declared inside a ROOT macro. -// -// NOTE: you muss do a casting to (void*) of the function that you pass to this -// constructor before invoking it in a macro, e.g. -// -// Double_t myfunction(Double_t *x, Double_t *par) -// { -// ... -// } -// -// MMcWeightEnergySpecCalc wcalc((void*)myfunction); -// -// tasklist.AddToList(&wcalc); -// -// otherwise ROOT will invoke the constructor McWeightEnergySpecCalc( -// const char* spectrum, const char *name, const char *title) -// -MMcWeightEnergySpecCalc::MMcWeightEnergySpecCalc(void* function, - const char *name, const char *title) -{ - fNewSpecIsPowLaw = kFALSE; - fNewSpectrum = new TF1("NewSpectrum",function,0,1,1); - - Init(name,title); -} - -// --------------------------------------------------------------------------- -// -// As before, but this is the constructor for real functions, i.e. it is called -// when invoked with the normal C++ compiler, i.e. not inside a ROOT macro. -// -MMcWeightEnergySpecCalc::MMcWeightEnergySpecCalc(Double_t (*function)(Double_t*x, Double_t* par), const Int_t npar, const char *name, const char *title) -{ - fNewSpecIsPowLaw = kFALSE; - fNewSpectrum = new TF1("NewSpectrum",function,0,1,1); - - Init(name,title); -} - - - -// ---------------------------------------------------------------------------- -// -// Destructor. Deletes the cloned fNewSpectrum. -// -MMcWeightEnergySpecCalc::~MMcWeightEnergySpecCalc() -{ - if (fNewSpectrum) - delete fNewSpectrum; -} - - - -// --------------------------------------------------------------------------- -// -// -Int_t MMcWeightEnergySpecCalc::PreProcess (MParList *pList) -{ - - fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); - if (!fMcEvt) - { - *fLog << err << dbginf << "MMcEvt not found... exit." << endl; - return kFALSE; - } - - fWeight = (MParameterD*)pList->FindCreateObj("MParameterD", "MWeight"); - if (!fWeight) - return kFALSE; - - return kTRUE; -} - - - -// ---------------------------------------------------------------------------- -// -// Executed each time a new root file is loaded -// We will need fCorsikaSlope and fE{Upp,Low}Lim to calculate the weights -// -Bool_t MMcWeightEnergySpecCalc::ReInit(MParList *plist) -{ - MMcRunHeader *runheader = (MMcRunHeader*)plist->FindObject("MMcRunHeader"); - if (!runheader) - { - *fLog << err << dbginf << "Error - MMcRunHeader not found... exit." << endl; - return kFALSE; - } - - MMcCorsikaRunHeader *corrunheader = (MMcCorsikaRunHeader*)plist->FindObject("MMcCorsikaRunHeader"); - if (!corrunheader) - { - *fLog << err << dbginf << "Error - MMcCorsikaRunHeader not found... exit." << endl; - return kFALSE; - } - - - fCorsikaSlope = (Double_t)corrunheader->GetSlopeSpec(); - fELowLim = (Double_t)corrunheader->GetELowLim(); - fEUppLim = (Double_t)corrunheader->GetEUppLim(); - - fTotalNumSimulatedShowers += runheader->GetNumSimulatedShowers(); - fAllEvtsTriggered |= runheader->GetAllEvtsTriggered(); - - - - *fLog << inf << "Slope of primaries' energy spectrum of Simulated showers: " << fCorsikaSlope << endl; - *fLog << inf << "Limits of energy range of Simulated showers: " - << fELowLim <<" - " << fEUppLim << endl; - *fLog << inf << "New Slope for Simulated showers: " << fNewSlope << endl; - *fLog << inf << "Total Number of Simulated showers: " - << fTotalNumSimulatedShowers << endl; - *fLog << inf << "Only triggered events avail: " << (fAllEvtsTriggered?"yes":"no") << endl; - - - - // - // Sanity checks to be sure that we won't divide by zero later on - // - if(fCorsikaSlope == -1. || fNewSlope == -1.) - { - *fLog << err << "The Slope of the power law must be different of -1... exit" << endl; - return kFALSE; - } - - - // - // Starting from fCorsikaSlope and fE{Upp,Low}Lim, calculate the integrals - // of both, the original Corsika spectrum and the new one. - // - // - // For the Corsika simulated spectrum (just a power law), we have: - // - fCorSpecInt = ( pow(fEUppLim,1+fCorsikaSlope) - pow(fELowLim,1+fCorsikaSlope) ) / ( 1+fCorsikaSlope ); - - - // - // For the new spectrum: - // - if (fNewSpecIsPowLaw) // just the integral of a power law - { - fNewSpecInt = ( pow(fEUppLim,1+fNewSlope) - pow(fELowLim,1+fNewSlope) )/ ( 1+fNewSlope ); - } - else - { - fNewSpectrum->SetRange(fELowLim, fEUppLim); - - // In this case we have to integrate the new spectrum numerically. We - // could do simply fNewSpectrum->Integral(fELowLim,fEUppLim), but ROOT - // fails integrating up to fEUppLim for a sharp cutoff spectrum - - // - // Trick to calculate the integral numerically (it works better than - // fNewSpectrum->Integral(fELowLim,fEUppLim) (although not perfectlly) - // - fNewSpectrum->SetNpx(1000); - TGraph gr(fNewSpectrum,"i"); - Int_t Npx = gr.GetN(); - Double_t* y = gr.GetY(); - - const Double_t integral = y[Npx-1]; - fNewSpecInt = integral; - } - - - return kTRUE; -} - - - -// ---------------------------------------------------------------------------- -// -// -Int_t MMcWeightEnergySpecCalc::Process() -{ - const Double_t energy = fMcEvt->GetEnergy(); - - const Double_t C = fCorSpecInt / fNewSpecInt; - Double_t weight = C; - - - if (fNewSpecIsPowLaw) - weight *= pow(energy,fNewSlope-fCorsikaSlope); - else - weight *= fNewSpectrum->Eval(energy) / pow(energy,fCorsikaSlope); - - - fWeight->SetVal(weight); - - return kTRUE; -} Index: trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.h =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.h (revision 7093) +++ (revision ) @@ -1,76 +1,0 @@ -#ifndef MARS_MMcWeightEnergySpecCalc -#define MARS_MWeigthEnergySlopeCalc - -#ifndef MARS_MTask -#include "MTask.h" -#endif - - -class MParList; -class MMcEvt; -class MParameterD; -class TF1; - - -class MMcWeightEnergySpecCalc : public MTask -{ - - private: - - const MMcEvt *fMcEvt; - MParameterD *fWeight; - - TF1* fNewSpectrum; // Function with the new spectrum - - Bool_t fNewSpecIsPowLaw; // Tells whether the new spectrum is introduced - //as a TF1 object or not - Double_t fCorsikaSlope; // Slope of energy spectrum generated with Corsika - Double_t fNewSlope; // Slope of the new spectrum (if it is a power law) - Double_t fELowLim; - Double_t fEUppLim; // Limits of energy range for generation - Double_t fCorSpecInt; // Value of the integrals of the Corsika and new - Double_t fNewSpecInt; //spectrum respectively, between fElow and fUppLim - - UInt_t fTotalNumSimulatedShowers; - Bool_t fAllEvtsTriggered; - - void Init(const char *name, const char *title); - Bool_t ReInit(MParList *plist); - - Int_t PreProcess(MParList *pList); - Int_t Process(); - - - public: - - MMcWeightEnergySpecCalc(const Float_t slope, - const char *name=NULL, const char *title=NULL); - - MMcWeightEnergySpecCalc(const TF1& spectrum, - const char *name=NULL, const char *title=NULL); - - MMcWeightEnergySpecCalc(const char* spectrum, - const char *name=NULL, const char *title=NULL); - - MMcWeightEnergySpecCalc(void *function, - const char *name=NULL, const char *title=NULL); - - MMcWeightEnergySpecCalc(Double_t (*function)(Double_t* x, Double_t* par), - const Int_t npar, const char *name=NULL, const char *title=NULL); - - - ~MMcWeightEnergySpecCalc(); - - - ClassDef(MMcWeightEnergySpecCalc, 0) // Task to convert the spectrum of the MC showers simulated with Corsika to a different one, by using weights -}; - -#endif - - - - - - - - Index: trunk/MagicSoft/Mars/sponde.rc =================================================================== --- trunk/MagicSoft/Mars/sponde.rc (revision 7093) +++ trunk/MagicSoft/Mars/sponde.rc (revision 7094) @@ -1,1 +1,3 @@ EstimateEnergy.Rule: (0.380075+(MPointingPos.fZd*MPointingPos.fZd*0.00109028))*pow(MHillas.fSize,0.892462) + +#MMcSpectrumWeight.NewSlope: -2.26