Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2473) +++ trunk/MagicSoft/Mars/Changelog (revision 2474) @@ -1,3 +1,26 @@ -*-*- END OF LINE -*-*- + 2003/11/05: Marcos Lopez + + * mmontecarlo/MMcWeightEnergySpecCalc.[h,cc] + - Now, if the new spectrum for the MC showers is a power law, we don't + convert it to a TF1 function. + - Changed the constructor for the case in which the new spectrum is + passed as a TF1 function. Now we pass the TF1 object by reference. + - Thanks to the suggestions of T. Bretz, added three more constructors to + give the possibility of passing the shape of the new spectrum in other + different ways. Now, if the new spectrum that you want for the MC + showers is different from a power law, you can especify its shape either + with a TF1 function, with a string (char*), or with a general C++ + function defined by your own. + - In function Reinit(): added a sanity check to prevent from dividing + by zero. + - In PreProcess(): removed an unnecessary sentence. + - Fixed a compiling error which appeared under gcc 3.3 + + * macros/weights.C + - addapted to show the new features introduced. + + + 2003/11/05: Thomas Bretz Index: trunk/MagicSoft/Mars/macros/weights.C =================================================================== --- trunk/MagicSoft/Mars/macros/weights.C (revision 2473) +++ trunk/MagicSoft/Mars/macros/weights.C (revision 2474) @@ -23,4 +23,5 @@ \* ======================================================================== */ + ////////////////////////////////////////////////////////////////////////////// // // @@ -31,4 +32,20 @@ ////////////////////////////////////////////////////////////////////////////// + +// -------------------------------------------------------------------------- +// +// +Double_t myfunction(Double_t *x, Double_t *par) +{ + Double_t xx = x[0]; + + return pow(xx,-2)*exp(-xx/20); +} + + + +// -------------------------------------------------------------------------- +// +// void weights(TString filename="/up1/data/Magic-MC/CameraAll/Gammas/zbin0/Gamma_zbin0_0_7_1000to1009_w0-4:4:2.root") { @@ -63,23 +80,37 @@ reader.EnableBranch("fImpact"); - // ------------------ + + // ------------------------------------------------------------- // // Option 1. Just change the slope of the MC power law spectrum // - //MMcWeightEnergySpecCalc wcalc(-2.0); + //MMcWeightEnergySpecCalc wcalc(-2.0); //<-- Uncomment me + // - // Option 2. A completely differente specturm + // Option 2. A completely differente specturm pass as a TF1 function // 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); + //TF1 spec("spectrum","pow(x,[0])*exp(-x/[1])"); //<-- Uncomment me + //spec->SetParameter(0,-2.0); //<-- Uncomment me + //spec->SetParameter(1,50); //<-- Uncomment me + //MMcWeightEnergySpecCalc wcalc(spec); //<-- Uncomment me + + // + // Option 3. A completely differente specturm pass as a cahr* + // + //char* func = "pow(x,-2)"; //<-- Uncomment me + //MMcWeightEnergySpecCalc wcalc(func); //<-- Uncomment me + + // + // Option 4. A completely differente specturm pass as a c++ function + // + MMcWeightEnergySpecCalc wcalc((void*)myfunction); //<-- Uncomment me + // + //------------------------------------------------------------- + MFillH hfill(&h1,"MMcEvt"); MFillH hfill2(&h2,"MMcEvt"); hfill2.SetWeight("MWeight"); - // - //------------------ tasklist.AddToList(&reader); @@ -115,4 +146,5 @@ hist1->DrawClone(); hist2->DrawClone("same"); + } Index: trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.cc (revision 2473) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.cc (revision 2474) @@ -25,14 +25,26 @@ ////////////////////////////////////////////////////////////////////////////// // -// MWeightEnergySlopeCalc +// 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 spectra (given as a TF1 object) +// 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) +// -Corsika spectrun: dN/dE = A * E^(a) // with a = fCorsikaSlope, and A = N/integral{E*de} from ELowLim to EUppLim // @@ -43,19 +55,27 @@ // 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)) +// 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: // ------ -// -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. -// -// +// -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 @@ -85,6 +105,20 @@ +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, @@ -92,41 +126,91 @@ { 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, + 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 = spectrum; - - fName = name ? name : "MMcWeightEnergySpecCalc"; - fTitle = title ? title : "Task to calculate weights to change the energy spectrum"; - - AddToBranchList("MMcEvt.fEnergy"); - - fAllEvtsTriggered = kFALSE; - fTotalNumSimulatedShowers = 0; -} - - + 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() { - delete fNewSpectrum; -} + if (fNewSpectrum) + delete fNewSpectrum; +} + @@ -139,25 +223,25 @@ fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt"); if (!fMcEvt) - { - *fLog << err << dbginf << "MMcEvt not found... exit." << endl; - return kFALSE; - } - - fWeight = (MWeight*)pList->FindObject("MWeight"); + { + *fLog << err << dbginf << "MMcEvt not found... exit." << endl; + return kFALSE; + } + + fWeight = (MWeight*)pList->FindCreateObj("MWeight"); if (!fWeight) - { - fWeight = (MWeight*)pList->FindCreateObj("MWeight"); - if (!fWeight) - return kFALSE; - } - + { + *fLog << err << dbginf << "MWeight not found... exit." << endl; + 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 +// We will need fCorsikaSlope and fE{Upp,Low}Lim to calculate the weights // Bool_t MMcWeightEnergySpecCalc::ReInit(MParList *plist) @@ -178,7 +262,7 @@ - fCorsikaSlope = corrunheader->GetSlopeSpec(); - fELowLim = corrunheader->GetELowLim(); - fEUppLim = corrunheader->GetEUppLim(); + fCorsikaSlope = (Double_t)corrunheader->GetSlopeSpec(); + fELowLim = (Double_t)corrunheader->GetELowLim(); + fEUppLim = (Double_t)corrunheader->GetEUppLim(); fTotalNumSimulatedShowers += runheader->GetNumSimulatedShowers(); @@ -195,8 +279,18 @@ *fLog << inf << "Only triggered events avail: " << (fAllEvtsTriggered?"yes":"no") << endl; - - - // - // Starting from fCorsikaSlope, and fE{Upp,Low}Lim, calculate the integrals + + + // + // 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. // @@ -204,19 +298,19 @@ // 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 ); + fCorSpecInt = ( pow(fEUppLim,1+fCorsikaSlope) - pow(fELowLim,1+fCorsikaSlope) ) / ( 1+fCorsikaSlope ); + // // 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 ); + fNewSpecInt = ( pow(fEUppLim,1+fNewSlope) - pow(fELowLim,1+fNewSlope) )/ ( 1+fNewSlope ); } else - { // In this case we have to integrate the new spectrum numerically. We + { + 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 @@ -226,10 +320,10 @@ // fNewSpectrum->Integral(fELowLim,fEUppLim) (although not perfectlly) // - fNewSpectrum->SetNpx(1e4); + fNewSpectrum->SetNpx(1000); TGraph gr(fNewSpectrum,"i"); Int_t Npx = gr.GetN(); Double_t* y = gr.GetY(); - Double_t integral = y[Npx-1]; + const Double_t integral = y[Npx-1]; fNewSpecInt = integral; } @@ -240,4 +334,5 @@ + // ---------------------------------------------------------------------------- // @@ -247,6 +342,14 @@ const Double_t energy = fMcEvt->GetEnergy(); - Double_t weight = fCorSpecInt/fNewSpecInt * fNewSpectrum->Eval(energy) / pow(energy,fCorsikaSlope); + const Double_t C = fCorSpecInt / fNewSpecInt; + Double_t weight; + + + if (fNewSpecIsPowLaw) + weight = C * pow(energy,fNewSlope-fCorsikaSlope); + else + weight = C * fNewSpectrum->Eval(energy) / pow(energy,fCorsikaSlope); + fWeight->SetWeight( weight ); @@ -271,2 +374,5 @@ + + + Index: trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.h =================================================================== --- trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.h (revision 2473) +++ trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.h (revision 2474) @@ -7,5 +7,4 @@ - class MParList; class MMcEvt; @@ -13,26 +12,28 @@ class TF1; + class MMcWeightEnergySpecCalc : public MTask { -private: - const MMcEvt *fMcEvt; - MWeight *fWeight; - TF1* fNewSpectrum; // Function with the new spectrum + private: - Bool_t fNewSpecIsPowLaw; // Tells whether the new spectrum is introduce as - //a TF1 object or not + const MMcEvt *fMcEvt; + MWeight *fWeight; - 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 + 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); @@ -41,12 +42,24 @@ -public: + public: + MMcWeightEnergySpecCalc(const Float_t slope, - const char *name=NULL, const char *title=NULL); + const char *name=NULL, const char *title=NULL); - MMcWeightEnergySpecCalc(TF1* spectrum, - 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