source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/detector/tools/src/OpticsResponseBuilder.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: 13.2 KB
Line 
1// $Id: OpticsResponseBuilder.cc 3005 2012-01-06 08:06:11Z fenu $
2// Author: ale   2005/01/19
3
4/*****************************************************************************
5 * ESAF: Euso Simulation and Analysis Framework                              *
6 *                                                                           *
7 *  Id: OpticsResponseBuilder                                                *
8 *  Package: <packagename>                                                   *
9 *  Coordinator: <coordinator>                                               *
10 *                                                                           *
11 *****************************************************************************/
12
13//_____________________________________________________________________________
14//
15// OpticsResponseBuilder
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 "OpticsResponseBuilder.hh"
27#include "ParentPhoton.hh"
28#include "Photon.hh"
29#include "EORSample.hh"
30#include "EORHeader.hh"
31#include "EsafRandom.hh"
32#include "OpticalSystem.hh"
33#include "IdealFocalSurface.hh"
34
35#include <TFile.h>
36#include <TStyle.h>
37#include <TH1.h>
38#include <TH2F.h>
39#include <TH3F.h>
40#include <TVector3.h>
41#include <TGraph.h>
42#include <TMultiGraph.h>
43#include <TCanvas.h>
44#include <TGaxis.h>
45#include <TLegend.h>
46#include <TRandom.h>
47#include <TROOT.h>
48#include <TMath.h>
49#include <TTree.h>
50#include <TBenchmark.h>
51#include <TClonesArray.h>
52#include "Ntrace_optics.hh"
53#include "NOpticalSystem.hh"
54
55using namespace sou;
56
57ClassImp(OpticsResponseBuilder)
58
59//_____________________________________________________________________________
60OpticsResponseBuilder::OpticsResponseBuilder() : fDBName(""), fDBFile(0),
61    fDBTree(0), fOpticalSystem(0), fIdealFocalSurface(0), fPsf(0), fPsfZoom(0),
62    fPsfRProj(0), fCentroid(0.,0.,0.), fPoints(0), fParent(0), fSample(0) {
63    //
64    // Constructor
65    //
66
67
68    // default parameters
69    fDiscRadius = ((NOpticalSystem*)fOpticalSystem)->GetOpticalSystemParam()->surface[1].r_lens;
70    fEusoRadius = 1325.*mm;
71    fRmax       = 30.;
72    fTriggDiam   = 5.;
73
74    // default settings
75    fThetaMin = 0*deg;
76    fThetaMax = 30*deg;
77    fLambdaMin = 300*nm;
78    fLambdaMax = 400*nm;
79
80    fBinsTheta = 30;
81    fBinsLambda = 100;
82    fPsfBins = 100;
83    fPsfZoomSide = 15.*mm;
84
85    fPoints = new TClonesArray("TVector3");
86    SetNrays(10000);
87
88    fParent   = new ParentPhoton;
89
90}
91
92//_____________________________________________________________________________
93OpticsResponseBuilder::~OpticsResponseBuilder() {
94    //
95    // Destructor
96    //
97
98    SafeDelete( fDBFile );
99    SafeDelete( fPsf );
100    SafeDelete( fPsfZoom );
101    SafeDelete( fPsfRProj );
102    SafeDelete( fPoints );
103    SafeDelete( fParent );
104    SafeDelete( fSample );
105}
106
107//_____________________________________________________________________________
108void OpticsResponseBuilder::SetNrays(Int_t n){
109    //
110    // Reset clear fPoints and reset its size
111    //
112   
113    fNrays = n;
114    fPoints->Clear();
115    fPoints->Expand(n);
116    fPoints->BypassStreamer(kFALSE);
117   
118}
119
120
121//______________________________________________________________________________
122void OpticsResponseBuilder::Clear() {
123    //
124    // Clears all temporary histograms
125    //
126
127    SafeDelete( fPsf );
128    SafeDelete( fPsfZoom );
129    SafeDelete( fPsfRProj );
130    fPoints->Clear();
131
132}
133
134//______________________________________________________________________________
135void OpticsResponseBuilder::OpenDB() {
136    //
137    // Helper function. Opens and initializes the database
138    //
139
140   
141    if ( fDBName.IsNull() ) {
142        MsgForm(EsafMsg::Warning, "Database file name is empty. File not opened.");
143        return;
144    }
145
146
147    fDBFile = new TFile(fDBName,"RECREATE");
148   
149    if ( fDBFile->IsZombie() ) {
150        MsgForm(EsafMsg::Warning, "File %s not opened.",fDBName.Data());
151        return;
152    }
153   
154    EORHeader *header = new EORHeader;
155
156    // instance database objects
157    fSample = new EORSample;
158
159    fDBTree = new TTree("ORDatabase","ORDatabase");
160    fDBTree->Branch("EORSample","EORSample",&fSample);
161
162    // fill the header
163    // simulation parameters
164    header->SetOpticsName(fOpticalSystem->IsA()->GetName());
165    header->SetIdealSurfaceName(fIdealFocalSurface->IsA()->GetName());
166
167    header->SetDiscRadius( fDiscRadius );
168    header->SetEusoRadius( fEusoRadius );
169    header->SetRmax( fRmax );
170    header->SetTriggDiam( fTriggDiam );
171   
172    // run parameters
173    header->SetBins(fBinsTheta,fBinsLambda);
174    header->SetNrays( fNrays );
175
176    // header will be deleted by TTree destructor
177    fDBTree->GetUserInfo()->Add(header);
178
179    MsgForm(EsafMsg::Info,"Database file %s opened.",fDBName.Data());
180}
181
182//______________________________________________________________________________
183void OpticsResponseBuilder::CloseDB() {
184    //
185    //
186    //
187   
188    if ( !fDBFile ) {
189        MsgForm(EsafMsg::Warning,"No database file opened.");
190        return;
191    }
192   
193    if ( fDBFile != fDBTree->GetCurrentFile() ) {
194        MsgForm(EsafMsg::Warning,"Mismatch between database file and current TTree file.");
195        Msg(EsafMsg::Warning) << "Maybe a file change has occurred. ";
196        Msg(EsafMsg::Warning) << "Only the latter will be saved." << MsgDispatch;
197       
198        fDBFile = fDBTree->GetCurrentFile(); 
199    }
200   
201    fDBTree->Write();
202    SafeDelete( fDBTree );
203    SafeDelete( fDBFile );
204    SafeDelete( fSample );
205}
206
207//_____________________________________________________________________________
208Bool_t OpticsResponseBuilder::ComputeEfficacy(Float_t theta, Float_t phi, Float_t lambda) {
209    //
210    // Simulate and collect a uniform beam of photons for the chosen optical
211    // system.
212    //
213
214
215    if ( !fOpticalSystem ) {
216        Msg(EsafMsg::Warning) << "No OpticalSystem defined. Exiting."<< MsgDispatch;
217        return kFALSE;
218    }
219
220    if ( !fIdealFocalSurface ) {
221        Msg(EsafMsg::Warning) << "No IdealFocalSurface defined. Exiting."<< MsgDispatch;
222        return kFALSE;
223    }
224
225   
226    Double_t pos[3] = {0,0,0};
227    Double_t spos[3] = {0,0,0};
228
229    Int_t ntrans(0);
230    Photon* ph;
231    TVector3 dir(TMath::Sin(theta)*TMath::Cos(phi), TMath::Sin(theta)*TMath::Sin(phi),TMath::Cos(theta));
232    for( Int_t i(0); i<fNrays; i++) {
233
234        // new photon position
235        Double_t x,y,z;
236        RandomPosOnLense(x,y,z);
237        TVector3 pos_back(x,y,z);
238//      pos_back.Dump();
239        pos_back = pos_back - 1*dir;
240//      pos_back.Dump();
241        pos[0] = pos_back.X();
242        pos[1] = pos_back.Y();
243        pos[2] = pos_back.Z();
244       
245        // update the fParent
246        fParent->Build(0, spos, pos, theta, phi, lambda, 0, Fluorescence);
247        ph = &(fParent->GetPhoton());
248
249        ph = fOpticalSystem->Transport(ph);
250//         printf("optics history %d\n",ph->history);
251        // reject dead photons and going downward
252        if (!ph || ph->dir.Z() < 0.) continue;
253
254        fIdealFocalSurface->HitPosition(ph);
255        if (!ph->hitIfs) continue;
256        TClonesArray &ref = *fPoints;
257        new(ref[ntrans++]) TVector3(ph->posOnIfs);
258    }
259   
260    Float_t pflux = fNrays/(TMath::Pi()*TMath::Power(fDiscRadius,2)*TMath::Cos(theta));//wrong normalization for SIDECUT!
261    fTotalEfficacy = ntrans/pflux;
262    if ( fSample ) {
263        // fill the sample with the current results
264        fSample->SetTheta( theta );
265        fSample->SetLambda( lambda );
266        fSample->SetTotalEfficacy( fTotalEfficacy );
267    }
268   
269    FillPsfHist();
270
271    if ( !fPsfRProj ) {
272        MsgForm(EsafMsg::Warning,"No photons in the point spread function");
273        return kTRUE;
274    }
275
276    Int_t bin = fPsfRProj->GetXaxis()->FindBin(fTriggDiam/2.);
277    fTriggEfficacy = fPsfRProj->Integral(1,bin)/pflux;
278    fPsf->Scale(1./pflux);
279    fPsfZoom->Scale(1./pflux);
280   
281    if ( fSample ) {
282        fSample->SetTriggEfficacy( fTriggEfficacy );
283        fSample->SetCentroid( fCentroid );
284        fSample->SetPsfHist( fPsfZoom );
285    }
286
287    return kTRUE;
288}
289
290//______________________________________________________________________________
291void OpticsResponseBuilder::RandomPosOnLense(Double_t &x, Double_t &y, Double_t &z) {
292    //This method returns a random position on a bottom lense of EUSO.
293    TRandom* rndm = EsafRandom::Get();
294    Double_t LenseCurvature = ((NOpticalSystem*)fOpticalSystem)->GetOpticalSystemParam()->surface[1].Rc;
295    Double_t Ycut = ((NOpticalSystem*)fOpticalSystem)->GetOpticalSystemParam()->surface[1].r_cut;
296    Double_t fDiskRadius = ((NOpticalSystem*)fOpticalSystem)->GetOpticalSystemParam()->surface[1].r_lens;
297//     Double_t zmax = fOpticalSystem->FirstLensTop();
298    Double_t zmin = fOpticalSystem->Bottom();
299    Double_t CosThetaMin = TMath::Sqrt(pow(LenseCurvature,2) - pow(fDiskRadius,2))/LenseCurvature;
300   
301    while (1) {
302      Double_t cosTheta = CosThetaMin + rndm->Rndm()*(1-CosThetaMin);
303      Double_t phi = TMath::TwoPi()*rndm->Rndm();
304      Double_t sinTheta = TMath::Sqrt(1-pow(cosTheta,2));
305      y = LenseCurvature*sinTheta*TMath::Sin(phi);
306      if (fabs(y) < Ycut) {
307        x = LenseCurvature*sinTheta*TMath::Cos(phi);
308        z = LenseCurvature*(1-cosTheta) + zmin;
309//      printf("%g %g %g %g %g %g\n",x,y,z,zmin,zmax,fDiskRadius);
310        break;
311      }
312    }
313}
314//______________________________________________________________________________
315void OpticsResponseBuilder::FillPsfHist() {
316    //
317    // Build histograms
318    //
319
320    Float_t xmin, ymin;
321    Float_t xmax, ymax; 
322    xmin = ymin =-fEusoRadius;
323    xmax = ymax = fEusoRadius;
324    Int_t npos = fPoints->GetEntriesFast();
325
326    if ( fPsf ) {
327        //TOCHECK: is it fast?
328        fPsf->SetBins(20*fPsfBins, xmin, xmax, 20*fPsfBins, ymin, ymax);
329        fPsf->Reset();
330     } else {   
331        fPsf = new TH2F("psf", "Psf", 20*fPsfBins, xmin, xmax, 20*fPsfBins, ymin, ymax);
332        fPsf->SetDirectory(0);
333    }
334
335    // search for the peak
336    for (Int_t i=0; i< npos; i++) {
337        TVector3 *pos = (TVector3*)fPoints->At(i);
338        fPsf->Fill(pos->X(),pos->Y());
339    }
340    Int_t bx, by, bz;
341    fPsf->GetMaximumBin(bx, by, bz);
342
343
344    Double_t peakX = fPsf->GetXaxis()->GetBinCenter( bx );
345    Double_t peakY = fPsf->GetYaxis()->GetBinCenter( by );
346
347    TVector3 fAverage(0,0,0);
348
349    Int_t nnear(0);
350    Float_t spacethreshold = 20.*mm;
351
352    // find the center of mass around the peak
353    for (Int_t i=0; i< npos; i++ ){
354        TVector3 *pos = (TVector3*)fPoints->At(i);   
355        Double_t dist = TMath::Sqrt(TMath::Power(pos->X()-peakX,2)
356                +TMath::Power(pos->Y()-peakY,2)); 
357        if ( dist > spacethreshold) continue;
358        fAverage  += *pos;
359        nnear++;
360    }
361
362    if (nnear == 0) 
363        return;
364
365    fAverage *= 1./nnear;
366
367    xmin = fAverage.X()-fPsfZoomSide;
368    xmax = fAverage.X()+fPsfZoomSide;
369    ymin = fAverage.Y()-fPsfZoomSide;
370    ymax = fAverage.Y()+fPsfZoomSide;
371
372    // build the projection
373    if ( fPsfRProj ) { 
374        fPsfRProj->SetBins( fPsfBins, 0, fPsfZoomSide );
375        fPsfRProj->Reset();
376    } else {   
377        fPsfRProj = new TH1F("PsfRProj", "Psf R-projection", fPsfBins, 0, fPsfZoomSide); 
378        fPsfRProj->SetDirectory(0);
379    }
380
381    // build the zoom
382    if ( fPsfZoom ) {
383        fPsfZoom->SetBins( fPsfBins, xmin, xmax, fPsfBins, ymin, ymax);
384        fPsfZoom->Reset();
385    } else {
386        fPsfZoom = new TH2F("PsfZoom", "Psf Zoom", 3*fPsfBins, xmin, xmax, 3*fPsfBins, ymin, ymax);
387        fPsfZoom->SetDirectory(0);
388    }
389       
390
391    fCentroid = fAverage;
392
393    Float_t r;
394    for (Int_t i=0; i< nnear; i++) {
395        TVector3 *pos = (TVector3*)fPoints->At(i);   
396        fPsfZoom->Fill(pos->X(), pos->Y());
397        r = (*pos - fAverage).Mag();
398        if( r < fRmax) 
399            fPsfRProj->Fill( r );
400    }
401
402    fPsf->GetXaxis()->SetTitle("mm");
403    fPsf->GetYaxis()->SetTitle("mm");
404    fPsf->GetYaxis()->SetTitleOffset(1.2);
405
406    fPsfRProj->GetXaxis()->SetTitle("mm");
407
408}
409
410
411//______________________________________________________________________________
412void OpticsResponseBuilder::BuildEfficacyTables() {
413    //
414    //
415    //
416   
417    OpenDB();
418
419    Float_t th, wl;
420    Float_t thstep = (fThetaMax-fThetaMin)/fBinsTheta;
421    Float_t wlstep = (fLambdaMax-fLambdaMin)/fBinsLambda;
422
423    if ( fDBTree ) {
424        EORHeader* header = (EORHeader*)fDBTree->GetUserInfo()->FindObject("EORHeader");
425        if ( header ) {
426            header->SetThetaMin(fThetaMin);
427            header->SetThetaMax(fThetaMax);
428            header->SetThetaIndexMin(0);
429            header->SetThetaIndexMax(fBinsTheta);
430            header->SetLambdaMin(fLambdaMin);
431            header->SetLambdaMax(fLambdaMax);
432            header->SetLambdaIndexMin(0);
433            header->SetLambdaIndexMax(fBinsLambda); 
434            header->SetPsfBins(fPsfBins);       
435        }
436    }
437
438    for (Int_t i(0); i <=fBinsTheta; i++) {
439        th = fThetaMin + thstep*i; 
440        Msg(EsafMsg::Info) << Form("Theta: %.2f deg ", th/deg) << MsgFlush;
441        for (Int_t j(0); j <=fBinsLambda; j++) {
442            wl = fLambdaMin + wlstep*j; 
443            if ( (j+1)*10%fBinsLambda == 0 )
444            Msg(EsafMsg::Info) << " ... "<<j*100/fBinsLambda << '%'<< MsgFlush;   
445
446            ComputeEfficacy(th, 0. ,wl);
447
448            if ( fSample ) {
449                fSample->SetThetaIndex( i );
450                fSample->SetLambdaIndex( j );
451            }
452
453            if ( fDBTree ) fDBTree->Fill();
454
455            if ( fSample ) fSample->Clear();
456            Clear();
457        }
458        Msg(EsafMsg::Info) << " -done. " << MsgDispatch; 
459    }
460    CloseDB();
461}
Note: See TracBrowser for help on using the repository browser.