Changes between Version 29 and Version 30 of DatabaseBasedAnalysis/Examples

Timestamp:

09/12/18 16:04:34 (7 years ago)

Author:

tbretz

Comment:

--

DatabaseBasedAnalysis/Examples

@@ -1286 +1286 @@
 
 Note that this query can easily run 30min or more!
+
+== Theta-Square Plot ==
+
+There are several ways ot prepare a theta-square plot. If you want to change the cuts later on, you will have to download all the data and do it offline. However, this is time- and bandwidth-hungry. If you know your cuts already, an easy easy solution is to let mysql create the plot for you and just use a local program to plot it. It significantly decreases bandwidth and space requirement (maybe from several GB to just a few byte) and it significantly increase CPU requirement (because data does not have to be converted to ascii representation and back and compressed for the transfer but it fully processed in memory).
+
+Assuming you produce a root file which contains three branches '''Bin''' (the number of the bin), '''Weight''' (-0.2 for background and 1 for signal) and '''Counts''' for the number of bin entries. You can for example run the following query using rootifysql (see DatabaseBasedAnalysis/rootifysql). Note that you might want to adapt the theta-sq cut your are going to use and the number of bins into which you sort the data up to the theta-square cut.
+
+{{{!Spoiler
+{{{#!sql
+WITH Table7 AS
+(
+    WITH Table6 AS
+    (
+        WITH Table5 AS
+        ( 
+            WITH Table4 AS
+            (
+                WITH Table3 AS
+                (
+                    WITH Table2 AS
+                    (
+                        WITH Table1 AS
+                        (
+                            WITH Table0 AS
+                            (
+                                SELECT -- 0
+                                    Weight,
+                                    Size,
+                                    NumUsedPixels,
+                                    NumIslands,
+                                    Leakage1,
+                                    MeanX,
+                                    MeanY,
+                                    CosDelta,
+                                    SinDelta,
+                                    M3Long,
+                                    SlopeLong,
+                                    Width/Length      AS WdivL,
+                                    PI()*Width*Length AS Area,
+                                    cosa*X - sina*Y   AS PX,
+                                    cosa*Y + sina*X   AS PY
+                                FROM RunInfo
+                                LEFT JOIN Images   USING (FileId)
+                                LEFT JOIN Position USING (FileId, EvtNumber)
+                                CROSS JOIN Wobble
+                                WHERE
+                                    FileId BETWEEN 131101000 AND 131107000
+                                AND
+                                    fSourceKey=5
+                                AND
+                                    fRunTypeKey=1
+                                AND
+                                    fZenithDistanceMax<35
+                                AND
+                                    fR750Cor>0.9e0*fR750Ref
+                                AND
+                                    NumUsedPixels>5.5
+                                AND
+                                    Leakage1<0.1
+
+                            ) -- AS Table0
+
+                            SELECT -- 1
+                                Weight, CosDelta, SinDelta, M3Long, SlopeLong, Leakage1, WdivL,
+                                MeanX - PX*1.02e0 AS DX,
+                                MeanY - PY*1.02e0 AS DY
+                            FROM
+                                Table0
+                            WHERE
+                                Area < LOG10(Size)*898-1535
+
+                        ) -- AS Table1
+
+                        SELECT -- 2
+                            Weight, CosDelta, SinDelta, DX, DY, M3Long, SlopeLong, Leakage1, WdivL,
+                            SQRT(DX*DX + DY*DY) AS Norm
+                        FROM
+                            Table1
+
+                    ) -- AS Table2
+
+                    SELECT -- 3
+                        Weight, M3Long, SlopeLong, Leakage1, WdivL, Norm,
+                        LEAST(GREATEST((CosDelta*DY - SinDelta*DX)/Norm, -1), 1) AS LX,
+                        SIGN(CosDelta*DX + SinDelta*DY) AS Sign
+                    FROM
+                        Table2
+
+                ) -- AS Table3
+
+                SELECT -- 4
+                    Weight, Leakage1, WdivL, LX,
+                    Norm          *0.0117193246260285378e0 AS Dist,
+                    M3Long   *Sign*0.0117193246260285378e0 AS M3L,
+                    SlopeLong*Sign/0.0117193246260285378e0 AS Slope
+                FROM
+                        Table3
+
+            ) -- AS Table4
+
+            SELECT -- 5
+                Weight, WdivL, Dist, LX, M3L, Slope,
+                1.299e0 + 0.0632e0*Slope + 1.67972e0*(1-1/(1+4.86232e0*Leakage1)) AS Xi
+                FROM
+                Table4
+
+        ) -- AS Table5
+
+        SELECT -- 6
+            Weight, Dist, LX,
+            IF (M3L<-0.07 || (Dist-0.5e0)*7.2e0-Slope<0, -Xi, Xi) * (1-WdivL) AS Disp
+        FROM
+            Table5
+
+    ) -- AS Table6
+
+    SELECT -- 7
+        Weight,
+        (Disp*Disp + Dist*Dist - 2*Disp*Dist*SQRT(1-LX*LX)) AS ThetaSq
+    FROM
+        Table6
+
+) -- AS Table8
+
+SELECT -- 9
+
+    /* The number 0.024 corresponds to the theta-sq cut you want to apply.  */
+    /* The number 3 is the number of bins into which the values smaller     */
+    /* the theta-sq cut are divide. We are basically calculating the bin-   */
+    /* number into which the event is sorted later. 0.024/3 is the bin-     */
+    /* width.                                                               */
+
+    FLOOR(ThetaSq/(0.024/3)) AS Bin,
+    Weight,
+    COUNT(*) AS `Count`
+FROM
+    Table7
+GROUP BY
+    Bin, Weight
+
+/* If no printed on the console, strictly speaking, this is not necessary */
+ORDER BY
+    Bin, Weight
+}}}
+}}}
+
+In a second step, you might want to produce a theta-square plot with your preferred language (here: root/C++). '''It is very important that you ensure that your plotting routine is consistent with the the histogram data produced by the former SQL query in terms of bin-width!''' Here is an example root macro to plot the data:
+
+{{{#!Spoiler
+{{{#!C++
+// ============= Significance according to Li&Ma ===============
+Double_t LiMa(Double_t s, Double_t b, Double_t alpha)
+{
+    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 ? -1 : TMath::Sqrt((l+m)*2);
+}
+
+void plot_thetasq()
+{
+    // ====================== Setup the plot ===================
+
+    double cut   = 0.024;      // Theta-square cut 
+    int    icut  = 3;          // Number of bins at which the theta-square cut should be displayed
+    double max   = 0.98;       // maximum to be displayed on the x-axis (in degree^2)
+    double width = cut/icut;   // width of a single bin (here: a third of the theta-square cut)
+
+    const char *title     = "FACT";
+    const char *watermark = "preliminary";
+
+    // ================= Initialize root chain =================
+
+    TChain chain("Result");
+    chain.AddFile("thetasq-hist.root");
+
+    // Define variables for all leaves to be accessed
+    // By definition rootifysql writes only doubles
+    double Bin, Weight, Count;
+
+    // Connect the variables to the cordesponding leaves
+    chain.SetBranchAddress("Bin",    &Bin);
+    chain.SetBranchAddress("Weight", &Weight);
+    chain.SetBranchAddress("Count",  &Count);
+
+    // ==================== Setting up Histogram ===============
+
+    int num = floor(max/width);
+
+    // Create a histogram for on- and off-data
+    TH1F hon ("on",  "", num, 0, num*width);
+    TH1F hoff("off", "", num, 0, num*width);
+
+    hon.Sumw2();
+    hoff.Sumw2();
+
+    // Loop over all events
+    for (int i=0; i<chain.GetEntries(); i++)
+    {
+        // read the i-th event from the file
+        chain.GetEntry(i);
+
+        if (Weight>0)
+        {
+            hon.SetBinContent((int)Bin+1, Count);
+            hon.SetBinError(  (int)Bin+1, sqrt(Count));
+        }
+        else
+        {
+            hoff.SetBinContent((int)Bin+1, -Weight*Count);
+            hoff.SetBinError(  (int)Bin+1, sqrt(-Weight*Count));
+        }
+    }
+
+    // ========================= Result ========================
+
+    double signal     = hon.Integral(1, icut);
+    double background = hoff.Integral(1, icut);
+
+    double significance = LiMa(signal, background*5, 0.2);
+
+    cout << "Signal: " << signal       << endl;
+    cout << "Backgr: " << background   << endl;
+    cout << "Signif: " << significance << endl;
+
+    // ======================= Plotting ========================
+
+    TCanvas *c = new TCanvas;
+
+    c->SetBorderMode(0);
+    c->SetFrameBorderMode(0);
+    c->SetFillColor(kWhite);
+
+    c->SetLeftMargin(0.12);
+    c->SetRightMargin(0.01);
+    c->SetBottomMargin(0.16);
+    c->SetTopMargin(0.18);
+
+    c->SetGridy();
+
+    hon.SetStats(kFALSE);
+    hon.SetMinimum(0);
+    hon.SetLineColor(kBlack);
+    hon.SetMarkerColor(kBlack);
+    hon.SetLineWidth(2);
+
+    hoff.GetXaxis()->SetLabelSize(0.06);
+    hoff.GetXaxis()->SetTitleSize(0.06);
+    hoff.GetXaxis()->SetTitleOffset(0.95);
+    hoff.GetYaxis()->SetLabelSize(0.06);
+    hoff.GetYaxis()->SetTitleSize(0.06);
+    hoff.GetYaxis()->CenterTitle();
+    hoff.SetXTitle("#vartheta^2");
+    hoff.SetYTitle("Counts");
+    hoff.SetTitleSize(0.07);
+    hoff.SetTitle("");
+
+    hoff.SetStats(kFALSE);
+    hoff.SetLineWidth(2);
+    hoff.SetLineColor(kBlack);
+    hoff.SetMarkerColor(kBlack);
+    hoff.SetFillColor(17);
+
+    hoff.SetMinimum(0);
+    hoff.SetMaximum(TMath::Max(hon.GetMaximum(), hoff.GetMaximum())*1.1);
+
+    hoff.DrawCopy("bar");
+    hon.DrawCopy("same");
+
+    // draw a preliminary tag
+    TLatex text;
+    text.SetTextColor(kWhite);
+    text.SetBit(TLatex::kTextNDC);
+    text.SetTextSize(0.07);
+    text.SetTextAngle(2.5);
+
+    TString wm(watermark);
+    text.DrawLatex(0.45, 0.2, wm);
+
+    // draw line showing cut
+    TLine line;
+    line.SetLineColor(14);
+    line.SetLineStyle(7);
+    line.DrawLine(icut*width, 0, icut*width, hoff.GetMaximum());
+
+    // Add a title above the plot
+    TPaveText *pave=new TPaveText(0.12, 0.83, 0.99, 0.975, "blNDC");
+    pave->SetBorderSize(1);
+    pave->SetLabel(title);
+
+    TText *ptxt = pave->AddText(" ");
+    ptxt->SetTextAlign(20);
+
+    ptxt = pave->AddText(Form("Significance  %.1f#sigma,  off-scale  0.2",
+                              significance));
+    ptxt->SetTextAlign(21);
+
+    ptxt = pave->AddText(Form("%.1f excess events,  %.1f background events",
+                              signal-background, background));
+    ptxt->SetTextAlign(21);
+
+    pave->SetBit(kCanDelete);
+    pave->Draw("br");
+}
+}}}
+}}}
+
+This is how the result will look like:
+
+[[Image(ThetaSquare.png)]]
+
+
+