Changes between Version 23 and Version 24 of DatabaseBasedAnalysis/Examples

Timestamp:

08/16/18 20:33:58 (6 years ago)

Author:

tbretz

Comment:

--

DatabaseBasedAnalysis/Examples

@@ -423 +423 @@
     hnew.DrawCopy();
     hold.DrawCopy("same");
+}
+}}}
+
+
+== Optim Disp (Leakage==0) ==
+
+{{{#!Spoiler
+{{{#!cpp
+#include <iostream>
+#include <iomanip>
+
+#include <TMath.h>
+#include <TH2.h>
+#include <TChain.h>
+#include <TCanvas.h>
+
+void optimdisp2()
+{
+    // Create chain for the tree Result
+    // This is just easier than using TFile/TTree
+    TChain c("Result");
+
+    // Add the input file to the
+    c.AddFile("simulation.root");
+
+    // Define variables for all leaves to be accessed
+    // By definition rootifysql writes only doubles
+    double FileId, EvtNumber, X, Y, MeanX, MeanY, Width, Length, CosDelta, SinDelta,
+        M3Long, SlopeLong, Leakage1, SlopeSpreadWeighted, Size,
+        ConcCore, ConcCOG, NumIslands, NumUsedPixels, Zd;
+
+    // Connect the variables to the cordesponding leaves
+    c.SetBranchAddress("FileId", &FileId);
+    c.SetBranchAddress("EvtNumber", &EvtNumber);
+    c.SetBranchAddress("X", &X);
+    c.SetBranchAddress("Y", &Y);
+    c.SetBranchAddress("MeanX", &MeanX);
+    c.SetBranchAddress("MeanY", &MeanY);
+    c.SetBranchAddress("Width", &Width);
+    c.SetBranchAddress("Length", &Length);
+    c.SetBranchAddress("CosDelta", &CosDelta);
+    c.SetBranchAddress("SinDelta", &SinDelta);
+    c.SetBranchAddress("M3Long", &M3Long);
+    c.SetBranchAddress("SlopeLong", &SlopeLong);
+    c.SetBranchAddress("Leakage1", &Leakage1);
+    c.SetBranchAddress("NumIslands", &NumIslands);
+    c.SetBranchAddress("NumUsedPixels", &NumUsedPixels);
+    c.SetBranchAddress("Size", &Size);
+
+    // Set some constants (they could be included in the database
+    // in the future)
+    double mm2deg = +0.0117193246260285378;
+
+    // -------------------- Source dependent parameter calculation -------------------
+
+    TH2F hist("", "",
+              26, 1.300-0.0025, 1.430-0.0025,
+              20, 0.040-0.0025, 0.140-0.0025);
+
+    for (float xi0=1.30; xi0<1.43; xi0+=0.005)
+        for (float slope0=0.04; slope0<0.14; slope0+=0.005)
+        {
+            int cnt_on = 0;
+
+            for (int i=0; i<c.GetEntries(); i++)
+            {
+                // read the i-th event from the file
+                c.GetEntry(i);
+
+                // First calculate all cuts to speedup the analysis
+                double area = TMath::Pi()*Width*Length;
+
+                // The abberation correction does increase also Width and Length by 1.02
+
+                bool cutq = NumIslands<3.5 && NumUsedPixels>5.5 && Leakage1==0;
+                if (!cutq)
+                    continue;
+
+                bool cut0 = area < log10(Size)*898-1535;
+                if (!cut0)
+                    continue;
+
+                int angle = 0;
+
+                // -------------------- Source dependent parameter calculation -------------------
+
+                double cr = cos(angle*TMath::DegToRad());
+                double sr = sin(angle*TMath::DegToRad());
+
+                double px = cr*X-sr*Y;
+                double py = cr*Y+sr*X;
+
+                double dx = MeanX - px*1.02;
+                double dy = MeanY - py*1.02;
+
+                double norm = sqrt(dx*dx + dy*dy);
+                double dist = norm*mm2deg;
+
+                double lx = min(max((CosDelta*dy - SinDelta*dx)/norm, -1.), 1.);
+                double ly = min(max((CosDelta*dx + SinDelta*dy)/norm, -1.), 1.);
+
+                double alpha = asin(lx);
+                double sgn   = TMath::Sign(1., ly);
+
+                // ------------------------------- Application ----------------------------------
+
+                double m3l   = M3Long*sgn*mm2deg;
+                double slope = SlopeLong*sgn/mm2deg;
+
+                // --------------------------------- Analysis -----------------------------------
+
+                double xi = xi0 + slope0 *slope;
+
+                double sign1 = m3l+0.07;
+                double sign2 = (dist-0.5)*7.2-slope;
+
+                double disp  = (sign1<0 || sign2<0 ? -xi : xi)*(1-Width/Length);
+
+                double thetasq = disp*disp + dist*dist - 2*disp*dist*sqrt(1-lx*lx);
+
+                if (thetasq<0.024)
+                    cnt_on++;
+            }
+
+            hist.Fill(xi0, slope0, cnt_on);
+
+            cout << xi0 << " " << slope0 << " " << cnt_on << endl;
+        }
+
+    TCanvas *canv = new TCanvas;
+    canv->SetTopMargin(0.01);
+
+    hist.SetStats(kFALSE);
+    hist.SetContour(100);
+    hist.SetXTitle("Xi_{0} / deg");
+    hist.SetYTitle("Slope_{0} #dot ns/deg^{2} ");
+    hist.GetYaxis()->SetTitleOffset(1.2);
+    hist.DrawCopy("colz cont2");
+
+}
+}}}
+}}}
+
+Results in
+
+[[Image(Result.png)]]
+
+== Optim Disp (Leakage>0) ==
+
+{{{#!Spoiler
+{{{#!cpp
+#include <iostream>
+#include <iomanip>
+
+#include <TMath.h>
+#include <TH2.h>
+#include <TChain.h>
+#include <TCanvas.h>
+
+void optimdisp3()
+{
+    // Create chain for the tree Result
+    // This is just easier than using TFile/TTree
+    TChain c("Result");
+
+    // Add the input file to the
+    c.AddFile("simulation.root");
+
+    // Define variables for all leaves to be accessed
+    // By definition rootifysql writes only doubles
+    double FileId, EvtNumber, X, Y, MeanX, MeanY, Width, Length, CosDelta, SinDelta,
+        M3Long, SlopeLong, Leakage1, SlopeSpreadWeighted, Size,
+        ConcCore, ConcCOG, NumIslands, NumUsedPixels, Zd;
+
+    // Connect the variables to the cordesponding leaves
+    c.SetBranchAddress("FileId", &FileId);
+    c.SetBranchAddress("EvtNumber", &EvtNumber);
+    c.SetBranchAddress("X", &X);
+    c.SetBranchAddress("Y", &Y);
+    c.SetBranchAddress("MeanX", &MeanX);
+    c.SetBranchAddress("MeanY", &MeanY);
+    c.SetBranchAddress("Width", &Width);
+    c.SetBranchAddress("Length", &Length);
+    c.SetBranchAddress("CosDelta", &CosDelta);
+    c.SetBranchAddress("SinDelta", &SinDelta);
+    c.SetBranchAddress("M3Long", &M3Long);
+    c.SetBranchAddress("SlopeLong", &SlopeLong);
+    c.SetBranchAddress("Leakage1", &Leakage1);
+    c.SetBranchAddress("NumIslands", &NumIslands);
+    c.SetBranchAddress("NumUsedPixels", &NumUsedPixels);
+    //c.SetBranchAddress("SlopeSpreadWeighted", &SlopeSpreadWeighted);
+    c.SetBranchAddress("Size", &Size);
+    c.SetBranchAddress("Zd", &Zd);
+    //c.SetBranchAddress("ConcCOG", &ConcCOG);
+    //c.SetBranchAddress("ConcCore", &ConcCore);
+
+    // Set some constants (they could be included in the database
+    // in the future)
+    double mm2deg = +0.0117193246260285378;
+    //double abberation = 1.02;
+
+    // -------------------- Source dependent parameter calculation -------------------
+
+    TH2F hist("", "",
+              32, 0.800-0.0250, 2.400-0.0250,
+              60, 3.000-0.0250, 6.000-0.0250);
+
+
+    for (float l0=0.8; l0<2.4; l0+=0.05)
+        for (float l1=3; l1<6; l1+=0.05)
+        {
+            int cnt_on = 0;
+
+            for (int i=0; i<c.GetEntries(); i++)
+            {
+                // read the i-th event from the file
+                c.GetEntry(i);
+
+                // First calculate all cuts to speedup the analysis
+                double area = TMath::Pi()*Width*Length;
+
+                // The abberation correction does increase also Width and Length by 1.02
+
+                bool cutq = NumIslands<3.5 && NumUsedPixels>5.5 && Leakage1>0;
+                if (!cutq)
+                    continue;
+
+                bool cut0 = area < log10(Size)*898-1535;
+                if (!cut0)
+                    continue;
+
+                int angle = 0;
+
+                // -------------------- Source dependent parameter calculation -------------------
+
+                double cr = cos(angle*TMath::DegToRad());
+                double sr = sin(angle*TMath::DegToRad());
+
+                double px = cr*X-sr*Y;
+                double py = cr*Y+sr*X;
+
+                double dx = MeanX - px*1.02;
+                double dy = MeanY - py*1.02;
+
+                double norm = sqrt(dx*dx + dy*dy);
+                double dist = norm*mm2deg;
+
+                double lx = min(max((CosDelta*dy - SinDelta*dx)/norm, -1.), 1.);
+                double ly = min(max((CosDelta*dx + SinDelta*dy)/norm, -1.), 1.);
+
+                double alpha = asin(lx);
+                double sgn   = TMath::Sign(1., ly);
+
+                // ------------------------------- Application ----------------------------------
+
+                double m3l   = M3Long*sgn*mm2deg;
+                double slope = SlopeLong*sgn/mm2deg;
+
+                // --------------------------------- Analysis -----------------------------------
+
+                double xi = 1.36 + 0.085*slope + l0 *(1-1/(1+l1*Leakage1));
+
+                // Optim: 1.36 + 0.085
+
+                double sign1 = m3l+0.07;
+                double sign2 = (dist-0.5)*7.2-slope;
+
+                double disp  = (sign1<0 || sign2<0 ? -xi : xi)*(1-Width/Length);
+
+                double thetasq = disp*disp + dist*dist - 2*disp*dist*sqrt(1-lx*lx);
+
+                if (thetasq<0.024)
+                    cnt_on++;
+            }
+
+            hist.Fill(l0, l1, cnt_on);
+
+            cout << l0 << " " << l1 << " " << cnt_on << endl;
+        }
+
+    TCanvas *canv = new TCanvas;
+    canv->SetTopMargin(0.01);
+
+    hist.SetStats(kFALSE);
+    hist.SetContour(100);
+    hist.SetXTitle("L_{0} / deg");
+    hist.SetYTitle("L_{1} / deg");
+    hist.GetYaxis()->SetTitleOffset(1.2);
+    hist.DrawCopy("colz cont2");
+
+}
+}}}
+}}}
+
+Results in
+
+[[Image(Result2.png)]]
+
+== Energy Optimization == 
+
+{{{#!cpp
+#include <iostream>
+#include <iomanip>
+
+#include <TMath.h>
+#include <TF1.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TProfile.h>
+#include <TChain.h>
+#include <TGraph.h>
+#include <TCanvas.h>
+#include <TStopwatch.h>
+
+void optimenergy()
+{
+    // Create chain for the tree Result
+    // This is just easier than using TFile/TTree
+    TChain c("Result");
+
+    // Add the input file to the
+    c.AddFile("simulation.root");
+
+    // Define variables for all leaves to be accessed
+    // By definition rootifysql writes only doubles
+    double X, Y, MeanX, MeanY, Width, Length, CosDelta, SinDelta,
+        M3Long, SlopeLong, Leakage1, SlopeSpreadWeighted, Size,
+        ConcCore, ConcCOG, NumIslands, NumUsedPixels, Zd, Energy;
+
+    // Connect the variables to the cordesponding leaves
+    //c.SetBranchAddress("FileId", &FileId);
+    //c.SetBranchAddress("EvtNumber", &EvtNumber);
+    c.SetBranchAddress("X", &X);
+    c.SetBranchAddress("Y", &Y);
+    c.SetBranchAddress("MeanX", &MeanX);
+    c.SetBranchAddress("MeanY", &MeanY);
+    c.SetBranchAddress("Width", &Width);
+    c.SetBranchAddress("Length", &Length);
+    c.SetBranchAddress("CosDelta", &CosDelta);
+    c.SetBranchAddress("SinDelta", &SinDelta);
+    c.SetBranchAddress("M3Long", &M3Long);
+    c.SetBranchAddress("SlopeLong", &SlopeLong);
+    c.SetBranchAddress("Leakage1", &Leakage1);
+    c.SetBranchAddress("NumIslands", &NumIslands);
+    c.SetBranchAddress("NumUsedPixels", &NumUsedPixels);
+    c.SetBranchAddress("Size", &Size);
+    c.SetBranchAddress("Zd", &Zd);
+    c.SetBranchAddress("Energy", &Energy);
+
+    // Set some constants (they could be included in the database
+    // in the future)
+    double mm2deg = +0.0117193246260285378;
+    //double abberation = 1.02;
+
+    // -------------------- Source dependent parameter calculation -------------------
+
+    TH2F h2svse ("H_SvsE",      "", 26*3, 2.2, 4.8, 3*30, 1.3, 4.3);
+    TH2F h2est  ("H_EstVsMC",   "", 26*3, 2.2, 4.8, 26*3, 2.2, 4.8);
+    TH2F h2bias ("H_BiasLog",   "", 26*3, 2.2, 4.8, 100, -1, 1);
+    TH2F h2lin  ("H_BiasLin",   "", 26*3, 2.2, 4.8, 100, -1, 2);
+
+    TH2F h2size ("H_ResSize",   "", 26*3, 2.2, 4.8, 100, -1, 1);
+    TH2F h2dist ("H_ResDist",   "", 50,  0, 2.5,    100, -1, 1);
+    TH2F h2slope("H_ResSlope",  "", 50, -10, 10,    100, -1, 1);
+    TH2F h2zd   ("H_ResZd",     "", 90, 0,   90,    100, -1, 1);
+
+    h2svse.SetStats(kFALSE);
+    h2bias.SetStats(kFALSE);
+    h2lin.SetStats(kFALSE);
+    h2est.SetStats(kFALSE);
+    h2size.SetStats(kFALSE);
+    h2dist.SetStats(kFALSE);
+    h2slope.SetStats(kFALSE);
+    h2zd.SetStats(kFALSE);
+
+    // Loop over all events
+    TStopwatch clock;
+        // Loop over all wobble positions in the camera
+    for (int i=0; i<c.GetEntries(); i++)
+    {
+        // read the i-th event from the file
+        c.GetEntry(i);
+
+        // First calculate all cuts to speedup the analysis
+        double area = TMath::Pi()*Width*Length;
+
+        // The abberation correction does increase also Width and Length by 1.02
+
+        int angle = 0;
+
+        // -------------------- Source dependent parameter calculation -------------------
+
+        double cr = cos(angle*TMath::DegToRad());
+        double sr = sin(angle*TMath::DegToRad());
+
+        double px = cr*X-sr*Y;
+        double py = cr*Y+sr*X;
+
+        double dx = MeanX - px*1.022;
+        double dy = MeanY - py*1.022;
+
+        double norm = sqrt(dx*dx + dy*dy);
+        double dist = norm*mm2deg;
+
+        double lx = min(max((CosDelta*dy - SinDelta*dx)/norm, -1.), 1.);
+        double ly = min(max((CosDelta*dx + SinDelta*dy)/norm, -1.), 1.);
+
+        double alpha = asin(lx);
+        double sgn   = TMath::Sign(1., ly);
+
+        // ------------------------------- Application ----------------------------------
+
+        double m3l   = M3Long*sgn*mm2deg;
+        double slope = SlopeLong*sgn/mm2deg;
+
+        // --------------------------------- Analysis -----------------------------------
+
+        //double xi = 1.34723 + 0.15214 *slope + 0.970704*(1-1/(1+8.89826*Leakage1));
+        double xi = 1.340   + 0.0755  *slope + 1.67972 *(1-1/(1+4.86232*Leakage1));
+
+        double sign1 = m3l+0.07;
+        double sign2 = (dist-0.5)*7.2-slope;
+
+        double disp  = (sign1<0 || sign2<0 ? -xi : xi)*(1-Width/Length);
+
+        double thetasq = disp*disp + dist*dist - 2*disp*dist*sqrt(1-lx*lx);
+
+        if (thetasq<0.024)
+            continue;
+
+        bool cutq = NumIslands<3.5 && NumUsedPixels>5.5 && Leakage1<0.1;
+        if (!cutq)
+            continue;
+
+        bool cut0 = area < log10(Size)*898-1535;
+        if (!cut0)
+            continue;
+
+        h2svse.Fill(log10(Energy), log10(Size));
+
+        double step = 0.8;
+        double dd = dist<step ? -(dist-step)*0.4 : (dist-step)*0.65;
+
+        double energy = pow(10, 0.33 + log10(Size)/1.1 + Zd*0.02 + dd + slope*0.01);
+
+        h2est.Fill(log10(Energy),   log10(energy));
+        h2bias.Fill(log10(energy), log10(energy)-log10(Energy));
+        h2lin.Fill(log10(energy),  (energy-Energy)/Energy);
+
+        h2size.Fill (log10(Size), log10(energy)-log10(Energy));
+        h2dist.Fill (dist,        log10(energy)-log10(Energy));
+        h2slope.Fill(slope,       log10(energy)-log10(Energy));
+        h2zd.Fill(   Zd,          log10(energy)-log10(Energy));
+    }
+
+    clock.Print();
+
+    TF1 fx("f", "x", -100, 100);
+    TF1 f0("f", "0", -100, 100);
+
+    TCanvas *canv = new TCanvas("Canvas", "Energy estimation", 800, 900);
+    canv->Divide(2,4);
+
+    canv->cd(1);
+    gPad->SetGridy();
+    h2size.DrawCopy("colz");
+    f0.DrawCopy("same");
+
+    canv->cd(3);
+    gPad->SetGridy();
+    h2dist.DrawCopy("colz");
+    f0.DrawCopy("same");
+
+    canv->cd(5);
+    gPad->SetGridy();
+    h2slope.DrawCopy("colz");
+    f0.DrawCopy("same");
+
+    canv->cd(7);
+    gPad->SetGridy();
+    h2zd.DrawCopy("colz");
+    f0.DrawCopy("same");
+
+    canv->cd(2);
+    gPad->SetGridy();
+    h2svse.DrawCopy("colz");
+    fx.DrawCopy("same");
+
+    canv->cd(4);
+    gPad->SetGridy();
+    h2est.DrawCopy("colz");
+    fx.DrawCopy("same");
+
+    canv->cd(6);
+    gPad->SetGridy();
+    h2bias.DrawCopy("colz");
+    f0.DrawCopy("same");
+
+
+    canv->cd(8);
+
+    TGraph grlin;
+    TGraph grlog;
+
+    for (int x=1; x<h2bias.GetNbinsX(); x++)
+    {
+        TH1D *p = h2bias.ProjectionY("_py", x, x);
+        if (p->GetRMS()>0)
+            grlog.SetPoint(grlog.GetN(), h2bias.GetXaxis()->GetBinCenter(x), pow(10, p->GetRMS())-1);
+        delete p;
+
+        p = h2lin.ProjectionY("_py", x, x);
+        if (p->GetRMS()>0)
+            grlin.SetPoint(grlin.GetN(), h2lin.GetXaxis()->GetBinCenter(x), p->GetRMS());
+        delete p;
+    }
+
+    grlog.SetMarkerColor(kBlue);
+    grlin.SetMinimum(0);
+    grlin.DrawClone("A*");
+    grlog.DrawClone("*");
 }
 }}}
@@ -573 +1094 @@
 
 Note that this query can easily run 30min or more!
-
-== Optim Disp (Leakage==0) ==
-
-{{{#!Spoiler
-{{{#!cpp
-#include <iostream>
-#include <iomanip>
-
-#include <TMath.h>
-#include <TH2.h>
-#include <TChain.h>
-#include <TCanvas.h>
-
-void optimdisp2()
-{
-    // Create chain for the tree Result
-    // This is just easier than using TFile/TTree
-    TChain c("Result");
-
-    // Add the input file to the
-    c.AddFile("simulation.root");
-
-    // Define variables for all leaves to be accessed
-    // By definition rootifysql writes only doubles
-    double FileId, EvtNumber, X, Y, MeanX, MeanY, Width, Length, CosDelta, SinDelta,
-        M3Long, SlopeLong, Leakage1, SlopeSpreadWeighted, Size,
-        ConcCore, ConcCOG, NumIslands, NumUsedPixels, Zd;
-
-    // Connect the variables to the cordesponding leaves
-    c.SetBranchAddress("FileId", &FileId);
-    c.SetBranchAddress("EvtNumber", &EvtNumber);
-    c.SetBranchAddress("X", &X);
-    c.SetBranchAddress("Y", &Y);
-    c.SetBranchAddress("MeanX", &MeanX);
-    c.SetBranchAddress("MeanY", &MeanY);
-    c.SetBranchAddress("Width", &Width);
-    c.SetBranchAddress("Length", &Length);
-    c.SetBranchAddress("CosDelta", &CosDelta);
-    c.SetBranchAddress("SinDelta", &SinDelta);
-    c.SetBranchAddress("M3Long", &M3Long);
-    c.SetBranchAddress("SlopeLong", &SlopeLong);
-    c.SetBranchAddress("Leakage1", &Leakage1);
-    c.SetBranchAddress("NumIslands", &NumIslands);
-    c.SetBranchAddress("NumUsedPixels", &NumUsedPixels);
-    c.SetBranchAddress("Size", &Size);
-
-    // Set some constants (they could be included in the database
-    // in the future)
-    double mm2deg = +0.0117193246260285378;
-
-    // -------------------- Source dependent parameter calculation -------------------
-
-    TH2F hist("", "",
-              26, 1.300-0.0025, 1.430-0.0025,
-              20, 0.040-0.0025, 0.140-0.0025);
-
-    for (float xi0=1.30; xi0<1.43; xi0+=0.005)
-        for (float slope0=0.04; slope0<0.14; slope0+=0.005)
-        {
-            int cnt_on = 0;
-
-            for (int i=0; i<c.GetEntries(); i++)
-            {
-                // read the i-th event from the file
-                c.GetEntry(i);
-
-                // First calculate all cuts to speedup the analysis
-                double area = TMath::Pi()*Width*Length;
-
-                // The abberation correction does increase also Width and Length by 1.02
-
-                bool cutq = NumIslands<3.5 && NumUsedPixels>5.5 && Leakage1==0;
-                if (!cutq)
-                    continue;
-
-                bool cut0 = area < log10(Size)*898-1535;
-                if (!cut0)
-                    continue;
-
-                int angle = 0;
-
-                // -------------------- Source dependent parameter calculation -------------------
-
-                double cr = cos(angle*TMath::DegToRad());
-                double sr = sin(angle*TMath::DegToRad());
-
-                double px = cr*X-sr*Y;
-                double py = cr*Y+sr*X;
-
-                double dx = MeanX - px*1.02;
-                double dy = MeanY - py*1.02;
-
-                double norm = sqrt(dx*dx + dy*dy);
-                double dist = norm*mm2deg;
-
-                double lx = min(max((CosDelta*dy - SinDelta*dx)/norm, -1.), 1.);
-                double ly = min(max((CosDelta*dx + SinDelta*dy)/norm, -1.), 1.);
-
-                double alpha = asin(lx);
-                double sgn   = TMath::Sign(1., ly);
-
-                // ------------------------------- Application ----------------------------------
-
-                double m3l   = M3Long*sgn*mm2deg;
-                double slope = SlopeLong*sgn/mm2deg;
-
-                // --------------------------------- Analysis -----------------------------------
-
-                double xi = xi0 + slope0 *slope;
-
-                double sign1 = m3l+0.07;
-                double sign2 = (dist-0.5)*7.2-slope;
-
-                double disp  = (sign1<0 || sign2<0 ? -xi : xi)*(1-Width/Length);
-
-                double thetasq = disp*disp + dist*dist - 2*disp*dist*sqrt(1-lx*lx);
-
-                if (thetasq<0.024)
-                    cnt_on++;
-            }
-
-            hist.Fill(xi0, slope0, cnt_on);
-
-            cout << xi0 << " " << slope0 << " " << cnt_on << endl;
-        }
-
-    TCanvas *canv = new TCanvas;
-    canv->SetTopMargin(0.01);
-
-    hist.SetStats(kFALSE);
-    hist.SetContour(100);
-    hist.SetXTitle("Xi_{0} / deg");
-    hist.SetYTitle("Slope_{0} #dot ns/deg^{2} ");
-    hist.GetYaxis()->SetTitleOffset(1.2);
-    hist.DrawCopy("colz cont2");
-
-}
-}}}
-}}}
-
-Results in
-
-[[Image(Result.png)]]
-
-== Optim Disp (Leakage>0) ==
-
-{{{#!Spoiler
-{{{#!cpp
-#include <iostream>
-#include <iomanip>
-
-#include <TMath.h>
-#include <TH2.h>
-#include <TChain.h>
-#include <TCanvas.h>
-
-void optimdisp3()
-{
-    // Create chain for the tree Result
-    // This is just easier than using TFile/TTree
-    TChain c("Result");
-
-    // Add the input file to the
-    c.AddFile("simulation.root");
-
-    // Define variables for all leaves to be accessed
-    // By definition rootifysql writes only doubles
-    double FileId, EvtNumber, X, Y, MeanX, MeanY, Width, Length, CosDelta, SinDelta,
-        M3Long, SlopeLong, Leakage1, SlopeSpreadWeighted, Size,
-        ConcCore, ConcCOG, NumIslands, NumUsedPixels, Zd;
-
-    // Connect the variables to the cordesponding leaves
-    c.SetBranchAddress("FileId", &FileId);
-    c.SetBranchAddress("EvtNumber", &EvtNumber);
-    c.SetBranchAddress("X", &X);
-    c.SetBranchAddress("Y", &Y);
-    c.SetBranchAddress("MeanX", &MeanX);
-    c.SetBranchAddress("MeanY", &MeanY);
-    c.SetBranchAddress("Width", &Width);
-    c.SetBranchAddress("Length", &Length);
-    c.SetBranchAddress("CosDelta", &CosDelta);
-    c.SetBranchAddress("SinDelta", &SinDelta);
-    c.SetBranchAddress("M3Long", &M3Long);
-    c.SetBranchAddress("SlopeLong", &SlopeLong);
-    c.SetBranchAddress("Leakage1", &Leakage1);
-    c.SetBranchAddress("NumIslands", &NumIslands);
-    c.SetBranchAddress("NumUsedPixels", &NumUsedPixels);
-    //c.SetBranchAddress("SlopeSpreadWeighted", &SlopeSpreadWeighted);
-    c.SetBranchAddress("Size", &Size);
-    c.SetBranchAddress("Zd", &Zd);
-    //c.SetBranchAddress("ConcCOG", &ConcCOG);
-    //c.SetBranchAddress("ConcCore", &ConcCore);
-
-    // Set some constants (they could be included in the database
-    // in the future)
-    double mm2deg = +0.0117193246260285378;
-    //double abberation = 1.02;
-
-    // -------------------- Source dependent parameter calculation -------------------
-
-    TH2F hist("", "",
-              32, 0.800-0.0250, 2.400-0.0250,
-              60, 3.000-0.0250, 6.000-0.0250);
-
-
-    for (float l0=0.8; l0<2.4; l0+=0.05)
-        for (float l1=3; l1<6; l1+=0.05)
-        {
-            int cnt_on = 0;
-
-            for (int i=0; i<c.GetEntries(); i++)
-            {
-                // read the i-th event from the file
-                c.GetEntry(i);
-
-                // First calculate all cuts to speedup the analysis
-                double area = TMath::Pi()*Width*Length;
-
-                // The abberation correction does increase also Width and Length by 1.02
-
-                bool cutq = NumIslands<3.5 && NumUsedPixels>5.5 && Leakage1>0;
-                if (!cutq)
-                    continue;
-
-                bool cut0 = area < log10(Size)*898-1535;
-                if (!cut0)
-                    continue;
-
-                int angle = 0;
-
-                // -------------------- Source dependent parameter calculation -------------------
-
-                double cr = cos(angle*TMath::DegToRad());
-                double sr = sin(angle*TMath::DegToRad());
-
-                double px = cr*X-sr*Y;
-                double py = cr*Y+sr*X;
-
-                double dx = MeanX - px*1.02;
-                double dy = MeanY - py*1.02;
-
-                double norm = sqrt(dx*dx + dy*dy);
-                double dist = norm*mm2deg;
-
-                double lx = min(max((CosDelta*dy - SinDelta*dx)/norm, -1.), 1.);
-                double ly = min(max((CosDelta*dx + SinDelta*dy)/norm, -1.), 1.);
-
-                double alpha = asin(lx);
-                double sgn   = TMath::Sign(1., ly);
-
-                // ------------------------------- Application ----------------------------------
-
-                double m3l   = M3Long*sgn*mm2deg;
-                double slope = SlopeLong*sgn/mm2deg;
-
-                // --------------------------------- Analysis -----------------------------------
-
-                double xi = 1.36 + 0.085*slope + l0 *(1-1/(1+l1*Leakage1));
-
-                // Optim: 1.36 + 0.085
-
-                double sign1 = m3l+0.07;
-                double sign2 = (dist-0.5)*7.2-slope;
-
-                double disp  = (sign1<0 || sign2<0 ? -xi : xi)*(1-Width/Length);
-
-                double thetasq = disp*disp + dist*dist - 2*disp*dist*sqrt(1-lx*lx);
-
-                if (thetasq<0.024)
-                    cnt_on++;
-            }
-
-            hist.Fill(l0, l1, cnt_on);
-
-            cout << l0 << " " << l1 << " " << cnt_on << endl;
-        }
-
-    TCanvas *canv = new TCanvas;
-    canv->SetTopMargin(0.01);
-
-    hist.SetStats(kFALSE);
-    hist.SetContour(100);
-    hist.SetXTitle("L_{0} / deg");
-    hist.SetYTitle("L_{1} / deg");
-    hist.GetYaxis()->SetTitleOffset(1.2);
-    hist.DrawCopy("colz cont2");
-
-}
-}}}
-}}}
-
-Results in
-
-[[Image(Result2.png)]]
-
-== Energy Optimization == 
-
-{{{#!cpp
-#include <iostream>
-#include <iomanip>
-
-#include <TMath.h>
-#include <TF1.h>
-#include <TH1.h>
-#include <TH2.h>
-#include <TProfile.h>
-#include <TChain.h>
-#include <TGraph.h>
-#include <TCanvas.h>
-#include <TStopwatch.h>
-
-void optimenergy()
-{
-    // Create chain for the tree Result
-    // This is just easier than using TFile/TTree
-    TChain c("Result");
-
-    // Add the input file to the
-    c.AddFile("simulation.root");
-
-    // Define variables for all leaves to be accessed
-    // By definition rootifysql writes only doubles
-    double X, Y, MeanX, MeanY, Width, Length, CosDelta, SinDelta,
-        M3Long, SlopeLong, Leakage1, SlopeSpreadWeighted, Size,
-        ConcCore, ConcCOG, NumIslands, NumUsedPixels, Zd, Energy;
-
-    // Connect the variables to the cordesponding leaves
-    //c.SetBranchAddress("FileId", &FileId);
-    //c.SetBranchAddress("EvtNumber", &EvtNumber);
-    c.SetBranchAddress("X", &X);
-    c.SetBranchAddress("Y", &Y);
-    c.SetBranchAddress("MeanX", &MeanX);
-    c.SetBranchAddress("MeanY", &MeanY);
-    c.SetBranchAddress("Width", &Width);
-    c.SetBranchAddress("Length", &Length);
-    c.SetBranchAddress("CosDelta", &CosDelta);
-    c.SetBranchAddress("SinDelta", &SinDelta);
-    c.SetBranchAddress("M3Long", &M3Long);
-    c.SetBranchAddress("SlopeLong", &SlopeLong);
-    c.SetBranchAddress("Leakage1", &Leakage1);
-    c.SetBranchAddress("NumIslands", &NumIslands);
-    c.SetBranchAddress("NumUsedPixels", &NumUsedPixels);
-    c.SetBranchAddress("Size", &Size);
-    c.SetBranchAddress("Zd", &Zd);
-    c.SetBranchAddress("Energy", &Energy);
-
-    // Set some constants (they could be included in the database
-    // in the future)
-    double mm2deg = +0.0117193246260285378;
-    //double abberation = 1.02;
-
-    // -------------------- Source dependent parameter calculation -------------------
-
-    TH2F h2svse ("H_SvsE",      "", 26*3, 2.2, 4.8, 3*30, 1.3, 4.3);
-    TH2F h2est  ("H_EstVsMC",   "", 26*3, 2.2, 4.8, 26*3, 2.2, 4.8);
-    TH2F h2bias ("H_BiasLog",   "", 26*3, 2.2, 4.8, 100, -1, 1);
-    TH2F h2lin  ("H_BiasLin",   "", 26*3, 2.2, 4.8, 100, -1, 2);
-
-    TH2F h2size ("H_ResSize",   "", 26*3, 2.2, 4.8, 100, -1, 1);
-    TH2F h2dist ("H_ResDist",   "", 50,  0, 2.5,    100, -1, 1);
-    TH2F h2slope("H_ResSlope",  "", 50, -10, 10,    100, -1, 1);
-    TH2F h2zd   ("H_ResZd",     "", 90, 0,   90,    100, -1, 1);
-
-    h2svse.SetStats(kFALSE);
-    h2bias.SetStats(kFALSE);
-    h2lin.SetStats(kFALSE);
-    h2est.SetStats(kFALSE);
-    h2size.SetStats(kFALSE);
-    h2dist.SetStats(kFALSE);
-    h2slope.SetStats(kFALSE);
-    h2zd.SetStats(kFALSE);
-
-    // Loop over all events
-    TStopwatch clock;
-        // Loop over all wobble positions in the camera
-    for (int i=0; i<c.GetEntries(); i++)
-    {
-        // read the i-th event from the file
-        c.GetEntry(i);
-
-        // First calculate all cuts to speedup the analysis
-        double area = TMath::Pi()*Width*Length;
-
-        // The abberation correction does increase also Width and Length by 1.02
-
-        int angle = 0;
-
-        // -------------------- Source dependent parameter calculation -------------------
-
-        double cr = cos(angle*TMath::DegToRad());
-        double sr = sin(angle*TMath::DegToRad());
-
-        double px = cr*X-sr*Y;
-        double py = cr*Y+sr*X;
-
-        double dx = MeanX - px*1.022;
-        double dy = MeanY - py*1.022;
-
-        double norm = sqrt(dx*dx + dy*dy);
-        double dist = norm*mm2deg;
-
-        double lx = min(max((CosDelta*dy - SinDelta*dx)/norm, -1.), 1.);
-        double ly = min(max((CosDelta*dx + SinDelta*dy)/norm, -1.), 1.);
-
-        double alpha = asin(lx);
-        double sgn   = TMath::Sign(1., ly);
-
-        // ------------------------------- Application ----------------------------------
-
-        double m3l   = M3Long*sgn*mm2deg;
-        double slope = SlopeLong*sgn/mm2deg;
-
-        // --------------------------------- Analysis -----------------------------------
-
-        //double xi = 1.34723 + 0.15214 *slope + 0.970704*(1-1/(1+8.89826*Leakage1));
-        double xi = 1.340   + 0.0755  *slope + 1.67972 *(1-1/(1+4.86232*Leakage1));
-
-        double sign1 = m3l+0.07;
-        double sign2 = (dist-0.5)*7.2-slope;
-
-        double disp  = (sign1<0 || sign2<0 ? -xi : xi)*(1-Width/Length);
-
-        double thetasq = disp*disp + dist*dist - 2*disp*dist*sqrt(1-lx*lx);
-
-        if (thetasq<0.024)
-            continue;
-
-        bool cutq = NumIslands<3.5 && NumUsedPixels>5.5 && Leakage1<0.1;
-        if (!cutq)
-            continue;
-
-        bool cut0 = area < log10(Size)*898-1535;
-        if (!cut0)
-            continue;
-
-        h2svse.Fill(log10(Energy), log10(Size));
-
-        double step = 0.8;
-        double dd = dist<step ? -(dist-step)*0.4 : (dist-step)*0.65;
-
-        double energy = pow(10, 0.33 + log10(Size)/1.1 + Zd*0.02 + dd + slope*0.01);
-
-        h2est.Fill(log10(Energy),   log10(energy));
-        h2bias.Fill(log10(energy), log10(energy)-log10(Energy));
-        h2lin.Fill(log10(energy),  (energy-Energy)/Energy);
-
-        h2size.Fill (log10(Size), log10(energy)-log10(Energy));
-        h2dist.Fill (dist,        log10(energy)-log10(Energy));
-        h2slope.Fill(slope,       log10(energy)-log10(Energy));
-        h2zd.Fill(   Zd,          log10(energy)-log10(Energy));
-    }
-
-    clock.Print();
-
-    TF1 fx("f", "x", -100, 100);
-    TF1 f0("f", "0", -100, 100);
-
-    TCanvas *canv = new TCanvas("Canvas", "Energy estimation", 800, 900);
-    canv->Divide(2,4);
-
-    canv->cd(1);
-    gPad->SetGridy();
-    h2size.DrawCopy("colz");
-    f0.DrawCopy("same");
-
-    canv->cd(3);
-    gPad->SetGridy();
-    h2dist.DrawCopy("colz");
-    f0.DrawCopy("same");
-
-    canv->cd(5);
-    gPad->SetGridy();
-    h2slope.DrawCopy("colz");
-    f0.DrawCopy("same");
-
-    canv->cd(7);
-    gPad->SetGridy();
-    h2zd.DrawCopy("colz");
-    f0.DrawCopy("same");
-
-    canv->cd(2);
-    gPad->SetGridy();
-    h2svse.DrawCopy("colz");
-    fx.DrawCopy("same");
-
-    canv->cd(4);
-    gPad->SetGridy();
-    h2est.DrawCopy("colz");
-    fx.DrawCopy("same");
-
-    canv->cd(6);
-    gPad->SetGridy();
-    h2bias.DrawCopy("colz");
-    f0.DrawCopy("same");
-
-
-    canv->cd(8);
-
-    TGraph grlin;
-    TGraph grlog;
-
-    for (int x=1; x<h2bias.GetNbinsX(); x++)
-    {
-        TH1D *p = h2bias.ProjectionY("_py", x, x);
-        if (p->GetRMS()>0)
-            grlog.SetPoint(grlog.GetN(), h2bias.GetXaxis()->GetBinCenter(x), pow(10, p->GetRMS())-1);
-        delete p;
-
-        p = h2lin.ProjectionY("_py", x, x);
-        if (p->GetRMS()>0)
-            grlin.SetPoint(grlin.GetN(), h2lin.GetXaxis()->GetBinCenter(x), p->GetRMS());
-        delete p;
-    }
-
-    grlog.SetMarkerColor(kBlue);
-    grlin.SetMinimum(0);
-    grlin.DrawClone("A*");
-    grlog.DrawClone("*");
-}
-}}}