Changes between Version 4 and Version 5 of DatabaseBasedAnalysis/Examples

Timestamp:

08/12/18 14:02:23 (6 years ago)

Author:

tbretz

Comment:

--

DatabaseBasedAnalysis/Examples

@@ +1 @@
+[[TOC]]
+
+
+== Optimize Area-cut ==
+
 Assume you have written a large data set of good-weather, low zenith-distance Crab data with the following query to a file. The following will give you a few example of how you can make use of that in your analysis. To speed up things or decrease the file size (715MB), you might want to add additional cuts depending on your analysis. If you want to apply a Theta-square cut, please refer to DatabaseBasedAnalysis#Dataretrieval.
 
@@ -33 +38 @@
 }}}
 
+{{{!#cpp
+#include <iostream>
+#include <iomanip>
+
+#include <TMath.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TProfile.h>
+#include <TChain.h>
+#include <TGraph.h>
+#include <TCanvas.h>
+
+Double_t LiMa(Double_t s, Double_t b, Double_t alpha=0.2)
+{
+    const Double_t sum = s+b;
+
+    if (s<0 || b<0 || alpha<=0)
+        return -1;
+
+    const Double_t l = s==0 ? 0 : s*TMath::Log(s/sum*(alpha+1)/alpha);
+    const Double_t m = b==0 ? 0 : b*TMath::Log(b/sum*(alpha+1)      );
+
+    return l+m<0 ? 0 : TMath::Sqrt((l+m)*2);
+}
+
+void optimarea()
+{
+    // 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("crab2018.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, Size,
+        ConcCore, ConcCOG, NumIslands, NumUsedPixels, Zd;
+
+    // Connect the variables to the cordesponding leaves
+    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 -------------------
+
+    // Create a histogram for on- and off-data
+    TH1F hold("old",  "", 100, 0, 1);
+    TH1F hnew("onew", "", 100, 0, 1);
+
+    TH2F h2area_on( "H_AvsS_on",   "", 24, 1.5, 4.5, 200, -2, 0);
+    TH2F h2area_of( "H_AvsS_of",   "", 24, 1.5, 4.5, 200, -2, 0);
+    TH2F h2area_sig("H_AvsS_sig",  "", 24, 1.5, 4.5, 200, -2, 0);
+
+    // Loop over all events
+    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;
+
+        bool cutq = NumIslands<3.5 && NumUsedPixels>5.5 && Leakage1<0.1;
+        if (!cutq)
+            continue;
+
+        // Loop over all wobble positions in the camera
+        for (int angle=0; angle<360; angle+=60)
+        {
+
+        // -------------------- 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.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_old*disp_old + dist*dist - 2*disp_old*dist*sqrt(1-lx*lx);
+
+            if (thetasq>0.024)
+                continue;
+
+            if (angle==0)
+                h2area_on.Fill(log10(Size), log10(area*mm2deg*mm2deg));
+            else
+                h2area_of.Fill(log10(Size), log10(area*mm2deg*mm2deg));
+        }
+    }
+
+    for (int x=1; x<=h2area_sig.GetNbinsX(); x++)
+        for (int y=1; y<=h2area_sig.GetNbinsY(); y++)
+        {
+            double on = h2area_on.Integral(x, x, 1, y);
+            double of = h2area_of.Integral(x, x, 1, y);
+
+            h2area_sig.SetBinContent(x, y, Error(on, of));
+        }
+
+    TGraph gmax;
+    for (int x=1; x<=h2area_sig.GetNbinsX(); x++)
+    {
+        TH1D *p = h2area_sig.ProjectionY("_py", x, x);
+        double max = p->GetXaxis()->GetBinCenter(p->GetMaximumBin());
+        delete p;
+
+        gmax.SetPoint(gmax.GetN(),
+                      h2area_sig.GetXaxis()->GetBinCenter(x),
+                      max);
+    }
+
+    clock.Print();
+
+    h2area_on.SetStats(kFALSE);
+    h2area_of.SetStats(kFALSE);
+    h2area_sig.SetStats(kFALSE);
+
+    TCanvas *canv = new TCanvas;
+    canv->Divide(2,2);
+
+    canv->cd(1);
+    h2area_on.DrawCopy("colz");
+    gmax.DrawClone("*");
+
+    canv->cd(2);
+    h2area_sig.DrawCopy("colz");
+    gmax.DrawClone("*");
+
+    canv->cd(3);
+    h2area_of.DrawCopy("colz");
+    gmax.DrawClone("*");
+
+    canv->cd(4);
+}
+}}}
+
 == Optimize Disp ==