source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/reconstruction/recoviewer/src/RecoTreePainter.cc @ 117

Last change on this file since 117 was 117, checked in by moretto, 11 years ago

ESAF version compilable on mac OS
File size: 61.7 KB
Line 
1// $Id: RecoTreePainter.cc 2927 2011-06-16 18:02:41Z mabl $
2// Author: Dmitry.Naumov   2006/02/03
3
4/*****************************************************************************
5 * ESAF: Euso Simulation and Analysis Framework                              *
6 *                                                                           *
7 *  Id: RecoTreePainter                                                      *
8 *  Package: RecoViewer                                                      *
9 *  Coordinator: Dmitry.Naumov                                               *
10 *                                                                           *
11 *****************************************************************************/
12
13//_____________________________________________________________________________
14//
15// RecoTreePainter
16//
17// <extensive class description>
18//
19//   Config file parameters
20//   ======================
21//
22//   <parameter name>: <parameter description>
23//   -Valid options: <available options>
24//
25
26#include "RecoTreePainter.hh"
27#include "RecoRootEvent.hh"
28#include "EGViewer.hh"
29#include "ESystemOfUnits.hh"
30
31#include <TEventList.h>
32#include <TH1F.h>
33#include <TH1D.h>
34#include <TH2F.h>
35#include <THStack.h>
36#include <TProfile.h>
37#include <TProfile2D.h>
38#include <TList.h>
39#include <TSystem.h>
40#include <TStyle.h>
41#include <TCanvas.h>
42#include <TTree.h>
43#include <TMath.h>
44#include <TF1.h>
45#include <TLegend.h>
46#include <TLegend.h>
47#include <TGraphErrors.h>
48#include <TText.h>
49#include <TPaveText.h>
50#include <TDirectory.h>
51
52using namespace TMath;
53
54ClassImp(RecoTreePainter)
55
56using namespace sou;
57//_____________________________________________________________________________
58RecoTreePainter::RecoTreePainter(TTree* tree): EsafMsgSource(), fEvent(NULL), fSaveEpsFiles(kFALSE) {
59    //
60    //
61    //
62    if (!tree) Error("RecoTreePainter","ctor must be initialized with valid TTree pointer to the recotree\n");
63    fHistoList = new TList;
64    fEvent = new RecoRootEvent;
65    tree->SetBranchAddress("events",&fEvent);
66    fNevents = tree->GetEntries();
67    fTree = tree;
68    Info("RecoTreePainter", "Initialized TTree with %d entries",fNevents);
69
70    fEnergyIndependentFromDir = kTRUE;
71    fUseAA1 = kFALSE;
72
73    //
74    left_edge["Energy"] = 13;
75    left_edge["Theta"]  = 0;
76    left_edge["Phi"]    = 0;
77
78    right_edge["Energy"] = 15;
79    right_edge["Theta"]  = 102;
80    right_edge["Phi"]    = 360;
81    fFilename = "EnergyPlots";
82    gStyle->SetCanvasColor(0);
83    fEpsDir = "./";
84
85    fEventList = new TEventList("recoEventList");
86}
87
88//_____________________________________________________________________________
89RecoTreePainter::~RecoTreePainter() {
90    //
91    // Destructor
92    //
93
94    fHistoList->Delete();
95    delete fHistoList;
96}
97
98//______________________________________________________________________________
99Long64_t RecoTreePainter::UpdateEventList() {
100    //
101    // Update the eventlist according the current cut
102    //
103
104    return fTree->Draw(">>recoEventList",fCut);
105 
106}
107//______________________________________________________________________________
108Long64_t RecoTreePainter::UpdateEventList( const char* selection ) {
109    //
110    // Update the eventlist according selection
111    //
112
113    fCut = selection;
114    return UpdateEventList();
115   
116}
117
118
119//______________________________________________________________________________
120Long64_t RecoTreePainter::GetEntry( Int_t i ) {
121
122
123    if ( i >= fEventList->GetN() ) {
124        MsgForm(EsafMsg::Info,"Entry id %d is greater than the maximum number of selected events (%d)",
125             i,fEventList->GetN());
126        return kFALSE;
127    }
128
129    Long64_t entry = fEventList->GetEntry(i);
130    return fTree->GetEntry(entry);
131}
132
133//_____________________________________________________________________________
134void RecoTreePainter::Loop() {
135    //
136    //
137    //
138    Int_t   progress    = 0;
139    UpdateEventList();
140
141    if ( !fEventList->GetN() ) return;
142
143    for (Int_t i(0); i<fEventList->GetN(); i++) {
144         GetEntry(i);
145         FillHistos();
146
147         if (100.*i/fEventList->GetN() >= progress) {
148         Msg(EsafMsg::Info).SetProgress(progress);
149         Msg(EsafMsg::Info) << "Loop():" << MsgCount;
150         progress+=10;
151          }
152    }
153   
154    FillAngularResolutions();
155   
156    Msg(EsafMsg::Info) << "Loop():" << MsgCount;
157    if (progress <= 100) Msg(EsafMsg::Info) << MsgDispatch;
158}
159
160//_____________________________________________________________________________
161void RecoTreePainter::FillHistos() {
162    //
163    //
164    //
165    MakeBaseDataMap();
166    FillSimHistos();
167    FillThetaHistos();
168    FillPhiHistos();
169    FillDirectionHistos();
170    FillEnergyHistos();
171}
172
173//______________________________________________________________________________
174void RecoTreePainter::FillSimHistos() {
175    //
176    //
177    //
178    map <string, Bool_t>::iterator it;
179    // BULK of simulated information
180    TH1 *hThetaSim = GetHisto("ThetaSimulated_1D", "Theta Simulated", "1D", 100, left_edge["Theta"], right_edge["Theta"],0,0,0);
181    hThetaSim->SetXTitle("#theta_{sim} [deg]");
182    hThetaSim->SetYTitle("Entries");
183    hThetaSim->Fill(fSim["ThetaSim"]*RadToDeg());
184    for (it=fCuts.begin(); it != fCuts.end(); it++) {
185         TString name  = "ThetaSimulated"  + it->first;
186         TString title = "Theta Simulated in " + it->first;
187         TH1* h = GetHisto(name+"_1D",title,"1D",100,left_edge["Theta"],right_edge["Theta"],0,0,0);
188     if (it->second) h->Fill(fSim["ThetaSim"]*RadToDeg());
189
190    }
191    TH1 *hPhiSim = GetHisto("PhiSimulated_1D", "Phi Simulated", "1D", 100, left_edge["Phi"], right_edge["Phi"],0,0,0);
192    hPhiSim->SetXTitle("#varphi_{sim} [deg]");
193    hPhiSim->SetYTitle("Entries");
194    hPhiSim->Fill(fSim["PhiSim"]*RadToDeg());
195    for (it=fCuts.begin(); it != fCuts.end(); it++) {
196         TString name  = "PhiSimulated"  + it->first;
197         TString title = "Phi Simulated in " + it->first;
198         TH1* h = GetHisto(name+"_1D",title,"1D",100,left_edge["Phi"],right_edge["Phi"],0,0,0);
199     if (it->second) h->Fill(fSim["PhiSim"]*RadToDeg());
200
201    } 
202    TH1* h3 = GetHisto("EnergySimulated_1D","Energy Simulated","1D",100,left_edge["Energy"],right_edge["Energy"],0,0,0);
203    h3->SetXTitle("Simulated Energy/EeV");
204    h3->SetYTitle("Entries");
205    h3->Fill(fSim["EnergySim"]);
206    for (it=fCuts.begin(); it != fCuts.end(); it++) {
207         TString name  = "EnergySimulated"  + it->first;
208         TString title = "EnergySimulated in " + it->first;
209         TH1* h = GetHisto(name+"_1D",title,"1D",100,left_edge["Energy"],right_edge["Energy"],0,0,0);
210     if (it->second) h->Fill(fSim["EnergySim"]);
211
212    }
213}
214
215//______________________________________________________________________________
216void RecoTreePainter::FillThetaHistos() {
217    //
218    //
219    //
220    map <string, Bool_t>::iterator it;
221   
222    if (!AcceptDirection()) return;
223
224    Int_t nBinsTheta = TMath::FloorNint(right_edge["Theta"]-left_edge["Theta"])-1;
225    Int_t nBinsPhi = TMath::FloorNint((right_edge["Phi"]-left_edge["Phi"])/2)-1;
226    Double_t errmax = 30;
227    Int_t errbins = 2*TMath::FloorNint(errmax)-1;
228
229    // Theta reco vs theta sim
230    TH1 *h1 = GetHisto("ThetaSim_vs_Rec_2D","Theta Sim vs Rec","2D",nBinsTheta,left_edge["Theta"],right_edge["Theta"],
231                                                                   nBinsTheta,left_edge["Theta"],right_edge["Theta"]);
232    h1->SetXTitle("#theta_{sim} [deg]");
233    h1->SetYTitle("#theta_{rec} [deg]");
234    h1->Fill(fSim["ThetaSim"]*RadToDeg(),fRec["ThetaRec"]*RadToDeg());
235    // theta error
236    h1 = GetHisto("ThetaError","ThetaError","1D",errbins,-errmax,errmax,0,0,0);
237    h1->SetXTitle("#Delta #theta [deg]");
238    h1->SetYTitle("Entries");
239    h1->Fill(fRec["ThetaError"]*RadToDeg());
240    // abs theta error
241    h1 = GetHisto("AbsThetaError","AbsThetaError","1D",180,0,errmax,0,0,0);
242    h1->SetXTitle("|#Delta #theta| [deg]");
243    h1->SetYTitle("Entries");
244    h1->Fill(Abs(fRec["ThetaError"]*RadToDeg()));
245   
246     // Fill Theta Error of vs Theta Simulated
247    h1 = GetHisto("ThetaError_vs_Theta_2D","ThetaError vs Theta","2D",nBinsTheta,left_edge["Theta"],right_edge["Theta"],
248                                            errbins,-errmax,errmax);
249    h1->SetXTitle("#theta_{sim} [deg]");
250    h1->SetYTitle("#Delta #theta [deg]");
251    h1->Fill(fSim["ThetaSim"]*RadToDeg(),fRec["ThetaError"]*RadToDeg());
252    // Fill Theta Error of vs Phi Simulated
253    h1 = GetHisto("ThetaError_vs_Phi_2D","ThetaError vs Phi","2D",nBinsPhi,left_edge["Phi"],right_edge["Phi"],
254                                            errbins,-errmax,errmax);
255    h1->SetXTitle("#varphi_{sim} [deg]");
256    h1->SetYTitle("#Delta #theta [deg]");
257    h1->Fill(fSim["PhiSim"]*RadToDeg(),fRec["ThetaError"]*RadToDeg());
258     // Fill Theta Error of vs signal/noise ratio
259    h1 = GetHisto("ThetaError_vs_SignalToNoise_2D","ThetaError vs S/N","2D",100,0,100,errbins,-errmax,errmax);
260    h1->SetXTitle("S/N");
261    h1->SetYTitle("#Delta #theta [deg]");
262    h1->Fill(fRec["SignalToNoise"],fRec["ThetaError"]*RadToDeg());
263     // Fill Theta Error vs Energy Simulated
264    h1 = GetHisto("ThetaError_vs_EnergySim_2D","ThetaError vs Energy","2D",100,left_edge["Energy"],right_edge["Energy"],errbins,-errmax,errmax);
265    h1->SetXTitle("Simulated Energy/EeV");
266    h1->SetYTitle("#Delta #theta [deg]");
267    h1->Fill(fSim["EnergySim"],fRec["ThetaError"]*RadToDeg());
268     // Fill Theta Error vs Distance in FOV
269    h1 = GetHisto("ThetaError_vs_Distance_2D","ThetaError vs Distance in FOV","2D",100,0,300,errbins,-errmax,errmax);
270    h1->SetXTitle("Distance in FOV [km]");
271    h1->SetYTitle("#Delta #theta [deg]");
272    h1->Fill(fSim["FOVDistance"],fRec["ThetaError"]*RadToDeg());
273
274    // histos with cuts
275    for (it=fCuts.begin(); it != fCuts.end(); it++) {
276        TString name  = "ThetaSim_vs_Rec"  + it->first;
277        TString title = "ThetaSim_vs_Rec " + it->first;
278        TH1* h = NULL;
279
280        // Fill Theta Simulated vs Theta Reconstructed
281        h = GetHisto(name+"_2D",title,"2D",nBinsTheta,left_edge["Theta"],right_edge["Theta"],
282                nBinsTheta,left_edge["Theta"],right_edge["Theta"]);
283        if (it->second) h->Fill(fSim["ThetaSim"]*RadToDeg(),fRec["ThetaRec"]*RadToDeg());
284        // theta error
285        name = "ThetaError" + it->first;
286        title = "ThetaError " + it->first;
287        h = GetHisto(name+"_1D",title,"1D",errbins,-errmax,errmax,0,0,0);
288        if (it->second) h->Fill(fRec["ThetaError"]*RadToDeg());
289        // absolute theta error
290        name = "AbsThetaError" + it->first;
291        title = "AbsThetaError " + it->first;
292        h = GetHisto(name+"_1D",title,"1D",90,0,errmax,0,0,0);
293        if (it->second) h->Fill(Abs(fRec["ThetaError"]*RadToDeg()));
294        // Fill Theta Error of vs Theta Simulated
295        name  = "ThetaError_vs_Theta" + it->first;
296        title = "ThetaError vs Theta"+ it->first;
297        h = GetHisto(name+"_2D",title,"2D",nBinsTheta,left_edge["Theta"],right_edge["Theta"],
298                errbins,-errmax,errmax);
299        if (it->second) h->Fill(fSim["ThetaSim"]*RadToDeg(),fRec["ThetaError"]*RadToDeg());
300        // Fill Theta Error of vs Phi Simulated
301        name  = "ThetaError_vs_Phi" + it->first;
302        title = "ThetaError vs Phi"+ it->first;
303        h = GetHisto(name+"_2D",title,"2D",nBinsTheta,left_edge["Phi"],right_edge["Phi"],
304                errbins,-errmax,errmax);
305        if (it->second) h->Fill(fSim["PhiSim"]*RadToDeg(),fRec["ThetaError"]*RadToDeg());
306        // Fill Theta Error of vs signal/noise ratio
307        name  = "ThetaError_vs_SignalToNoise" + it->first;
308        title = "ThetaError vs S/N"+ it->first;
309        h = GetHisto(name+"_2D",title,"2D",100,0,100,errbins,-errmax,errmax);
310        if (it->second) h->Fill(fRec["SignalToNoise"],fRec["ThetaError"]*RadToDeg());
311        // Fill Theta Error vs Energy Simulated
312        name  = "ThetaError_vs_EnergySim" + it->first;
313        title = "ThetaError vs Energy"+ it->first;
314        h = GetHisto(name+"_2D",title,"2D",100,left_edge["Energy"],right_edge["Energy"],errbins,-errmax,errmax);
315        if (it->second) h->Fill(fSim["EnergySim"],fRec["ThetaError"]*RadToDeg());
316        // Fill Theta Error vs Distance in FOV
317        name = "ThetaError_vs_Distance" + it->first;
318        title = "ThetaError vs Distance in FOV"+ it->first;
319        h = GetHisto(name+"_2D",title,"2D",100,0,300,errbins,-errmax,errmax);
320        if (it->second) h->Fill(fSim["FOVDistance"],fRec["ThetaError"]*RadToDeg());
321    }
322}
323
324//______________________________________________________________________________
325void RecoTreePainter::FillPhiHistos() {
326    //
327    //
328    //
329    map <string, Bool_t>::iterator it;
330   
331    if (!AcceptDirection()) return;
332
333    Int_t nBinsTheta = TMath::FloorNint(right_edge["Theta"]-left_edge["Theta"])-1;
334    Int_t nBinsPhi = TMath::FloorNint((right_edge["Phi"]-left_edge["Phi"])/2)-1;
335    Double_t errmax = 180;
336    Int_t errbins = 2*TMath::FloorNint(errmax/2)-1;
337   
338    TH1* h1 = GetHisto("PhiSim_vs_Rec_2D","Phi Sim vs Rec","2D",nBinsPhi,left_edge["Phi"],right_edge["Phi"],
339                                                                   nBinsPhi,left_edge["Phi"],right_edge["Phi"]);
340    h1->SetXTitle("#varphi_{sim} [deg]");
341    h1->SetYTitle("#varphi_{rec} [deg]");
342    h1->Fill(fSim["PhiSim"]*RadToDeg(),fRec["PhiRec"]*RadToDeg());
343    // phi error
344    h1 = GetHisto("PhiError","PhiError","1D",errbins,-errmax,errmax,0,0,0);
345    h1->SetXTitle("#Delta #varphi [deg]");
346    h1->SetYTitle("Entries");
347    h1->Fill(fRec["PhiError"]*RadToDeg());
348    // abs phi error
349    h1 = GetHisto("AbsPhiError","AbsPhiError","1D",180,0,errmax,0,0,0);
350    h1->SetXTitle("|#Delta #varphi| [deg]");
351    h1->SetYTitle("Entries");
352    h1->Fill(Abs(fRec["PhiError"]*RadToDeg()));
353
354     // Fill phi Error of vs Theta Simulated
355    h1 = GetHisto("PhiError_vs_Theta_2D","PhiError vs Theta","2D",nBinsTheta,left_edge["Theta"],right_edge["Theta"],
356                                            errbins,-errmax,errmax);
357    h1->SetXTitle("#theta_{sim} [deg]");
358    h1->SetYTitle("#Delta #varphi [deg]");
359    h1->Fill(fSim["ThetaSim"]*RadToDeg(),fRec["PhiError"]*RadToDeg());
360    // Fill phi Error of vs Phi Simulated
361    h1 = GetHisto("PhiError_vs_Phi_2D","PhiError vs Phi","2D",nBinsPhi,left_edge["Phi"],right_edge["Phi"],
362                                            errbins,-errmax,errmax);
363    h1->SetXTitle("#varphi_{sim} [deg]");
364    h1->SetYTitle("#Delta #varphi [deg]");
365    h1->Fill(fSim["PhiSim"]*RadToDeg(),fRec["PhiError"]*RadToDeg());
366     // Fill Phi Error of vs signal/noise ratio
367    h1 = GetHisto("PhiError_vs_SignalToNoise_2D","PhiError vs S/N","2D",100,0,100,errbins,-errmax,errmax);
368    h1->SetXTitle("S/N");
369    h1->SetYTitle("#Delta #varphi [deg]");
370    h1->Fill(fRec["SignalToNoise"],fRec["PhiError"]*RadToDeg());
371     // Fill phi Error vs Energy Simulated
372    h1 = GetHisto("PhiError_vs_EnergySim_2D","PhiError vs Energy","2D",100,left_edge["Energy"],right_edge["Energy"],errbins,-errmax,errmax);
373    h1->SetXTitle("Simulated Energy/EeV");
374    h1->SetYTitle("#Delta #varphi [deg]");
375    h1->Fill(fSim["EnergySim"],fRec["PhiError"]*RadToDeg());
376     // Fill phi Error vs Distance in FOV
377    h1 = GetHisto("PhiError_vs_Distance_2D","PhiError vs Distance in FOV","2D",100,0,300,60,-180,180);
378    h1->SetXTitle("Distance in FOV [km]");
379    h1->SetYTitle("#Delta #varphi [deg]");
380    h1->Fill(fSim["FOVDistance"],fRec["PhiError"]*RadToDeg());
381
382    // histos with cuts
383
384    for (it=fCuts.begin(); it != fCuts.end(); it++) {
385         TString name  = "PhiSim_vs_Rec"  + it->first;
386         TString title = "PhiSim_vs_Rec " + it->first;
387         TH1* h = NULL;
388         
389     // Fill Phi Simulated vs Phi Reconstructed
390     h = GetHisto(name+"_2D",title,"2D",nBinsPhi,left_edge["Phi"],right_edge["Phi"],
391                                            nBinsPhi,left_edge["Phi"],right_edge["Phi"]);
392         if (it->second) h->Fill(fSim["PhiSim"]*RadToDeg(),fRec["PhiRec"]*RadToDeg());
393    // phi error
394        name = "PhiError" + it->first;
395        title = "PhiError " + it->first;
396        h = GetHisto(name+"_1D",title,"1D",errbins,-errmax,errmax,0,0,0);
397        if (it->second) h->Fill(fRec["PhiError"]*RadToDeg());
398    // absolute phi error
399        name = "AbsPhiError" + it->first;
400        title = "AbsPhiError" + it->first;
401        h = GetHisto(name+"_1D",title,"1D",90,0,errmax,0,0,0);
402        if (it->second) h->Fill(Abs(fRec["PhiError"]*RadToDeg()));
403     // Fill Phi Error of vs Theta Simulated
404     name  = "PhiError_vs_Theta" + it->first;
405     title = "PhiError vs Theta"+ it->first;
406     h = GetHisto(name+"_2D",title,"2D",nBinsTheta,left_edge["Theta"],right_edge["Theta"],
407                                            errbins,-errmax,errmax);
408         if (it->second) h->Fill(fSim["ThetaSim"]*RadToDeg(),fRec["PhiError"]*RadToDeg());
409     // Fill Phi Error of vs Phi Simulated
410     name  = "PhiError_vs_Phi" + it->first;
411     title = "PhiError vs Phi"+ it->first;
412     h = GetHisto(name+"_2D",title,"2D",nBinsPhi,left_edge["Phi"],right_edge["Phi"],
413                                            errbins,-errmax,errmax);
414         if (it->second) h->Fill(fSim["PhiSim"]*RadToDeg(),fRec["PhiError"]*RadToDeg());
415     // Fill Phi Error of vs signal/noise ratio
416     name  = "PhiError_vs_SignalToNoise" + it->first;
417     title = "PhiError vs S/N"+ it->first;
418     h = GetHisto(name+"_2D",title,"2D",100,0,100,errbins,-errmax,errmax);
419         if (it->second) h->Fill(fRec["SignalToNoise"],fRec["PhiError"]*RadToDeg());
420     // Fill Phi Error vs Energy Simulated
421     name  = "PhiError_vs_EnergySim" + it->first;
422     title = "PhiError vs Energy"+ it->first;
423     h = GetHisto(name+"_2D",title,"2D",100,left_edge["Energy"],right_edge["Energy"],errbins,-errmax,errmax);
424         if (it->second) h->Fill(fSim["EnergySim"],fRec["PhiError"]*RadToDeg());
425     // Fill Phi Error vs Distance in FOV
426     name = "PhiError_vs_Distance" + it->first;
427     title = "PhiError vs Distance in FOV"+ it->first;
428     h = GetHisto(name+"_2D",title,"2D",100,0,300,errbins,-errmax,errmax);
429         if (it->second) h->Fill(fSim["FOVDistance"],fRec["PhiError"]*RadToDeg());
430
431   }
432}
433
434//______________________________________________________________________________
435void RecoTreePainter::FillAngularResolutions() {
436    //
437    //
438    //
439   
440    map <string, Bool_t>::iterator it;
441   
442    // search theta68
443    TH1 *h1 = GetHisto("AbsThetaError","","1D",0,0,0,0,0,0);
444    Double_t integral = h1->Integral();
445    Double_t theta68 = -1.;
446    Int_t bin68 = -1;
447    for( Int_t i(0); i<h1->GetNbinsX(); i++ ) {
448        if ( (h1->Integral(1,i+1))/integral >= 0.68 ) {
449            bin68 = i+1;
450            theta68 = h1->GetBinCenter(i+1);
451            break;
452        }
453    }
454    Double_t th68max = h1->GetBinLowEdge(bin68)+h1->GetBinWidth(bin68);
455    TH1 *h68 = GetHisto("Theta68","Theta68","1D",bin68,0,th68max,0,0,0);
456    for( Int_t i(0); i<bin68; i++ ) h68->SetBinContent(i+1,h1->GetBinContent(i+1));
457
458    for (it=fCuts.begin(); it != fCuts.end(); it++) {
459        TString name = "AbsThetaError" + it->first;
460        TString title = "AbsThetaError " + it->first;
461        TH1 *h = GetHisto(name+"_1D",title,"1D",0,0,0,0,0,0);
462    // search theta68
463    integral = h->Integral();
464    theta68 = -1.;
465    bin68 = -1;
466    for( Int_t i(0); i<h->GetNbinsX(); i++ ) {
467        if ( (h->Integral(1,i+1))/integral >= 0.68 ) {
468            bin68 = i+1;
469            theta68 = h->GetBinCenter(i+1);
470            break;
471        }
472    }
473    th68max = h->GetBinLowEdge(bin68)+h->GetBinWidth(bin68);
474    name = "Theta68" + it->first;
475    title = "Theta68 " + it->first;
476    TH1 *h68 = GetHisto(name+"_1D",title,"1D",bin68,0,th68max,0,0,0);
477    for( Int_t i(0); i<bin68; i++ ) h68->SetBinContent(i+1,h->GetBinContent(i+1));
478    }
479   
480    // search phi68
481    h1 = GetHisto("AbsPhiError","AbsPhiError","1D",0,0,0,0,0,0);
482    integral = h1->Integral();
483    Double_t phi68 = -1.;
484    bin68 = -1;
485    for( Int_t i(0); i<h1->GetNbinsX(); i++ ) {
486        if ( (h1->Integral(1,i+1))/integral >= 0.68 ) {
487            bin68 = i+1;
488            phi68 = h1->GetBinCenter(i+1);
489            break;
490        }
491    }
492    Double_t phi68max = h1->GetBinLowEdge(bin68)+h1->GetBinWidth(bin68);
493    h68 = GetHisto("Phi68","Phi68","1D",bin68,0,phi68max,0,0,0);
494    for( Int_t i(0); i<bin68; i++ ) h68->SetBinContent(i+1,h1->GetBinContent(i+1));
495
496    for (it=fCuts.begin(); it != fCuts.end(); it++) {
497        TString name = "AbsPhiError" + it->first;
498        TString title = "AbsPhiError" + it->first;
499        TH1 *h = GetHisto(name+"_1D",title,"1D",0,0,0,0,0,0);
500    // search phi68
501    integral = h->Integral();
502    phi68 = -1.;
503    bin68 = -1;
504    for( Int_t i(0); i<h->GetNbinsX(); i++ ) {
505        if ( (h->Integral(1,i+1))/integral >= 0.68 ) {
506            bin68 = i+1;
507            phi68 = h->GetBinCenter(i+1);
508            break;
509        }
510    }
511    phi68max = h->GetBinLowEdge(bin68)+h->GetBinWidth(bin68);
512    name = "Phi68" + it->first;
513    title = "Phi68 " + it->first;
514    TH1 *h68 = GetHisto(name+"_1D",title,"1D",bin68,0,phi68max,0,0,0);
515    for( Int_t i(0); i<bin68; i++ ) h68->SetBinContent(i+1,h->GetBinContent(i+1));
516    }
517
518    // calculate the resolution in the direction
519    for (it=fCuts.begin(); it != fCuts.end(); it++) {
520        TString name = "DirError" + it->first;
521        TString title = "DirError" + it->first;
522        TH1 *h = GetHisto(name+"_1D",title,"1D",0,0,0,0,0,0);
523    // search psi68
524    integral = h->Integral();
525    phi68 = -1.;
526    bin68 = -1;
527    for( Int_t i(0); i<h->GetNbinsX(); i++ ) {
528        if ( (h->Integral(1,i+1))/integral >= 0.68 ) {
529            bin68 = i+1;
530            phi68 = h->GetBinCenter(i+1);
531            break;
532        }
533    }
534    phi68max = h->GetBinLowEdge(bin68)+h->GetBinWidth(bin68);
535    name = "Psi68" + it->first;
536    title = "Psi68 " + it->first;
537    TH1 *h68 = GetHisto(name+"_1D",title,"1D",bin68,0,phi68max,0,0,0);
538    for( Int_t i(0); i<bin68; i++ ) h68->SetBinContent(i+1,h->GetBinContent(i+1));
539    }
540
541}
542
543//______________________________________________________________________________
544void RecoTreePainter::FillDirectionHistos() {
545    //
546    //
547    //
548    map <string, Bool_t>::iterator it;
549   
550    if (!AcceptDirection()) return;
551
552    for (it=fCuts.begin(); it != fCuts.end(); it++) {
553         TString name  = "DirError_vs_Phi"  + it->first;
554         TString title = "DirError_vs_Phi " + it->first;
555         TH1* h = NULL;
556     // Fill Direction Error vs Phi Simulated
557     h = GetHisto(name+"_2D",title,"2D",100,left_edge["Phi"],right_edge["Phi"],
558                                            100,0,180);
559         if (it->second) h->Fill(fSim["PhiSim"]*RadToDeg(),fRec["DirectionError"]*RadToDeg());
560     // Fill Dir Error
561     name  = "DirError" + it->first;
562     title = "DirError "+ it->first;
563     h = GetHisto(name+"_1D",title,"1D",360,0,180,0,0,0);
564     if (it->second) h->Fill(fRec["DirectionError"]*RadToDeg());
565     // Fill Dir Error of vs Theta Simulated
566     name  = "DirError_vs_Theta" + it->first;
567     title = "DirError vs Theta"+ it->first;
568     h = GetHisto(name+"_2D",title,"2D",100,left_edge["Theta"],right_edge["Theta"],
569                                            100,0,180);
570         if (it->second) h->Fill(fSim["ThetaSim"]*RadToDeg(),fRec["DirectionError"]*RadToDeg());
571     // Fill Theta Error of vs signal/noise ratio
572     name  = "DirError_vs_SignalToNoise" + it->first;
573     title = "DirError vs S/N"+ it->first;
574     h = GetHisto(name+"_2D",title,"2D",100,0,100,100,0,180);
575         if (it->second) h->Fill(fRec["SignalToNoise"],fRec["DirectionError"]*RadToDeg());
576     // Fill Theta Error vs Energy Simulated
577     name  = "DirError_vs_EnergySim" + it->first;
578     title = "DirError vs Energy"+ it->first;
579     h = GetHisto(name+"_2D",title,"2D",100,left_edge["Energy"],right_edge["Energy"],100,0,180);
580         if (it->second) h->Fill(fSim["EnergySim"],fRec["DirectionError"]*RadToDeg());
581   }
582
583}
584
585//______________________________________________________________________________
586void RecoTreePainter::FillEnergyHistos() {
587    //
588    //
589    //
590   
591    map <string, Bool_t>::iterator it;
592   
593    if (!AcceptEnergy()) return;
594
595    TH1* h1 = GetHisto("EnergySim_vs_Rec_2D","Energy Sim vs Rec","2D",100,left_edge["Energy"],right_edge["Energy"],
596                                                                   100,left_edge["Energy"],right_edge["Energy"]);
597    h1->SetXTitle("Simulated Energy/EeV");
598    h1->SetYTitle("Reconstructed Energy/EeV");
599    h1->Fill(fSim["EnergySim"],fRec["EnergyRec"]);
600
601    for (it=fCuts.begin(); it != fCuts.end(); it++) {
602        TString name  = "EnergySim_vs_Rec"  + it->first;
603        TString title = "EnergySim_vs_Rec " + it->first;
604        TH1* h = NULL;
605       
606        // Fill Energy Simulated vs Energy Reconstructed
607        h = GetHisto(name+"_2D",title,"2D",100,left_edge["Energy"],right_edge["Energy"],
608                100,left_edge["Energy"],right_edge["Energy"]);
609        if (it->second) h->Fill(fSim["EnergySim"],fRec["EnergyRec"]);
610       
611        // Fill Relative Error of vs Energy Simulated
612        name  = "EnergyErrorRelative_vs_Energy" + it->first;
613        title = "EnergyError Relative vs Energy"+ it->first;
614        h = GetHisto(name+"_2D",title,"2D",100,left_edge["Energy"],right_edge["Energy"],
615                100,-2,2);
616        if (it->second) h->Fill(fSim["EnergySim"],fRec["EnergyErrorRelative"]);
617       
618        // Fill Relative Error of vs signal/noise ratio
619        name  = "EnergyErrorRelative_vs_SignalToNoise" + it->first;
620        title = "EnergyError Relative vs S/N"+ it->first;
621        h = GetHisto(name+"_2D",title,"2D",100,0,100,100,-2,2);
622        if (it->second) h->Fill(fRec["SignalToNoise"],fRec["EnergyErrorRelative"]);
623       
624        // Fill signal/noise ratio vs Energy Simulated
625        name  = "SignalToNoise_vs_EnergySim" + it->first;
626        title = "S/N vs Energy"+ it->first;
627        h = GetHisto(name+"_2D",title,"2D",100,left_edge["Energy"],right_edge["Energy"],100,0,100);
628        if (it->second) h->Fill(fSim["EnergySim"],fRec["SignalToNoise"]);
629       
630        // Fill simulated Theta angle for this theta bin
631        name  = "SimulatedTheta"+it->first;
632        title = "Simulated Theta in "+it->first;
633        h = GetHisto(name,title,"1D",100,left_edge["Theta"],right_edge["Theta"],0,0,0);
634        if (it->second) h->Fill(fSim["ThetaSim"]*RadToDeg());
635       
636        // Fill simulated Phi angle for this theta bin
637        name  = "SimulatedPhi"+it->first;
638        title = "Simulated Phi in "+it->first;
639        h = GetHisto(name,title,"1D",100,left_edge["Phi"],right_edge["Phi"],0,0,0);
640        if (it->second) h->Fill(fSim["PhiSim"]*RadToDeg());
641       
642        // Fill reconstructed Energy  for this theta bin
643        name  = "ReconstructedEnergy:"+it->first;
644        title = "Reconstructed Energy in "+it->first;
645        h = GetHisto(name,title,"1D",100,left_edge["Energy"],right_edge["Energy"],0,0,0);
646        if (it->second) h->Fill(fRec["EnergyRec"]);
647       
648        // Fill simulated Energy  for this theta bin
649        name  = "AcceptedEventsRatio:"+it->first;
650        title = "Fraction of accepted events in "+it->first;
651        h = GetHisto(name,title,"1D",100,left_edge["Energy"],right_edge["Energy"],0,0,0);
652        //   if (it->second && TMath::Abs(fRec["EnergyErrorRelative"])<0.3) h->Fill(fSim["EnergySim"]);
653        if (it->second ) h->Fill(fSim["EnergySim"]);
654       
655        // Fill statistical expectation for the resolution
656        name  = "StatisticalExpectation"+it->first;
657        title = "Statistical Expectation in "+it->first;
658        h = GetHisto(name,title,"1D",100,0,100,0,0,0);
659        if (it->second) h->Fill(fRec["SignalToNoise"],1./TMath::Sqrt(fRec["Signal"]));
660    }
661
662
663    TProfile2D* h2 = (TProfile2D*)GetHisto("EnergyError_prof2","Relative Energy Error","prof2",10,left_edge["Energy"],right_edge["Energy"],
664                                                         10,left_edge["Theta"],right_edge["Theta"]);
665    h2->SetXTitle("Simulated Energy/EeV");
666    h2->SetYTitle("Simulated Theta/deg");
667    h2->Fill(fSim["EnergySim"],fSim["ThetaSim"]*TMath::RadToDeg(),fRec["EnergyErrorRelative"]);
668
669    TH2F* h21 = (TH2F*)GetHisto("EnergyRatio_2D","Reconstructed/Simulated","2D",100,left_edge["Energy"],right_edge["Energy"],100,0,3);
670    h21->SetXTitle("Simulated Energy/EeV");
671    h21->SetYTitle("Reconstructed/Simulated Energy");
672    h21->Fill(fSim["EnergySim"],fRec["EnergyRec"]/fSim["EnergySim"]);
673
674    TH1 *h4 = GetHisto("Altitude_vs_Theta","Altitude vs Theta","2D",100,left_edge["Theta"],right_edge["Theta"],100,0,25);
675    h4->Fill(fSim["ThetaSim"]*RadToDeg(),fRec["Altitude"]);
676
677    TH1 *h5 = GetHisto("Quality_vs_Energy","Quality vs Energy","2D",100,left_edge["Energy"],right_edge["Energy"],100,0,25);
678    h5->Fill(fSim["EnergySim"],fRec["EnergyQuality"]);
679}
680//_____________________________________________________________________________
681Bool_t RecoTreePainter::AcceptDirection() {
682     //
683     //
684     //
685     if ( fEvent->GetRecoTrackDirection2().GetQuality() < 0 )                                      return kFALSE;
686     //if ( Sqrt(Power(fEvent->GetTruth().GetMaxPos().X(),2)+Power(fEvent->GetTruth().GetMaxPos().X(),2)) > 200 ) return kFALSE;
687     //if ( fEvent->GetTruth().GetEnergy()<1.e14 ) return kFALSE;
688     return kTRUE;
689}
690
691//_____________________________________________________________________________
692Bool_t RecoTreePainter::AcceptEnergy() {
693    //
694    //
695    //
696     if ( !(AcceptDirection() || fEnergyIndependentFromDir ) ) return kFALSE;
697     if (fEvent->GetRecoEnergy().GetEnergy() == 5.e14)                                            return kFALSE;
698     if (fEvent->GetRecoEnergy().GetQuality() < 0 )     return kFALSE;
699//     if (fEvent->GetRecoEnergy().GetQuality() < 0 || fEvent->GetRecoEnergy().GetQuality()> 30)     return kFALSE;
700     if (fEvent->GetRecoEnergy().GetEnergy() <= 0 ) return kFALSE;
701     return kTRUE;
702}
703
704//_____________________________________________________________________________
705void RecoTreePainter::MakeBaseDataMap() {
706    //
707    //
708    //
709    fSim["EnergySim"] = (Double_t)Log10(fEvent->GetTruth().GetEnergy());
710    fSim["ThetaSim"]  = (Double_t)fEvent->GetTruth().GetTheta();
711    fSim["PhiSim"]    = (Double_t)fEvent->GetTruth().GetPhi();
712    fSim["FOVDistance"] = (Double_t)Sqrt(Power(fEvent->GetTruth().GetMaxPos().X(),2)+Power(fEvent->GetTruth().GetMaxPos().Y(),2));
713
714    if (!fUseAA1) {
715        fRec["ThetaRec"]         = (Double_t)fEvent->GetRecoTrackDirection2().GetTheta()*DegToRad();
716        fRec["PhiRec"]           = (Double_t)fEvent->GetRecoTrackDirection2().GetPhi()*DegToRad();
717        fRec["ThetaError"]       = (Double_t)fEvent->GetRecoTrackDirection2().GetErrorTheta();
718        fRec["PhiError"]         = (Double_t)fEvent->GetRecoTrackDirection2().GetErrorPhi();
719        fRec["DirectionError"]   = (Double_t)fEvent->GetRecoTrackDirection2().GetErrorDirection()*DegToRad();
720    } else {
721        fRec["ThetaRec"]         = (Double_t)fEvent->GetRecoTrackDirection2().GetThetaAA1()*DegToRad();
722        fRec["PhiRec"]           = (Double_t)fEvent->GetRecoTrackDirection2().GetPhiAA1()*DegToRad();
723        fRec["ThetaError"]       = (Double_t)(fEvent->GetRecoTrackDirection2().GetThetaAA1() - fEvent->GetTruth().GetTheta()*RadToDeg());
724        fRec["PhiError"]         = (Double_t)(fEvent->GetRecoTrackDirection2().GetPhiAA1() - fEvent->GetTruth().GetPhi()*RadToDeg());
725        fRec["DirectionError"]   = (Double_t)fEvent->GetRecoTrackDirection2().GetErrorDirectionAA1()*DegToRad();
726    }
727    if (fRec["ThetaError"]>180) fRec["ThetaError"] -= 360.;
728    if (fRec["ThetaError"]<-180) fRec["ThetaError"] += 360.;
729    fRec["ThetaError"] *= DegToRad();
730    if (fRec["PhiError"]>180) fRec["PhiError"] -= 360.;
731    if (fRec["PhiError"]<-180) fRec["PhiError"] += 360.;
732    fRec["PhiError"] *= DegToRad();
733    fRec["TDPError"]         = (Double_t)fEvent->GetRecoTrackDirection2().GetErrorTDP()*DegToRad();
734
735    fRec["EnergyRec"]           = (Double_t)TMath::Log10(fEvent->GetRecoEnergy().GetEnergy());
736    fRec["EnergyError"]         = (Double_t)fEvent->GetRecoEnergy().GetEnergyMCError();
737    fRec["EnergyErrorRelative"] = (Double_t)fEvent->GetRecoEnergy().GetEnergyMCError()/fEvent->GetTruth().GetEnergy();
738    fRec["EnergyQuality"]       = (Double_t)fEvent->GetRecoEnergy().GetQuality();
739    fRec["SignalToNoise"]       = (Double_t)fEvent->GetRecoEnergy().GetSignalCluster()/TMath::Sqrt(fEvent->GetRecoEnergy().GetNoiseCluster());
740    fRec["Signal"]              = (Double_t)fEvent->GetRecoEnergy().GetSignalCluster();
741    fRec["Altitude"]            = (Double_t)fEvent->GetRecoEnergy().GetAltitude()/km;
742
743    // seting up the cuts
744    fCuts["0.lt.Theta.lt.15"]  = 0*DegToRad()  < fEvent->GetTruth().GetTheta() && fEvent->GetTruth().GetTheta() <= 15*DegToRad();
745    fCuts["15.lt.Theta.lt.30"] = 15*DegToRad() < fEvent->GetTruth().GetTheta() && fEvent->GetTruth().GetTheta() <= 30*DegToRad();
746    fCuts["30.lt.Theta.lt.45"] = 30*DegToRad() < fEvent->GetTruth().GetTheta() && fEvent->GetTruth().GetTheta() <= 45*DegToRad();
747    fCuts["45.lt.Theta.lt.60"] = 45*DegToRad() < fEvent->GetTruth().GetTheta() && fEvent->GetTruth().GetTheta() <= 60*DegToRad();
748    fCuts["60.lt.Theta.lt.75"] = 60*DegToRad() < fEvent->GetTruth().GetTheta() && fEvent->GetTruth().GetTheta() <= 75*DegToRad();
749    fCuts["75.lt.Theta.lt.102"] = 75*DegToRad() < fEvent->GetTruth().GetTheta() && fEvent->GetTruth().GetTheta() <= 102*DegToRad();
750}
751
752//_____________________________________________________________________________
753void RecoTreePainter::DrawAll() {
754    //
755    //
756    //
757    DrawTheta();
758    DrawPhi();
759    //DrawDirection();
760    DrawEnergy();
761}
762
763//_____________________________________________________________________________
764void RecoTreePainter::DrawTheta() {
765    //
766    //
767    //
768    fThetaViewer = new EGViewer();
769    fThetaViewer->SetName("ThetaReconstuction");
770    fThetaViewer->SetTitle("Theta Reconstuction");
771    fThetaViewer->SetWindowName("Theta Reconstuction");
772
773    map <string, Bool_t>::iterator it;
774    TString cname = "Theta ";
775    TCanvas *c1 = fThetaViewer->FindCanvas(cname);
776    if(!c1 )  c1 = fThetaViewer->AddTab(cname,cname);
777    c1->Clear();
778    c1->cd();
779    c1->Divide(3,3);
780    // draw simulated theta distribution: all events and superimposed in the given bin
781    c1->cd(1);
782    TH1* h1 = GetHisto("ThetaSimulated_1D","Theta Simulated","1D",0,0,0,0,0,0);
783    TH1* h2 = GetHisto("ThetaError","","1D",0,0,0,0,0,0);
784    h1->Draw();
785    TPaveText *summary = new TPaveText(0.65,0.75,0.98,0.98,"brNDC");
786    summary->SetBorderSize(1);
787    summary->SetTextAlign(12);
788    summary->AddText("Events");
789    summary->AddText(Form("%d simulated",(Int_t)h1->GetEntries()));
790    summary->AddText(Form("%d reconstructed",(Int_t)h2->GetEntries()));
791    summary->Draw();
792
793
794
795
796    // draw simulated theta vs reconstructed theta
797    c1->cd(2);
798    TH1* h3 = GetHisto("ThetaSim_vs_Rec_2D","","2D",0,0,0,0,0,0);
799    h3->Draw("colz");
800    // draw theta error
801    c1->cd(4);
802    h2->Draw();
803    // draw theta abs error
804    TVirtualPad *p51 = c1->cd(5);
805    p51->SetLogy();
806    TH1* h21 = GetHisto("AbsThetaError","","1D",0,0,0,0,0,0);
807    h21->Draw();
808    TH1* h22 = GetHisto("Theta68","","1D",0,0,0,0,0,0);
809    h22->SetFillColor(kRed);
810    h22->Draw("same");
811
812    // draw Theta Error of vs Theta Simulated
813    c1->cd(3);
814    TH1* h4 = GetHisto("ThetaError_vs_Theta_2D","","2D",0,0,0,0,0,0);
815    h4->Draw("colz");
816    // draw Theta Error of vs Theta Simulated
817    c1->cd(6);
818    TH1* h5 = GetHisto("ThetaError_vs_Phi_2D","","2D",0,0,0,0,0,0);
819    h5->Draw("colz");
820    // draw Theta Error of vs S/N
821    c1->cd(7);
822    TH1* h6 = GetHisto("ThetaError_vs_SignalToNoise_2D","","2D",0,0,0,0,0,0);
823    h6->Draw("colz");
824    // draw Theta Error of vs Energy
825    c1->cd(8);
826    TH1* h7 = GetHisto("ThetaError_vs_EnergySim_2D","","2D",0,0,0,0,0,0);
827    h7->Draw("colz");
828    // draw Theta Error of vs Energy
829    c1->cd(9);
830    TH1* h8 = GetHisto("ThetaError_vs_Distance_2D","","2D",0,0,0,0,0,0);
831    h8->Draw("colz");
832
833    Int_t cuts_size = fCuts.size();
834    Double_t effvstheta[(const Int_t)cuts_size]; 
835    Int_t i(0);
836    for (it=fCuts.begin(); it != fCuts.end(); it++) {
837        TCanvas *c = fThetaViewer->FindCanvas(cname+(it->first).c_str());
838        if(!c )  c = fThetaViewer->AddTab(cname+(it->first).c_str(),cname+(it->first).c_str());
839        c->Clear();
840        c->cd();
841        c->Divide(3,3);
842        TString name, title;
843        // draw simulated theta distribution: all events and superimposed in the given bin
844        c->cd(1);
845        TH1* h1 = GetHisto("ThetaSimulated_1D","Theta Simulated","1D",0,0,0,0,0,0);
846        h1->DrawCopy();
847
848        name = "ThetaSimulated"  + it->first;
849        TH1* hTot = GetHisto(name+"_1D","","1D",0,0,0,0,0,0);
850        TH1* hCopyTot = hTot->DrawCopy("same");
851        hCopyTot->SetFillColor(kBlue);
852        hCopyTot->SetLineColor(kBlue);
853       
854        name  = "SimulatedTheta"+it->first;
855        TH1* h2 = GetHisto(name,"","1D",0,0,0,0,0,0);
856        TH1* hCopyRe = h2->DrawCopy("same");
857        hCopyRe->SetXTitle(" #theta_{sim} [deg]");
858        hCopyRe->SetFillColor(kRed);
859        hCopyRe->SetLineColor(kRed);
860        name = "ThetaError" + it->first;
861        TH1 *hT = GetHisto(name+"_1D","","1D",0,0,0,0,0,0);
862       
863        TLegend *leg = new TLegend(0.65,0.75,0.98,0.98,"Events","brNDC");
864        leg->SetBorderSize(1);
865        leg->SetTextAlign(12);
866        leg->AddEntry(hCopyTot,Form("%d simulated",(Int_t)hCopyTot->GetEntries()),"F");
867        leg->AddEntry(hCopyRe,Form("%d reconstructed",(Int_t)hT->GetEntries()),"F");
868        leg->SetBit(kCanDelete);
869        leg->Draw();
870
871        //reco eff
872        effvstheta[i++] = (Double_t)hT->GetEntries()/hCopyTot->GetEntries(); 
873
874        // draw theta error
875        c->cd(4);
876
877        hT->Draw();
878        hT->SetXTitle("#Delta #theta [deg]");
879        hT->SetYTitle("Entries");
880
881        // draw abs theta error
882        TVirtualPad *p5 = c->cd(5);
883        p5->SetLogy();
884        name = "AbsThetaError" + it->first;
885        TH1 *hTabs = GetHisto(name+"_1D","","1D",0,0,0,0,0,0);
886        hTabs->Draw();
887        hTabs->SetXTitle("|#Delta #theta| [deg]");
888        hTabs->SetYTitle("Entries");
889        name = "Theta68" + it->first;
890        TH1 *h68 = GetHisto(name+"_1D","","1D",0,0,0,0,0,0);
891        h68->SetFillColor(kRed);
892        h68->Draw("same");
893
894        // draw simulated theta vs reconstructed theta
895        c->cd(2);
896        name  = "ThetaSim_vs_Rec"  + it->first;
897        TH1* h3 = GetHisto(name+"_2D","","2D",0,0,0,0,0,0);
898        h3->Draw("colz");
899        h3->SetXTitle("#theta_{sim} [deg]");
900        h3->SetYTitle("#theta_{rec} [deg]");
901        // draw Theta Error of vs Theta Simulated
902        c->cd(3);
903        name  = "ThetaError_vs_Theta" + it->first;
904        TH1* h4 = GetHisto(name+"_2D","","2D",0,0,0,0,0,0);
905        h4->SetXTitle("#theta_{sim} [deg]");
906        h4->SetYTitle("#Delta #theta [deg]");
907        h4->Draw("colz");
908        // draw Theta Error of vs Theta Simulated
909        c->cd(6);
910        name  = "ThetaError_vs_Phi" + it->first;
911        TH1* h5 = GetHisto(name+"_2D","","2D",0,0,0,0,0,0);
912        h5->SetXTitle("#varphi_{sim} [deg]");
913        h5->SetYTitle("#Delta #theta [deg]");
914        h5->Draw("colz");
915        // draw Theta Error of vs S/N
916        c->cd(7);
917        name  = "ThetaError_vs_SignalToNoise" + it->first;
918        TH1* h6 = GetHisto(name+"_2D","","2D",0,0,0,0,0,0);
919        h6->SetXTitle("S/N");
920        h6->SetYTitle("#Delta #theta [deg]");
921        h6->Draw("colz");
922        // draw Theta Error of vs Energy
923        c->cd(8);
924        name  = "ThetaError_vs_EnergySim" + it->first;
925        TH1* h7 = GetHisto(name+"_2D","","2D",0,0,0,0,0,0);
926        h7->SetXTitle("Simulated log_{10}Energy");
927        h7->SetYTitle("#Delta #theta [deg]");
928        h7->Draw("colz");
929        // draw Theta Error of vs Energy
930        c->cd(9);
931        name  = "ThetaError_vs_Distance" + it->first;
932        TH1* h8 = GetHisto(name+"_2D","","2D",0,0,0,0,0,0);
933        h8->SetXTitle("Distance in FOV [km]");
934        h8->SetYTitle("#Delta #theta [deg]");
935        h8->Draw("colz");
936
937    }
938
939    cname = "Summary";
940    TCanvas *c2 = fThetaViewer->FindCanvas(cname);
941    if(!c2 )  c2 = fThetaViewer->AddTab(cname,cname);
942    c2->Clear();
943    c2->cd();
944    c2->Divide(2,2);
945
946    Double_t theta[] = {0.,15.,30.,45.,60.,75.};
947    Double_t theta68[(const Int_t)cuts_size];
948    Double_t psi68[(const Int_t)cuts_size];
949    i = 0;
950    for (it=fCuts.begin(); it != fCuts.end(); it++) {
951        TString hname = "Theta68" + it->first;
952        TH1* histo = GetHisto(hname+"_1D","","1D",0,0,0,0,0,0);
953        Int_t nbins = histo->GetNbinsX();
954        theta68[i] = histo->GetBinCenter(nbins);
955        hname = "Psi68" + it->first;
956        histo = GetHisto(hname+"_1D","","1D",0,0,0,0,0,0);
957        nbins = histo->GetNbinsX();
958        psi68[i] = histo->GetBinCenter(nbins);
959        i++;
960    }
961
962    c2->cd(1);
963    TGraph *dirRes = new TGraph(6,theta,psi68);
964    dirRes->SetName("DirectionResolution");
965    dirRes->SetTitle("Direction Resolution");
966    dirRes->GetXaxis()->SetTitle("#theta [deg]");
967    dirRes->GetYaxis()->SetTitle("#psi_{68} [deg]");
968    dirRes->SetMarkerStyle(22);
969    dirRes->SetMarkerColor(kRed);
970    dirRes->Draw("AP");
971
972    c2->cd(2);
973    TGraph *thRes = new TGraph(6,theta,theta68);
974    thRes->SetName("ThetaResolution");
975    thRes->SetTitle("Theta Resolution");
976    thRes->GetXaxis()->SetTitle("#theta [deg]");
977    thRes->GetYaxis()->SetTitle("#theta_{68} [deg]");
978    thRes->SetMarkerStyle(21);
979    thRes->SetMarkerColor(kBlue);
980    thRes->Draw("AP");
981
982    // reconstruction efficiency vs theta
983    c2->cd(3);
984    TGraph *effTheta = new TGraph(6,theta,effvstheta);
985    effTheta->SetName("RecoEfficiency");
986    effTheta->SetTitle("Reconstruction Efficiency vs #theta");
987    effTheta->GetXaxis()->SetTitle("#theta [deg]");
988    effTheta->GetYaxis()->SetTitle("Eff.");
989    effTheta->SetMarkerStyle(22);
990    effTheta->SetMarkerColor(kBlue);
991    effTheta->Draw("AP");
992
993}
994
995//_____________________________________________________________________________
996void RecoTreePainter::DrawPhi() {
997    //
998    //
999    //
1000    fPhiViewer = new EGViewer();
1001    fPhiViewer->SetName("PhiReconstuction");
1002    fPhiViewer->SetTitle("Phi Reconstuction");
1003    fPhiViewer->SetWindowName("Phi Reconstuction");
1004
1005    map <string, Bool_t>::iterator it;
1006    TString cname = "Phi ";
1007    TCanvas *c1 = fPhiViewer->FindCanvas(cname);
1008    if(!c1 )  c1 = fPhiViewer->AddTab(cname,cname);
1009    c1->Clear();
1010    c1->cd();
1011    c1->Divide(3,3);
1012    // draw simulated phi distribution: all events and superimposed in the given bin
1013    c1->cd(1);
1014    TH1* h1 = GetHisto("PhiSimulated_1D","Phi Simulated","1D",0,0,0,0,0,0);
1015    h1->Draw();
1016    TH1* h2 = GetHisto("PhiError","","1D",0,0,0,0,0,0);
1017    TPaveText *summary = new TPaveText(0.65,0.75,0.98,0.98,"brNDC");
1018    summary->SetBorderSize(1);
1019    summary->SetTextAlign(12);
1020    summary->AddText("Events");
1021    summary->AddText(Form("%d simulated",(Int_t)h1->GetEntries()));
1022    summary->AddText(Form("%d reconstructed",(Int_t)h2->GetEntries()));
1023    summary->Draw();
1024    // draw simulated phi vs reconstructed theta
1025    c1->cd(2);
1026    TH1* h3 = GetHisto("PhiSim_vs_Rec_2D","","2D",0,0,0,0,0,0);
1027    h3->Draw("colz");
1028    // draw phi error
1029    c1->cd(4);
1030    h2->Draw();
1031    // draw phi abs error
1032    TVirtualPad *p51 = c1->cd(5);
1033    p51->SetLogy();
1034    TH1* h21 = GetHisto("AbsPhiError","","1D",0,0,0,0,0,0);
1035    h21->Draw();
1036    TH1* h22 = GetHisto("Phi68","","1D",0,0,0,0,0,0);
1037    h22->SetFillColor(kRed);
1038    h22->Draw("same");
1039
1040    // draw phi Error of vs Theta Simulated
1041    c1->cd(3);
1042    TH1* h4 = GetHisto("PhiError_vs_Theta_2D","","2D",0,0,0,0,0,0);
1043    h4->Draw("colz");
1044    // draw phi Error of vs Phi Simulated
1045    c1->cd(6);
1046    TH1* h5 = GetHisto("PhiError_vs_Phi_2D","","2D",0,0,0,0,0,0);
1047    h5->Draw("colz");
1048    // draw phi Error of vs S/N
1049    c1->cd(7);
1050    TH1* h6 = GetHisto("PhiError_vs_SignalToNoise_2D","","2D",0,0,0,0,0,0);
1051    h6->Draw("colz");
1052    // draw Phi Error of vs Energy
1053    c1->cd(8);
1054    TH1* h7 = GetHisto("PhiError_vs_EnergySim_2D","","2D",0,0,0,0,0,0);
1055    h7->Draw("colz");
1056    // draw Phi Error of vs Energy
1057    c1->cd(9);
1058    TH1* h8 = GetHisto("PhiError_vs_Distance_2D","","2D",0,0,0,0,0,0);
1059    h8->Draw("colz");
1060
1061    for (it=fCuts.begin(); it != fCuts.end(); it++) {
1062        TCanvas *c = fPhiViewer->FindCanvas((it->first).c_str());
1063        if(!c )  c = fPhiViewer->AddTab(cname+(it->first).c_str(),cname+(it->first).c_str());
1064        c->Clear();
1065        c->cd();
1066        c->Divide(3,3);
1067        TString name, title;
1068
1069        // draw simulated phi distribution: all events and superimposed in the given bin
1070        c->cd(1);
1071        TH1* h1 = GetHisto("ThetaSimulated_1D","Theta Simulated","1D",0,0,0,0,0,0);
1072        h1->DrawCopy();
1073
1074        name = "ThetaSimulated"  + it->first;
1075        TH1* hTot = GetHisto(name+"_1D","","1D",0,0,0,0,0,0);
1076        TH1* hCopyTot = hTot->DrawCopy("same");
1077        hCopyTot->SetFillColor(kBlue);
1078        hCopyTot->SetLineColor(kBlue);
1079
1080        name  = "SimulatedTheta"+it->first;
1081        TH1* h2 = GetHisto(name,"","1D",0,0,0,0,0,0);
1082        TH1* hCopyRe = h2->DrawCopy("same");
1083        hCopyRe->SetXTitle(" #theta_{sim} [deg]");
1084        hCopyRe->SetFillColor(kRed);
1085        hCopyRe->SetLineColor(kRed); 
1086
1087        TLegend *leg = new TLegend(0.65,0.75,0.98,0.98,"Events","brNDC");
1088        leg->SetBorderSize(1);
1089        leg->SetTextAlign(12);
1090        leg->AddEntry(hCopyTot,Form("%d simulated",(Int_t)hCopyTot->GetEntries()),"F");
1091        leg->AddEntry(hCopyRe,Form("%d reconstructed",(Int_t)hCopyRe->GetEntries()),"F");
1092        leg->SetBit(kCanDelete);
1093        leg->Draw();   
1094
1095        // draw simulated phi vs reconstructed phi
1096        c->cd(2);
1097        name  = "PhiSim_vs_Rec"  + it->first;
1098        TH1* h3 = GetHisto(name+"_2D","","2D",0,0,0,0,0,0);
1099        h3->Draw("colz");
1100        h3->SetXTitle("#varphi_{sim} [deg]");
1101        h3->SetYTitle("#varphi_{rec} [deg]");
1102       
1103        // draw theta error
1104        c->cd(4);
1105        name = "PhiError" + it->first;
1106        TH1 *hT = GetHisto(name+"_1D","","1D",0,0,0,0,0,0);
1107        hT->Draw();
1108        hT->SetXTitle("#Delta #varphi [deg]");
1109        hT->SetYTitle("Entries");
1110       
1111        // draw abs theta error
1112        TVirtualPad *p5 = c->cd(5);
1113        p5->SetLogy();
1114        name = "AbsPhiError" + it->first;
1115        TH1 *hTabs = GetHisto(name+"_1D","","1D",0,0,0,0,0,0);
1116        hTabs->Draw();
1117        hTabs->SetXTitle("|#Delta #varphi| [deg]");
1118        hTabs->SetYTitle("Entries");
1119        name = "Phi68" + it->first;
1120        TH1 *h68 = GetHisto(name+"_1D","","1D",0,0,0,0,0,0);
1121        h68->SetFillColor(kRed);
1122        h68->Draw("same");
1123       
1124        // draw phi Error of vs Theta Simulated
1125        c->cd(3);
1126        name  = "PhiError_vs_Theta" + it->first;
1127        TH1* h4 = GetHisto(name+"_2D","","2D",0,0,0,0,0,0);
1128        h4->SetXTitle("#theta_{sim} [deg]");
1129        h4->SetYTitle("#Delta #varphi [deg]");
1130        h4->Draw("colz");
1131       
1132        // draw phi Error of vs Theta Simulated
1133        c->cd(6);
1134        name  = "PhiError_vs_Phi" + it->first;
1135        TH1* h7 = GetHisto(name+"_2D","","2D",0,0,0,0,0,0);
1136        h7->SetXTitle("#varphi_{sim} [deg]");
1137        h7->SetYTitle("#Delta #varphi [deg]");
1138        h7->Draw("colz");
1139       
1140        // draw Phi Error of vs S/N
1141        c->cd(7);
1142        name  = "PhiError_vs_SignalToNoise" + it->first;
1143        TH1* h5 = GetHisto(name+"_2D","","2D",0,0,0,0,0,0);
1144        h5->SetXTitle("S/N");
1145        h5->SetYTitle("#Delta #varphi [deg]");
1146        h5->Draw("colz");
1147       
1148        // draw Theta Error of vs Energy
1149        c->cd(8);
1150        name  = "PhiError_vs_EnergySim" + it->first;
1151        TH1* h6 = GetHisto(name+"_2D","","2D",0,0,0,0,0,0);
1152        h6->SetXTitle("Simulated log_{10}Energy");        h6->SetYTitle("#Delta #varphi [deg]");
1153        h6->Draw("colz");
1154       
1155        // draw Theta Error of vs Energy
1156        c->cd(9);
1157        name  = "PhiError_vs_Distance" + it->first;
1158        TH1* h8 = GetHisto(name+"_2D","","2D",0,0,0,0,0,0);
1159        h8->SetXTitle("Distance in FOV [km]");
1160        h8->SetYTitle("#Delta #theta [deg]");
1161        h8->Draw("colz");
1162    }
1163    cname = "Summary";
1164    TCanvas *c2 = fPhiViewer->FindCanvas(cname);
1165    if(!c2 )  c2 = fPhiViewer->AddTab(cname,cname);
1166    c2->Clear();
1167    c2->cd();
1168    c2->Divide(2,2);
1169
1170    Int_t cuts_size = fCuts.size();
1171    Double_t theta[] = {0.,15.,30.,45.,60.,75.};
1172    Double_t phi68[(const Int_t)cuts_size];
1173    Double_t psi68[(const Int_t)cuts_size];
1174    Int_t i = 0;
1175    for (it=fCuts.begin(); it != fCuts.end(); it++) {
1176        TString hname = "Phi68" + it->first;
1177        TH1* histo = GetHisto(hname+"_1D","","1D",0,0,0,0,0,0);
1178        Int_t nbins = histo->GetNbinsX();
1179        phi68[i] = histo->GetBinCenter(nbins);
1180        hname = "Psi68" + it->first;
1181        histo = GetHisto(hname+"_1D","","1D",0,0,0,0,0,0);
1182        nbins = histo->GetNbinsX();
1183        psi68[i] = histo->GetBinCenter(nbins);
1184        i++;
1185    }
1186   
1187    c2->cd(1);
1188    TGraph *dirRes = new TGraph(6,theta,psi68);
1189    dirRes->SetName("DirectionResolution");
1190    dirRes->SetTitle("Direction Resolution");
1191    dirRes->GetXaxis()->SetTitle("#theta [deg]");
1192    dirRes->GetYaxis()->SetTitle("#psi_{68} [deg]");
1193    dirRes->SetMarkerStyle(22);
1194    dirRes->SetMarkerColor(kRed);
1195    dirRes->Draw("AP");
1196
1197    c2->cd(2);   
1198    TGraph *phiRes = new TGraph(6,theta,phi68);
1199    phiRes->SetName("PhiResolution");
1200    phiRes->SetTitle("Phi Resolution");
1201    phiRes->GetXaxis()->SetTitle("#theta [deg]");
1202    phiRes->GetYaxis()->SetTitle("#varphi_{68} [deg]");
1203    phiRes->SetMarkerStyle(21);
1204    phiRes->SetMarkerColor(kBlue);
1205    phiRes->Draw("AP");
1206
1207}
1208
1209//_____________________________________________________________________________
1210void RecoTreePainter::DrawEnergy() {
1211    //
1212    //
1213    //
1214    fEnergyViewer = new EGViewer();
1215    fEnergyViewer->SetName("EnergyReconstuction");
1216    fEnergyViewer->SetTitle("Energy Reconstuction");
1217    fEnergyViewer->SetWindowName("Energy Reconstuction");
1218
1219    map <string, Bool_t>::iterator it;
1220    for (it=fCuts.begin(); it != fCuts.end(); it++) {
1221        TCanvas *c = fEnergyViewer->FindCanvas((it->first).c_str());
1222        if(!c )  c = fEnergyViewer->AddTab((it->first).c_str(),(it->first).c_str());
1223        c->Clear();
1224        c->cd();
1225        c->Divide(3,3);
1226        TString name, title;
1227
1228        //
1229        // draw simulated theta distribution: all events and superimposed in
1230        // the given bin
1231        //
1232        c->cd(1);
1233        TH1* h1 = GetHisto("ThetaSimulated_1D","Theta Simulated","1D",0,0,0,0,0,0);
1234        h1->DrawCopy();
1235       
1236        name = "ThetaSimulated"  + it->first;
1237        TH1* hTot = GetHisto(name+"_1D","","1D",0,0,0,0,0,0);       
1238        TH1* hCopyTot = hTot->DrawCopy("same");
1239        hCopyTot->SetFillColor(kBlue);
1240        hCopyTot->SetLineColor(kBlue);
1241
1242
1243        name  = "SimulatedTheta"+it->first;
1244        TH1* h2 = GetHisto(name,"","1D",0,0,0,0,0,0);
1245        TH1* hCopyRe = h2->DrawCopy("same");
1246        hCopyRe->SetXTitle("Simulated #{Theta}/deg");
1247        hCopyRe->SetFillColor(kRed);
1248        hCopyRe->SetLineColor(kRed);
1249
1250        name  = "EnergySim_vs_Rec"  + it->first;
1251        TH1* h3 = GetHisto(name+"_2D","","2D",0,0,0,0,0,0);
1252        TLegend *legRec = new TLegend(0.65,0.75,0.98,0.98,"Entries","brNDC");
1253        legRec->SetBorderSize(1);
1254        legRec->SetTextAlign(12);
1255        legRec->AddEntry(hCopyTot,Form("%d simulated",(Int_t)hCopyTot->GetEntries()),"F");
1256        legRec->AddEntry(hCopyRe,Form("%d reconstructed",(Int_t)h3->GetEntries()),"F");
1257        legRec->SetBit(kCanDelete);
1258        legRec->Draw();
1259
1260        //
1261        // draw simulated energy vs reconstructed energy
1262        //
1263        c->cd(2);
1264
1265        h3->Draw("colz");
1266        h3->SetXTitle("Simulated log_{10}Energy");
1267        h3->SetYTitle("Reconstructed log_{10}Energy");
1268
1269
1270        //
1271        // draw Relative Error of vs Energy Simulated
1272        //
1273        c->cd(5);
1274        name  = "EnergyErrorRelative_vs_Energy" + it->first;
1275        TH1* h4 = GetHisto(name+"_2D","","2D",0,0,0,0,0,0);
1276        h4->SetXTitle("Simulated log_{10}Energy");
1277        h4->SetYTitle("#Delta E/E");
1278        h4->Draw("colz");
1279
1280        //
1281        // draw simulated energy superimposed by reconstructed energy
1282        //
1283        c->cd(4);
1284        name = "EnergySimulated"  + it->first;
1285        TH1* h5 = GetHisto(name+"_1D","","1D",0,0,0,0,0,0);
1286        h5->Draw();
1287        h5->SetLineWidth(2);
1288        h5->SetMinimum(0);
1289        name  = "ReconstructedEnergy:"+it->first;
1290        TH1* h6 = GetHisto(name,"","1D",0,0,0,0,0,0);
1291        h6->SetMarkerStyle(29);
1292        h6->SetMarkerSize(0.3);
1293        h6->Draw("EPsame");
1294        h6->SetXTitle("log_{10}Energy");
1295        h6->SetLineColor(kRed);
1296        h5->Rebin(4);
1297        h6->Rebin(4);
1298
1299        //
1300        // draw efficiency as a function of energy
1301        //
1302        c->cd(3);
1303        name = "AcceptedEventsRatio:"+it->first;
1304        TH1* h71 = GetHisto(name,"","1D",0,0,0,0,0,0);
1305        h71->Rebin(4);
1306        h71->Divide(h5);
1307        h71->SetMarkerStyle(29);
1308        h71->SetMarkerColor(kBlue);
1309        h71->Draw("P");
1310        h71->SetXTitle("log_{10}Energy");
1311        h71->SetLineColor(kBlue);
1312        h71->SetMaximum(1.1);
1313
1314
1315        //
1316        // draw Relative Error of vs Signal/Noise ratio
1317        //
1318        c->cd(6);
1319        name  = "EnergyErrorRelative_vs_SignalToNoise" + it->first;
1320        title = "EnergyError Relative vs S/N"+ it->first;
1321        TH2D* h8 = (TH2D*)GetHisto(name+"_2D",title,"2D",0,0,0,0,0,0);
1322        h8->Draw("colz");
1323        h8->SetXTitle("signal/noise");
1324        h8->SetYTitle("#Delta E/E");
1325
1326        //
1327        // draw Signal/Noise ratio as a function of Energy simulated
1328        //
1329        c->cd(7);
1330        name  = "SignalToNoise_vs_EnergySim" + it->first;
1331        title = "S/N vs Energy"+ it->first;
1332        TH2D* h9 = (TH2D*)GetHisto(name+"_2D",title,"2D",0,0,0,0,0,0);
1333        h9->Draw("colz");
1334        h9->SetXTitle("log_{10}Energy");
1335        h9->SetYTitle("signal/noise");
1336
1337        //
1338        // draw sigma and mean of the energy resolution as a function of
1339        // signal/noise
1340        //
1341        c->cd(8);
1342        name  = "EnergyErrorRelative_vs_SignalToNoise" + it->first;
1343        TH2D* h10 = (TH2D*)h8->Clone();
1344        h10->RebinX(4);
1345        h10->FitSlicesY();
1346        TH1D* h8_1 = (TH1D*)gDirectory->Get(name+"_2D_1");
1347        TH1D* h8_2 = (TH1D*)gDirectory->Get(name+"_2D_2");
1348        h8_1->SetMarkerColor(kBlue);
1349        h8_1->SetMarkerStyle(23);
1350        h8_1->SetMarkerSize(1.0);
1351        h8_1->SetXTitle("Signal/Noise");
1352        h8_1->SetYTitle("");
1353        h8_1->SetTitle("Energy resolution");
1354        h8_2->SetMarkerColor(kRed);
1355        h8_2->SetMarkerStyle(29);
1356        h8_2->SetMarkerSize(1.0);
1357
1358        TLegend *leg = new TLegend(0.68,0.75,0.98,0.95);
1359        leg->SetBit(kCanDelete);
1360        leg->SetTextFont(72);
1361
1362        if (h8_1) {
1363            h8_1->Draw("PE");
1364            h8_1->SetMaximum(2);
1365            h8_1->SetMinimum(-1);
1366            h8_1->SetStats(kFALSE);
1367            leg->AddEntry(h8_1,"#LTE#GT/E","p");
1368        }
1369        if (h8_2) {
1370            h8_2->Draw("PEsame");
1371            h8_2->SetStats(kFALSE);
1372            leg->AddEntry(h8_2,"#Delta E/E","p");
1373        }
1374
1375        leg->SetBorderSize(1);
1376        leg->Draw();
1377
1378        /*   name  = "StatisticalExpectation"+it->first;
1379             TH1* h11 = GetHisto(name,title,"1D",100,0,100,0,0,0);
1380             h11->Draw("same");
1381             h11->SetLineColor(kBlack);
1382         */
1383
1384        //
1385        // draw sigma and mean of the energy resolution as a function of energy
1386        //
1387        c->cd(9);
1388        name  = "EnergyErrorRelative_vs_Energy" + it->first;
1389        TH2D* h12 = (TH2D*)GetHisto(name+"_2D",title,"2D",100,left_edge["Energy"],right_edge["Energy"], 100,-2,2);
1390        TH2D* h13 = (TH2D*)h12->Clone();
1391        h13->RebinX(10);
1392        h13->FitSlicesY();
1393        TH1D* h13_1 = (TH1D*)gDirectory->Get(name+"_2D_1");
1394        TH1D* h13_2 = (TH1D*)gDirectory->Get(name+"_2D_2");
1395        h13_1->SetMarkerColor(kBlue);
1396        h13_1->SetMarkerStyle(23);
1397        h13_1->SetXTitle("log_{10}Energy");
1398        h13_1->SetYTitle("");
1399        h13_1->SetTitle("Energy resolution");
1400
1401        h13_2->SetMarkerColor(kRed);
1402        h13_2->SetMarkerStyle(29);
1403
1404        leg = new TLegend(0.68,0.75,0.98,0.95);
1405        leg->SetBit(kCanDelete);
1406        leg->SetTextFont(72);
1407
1408        if (h13_1) {
1409            h13_1->Draw("PE");
1410            h13_1->SetMaximum(2);
1411            h13_1->SetMinimum(-1);
1412            h13_1->SetStats(kFALSE);
1413            leg->AddEntry(h13_1,"#LTE#GT/E","p");
1414        }
1415        if (h13_2) {
1416            h13_2->Draw("PEsame");
1417            h13_2->SetStats(kFALSE);
1418            leg->AddEntry(h13_2,"#Delta E/E","p");
1419        }
1420
1421        leg->SetBorderSize(1);
1422        leg->Draw();
1423
1424    }
1425    TCanvas *c = fEnergyViewer->FindCanvas("c_summary");
1426    if(!c )  c = fEnergyViewer->AddTab("c_summary","Summary");
1427    c->Clear();
1428    c->cd();
1429    c->Divide(3,2);
1430    // draw sigma of energy resolution for all theta angles as a function of energy
1431    c->cd(1);
1432    TLegend *legend1 = new TLegend(0.7,0.6,0.95,0.98);
1433    legend1->SetTextFont(72);
1434    legend1->SetBit(kCanDelete);
1435    //    legend1->SetTextSize(0.02);
1436    Int_t cut(0);
1437
1438    THStack* eRelErr = new THStack("EnergyErrorRelative_vs_Energy","#sigma of Energy Resolution");
1439    eRelErr->SetBit(kCanDelete);
1440
1441    TH1* hCopy = NULL;
1442    for (it=fCuts.begin(); it != fCuts.end(); it++) {
1443        TString name  = "EnergyErrorRelative_vs_Energy" + it->first;
1444        TH1D* h13_2 = (TH1D*)gDirectory->Get(name+"_2D_2");
1445        cut++;
1446        hCopy = (TH1*)h13_2->Clone();
1447        hCopy->SetTitle("#sigma of Energy Resolution");
1448        hCopy->SetMarkerSize(0.5);
1449        hCopy->SetMarkerStyle(20+cut);
1450        hCopy->SetMarkerColor(cut+1);
1451        hCopy->SetStats(0);
1452        legend1->AddEntry(hCopy,TString(it->first).Data(),"p");
1453        eRelErr->Add(hCopy,"P");
1454    }
1455    eRelErr->Draw("nostack");
1456    eRelErr->GetXaxis()->SetTitle(hCopy->GetXaxis()->GetTitle());
1457    eRelErr->GetYaxis()->SetTitle(hCopy->GetYaxis()->GetTitle());
1458    legend1->SetBorderSize(1);
1459    legend1->Draw();
1460
1461    //
1462    // draw sigma of energy resolution for all theta angles as a function of signal/noise
1463    //
1464    c->cd(2);
1465    TLegend *legend2 = new TLegend(0.7,0.55,0.95,0.98);
1466    legend2->SetTextFont(72);
1467    legend2->SetBit(kCanDelete);
1468    //    legend1->SetTextSize(0.02);
1469    cut = 0;
1470    eRelErr = new THStack("EnergyErrorRelative_vs_SignalToNoise","#sigma of Energy Resolution");
1471    eRelErr->SetBit(kCanDelete);
1472    for (it=fCuts.begin(); it != fCuts.end(); it++) {
1473        TString name  = "EnergyErrorRelative_vs_SignalToNoise" + it->first;
1474        TH1D* h13_2 = (TH1D*)gDirectory->Get(name+"_2D_2");
1475        cut++;
1476
1477        hCopy = (TH1*)h13_2->Clone();
1478        hCopy->SetTitle("#sigma of Energy Resolution");
1479        hCopy->SetMarkerSize(0.5);
1480        hCopy->SetStats(kFALSE);
1481        hCopy->SetMarkerStyle(20+cut);
1482        hCopy->SetMarkerColor(cut+1);
1483        legend2->AddEntry(hCopy,TString(it->first).Data(),"p");
1484        eRelErr->Add(hCopy);
1485    }
1486    eRelErr->Draw("nostack");
1487    eRelErr->GetXaxis()->SetTitle(hCopy->GetXaxis()->GetTitle());
1488    eRelErr->GetYaxis()->SetTitle(hCopy->GetYaxis()->GetTitle());
1489    legend2->SetBorderSize(1);
1490    legend2->Draw();
1491
1492    //
1493    // draw efficiencies as a function of Energy for various bins in Theta
1494    //
1495    c->cd(3);
1496    TLegend *legend3 = new TLegend(0.7,0.6,0.95,0.98);
1497    legend3->SetTextFont(72);
1498    legend3->SetBit(kCanDelete);
1499
1500    eRelErr = new THStack("AcceptedEventsRatio","Fraction of accepted events");
1501    eRelErr->SetBit(kCanDelete);
1502
1503    cut = 0;
1504    for (it=fCuts.begin(); it != fCuts.end(); it++) {
1505        TString name = "AcceptedEventsRatio:"+it->first;
1506        TH1* h = GetHisto(name,"","1D",0,0,0,0,0,0);
1507        cut++;
1508        hCopy = (TH1*)h->Clone();
1509        hCopy->SetMarkerSize(0.5);
1510        hCopy->SetTitle("Fraction of accepted events");
1511        hCopy->SetStats(kFALSE);
1512        hCopy->SetMarkerStyle(20+cut);
1513        hCopy->SetMarkerColor(cut+1);
1514        legend3->AddEntry(hCopy,TString(it->first).Data(),"p");
1515        eRelErr->Add(hCopy,"P");
1516    }
1517    eRelErr->Draw("nostack");
1518    eRelErr->GetXaxis()->SetTitle(hCopy->GetXaxis()->GetTitle());
1519    eRelErr->GetYaxis()->SetTitle(hCopy->GetYaxis()->GetTitle());
1520    eRelErr->SetMaximum(eRelErr->GetMaximum("nostack"));
1521    eRelErr->SetMinimum(eRelErr->GetMinimum("nostack"));
1522    legend3->SetBorderSize(1);
1523    legend3->Draw();
1524
1525    //
1526    // draw reconstructed altitude as a function of Theta
1527    //
1528    c->cd(4);
1529    TH2D* h;
1530    h = (TH2D*)GetHisto("Altitude_vs_Theta","Altitude vs Theta","2D",0,0,0,0,0,0);
1531    h->SetYTitle("Altitude [km]");
1532    h->SetXTitle("#Theta [deg]");
1533    TH1* hProf;
1534    hProf = h->ProfileX("_pfx",-1,-1,"s");
1535    hProf->SetYTitle("Altitude [km]");
1536    hProf->SetBit(kCanDelete);
1537    hProf->Draw();
1538
1539    c->cd(5);
1540    h = (TH2D*)GetHisto("Quality_vs_Energy","","2D",0,0,0,0,0,0);
1541    h->SetYTitle("Quality");
1542    h->SetXTitle("#Theta [deg]");
1543    h->Rebin2D(5,5);
1544    h->Draw("box");
1545
1546    //
1547    // summary: reconstructed events vs simulated events
1548    //
1549    TH1* h3 = GetHisto("EnergySimulated_1D","","1D",0,0,0,0,0,0);
1550    c->cd(6);
1551    TPaveText *summary = new TPaveText(0.0,0.0,1.,1.,"brNDC");
1552    summary->SetTextAlign(12);
1553    summary->SetFillStyle(0);
1554    summary->SetBorderSize(0);
1555
1556    summary->AddText(Form("%d simulated events",(Int_t)h3->GetEntries()));
1557    summary->AddText(Form("%d reconstructed events",(Int_t)h->GetEntries()));
1558    summary->Draw();
1559
1560
1561}
1562
1563
1564//_____________________________________________________________________________
1565TH1* RecoTreePainter::GetHisto(const char *name, const char *title, Option_t *opt, Int_t xbins, Float_t xmin, Float_t xmax,
1566                                                                                   Int_t ybins, Float_t ymin, Float_t ymax) {
1567    //
1568    //
1569    //
1570    TH1* th = (TH1*)fHistoList->FindObject(name);
1571    TString option(opt);
1572    if ( !title ) title = name;
1573
1574    if (!th) {
1575        if      (option == "1D")    th = new TH1F( name, title, xbins, xmin, xmax);
1576        else if (option == "2D")    th = new TH2F( name, title, xbins, xmin, xmax, ybins, ymin, ymax);
1577        else if (option == "prof")  th = new TProfile( name, title, xbins, xmin, xmax);
1578        else if (option == "prof2") th = new TProfile2D( name, title, xbins, xmin, xmax,ybins, ymin, ymax,-2,2);
1579        else Error("GetHisto", "Option %s is not supported\n",opt);
1580        fHistoList->Add(th);
1581    }
1582
1583    th->SetStats(kFALSE);
1584    return th;
1585}
1586
1587//______________________________________________________________________________
1588void RecoTreePainter::SaveAsPs( Option_t* option, const char* prefix ) {
1589    TString opt(option);
1590    TString pfx = prefix;
1591    if ( !pfx.IsNull() ) pfx += "_";
1592
1593    if ((opt.Contains("A") || opt.Contains("T")) && fThetaViewer )
1594        fThetaViewer->SaveAsPS(pfx+"ThetaPlots.ps");
1595    if ((opt.Contains("A") || opt.Contains("P")) && fPhiViewer )
1596        fPhiViewer->SaveAsPS(pfx+"PhiPlots.ps");
1597    if ((opt.Contains("A") || opt.Contains("E")) && fEnergyViewer )
1598        fEnergyViewer->SaveAsPS(pfx+"EnergyPlots.ps");
1599}
Note: See TracBrowser for help on using the repository browser.