Changes between Version 1 and Version 2 of DatabaseBasedAnalysis/Examples

Timestamp:

08/12/18 13:53:59 (7 years ago)

Author:

tbretz

Comment:

--

DatabaseBasedAnalysis/Examples

@@ -32 +32 @@
    fR750Cor>0.9*fR750Ref
 }}}
+
+== Optimize Disp ==
+
+The following root-macro (to be compiled!)
+
+{{{!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>
+#include <TStopwatch.h>
+
+void optimdisp()
+{
+    // 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, 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);
+    c.SetBranchAddress("Zd", &Zd);
+
+    // 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     h2slope("H_VsSlope",   "", 75, -8,     7,    100, -2.5, 1.5);
+    TProfile p2slope("P_VsSlope",   "", 75, -8,     7);
+    TH2F     h2leak( "H_VsLeakage", "", 75,  0,     0.15, 100, -2.5, 1.5);
+    TH2F     h2m3l(  "H_M3long",    "", 75, -0.2,   0.6,  100, -2.5, 1.5);
+    TH2F     h2zd(   "H_Zd",        "", 30, 30,    60,    100, -2.5, 1.5);
+    TH2F     h2ni(   "H_NumIsl",    "", 10,  0.5,  10.5,  100, -2.5, 1.5);
+    TH2F     h2np(   "H_NumPix",    "", 10,  0.5, 100.5,  100, -2.5, 1.5);
+    TH2F     h2size( "H_Size",      "", 30,  1.5,   4.5,  100, -2.5, 1.5);
+
+    // 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;
+
+        bool cutq = NumIslands<3.5 && NumUsedPixels>5.5 && Leakage1<0.1;
+        if (!cutq)
+            continue;
+
+        bool cut0 = area < log10(Size)*898-1535;
+        if (!cut0)
+            continue;
+
+        // -------------------- Source dependent parameter calculation -------------------
+
+        int angle = 0;
+
+        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_old = 1.39252 + 0.154247*slope + 1.67972 *(1-1/(1+4.86232*Leakage1));
+        double xi_new = 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_old  = (sign1<0 || sign2<0 ? -xi_old : xi_old)*(1-Width/Length);
+        double disp_new  = (sign1<0 || sign2<0 ? -xi_new : xi_new)*(1-Width/Length);
+
+        double thetasq_old = disp_old*disp_old + dist*dist - 2*disp_old*dist*sqrt(1-lx*lx);
+        double thetasq_new = disp_new*disp_new + dist*dist - 2*disp_new*dist*sqrt(1-lx*lx);
+
+        // Fill the on- and off-histograms
+        hold.Fill(thetasq_old);
+        hnew.Fill(thetasq_new);
+
+        double residual = xi_new-dist/(1-Width/Length);
+
+        h2slope.Fill(slope,        residual);
+        p2slope.Fill(slope,        residual);
+
+        h2leak.Fill(Leakage1,      residual);
+        h2m3l.Fill( m3l,           residual);
+        h2zd.Fill(  Zd,            residual);
+        h2ni.Fill(  NumIslands,    residual);
+        h2np.Fill(  NumUsedPixels, residual);
+        h2size.Fill(log10(Size),   residual);
+    }
+
+    cout << hnew.GetBinContent(1) << endl;
+
+    clock.Print();
+
+    TCanvas *canv = new TCanvas;
+    canv->Divide(2,2);
+
+    canv->cd(1);
+    gPad->SetGridy();
+    h2slope.DrawCopy("colz");
+    p2slope.DrawCopy("same");
+
+    canv->cd(2);
+    gPad->SetGridy();
+    h2leak.DrawCopy("colz");
+
+    canv->cd(3);
+    gPad->SetGridy();
+    //h2m3l.DrawCopy("colz");
+    h2zd.DrawCopy("colz");
+    //h2ni.DrawCopy("colz");
+    //h2np.DrawCopy("colz");
+    //h2size.DrawCopy("colz");
+
+    canv->cd(4);
+    // Plot the result
+    hold.SetLineColor(kRed);
+    hnew.SetMinimum(0);
+    hnew.DrawCopy();
+    hold.DrawCopy("same");
+}
+}}}