Changeset 2031


Ignore:
Timestamp:
04/28/03 16:51:29 (21 years ago)
Author:
moralejo
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r2030 r2031  
    11                                                 -*-*- END OF LINE -*-*-
     2 2003/04/28: Abelardo Moralejo
     3
     4   * mmontecarlo/MMcEnergyEst.[h,cc]
     5     - Lots of fixes after Thomas suggestions. Now cuts are not part
     6       of the class, but introduced via a new MFilter* member. Changed
     7       all Char_t* for TString. Changed own TMiniut pointer by gMinuit.
     8       Removed couts and used fLog instead. Function fcn is no longer
     9       declared external.
     10
    211 2003/04/27: Thomas Bretz
    312

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

    r1982 r2031  
    3030//                                                                         //
    3131// Class for otimizing the parameters of the energy estimator              //
     32//                                                                         //
     33// FIXME: the class must be made more flexible, allowing for different     //
     34// parametrizations to be used. Also a new class is needed, a container    //
     35// for the parameters of the energy estimator.                             //
     36// FIXME: the starting values of the parameters are now fixed.             //
    3237//                                                                         //
    3338//                                                                         //
    … …  
    5560#include "MLogManip.h"
    5661
    57 extern void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
    58 TMinuit *minuit;
     62
     63//------------------------------------------------------------------------
     64//
     65// fcn calculates the function to be minimized (using TMinuit::Migrad)
     66//
     67void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
     68{
     69  MEvtLoop *evtloop = (MEvtLoop*)gMinuit->GetObjectFit();
     70 
     71  MTaskList *tlist  = (MTaskList*)evtloop->GetParList()->FindObject("MTaskList"); // GetTaskList();
     72
     73  MChisqEval      *eval = (MChisqEval*)     tlist->FindObject("MChisqEval");
     74  MEnergyEstParam *eest = (MEnergyEstParam*)tlist->FindObject("MEnergyEstParam");
     75
     76  eest->SetCoeff(TArrayD(eest->GetNumCoeff(), par));
     77
     78  evtloop->Eventloop();
     79
     80  f = eval->GetChisq();
     81}
    5982
    6083ClassImp(MMcEnergyEst);
    … …  
    7194    SetHillasName("MHillas");
    7295    SetHillasSrcName("MHillasSrc");
    73     SetHadronnessName("MHadronness");
    7496}
    7597
    … …  
    84106  MFEventSelector selector;
    85107  selector.SetNumSelectEvts(fNevents);
    86   cout << "Read events from file '" << fInFile << "'" << endl;   
     108  *fLog << inf << "Read events from file '" << fInFile << "'" << endl;   
    87109
    88110  MReadTree read("Events");
    … …  
    91113  read.SetSelector(&selector);
    92114
    93   MFCT1SelFinal filterhadrons;
    94   filterhadrons.SetHadronnessName(fHadronnessName);
    95   filterhadrons.SetCuts(fMaxHadronness, fMaxAlpha, fMaxDist);
    96   filterhadrons.SetInverted();
    97 
    98   cout << "Define columns of matrix" << endl;
     115  *fLog << inf << "Define columns of matrix" << endl;
    99116  MHMatrix matrix;
    100117  Int_t colenergy = matrix.AddColumn("MMcEvt.fEnergy");
    … …  
    103120  if (colenergy < 0  ||  colimpact < 0)
    104121  {
    105     cout << "colenergy, colimpact = " << colenergy << ",  "
    106          << colimpact << endl;
     122    *fLog << err << dbginf << "colenergy, colimpact = " << colenergy << ",  "
     123          << colimpact << endl;
    107124    return;
    108125  }
    … …  
    112129  eest.InitMapping(&matrix);
    113130 
    114   cout << "--------------------------------------" << endl;
    115   cout << "Fill events into the matrix" << endl;
    116   if ( !matrix.Fill(&parlist, &read, &filterhadrons) )
     131  *fLog << inf << "--------------------------------------" << endl;
     132  *fLog << inf << "Fill events into the matrix" << endl;
     133  if ( !matrix.Fill(&parlist, &read, fEventFilter) )
    117134    return;
    118   cout << "Matrix was filled with " << matrix.GetNumRows()
    119       << " events" << endl; 
     135  *fLog << inf << "Matrix was filled with " << matrix.GetNumRows()
     136        << inf << " events" << endl; 
    120137
    121138  //-----------------------------------------------------------------------
    … …  
    123140  // Setup the eventloop which will be executed in the fcn of MINUIT
    124141  //
    125   cout << "--------------------------------------" << endl;
    126   cout << "Setup event loop to be used in 'fcn'" << endl;
     142  *fLog << inf << "--------------------------------------" << endl;
     143  *fLog << inf << "Setup event loop to be used in 'fcn'" << endl;
    127144
    128145  MTaskList tasklist;
    … …  
    148165
    149166
    150   cout << "Event loop was setup" << endl;
     167  *fLog << inf << "Event loop was setup" << endl;
    151168  MEvtLoop evtloop;
    152169  evtloop.SetParList(&parlist);
    … …  
    159176  // Be careful: This is not thread safe
    160177  //
    161   cout << "--------------------------------------" << endl;
    162   cout << "Start minimization" << endl;
    163 
    164   minuit = new TMinuit(12);
    165   minuit->SetPrintLevel(-1);
    166 
    167   minuit->SetFCN(fcn);
    168   minuit->SetObjectFit(&evtloop);
    169 
    170   // Ready for: minuit->mnexcm("SET ERR", arglist, 1, ierflg)
    171 
    172   if (minuit->SetErrorDef(1))
    173     {
    174       cout << "SetErrorDef failed." << endl;
     178  *fLog << inf << "--------------------------------------" << endl;
     179  *fLog << inf << "Start minimization" << endl;
     180
     181  gMinuit = new TMinuit(12);
     182  gMinuit->SetPrintLevel(-1);
     183
     184  gMinuit->SetFCN(fcn);
     185  gMinuit->SetObjectFit(&evtloop);
     186
     187  // Ready for: gMinuit->mnexcm("SET ERR", arglist, 1, ierflg)
     188
     189  if (gMinuit->SetErrorDef(1))
     190    {
     191      *fLog << err << dbginf << "SetErrorDef failed." << endl;
    175192      return;
    176193    }
    … …  
    211228      Double_t limup = 0;
    212229
    213       Bool_t rc = minuit->DefineParameter(i, name, vinit, step, limlo, limup);
     230      Bool_t rc = gMinuit->DefineParameter(i, name, vinit, step, limlo, limup);
    214231      if (!rc)
    215232        continue;
    216233
    217       cout << "Error in defining parameter #" << i << endl;
     234      *fLog << err << dbginf << "Error in defining parameter #" << i << endl;
    218235      return;
    219236    }
    … …  
    230247      Double_t limup = 0;
    231248
    232       Bool_t rc = minuit->DefineParameter(i+fA.GetSize(), name, vinit, step, limlo, limup);
     249      Bool_t rc = gMinuit->DefineParameter(i+fA.GetSize(), name, vinit, step, limlo, limup);
    233250      if (!rc)
    234251        continue;
    235252
    236       cout << "Error in defining parameter #" << i+fA.GetSize() << endl;
     253      *fLog << err << dbginf << "Error in defining parameter #" << i+fA.GetSize() << endl;
    237254      return;
    238255    }
    … …  
    244261
    245262  gLog.SetNullOutput(kTRUE);
    246   Bool_t rc = minuit->Migrad();
     263  Bool_t rc = gMinuit->Migrad();
    247264  gLog.SetNullOutput(kFALSE);
    248265 
    249266  if (rc)
    250267    {
    251       cout << "Migrad failed." << endl;
     268      *fLog << err << dbginf << "Migrad failed." << endl;
    252269      return;
    253270    }
    254271
    255   cout << endl;
     272  *fLog << inf << endl;
    256273  clock.Stop();
    257274  clock.Print();
    258   cout << endl;
    259 
    260   cout << "Resulting Chisq: " << minuit->fAmin << endl;
     275  *fLog << inf << endl;
     276
     277  *fLog << inf << "Resulting Chisq: " << gMinuit->fAmin << endl;
    261278
    262279  //
    … …  
    266283    {
    267284      Double_t x1, x2;
    268       minuit->GetParameter(i,x1,x2);
     285      gMinuit->GetParameter(i,x1,x2);
    269286      fA[i] = x1;
    270287    }
    … …  
    272289    {
    273290      Double_t x1, x2;
    274       minuit->GetParameter(i,x1,x2);
     291      gMinuit->GetParameter(i,x1,x2);
    275292      fB[i-fA.GetSize()] = x1;
    276293    }
    … …  
    278295  //    eest.Print();
    279296  eest.StopMapping();
    280   cout << "Minimization finished" << endl;
     297  *fLog << inf << "Minimization finished" << endl;
    281298
    282299}
    … …  
    289306{
    290307  for (Int_t i = 0; i < fA.GetSize(); i++)
    291     cout << "Parameter " << i << ": " << fA[i] << endl;
     308    *fLog << inf << "Parameter " << i << ": " << fA[i] << endl;
    292309
    293310  for (Int_t i = fA.GetSize(); i < fA.GetSize()+fB.GetSize(); i++)
    294     cout << "Parameter " << i << ": " << fB[i-fA.GetSize()] << endl;
     311    *fLog << inf << "Parameter " << i << ": " << fB[i-fA.GetSize()] << endl;
    295312
    296313  /*
    … …  
    298315    Double_t amin, edm, errdef;
    299316    Int_t nvpar, nparx, icstat;
    300     minuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
    301     minuit->mnprin(3, amin);
     317    gMinuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
     318    gMinuit->mnprin(3, amin);
    302319  */
    303320
    304321}
    305322
    306 //------------------------------------------------------------------------
    307 //
    308 // fcn calculates the function to be minimized (using TMinuit::Migrad)
    309 //
    310 void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
    311 {
    312   MEvtLoop *evtloop = (MEvtLoop*)minuit->GetObjectFit();
    313  
    314   MTaskList *tlist  = (MTaskList*)evtloop->GetParList()->FindObject("MTaskList"); // GetTaskList();
    315 
    316   MChisqEval      *eval = (MChisqEval*)     tlist->FindObject("MChisqEval");
    317   MEnergyEstParam *eest = (MEnergyEstParam*)tlist->FindObject("MEnergyEstParam");
    318 
    319   eest->SetCoeff(TArrayD(eest->GetNumCoeff(), par));
    320 
    321   evtloop->Eventloop();
    322 
    323   f = eval->GetChisq();
    324 }
    325 
    326 
    327 
     323

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

    r1964 r2031  
    1010#endif
    1111
     12#include "MFilter.h"
     13
    1214class MMcEnergyEst : public MParContainer
    1315{
    1416private:
    1517
    16   Char_t *fInFile, *fOutFile;
    17   Char_t *fHillasName, *fHillasSrcName, *fHadronnessName;
     18  TString fInFile, fOutFile;
     19  TString fHillasName;
     20  TString fHillasSrcName;
    1821  Int_t   fNevents;
    19   Float_t fMaxHadronness, fMaxAlpha, fMaxDist;
     22
     23  MFilter *fEventFilter; //!
    2024
    2125  TArrayD fA;
    … …  
    2529  MMcEnergyEst(const char *name=NULL, const char *title=NULL);
    2630
    27   void SetInFile(Char_t *name)         {fInFile = name;}
    28   void SetOutFile(Char_t *name)        {fOutFile = name;}
    29   void SetHillasName(Char_t *name)     {fHillasName = name;}
    30   void SetHillasSrcName(Char_t *name)  {fHillasSrcName = name;}
    31   void SetHadronnessName(Char_t *name) {fHadronnessName = name;}
    32   void SetNevents(Int_t n)             {fNevents = n;}
    33   void SetMaxHadronness(Float_t x)     {fMaxHadronness = x;}
    34   void SetMaxAlpha(Float_t x)          {fMaxAlpha = x;}
    35   void SetMaxDist(Float_t x)           {fMaxDist = x;}
     31  void SetInFile(const TString &name)         {fInFile = name;}
     32  void SetOutFile(const TString &name)        {fOutFile = name;}
     33  void SetHillasName(const TString &name)     {fHillasName = name;}
     34  void SetHillasSrcName(const TString &name)  {fHillasSrcName = name;}
     35  void SetEventFilter(MFilter *filter)        {fEventFilter = filter;}
     36  void SetNevents(Int_t n)                    {fNevents = n;}
    3637
    37   Char_t* GetInFile()         {return fInFile;}
    38   Char_t* GetOutFile()        {return fOutFile;}
    39   Char_t* GetHillasName()     {return fHillasName;}
    40   Char_t* GetHillasSrcName()  {return fHillasSrcName;}
    41   Char_t* GetHadronnessName() {return fHadronnessName;}
    42   Int_t   GetNevents()        {return fNevents;}
    43   Float_t GetMaxHadronness()  {return fMaxHadronness;}
    44   Float_t GetMaxAlpha()       {return fMaxAlpha;}
    45   Float_t GetMaxDist()        {return fMaxDist;}
     38  TString GetInFile()         const {return fInFile;}
     39  TString GetOutFile()        const {return fOutFile;}
     40  TString GetHillasName()     const {return fHillasName;}
     41  TString GetHillasSrcName()  const {return fHillasSrcName;}
     42  Int_t   GetNevents()        const {return fNevents;}
    4643
    4744  Double_t GetCoeff(Int_t i) { return i<fA.GetSize()? fA[i] : fB[i-fA.GetSize()]; }
Note: See TracChangeset for help on using the changeset viewer.