Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2447) +++ trunk/MagicSoft/Mars/Changelog (revision 2448) @@ -1,3 +1,26 @@ -*-*- END OF LINE -*-*- + 2003/10/31: Marcos Lopez + + * mhist/MFill.cc: + - Fixed a bug in function PreProcess(MParList *pList). Inside the + conditional sentence "if (!fWeight && !fWeightName.IsNull())", + fWeight never pointed to the object MWeight recoverd from the + parameter list. + + * mhistmc/MHMcEnergyImpact.cc + - In the Fill function, pass the weight to the histogram fHist. + + * mmontecarlo/MMcWeightEnergySpecCalc.[h,cc] + - Added new class for changing the energy spectrum of the showers + simulated with Corsika to a different one, be using weights + + * mmontecarlo/Makefile, MonteCarloLinkDef.h + - Added the new class. + + * macros/weights.C + - Added macro showing how to transform the spectrum of the MC showers. + + + 2003/10/31: Thomas Bretz Index: trunk/MagicSoft/Mars/macros/weights.C =================================================================== --- trunk/MagicSoft/Mars/macros/weights.C (revision 2448) +++ trunk/MagicSoft/Mars/macros/weights.C (revision 2448) @@ -0,0 +1,130 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// // +// This macro shows how to use the class MMcWeightEnergySpecCalc // +// to convert the energy spectrum of the MC showers generated with Corsika, // +// to a different one. // +// // +////////////////////////////////////////////////////////////////////////////// + +void weights(TString filename="/up1/data/Magic-MC/CameraAll/Gammas/zbin0/Gamma_zbin0_0_7_1000to1009_w0-4:4:2.root") +{ + + // + // PartList + // + MParList parlist; + MTaskList tasklist; + + MHMcEnergyImpact h1("h1"); + MHMcEnergyImpact h2("h2"); + parlist.AddToList(&h1); + parlist.AddToList(&h2); + + MBinning binsenergy("BinningEnergy"); + binsenergy.SetEdgesLog(100, 1, 1e5); + parlist.AddToList(&binsenergy); + + MBinning binsimpact("BinningImpact"); + binsimpact.SetEdges(100, 0, 450); + parlist.AddToList(&binsimpact); + + parlist.AddToList(&tasklist); + + + // + // TaskList + // + MReadMarsFile reader("Events", filename); + reader.EnableBranch("fEnergy"); + reader.EnableBranch("fImpact"); + + // ------------------ + // + // Option 1. Just change the slope of the MC power law spectrum + // + MMcWeightEnergySpecCalc wcalc(-2.0); + // + // Option 2. A completely differente specturm + // e.g. spectrum with exponential cutoff + // + TF1* spec = new TF1("spectrum","pow(x,[0])*exp(-x/[1])"); + spec->SetParameter(0,-2.0); // Spectral index + spec->SetParameter(1,50); // Spectral index + MMcWeightEnergySpecCalc wcalc(spec); + + MFillH hfill(&h1,"MMcEvt"); + MFillH hfill2(&h2,"MMcEvt"); + hfill2.SetWeight("MWeight"); + // + //------------------ + + tasklist.AddToList(&reader); + tasklist.AddToList(&wcalc); + tasklist.AddToList(&hfill); + tasklist.AddToList(&hfill2); + + + // + // EventLoop + // + MEvtLoop magic; + magic.SetParList(&parlist); + + if (!magic.Eventloop()) + return; + + tasklist.PrintStatistics(); + parlist.Print(); + + + // + // Draw the Results + // + TCanvas *c = new TCanvas(); + c->SetLogy(); + c->SetLogx(); + + TH1D* hist1 = (h1->GetHist())->ProjectionX(); + TH1D* hist2 = (h2->GetHist())->ProjectionX(); + hist2->SetLineColor(2); + + hist1->DrawClone(); + hist2->DrawClone("same"); +} + + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/mhist/MFillH.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MFillH.cc (revision 2447) +++ trunk/MagicSoft/Mars/mhist/MFillH.cc (revision 2448) @@ -383,6 +383,10 @@ { fWeight = (MWeight*)pList->FindObject(fWeightName, "MWeight"); - *fLog << err << fWeightName << " [MWeight] not found... aborting." << endl; - return kFALSE; + + if (!fWeight) + { + *fLog << err << fWeightName << " [MWeight] not found... aborting." << endl; + return kFALSE; + } } Index: trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.cc (revision 2448) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.cc (revision 2448) @@ -0,0 +1,272 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MWeightEnergySlopeCalc +// +// 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 spectra (given as a TF1 object) +// +// 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) +// +// (The factor B/A is introduced in order both the original and new spectrum +// have the same area (i.e. number of showers)) +// +// Note: +// ------ +// -In the calculations, the new spectrum is treated internally as a TF1 +// object, to give to the user the option of applying a general spectrum. +// Is the new spectrum is just a power law (i.e. the user only specify the +// slope),it is converted to a TF1 object for consistency. +// -With the original corsika spectrum, as it is always a power law, not TF1 +// object is used to represent the spectrum. +// +// +// Input Containers: +// MMcEvt, MMcRunHeader, MMcCorsikaRunHeader +// +// Output Container: +// MWeight +// +////////////////////////////////////////////////////////////////////////////// + +#include "MMcWeightEnergySpecCalc.h" + +#include "MParList.h" +#include "MLog.h" +#include "MLogManip.h" +#include "MMcEvt.hxx" +#include "MMcRunHeader.hxx" +#include "MMcCorsikaRunHeader.h" +#include "MWeight.h" + +#include "TF1.h" +#include "TGraph.h" + +ClassImp(MMcWeightEnergySpecCalc); + +using namespace std; + + + +// --------------------------------------------------------------------------- +// +// +MMcWeightEnergySpecCalc::MMcWeightEnergySpecCalc(Float_t slope, + const char *name, const char *title) +{ + fNewSpecIsPowLaw = kTRUE; + + fNewSlope = slope; + + fNewSpectrum = new TF1("NewSpectrum","pow(x,[0])"); + fNewSpectrum->SetParameter(0,fNewSlope); + + + fName = name ? name : "MMcWeightEnergySpecCalc"; + fTitle = title ? title : "Task to calculate weights to change the energy spectrum"; + + AddToBranchList("MMcEvt.fEnergy"); + + fAllEvtsTriggered = kFALSE; + fTotalNumSimulatedShowers = 0; +} + +MMcWeightEnergySpecCalc::MMcWeightEnergySpecCalc(TF1* spectrum, + const char *name, const char *title) +{ + fNewSpecIsPowLaw = kFALSE; + + fNewSpectrum = spectrum; + + fName = name ? name : "MMcWeightEnergySpecCalc"; + fTitle = title ? title : "Task to calculate weights to change the energy spectrum"; + + AddToBranchList("MMcEvt.fEnergy"); + + fAllEvtsTriggered = kFALSE; + fTotalNumSimulatedShowers = 0; +} + + +MMcWeightEnergySpecCalc::~MMcWeightEnergySpecCalc() +{ + 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 = (MWeight*)pList->FindObject("MWeight"); + if (!fWeight) + { + fWeight = (MWeight*)pList->FindCreateObj("MWeight"); + if (!fWeight) + return kFALSE; + } + + return kTRUE; +} + + +// ---------------------------------------------------------------------------- +// +// Executed each time a new root file is loaded +// We will need the 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 = corrunheader->GetSlopeSpec(); + fELowLim = corrunheader->GetELowLim(); + fEUppLim = 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; + + + + // + // 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: + // + Double_t a = fCorsikaSlope; + fCorSpecInt = ( pow(fEUppLim,1+a) - pow(fELowLim,1+a) ) / ( 1+a ); + + // + // For the new spectrum: + // + fNewSpectrum->SetRange(fELowLim, fEUppLim); + + if (fNewSpecIsPowLaw) // just the integral of a power law + { + Double_t b = fNewSlope; + fNewSpecInt = ( pow(fEUppLim,1+b) - pow(fELowLim,1+b) ) / ( 1+b ); + } + else + { // 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(1e4); + TGraph gr(fNewSpectrum,"i"); + Int_t Npx = gr.GetN(); + Double_t* y = gr.GetY(); + + Double_t integral = y[Npx-1]; + fNewSpecInt = integral; + } + + + return kTRUE; +} + + +// ---------------------------------------------------------------------------- +// +// +Int_t MMcWeightEnergySpecCalc::Process() +{ + const Double_t energy = fMcEvt->GetEnergy(); + + Double_t weight = fCorSpecInt/fNewSpecInt * fNewSpectrum->Eval(energy) / pow(energy,fCorsikaSlope); + + fWeight->SetWeight( weight ); + + return kTRUE; +} + + + + + + + + + + + + + + + + + + Index: trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.h =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.h (revision 2448) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.h (revision 2448) @@ -0,0 +1,63 @@ +#ifndef MARS_MMcWeightEnergySpecCalc +#define MARS_MWeigthEnergySlopeCalc + +#ifndef MARS_MTask +#include "MTask.h" +#endif + + + +class MParList; +class MMcEvt; +class MWeight; +class TF1; + +class MMcWeightEnergySpecCalc : public MTask +{ +private: + const MMcEvt *fMcEvt; + MWeight *fWeight; + + TF1* fNewSpectrum; // Function with the new spectrum + + Bool_t fNewSpecIsPowLaw; // Tells whether the new spectrum is introduce as + //a TF1 object or not + + UInt_t fTotalNumSimulatedShowers; + Bool_t fAllEvtsTriggered; + Float_t fCorsikaSlope; // Slope of energy spectrum generated with Corsika + Float_t fNewSlope; // Slope of the new spectrum (if it is a power law) + Float_t fELowLim; + Float_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 + + + + 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(TF1* spectrum, + 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/mmontecarlo/Makefile =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/Makefile (revision 2447) +++ trunk/MagicSoft/Mars/mmontecarlo/Makefile (revision 2448) @@ -32,5 +32,6 @@ MMcTimeGenerate.cc \ MMcTriggerRateCalc.cc \ - MMcEnergyEst.cc + MMcEnergyEst.cc \ + MMcWeightEnergySpecCalc.cc SRCS = $(SRCFILES) Index: trunk/MagicSoft/Mars/mmontecarlo/MonteCarloLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MonteCarloLinkDef.h (revision 2447) +++ trunk/MagicSoft/Mars/mmontecarlo/MonteCarloLinkDef.h (revision 2448) @@ -10,3 +10,5 @@ #pragma link C++ class MMcTriggerRateCalc+; #pragma link C++ class MMcEnergyEst+; +#pragma link C++ class MMcWeightEnergySpecCalc+; + #endif