Changeset 2474 for trunk/MagicSoft


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

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r2473 r2474  
    11                                                 -*-*- END OF LINE -*-*-
     2 2003/11/05: Marcos Lopez
     3 
     4   * mmontecarlo/MMcWeightEnergySpecCalc.[h,cc]
     5     - Now, if the new spectrum for the MC showers is a power law, we don't
     6       convert it to a TF1 function.
     7     - Changed the constructor for the case in which the new spectrum is
     8       passed as a TF1 function. Now we pass the TF1 object by reference.
     9     - Thanks to the suggestions of T. Bretz, added three more constructors to
     10       give the possibility of passing the shape of the new spectrum in other
     11       different ways. Now, if the new spectrum that you want for the MC
     12       showers is different from a power law, you can especify its shape either
     13       with a TF1 function, with a string (char*), or with a general C++
     14       function defined by your own.
     15     - In function Reinit(): added a sanity check to prevent from dividing
     16       by zero.
     17     - In PreProcess(): removed an unnecessary sentence.
     18     - Fixed a compiling error which appeared under gcc 3.3
     19       
     20   * macros/weights.C
     21     - addapted to show the new features introduced.
     22
     23
     24
    225 2003/11/05: Thomas Bretz
    326 

  • trunk/MagicSoft/Mars/macros/weights.C

    r2449 r2474  
    2323\* ======================================================================== */
    2424
     25
    2526//////////////////////////////////////////////////////////////////////////////
    2627//                                                                          //
    … …  
    3132//////////////////////////////////////////////////////////////////////////////
    3233
     34
     35// --------------------------------------------------------------------------
     36//
     37//
     38Double_t myfunction(Double_t *x, Double_t *par)
     39{
     40  Double_t xx = x[0];
     41
     42  return pow(xx,-2)*exp(-xx/20); 
     43}
     44
     45
     46
     47// --------------------------------------------------------------------------
     48//
     49//
    3350void weights(TString filename="/up1/data/Magic-MC/CameraAll/Gammas/zbin0/Gamma_zbin0_0_7_1000to1009_w0-4:4:2.root")
    3451{
    … …  
    6380    reader.EnableBranch("fImpact");
    6481
    65     // ------------------
     82
     83    // -------------------------------------------------------------
    6684    //
    6785    // Option 1. Just change the slope of the MC power law spectrum
    6886    //
    69     //MMcWeightEnergySpecCalc wcalc(-2.0);
     87    //MMcWeightEnergySpecCalc wcalc(-2.0);                //<-- Uncomment me
     88
    7089    //
    71     // Option 2. A completely differente specturm
     90    // Option 2. A completely differente specturm pass as a TF1 function
    7291    //           e.g. spectrum with exponential cutoff
    7392    //
    74     TF1* spec = new TF1("spectrum","pow(x,[0])*exp(-x/[1])");
    75     spec->SetParameter(0,-2.0); // Spectral index
    76     spec->SetParameter(1,50); // Spectral index
    77     MMcWeightEnergySpecCalc wcalc(spec);
     93    //TF1 spec("spectrum","pow(x,[0])*exp(-x/[1])");      //<-- Uncomment me
     94    //spec->SetParameter(0,-2.0);                         //<-- Uncomment me
     95    //spec->SetParameter(1,50);                           //<-- Uncomment me
     96    //MMcWeightEnergySpecCalc wcalc(spec);                //<-- Uncomment me
     97 
     98    //
     99    // Option 3. A completely differente specturm pass as a cahr*
     100    //           
     101    //char* func = "pow(x,-2)";                           //<-- Uncomment me
     102    //MMcWeightEnergySpecCalc wcalc(func);                //<-- Uncomment me
     103
     104    //
     105    // Option 4. A completely differente specturm pass as a c++ function
     106    //     
     107    MMcWeightEnergySpecCalc wcalc((void*)myfunction);   //<-- Uncomment me
     108    //
     109    //-------------------------------------------------------------
     110
    78111
    79112    MFillH hfill(&h1,"MMcEvt");
    80113    MFillH hfill2(&h2,"MMcEvt");
    81114    hfill2.SetWeight("MWeight");
    82     //
    83     //------------------
    84115
    85116    tasklist.AddToList(&reader);
    … …  
    115146    hist1->DrawClone();
    116147    hist2->DrawClone("same");   
     148
    117149}
    118150

  • 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.