source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/common/eventviewer/src/EShowerPainter.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: 11.1 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: EShowerPainter.cc 2795 2008-02-11 12:50:20Z naumov $
3// Naumov Dmitry created Jun, 28 2004
4
5#include "EShowerPainter.hh"
6#include "ESystemOfUnits.hh"
7#include "EShower.hh"
8#include "EShowerStep.hh"
9#include "Etypes.hh"
10
11#include "Riostream.h"
12
13#include <TGeoBBox.h>
14#include <TGeoVolume.h>
15#include <TGeoMatrix.h>
16#include <TVector3.h>
17#include <TGeoManager.h>
18#include <TGeoTube.h>
19#include <TPad.h>
20#include <TTimer.h>
21#include <TSeqCollection.h>
22#include <TH2F.h>
23#include <TH3F.h>
24#include <TMath.h>
25#include <TList.h>
26
27using namespace TMath;
28using namespace sou;
29
30ClassImp(EShowerPainter)
31
32//_____________________________________________________________________________
33EShowerPainter::EShowerPainter( EShower *sh ) {
34    //
35    // Constructor
36    //
37
38    fShower = sh;
39
40    fNumFrames = 100;
41    fFrame = -1;
42    fWorldBuild = kFALSE;
43    fBoxDefined = kFALSE;
44    fHistoContainer = new TList();
45    fThreshold = 1.e-3;
46
47    if ( !fShower || fShower->GetNumSteps() <= 0 ) {
48        Info("Build()","Empty Shower ( NumSteps = 0 ). Painter made zombie.");
49        MakeZombie();
50        return;
51    }
52
53}
54
55//_____________________________________________________________________________
56EShowerPainter::~EShowerPainter() {
57    //
58    // Destructor
59    //
60    fHistoContainer->Delete();
61
62}
63//______________________________________________________________________________
64void EShowerPainter::DefineBox() {
65
66    if (IsBoxDefined()) return;
67    Int_t n = fShower->GetNumSteps();
68    Double_t Xmin(kHuge), Xmax(-kHuge), Ymin(kHuge), Ymax(-kHuge), Zmin(kHuge), Zmax(-kHuge), Nmax(1);
69   
70    for (Int_t i=0; i<n; i++) {
71        if (Nmax < fShower->GetStep(i)->GetNumElectrons() )
72            Nmax = fShower->GetStep(i)->GetNumElectrons();
73    }   
74   
75    for (Int_t i=0; i<n; i++) {
76        Double_t xmean = (fShower->GetStep(i)->GetPosi().X() + fShower->GetStep(i)->GetPosf().X())/2/km;
77        Double_t ymean = (fShower->GetStep(i)->GetPosi().Y() + fShower->GetStep(i)->GetPosf().Y())/2/km;
78        Double_t zmean = (fShower->GetStep(i)->GetPosi().Z() + fShower->GetStep(i)->GetPosf().Z())/2/km;
79        Double_t ne    =  fShower->GetStep(i)->GetNumElectrons();
80        if (ne/Nmax > fThreshold) {
81            Xmin = (Xmin < xmean ? Xmin : xmean); 
82            Xmax = (Xmax > xmean ? Xmax : xmean); 
83            Ymin = (Ymin < ymean ? Ymin : ymean); 
84            Ymax = (Ymax > ymean ? Ymax : ymean); 
85            Zmin = (Zmin < zmean ? Zmin : zmean); 
86            Zmax = (Zmax > zmean ? Zmax : zmean); 
87        }
88    }
89
90    fXmin = Xmin;
91    fXmax = Xmax;
92    fYmin = Ymin;
93    fYmax = Ymax;
94    fZmin = Zmin;
95    fZmax = Zmax;
96    fNmax = Nmax;
97    fBoxDefined = kTRUE;
98}
99//______________________________________________________________________________
100void EShowerPainter::Build() {
101    //
102    // Find the size of the box that contains the shower
103    //
104
105    if ( IsZombie() ) return;
106
107    DefineBox();
108    Int_t n = fShower->GetNumSteps();
109
110   
111    fDx = 1.;
112    fDy = 1.;
113    fDz = (fShower->GetStep(0)->GetPosi() - fShower->GetStep(n-1)->GetPosf()).Mag()/2/km;
114   
115    TVector3 initPos  = (fShower->GetStep(0)->GetPosi() + fShower->GetStep(0)->GetPosf())*0.5;
116    TVector3 lastPos  = (fShower->GetStep(n-1)->GetPosi() + fShower->GetStep(n-1)->GetPosf())*0.5;
117    fCenter = (initPos + lastPos)*0.5*(1/km);
118    // build the shower track
119    Double_t kScale = 1./fNmax;
120    fShowerVolume = gGeoManager->MakeBox("fShowerVolume",fVacuum,fDx,fDy,fDz);
121    Double_t dz_shift = 0;
122
123    for (Int_t i=0; i<n; i++) {
124        TVector3 stepCenter = (fShower->GetStep(i)->GetPosi() + fShower->GetStep(i)->GetPosf())*0.5;
125        double dz = (fShower->GetStep(i)->GetPosi() - fShower->GetStep(i)->GetPosf()).Mag()/2/km;
126        double rmax = fShower->GetStep(i)->GetNumElectrons()*kScale;
127
128        if (fShower->GetStep(i)->GetNumElectrons()/fNmax > fThreshold) {
129            rmax = rmax < 0.001 ? 0.001 : rmax;
130            TString name = Form("tube_%d",i);
131            dz_shift = (stepCenter-initPos).Mag()/km-fDz; 
132            TGeoVolume *s = gGeoManager->MakeTube(name, fVacuum, 0.9*rmax, rmax, dz);
133            s->SetLineColor(50+(Int_t)(50.*log(fShower->GetStep(i)->GetNumElectrons()/fNmax)));
134            fShowerVolume->AddNode(s,1, new TGeoTranslation(0,0,-dz_shift));
135        }
136    }
137
138}
139
140//______________________________________________________________________________
141void EShowerPainter::BuildWorld() {
142    //
143    // Build World based on the shower geometry
144    //
145
146    if ( IsZombie() ) return;
147
148    if( IsWorldBuilt() ) return;
149
150    new TGeoManager("ShowerViewer","Shower Viewer");
151    TGeoMaterial *matVacuum = new TGeoMaterial("Al", 26.98,13,2.7);
152    fVacuum = new TGeoMedium("Al", 1, matVacuum);
153    fTOP = gGeoManager->MakeBox("fTOP",fVacuum,230,230,20);
154    gGeoManager->SetTopVolume(fTOP);   
155    Build();
156    Double_t OffSet = 1.1;
157    Double_t x1(-5),x2(5),y1(-5),y2(5); // make 10x10x30 kms box for default
158    if (fXmin < x1) x1 = fXmin*OffSet;
159    if (fXmax > x2) x2 = fXmax*OffSet;
160    if (fYmin < y1) y1 = fYmin*OffSet;
161    if (fYmax > y2) y2 = fYmax*OffSet;       
162    fFoVVolume = gGeoManager->MakeBox("fFoVVolume",fVacuum,(x2-x1)/2,(y2-y1)/2,15);
163    Double_t fTheta = fShower->GetTheta()*RadToDeg();
164    Double_t fPhi   = fShower->GetPhi()*RadToDeg();
165    TGeoTranslation tr;
166    TGeoRotation rot;
167    tr.SetTranslation(fCenter.X(),fCenter.Y(),fCenter.Z());
168    rot.SetAngles(fTheta,Min(fPhi,360-fPhi),0);
169    TGeoHMatrix h = tr*rot; 
170    fShowerVolume->SetLineColor(kGreen);
171    fFoVVolume->SetLineColor(kBlue);
172    fTOP->AddNode(fShowerVolume,1,new TGeoHMatrix(h));
173    fTOP->AddNode(fFoVVolume,1, new TGeoTranslation((x1+x2)/2,(y1+y2)/2,15));
174    gGeoManager->CloseGeometry();
175}
176
177//______________________________________________________________________________
178void EShowerPainter::UpdateVisibility( Option_t *opt ) {
179    //
180    //  Fast update ov fisibility status of the objects
181    //
182
183    if ( !gGeoManager ) return;
184
185    TGeoVolume *top;
186    if ( !(top = gGeoManager->GetTopVolume()) ) return;
187    TGeoVolume *shower = top->FindNode("fShowerVolume_1")->GetVolume();
188    if ( !shower )                              return;
189
190    Int_t n;
191    if ( fFrame == -1)
192       n = shower->GetNdaughters();
193    else
194       n = (Int_t)(((Float_t)shower->GetNdaughters()/fNumFrames)*fFrame);
195    shower->SetVisibility(kFALSE);
196    shower->VisibleDaughters(kFALSE);
197    for (Int_t i=0; i<shower->GetNdaughters(); i++) {
198        Bool_t vis = ( i <= n );
199        if ( shower->GetNode(i)->GetVolume()->IsVisible() != vis ) {
200            shower->GetNode(i)->GetVolume()->SetVisibility(vis);
201        }
202    }
203    shower->SetVisibility(kTRUE);   
204    shower->VisibleDaughters(kTRUE);
205    return;
206}
207
208//______________________________________________________________________________
209void EShowerPainter::Draw( Option_t* ) {
210    //
211    // Draw the shower in its world volume
212    //
213    BuildWorld();
214    UpdateVisibility();
215
216    if ( !IsZombie() ) {
217           fTOP->Draw();
218    }
219}
220
221//______________________________________________________________________________
222void EShowerPainter::DrawXY() {
223    // Draw the shower projection on XY Earth surface     
224   
225    // build histograms
226    TH2F *th0 = GetShowerHisto2D("ShowerXY", "XY projection");
227     for (Int_t i=0; i<fShower->GetNumSteps(); i++) {
228        double xmean = (fShower->GetStep(i)->GetPosi().X() + fShower->GetStep(i)->GetPosf().X())/2/km;
229        double ymean = (fShower->GetStep(i)->GetPosi().Y() + fShower->GetStep(i)->GetPosf().Y())/2/km;
230        double Ne    = fShower->GetStep(i)->GetNumElectrons();
231        th0->Fill(xmean,ymean,Ne);
232    }
233    th0->Draw("colz");
234}
235
236//______________________________________________________________________________
237void EShowerPainter::DrawXYZ() {
238    // Draw the shower track in XYZ
239   
240    // build histograms
241    TH3F *th0 = GetShowerHisto3D("ShowerXYZ", "XYZ development");
242     for (Int_t i=0; i<fShower->GetNumSteps(); i++) {
243        double xmean = (fShower->GetStep(i)->GetPosi().X() + fShower->GetStep(i)->GetPosf().X())/2/km;
244        double ymean = (fShower->GetStep(i)->GetPosi().Y() + fShower->GetStep(i)->GetPosf().Y())/2/km;
245        double zmean = (fShower->GetStep(i)->GetPosi().Z() + fShower->GetStep(i)->GetPosf().Z())/2/km;
246        double Ne    = fShower->GetStep(i)->GetNumElectrons();
247        th0->Fill(xmean,ymean,zmean,Ne/fNmax);
248    }
249    th0->Draw("box");
250}
251
252//_____________________________________________________________________________
253void EShowerPainter::Animate() {
254    //
255    // Move to next animation step and increase the counter
256    //
257   
258    if ( IsEnded() ) {
259        // leave the last frame to be painted
260        fFrame = fNumFrames-1;
261        Stop();
262        return;
263    } 
264
265    NextFrame();
266}
267
268//_____________________________________________________________________________
269void EShowerPainter::NextFrame() {
270    //
271    // Move to the next animation frame
272    //
273   
274    if ( !gGeoManager ) return;
275
276    TGeoVolume *top = gGeoManager->GetTopVolume();
277    TGeoVolume *shower = top->FindNode("fShowerVolume_1")->GetVolume();
278   
279    Int_t n = (Int_t)(((Float_t)shower->GetNdaughters()/fNumFrames)*fFrame);
280
281    for (Int_t i=0; i<shower->GetNdaughters(); i++) {
282        Bool_t vis = ( i <= n );
283        if ( shower->GetNode(i)->GetVolume()->IsVisible() != vis ) {
284             shower->GetNode(i)->GetVolume()->SetVisibility(vis);
285        }
286    }
287    fFrame++;
288}
289
290//_____________________________________________________________________________
291TH2F *EShowerPainter::GetShowerHisto2D(const char *name, const char *title) {
292    // Draw XY shower track projection. X and Y of the histogram have always center 
293    // at zero for a better display
294    //
295    DefineBox();
296
297    if ( title == 0 ) title = name;
298
299   
300    TH2F* th = (TH2F*)fHistoContainer->FindObject(name);
301   
302    if (!th) {
303        Double_t OffSet = 1.1;
304        Double_t x1(-5),x2(5),y1(-5),y2(5); // make 5 kms box for default
305        if (fXmin < x1) x1 = fXmin*OffSet;
306        if (fXmax > x2) x2 = fXmax*OffSet;
307        if (fYmin < y1) y1 = fYmin*OffSet;
308        if (fYmax > y2) y2 = fYmax*OffSet;
309        // define number of steps
310        Int_t nx = (Int_t)(x2-x1);
311        Int_t ny = (Int_t)(y2-y1);
312       
313        th = new TH2F( name, title,nx,x1,x2,ny,y1,y2);
314        fHistoContainer->Add(th);
315        th->SetStats(0);
316        th->SetXTitle("X [km]");
317        th->SetYTitle("Y [km]");
318    }
319
320    return th;
321}
322//_____________________________________________________________________________
323TH3F *EShowerPainter::GetShowerHisto3D(const char *name, const char *title) {
324    //
325    //
326    //
327    DefineBox();
328
329    if ( title == 0 ) title = name;
330       
331    TH3F* th = (TH3F*)fHistoContainer->FindObject(name);
332   
333    if (!th) {
334        Double_t OffSet = 1.1;
335        Double_t x1(-5),x2(5),y1(-5),y2(5),z1(0),z2(30); // make 10x10x30 kms box for default
336        if (fXmin < x1) x1 = fXmin*OffSet;
337        if (fXmax > x2) x2 = fXmax*OffSet;
338        if (fYmin < y1) y1 = fYmin*OffSet;
339        if (fYmax > y2) y2 = fYmax*OffSet;       
340        th = new TH3F( name, title, 30,x1,x2,30,y1,y2,30,z1,z2); 
341        fHistoContainer->Add(th);
342        th->SetStats(0);
343        th->SetXTitle("X [km]");
344        th->SetYTitle("Y [km]");
345        th->SetZTitle("Z [km]");
346    }
347
348    return th;
349}
Note: See TracBrowser for help on using the repository browser.