Changeset 6924

Timestamp:

04/11/05 16:15:08 (20 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• manalysis/MMatrixLoop.cc (modified) (3 diffs)
• manalysis/MMatrixLoop.h (modified) (3 diffs)
• mhflux/MHEnergyEst.cc (modified) (12 diffs)
• mjobs/MJOptimize.cc (modified) (4 diffs)
• mjobs/MJOptimize.h (modified) (2 diffs)

trunk/MagicSoft/Mars/manalysis/MMatrixLoop.cc

@@ -46 +46 @@
 //  through.
 //
-MMatrixLoop::MMatrixLoop(MHMatrix *mat, const char *name, const char *title) : fMatrix(mat)
+MMatrixLoop::MMatrixLoop(MHMatrix *mat, const char *name, const char *title) : fMatrix(mat), fOperationMode(kDefault)
 {
     fName  = name  ? name  : gsDefName.Data();
@@ -75 +75 @@
 Int_t MMatrixLoop::PreProcess(MParList *plist)
 {
-    fNumRow = 0;
+    fNumRow = fOperationMode==kOdd ? 1 : 0;
 
     return fMatrix ? kTRUE : kFALSE;
@@ -87 +87 @@
 Int_t MMatrixLoop::Process()
 {
-    return fMatrix->SetNumRow(fNumRow++);
+    const Int_t rc = fMatrix->SetNumRow(fNumRow);
+    fNumRow += fOperationMode==kDefault ? 1 : 2;
+    return rc;
 }

trunk/MagicSoft/Mars/manalysis/MMatrixLoop.h

@@ -10 +10 @@
 class MMatrixLoop : public MRead
 {
+public:
+    enum OperationMode_t {
+        kDefault,
+        kEven,
+        kOdd
+    };
 private:
     // MMatrixLoop
@@ -17 +23 @@
     MHMatrix *fMatrix;
     Int_t     fNumRow;    //! Number of dimensions of histogram
+
+    Byte_t fOperationMode;
 
     // MRead
@@ -33 +41 @@
     MMatrixLoop(MHMatrix *mat, const char *name=NULL, const char *title=NULL);
 
+    void SetOperationMode(OperationMode_t mode) { fOperationMode = mode; }
+
     ClassDef(MMatrixLoop, 0) // Task 'reading' events from a MHMatrix
 };

trunk/MagicSoft/Mars/mhflux/MHEnergyEst.cc

@@ -78 +78 @@
     fHResolution.SetDirectory(NULL);
     fHResolution.SetName("EnergyRes");
-    fHResolution.SetTitle("Histogram in \\frac{\\Delta E}{E_{mc}} vs E_{est} and E_{mc}");
+    fHResolution.SetTitle("Histogram in \\frac{\\Delta lg(E)}{lg(E_{mc})} vs E_{est} and E_{mc}");
     fHResolution.SetXTitle("E_{est} [GeV]");
     fHResolution.SetYTitle("E_{mc} [GeV]");
-    fHResolution.SetZTitle("\\frac{E_{est} - E_{mc}}{E_{mc}}");
+    fHResolution.SetZTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}");
 
     fHImpact.SetDirectory(NULL);
     fHImpact.SetName("Impact");
-    fHImpact.SetTitle("\\frac{\\Delta E}{E} vs Impact parameter");
+    fHImpact.SetTitle("\\frac{\\Delta lg(E)}{lg(E)} vs Impact parameter");
     fHImpact.SetXTitle("Impact parameter [m]");
-    fHImpact.SetYTitle("\\frac{E_{est} - E_{mc}}{E_{mc}}");
+    fHImpact.SetYTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}");
 
     fHEnergy.Sumw2();
@@ -97 +97 @@
     binst.SetEdgesCos(50, 0, 60);
     binsi.SetEdges(10, 0, 400);
-    binsr.SetEdges(10, 0, 1);
+    binsr.SetEdges(50, -1.3, 1.3);
 
     SetBinning(&fHEnergy,     &binse, &binse, &binst);
@@ -164 +164 @@
     const Double_t imp   = fMatrix ? GetVal(1) : fMcEvt->GetImpact()/100;
     const Double_t theta = fMatrix ? GetVal(2) : fMcEvt->GetTelescopeTheta()*TMath::RadToDeg();
-    const Double_t res   = (eest-etru)/etru;
+    const Double_t res   = (log10(eest)-log10(etru));///log10(etru);
+    const Double_t resE   = (eest-etru)/etru;
 
     fHEnergy.Fill(eest, etru, theta, w);
-    fHResolution.Fill(eest, etru, TMath::Abs(res), w);
-    fHImpact.Fill(imp, TMath::Abs(res), w);
-
-    fChisq += TMath::Abs(res);//*res;
+    fHResolution.Fill(eest, etru, resE, w);
+    fHImpact.Fill(imp, resE, w);
+
+    fChisq += res*res;
     fBias  += res;
 
@@ -181 +182 @@
     fBias  /= GetNumExecutions();
 
-    Double_t res = fChisq; //TMath::Sqrt(fChisq - fBias*fBias);
-
-    fResult->SetVal(TMath::IsNaN(res)?0:res);/// GetNumExecutions());
-
-    *fLog << all << "Mean Energy Resoltuion: " << Form("%.1f%%", fResult->GetVal()*100) << endl;
-    *fLog << all << "Energy Bias at:         " << Form("%.1f%%", fBias*100) << endl;
+    Double_t sigma = TMath::Sqrt(fChisq);//+TMath::Abs(fBias);
+    fResult->SetVal(sigma);
+
+    *fLog << all << "Mean log10(Energy) Resoltion: " << Form("%.1f%%", TMath::Sqrt(fChisq-fBias*fBias)*100) << endl;
+    *fLog << all << "Mean log10(Energy) Bias:      " << Form("%.1f%%", fBias*100) << endl;
 
     return kTRUE;
@@ -239 +239 @@
             }
 
-            //pad->GetPad(1)->GetPad(2)->cd(2);
+            pad->GetPad(1)->GetPad(2)->cd(2);
+            /*=*/fHResolution.ProjectionZ("Resolution");
             ///*h =*/ fHImpact.ProjectionX("Impact", -1, 9999, "e");
         }
@@ -294 +295 @@
     gPad->SetBorderMode(0);
     gPad->SetLogx();
+    gPad->SetGridx();
+    gPad->SetGridy();
 
     //gROOT->GetListOfCleanups()->Add(gPad); // WHY?
@@ -301 +304 @@
     h2->SetBit(kCanDelete);
     h2->SetFillColor(kBlue);
+    h2->SetLineColor(kRed);
 
     TH1D *h1 = h2->ProfileX(Form("Prof%s", h.GetName()), -1, 9999, "s");
@@ -310 +314 @@
     h1->SetStats(kFALSE);
 
- 
-    //h1->Draw("E3");
-    h2->Draw();
+    h2->Draw("");
+    h1->Draw("E3same");
     h1->Draw("Chistsame");
 
@@ -346 +349 @@
     gPad->SetBorderMode(0);
 
+    gPad->SetGridx();
+    gPad->SetGridy();
     gPad->SetLogx();
     h = (TH1D*)fHEnergy.Project3D("ey");
@@ -376 +381 @@
     pad3->Divide(2, 1, 1e-10, 1e-10);
     pad3->cd(1);
-    gPad->SetBorderMode(0);/*
+    gPad->SetBorderMode(0);
+    gPad->SetGridx();
+    gPad->SetGridy();
+    h = fHEnergy.Project3D("ez");
+    h->SetTitle("Zenith Angle Distribution");
+    h->SetBit(TH1::kNoStats);
+    h->SetDirectory(NULL);
+    h->SetXTitle("\\Theta [\\circ]");
+    h->SetBit(kCanDelete);
+    h->Draw();
+
+    pad3->cd(2);
+    gPad->SetBorderMode(0);
+    gPad->SetGridx();
+    gPad->SetGridy();
+    /*
     h = fHImpact.ProjectionX("Impact", -1, 9999, "e");
     h->SetBit(TH1::kNoStats);
@@ -384 +404 @@
     h->SetBit(kCanDelete);
     h->Draw();*/
-    h = fHEnergy.Project3D("ez");
-    h->SetTitle("Distribution of Theta");
-    h->SetBit(TH1::kNoStats);
+    // ----------------------------------------
+    h = fHResolution.ProjectionZ("Resolution");
+    //h->SetXTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}");
+    h->SetYTitle("Counts");
+    h->SetTitle("Distribution of \\Delta lg(E) / lg(E_{mc})");
     h->SetDirectory(NULL);
-    h->SetXTitle("\\Theta [\\circ]");
     h->SetBit(kCanDelete);
-    h->Draw();
-
-    pad3->cd(2);
-    gPad->SetBorderMode(0);
+    h->GetXaxis()->SetRangeUser(-1.3, 1.3);
+    h->Draw("");
+    // ----------------------------------------
 
     pad->cd(2);
@@ -413 +433 @@
     h = MakePlot(fHResolution, "zy");
     h->SetXTitle("E_{mc} [GeV]");
-    h->SetYTitle("\\frac{E_{est} - E_{mc}}{E_{mc}}");
-    h->SetTitle("Energy resolution \\Delta E / E vs Monte Carlo energy E_{mc}");
-    h->SetMinimum(0);
-    h->SetMaximum(1);
+    h->SetYTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}");
+    h->SetTitle("Energy resolution \\Delta lg(E) / lg(E) vs Monte Carlo energy E_{mc}");
+    h->SetMinimum(-1.3);
+    h->SetMaximum(1.3);
 
     pad2->cd(3);
     h = MakePlot(fHResolution, "zx");
     h->SetXTitle("E_{est} [GeV]");
-    h->SetYTitle("\\frac{E_{est} - E_{mc}}{E_{mc}}");
-    h->SetTitle("Energy Resolution \\Delta E / E vs estimated Energy E_{est}");
-    h->SetMinimum(0);
-    h->SetMaximum(1);
+    h->SetYTitle("\\frac{lg(E_{est}) - lg(E_{mc})}{lg(E_{mc})}");
+    h->SetTitle("Energy Resolution \\Delta lg(E) / lg(E) vs estimated Energy E_{est}");
+    h->SetMinimum(-1.3);
+    h->SetMaximum(1.3);
 }
 

trunk/MagicSoft/Mars/mjobs/MJOptimize.cc

@@ -228 +228 @@
 }
 
-MJOptimize::MJOptimize() : fDebug(-1), fNumEvents(0), fType(kSimplex), fNumMaxCalls(0), fTolerance(0)
+MJOptimize::MJOptimize() : fDebug(-1), fNumEvents(0), fType(kSimplex), fNumMaxCalls(0), fTolerance(0), fTestTrain(kFALSE)
 {
     fRules.SetOwner();
@@ -642 +642 @@
     }
 
+    MMatrixLoop *loop = dynamic_cast<MMatrixLoop*>(parlist.FindTask("MRead"));
+
     TString txt("Starting ");
     switch (fType)
@@ -666 +668 @@
     *fLog << inf << "Number of Parameters: " << fParameters.GetSize() << endl;
 
+    // In case the reader is the matrix loop and testrain is enabled
+    // switch on even mode...
+    if (loop && fTestTrain)
+        loop->SetOperationMode(MMatrixLoop::kEven);
+
     if (!Optimize(evtloop))
         return kFALSE;
@@ -672 +679 @@
 
     fEvtLoop->SetDisplay(fDisplay);
-    return Fcn(fParameters);
+    if (!Fcn(fParameters))
+        return kFALSE;
+
+    // In case the reader is the matrix loop and testrain is enabled
+    // switch on odd mode...
+    if (!loop || !fTestTrain)
+        return kTRUE;
+
+    loop->SetOperationMode(MMatrixLoop::kOdd);
+
     // Done already in Fcn
     // list.SetVariables(fParameters);
+    return Fcn(fParameters);
 }
 

trunk/MagicSoft/Mars/mjobs/MJOptimize.h

@@ -84 +84 @@
     UInt_t  fNumMaxCalls;
     Float_t fTolerance;
+    Bool_t  fTestTrain;
 
     Bool_t Optimize(MEvtLoop &evtloop);
@@ -108 +109 @@
     void SetNumMaxCalls(UInt_t num=0) { fNumMaxCalls=num; }
     void SetTolerance(Float_t tol=0)  { fTolerance=tol; }
+    void EnableTestTrain(Bool_t b=kTRUE) { fTestTrain=b; }
 
     // Parameter access