/* ======================================================================== *\ ! ! * ! * 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 // ////////////////////////////////////////////////////////////////////////////// #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; 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 = (MWeight*)pList->FindCreateObj("MWeight"); if (!fWeight) { *fLog << err << dbginf << "MWeight not found... exit." << endl; 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; if (fNewSpecIsPowLaw) weight = C * pow(energy,fNewSlope-fCorsikaSlope); else weight = C * fNewSpectrum->Eval(energy) / pow(energy,fCorsikaSlope); fWeight->SetWeight( weight ); return kTRUE; }