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 Index: trunk/MagicSoft/include-Classes/MMcFormat/MMcEvt.cxx =================================================================== --- trunk/MagicSoft/include-Classes/MMcFormat/MMcEvt.cxx (revision 7093) +++ trunk/MagicSoft/include-Classes/MMcFormat/MMcEvt.cxx (revision 7094) @@ -1,10 +1,30 @@ -#include "MMcEvt.hxx" - -#include "MLog.h" -#include "MLogManip.h" - -//========== -// MMcEvt -// +/* ======================================================================== *\ +! +! * +! * 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): +! +! Copyright: MAGIC Software Development, 2000-2005 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MMcEvt +// // This class handles and contains the MonteCarlo information // with which the events have been generated @@ -19,6 +39,20 @@ // So, fTelescopeTheta=90, fTelescopePhi = 0 means the telescope is // pointing horizontally towards South. For an explanation, see also -// TDAS 02-11. -// +// TDAS 02-11. +// +// Version 4: +// - Added member fFadcTimeJitter +// +// Version 5: +// - removed fPartId, fEnergy, fImpact, fTelescopeTheta, fTelescopePhi +// - derives now from MMcEvtBasic which contains all these values +// - moved ParticleId_t to base class MMcEvtBasic +// +///////////////////////////////////////////////////////////////////////////// +#include "MMcEvt.hxx" + +#include "MLog.h" +#include "MLogManip.h" + ClassImp(MMcEvt); @@ -26,9 +60,10 @@ +// -------------------------------------------------------------------------- +// +// Default constructor. Calls Clear() +// MMcEvt::MMcEvt() { - // - // default constructor - // set all values to zero fName = "MMcEvt"; fTitle = "Event info from Monte Carlo"; @@ -37,107 +72,40 @@ } -MMcEvt::MMcEvt( UInt_t fEvtNum, - UShort_t usPId, - Float_t fEner, - Float_t fThi0, - Float_t fFirTar, - Float_t fzFirInt, - Float_t fThet, - Float_t fPhii, - Float_t fCorD, - Float_t fCorX, - Float_t fCorY, - Float_t fImpa, - Float_t fTPhii, - Float_t fTThet, - Float_t fTFirst, - Float_t fTLast, - Float_t fL_Nmax, - Float_t fL_t0, - Float_t fL_tmax, - Float_t fL_a, - Float_t fL_b, - Float_t fL_c, - Float_t fL_chi2, - UInt_t uiPin, - UInt_t uiPat, - UInt_t uiPre, - UInt_t uiPco, - UInt_t uiPelS, - UInt_t uiPelC, - Float_t elec, - Float_t muon, - Float_t other, - Float_t fadc_jitter) { - +// -------------------------------------------------------------------------- +// +// Constructor. Use this to set all data members +// +// THIS FUNCTION IS FOR THE SIMULATION OLNY. +// DON'T USE THIS MEMBERFUNCTION IN THE ANALYSIS. +// +MMcEvt::MMcEvt(UInt_t fEvtNum, ParticleId_t usPId, Float_t fEner, + Float_t fThi0, Float_t fFirTar, Float_t fzFirInt, + Float_t fThet, Float_t fPhii, Float_t fCorD, + Float_t fCorX, Float_t fCorY, Float_t fImpa, + Float_t fTPhii, Float_t fTThet, Float_t fTFirst, + Float_t fTLast, Float_t fL_Nmax, Float_t fL_t0, + Float_t fL_tmax, Float_t fL_a, Float_t fL_b, + Float_t fL_c, Float_t fL_chi2, UInt_t uiPin, + UInt_t uiPat, UInt_t uiPre, UInt_t uiPco, + UInt_t uiPelS, UInt_t uiPelC, Float_t elec, + Float_t muon, Float_t other, Float_t fadc_jitter) +{ fName = "MMcEvt"; fTitle = "Event info from Monte Carlo"; - // - // constuctor II - // - // All datamembers are parameters. - // - // Don't use this memberfunction in analysis - // - - fEvtNumber = fEvtNum; - fPartId = usPId ; - fEnergy = fEner ; - fThick0 = fThi0; - fFirstTarget = fFirTar; - fZFirstInteraction = fzFirInt; - - fTheta = fThet ; - fPhi = fPhii ; - - fCoreD = fCorD ; - fCoreX = fCorX ; - fCoreY = fCorY ; - fImpact = fImpa ; - - fTelescopePhi = fTPhii; - fTelescopeTheta = fTThet; - fTimeFirst = fTFirst; - fTimeLast = fTLast; - fLongiNmax = fL_Nmax; - fLongit0 = fL_t0; - fLongitmax = fL_tmax; - fLongia = fL_a; - fLongib = fL_b; - fLongic = fL_c; - fLongichi2 = fL_chi2; - - - fPhotIni = uiPin ; - fPassPhotAtm = uiPat ; - fPassPhotRef = uiPre ; - fPassPhotCone = uiPco ; - fPhotElfromShower = uiPelS ; - fPhotElinCamera = uiPelC ; - - fElecCphFraction=elec; - fMuonCphFraction=muon; - fOtherCphFraction=other; - - fFadcTimeJitter = fadc_jitter; -} - - - -MMcEvt::~MMcEvt() { - // - // default destructor - // -} - - - - + + Fill(fEvtNum, usPId, fEner, fThi0, fFirTar, fzFirInt, fThet, + fPhii, fCorD, fCorX, fCorY, fImpa, fTPhii, fTThet, fTFirst, + fTLast, fL_Nmax, fL_t0, fL_tmax, fL_a, fL_b, fL_c, fL_chi2, + uiPin, uiPat, uiPre, uiPco, uiPelS, uiPelC, elec, muon, other, + fadc_jitter); +} + +// -------------------------------------------------------------------------- +// +// reset all values to values as nonsense as possible +// void MMcEvt::Clear(Option_t *opt) { - // - // reset all values to values as nonsense as possible - // - fPartId = 0; + fPartId = kUNDEFINED; fEnergy = -1; @@ -162,91 +130,68 @@ } -void MMcEvt::Fill( UInt_t fEvtNum, - UShort_t usPId, - Float_t fEner, - Float_t fThi0, - Float_t fFirTar, - Float_t fzFirInt, - Float_t fThet, - Float_t fPhii, - Float_t fCorD, - Float_t fCorX, - Float_t fCorY, - Float_t fImpa, - Float_t fTPhii, - Float_t fTThet, - Float_t fTFirst, - Float_t fTLast, - Float_t fL_Nmax, - Float_t fL_t0, - Float_t fL_tmax, - Float_t fL_a, - Float_t fL_b, - Float_t fL_c, - Float_t fL_chi2, - UInt_t uiPin, - UInt_t uiPat, - UInt_t uiPre, - UInt_t uiPco, - UInt_t uiPelS, - UInt_t uiPelC, - Float_t elec, - Float_t muon, - Float_t other, - Float_t fadc_jitter) { - // - // All datamembers are filled with the correspondin parameters. - // - // Don't use this memberfunction in analysis - // - - fEvtNumber = fEvtNum; - fPartId = usPId ; - fEnergy = fEner ; - fThick0 = fThi0; - fFirstTarget = fFirTar; - fZFirstInteraction = fzFirInt; - - fTheta = fThet ; - fPhi = fPhii ; - - fCoreD = fCorD ; - fCoreX = fCorX ; - fCoreY = fCorY ; - fImpact = fImpa ; - - fTelescopePhi = fTPhii; - fTelescopeTheta = fTThet; - fTimeFirst = fTFirst; - fTimeLast = fTLast; - fLongiNmax = fL_Nmax; - fLongit0 = fL_t0; - fLongitmax = fL_tmax; - fLongia = fL_a; - fLongib = fL_b; - fLongic = fL_c; - fLongichi2 = fL_chi2; - - fPhotIni = uiPin ; - fPassPhotAtm = fPhotIni-uiPat ; - fPassPhotRef = fPassPhotAtm-uiPre ; - fPassPhotCone = uiPco ; - fPhotElfromShower = uiPelS ; - fPhotElinCamera = uiPelC ; - - fElecCphFraction=elec; - fMuonCphFraction=muon; - fOtherCphFraction=other; - - fFadcTimeJitter = fadc_jitter; -} - -/* -void MMcEvt::AsciiWrite(ofstream &fout) const -{ - fout << fEnergy << " "; - fout << fTheta ; -} -*/ +// -------------------------------------------------------------------------- +// +// Use this to set all data members +// +// THIS FUNCTION IS FOR THE SIMULATION OLNY. +// DON'T USE THIS MEMBERFUNCTION IN THE ANALYSIS. +// +void MMcEvt::Fill( UInt_t fEvtNum, ParticleId_t usPId, Float_t fEner, + Float_t fThi0, Float_t fFirTar, Float_t fzFirInt, + Float_t fThet, Float_t fPhii, Float_t fCorD, + Float_t fCorX, Float_t fCorY, Float_t fImpa, + Float_t fTPhii, Float_t fTThet, Float_t fTFirst, + Float_t fTLast, Float_t fL_Nmax, Float_t fL_t0, + Float_t fL_tmax, Float_t fL_a, Float_t fL_b, + Float_t fL_c, Float_t fL_chi2, UInt_t uiPin, + UInt_t uiPat, UInt_t uiPre, UInt_t uiPco, + UInt_t uiPelS, UInt_t uiPelC, Float_t elec, + Float_t muon, Float_t other, Float_t fadc_jitter) +{ + // + // All datamembers are filled with the correspondin parameters. + // + // Don't use this memberfunction in analysis + // + fEvtNumber = fEvtNum; + fPartId = usPId ; + fEnergy = fEner ; + fThick0 = fThi0; + fFirstTarget = fFirTar; + fZFirstInteraction = fzFirInt; + + fTheta = fThet ; + fPhi = fPhii ; + + fCoreD = fCorD ; + fCoreX = fCorX ; + fCoreY = fCorY ; + fImpact = fImpa ; + + fTelescopePhi = fTPhii; + fTelescopeTheta = fTThet; + fTimeFirst = fTFirst; + fTimeLast = fTLast; + fLongiNmax = fL_Nmax; + fLongit0 = fL_t0; + fLongitmax = fL_tmax; + fLongia = fL_a; + fLongib = fL_b; + fLongic = fL_c; + fLongichi2 = fL_chi2; + + fPhotIni = uiPin ; + fPassPhotAtm = fPhotIni-uiPat ; + fPassPhotRef = fPassPhotAtm-uiPre ; + fPassPhotCone = uiPco ; + fPhotElfromShower = uiPelS ; + fPhotElinCamera = uiPelC ; + + fElecCphFraction=elec; + fMuonCphFraction=muon; + fOtherCphFraction=other; + + fFadcTimeJitter = fadc_jitter; +} // -------------------------------------------------------------------------- @@ -260,47 +205,11 @@ void MMcEvt::Print(Option_t *opt) const { - // - // print out the data member on screen - // + MMcEvtBasic::Print(opt); + TString str(opt); if (str.IsNull()) { - *fLog << all << endl; - *fLog << "Monte Carlo output:" << endl; - *fLog << " Particle Id: "; - switch(fPartId) - { - case kGAMMA: - *fLog << "Gamma" << endl; - break; - case kPROTON: - *fLog << "Proton" << endl; - break; - case kHELIUM: - *fLog << "Helium" << endl; - break; - } - *fLog << " Energy: " << fEnergy << "GeV" << endl; - *fLog << " Impactpar.: " << fImpact/100 << "m" << endl; *fLog << " Photoelectrons: " << fPhotElfromShower << endl; - *fLog << endl; return; } - if (str.Contains("id", TString::kIgnoreCase)) - switch(fPartId) - { - case kGAMMA: - *fLog << "Particle: Gamma" << endl; - break; - case kPROTON: - *fLog << "Particle: Proton" << endl; - break; - case kHELIUM: - *fLog << "Particle: Helium" << endl; - break; - } - if (str.Contains("energy", TString::kIgnoreCase)) - *fLog << "Energy: " << fEnergy << "GeV" << endl; - if (str.Contains("impact", TString::kIgnoreCase)) - *fLog << "Impact: " << fImpact << "cm" << endl; -} +} Index: trunk/MagicSoft/include-Classes/MMcFormat/MMcEvt.hxx =================================================================== --- trunk/MagicSoft/include-Classes/MMcFormat/MMcEvt.hxx (revision 7093) +++ trunk/MagicSoft/include-Classes/MMcFormat/MMcEvt.hxx (revision 7094) @@ -1,244 +1,117 @@ -#ifndef __MMcEvt__ -#define __MMcEvt__ +#ifndef MARS_MMcEvt +#define MARS_MMcEvt -#ifndef MARS_MParContainer -#include "MParContainer.h" +#ifndef MARS_MMcEvtBasic +#include "MMcEvtBasic.h" #endif -// -// Version 4: -// Added member fFadcTimeJitter -// -class MMcEvt : public MParContainer +class MMcEvt : public MMcEvtBasic { +private: + UInt_t fEvtNumber; + Float_t fThick0; // [g/cm2] + Float_t fFirstTarget; // [] + Float_t fZFirstInteraction; // [cm] + + Float_t fTheta; // [rad] Theta angle of event + Float_t fPhi; // [rad] Phi angle of event (see class description) + + Float_t fCoreD; // [cm] Core d pos + Float_t fCoreX; // [cm] Core x pos + Float_t fCoreY; // [cm] Core y pos + + // Up to here, the info from the CORSIKA event header. + + // Time of first and last photon: + Float_t fTimeFirst; // [ns] + Float_t fTimeLast; // [ns] + + // 6 parameters and chi2 of the NKG fit to the longitudinal + // particle distribution. See CORSIKA manual for explanation, + // section 4.42 "Longitudinal shower development": + // + Float_t fLongiNmax; // [particles] + Float_t fLongit0; // [g/cm2] + Float_t fLongitmax; // [g/cm2] + Float_t fLongia; // [g/cm2] + Float_t fLongib; // [] + Float_t fLongic; // [cm2/g] + Float_t fLongichi2; + + UInt_t fPhotIni; // [ph] Initial number of photons + UInt_t fPassPhotAtm; // [ph] Passed atmosphere + UInt_t fPassPhotRef; // [ph] Passed reflector(reflectivity + effective area) + UInt_t fPassPhotCone; // [ph] Within any valid pixel, before plexiglas + UInt_t fPhotElfromShower; // [phe] Passed qe, coming from the shower + UInt_t fPhotElinCamera; // [phe] usPhotElfromShower + mean of phe from NSB + + // Now follow the fraction of photons reaching the camera produced by + // electrons, muons and other particles respectively: + + Float_t fElecCphFraction; + Float_t fMuonCphFraction; + Float_t fOtherCphFraction; + + Float_t fFadcTimeJitter; + public: - // - // ParticleId for Monte Carlo simulation - // - enum ParticleId_t - { - kGAMMA = 1, - kPOSITRON = 2, - kELECTRON = 3, - kANTIMUON = 5, - kMUON = 6, - kPI0 = 7, - kNEUTRON = 13, - kPROTON = 14, - kHELIUM = 402, - kOXYGEN = 1608, - kIRON = 5626 - }; + MMcEvt(); + MMcEvt(UInt_t, ParticleId_t, Float_t, Float_t, Float_t, + Float_t, Float_t, Float_t, Float_t, Float_t, Float_t, + Float_t, Float_t, Float_t, Float_t, Float_t, Float_t, + Float_t, Float_t, Float_t, Float_t, Float_t, Float_t, + UInt_t, UInt_t, UInt_t, UInt_t, UInt_t, UInt_t, + Float_t, Float_t, Float_t, Float_t) ; -private: - UInt_t fEvtNumber; - UShort_t fPartId; // Type of particle - Float_t fEnergy; // [GeV] Energy - Float_t fThick0; // [g/cm2] - Float_t fFirstTarget; // [] - Float_t fZFirstInteraction; // [cm] + // Getter + UInt_t GetEvtNumber() const { return fEvtNumber; } //Get Event Number + Float_t GetTheta() const { return fTheta; } //Get Theta angle + Float_t GetPhi() const { return fPhi ; } //Get Phi angle - Float_t fTheta; // [rad] Theta angle of event - Float_t fPhi; // [rad] Phi angle of event (see class description) + Float_t GetCoreX() const { return fCoreX; } //Get Core x pos + Float_t GetCoreY() const { return fCoreY; } //Get Core y pos - Float_t fCoreD; // [cm] Core d pos - Float_t fCoreX; // [cm] Core x pos - Float_t fCoreY; // [cm] Core y pos - Float_t fImpact; // [cm] impact parameter + UInt_t GetPhotIni() const { return fPhotIni; } //Get Initial photons + UInt_t GetPassPhotAtm() const { return fPassPhotAtm;} //Get Passed atmosphere + UInt_t GetPassPhotRef() const { return fPassPhotRef; } //Get Passed reflector + UInt_t GetPassPhotCone() const { return fPassPhotCone; } //Get Passed glas + UInt_t GetPhotElfromShower() const { return fPhotElfromShower; } //Get Passed qe from shower + UInt_t GetPhotElinCamera() const { return fPhotElinCamera; } //Get Passed qe total + Float_t GetZFirstInteraction() const { return fZFirstInteraction; } - // Up to here, the info from the CORSIKA event header. + Float_t GetOtherCphFraction() const { return fOtherCphFraction; } - // Telescope orientation: - Float_t fTelescopePhi; // [rad] (see class description) - Float_t fTelescopeTheta; // [rad] + Float_t GetLongiNmax() const { return fLongiNmax; } + Float_t GetLongia() const { return fLongia; } + Float_t GetLongib() const { return fLongib; } + Float_t GetLongic() const { return fLongic; } + Float_t GetLongichi2() const { return fLongichi2; } + Float_t GetLongit0() const { return fLongit0; } + Float_t GetLongitmax() const { return fLongitmax; } - // Time of first and last photon: - Float_t fTimeFirst; // [ns] - Float_t fTimeLast; // [ns] + Float_t GetFadcTimeJitter() const { return fFadcTimeJitter; } - // 6 parameters and chi2 of the NKG fit to the longitudinal - // particle distribution. See CORSIKA manual for explanation, - // section 4.42 "Longitudinal shower development": - // - Float_t fLongiNmax; // [particles] - Float_t fLongit0; // [g/cm2] - Float_t fLongitmax; // [g/cm2] - Float_t fLongia; // [g/cm2] - Float_t fLongib; // [] - Float_t fLongic; // [cm2/g] - Float_t fLongichi2; + Float_t GetMuonCphFraction() const { return fMuonCphFraction; } - UInt_t fPhotIni; // [ph] Initial number of photons - UInt_t fPassPhotAtm; // [ph] Passed atmosphere - UInt_t fPassPhotRef; // [ph] Passed reflector(reflectivity + effective area) - UInt_t fPassPhotCone; // [ph] Within any valid pixel, before plexiglas - UInt_t fPhotElfromShower; // [phe] Passed qe, coming from the shower - UInt_t fPhotElinCamera; // [phe] usPhotElfromShower + mean of phe - // from NSB + // Setter + void SetTheta(Float_t Theta) { fTheta=Theta; } //Set Theta angle + void SetPhi(Float_t Phi) { fPhi=Phi; } //Set Phi angle + void SetCoreD(Float_t CoreD) { fCoreD=CoreD; } //Set Core d pos + void SetCoreX(Float_t CoreX) { fCoreX=CoreX; } //Set Core x pos + void SetCoreY(Float_t CoreY) { fCoreY=CoreY; } //Set Core y pos - // Now follow the fraction of photons reaching the camera produced by - // electrons, muons and other particles respectively: + void Fill( UInt_t, ParticleId_t, Float_t, Float_t, Float_t, + Float_t, Float_t, Float_t, Float_t, Float_t, Float_t, + Float_t, Float_t, Float_t, Float_t, Float_t, Float_t, + Float_t, Float_t, Float_t, Float_t, Float_t, Float_t, + UInt_t, UInt_t, UInt_t, UInt_t, UInt_t, UInt_t, + Float_t, Float_t, Float_t, Float_t) ; - Float_t fElecCphFraction; - Float_t fMuonCphFraction; - Float_t fOtherCphFraction; - - Float_t fFadcTimeJitter; + // TObject + void Print(Option_t *opt=NULL) const; + void Clear(Option_t *opt=NULL); - public: - MMcEvt() ; - - MMcEvt( UInt_t, UShort_t, Float_t, Float_t, Float_t, - Float_t, Float_t, Float_t, Float_t, Float_t, Float_t, - Float_t, Float_t, Float_t, Float_t, Float_t, Float_t, - Float_t, Float_t, Float_t, Float_t, Float_t, Float_t, - UInt_t, UInt_t, UInt_t, UInt_t, UInt_t, UInt_t, - Float_t, Float_t, Float_t, Float_t) ; - - ~MMcEvt(); - - void Clear(Option_t *opt=NULL); - - void Fill( UInt_t, UShort_t, Float_t, Float_t, Float_t, - Float_t, Float_t, Float_t, Float_t, Float_t, Float_t, - Float_t, Float_t, Float_t, Float_t, Float_t, Float_t, - Float_t, Float_t, Float_t, Float_t, Float_t, Float_t, - UInt_t, UInt_t, UInt_t, UInt_t, UInt_t, UInt_t, - Float_t, Float_t, Float_t, Float_t) ; - - //virtual void AsciiWrite(ofstream &fout) const; - - void Print(Option_t *opt=NULL) const; - - UInt_t GetEvtNumber() const { return fEvtNumber; } //Get Event Number - Short_t GetPartId() const { return fPartId; } //Get Type of particle - Float_t GetEnergy() const { return fEnergy; } //Get Energy - - Float_t GetTheta() const { return fTheta; } //Get Theta angle - Float_t GetPhi() const { return fPhi ; } //Get Phi angle - -/* Float_t GetCoreD() { return fCoreD; } //Get Core d pos */ - Float_t GetCoreX() { return fCoreX; } //Get Core x pos - Float_t GetCoreY() { return fCoreY; } //Get Core y pos - Float_t GetImpact() const { return fImpact;} //Get impact parameter - - UInt_t GetPhotIni() { return fPhotIni; } //Get Initial photons - UInt_t GetPassPhotAtm() { return fPassPhotAtm;} //Get Passed atmosphere - UInt_t GetPassPhotRef() { return fPassPhotRef; } //Get Passed reflector - UInt_t GetPassPhotCone() { return fPassPhotCone; } //Get Passed glas - UInt_t GetPhotElfromShower() { return fPhotElfromShower; } //Get Passed qe from shower - UInt_t GetPhotElinCamera() { return fPhotElinCamera; } //Get Passed qe total - Float_t GetZFirstInteraction() const { return fZFirstInteraction; } - - Float_t GetTelescopePhi() const { return fTelescopePhi; } - Float_t GetTelescopeTheta() const { return fTelescopeTheta; } - Float_t GetOtherCphFraction() const { return fOtherCphFraction; } - - Float_t GetLongiNmax() const { return fLongiNmax; } - Float_t GetLongia() const { return fLongia; } - Float_t GetLongib() const { return fLongib; } - Float_t GetLongic() const { return fLongic; } - Float_t GetLongichi2() const { return fLongichi2; } - Float_t GetLongit0() const { return fLongit0; } - Float_t GetLongitmax() const { return fLongitmax; } - - Float_t GetFadcTimeJitter() const { return fFadcTimeJitter; } - - Float_t GetMuonCphFraction() const { return fMuonCphFraction; } - - void SetPartId(Short_t PartId) - {fPartId=PartId;} //Set Type of particle - - TString GetParticleName() const - { - switch (fPartId) - { - case kGAMMA: return "Gamma"; - case kPOSITRON: return "Positron"; - case kELECTRON: return "Electron"; - case kANTIMUON: return "Anti-Muon"; - case kMUON: return "Muon"; - case kPI0: return "Pi-0"; - case kNEUTRON: return "Neutron"; - case kPROTON: return "Proton"; - case kHELIUM: return "Helium"; - case kOXYGEN: return "Oxygen"; - case kIRON: return "Iron"; - } - - return Form("Id:%d", fPartId); - } - TString GetParticleSymbol() const - { - switch (fPartId) - { - case kGAMMA: return "\\gamma"; - case kPOSITRON: return "e^{+}"; - case kELECTRON: return "e^{-}"; - case kANTIMUON: return "\\mu^{+}"; - case kMUON: return "\\mu^{-}"; - case kPI0: return "\\pi^{0}"; - case kNEUTRON: return "n"; - case kPROTON: return "p"; - case kHELIUM: return "He"; - case kOXYGEN: return "O"; - case kIRON: return "Fe"; - } - - return Form("Id:%d", fPartId); - } - TString GetEnergyStr() const - { - if (fEnergy>1000) - return Form("%.1fTeV", fEnergy/1000); - - if (fEnergy>10) - return Form("%dGeV", (Int_t)(fEnergy+.5)); - - if (fEnergy>1) - return Form("%.1fGeV", fEnergy); - - return Form("%dMeV", (Int_t)(fEnergy*1000+.5)); - } - - void SetEnergy(Float_t Energy) - { fEnergy=Energy; } //Set Energy - - void SetTheta(Float_t Theta) - { fTheta=Theta; } //Set Theta angle - - void SetPhi(Float_t Phi) - { fPhi=Phi; } //Set Phi angle - - void SetCoreD(Float_t CoreD) - { fCoreD=CoreD; } //Set Core d pos - - void SetCoreX(Float_t CoreX) - { fCoreX=CoreX; } //Set Core x pos - - void SetCoreY(Float_t CoreY ) - { fCoreY=CoreY; } //Set Core y pos - - void SetImpact(Float_t Impact) - { fImpact=Impact;} //Set impact parameter - - // DO NOT USE THIS IS A WORKAROUND! - void SetTelescopeTheta(Float_t Theta) { fTelescopeTheta=Theta; } - void SetTelescopePhi(Float_t Phi) { fTelescopePhi=Phi; } - - -/* void SetPhotIni(Short_t PhotIni) */ -/* { fPhotIni=PhotIni; } //Set Initial photons */ -/* void SetPassPhotAtm(Short_t PassPhotAtm) */ -/* { fPassPhotAtm=PassPhotAtm;} //Set Passed atmosphere */ -/* void SetPassPhotRef(Short_t PassPhotRef) */ -/* { fPassPhotRef=PassPhotRef ; } //Set Passed reflector */ -/* void SetPassPhotCone(Short_t PhotCon) */ -/* { fPassPhotCone=PhotCon; } //Set Passed glas */ - - - ClassDef(MMcEvt, 4) //Stores Montecarlo Information of one event (eg. the energy) - + ClassDef(MMcEvt, 5) //Stores Montecarlo Information of one event (eg. the energy) }; Index: trunk/MagicSoft/include-Classes/MMcFormat/MMcEvtBasic.cc =================================================================== --- trunk/MagicSoft/include-Classes/MMcFormat/MMcEvtBasic.cc (revision 7093) +++ trunk/MagicSoft/include-Classes/MMcFormat/MMcEvtBasic.cc (revision 7094) @@ -40,6 +40,14 @@ // South. // +// +// Version 1: +// New container to keep the very basic informations on the +// original MC events produced by Corsika +// +// Version 2: +// - added typedef for ParticleId_t from MMcEvt +// - replaced MMcEvt::ParticleId_t by ParticleId_t +// ///////////////////////////////////////////////////////////////////////////// - #include "MMcEvtBasic.h" @@ -52,55 +60,38 @@ // -------------------------------------------------------------------------- -// Default constructor set all values to zero +// +// Default constructor. Calls Clear() // MMcEvtBasic::MMcEvtBasic() -{ - fName = "MMcEvtBasic"; - fTitle = "Basic event info from Monte Carlo"; +{ + fName = "MMcEvtBasic"; + fTitle = "Basic event info from Monte Carlo"; - Clear(); + Clear(); } // -------------------------------------------------------------------------- -// constuctor II // -// All datamembers are parameters. +// Constructor. Use this to set all data members // -MMcEvtBasic::MMcEvtBasic(MMcEvt::ParticleId_t usPId, - Float_t fEner, - Float_t fImpa, - Float_t fTPhii, - Float_t fTThet) +// THIS FUNCTION IS FOR THE SIMULATION OLNY. +// DON'T USE THIS MEMBERFUNCTION IN THE ANALYSIS. +// +MMcEvtBasic::MMcEvtBasic(ParticleId_t usPId, Float_t fEner, + Float_t fImpa, Float_t fTPhii, Float_t fTThet) { - fName = "MMcEvtBasic"; - fTitle = "Basic event info from Monte Carlo"; + fName = "MMcEvtBasic"; + fTitle = "Basic event info from Monte Carlo"; - fPartId = usPId; - fEnergy = fEner; - fImpact = fImpa; - fTelescopePhi = fTPhii; - fTelescopeTheta = fTThet; + Fill(usPId, fEner, fImpa, fTPhii, fTThet); } - - -// -------------------------------------------------------------------------- -// default destructor -// -MMcEvtBasic::~MMcEvtBasic() { -} - // -------------------------------------------------------------------------- // -// Reset all values. We have no "error" value for fPartId, -// we just set kGAMMA +// Reset all values: Fill(kUNDEFINED, -1, -1, 0, 0) // void MMcEvtBasic::Clear(Option_t *opt) { - fPartId = MMcEvt::kGAMMA; - fEnergy = -1.; - fImpact = -1.; - fTelescopeTheta = -1.; - fTelescopePhi = -1.; + Fill(kUNDEFINED, -1, -1, 0, 0); } @@ -109,15 +100,14 @@ // Fill all data members // -void MMcEvtBasic::Fill(MMcEvt::ParticleId_t usPId, - Float_t fEner, - Float_t fImpa, - Float_t fTPhii, - Float_t fTThet) +void MMcEvtBasic::Fill(ParticleId_t usPId, Float_t fEner, + Float_t fImpa, Float_t fTPhii, Float_t fTThet) { - fPartId = usPId; - fEnergy = fEner; - fImpact = fImpa; - fTelescopePhi = fTPhii; - fTelescopeTheta = fTThet; + fPartId = usPId; + + fEnergy = fEner; + fImpact = fImpa; + + fTelescopePhi = fTPhii; + fTelescopeTheta = fTThet; } @@ -131,50 +121,23 @@ void MMcEvtBasic::Print(Option_t *opt) const { - // - // print out the data member on screen - // - TString str(opt); - if (str.IsNull()) + // + // print out the data member on screen + // + TString str(opt); + if (str.IsNull()) { - *fLog << all << endl; - *fLog << "Monte Carlo output:" << endl; - *fLog << " Particle Id: "; - switch(fPartId) - { - case MMcEvt::kGAMMA: - *fLog << "Gamma" << endl; - break; - case MMcEvt::kPROTON: - *fLog << "Proton" << endl; - break; - case MMcEvt::kHELIUM: - *fLog << "Helium" << endl; - break; - default: - break; - } - *fLog << " Energy: " << fEnergy << "GeV" << endl; - *fLog << " Impactpar.: " << fImpact/100 << "m" << endl; - *fLog << endl; - return; + *fLog << all << endl; + *fLog << "Monte Carlo output:" << endl; + *fLog << " Particle Id: " << GetParticleName() << endl; + *fLog << " Energy: " << fEnergy << "GeV" << endl; + *fLog << " Impactparam.: " << fImpact/100 << "m" << endl; + *fLog << endl; + return; } - if (str.Contains("id", TString::kIgnoreCase)) - switch(fPartId) - { - case MMcEvt::kGAMMA: - *fLog << "Particle: Gamma" << endl; - break; - case MMcEvt::kPROTON: - *fLog << "Particle: Proton" << endl; - break; - case MMcEvt::kHELIUM: - *fLog << "Particle: Helium" << endl; - break; - default: - break; - } - if (str.Contains("energy", TString::kIgnoreCase)) - *fLog << "Energy: " << fEnergy << "GeV" << endl; - if (str.Contains("impact", TString::kIgnoreCase)) - *fLog << "Impact: " << fImpact << "cm" << endl; + if (str.Contains("id", TString::kIgnoreCase)) + *fLog << "Particle: " << GetParticleName() << endl; + if (str.Contains("energy", TString::kIgnoreCase)) + *fLog << "Energy: " << fEnergy << "GeV" << endl; + if (str.Contains("impact", TString::kIgnoreCase)) + *fLog << "Impact: " << fImpact << "cm" << endl; } Index: trunk/MagicSoft/include-Classes/MMcFormat/MMcEvtBasic.h =================================================================== --- trunk/MagicSoft/include-Classes/MMcFormat/MMcEvtBasic.h (revision 7093) +++ trunk/MagicSoft/include-Classes/MMcFormat/MMcEvtBasic.h (revision 7094) @@ -1,4 +1,4 @@ -#ifndef __MMcEvtBasic__ -#define __MMcEvtBasic__ +#ifndef MARS_MMcEvtBasic +#define MARS_MMcEvtBasic #ifndef MARS_MParContainer @@ -6,41 +6,44 @@ #endif -#include "MMcEvt.hxx" - -// -// Version 1: -// New container to keep the very basic informations on the -// original MC events produced by Corsika -// class MMcEvtBasic : public MParContainer { +public: + enum ParticleId_t + { + kUNDEFINED = -1, + kGAMMA = 1, + kPOSITRON = 2, + kELECTRON = 3, + kANTIMUON = 5, + kMUON = 6, + kPI0 = 7, + kNEUTRON = 13, + kPROTON = 14, + kHELIUM = 402, + kOXYGEN = 1608, + kIRON = 5626 + }; private: - MMcEvt::ParticleId_t fPartId; // Type of particle - Float_t fEnergy; // [GeV] Energy - Float_t fImpact; // [cm] impact parameter + ParticleId_t fPartId; // Type of particle + Float_t fEnergy; // [GeV] Energy + Float_t fImpact; // [cm] impact parameter // Telescope orientation (see TDAS 02-11 regarding the // precise meaning of these angles): - Float_t fTelescopePhi; // [rad] - Float_t fTelescopeTheta; // [rad] - + Float_t fTelescopePhi; // [rad] + Float_t fTelescopeTheta; // [rad] public: MMcEvtBasic(); - - MMcEvtBasic(MMcEvt::ParticleId_t, Float_t, Float_t, Float_t, Float_t); - ~MMcEvtBasic(); + MMcEvtBasic(ParticleId_t, Float_t, Float_t, Float_t, Float_t); - void Clear(Option_t *opt=NULL); + // Getter + ParticleId_t GetPartId() const { return fPartId; } - void Fill(MMcEvt::ParticleId_t, Float_t, Float_t, Float_t, Float_t); - - void Print(Option_t *opt=NULL) const; - - MMcEvt::ParticleId_t GetPartId() const { return fPartId; } Float_t GetEnergy() const { return fEnergy; } Float_t GetImpact() const { return fImpact; } + Float_t GetTelescopePhi() const { return fTelescopePhi; } Float_t GetTelescopeTheta() const { return fTelescopeTheta; } @@ -50,15 +53,16 @@ switch (fPartId) { - case MMcEvt::kGAMMA: return "Gamma"; - case MMcEvt::kPOSITRON: return "Positron"; - case MMcEvt::kELECTRON: return "Electron"; - case MMcEvt::kANTIMUON: return "Anti-Muon"; - case MMcEvt::kMUON: return "Muon"; - case MMcEvt::kPI0: return "Pi-0"; - case MMcEvt::kNEUTRON: return "Neutron"; - case MMcEvt::kPROTON: return "Proton"; - case MMcEvt::kHELIUM: return "Helium"; - case MMcEvt::kOXYGEN: return "Oxygen"; - case MMcEvt::kIRON: return "Iron"; + case kUNDEFINED:return "Undefined"; + case kGAMMA: return "Gamma"; + case kPOSITRON: return "Positron"; + case kELECTRON: return "Electron"; + case kANTIMUON: return "Anti-Muon"; + case kMUON: return "Muon"; + case kPI0: return "Pi-0"; + case kNEUTRON: return "Neutron"; + case kPROTON: return "Proton"; + case kHELIUM: return "Helium"; + case kOXYGEN: return "Oxygen"; + case kIRON: return "Iron"; } @@ -70,15 +74,16 @@ switch (fPartId) { - case MMcEvt::kGAMMA: return "\\gamma"; - case MMcEvt::kPOSITRON: return "e^{+}"; - case MMcEvt::kELECTRON: return "e^{-}"; - case MMcEvt::kANTIMUON: return "\\mu^{+}"; - case MMcEvt::kMUON: return "\\mu^{-}"; - case MMcEvt::kPI0: return "\\pi^{0}"; - case MMcEvt::kNEUTRON: return "n"; - case MMcEvt::kPROTON: return "p"; - case MMcEvt::kHELIUM: return "He"; - case MMcEvt::kOXYGEN: return "O"; - case MMcEvt::kIRON: return "Fe"; + case kUNDEFINED:return "N/A"; + case kGAMMA: return "\\gamma"; + case kPOSITRON: return "e^{+}"; + case kELECTRON: return "e^{-}"; + case kANTIMUON: return "\\mu^{+}"; + case kMUON: return "\\mu^{-}"; + case kPI0: return "\\pi^{0}"; + case kNEUTRON: return "n"; + case kPROTON: return "p"; + case kHELIUM: return "He"; + case kOXYGEN: return "O"; + case kIRON: return "Fe"; } @@ -101,18 +106,19 @@ - void SetPartId(MMcEvt::ParticleId_t id) - { fPartId = id; } - - void SetEnergy(Float_t Energy) - { fEnergy=Energy; } //Set Energy - - void SetImpact(Float_t Impact) - { fImpact=Impact;} //Set impact parameter + // Setter + void SetPartId(ParticleId_t id) { fPartId = id; } + void SetEnergy(Float_t Energy) { fEnergy=Energy; } //Set Energy + void SetImpact(Float_t Impact) { fImpact=Impact;} //Set impact parameter void SetTelescopeTheta(Float_t Theta) { fTelescopeTheta=Theta; } - void SetTelescopePhi (Float_t Phi) { fTelescopePhi=Phi; } - ClassDef(MMcEvtBasic, 1) //Stores Basic Montecarlo Information of one event + void Fill(ParticleId_t, Float_t, Float_t, Float_t, Float_t); + + // TObject + void Clear(Option_t *opt=NULL); + void Print(Option_t *opt=NULL) const; + + ClassDef(MMcEvtBasic, 2) //Stores Basic Montecarlo Information of one event };