Changeset 7409


Ignore:
Timestamp:
11/18/05 17:15:30 (19 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft
Files:
16 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r7408 r7409  
    2222   * mhflux/MMcSpectrumWeight.cc:
    2323     - fixed a problem with using X more than once in the formula
     24
     25   * ganymed.cc:
     26     - improved UsageInfo
     27
     28   * macros/optim/optimdisp.C, macros/optim/optimenergy.C:
     29     - added some  more examples
     30
     31   * mbadpixels/MBadPixelsCalc.cc:
     32     - updated authors
     33
     34   * mbase/BaseIncl.h:
     35     - added TArrayD
     36
     37   * mbase/MLogHtml.cc:
     38     - added some includes suggested by Ching-Cheng
     39
     40   * mbase/MMath.[h,cc]:
     41     - added some function for the analytical result of special fits
     42
     43   * mfilter/MFEnergySlope.[h,cc]:
     44     - some updated to make it better fit into Mars
     45     - upadted the user interface
     46
     47   * mhflux/MHEnergyEst.[h,cc]:
     48     - updated to let everything fit what is commonly used. This
     49       was just discussed with Abelardo
     50     - some updates to the plots
     51
     52   * mjobs/MDataSet.cc:
     53     - added some includes suggested by Ching-Cheng
     54     - upadted some output
     55     - remove whitespaces from read TString
     56
     57   * mmc/MMcCorsikaRunHeader.h:
     58     - added new Getter
     59
     60   * mranforest/MRanForest.cc:
     61     - some updates to Grow-output
    2462
    2563

  • trunk/MagicSoft/Mars/NEWS

    r7408 r7409  
    22
    33 *** Version  <cvs>
     4
     5   - general: Updated MMath with new functions to calculate the results of
     6     a exponential, logarithmic and powerlaw fits analytically.
    47
    58   - general: Updated some macros with comments:
    … …  
    4447
    4548   - sponde: Added a plot E^2*dN/dE
     49
     50   - sponde: The energy estimator plot should now show values like
     51     they are commonly used.
     52
     53   - sponde: Now MMcSpectrumWeight also excepts formulas with two X
     54     (a powerlaw with cutoff didn't work before)
    4655
    4756

  • trunk/MagicSoft/Mars/ganymed.cc

    r7380 r7409  
    6161    gLog << "   -q                        Quit when job is finished" << endl;
    6262    gLog << "   -f                        Force overwrite of existing files" << endl;
    63     gLog << "   --n=[n]                   Analysis number" << endl;
     63    gLog << "   --n=number                Analysis number (required if not in dataset file)" << endl;
    6464    gLog << "   --out=path                Path to write the all output to [def=local path]" << endl;
    6565    gLog << "   --ind=path                Path to data/star files [default=datacenter path]" << endl;

  • trunk/MagicSoft/Mars/macros/optim/optimdisp.C

    r7401 r7409  
    4242     opt.AddPreCut(&cuts);
    4343
     44     -------------------- Energy Slope --------------------
     45     MFEnergySlope slope(-2.8);
     46     opt.AddPreCut(&slope);
     47
    4448     -------------------- Other cuts ----------------------
    4549     opt.AddPreCut("MNewImagePar.fLeakage1<0.0001");

  • trunk/MagicSoft/Mars/macros/optim/optimenergy.C

    r7401 r7409  
    3838     opt.AddPreCut(&cuts);
    3939
     40     -------------------- Energy Slope --------------------
     41     MFEnergySlope slope(-2.8);
     42     opt.AddPreCut(&slope);
     43
    4044     -------------------- Other cuts ----------------------
    4145     opt.AddPreCut("MPointingPos.fZd<7");

  • trunk/MagicSoft/Mars/mbadpixels/MBadPixelsCalc.cc

    r7350 r7409  
    1717!
    1818!   Author(s): Thomas Bretz, 02/2004 <mailto:tbretz@astro.uni.wuerzburg.de>
    19 !
    20 !   Copyright: MAGIC Software Development, 2000-2004
     19!   Author(s): Stefan Ruegamer, 08/2005 <mailto:snruegam@astro.uni.wuerzburg.de>
     20!
     21!   Copyright: MAGIC Software Development, 2000-2005
    2122!
    2223!

  • trunk/MagicSoft/Mars/mbase/BaseIncl.h

    r7181 r7409  
    11#ifndef __CINT__
    22
     3#include <TArrayD.h>
    34#include <TVector3.h>
    45
    5 /*
    6 #include <fstream.h>
    7 
    8 #include <TFile.h>
    9 #include <TTree.h>
    10 
    11 #include <TGListBox.h>
    12 */
    13 
    146#endif // __CINT__

  • trunk/MagicSoft/Mars/mbase/MLogHtml.cc

    r6870 r7409  
    3030//////////////////////////////////////////////////////////////////////////////
    3131#include "MLogHtml.h"
     32
     33#include <string.h>  // necessary for Fedora core 2 with kernel 2.6.9-1.667 #1 and gcc 3.4.2
     34#include <errno.h>   // necessary for Fedora core 2 with kernel 2.6.9-1.667 #1 and gcc 3.4.2
    3235
    3336#include <fstream>  // ofstream

  • trunk/MagicSoft/Mars/mbase/MMath.cc

    r7386 r7409  
    3636#endif
    3737
     38#ifndef ROOT_TArrayD
     39#include <TArrayD.h>
     40#endif
    3841// --------------------------------------------------------------------------
    3942//
    … …  
    248251    return InterpolParabLin(vx0, vy, TMath::Cos(x));
    249252}
     253
     254// --------------------------------------------------------------------------
     255//
     256// Analytically calculated result of a least square fit of:
     257//    y = A*e^(B*x)
     258// Equal weights
     259//
     260// It returns TArrayD(2) = { A, B };
     261//
     262// see: http://mathworld.wolfram.com/LeastSquaresFittingExponential.html
     263//
     264TArrayD MMath::LeastSqFitExpW1(Int_t n, Double_t *x, Double_t *y)
     265{
     266    Double_t sumxsqy  = 0;
     267    Double_t sumylny  = 0;
     268    Double_t sumxy    = 0;
     269    Double_t sumy     = 0;
     270    Double_t sumxylny = 0;
     271    for (int i=0; i<n; i++)
     272    {
     273        sumylny  += y[i]*TMath::Log(y[i]);
     274        sumxy    += x[i]*y[i];
     275        sumxsqy  += x[i]*x[i]*y[i];
     276        sumxylny += x[i]*y[i]*TMath::Log(y[i]);
     277        sumy     += y[i];
     278    }
     279
     280    const Double_t dev = sumy*sumxsqy - sumxy*sumxy;
     281
     282    const Double_t a = (sumxsqy*sumylny - sumxy*sumxylny)/dev;
     283    const Double_t b = (sumy*sumxylny - sumxy*sumylny)/dev;
     284
     285    TArrayD rc(2);
     286    rc[0] = TMath::Exp(a);
     287    rc[1] = b;
     288    return rc;
     289}
     290
     291// --------------------------------------------------------------------------
     292//
     293// Analytically calculated result of a least square fit of:
     294//    y = A*e^(B*x)
     295// Greater weights to smaller values
     296//
     297// It returns TArrayD(2) = { A, B };
     298//
     299// see: http://mathworld.wolfram.com/LeastSquaresFittingExponential.html
     300//
     301TArrayD MMath::LeastSqFitExp(Int_t n, Double_t *x, Double_t *y)
     302{
     303    // -------- Greater weights to smaller values ---------
     304    Double_t sumlny  = 0;
     305    Double_t sumxlny = 0;
     306    Double_t sumxsq  = 0;
     307    Double_t sumx    = 0;
     308    for (int i=0; i<n; i++)
     309    {
     310        sumlny  += TMath::Log(y[i]);
     311        sumxlny += x[i]*TMath::Log(y[i]);
     312
     313        sumxsq  += x[i]*x[i];
     314        sumx    += x[i];
     315    }
     316
     317    const Double_t dev = n*sumxsq-sumx*sumx;
     318
     319    const Double_t a = (sumlny*sumxsq - sumx*sumxlny)/dev;
     320    const Double_t b = (n*sumxlny - sumx*sumlny)/dev;
     321
     322    TArrayD rc(2);
     323    rc[0] = TMath::Exp(a);
     324    rc[1] = b;
     325    return rc;
     326}
     327
     328// --------------------------------------------------------------------------
     329//
     330// Analytically calculated result of a least square fit of:
     331//    y = A+B*ln(x)
     332//
     333// It returns TArrayD(2) = { A, B };
     334//
     335// see: http://mathworld.wolfram.com/LeastSquaresFittingLogarithmic.html
     336//
     337TArrayD LeastSqFitLog(Int_t n, Double_t *x, Double_t *y)
     338{
     339    Double_t sumylnx  = 0;
     340    Double_t sumy     = 0;
     341    Double_t sumlnx   = 0;
     342    Double_t sumlnxsq = 0;
     343    for (int i=0; i<n; i++)
     344    {
     345        sumylnx  += y[i]*TMath::Log(x[i]);
     346        sumy     += y[i];
     347        sumlnx   += TMath::Log(x[i]);
     348        sumlnxsq += TMath::Log(x[i])*TMath::Log(x[i]);
     349    }
     350
     351    const Double_t b = (n*sumylnx-sumy*sumlnx)/(n*sumlnxsq-sumlnx*sumlnx);
     352    const Double_t a = (sumy-b*sumlnx)/n;
     353
     354    TArrayD rc(2);
     355    rc[0] = a;
     356    rc[1] = b;
     357    return rc;
     358}
     359
     360// --------------------------------------------------------------------------
     361//
     362// Analytically calculated result of a least square fit of:
     363//    y = A*x^B
     364//
     365// It returns TArrayD(2) = { A, B };
     366//
     367// see: http://mathworld.wolfram.com/LeastSquaresFittingPowerLaw.html
     368//
     369TArrayD LeastSqFitPowerLaw(Int_t n, Double_t *x, Double_t *y)
     370{
     371    Double_t sumlnxlny  = 0;
     372    Double_t sumlnx   = 0;
     373    Double_t sumlny    = 0;
     374    Double_t sumlnxsq   = 0;
     375    for (int i=0; i<n; i++)
     376    {
     377        sumlnxlny  += TMath::Log(x[i])*TMath::Log(y[i]);
     378        sumlnx     += TMath::Log(x[i]);
     379        sumlny     += TMath::Log(y[i]);
     380        sumlnxsq   += TMath::Log(x[i])*TMath::Log(x[i]);
     381    }
     382
     383    const Double_t b = (n*sumlnxlny-sumlnx*sumlny)/(n*sumlnxsq-sumlnx*sumlnx);
     384    const Double_t a = (sumlny-b*sumlnx)/n;
     385
     386    TArrayD rc(2);
     387    rc[0] = TMath::Exp(a);
     388    rc[1] = b;
     389    return rc;
     390}

  • trunk/MagicSoft/Mars/mbase/MMath.h

    r7384 r7409  
    77
    88class TVector3;
     9class TArrayD;
    910
    1011namespace MMath
    … …  
    3132    Double_t InterpolParabCos(const TVector3 &vx, const TVector3 &vy, Double_t x);
    3233
     34    TArrayD LeastSqFitExpW1(Int_t n, Double_t *x, Double_t *y);
     35    TArrayD LeastSqFitExp(Int_t n, Double_t *x, Double_t *y);
     36    TArrayD LeastSqFitLog(Int_t n, Double_t *x, Double_t *y);
     37    TArrayD LeastSqFitPowerLaw(Int_t n, Double_t *x, Double_t *y);
     38
    3339    inline Double_t Sgn(Double_t d) { return d<0 ? -1 : 1; }
    3440}

  • trunk/MagicSoft/Mars/mfilter/MFEnergySlope.cc

    r2206 r7409  
    1818!   Author(s): Antonio Stamerra  02/2003 <mailto:antonio.stamerra@pi.infn.it>
    1919!
    20 !   Copyright: MAGIC Software Development, 2000-2003
     20!   Copyright: MAGIC Software Development, 2000-2005
    2121!
    2222!
    … …  
    6565// --------------------------------------------------------------------------
    6666//
     67//  Default Constructor
     68//
     69MFEnergySlope::MFEnergySlope(const char *name, const char *title):
     70    fNewSlope(-1), fMcMinEnergy(-1.), fMcMaxEnergy(-1.)
     71{
     72    fName  = name  ? name  : "MFEnergySlope";
     73    fTitle = title ? title : "Filter to select energy with given slope";
     74}
     75
     76// --------------------------------------------------------------------------
     77//
    6778//     Constructor
    6879//
    69 MFEnergySlope::MFEnergySlope(const char *name, const char *title):
    70   fNumSelectedEvts(0), fNewSlope(-1), fMcMinEnergy(-1.), fMcMaxEnergy(-1.)
    71 {
    72   //    fContName = cname;
    73   fName  = name  ? name  : "MFEnergySlope";
    74   fTitle = title ? title : "Filter to select energy with given slope";
     80MFEnergySlope::MFEnergySlope(Float_t slope, const char *name, const char *title):
     81    fNewSlope(TMath::Abs(slope)), fMcMinEnergy(-1.), fMcMaxEnergy(-1.)
     82{
     83    fName  = name  ? name  : "MFEnergySlope";
     84    fTitle = title ? title : "Filter to select energy with given slope";
    7585}
    7686
    … …  
    8494Int_t MFEnergySlope::PreProcess(MParList *pList)
    8595{
    86  
     96    if (fNewSlope<0)
     97    {
     98        *fLog << err << "New slope still undefined... aborting." << endl;
     99        return kFALSE;
     100    }
     101
     102    fEvt = (MMcEvt*)pList->FindObject("MMcEvt");
     103    if (!fEvt)
     104    {
     105        *fLog << err << "MMcEvt not found... aborting." << endl;
     106        return kFALSE;
     107    }
     108
     109    return kTRUE;
     110}
     111
     112// --------------------------------------------------------------------------
     113//
     114//   Preprocess
     115// 
     116//  MC slope and min/max energy are read
     117//  Normalization factor is computed
     118//
     119Bool_t MFEnergySlope::ReInit(MParList *pList)
     120{
    87121    MMcCorsikaRunHeader *runheader = (MMcCorsikaRunHeader*)pList->FindObject("MMcCorsikaRunHeader");
    88 
    89122    if (!runheader)
    90       {
    91         *fLog << err << dbginf << fName << " [MMcCorsikaRunHeader] not found... aborting." << endl;
    92         return kFALSE;
    93       }
     123    {
     124        *fLog << err << "MMcCorsikaRunHeader not found... aborting." << endl;
     125        return kFALSE;
     126    }
     127
    94128    //
    95     // Read info from the MC sample (it must be generated with 
     129    // Read info from the MC sample (it must be generated with
    96130    //   reflector ver.0.6 and camera ver. 0.6)
    97131    //
    98     fMcSlope = runheader->GetSlopeSpec();   
     132    fMcSlope = runheader->GetSlopeSpec();
    99133    if (fMcMinEnergy<0)
    100       fMcMinEnergy = runheader->GetELowLim();
     134        fMcMinEnergy = runheader->GetELowLim();
    101135    if (fMcMinEnergy<0)
    102       fMcMaxEnergy = runheader->GetEUppLim();
    103 
    104     *fLog << inf;
    105     *fLog << "MFEnergySlope::PreProcess: fetched MC info:" << endl;
    106     *fLog << "  E Slope:     " << fMcSlope << endl;
    107     *fLog << "  E Min:       " << fMcMinEnergy << endl;
    108     *fLog << "  E Max:       " << fMcMaxEnergy << endl;
    109     *fLog << "  New E Slope: " << fNewSlope << endl;
    110    
     136        fMcMaxEnergy = runheader->GetEUppLim();
     137
    111138    // Slope is used with positive values in the code
    112139    if (fNewSlope < 0)
    113       fNewSlope *= -1;
     140        fNewSlope *= -1;
    114141    if (fMcSlope < 0)
    115       fMcSlope *= -1;
    116 
    117 
    118   // Set normalization on energy 
    119     fN0 = pow(fNewSlope>fMcSlope?fMcMinEnergy:fMcMaxEnergy,fNewSlope-fMcSlope);
    120 
    121   *fLog << "Normalization factor:" <<fN0 << endl;
    122 
    123   //---
    124     fEvt = (MMcEvt*)pList->FindObject("MMcEvt");
    125     if (!fEvt)
    126       {
    127         *fLog << err << dbginf << fName << " [MMcEvt] not found... aborting." << endl;
    128         return kFALSE;
    129       }
     142        fMcSlope *= -1;
     143
     144    // Calculate normalization factor
     145    fN0 = TMath::Power(fNewSlope>fMcSlope ? fMcMinEnergy : fMcMaxEnergy, fNewSlope-fMcSlope);
     146
     147    // Print some infos
     148    *fLog << inf;
     149    *fLog << "Fetched MC info:" << endl;
     150    *fLog << "  Change E Slope " << fMcSlope << " (" << fMcMinEnergy << " < E < ";
     151    *fLog << fMcMaxEnergy << ") to " << fNewSlope << endl;
     152    *fLog << "  Norm factor: " << fN0 << endl;
    130153
    131154    return kTRUE;
    … …  
    144167Int_t MFEnergySlope::Process()
    145168{
    146   fResult = kTRUE;
    147 
    148   // Energy slopes are the same: skip it
    149   if (fNewSlope == fMcSlope)
     169    fResult = kTRUE;
     170
     171    // Energy slopes are the same: skip it
     172    if (fNewSlope == fMcSlope)
     173        return kTRUE;
     174
     175    //  The value of the normalized spectrum is compared with a
     176    //   random value in [0,1];
     177    //   if the latter is higher the event is accepted
     178    const Float_t energy = fEvt->GetEnergy();
     179
     180    /*
     181     //
     182     // If energy bounds different from MC ones are selected, then
     183     // events outside these bounds are rejected, as their slope has
     184     // not been changed.
     185     //
     186     if (energy > fMcMaxEnergy || energy < fMcMinEnergy)
     187     {
     188        fResult = kFALSE;
     189        return kTRUE;
     190     }
     191     */
     192
     193    const Float_t Nexp = fN0 * pow(energy,fMcSlope-fNewSlope);
     194    const Float_t Nrnd = gRandom->Uniform();
     195
     196    fResult = Nexp > Nrnd;
     197
     198    //if (fResult)
     199    //    fNumSelectedEvts++;
     200
    150201    return kTRUE;
    151  
    152   //  The value of the normalized spectrum is compared with a
    153   //   random value in [0,1];
    154   //   if the latter is higher the event is accepted
    155   const Float_t energy = fEvt->GetEnergy();
    156 
    157   /*
    158   //
    159   // If energy bounds different from MC ones are selected, then
    160   // events outside these bounds are rejected, as their slope has
    161   // not been changed.
    162   //
    163   if (energy > fMcMaxEnergy || energy < fMcMinEnergy)
    164     {
    165       fResult = kFALSE;
    166       return kTRUE;
    167     }
    168   */
    169 
    170   const Float_t Nexp = fN0 * pow(energy,fMcSlope-fNewSlope);
    171   const Float_t Nrnd = gRandom->Uniform();
    172 
    173   fResult = Nexp > Nrnd;
    174 
    175   if (!fResult)
    176       return kTRUE;
    177 
    178   fNumSelectedEvts++;
    179   return kTRUE;
    180 }
    181 
     202}
     203
     204
     205// --------------------------------------------------------------------------
     206//
     207//   MFEnergySlope.NewSlope: -2.8
     208//
     209Int_t MFEnergySlope::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
     210{
     211    Bool_t rc = kFALSE;
     212    if (IsEnvDefined(env, prefix, "NewSlope", print))
     213    {
     214        rc = kTRUE;
     215        SetNewSlope(GetEnvValue(env, prefix, "NewSlope", fNewSlope));
     216    }
     217    return rc;
     218}

  • trunk/MagicSoft/Mars/mfilter/MFEnergySlope.h

    r2206 r7409  
    11#ifndef MARS_MFEnergySlope
    22#define MARS_MFEnergySlope
    3 /////////////////////////////////////////////////////////////////////////////
    4 //                                                                         //
    5 // MFEnergySlope                                                           //
    6 //                                                                         //
    7 /////////////////////////////////////////////////////////////////////////////
    83
    94#ifndef MARS_MFilter
    … …  
    116#endif
    127
     8class MMcEvt;
    139class MParList;
    14 class MMcEvt;
    1510class MMcCorsikaRunHeader;
    1611
    … …  
    1813{
    1914private:
    20     Int_t fNumSelectedEvts; // counter for number of selected events
     15    //Int_t fNumSelectedEvts; // counter for number of selected events
    2116
    22     MMcEvt *fEvt;           // Events used to determin energy slope
     17    MMcEvt *fEvt;           //! Events used to determin energy slope
    2318
    24     Bool_t fResult;         // Result returned by IsExpressionTrue
    2519    Float_t fNewSlope;      // New slope set by user
     20    Float_t fMcSlope;       //! Original energy slope from MC data
    2621
    27     Float_t fMcSlope;       // Original energy slope from MC data
    28     Float_t fMcMinEnergy;   // Starting energy of MC data
    29     Float_t fMcMaxEnergy;   // Ending energy of MC data
     22    Float_t fMcMinEnergy;   //! Starting energy of MC data
     23    Float_t fMcMaxEnergy;   //! Ending energy of MC data
    3024
    31     Float_t fN0;            // Normalization factor
     25    Float_t fN0;            //! Normalization factor
    3226
    33     Int_t PreProcess(MParList *pList);
    34     Int_t Process();
     27    Bool_t  fResult;        //! Result returned by IsExpressionTrue
     28
     29    // MTask
     30    Int_t  PreProcess(MParList *pList);
     31    Bool_t ReInit(MParList *pList);
     32    Int_t  Process();
     33
     34    // MFilter
     35    Bool_t IsExpressionTrue() const { return fResult; }
    3536
    3637public:
    3738    MFEnergySlope(const char *name=NULL, const char *title=NULL);
     39    MFEnergySlope(Float_t slope, const char *name=NULL, const char *title=NULL);
    3840
    39     Bool_t IsExpressionTrue() const { return fResult; }
     41    // Setter
     42    void SetNewSlope(Float_t f)    { fNewSlope = TMath::Abs(f); }
     43    void SetMcSlope(Float_t f)     { fMcSlope  = TMath::Abs(f); }
    4044
    41     // Slope is used with positive values in the code
    42     void SetNewSlope(Float_t f) {fNewSlope = TMath::Abs(f);}
    43     void SetMcSlope(Float_t f) {fMcSlope = TMath::Abs(f);}
     45    void SetMcMinEnergy(Float_t f) { fMcMinEnergy = f; }
     46    void SetMcMaxEnergy(Float_t f) { fMcMaxEnergy = f; }
    4447
    45     void SetMcMinEnergy(Float_t f) {fMcMinEnergy = f;}
    46     void SetMcMaxEnergy(Float_t f) {fMcMaxEnergy = f;}
     48    // MParContainer
     49    Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
    4750
    4851    ClassDef(MFEnergySlope, 0) // A Filter to select events with a given energy slope

  • trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc

    r7388 r7409  
    3535#include "MHEnergyEst.h"
    3636
     37#include <TF1.h>
    3738#include <TLine.h>
    3839#include <TCanvas.h>
    … …  
    7980    fHResolution.SetDirectory(NULL);
    8081    fHResolution.SetName("EnergyRes");
    81     fHResolution.SetTitle("Histogram in \\Delta lg(E) vs E_{est} and E_{mc}");
     82    fHResolution.SetTitle("Histogram in \\Delta E/E vs E_{est} and E_{mc}");
    8283    fHResolution.SetXTitle("E_{est} [GeV]");
    8384    fHResolution.SetYTitle("E_{mc} [GeV]");
    84     fHResolution.SetZTitle("lg(E_{est}) - lg(E_{mc})");
     85    fHResolution.SetZTitle("E_{est}/E_{mc}-1");
    8586
    8687    fHImpact.SetDirectory(NULL);
    8788    fHImpact.SetName("Impact");
    88     fHImpact.SetTitle("\\Delta lg(E) vs Impact parameter");
     89    fHImpact.SetTitle("\\Delta E/E vs Impact parameter");
    8990    fHImpact.SetXTitle("Impact parameter [m]");
    90     fHImpact.SetYTitle("lg(E_{est}) - lg(E_{mc})");
     91    fHImpact.SetYTitle("E_{est}/E_{mc}-1");
    9192
    9293    fHEnergy.Sumw2();
    … …  
    9596
    9697    MBinning binsi, binse, binst, binsr;
    97     binse.SetEdgesLog(15, 10, 1000000);
     98    binse.SetEdgesLog(25, 10, 1000000);
    9899    binst.SetEdgesASin(51, -0.005, 0.505);
    99100    binsi.SetEdges(10, 0, 400);
    100     binsr.SetEdges(50, -1.3, 1.3);
     101    binsr.SetEdges(75, -1.75, 1.75);
    101102
    102103    SetBinning(&fHEnergy,     &binse, &binse, &binst);
    … …  
    149150    fChisq = 0;
    150151    fBias  = 0;
     152
    151153    fHEnergy.Reset();
    152154    fHImpact.Reset();
    … …  
    155157    return kTRUE;
    156158}
    157 /*
    158 #include <TSpline.h>
    159 Double_t x[] = {33, 75, 153, 343, 745, 1617, 3509, 7175};
    160 Double_t y[] = {2,  24, 302, 333, 132,   53,   11,    1};
    161 Int_t n = 8;
    162 TSpline3 spline("", x, y, n);
    163 */
     159
    164160// --------------------------------------------------------------------------
    165161//
    … …  
    172168    const Double_t imp   = fMatrix ? GetVal(1) : fMcEvt->GetImpact()/100;
    173169    const Double_t theta = fMatrix ? GetVal(2) : fMcEvt->GetTelescopeTheta()*TMath::RadToDeg();
    174     const Double_t res   = log10(eest)-log10(etru);///log10(etru);
    175170    const Double_t resE  = (eest-etru)/etru;
    176171
    … …  
    179174    fHImpact.Fill(imp, resE, w);
    180175
    181     //if (etru>125 && etru<750)
    182     //    return kCONTINUE;
    183 
    184     // lg(N)  = 4.3*lg(E)-6
    185     // lg(N0) = 4.3*2.2-6
    186 
    187     const Double_t n  = 2.;///spline.Eval(etru); //pow(10, -4.3*(log10(etru)-2.2));
    188     if (n<=0)
    189         return kCONTINUE;
    190 
    191     fChisq += log10(etru)>0 ? res*res*n/2   : res*res;
    192     fBias  += log10(etru)>0 ? res*sqrt(n/2) : res;
     176    // For the fit we use a different quantity
     177    //const Double_t res = TMath::Log10(eest/etru);
     178    const Double_t res = eest-etru;
     179
     180    fChisq += res*res;
     181    fBias  += res;
    193182
    194183    return kTRUE;
    195184}
    196 /*
    197 void MHEnergyEst::CalcChisq(Double_t &chisq, Double_t &prob) const
    198 {
    199     TH1D *h1 = (TH1D*)fHEnergy.Project3D("dum2_ex");
    200     TH1D *h2 = (TH1D*)fHEnergy.Project3D("dum2_ey");
    201 
    202     Int_t df = 0;
    203     chisq    = 0;
    204     prob     = 0;
    205 
    206     for (Int_t i=1; i<=h1->GetNbinsX(); i++)
    207     {
    208         if (h1->GetBinContent(i)>0 && h2->GetBinContent(i)>0)
    209         {
    210             const Double_t bin1 = log10(h1->GetBinContent(i));
    211             const Double_t bin2 = log10(h2->GetBinContent(i));
    212 
    213             const Double_t temp = bin1-bin2;
    214 
    215             chisq += temp*temp/(bin1+bin2);
    216             df ++;
    217         }
    218     }
    219  
    220     prob = TMath::Prob(0.5*chisq, Int_t(0.5*df));
    221 
    222     chisq /= df;
    223  
    224     delete h1;
    225     delete h2;
    226 }
    227 */
     185
    228186Bool_t MHEnergyEst::Finalize()
    229187{
    230     //Double_t chisq, prob;
    231     //CalcChisq(chisq, prob);
    232  
    233188    fChisq /= GetNumExecutions();
    234189    fBias  /= GetNumExecutions();
    235190
    236     fResult->SetVal(TMath::Sqrt(fChisq));
    237 
    238 /*
    239     Double_t sigma = TMath::Sqrt(fChisq);//+TMath::Abs(fBias);
    240     fResult->SetVal(sigma);
    241 */
     191    fResult->SetVal(fChisq);
     192
    242193    Print();
    243194
    … …  
    247198void MHEnergyEst::Print(Option_t *o) const
    248199{
    249     const Double_t res  =   TMath::Sqrt(fChisq-fBias*fBias);
    250     const Double_t resp =   TMath::Power(10,  res)-1;
    251     const Double_t resm = 1-TMath::Power(10, -res);
    252 
     200    //    const Double_t resp =   TMath::Power(10,  res)-1;
     201    //    const Double_t resm = 1-TMath::Power(10, -res);
     202
     203    const Double_t res = TMath::Sqrt(fChisq-fBias*fBias);
    253204    if (TMath::IsNaN(res))
    254205    {
    … …  
    257208    }
    258209
    259     *fLog << all;
    260     *fLog << "Mean log10(Energy) Resoltion: +/- " << Form("%4.2f", res) << endl;
    261     *fLog << "Mean log10(Energy) Bias:         "  << Form("%+4.2f", fBias) << endl;
    262     *fLog << "Asymmetric Energy  Resoltion:   +"  << Form("%4.1f%%", resp*100) << " -" << Form("%4.1f%%", resm*100) << endl;
     210    TH1D *h = (TH1D*)fHResolution.ProjectionZ("Resolution");
     211    h->Fit("gaus", "Q0");
     212
     213    TF1 *f = h->GetFunction("gaus");
     214
     215    *fLog << all << "F=" << fChisq << endl;
     216    *fLog << "Results from Histogram:" << endl;
     217    *fLog << " Mean  of Delta E/E: " << Form("%+4.2f%%", 100*h->GetMean()) << endl;
     218    *fLog << " RMS   of Delta E/E: " << Form("%4.2f%%",  100*h->GetRMS()) << endl;
     219    *fLog << "Results from Histogram-Fit:" << endl;
     220    *fLog << " Bias  of Delta E/E: " << Form("%+4.2f%%", 100*f->GetParameter(1)) << endl;
     221    *fLog << " Sigma of Delta E/E: " << Form("%4.2f%%",  100*f->GetParameter(2)) << endl;
     222    *fLog << " Res   of Delta E/E: " << Form("%4.2f%%",  100*TMath::Hypot(f->GetParameter(1), f->GetParameter(2))) << endl;
     223
     224    delete h;
    263225}
    264226
    … …  
    411373    h1->SetBit(kCanDelete);
    412374    h1->SetLineWidth(2);
    413     h1->SetLineColor(kRed);
     375    h1->SetLineColor(kBlue);
    414376    h1->SetFillStyle(4000);
    415377    h1->SetStats(kFALSE);
    416378
    417379    h2->Draw("");
    418     h1->Draw("E3 hist C same");
     380    h1->Draw("E0 hist C same");
    419381//    h1->Draw("Chistsame");
    420382
    … …  
    509471    //h->SetXTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}");
    510472    h->SetYTitle("Counts");
    511     h->SetTitle("Distribution of \\Delta lg(E)");
     473    h->SetTitle("Distribution of \\Delta E/E");
    512474    h->SetDirectory(NULL);
    513475    h->SetBit(kCanDelete);
    … …  
    535497    h = MakePlot(fHResolution, "zy");
    536498    h->SetXTitle("E_{mc} [GeV]");
    537     h->SetYTitle("lg(E_{est}) - lg(E_{mc})");
    538     h->SetTitle("Energy resolution \\Delta lg(E) vs Monte Carlo energy E_{mc}");
     499    h->SetYTitle("(E_{est}/E_{mc}-1");
     500    h->SetTitle("Energy resolution \\Delta E/E vs Monte Carlo Energy E_{mc}");
    539501    h->SetMinimum(-1.3);
    540502    h->SetMaximum(1.3);
    … …  
    543505    h = MakePlot(fHResolution, "zx");
    544506    h->SetXTitle("E_{est} [GeV]");
    545     h->SetYTitle("lg(E_{est}) - lg(E_{mc})");
    546     h->SetTitle("Energy Resolution \\Delta lg(E) vs estimated Energy E_{est}");
     507    h->SetYTitle("(E_{est}/E_{mc}-1");
     508    h->SetTitle("Energy resolution \\Delta E/E vs estimated Energy E_{est}");
    547509    h->SetMinimum(-1.3);
    548510    h->SetMaximum(1.3);

  • trunk/MagicSoft/Mars/mjobs/MDataSet.cc

    r7389 r7409  
    7474#include "MDataSet.h"
    7575
     76#include <string.h>  // necessary for Fedora core 2 with kernel 2.6.9-1.667 #1 and gcc 3.4.2
     77#include <errno.h>   // necessary for Fedora core 2 with kernel 2.6.9-1.667 #1 and gcc 3.4.2
     78
    7679#include <stdlib.h>
    7780#include <fstream>
    … …  
    222225    fIsWobbleMode = env.GetValue("WobbleMode", kFALSE);
    223226    fComment      = env.GetValue("Comment",    "");
     227
     228    fNameSource = fNameSource.Strip(TString::kBoth);
     229    fCatalog    = fCatalog.Strip(TString::kBoth);
    224230}
    225231
    … …  
    246252    if (!IsValid())
    247253    {
    248         gLog << "Dataset: " << fName << " <invalid>" << endl;
     254        gLog << "Dataset: " << fName << " <invalid - no analysis number available>" << endl;
    249255        return;
    250256    }

  • trunk/MagicSoft/Mars/mranforest/MRanForest.cc

    r7398 r7409  
    419419
    420420        if(calcResolution)
    421             *fLog << "no. of tree    no. of nodes    resolution in % (from oob-data -> overest. of error)" << endl;
     421            *fLog << "no. of tree  no. of nodes  resolution in % (from oob-data -> overest. of error)" << endl;
    422422        else
    423             *fLog << "no. of tree    no. of nodes    rms in % (from oob-data -> overest. of error)" << endl;
     423            *fLog << "no. of tree  no. of nodes  rms in % (from oob-data -> overest. of error)" << endl;
    424424                     //        12345678901234567890123456789012345678901234567890
    425425    }
    … …  
    509509    // give running output
    510510    *fLog << inf << setw(5)  << fTreeNo;
    511     *fLog << inf << setw(20) << fRanTree->GetNumEndNodes();
    512     *fLog << inf << Form("%20.2f", ferr*100.);
     511    *fLog << inf << setw(18) << fRanTree->GetNumEndNodes();
     512    *fLog << inf << Form("%18.2f", ferr*100.);
    513513    *fLog << inf << endl;
    514514

  • trunk/MagicSoft/include-Classes/MMcFormat/MMcCorsikaRunHeader.h

    r6912 r7409  
    9797    Float_t GetSlopeSpec() const { return fSlopeSpec; }
    9898    Float_t GetWobbleMode() const { return fWobbleMode; }
     99    Float_t GetCorsikaVersion() const { return fCorsikaVersion; }
    99100
    100101    Int_t GetNumCT() const { return fNumCT; }
Note: See TracChangeset for help on using the changeset viewer.