Changeset 2474 for trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.cc

Timestamp:

11/05/03 14:47:52 (21 years ago)

Author:

marcos

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.cc (modified) (12 diffs)

trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.cc

@@ -25 +25 @@
 //////////////////////////////////////////////////////////////////////////////
 //
-//  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 +55 @@
 //  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 +105 @@
 
 
+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 +126 @@
 {
     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 +223 @@
     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 +262 @@
 
 
-    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 +279 @@
     *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 +298 @@
     // 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 +320 @@
         // 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 +334 @@
 
 
+
 // ----------------------------------------------------------------------------
 //
@@ -247 +342 @@
     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 +374 @@
 
 
+
+
+