Ignore:
Timestamp:
11/05/03 14:47:52 (21 years ago)
Author:
marcos
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mmontecarlo
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.cc

    r2448 r2474  
    2525//////////////////////////////////////////////////////////////////////////////
    2626//
    27 //  MWeightEnergySlopeCalc
     27//  MMcWeightEnergySlopeCalc
    2828//               
    2929//  Change the spectrum of the MC showers simulated with Corsika (a power law)
    3030//  to a new one, which can be either, again a power law but with a different
    31 //  spectral index, or a generalizeed spectra (given as a TF1 object)
     31//  spectral index, or a generalizeed spectrum. The new spectrum can be
     32//  pass to this class in different ways:
     33//    1. Is the new spectrum will be a power law, just introduce the slope
     34//       of this power law.
     35//    2. Is the new spectrum will have a general shape, different options are
     36//       available:
     37//       a) The new spectrum is pass as a TF1 function
     38//       b) The new spectrum is pass as a char*
     39//       c) The new spectrum is pass as a "interpreted function", i.e., a
     40//          function defined inside a ROOT macro, which will be invoked by the
     41//          ROOT Cint itself. This is the case when we use ROOT macros.
     42//       d) The new spectrum is pass as a "real function", i.e., a
     43//          function defined inside normal c++ file.
    3244//
    3345//  Method:
    3446//  ------
    3547//                                 
    36 //  -Corsika spectrun: dN/dE = A * E^(-a)
     48//  -Corsika spectrun: dN/dE = A * E^(a) 
    3749//    with a = fCorsikaSlope, and A = N/integral{E*de} from ELowLim to EUppLim
    3850//
    … …  
    4355//  a weight to each event, given by:
    4456//
    45 //     W(E) = B/A * g(E)/E^(-a)
    46 //
    47 //  (The factor B/A is introduced in order both the original and new spectrum
    48 //  have the same area (i.e. number of showers))
     57//     W(E) = B/A * g(E)/E^(a)
     58//
     59//  In the case the new spectrum is simply a power law: dN/dE = B * E^(b), we
     60//  have:
     61//
     62//     W(E) = B/A * E^(b-a)
     63//
     64//  (The factor B/A is used in order both the original and new spectrum have
     65//   the same area (i.e. in order they represent the same number of showers))
    4966//
    5067//  Note:
    5168//  ------
    52 //   -In the calculations, the new spectrum is treated internally as a TF1
    53 //    object, to give to the user the option of applying a general spectrum.
    54 //    Is the new spectrum is just a power law (i.e. the user only specify the
    55 //    slope),it is converted to a TF1 object for consistency.
    56 //   -With the original corsika spectrum, as it is always a power law, not TF1
    57 //    object is used to represent the spectrum.
    58 //
    59 //
     69//   -If the the new spectrum is just a power law (i.e. the user only specify
     70//     the slope), the needed calculations (such as the integral of the
     71//     spectrum) are done analytically. But if the new spectrum is given as a
     72//     TF1 object, the calculations is done numerically.
     73//
     74//  ToDo:
     75//  -----
     76//   -Give to the user also the possibility to specify the integral of the
     77//    spectrum as another TF1 object (or C++ function)
     78//
     79//
    6080//  Input Containers:                         
    6181//   MMcEvt, MMcRunHeader, MMcCorsikaRunHeader
    … …  
    85105
    86106
     107void MMcWeightEnergySpecCalc::Init(const char *name, const char *title)
     108{
     109
     110    fName  = name  ? name  : "MMcWeightEnergySpecCalc";
     111    fTitle = title ? title : "Task to calculate weights to change the energy spectrum";
     112   
     113    AddToBranchList("MMcEvt.fEnergy");
     114
     115    fAllEvtsTriggered         =  kFALSE;
     116    fTotalNumSimulatedShowers =  0;
     117}
     118
     119
    87120// ---------------------------------------------------------------------------
    88121//
     122// Constructor. The new spectrum will be just a power law.
    89123//
    90124MMcWeightEnergySpecCalc::MMcWeightEnergySpecCalc(Float_t slope,
    … …  
    92126{
    93127    fNewSpecIsPowLaw = kTRUE;
    94 
    95128    fNewSlope = slope;
    96129
    97     fNewSpectrum = new TF1("NewSpectrum","pow(x,[0])");
    98     fNewSpectrum->SetParameter(0,fNewSlope);
    99 
    100 
    101     fName  = name  ? name  : "MMcWeightEnergySpecCalc";
    102     fTitle = title ? title : "Task to calculate weights to change the energy spectrum";
    103 
    104     AddToBranchList("MMcEvt.fEnergy");
    105 
    106     fAllEvtsTriggered         =  kFALSE;
    107     fTotalNumSimulatedShowers =  0;
    108 }
    109 
    110 MMcWeightEnergySpecCalc::MMcWeightEnergySpecCalc(TF1* spectrum,
     130    fNewSpectrum = NULL;
     131
     132    Init(name,title);
     133}
     134
     135//
     136// Constructor. The new spectrum will have a general shape, given by the user
     137// as a TF1 function.
     138//
     139MMcWeightEnergySpecCalc::MMcWeightEnergySpecCalc(const TF1& spectrum,
    111140                                          const char *name, const char *title)
    112141{
    113142    fNewSpecIsPowLaw = kFALSE;
    114 
    115     fNewSpectrum = spectrum;
    116 
    117     fName  = name  ? name  : "MMcWeightEnergySpecCalc";
    118     fTitle = title ? title : "Task to calculate weights to change the energy spectrum";
    119 
    120     AddToBranchList("MMcEvt.fEnergy");
    121 
    122     fAllEvtsTriggered         =  kFALSE;
    123     fTotalNumSimulatedShowers =  0;
    124 }
    125 
    126 
     143    fNewSpectrum = (TF1*)spectrum.Clone();
     144
     145    Init(name,title);
     146}
     147
     148//
     149// As before, but the function which represent the new spectrum is given as
     150// a char* . Starting from it, we build a TF1 function
     151//
     152MMcWeightEnergySpecCalc::MMcWeightEnergySpecCalc(const char* spectrum,
     153                                          const char *name, const char *title)
     154{
     155    fNewSpecIsPowLaw = kFALSE;
     156    fNewSpectrum = new TF1("NewSpectrum",spectrum);
     157
     158    Init(name,title);
     159}
     160
     161//
     162// As before, but the new spectrum is given as a intrepreted C++ function.
     163// Starting from it we build a TF1 function.
     164// This constructor is called for interpreted functions by CINT, i.e., when
     165// the functions are declared inside a ROOT macro.
     166//
     167// NOTE: you muss do a casting to (void*) of the function that you pass to this
     168//       constructor before invoking it in a macro, e.g.
     169//
     170//       Double_t myfunction(Double_t *x, Double_t *par)
     171//         {
     172//           ...
     173//         }
     174//
     175//        MMcWeightEnergySpecCalc wcalc((void*)myfunction);
     176//
     177//        tasklist.AddToList(&wcalc);
     178//
     179//        otherwise ROOT will invoke the constructor McWeightEnergySpecCalc(
     180//         const char* spectrum, const char *name, const char *title)
     181//
     182MMcWeightEnergySpecCalc::MMcWeightEnergySpecCalc(void* function,
     183                                          const char *name, const char *title)
     184{
     185    fNewSpecIsPowLaw = kFALSE;
     186    fNewSpectrum = new TF1("NewSpectrum",function,0,1,1);
     187
     188    Init(name,title);
     189}
     190
     191//
     192// As before, but this is the constructor for real functions, i.e. it is called
     193// when invoked with the normal C++ compiler, i.e. not inside a ROOT macro.
     194//
     195MMcWeightEnergySpecCalc::MMcWeightEnergySpecCalc(Double_t (*function)(Double_t*x, Double_t* par), const Int_t npar,  const char *name, const char *title)
     196{
     197    fNewSpecIsPowLaw = kFALSE;
     198    fNewSpectrum = new TF1("NewSpectrum",function,0,1,1);
     199
     200    Init(name,title);
     201}
     202
     203
     204
     205// ----------------------------------------------------------------------------
     206//
     207// Destructor. Deletes the cloned fNewSpectrum.
     208//
    127209MMcWeightEnergySpecCalc::~MMcWeightEnergySpecCalc()
    128210{
    129   delete fNewSpectrum;
    130 }
     211  if (fNewSpectrum)
     212    delete fNewSpectrum;
     213}
     214
    131215
    132216
    … …  
    139223    fMcEvt = (MMcEvt*)pList->FindObject("MMcEvt");
    140224    if (!fMcEvt)
    141     {
    142         *fLog << err << dbginf << "MMcEvt not found... exit." << endl;
    143         return kFALSE;
    144     }
    145 
    146     fWeight = (MWeight*)pList->FindObject("MWeight");
     225      {
     226        *fLog << err << dbginf << "MMcEvt not found... exit." << endl;
     227        return kFALSE;
     228      }
     229
     230    fWeight = (MWeight*)pList->FindCreateObj("MWeight");
    147231    if (!fWeight)
    148     {
    149       fWeight = (MWeight*)pList->FindCreateObj("MWeight");
    150         if (!fWeight)
    151             return kFALSE;
    152     }
    153 
     232      {
     233        *fLog << err << dbginf << "MWeight not found... exit." << endl;
     234        return kFALSE;
     235      }
     236   
    154237    return kTRUE;
    155238}
    156239
    157240
     241
    158242// ----------------------------------------------------------------------------
    159243//
    160244// Executed each time a new root file is loaded
    161 // We will need the fCorsikaSlope, and fE{Upp,Low}Lim to calculate the weights
     245// We will need fCorsikaSlope and fE{Upp,Low}Lim to calculate the weights
    162246//
    163247Bool_t MMcWeightEnergySpecCalc::ReInit(MParList *plist)
    … …  
    178262
    179263
    180     fCorsikaSlope = corrunheader->GetSlopeSpec();
    181     fELowLim      = corrunheader->GetELowLim();
    182     fEUppLim      = corrunheader->GetEUppLim();
     264    fCorsikaSlope = (Double_t)corrunheader->GetSlopeSpec();
     265    fELowLim      = (Double_t)corrunheader->GetELowLim();
     266    fEUppLim      = (Double_t)corrunheader->GetEUppLim();
    183267
    184268    fTotalNumSimulatedShowers += runheader->GetNumSimulatedShowers();   
    … …  
    195279    *fLog << inf << "Only triggered events avail: " << (fAllEvtsTriggered?"yes":"no") << endl;
    196280
    197    
    198 
    199     //
    200     // Starting from fCorsikaSlope, and fE{Upp,Low}Lim, calculate the integrals
     281
     282
     283    //
     284    // Sanity checks to be sure that we won't divide by zero later on
     285    //
     286    if(fCorsikaSlope == -1. || fNewSlope == -1.)
     287      {
     288        *fLog << err << "The Slope of the power law must be different of -1... exit" << endl;
     289        return kFALSE;
     290      }
     291
     292
     293    //
     294    // Starting from fCorsikaSlope and fE{Upp,Low}Lim, calculate the integrals
    201295    // of both, the original Corsika spectrum and the new one.
    202296    //
    … …  
    204298    // For the Corsika simulated spectrum (just a power law), we have:
    205299    //
    206     Double_t a = fCorsikaSlope;
    207     fCorSpecInt = ( pow(fEUppLim,1+a) - pow(fELowLim,1+a) ) / ( 1+a );
     300    fCorSpecInt = ( pow(fEUppLim,1+fCorsikaSlope) - pow(fELowLim,1+fCorsikaSlope) ) / ( 1+fCorsikaSlope );
    208301 
     302 
    209303    //
    210304    // For the new spectrum:
    211305    //
    212     fNewSpectrum->SetRange(fELowLim, fEUppLim);
    213    
    214306    if (fNewSpecIsPowLaw)     // just the integral of a power law
    215307      {
    216         Double_t b = fNewSlope;
    217         fNewSpecInt = ( pow(fEUppLim,1+b) - pow(fELowLim,1+b) ) / ( 1+b );
     308        fNewSpecInt = ( pow(fEUppLim,1+fNewSlope) - pow(fELowLim,1+fNewSlope) )/ ( 1+fNewSlope );
    218309      }
    219310    else
    220       { // In this case we have to integrate the new spectrum numerically. We
     311      {
     312        fNewSpectrum->SetRange(fELowLim, fEUppLim);
     313
     314        // In this case we have to integrate the new spectrum numerically. We
    221315        // could do simply fNewSpectrum->Integral(fELowLim,fEUppLim), but ROOT
    222316        // fails integrating up to fEUppLim for a sharp cutoff spectrum
    … …  
    226320        // fNewSpectrum->Integral(fELowLim,fEUppLim) (although not perfectlly)
    227321        //
    228         fNewSpectrum->SetNpx(1e4);
     322        fNewSpectrum->SetNpx(1000);
    229323        TGraph gr(fNewSpectrum,"i");
    230324        Int_t Npx = gr.GetN();
    231325        Double_t* y = gr.GetY();
    232326
    233         Double_t integral = y[Npx-1];
     327        const Double_t integral = y[Npx-1];
    234328        fNewSpecInt = integral;
    235329    }
    … …  
    240334
    241335
     336
    242337// ----------------------------------------------------------------------------
    243338//
    … …  
    247342    const Double_t energy = fMcEvt->GetEnergy();
    248343
    249     Double_t weight = fCorSpecInt/fNewSpecInt * fNewSpectrum->Eval(energy) / pow(energy,fCorsikaSlope);
     344    const Double_t C = fCorSpecInt / fNewSpecInt;
     345    Double_t weight;
     346
     347
     348    if (fNewSpecIsPowLaw)   
     349      weight = C * pow(energy,fNewSlope-fCorsikaSlope);
     350    else
     351      weight = C * fNewSpectrum->Eval(energy) / pow(energy,fCorsikaSlope);
    250352     
     353
    251354    fWeight->SetWeight( weight );
    252355
    … …  
    271374
    272375
     376
     377
     378

  • trunk/MagicSoft/Mars/mmontecarlo/MMcWeightEnergySpecCalc.h

    r2448 r2474  
    77
    88
    9 
    109class MParList;
    1110class MMcEvt;
    … …  
    1312class TF1;
    1413
     14
    1515class MMcWeightEnergySpecCalc : public MTask
    1616{
    17 private:
    18     const MMcEvt  *fMcEvt;
    19     MWeight *fWeight;
    2017
    21     TF1* fNewSpectrum;       // Function with the new spectrum
     18 private:
    2219
    23     Bool_t fNewSpecIsPowLaw; // Tells whether the new spectrum is introduce as
    24                              //a TF1 object or not
     20    const MMcEvt *fMcEvt;
     21    MWeight  *fWeight;
    2522
    26     UInt_t fTotalNumSimulatedShowers;
    27     Bool_t fAllEvtsTriggered;
    28     Float_t fCorsikaSlope; // Slope of energy spectrum generated with Corsika
    29     Float_t fNewSlope;     // Slope of the new spectrum (if it is a power law)
    30     Float_t fELowLim;
    31     Float_t fEUppLim;      // Limits of energy range for generation
    32     Double_t fCorSpecInt;  // Value of the integrals of the Corsika and new
    33     Double_t fNewSpecInt;  //spectrum respectively, between fElow and fUppLim
     23    TF1*     fNewSpectrum;     // Function with the new spectrum
    3424
     25    Bool_t   fNewSpecIsPowLaw; // Tells whether the new spectrum is introduced
     26                               //as a TF1 object or not
     27    Double_t fCorsikaSlope; // Slope of energy spectrum generated with Corsika
     28    Double_t fNewSlope;     // Slope of the new spectrum (if it is a power law)
     29    Double_t fELowLim;     
     30    Double_t fEUppLim;      // Limits of energy range for generation
     31    Double_t fCorSpecInt;   // Value of the integrals of the Corsika and new
     32    Double_t fNewSpecInt;   //spectrum respectively, between fElow and fUppLim
    3533
     34    UInt_t   fTotalNumSimulatedShowers;
     35    Bool_t   fAllEvtsTriggered;
    3636
     37    void   Init(const char *name, const char *title);
    3738    Bool_t ReInit(MParList *plist);
    3839
    … …  
    4142   
    4243
    43 public:
     44 public:
     45
    4446    MMcWeightEnergySpecCalc(const Float_t slope,
    45                           const char *name=NULL, const char *title=NULL);
     47                            const char *name=NULL, const char *title=NULL);
    4648
    47     MMcWeightEnergySpecCalc(TF1* spectrum,
    48                           const char *name=NULL, const char *title=NULL);
     49    MMcWeightEnergySpecCalc(const TF1& spectrum,
     50                            const char *name=NULL, const char *title=NULL);
     51
     52    MMcWeightEnergySpecCalc(const char* spectrum,
     53                            const char *name=NULL, const char *title=NULL);
     54   
     55    MMcWeightEnergySpecCalc(void *function,
     56                            const char *name=NULL, const char *title=NULL);
     57 
     58    MMcWeightEnergySpecCalc(Double_t (*function)(Double_t* x, Double_t* par),
     59              const Int_t npar, const char *name=NULL, const char *title=NULL);
     60   
    4961
    5062    ~MMcWeightEnergySpecCalc();
     63
    5164
    5265    ClassDef(MMcWeightEnergySpecCalc, 0) // Task to convert the spectrum of the MC showers simulated with Corsika to a different one, by using weights
Note: See TracChangeset for help on using the changeset viewer.