source: JEM-EUSO/esaf_cc_at_lal/packages/reconstruction/modules/shower/energy/src/EnergyViewer.cc @ 114

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

actual version of ESAF at CCin2p3
File size: 27.8 KB
Line 
1// $Id: EnergyViewer.cc 2946 2011-06-25 13:37:26Z mabl $
2// Author: Dmitry.Naumov   2005/09/25
3
4/*****************************************************************************
5 * ESAF: Euso Simulation and Analysis Framework                              *
6 *                                                                           *
7 *  Id: EnergyViewer                                                         *
8 *  Package:  Energy                                                         *
9 *  Coordinator: Dmitry.Naumov                                               *
10 *                                                                           *
11 *****************************************************************************/
12
13//_____________________________________________________________________________
14//
15// EnergyViewer
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#ifdef USE_VISUALIZATION
27
28#include "EnergyViewer.hh"
29#include "EnergyModule.hh"
30#include "RecoPixelData.hh"
31#include "RecoEvent.hh"
32#include "Etypes.hh"
33#include "RecoShowerTrack.hh"
34#include "EPmtGeo.hh"
35#include "EEvent.hh"
36#include "EShower.hh"
37#include "EAtmosphere.hh"
38#include "ETruth.hh"
39#include "EAtmosphereHistoPainter.hh"
40#include "EShowerHistoPainter.hh"
41#include "EConst.hh"
42#include "EGeometry.hh"
43
44#include <TCanvas.h>
45#include <TPad.h>
46#include <TGraph.h>
47#include <TPolyLine.h>
48#include <TPolyLine3D.h>
49#include <TPolyMarker.h>
50#include <TPolyMarker3D.h>
51#include <TStyle.h>
52#include <TColor.h>
53#include <TText.h>
54#include <TF1.h>
55#include <TView.h>
56#include <TProfile.h>
57#include <TSystem.h>
58#include <TH3F.h>
59
60using namespace TMath;
61using namespace sou;
62
63ClassImp(EnergyViewer)
64 
65//_____________________________________________________________________________
66EnergyViewer::EnergyViewer(EnergyModule* fEM) : EsafMsgSource() {
67  //
68  // Constructor
69  //
70    fCanvas           = NULL;
71    fEnergyModule     = fEM;
72    fColNum           = 200;
73    fColorDefined     = kFALSE;
74    fFSViewBuild      = kFALSE; 
75    fDefinePixelsView = kFALSE;
76
77    DefineColors();
78    gStyle->SetCanvasColor(kWhite);
79}
80
81//_____________________________________________________________________________
82EnergyViewer::~EnergyViewer() {
83    //
84    // Destructor
85    //
86    SafeDelete(fColor);
87    SafeDelete(fCanvas);
88}
89
90
91//_____________________________________________________________________________
92void EnergyViewer::InitFilePS() {
93    //
94    //
95    //
96    fCanvas = new TCanvas("EnergyModule");
97    fPsFile = Form("output/debug.run%d.event%d.ps",
98                    fEnergyModule->fEvent->GetHeader().GetRun(),
99                    fEnergyModule->fEvent->GetHeader().GetNum());
100    fCanvas->Print(fPsFile+"[");
101}
102
103//_____________________________________________________________________________
104void EnergyViewer::CloseFilePS() {
105    //
106    //
107    //
108
109    if (!fCanvas) return;
110 
111    fCanvas->Print(fPsFile+"]");
112    gSystem->Exec(Form("rm -rf %s.gz ; gzip -9 %s", fPsFile.Data(), fPsFile.Data()));
113    Msg(EsafMsg::Info) << fPsFile << " is gzipped" << MsgDispatch;
114    fCanvas->Clear();
115}
116
117//_____________________________________________________________________________
118void EnergyViewer::DisplayFSTrack() {
119    //
120    //
121    //
122
123    if (!fCanvas) return;
124
125    fEnergyModule->fGraphTheta->GetFunction("pol2")->SetLineColor(kMagenta);
126    fEnergyModule->fGraphPhi->GetFunction("pol2")->SetLineColor(kMagenta);
127    fEnergyModule->fGraphX->GetFunction("pol1")->SetLineColor(kMagenta);
128    fEnergyModule->fGraphY->GetFunction("pol1")->SetLineColor(kMagenta);
129    fEnergyModule->fGraphXY->GetFunction("pol1")->SetLineColor(kMagenta);
130
131    fCanvas->Clear(); 
132    fCanvas->Divide(1,2); fCanvas->cd(1);
133    SetHistoStyle(fEnergyModule->fHistoTheta);
134    fEnergyModule->fHistoTheta->Draw();
135    fEnergyModule->fGraphTheta->Draw("X");
136    fEnergyModule->fHistoTheta->GetYaxis()->SetTitle("#Theta, [deg]");
137    fEnergyModule->fHistoTheta->GetXaxis()->SetTitle("Gtu Hit");
138   
139    fCanvas->cd(2);
140    SetHistoStyle(fEnergyModule->fHistoPhi);
141    fEnergyModule->fHistoPhi->Draw();
142    fEnergyModule->fGraphPhi->Draw("X");
143    fEnergyModule->fHistoPhi->GetYaxis()->SetTitle("#Phi, [deg]");
144    fEnergyModule->fHistoPhi->GetXaxis()->SetTitle("Gtu Hit");
145   
146    fCanvas->Print(fPsFile);
147   
148    fCanvas->Clear();
149    fCanvas->Divide(1,2); fCanvas->cd(1);
150    SetHistoStyle(fEnergyModule->fHistoX);
151    fEnergyModule->fHistoX->Draw();
152    fEnergyModule->fGraphX->Draw("X");
153    fEnergyModule->fHistoX->GetYaxis()->SetTitle("X, [mm]");
154    fEnergyModule->fHistoX->GetXaxis()->SetTitle("Gtu Hit");
155   
156    fCanvas->cd(2);
157    SetHistoStyle(fEnergyModule->fHistoY);
158    fEnergyModule->fHistoY->Draw();
159    fEnergyModule->fGraphY->Draw("X");
160    fEnergyModule->fHistoY->GetYaxis()->SetTitle("Y, [mm]");
161    fEnergyModule->fHistoY->GetXaxis()->SetTitle("Gtu Hit");
162
163    fCanvas->Print(fPsFile);
164
165    fCanvas->Clear();
166    SetHistoStyle(fEnergyModule->fHistoXY);
167    fEnergyModule->fHistoXY->SetMarkerColor(kGreen);
168    fEnergyModule->fHistoXY->SetMarkerStyle(kMultiply);
169    fEnergyModule->fHistoXY->SetMarkerSize(1.e-3*Min(gPad->GetWh(),gPad->GetWw()));
170    fEnergyModule->fHistoXY->Draw();
171    fEnergyModule->fGraphXY->Draw("X");
172    fCanvas->Print(fPsFile);
173    fCanvas->Clear();
174}
175
176
177//_____________________________________________________________________________
178void EnergyViewer::DisplayRecoveryDebug(TH1* Psf, Float_t IntegratedEfficacy, Float_t time) {
179    //
180    //
181    //
182
183    if (!fCanvas) return;
184    if (ceil(time) - time) return;
185
186    DefinePixelsView();
187
188    Pixels_t OneGtuPixels;
189    map <Int_t, TPolyLine>::iterator itPmt;
190    Pixels_t::iterator itPix;
191
192    TPaveText *cl = NULL;
193    TText *text = new TText();
194
195    fCanvas->Clear();
196
197    if (gPad) {
198        Float_t OffSet = 5; 
199        gPad->Range(fXmin-OffSet,fYmin-OffSet,fXmax+OffSet,fYmax+OffSet);
200    }
201
202    for (itPmt = fMyPixelLines.begin(); itPmt != fMyPixelLines.end(); itPmt++)
203         itPmt->second.Draw();
204
205    fEnergyModule->fHistoXY->SetStats(kFALSE);
206    fEnergyModule->fHistoXY->Draw("same");
207    Psf->Draw("samecolz");
208    if (fEnergyModule->fGraphXY) fEnergyModule->fGraphXY->Draw("X");
209
210    OneGtuPixels.clear();
211    OneGtuPixels = fEnergyModule->fGtuMyPixels[Int_t(time)];
212
213    for (itPix = OneGtuPixels.begin(); itPix != OneGtuPixels.end(); itPix++) {
214         Int_t Count = (*itPix)->GetCounts();
215         Int_t PixUid = (*itPix)->GetPixelId();
216         if (Count > 0) {
217             fMyPixelLines[PixUid].SetFillColor(fPalette[Count]);
218             fMyPixelLines[PixUid].Draw("f");
219             fMyPixelLines[PixUid].Draw();
220             TVector3 PixelCenter = fEnergyModule->fRunpars->PixelCenter(PixUid);
221             Float_t x = PixelCenter.X()-0.2*fEnergyModule->fRunpars->GetPmtData().GetPadSide();
222             Float_t y = PixelCenter.Y()-0.2*fEnergyModule->fRunpars->GetPmtData().GetPadSide();
223             Float_t padMin = Min(gPad->GetWh(),gPad->GetWw());
224             text->SetTextSize(5.e-5*padMin);
225             text->SetTextColor(5);
226             text->DrawText(x,y,Form("%d",Count));
227          }
228    }
229
230    TString pavetext = Form("Efficacies: total = %f mm^{2}, integrated %f mm^{2} at %3.2f ms",Psf->Integral(), IntegratedEfficacy, time);
231
232    DrawLabel(cl,pavetext);
233
234    fCanvas->Print(fPsFile);
235
236    delete cl;
237    delete text;
238}
239
240//_____________________________________________________________________________
241void EnergyViewer::DisplayRecovery() {
242    //
243    //
244    //
245    if (!fCanvas) return;
246    fCanvas->Clear();
247    fCanvas->Divide(2,2); 
248    fCanvas->cd(1);
249
250    fEnergyModule->fHistoOverlap->Draw("histo");
251    fEnergyModule->fHistoOverlap->GetYaxis()->SetTitle("Psf Integrated over Gtu Overlap with pixels, [mm^{2}]");
252    fEnergyModule->fHistoOverlap->GetXaxis()->SetTitle("Gtu Hit");
253   
254    fCanvas->cd(2);
255    fEnergyModule->fHistoIdealOverlap->Draw("histo");
256    fEnergyModule->fHistoIdealOverlap->GetYaxis()->SetTitle("Psf Ideal Overlap with pixels, [mm^{2}]");
257    fEnergyModule->fHistoIdealOverlap->GetXaxis()->SetTitle("Gtu Hit");
258//     map <Int_t, TPolyLine>::iterator itPix;
259//     fCanvas->cd(3);
260//     for (itPix = fMyPixelLines.begin(); itPix != fMyPixelLines.end(); itPix++) {
261//          itPix->second.Draw("f");
262//          itPix->second.Draw();
263//     }
264//     if (fEnergyModule->fGraphXY) fEnergyModule->fGraphXY->Draw("X");
265   
266    fCanvas->Print(fPsFile,"Preview");
267}
268
269//_____________________________________________________________________________
270void EnergyViewer::DefineColors() {
271
272    if (fColorDefined) return;
273
274    Int_t Red, Green, Blue;
275
276    fPalette[0] = 1;
277
278    for (Int_t iColor = 1; iColor < 1000; iColor++) {
279       if (iColor <= fColNum) {
280
281      if (iColor < 0.3*fColNum) {Blue = 1; Green = Red = 0;}
282      else if (iColor > 0.7*fColNum) {Red = 1; Green = Blue = 0;}
283      else {Green = 1; Red = Blue = 0;}
284     
285      if(!gROOT->GetColor(300+iColor))
286        fColor = new TColor(300+iColor,
287                            Red*(0.6+0.4*iColor/fColNum),
288                            Green*(0.6+0.5*iColor/fColNum),
289                            Blue*(0.6+1.*iColor/fColNum),"");
290      else {
291        fColor = gROOT->GetColor(300+iColor);
292        fColor->SetRGB(Red*(0.6+0.4*iColor/fColNum),
293                       Green*(0.6+0.5*iColor/fColNum),
294                       Blue*(0.6+1.*iColor/fColNum));
295      }
296      fPalette[iColor] = 300 + iColor;
297    }
298    else 
299      fPalette[iColor] = 300 + fColNum;
300  }
301
302    fColorDefined = kTRUE;
303}
304
305//_____________________________________________________________________________
306void EnergyViewer::DefinePixelsView() {
307
308    if (!fCanvas || fDefinePixelsView) return;
309
310    PMTs_t::iterator itPmt;
311
312    fMyPixelLines.clear();
313    fPixelLines.clear(); fPixelLines3D.clear(); 
314    fPmtLines.clear();   fPmtLines3D.clear();
315
316    fXmin =  fEnergyModule->fXmin;
317    fXmax =  fEnergyModule->fXmax;
318    fYmin =  fEnergyModule->fYmin;
319    fYmax =  fEnergyModule->fYmax;
320    fZmin =  kHuge;
321    fZmax = -kHuge;
322 
323    Float_t x[5], y[5], p[15];
324
325    // make list of lines to draw around each pixels with signal
326
327    for (itPmt = fEnergyModule->fMyPixels.begin(); itPmt != fEnergyModule->fMyPixels.end(); itPmt++) {
328      for(Int_t iCorner = 0; iCorner < 4; iCorner++) {
329        TVector3 vPix = fEnergyModule->fRunpars->PixelCorner(*itPmt,iCorner);
330        x[iCorner] = vPix.X(); y[iCorner] = vPix.Y();
331      }
332      x[4] = x[0]; y[4] = y[0];
333      fMyPixelLines[*itPmt].SetPolyLine(5,x,y);
334    }
335   
336    // make list of lines to draw around each PMT
337 
338    for (itPmt = fEnergyModule->fPmts.begin(); itPmt != fEnergyModule->fPmts.end(); itPmt++) {
339      for(Int_t iCorner(0); iCorner < 4; iCorner++) {
340        TVector3 vPmt = fEnergyModule->fRunpars->PmtCorner(*itPmt,iCorner);
341
342        p[3*iCorner]   =  x[iCorner] =  vPmt.X();
343        p[3*iCorner+1] =  y[iCorner] =  vPmt.Y();
344        p[3*iCorner+2] = -vPmt.Z();
345
346        if (fXmin >  vPmt.X()) fXmin =  vPmt.X();
347        if (fYmin >  vPmt.Y()) fYmin =  vPmt.Y();
348        if (fZmin > -vPmt.Z()) fZmin = -vPmt.Z();
349        if (fXmax <  vPmt.X()) fXmax =  vPmt.X();
350        if (fYmax <  vPmt.Y()) fYmax =  vPmt.Y();
351        if (fZmax < -vPmt.Z()) fZmax = -vPmt.Z();
352      }
353
354      p[12] = x[4] = x[0]; p[13] = y[4] = y[0]; p[14] = p[2];
355
356      fPmtLines[*itPmt].SetPolyLine(5,x,y);
357      fPmtLines[*itPmt].SetLineColor(20);
358
359      fPmtLines3D[*itPmt].SetPolyLine(5,p);
360      fPmtLines3D[*itPmt].SetLineColor(kBlack);
361    }
362
363    // make list of lines to draw around each pixels
364
365    for (itPmt = fEnergyModule->fPixelsBgr.begin(); itPmt != fEnergyModule->fPixelsBgr.end(); itPmt++) {
366      for(Int_t iCorner(0); iCorner < 4; iCorner++) {
367        TVector3 vPmt = fEnergyModule->fRunpars->PixelCorner(*itPmt,iCorner);
368        p[3*iCorner]   =  x[iCorner] =  vPmt.X();
369        p[3*iCorner+1] =  y[iCorner] =  vPmt.Y();
370        p[3*iCorner+2] = -vPmt.Z();
371      }
372     
373      p[12] = x[4] = x[0]; p[13] = y[4] = y[0]; p[14] = p[2];
374     
375      fPixelLines[*itPmt].SetPolyLine(5,x,y);
376      fPixelLines[*itPmt].SetLineColor(10);
377      fPixelLines[*itPmt].SetFillColor(fPalette[0]);
378   
379      fPixelLines3D[*itPmt].SetPolyLine(5,p);
380      fPixelLines3D[*itPmt].SetLineColor(kWhite);
381    }
382    fDefinePixelsView = kTRUE;
383}
384
385//_____________________________________________________________________________
386void EnergyViewer::DefineFSView() {
387
388  if (!fCanvas || fFSViewBuild) return;
389
390  Float_t x[5], y[5], p[15];
391  Float_t PhiMin, PhiMax, PhiMin1, PhiMax1;
392
393  PMTs_t::iterator itPmt;
394
395  fPmtLinesFS3D.clear();
396
397  fGXmin = fGYmin = fGZmin = PhiMin = PhiMin1 =  kHuge;
398  fGXmax = fGYmax = fGZmax = PhiMax = PhiMax1 = -kHuge;
399
400  for (itPmt = fEnergyModule->fPmts.begin(); itPmt != fEnergyModule->fPmts.end(); itPmt++) {
401    TVector3 vPmt = fEnergyModule->fRunpars->GetPmtGeo(*itPmt)->GetCenter();
402    Float_t Phi = vPmt.Phi();
403    if (Phi < PhiMin) PhiMin = Phi;
404    if (Phi > PhiMax) PhiMax = Phi;
405
406    if (Phi < 0.) Phi += 2*Pi();
407
408    if (Phi < PhiMin1) PhiMin1 = Phi;
409    if (Phi > PhiMax1) PhiMax1 = Phi;
410  }
411
412  Int_t Case = 0;
413
414  if (PhiMax > 0.7*Pi() && PhiMin < -0.7*Pi()) {
415    Case = 1;
416    PhiMin = PhiMin1; PhiMax = PhiMax1;
417  }
418
419  for(Int_t i(1); i <= fEnergyModule->fRunpars->GetNumPmts(); i++) {
420    TVector3 vPmt = fEnergyModule->fRunpars->GetPmtGeo(i)->GetCenter();
421    Float_t phi = vPmt.Phi();
422    if (Case && vPmt.Phi() < 0) phi = vPmt.Phi() + 2*Pi();
423    if (phi > PhiMin - 0.05*Pi() && phi < PhiMax + 0.05*Pi()) {
424      for(Int_t iCorner(0); iCorner < 4; iCorner++) {
425        TVector3 vPmt = fEnergyModule->fRunpars->PmtCorner(i,iCorner);
426       
427        p[3*iCorner]   =  x[iCorner] =  vPmt.X();
428        p[3*iCorner+1] =  y[iCorner] =  vPmt.Y();
429        p[3*iCorner+2] = -vPmt.Z();
430       
431        if (fGXmin >  vPmt.X()) fGXmin =  vPmt.X();
432        if (fGYmin >  vPmt.Y()) fGYmin =  vPmt.Y();
433        if (fGZmin > -vPmt.Z()) fGZmin = -vPmt.Z();
434        if (fGXmax <  vPmt.X()) fGXmax =  vPmt.X();
435        if (fGYmax <  vPmt.Y()) fGYmax =  vPmt.Y();
436        if (fGZmax < -vPmt.Z()) fGZmax = -vPmt.Z();
437      }
438     
439      p[12] = x[4] = x[0]; p[13] = y[4] = y[0]; p[14] = p[2];
440     
441      fPmtLinesFS3D[i].SetPolyLine(5,p);
442      fPmtLinesFS3D[i].SetLineColor(kBlack);
443    }
444  }
445
446  fFSViewBuild = kTRUE;
447}
448
449//_____________________________________________________________________________
450void EnergyViewer::DisplayPixels() {
451
452  if (!fCanvas) return;
453  DefinePixelsView();
454
455  Int_t Count, PixUid;
456  Pixels_t OneGtuPixels;
457  map <Int_t,Int_t> PixelsSignal;
458  Pixels_t::iterator itPix;
459  GtuPixels_t::iterator itGtu;
460  map <Int_t, TPolyLine>::iterator itPmt;
461
462  PixelsSignal.clear();
463  fCanvas->Clear();
464
465  if (gPad) {
466    Float_t OffSet = 5; 
467    gPad->Range(fXmin-OffSet,fYmin-OffSet,fXmax+OffSet,fYmax+OffSet);
468  }
469
470  // Draw PMTs
471
472  for (itPmt = fPmtLines.begin(); itPmt != fPmtLines.end(); itPmt++) {
473    itPmt->second.Draw("f");
474    itPmt->second.Draw();
475  }
476
477  // Integrate events over all pixels
478
479  for (itPix = fEnergyModule->fActivePixels.begin(); itPix != fEnergyModule->fActivePixels.end(); itPix++){
480    PixelsSignal[(*itPix)->GetPixelId()] += (*itPix)->GetCounts();
481  }
482
483  // Draw integral signals
484
485  for (itPmt = fPixelLines.begin(); itPmt != fPixelLines.end(); itPmt++) {
486    Count = PixelsSignal[itPmt->first];
487    itPmt->second.SetFillColor(fPalette[Count]);
488    itPmt->second.Draw("f");
489    itPmt->second.Draw();
490  }
491
492  if (fEnergyModule->fGraphXY) fEnergyModule->fGraphXY->Draw("X");
493  fEnergyModule->fHistoXY->SetStats(0);
494  fEnergyModule->fHistoXY->Draw("SAME");
495
496  // show the legend on top of the page
497 
498  TPaveText *cl = NULL;
499  TString pavetext = Form("Integral signal over all GTUs in [%d,%d] microsecond interval",fEnergyModule->fGtuMin,fEnergyModule->fGtuMax);
500
501  DrawLabel(cl,pavetext);
502
503  fCanvas->Print(fPsFile);
504
505  SafeDelete(cl);
506
507  // Restore default color for inactive pixels
508
509  for (itPmt = fPixelLines.begin(); itPmt != fPixelLines.end(); itPmt++) 
510    itPmt->second.SetFillColor(fPalette[0]);
511
512  if (fEnergyModule->fDebugPs >= 2) {
513
514    // Loop for all Gtu and fill the map of all active pixels
515 
516    for (itGtu = fEnergyModule->fGtuMyPixels.begin(); itGtu != fEnergyModule->fGtuMyPixels.end(); itGtu++) {
517
518      fCanvas->Clear();
519
520      // Draw PMTs
521
522      for (itPmt = fPmtLines.begin(); itPmt != fPmtLines.end(); itPmt++) {
523        itPmt->second.Draw("f");
524        itPmt->second.Draw();
525      }
526
527      // Define colors and signals in each pixels per one Gtu
528
529      OneGtuPixels.clear(); PixelsSignal.clear();
530      OneGtuPixels = itGtu->second;
531
532      for (itPix = OneGtuPixels.begin(); itPix != OneGtuPixels.end(); itPix++) {
533        Count = (*itPix)->GetCounts();
534        PixUid = (*itPix)->GetPixelId();
535        fPixelLines[PixUid].SetFillColor(fPalette[Count]);
536        fPixelLines[PixUid].SetLineColor(kRed);
537        PixelsSignal[PixUid] = Count;
538      }
539
540      // Draw Pixels
541
542      TText *text = new TText();
543
544      for (itPmt = fPixelLines.begin(); itPmt != fPixelLines.end(); itPmt++) {
545        itPmt->second.Draw("f");
546        itPmt->second.Draw();
547        Count = PixelsSignal[itPmt->first];
548        if (Count > 0) {
549          TVector3 PixelCenter = fEnergyModule->fRunpars->PixelCenter(itPmt->first);
550          Float_t x = PixelCenter.X()-0.2*fEnergyModule->fRunpars->GetPmtData().GetPadSide();
551          Float_t y = PixelCenter.Y()-0.2*fEnergyModule->fRunpars->GetPmtData().GetPadSide();
552          Float_t padMin = Min(gPad->GetWh(),gPad->GetWw());
553          text->SetTextSize(5.e-5*padMin);
554          text->SetTextColor(5);
555          text->DrawText(x,y,Form("%d",Count));
556        }
557      }
558
559      if (fEnergyModule->fGraphXY) fEnergyModule->fGraphXY->Draw("X");
560
561      TString pavetext = Form("Signal at GTU = %d",itGtu->first);
562      DrawLabel(cl,pavetext);
563
564      fCanvas->Print(fPsFile);
565
566      SafeDelete(text);
567      SafeDelete(cl);
568
569      // Restore default color for inactive pixels
570
571      for (itPmt = fPixelLines.begin(); itPmt != fPixelLines.end(); itPmt++) {
572        itPmt->second.SetFillColor(fPalette[0]);
573        itPmt->second.SetLineColor(kWhite);
574      }
575    }
576  }
577}
578
579//_____________________________________________________________________________
580void EnergyViewer::DisplayPixels3D(TString DrawSurf) {
581
582    if (!fCanvas) return;
583    Int_t Count;
584    Float_t Xmin, Xmax, Ymin, Ymax, Zmin, Zmax;
585    map <Int_t,Int_t> PixelsSignal;
586    map<Int_t,TPolyLine3D> Pmts3D;
587    Pixels_t::iterator itPix;
588    map <Int_t, TPolyLine3D>::iterator itPmt;
589
590    fCanvas->Clear();
591    PixelsSignal.clear();
592    Pmts3D.clear();
593   
594    if (DrawSurf == "FocalSurf") {
595      DefineFSView();
596      Xmin = fGXmin; Xmax = fGXmax; Ymin = fGYmin; Ymax = fGYmax;
597      Zmin = fGZmin; Zmax = fGZmax;
598      Pmts3D = fPmtLinesFS3D;
599    }
600    else if (DrawSurf == "ClusterSurf") {
601      DefinePixelsView();
602      Xmin = fXmin; Xmax = fXmax; Ymin = fYmin; Ymax = fYmax;
603      Zmin = fZmin; Zmax = fZmax;
604      Pmts3D = fPmtLines3D;
605    } else {
606        throw invalid_argument("DisplayPixels3D expects either "
607                "'FocalSurf' or 'ClusterSurf' as parameter.");
608    }
609
610    TView *view = NULL;
611#if ( ROOT_VERSION_CODE <= ROOT_VERSION(5,16,00) )
612        view = new TView(11);
613#else
614        view =TView::CreateView(1);
615#endif
616
617    view->ShowAxis();
618    view->SetRange(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
619
620    Float_t xm = 0.5*(Xmax+Xmin);
621    Float_t ym = 0.5*(Ymax+Ymin);
622
623    if (xm > 0 && ym > 0) view->RotateView(-135,60);
624    else if (xm < 0 && ym < 0) view->RotateView(45,60);
625    else if (xm < 0 && ym > 0) view->RotateView(-45,60);
626    else if (xm > 0 && ym < 0) view->RotateView(135,60);
627
628    // Draw PMTs
629
630    for (itPmt = Pmts3D.begin(); itPmt != Pmts3D.end(); itPmt++)
631      itPmt->second.Draw();
632
633    // Calculate number of counts for each pixel
634
635    for (itPix = fEnergyModule->fActivePixels.begin(); 
636         itPix != fEnergyModule->fActivePixels.end(); itPix++)
637      PixelsSignal[(*itPix)->GetPixelId()] += (*itPix)->GetCounts();
638
639    // Draw integral signals
640
641    for (itPmt = fPixelLines3D.begin(); itPmt != fPixelLines3D.end(); itPmt++) {
642      Count = PixelsSignal[itPmt->first];
643      if (Count) {
644        itPmt->second.SetLineColor(fPalette[Count]);
645        itPmt->second.Draw();
646      }
647    }
648
649    // show the legend on top of the page
650
651    TPaveText *cl = NULL;
652    TString pavetext = Form("Integral signal over all GTUs in [%d,%d] microsecond interval",fEnergyModule->fGtuMin,fEnergyModule->fGtuMax);
653   
654    DrawLabel(cl,pavetext);
655
656    fCanvas->Print(fPsFile);
657
658    // Restore default color for inactive pixels
659
660    for (itPmt = fPixelLines3D.begin(); itPmt != fPixelLines3D.end(); itPmt++)
661      itPmt->second.SetLineColor(kWhite);
662
663    SafeDelete(cl);
664    SafeDelete(view);
665}
666
667//_____________________________________________________________________________
668void EnergyViewer::DisplayCherenkov(TH1F *htime, Int_t nPoints, Float_t x[], Float_t y[]) {
669
670  if (!fCanvas) return;
671  fCanvas->Clear();
672
673  TList *functions = htime->GetListOfFunctions();                       
674  TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");
675
676  htime->Draw("histo");
677
678  TGraph *gr1 = new TGraph(nPoints,x,y);
679  gr1->SetMarkerColor(kBlue);
680  gr1->SetMarkerStyle(21);
681  gr1->Draw("P");
682  if (pm) pm->Draw();
683     
684  // fCanvas->Print(fPsFileTime+"(");
685  fCanvas->Print(fPsFile);
686  SafeDelete(gr1);
687}
688
689
690//_____________________________________________________________________________
691void EnergyViewer::DisplayTimeFit(TH1F *htime, TF1 *func, TF1 *bg)
692{
693    if (!fCanvas) return;
694    fCanvas->Clear();
695
696    htime->SetStats(0);
697    htime->Draw("PE");
698    if (func) func->Draw("same");
699    if (bg)   bg->Draw("same");
700
701    TPaveText *cl = new TPaveText(0.7,0.85,1.2,1.,"brNDC");
702    cl->SetTextFont(72); cl->SetTextSize(0.02);
703    cl->SetFillColor(0); cl->SetTextAlign(12);
704
705    // Event Information
706
707    cl->AddText("Truth Information");
708    cl->AddText(Form("Energy = %.2E MeV",
709                     fEnergyModule->fEvent->GetHeader().GetTrueEnergy()));
710    cl->AddText(Form("#Theta = %.2f deg",
711                     fEnergyModule->fEvent->GetHeader().GetTrueTheta()*RadToDeg()));
712    cl->AddText(Form("#phi   = %.2f deg",
713                     fEnergyModule->fEvent->GetHeader().GetTruePhi()*RadToDeg()));
714    cl->AddText("Reco Information");
715    if (!IsNaN(fEnergyModule->fEnergy)) 
716      cl->AddText(Form("Energy = %.2E #pm %2.E MeV",
717                       fEnergyModule->fEnergy,fEnergyModule->fEnergyError));
718    if (!IsNaN(fEnergyModule->fTheta)) 
719      cl->AddText(Form("#Theta = %.2f deg",fEnergyModule->fTheta*RadToDeg()));
720    if (!IsNaN(fEnergyModule->fPhi))
721      cl->AddText(Form("#phi   = %.2f deg",fEnergyModule->fPhi*RadToDeg()));
722    cl->Draw();
723
724    fCanvas->Print(fPsFile);
725
726    SafeDelete(cl);
727
728}
729
730
731//_____________________________________________________________________________
732void EnergyViewer::DrawHisto(TH1F *histo) {   
733    //
734    //
735    //
736    if (!fCanvas) return;
737    fCanvas->Clear();
738
739    histo->Draw();
740    fCanvas->Print(fPsFile);
741}
742
743
744//_____________________________________________________________________________
745void EnergyViewer::Clean() {
746
747    fMyPixelLines.clear();
748    fPixelLines.clear(); fPixelLines3D.clear();
749    fPmtLines.clear();   fPmtLines3D.clear();
750    fPmtLinesFS3D.clear();
751
752    fFSViewBuild      = kFALSE; 
753    fDefinePixelsView = kFALSE;
754}
755 
756//_____________________________________________________________________________
757void EnergyViewer::SetHistoStyle(TH2F *histo) {
758  histo->GetXaxis()->SetLabelFont(63);
759  histo->GetXaxis()->SetLabelSize(15);
760  histo->GetYaxis()->SetLabelFont(63);
761  histo->GetYaxis()->SetLabelSize(15);
762 
763  histo->GetXaxis()->SetTitleFont(73);
764  histo->GetXaxis()->SetTitleSize(15);
765  histo->GetYaxis()->SetTitleFont(73);
766  histo->GetYaxis()->SetTitleSize(15);
767 
768  histo->SetTitleOffset(1.3,"X");
769  histo->SetTitleOffset(1.3,"Y");
770
771  histo->GetXaxis()->SetNdivisions(507,kTRUE);
772  histo->GetYaxis()->SetNdivisions(505,kTRUE);
773
774  histo->SetMarkerStyle(kFullDotMedium);
775  histo->SetMarkerColor(kBlue);
776}
777
778
779//_____________________________________________________________________________
780void EnergyViewer::DrawLabel(TPaveText *cl, TString pavetext) {
781  cl = new TPaveText(0.1,0.96,0.9,1.,"brNDC");
782  cl->SetTextFont(72);
783  cl->SetTextSize(0.03);
784  cl->SetFillColor(kGreen);
785  cl->SetTextAlign(22);
786  cl->SetTextColor(kBlue);
787  cl->AddText(pavetext);
788  cl->Draw();
789}
790
791//_____________________________________________________________________________
792void EnergyViewer::DisplayShowerInfo() {
793    //
794    //
795    //
796    if (!fCanvas) return;
797    fCanvas->Clear();
798    EEvent          *SimuEvent  = fEnergyModule->fEvent->GetSimuEvent();
799    if (!SimuEvent) return;
800    EAtmosphere     *atmosphere = SimuEvent->GetAtmosphere();
801    EShower         *shower     = SimuEvent->GetShower();
802    ETruth          *Truth      = SimuEvent->GetTruth();
803    ERunParameters  *runpars    = SimuEvent->GetRunPars();
804    const EGeometry *detgeom    = SimuEvent->GetGeometry();
805    Double_t GTUlength = runpars->GetPmtData().GetGtuLength();
806
807    if (!atmosphere || !shower || !Truth || !runpars) return;
808
809    TVector3 DetectorPos = detgeom->GetPos();
810    EAtmosphereHistoPainter *gAtmosphereHistoPainter = 
811                             new EAtmosphereHistoPainter( atmosphere,detgeom,GTUlength,Truth,shower);
812
813    fCanvas->Divide(2,2);
814
815    fCanvas->cd(1);
816    gAtmosphereHistoPainter->Draw("alt_Nph_F");
817    fCanvas->cd(2);
818    gAtmosphereHistoPainter->Draw("alt_Yield_F");
819    fCanvas->cd(3);
820    gAtmosphereHistoPainter->Draw("alt_tottrans_details");
821    fCanvas->cd(4);
822    gAtmosphereHistoPainter->Draw("alt_fluo_trans");
823    fCanvas->Print(fPsFile);
824    SafeDelete(gAtmosphereHistoPainter);
825
826    // check reconstructed shower
827    fCanvas->Clear();
828    fCanvas->Divide(2,2);
829
830   
831    TList *list = (TList*)fEnergyModule->fShowerTrack->GetDebugList();
832    if (list->GetEntries()==0) fEnergyModule->fShowerTrack->DoDebugHistos();
833    fCanvas->cd(1);
834    TH1* h1 = (TH1*)list->At(9);
835    h1->Draw();
836    h1->SetStats(0);
837    h1->SetXTitle("altitude, [km]");
838    h1->SetYTitle("Fluorescence yield");
839
840    fCanvas->cd(2);
841    TH1* h2 = (TH1*)list->At(11);
842    h2->Draw();
843    h2->SetStats(0);
844    h2->SetXTitle("altitude, [km]");
845    h2->SetYTitle("Fluorescence {#gamma} at EP");
846
847    fCanvas->cd(3);
848    TH1* h3 = (TH1*)list->At(10);
849    h3->Draw();
850    h3->SetStats(0);
851    h3->SetXTitle("altitude, [km]");
852    h3->SetYTitle("Total transmission");
853
854    fCanvas->cd(4);
855    TPolyMarker3D *track = (TPolyMarker3D*)list->At(12);
856    Double_t Xmin =  fEnergyModule->fShowerTrack->GetXmin();
857    Double_t Ymin =  fEnergyModule->fShowerTrack->GetYmin();
858    Double_t Zmin =  fEnergyModule->fShowerTrack->GetZmin();
859    Double_t Xmax =  fEnergyModule->fShowerTrack->GetXmax();
860    Double_t Ymax =  fEnergyModule->fShowerTrack->GetYmax();
861    Double_t Zmax =  fEnergyModule->fShowerTrack->GetZmax();
862
863    TH3F* h3D = new TH3F("h3D","",2,Xmin,Xmax,2,Ymin,Ymax,2,Zmin,Zmax);
864    h3D->SetStats(0);
865    h3D->Draw();
866    track->SetMarkerColor(kRed);
867    track->SetMarkerStyle(8);
868    track->SetMarkerSize(0.1);
869    track->Draw();
870    fCanvas->Print(fPsFile);
871    SafeDelete(h3D);
872}
873
874//_____________________________________________________________________________
875void EnergyViewer::CompareRecoVersusSimu() {
876    //
877    // This method compares reconstructed information to that simulated on every step of the chain to help to find
878    // problems in the reconstruction or simulation (if any)
879    // 1.Shower (development vs time, Theta)
880    // 2.Fluorescence and Cherenkov production
881    // 3.Atmosphere response
882    // 4.Photons on entrance pupil
883    // 5.Photons on focal surface (or optical adaptor?)
884    // 6.p.e.
885
886    CompareRecoVersusSimu_Shower();
887
888}
889//_____________________________________________________________________________
890void EnergyViewer::CompareRecoVersusSimu_Shower() {
891    TList *list = (TList*)fEnergyModule->fShowerTrack->GetDebugList();
892    if (list->GetEntries()==0) fEnergyModule->fShowerTrack->DoDebugHistos();
893
894    if (!fCanvas) return;
895    fCanvas->Clear();
896
897    fCanvas->Divide(1,2);
898    fCanvas->cd(1);
899    fEnergyModule->fShowerTrack->Debug(fCanvas,0);
900    fCanvas->cd(2);
901    EEvent          *SimuEvent  = fEnergyModule->fEvent->GetSimuEvent();
902    if (!SimuEvent) return;
903    EShower         *shower   = SimuEvent->GetShower();
904    const EGeometry *geometry = SimuEvent->GetGeometry();
905    EDetector *detector = SimuEvent->GetDetector();
906    EShowerHistoPainter *painter = new EShowerHistoPainter(shower,geometry,detector);
907    painter->GetGtuHisto("Shower_Ne_GTU")->Draw();
908    fCanvas->Print(fPsFile);
909
910}
911
912//_____________________________________________________________________________
913void EnergyViewer::DisplayAltitude(TH1F* histo, TF1* fit) {
914    //
915    //
916    //
917    if (!fCanvas) return;
918    fCanvas->Clear();
919    histo->Draw("E");
920    fit->Draw("same");
921    histo->SetXTitle("time, [#mu sec]");
922    fCanvas->Print(fPsFile);
923}
924
925#endif
Note: See TracBrowser for help on using the repository browser.