Changeset 2031 for trunk/MagicSoft/Mars/mmontecarlo

Timestamp:

04/28/03 16:51:29 (22 years ago)

Author:

moralejo

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mmontecarlo

Files:

• MMcEnergyEst.cc (modified) (18 diffs)
• MMcEnergyEst.h (modified) (2 diffs)

trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.cc

@@ -30 +30 @@
 //                                                                         //
 // Class for otimizing the parameters of the energy estimator              //
+//                                                                         //
+// FIXME: the class must be made more flexible, allowing for different     //
+// parametrizations to be used. Also a new class is needed, a container    //
+// for the parameters of the energy estimator.                             //
+// FIXME: the starting values of the parameters are now fixed.             //
 //                                                                         //
 //                                                                         //
@@ -55 +60 @@
 #include "MLogManip.h"
 
-extern void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
-TMinuit *minuit;
+
+//------------------------------------------------------------------------
+//
+// fcn calculates the function to be minimized (using TMinuit::Migrad)
+//
+void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
+{
+  MEvtLoop *evtloop = (MEvtLoop*)gMinuit->GetObjectFit();
+ 
+  MTaskList *tlist  = (MTaskList*)evtloop->GetParList()->FindObject("MTaskList"); // GetTaskList();
+
+  MChisqEval      *eval = (MChisqEval*)     tlist->FindObject("MChisqEval");
+  MEnergyEstParam *eest = (MEnergyEstParam*)tlist->FindObject("MEnergyEstParam");
+
+  eest->SetCoeff(TArrayD(eest->GetNumCoeff(), par));
+
+  evtloop->Eventloop();
+
+  f = eval->GetChisq();
+}
 
 ClassImp(MMcEnergyEst);
@@ -71 +94 @@
     SetHillasName("MHillas");
     SetHillasSrcName("MHillasSrc");
-    SetHadronnessName("MHadronness");
 }
 
@@ -84 +106 @@
   MFEventSelector selector;
   selector.SetNumSelectEvts(fNevents);
-  cout << "Read events from file '" << fInFile << "'" << endl;    
+  *fLog << inf << "Read events from file '" << fInFile << "'" << endl;    
 
   MReadTree read("Events");
@@ -91 +113 @@
   read.SetSelector(&selector);
 
-  MFCT1SelFinal filterhadrons;
-  filterhadrons.SetHadronnessName(fHadronnessName);
-  filterhadrons.SetCuts(fMaxHadronness, fMaxAlpha, fMaxDist);
-  filterhadrons.SetInverted();
-
-  cout << "Define columns of matrix" << endl;
+  *fLog << inf << "Define columns of matrix" << endl;
   MHMatrix matrix;
   Int_t colenergy = matrix.AddColumn("MMcEvt.fEnergy");
@@ -103 +120 @@
   if (colenergy < 0  ||  colimpact < 0)
   {
-    cout << "colenergy, colimpact = " << colenergy << ",  " 
-         << colimpact << endl;
+    *fLog << err << dbginf << "colenergy, colimpact = " << colenergy << ",  " 
+          << colimpact << endl;
     return;
   }
@@ -112 +129 @@
   eest.InitMapping(&matrix);
   
-  cout << "--------------------------------------" << endl;
-  cout << "Fill events into the matrix" << endl;
-  if ( !matrix.Fill(&parlist, &read, &filterhadrons) )
+  *fLog << inf << "--------------------------------------" << endl;
+  *fLog << inf << "Fill events into the matrix" << endl;
+  if ( !matrix.Fill(&parlist, &read, fEventFilter) )
     return;
-  cout << "Matrix was filled with " << matrix.GetNumRows() 
-       << " events" << endl;  
+  *fLog << inf << "Matrix was filled with " << matrix.GetNumRows() 
+        << inf << " events" << endl;  
 
   //-----------------------------------------------------------------------
@@ -123 +140 @@
   // Setup the eventloop which will be executed in the fcn of MINUIT 
   //
-  cout << "--------------------------------------" << endl;
-  cout << "Setup event loop to be used in 'fcn'" << endl;
+  *fLog << inf << "--------------------------------------" << endl;
+  *fLog << inf << "Setup event loop to be used in 'fcn'" << endl;
 
   MTaskList tasklist;
@@ -148 +165 @@
 
 
-  cout << "Event loop was setup" << endl;
+  *fLog << inf << "Event loop was setup" << endl;
   MEvtLoop evtloop;
   evtloop.SetParList(&parlist);
@@ -159 +176 @@
   // Be careful: This is not thread safe
   //
-  cout << "--------------------------------------" << endl;
-  cout << "Start minimization" << endl;
-
-  minuit = new TMinuit(12);
-  minuit->SetPrintLevel(-1);
-
-  minuit->SetFCN(fcn);
-  minuit->SetObjectFit(&evtloop);
-
-  // Ready for: minuit->mnexcm("SET ERR", arglist, 1, ierflg)
-
-  if (minuit->SetErrorDef(1))
-    {
-      cout << "SetErrorDef failed." << endl;
+  *fLog << inf << "--------------------------------------" << endl;
+  *fLog << inf << "Start minimization" << endl;
+
+  gMinuit = new TMinuit(12);
+  gMinuit->SetPrintLevel(-1);
+
+  gMinuit->SetFCN(fcn);
+  gMinuit->SetObjectFit(&evtloop);
+
+  // Ready for: gMinuit->mnexcm("SET ERR", arglist, 1, ierflg)
+
+  if (gMinuit->SetErrorDef(1))
+    {
+      *fLog << err << dbginf << "SetErrorDef failed." << endl;
       return;
     }
@@ -211 +228 @@
       Double_t limup = 0; 
 
-      Bool_t rc = minuit->DefineParameter(i, name, vinit, step, limlo, limup);
+      Bool_t rc = gMinuit->DefineParameter(i, name, vinit, step, limlo, limup);
       if (!rc)
         continue;
 
-      cout << "Error in defining parameter #" << i << endl;
+      *fLog << err << dbginf << "Error in defining parameter #" << i << endl;
       return;
     }
@@ -230 +247 @@
       Double_t limup = 0;
 
-      Bool_t rc = minuit->DefineParameter(i+fA.GetSize(), name, vinit, step, limlo, limup);
+      Bool_t rc = gMinuit->DefineParameter(i+fA.GetSize(), name, vinit, step, limlo, limup);
       if (!rc)
         continue;
 
-      cout << "Error in defining parameter #" << i+fA.GetSize() << endl;
+      *fLog << err << dbginf << "Error in defining parameter #" << i+fA.GetSize() << endl;
       return;
     }
@@ -244 +261 @@
 
   gLog.SetNullOutput(kTRUE);
-  Bool_t rc = minuit->Migrad();
+  Bool_t rc = gMinuit->Migrad();
   gLog.SetNullOutput(kFALSE);
   
   if (rc)
     {
-      cout << "Migrad failed." << endl;
+      *fLog << err << dbginf << "Migrad failed." << endl;
       return;
     }
 
-  cout << endl;
+  *fLog << inf << endl;
   clock.Stop();
   clock.Print();
-  cout << endl;
-
-  cout << "Resulting Chisq: " << minuit->fAmin << endl;
+  *fLog << inf << endl;
+
+  *fLog << inf << "Resulting Chisq: " << gMinuit->fAmin << endl;
 
   //
@@ -266 +283 @@
     {
       Double_t x1, x2;
-      minuit->GetParameter(i,x1,x2);
+      gMinuit->GetParameter(i,x1,x2);
       fA[i] = x1;
     }
@@ -272 +289 @@
     {
       Double_t x1, x2;
-      minuit->GetParameter(i,x1,x2);
+      gMinuit->GetParameter(i,x1,x2);
       fB[i-fA.GetSize()] = x1;
     }
@@ -278 +295 @@
   //    eest.Print();
   eest.StopMapping();
-  cout << "Minimization finished" << endl;
+  *fLog << inf << "Minimization finished" << endl;
 
 }
@@ -289 +306 @@
 {
   for (Int_t i = 0; i < fA.GetSize(); i++)
-    cout << "Parameter " << i << ": " << fA[i] << endl; 
+    *fLog << inf << "Parameter " << i << ": " << fA[i] << endl; 
 
   for (Int_t i = fA.GetSize(); i < fA.GetSize()+fB.GetSize(); i++)
-    cout << "Parameter " << i << ": " << fB[i-fA.GetSize()] << endl; 
+    *fLog << inf << "Parameter " << i << ": " << fB[i-fA.GetSize()] << endl; 
 
   /*
@@ -298 +315 @@
     Double_t amin, edm, errdef;
     Int_t nvpar, nparx, icstat;
-    minuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
-    minuit->mnprin(3, amin);
+    gMinuit->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
+    gMinuit->mnprin(3, amin);
   */
 
 }
 
-//------------------------------------------------------------------------
-//
-// fcn calculates the function to be minimized (using TMinuit::Migrad)
-//
-void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
-{
-  MEvtLoop *evtloop = (MEvtLoop*)minuit->GetObjectFit();
- 
-  MTaskList *tlist  = (MTaskList*)evtloop->GetParList()->FindObject("MTaskList"); // GetTaskList();
-
-  MChisqEval      *eval = (MChisqEval*)     tlist->FindObject("MChisqEval");
-  MEnergyEstParam *eest = (MEnergyEstParam*)tlist->FindObject("MEnergyEstParam");
-
-  eest->SetCoeff(TArrayD(eest->GetNumCoeff(), par));
-
-  evtloop->Eventloop();
-
-  f = eval->GetChisq();
-}
-
-
-
+

trunk/MagicSoft/Mars/mmontecarlo/MMcEnergyEst.h

@@ -10 +10 @@
 #endif
 
+#include "MFilter.h"
+
 class MMcEnergyEst : public MParContainer
 {
 private:
 
-  Char_t *fInFile, *fOutFile;
-  Char_t *fHillasName, *fHillasSrcName, *fHadronnessName;
+  TString fInFile, fOutFile;
+  TString fHillasName;
+  TString fHillasSrcName;
   Int_t   fNevents;
-  Float_t fMaxHadronness, fMaxAlpha, fMaxDist;
+
+  MFilter *fEventFilter; //!
 
   TArrayD fA;
@@ -25 +29 @@
   MMcEnergyEst(const char *name=NULL, const char *title=NULL);
 
-  void SetInFile(Char_t *name)         {fInFile = name;}
-  void SetOutFile(Char_t *name)        {fOutFile = name;}
-  void SetHillasName(Char_t *name)     {fHillasName = name;}
-  void SetHillasSrcName(Char_t *name)  {fHillasSrcName = name;}
-  void SetHadronnessName(Char_t *name) {fHadronnessName = name;}
-  void SetNevents(Int_t n)             {fNevents = n;}
-  void SetMaxHadronness(Float_t x)     {fMaxHadronness = x;}
-  void SetMaxAlpha(Float_t x)          {fMaxAlpha = x;}
-  void SetMaxDist(Float_t x)           {fMaxDist = x;}
+  void SetInFile(const TString &name)         {fInFile = name;}
+  void SetOutFile(const TString &name)        {fOutFile = name;}
+  void SetHillasName(const TString &name)     {fHillasName = name;}
+  void SetHillasSrcName(const TString &name)  {fHillasSrcName = name;}
+  void SetEventFilter(MFilter *filter)        {fEventFilter = filter;}
+  void SetNevents(Int_t n)                    {fNevents = n;}
 
-  Char_t* GetInFile()         {return fInFile;}
-  Char_t* GetOutFile()        {return fOutFile;}
-  Char_t* GetHillasName()     {return fHillasName;}
-  Char_t* GetHillasSrcName()  {return fHillasSrcName;}
-  Char_t* GetHadronnessName() {return fHadronnessName;}
-  Int_t   GetNevents()        {return fNevents;}
-  Float_t GetMaxHadronness()  {return fMaxHadronness;}
-  Float_t GetMaxAlpha()       {return fMaxAlpha;}
-  Float_t GetMaxDist()        {return fMaxDist;}
+  TString GetInFile()         const {return fInFile;}
+  TString GetOutFile()        const {return fOutFile;}
+  TString GetHillasName()     const {return fHillasName;}
+  TString GetHillasSrcName()  const {return fHillasSrcName;}
+  Int_t   GetNevents()        const {return fNevents;}
 
   Double_t GetCoeff(Int_t i) { return i<fA.GetSize()? fA[i] : fB[i-fA.GetSize()]; }