wiki:DatabaseBasedAnalysis/Examples

Version 2 (modified by tbretz, 6 years ago) ( diff )

--

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.

SELECT
   Events.MeanX,
   Events.MeanY,
   Events.Width,
   Events.Length,
   Events.CosDelta,
   Events.SinDelta,
   Events.M3Long,
   Events.SlopeLong,
   Events.Leakage1,
   Events.NumIslands,
   Events.NumUsedPixels,
   Events.Size,
   Position.X,
   Position.Y
FROM
   Events
LEFT JOIN Position USING (FileId, EvtNumber)
LEFT JOIN RunInfo  USING (FileId)
WHERE
   fSourceKey=5
AND
   fRunTypeKey=1
AND
   FileId BETWEEN 170800000 AND 180800000
AND
   fZenithDistanceMax<35
AND
   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 double disp_new = (sign1<0
sign2<0 ? -xi_old : xi_old)*(1-Width/Length);
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");

} }}}

Attachments (4)

Download all attachments as: .zip

Note: See TracWiki for help on using the wiki.