Version 3 (modified by 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!)
#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"); }
Attachments (4)
- Result.png (52.2 KB ) - added by 6 years ago.
- Result2.png (79.3 KB ) - added by 6 years ago.
- Result4.png (63.3 KB ) - added by 6 years ago.
- ThetaSquare.png (15.9 KB ) - added by 6 years ago.
Download all attachments as: .zip
Note:
See TracWiki
for help on using the wiki.