Changes between Version 25 and Version 26 of DatabaseBasedAnalysis/Examples


Ignore:
Timestamp:
08/17/18 19:01:30 (7 years ago)
Author:
tbretz
Comment:

--

Legend:

Unmodified
Added
Removed
Modified
  • DatabaseBasedAnalysis/Examples

    v25 v26  
    720720
    721721[[Image(Result2.png)]]
     722
     723== Optimize Disp (on Crab Data) ==
     724
     725This is the same analysis using the Excess-Error on Crab Data.
     726
     727{{{#!Spoiler
     728{{{#!cpp
     729#include <iostream>
     730#include <iomanip>
     731
     732#include <TMath.h>
     733#include <TH2.h>
     734#include <TChain.h>
     735#include <TCanvas.h>
     736#include <TGraph.h>
     737#include <TStopwatch.h>
     738
     739Double_t LiMa(Double_t s, Double_t b, Double_t alpha=0.2)
     740{
     741    const Double_t sum = s+b;
     742
     743    if (s<0 || b<0 || alpha<=0)
     744        return -1;
     745
     746    const Double_t l = s==0 ? 0 : s*TMath::Log(s/sum*(alpha+1)/alpha);
     747    const Double_t m = b==0 ? 0 : b*TMath::Log(b/sum*(alpha+1)      );
     748
     749    return l+m<0 ? 0 : TMath::Sqrt((l+m)*2);
     750}
     751
     752Double_t Error(Double_t s, Double_t b, Double_t alpha=0.2)
     753{
     754    const Double_t sN = s + alpha*alpha*b;
     755
     756    const Double_t error = sN<0 ? 0 : TMath::Sqrt(sN);
     757    if (error==0)
     758        return 0;
     759
     760    const Double_t Ns = s - alpha*b;
     761
     762    return Ns<0 ? 0 : Ns/error;
     763}
     764
     765
     766double ganymed2(double xi0, double slope0)
     767{
     768    // Create chain for the tree Result
     769    // This is just easier than using TFile/TTree
     770    TChain c("Result");
     771
     772    // Add the input file to the
     773    c.AddFile("crab-data-only.root");
     774
     775    // Define variables for all leaves to be accessed
     776    // By definition rootifysql writes only doubles
     777    double FileId, EvtNumber, X, Y, MeanX, MeanY, Width, Length, CosDelta, SinDelta,
     778        M3Long, SlopeLong, Leakage1, SlopeSpreadWeighted, Size,
     779        ConcCore, ConcCOG, NumIslands, NumUsedPixels;
     780
     781    // Connect the variables to the cordesponding leaves
     782    c.SetBranchAddress("FileId", &FileId);
     783    c.SetBranchAddress("EvtNumber", &EvtNumber);
     784    c.SetBranchAddress("X", &X);
     785    c.SetBranchAddress("Y", &Y);
     786    c.SetBranchAddress("MeanX", &MeanX);
     787    c.SetBranchAddress("MeanY", &MeanY);
     788    c.SetBranchAddress("Width", &Width);
     789    c.SetBranchAddress("Length", &Length);
     790    c.SetBranchAddress("CosDelta", &CosDelta);
     791    c.SetBranchAddress("SinDelta", &SinDelta);
     792    c.SetBranchAddress("M3Long", &M3Long);
     793    c.SetBranchAddress("SlopeLong", &SlopeLong);
     794    c.SetBranchAddress("Leakage1", &Leakage1);
     795    c.SetBranchAddress("NumIslands", &NumIslands);
     796    c.SetBranchAddress("NumUsedPixels", &NumUsedPixels);
     797    //c.SetBranchAddress("SlopeSpreadWeighted", &SlopeSpreadWeighted);
     798    c.SetBranchAddress("Size", &Size);
     799    //c.SetBranchAddress("ConcCOG", &ConcCOG);
     800    //c.SetBranchAddress("ConcCore", &ConcCore);
     801
     802    // Set some constants (they could be included in the database
     803    // in the future)
     804    double mm2deg = +0.0117193246260285378;
     805    //double abberation = 1.02;
     806
     807    // -------------------- Source dependent parameter calculation -------------------
     808
     809    // Create a histogram for on- and off-data
     810    TH1F hon("on",   "", 55, 0, 1);
     811    TH1F hoff("off", "", 55, 0, 1);
     812
     813    long cnt_on  = 0;
     814    long cnt_off = 0;
     815
     816    for (int i=0; i<c.GetEntries(); i++)
     817    {
     818        // read the i-th event from the file
     819        c.GetEntry(i);
     820
     821        // First calculate all cuts to speedup the analysis
     822        double area = TMath::Pi()*Width*Length;
     823
     824        bool cutq = NumIslands<3.5 && NumUsedPixels>5.5 && Leakage1==0;
     825        if (!cutq)
     826            continue;
     827
     828        bool cut0 = area < log10(Size)*898-1535;
     829        if (!cut0)
     830            continue;
     831
     832        for (int angle=0; angle<360; angle+=60)
     833        {
     834            // -------------------- Source dependent parameter calculation -------------------
     835
     836            double cr = cos(angle*TMath::DegToRad());
     837            double sr = sin(angle*TMath::DegToRad());
     838
     839            double px = cr*X-sr*Y;
     840            double py = cr*Y+sr*X;
     841
     842            //                   blau             rot
     843            double dx = MeanX - px*1.02;
     844            double dy = MeanY - py*1.02;
     845
     846            double norm = sqrt(dx*dx + dy*dy);
     847            double dist = norm*mm2deg;
     848
     849            double lx = min(max((CosDelta*dy - SinDelta*dx)/norm, -1.), 1.);
     850            double ly = min(max((CosDelta*dx + SinDelta*dy)/norm, -1.), 1.);
     851
     852            double alpha = asin(lx);
     853            double sgn   = TMath::Sign(1., ly);
     854
     855            // ------------------------------- Application ----------------------------------
     856
     857            double m3l   = M3Long*sgn*mm2deg;
     858            double slope = SlopeLong*sgn/mm2deg;
     859
     860            // --------------------------------- Analysis -----------------------------------
     861
     862            double xi = xi0 + slope0  *slope;// + 1.67972 *(1-1/(1+4.86232*Leakage1));
     863
     864            double sign1 = m3l+0.07;
     865            double sign2 = (dist-0.5)*7.2-slope;
     866
     867            double disp  = (sign1<0 || sign2<0 ? -xi : xi)*(1-Width/Length);
     868
     869            //double thetasq = (disp*disp + dist*dist - 2*disp*dist*cos(alpha))*mm2deg*mm2deg;
     870            double thetasq = disp*disp + dist*dist - 2*disp*dist*sqrt(1-lx*lx);
     871
     872            // Fill the on- and off-histograms
     873            if (thetasq<0.024)
     874                angle==0 ? cnt_on++ : cnt_off++;
     875        }
     876    }
     877
     878    cout << cnt_on << "\t" << cnt_off << "\t" << cnt_on-cnt_off*0.2 << endl;
     879    return Error(cnt_on, cnt_off);
     880}
     881
     882
     883void ganymed_optimdisp()
     884{
     885    TH2F hist("", "",
     886              21, 1.200-0.0050, 1.410-0.0050,
     887              13, 0.010-0.0025, 0.140-0.0025);
     888
     889    for (float xi0=1.20; xi0<1.41; xi0+=0.01)
     890        for (float slope0=0.01; slope0<0.14; slope0+=0.01)
     891        {
     892            cout << xi0 << "\t" << slope0 << "\t";
     893            hist.Fill(xi0, slope0, ganymed2(xi0, slope0));
     894        }
     895
     896    TCanvas *canv = new TCanvas;
     897    canv->SetTopMargin(0.01);
     898
     899    hist.SetStats(kFALSE);
     900    hist.SetContour(100);
     901    hist.SetXTitle("Xi_{0} / deg");
     902    hist.SetYTitle("Slope_{0} #dot ns/deg^{2} ");
     903    hist.GetYaxis()->SetTitleOffset(1.2);
     904    hist.DrawCopy("colz cont2");
     905}
     906}}}
     907}}}
     908
     909Which gives the following result
     910
     911[[Image(Result3.png)]]
     912
     913
    722914
    723915== Energy Optimization ==