Changes between Version 18 and Version 19 of DatabaseBasedAnalysis/Examples

Timestamp:

08/15/18 15:46:04 (7 years ago)

Author:

tbretz

Comment:

--

DatabaseBasedAnalysis/Examples

@@ -573 +573 @@
 
 Note that this query can easily run 30min or more!
+
+== 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("*");
+}
+}}}