Changeset 1840 for trunk/MagicSoft/Mars/macros/estfit.C

Timestamp:

03/19/03 16:29:30 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/macros/estfit.C (modified) (5 diffs)

trunk/MagicSoft/Mars/macros/estfit.C

@@ -1 @@
-#include <fstream.h>
+#include <TStopwatch.h>
+#include <TMinuit.h>
 
-#include "MTask.h"
 #include "MParList.h"
+#include "MTaskList.h"
+#include "MEvtLoop.h"
 
-#include "MDataChain.h"
-
-class MChisqEval : public MTask
-{
-private:
-    static const TString gsDefName;
-    static const TString gsDefTitle;
-
-    Double_t fChisq; //! Evaluated chi square
-
-    MData   *fData0; // Data Member one (monte carlo data or chisq function)
-    MData   *fData1; // Data Member two (measured data)
-
-    // --------------------------------------------------------------------------
-    //
-    // Implementation of SavePrimitive. Used to write the call to a constructor
-    // to a macro. In the original root implementation it is used to write
-    // gui elements to a macro-file.
-    //
-    void StreamPrimitive(ofstream &out) const
-    {
-        out << "   MChisqEval " << GetUniqueName() << ";";
-        if (fData0)
-            out << "   " << GetUniqueName() << ".SetY1(\"" << fData0->GetRule() << "\");" << endl;
-        if (fData1)
-            out << "   " << GetUniqueName() << ".SetY1(\"" << fData1->GetRule() << "\");" << endl;
-    }
-
-    enum { kIsOwner = BIT(14) };
-
-public:
-    MChisqEval(const char *name=NULL, const char *title=NULL) : fData0(NULL), fData1(NULL)// : fMatrix(mat), fColumn(col), fEvalE(evale)
-    {
-        fName  = name  ? name  : gsDefName.Data();
-        fTitle = title ? title : gsDefTitle.Data();
-    }
-
-    MChisqEval(MData *y1, const char *name=NULL, const char *title=NULL) : fData0(NULL), fData1(NULL)// : fMatrix(mat), fColumn(col), fEvalE(evale)
-    {
-        fName  = name  ? name  : gsDefName.Data();
-        fTitle = title ? title : gsDefTitle.Data();
-        SetY1(y1);
-    }
-    MChisqEval(MData *y1, MData *y2, const char *name=NULL, const char *title=NULL) : fData0(NULL), fData1(NULL)// : fMatrix(mat), fColumn(col), fEvalE(evale)
-    {
-        fName  = name  ? name  : gsDefName.Data();
-        fTitle = title ? title : gsDefTitle.Data();
-        SetY1(y1);
-        SetY2(y2);
-    }
-    ~MChisqEval()
-    {
-        if (fData0 && (fData0->TestBit(kCanDelete) || TestBit(kIsOwner)))
-            delete fData0;
-
-        if (fData1 && (fData1->TestBit(kCanDelete) || TestBit(kIsOwner)))
-            delete fData1;
-    }
-
-    void SetY1(MData *data)
-    {
-        // Set MC value
-        if (fData0 && (fData0->TestBit(kCanDelete) || TestBit(kIsOwner)))
-            delete fData0;
-        fData0 = data;
-        fData0->SetBit(kCanDelete);
-        AddToBranchList(fData0->GetDataMember());
-    }
-    void SetY2(MData *data)
-    {
-        // Set measured/estimated value
-        if (fData1 && (fData1->TestBit(kCanDelete) || TestBit(kIsOwner)))
-            delete fData1;
-        fData1 = data;
-        fData1->SetBit(kCanDelete);
-        AddToBranchList(fData1->GetDataMember());
-    }
-    void SetY1(const TString data) { SetY1(new MDataChain(data)); }
-    void SetY2(const TString data) { SetY2(new MDataChain(data)); }
-
-    void SetOwner(Bool_t o=kTRUE) { o ? SetBit(kIsOwner) : ResetBit(kIsOwner); }
-
-    Bool_t PreProcess(MParList *plist)
-    {
-        fChisq = 0;
-
-        if (!fData0)
-            return kFALSE;
-
-        if (!fData0->PreProcess(plist))
-            return kFALSE;
-
-        if (fData1)
-            if (!fData1->PreProcess(plist))
-                return kFALSE;
-
-        return kTRUE;
-    }
-
-    Bool_t Process()
-    {
-        const Double_t y1 = fData0->GetValue();
-        const Double_t y2 = fData1 ? fData1->GetValue() : 0;
-
-        const Double_t dy  = y2-y1;
-        const Double_t err = fData1 ? y1*y1 : 1;
-
-        fChisq += dy*dy/err;
-        return kTRUE;
-    }
-
-    Bool_t PostProcess()
-    {
-        fChisq /= GetNumExecutions();
-        return kTRUE;
-    }
-
-    Double_t GetChisq() const { return fChisq; }
- 
-    ClassDef(MChisqEval, 0)
-};
-
-ClassImp(MChisqEval);
-
-const TString MChisqEval::gsDefName  = "MChisqEval";
-const TString MChisqEval::gsDefTitle = "Evaluate a chisq";
+#include "MReadTree.h"
+#include "MHMatrix.h"
+#include "MChisqEval.h"
+#include "MMatrixLoop.h"
+#include "MDataMember.h"
+#include "MDataElement.h"
+#include "MEnergyEstParam.h"
 
 // --------------------------------------------------------------------------
-
-#include <TMinuit.h>
-
-#include "MEvtLoop.h"
-#include "MTaskList.h"
-#include "MEnergyEstParam.h"
 
 void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
@@ -156 +33 @@
 
 // --------------------------------------------------------------------------
-#include <TStopwatch.h>
-
-#include "MHMatrix.h"
-#include "MEnergyEst.h"
-
-#include "MDataMember.h"
-#include "MDataElement.h"
-
-#include "MReadTree.h"
-#include "MMatrixLoop.h"
 
 void estfit(Bool_t evalenergy=kFALSE)
@@ -177 +44 @@
     Int_t col = matrix.AddColumn(evalenergy?"MMcEvt.fEnergy":"MMcEvt.fImpact");
 
-    MEnergyEstParam eest;
-    eest.Add("MHillasSrc");
+    MEnergyEstParam eest("Hillas");
+    eest.Add("HillasSrc");
     eest.InitMapping(&matrix);
 
-    MReadTree read("Events", "gammas.root");
+    MReadTree read("Events", "~/ct1/MC_ON2.root");
     read.DisableAutoScheme();
 
@@ -211 +78 @@
     //
     TMinuit minuit(12);
+    minuit.SetPrintLevel(-1);
     minuit.SetFCN(fcn);
 
@@ -224 +92 @@
     //
     TArrayD fA(5);
-    fA[0] =  -2539; // [cm]
-    fA[1] =    900; // [cm]
-    fA[2] =   17.5; // [cm]
-    fA[3] =  1;//    4;
-    fA[4] =   58.3;
+    fA[0] = -3907.74; //4916.4;    //-2414.75;
+    fA[1] = 1162.3;   //149.549;   // 1134.28;
+    fA[2] = 199.351;  //-558.209;  // 132.932;
+    fA[3] = 0.403192; //0.270725;  //0.292845;
+    fA[4] = 121.921;  //107.001;   // 107.001;
 
     TArrayD fB(7);
-    fB[0] =    -8.64; // [GeV]
-    fB[1] =   -0.069; // [GeV]
-    fB[2] =  0.00066; // [GeV]
-    fB[3] =    0.033; // [GeV]
-    fB[4] = 0.000226; // [GeV]
-    fB[5] =  4.14e-8; // [GeV]
-    fB[6] =    -0.06;
+    fB[0] = -100;
+    fB[1] = 10;
+    fB[2] = 0.1;
+    fB[3] = 10;
+    fB[4] = -1;
+    fB[5] = -0.1;
+    fB[6] = -0.1;
 
     // Set starting values and step sizes for parameters