source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/detector/tools/src/PixelMapBuilder.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: 29.4 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: PixelMapBuilder.cc 2927 2011-06-16 18:02:41Z mabl $
3// A.Thea created Apr, 15 2004
4
5#include "PixelMapBuilder.hh"
6#include "EusoDetector.hh"
7#include "DetectorTransportManager.hh"
8#include "EusoElectronics.hh"
9//#include "OAPxBunch.hh"
10#include "Photon.hh"
11#include "ParentPhoton.hh"
12#include "DetectorGeometry.hh"
13#include <TTree.h>
14#include <TChain.h>
15#include "TTimeStamp.h"
16#include "TString.h"
17#include <TCanvas.h>
18#include <TF1.h>
19#include <TText.h>
20#include <TEllipse.h>
21
22#include <TPolyMarker.h>
23#include <TSystem.h>
24#include <TGClient.h>
25#include <TGProgressBar.h>
26#include <TMutex.h>
27#include "EGViewer.hh"
28#include <TMath.h>
29
30using namespace TMath;
31
32//______________________________________________________________________________
33//
34// PixelMapBuilder
35//
36// Works in 2 modes
37//
38// Create the photon's file for the anaysis
39//
40// fNphotons
41// fDetector
42// fTransporter
43// fElectronics
44// fHitsRootFile
45// fHitsTree
46// fFoV
47// fHit
48// fHeader
49//
50//
51
52ClassImp(PixelMapBuilder)
53
54//______________________________________________________________________________
55PixelMapBuilder::PixelMapBuilder() : EsafConfigurable(), fHitsRootFile(0), fHitsTree(0), fHitsChain(0),
56    fFitHistsFile(0), fFitHistsTree(0), fFitUID(0), fH1FitTheta(0),
57    fH1FitPhi(0), fH2FitThetaPhi(0), fH2ThetaPhi(0), fDisplay(0) {
58    //
59    // ctor
60
61    fMaxBufSize = 5000000/(sizeof(Float_t)*2);
62    fBufferSize = 250 /*Mb*/ *1024*1024;
63    fNsigma       = 1.5;
64
65    fHitsThreshold = 0;
66    fNbins        = 100;
67
68    fBufIndex = kBufEmpty;
69
70    //    fFullFoV = Pi()/5.;
71    //    fFoV = Pi()/5.;
72    fFullFoV = Conf()->GetNum("DiffusePhotonsOnPupil.fThetaFOVMax")*TMath::DegToRad();
73    fFoV = fFullFoV;
74//    fBufferSize = 1;
75}
76
77//______________________________________________________________________________
78PixelMapBuilder::~PixelMapBuilder() {
79    //
80    // Destructor
81    //
82}
83
84
85//______________________________________________________________________________
86Bool_t PixelMapBuilder::AddToChain( const char* name ) {
87
88    if ( fHitsChain ) fHitsChain->Add(name,0);
89    return kTRUE;
90}
91
92//______________________________________________________________________________
93Bool_t PixelMapBuilder::OpenChain( const char* name ) {
94
95    fHitsChain = new TChain("oapxtree");
96
97    AddToChain(name);
98    fHitsChain->GetEntries();
99
100    fHitsChain->SetBranchAddress("Hit", &fHit);
101    fHitsChain->GetEntry();
102    fHeader = (Header*)fHitsChain->GetTree()->GetUserInfo()->FindObject("PixelMapBuilder::Header")->Clone();
103
104    fHitsTree = fHitsChain;
105    return kTRUE;
106}
107
108//______________________________________________________________________________
109Bool_t PixelMapBuilder::OpenChain( const char** name, Int_t n ) {
110
111    fHitsChain = new TChain("oapxtree");
112
113    for(Int_t i(0); i<n; i++) {
114        AddToChain(name[i]);
115    }
116    fHitsChain->GetEntries();
117
118    fHitsChain->SetBranchAddress("Hit", &fHit);
119    fHitsChain->GetEntry();
120    fHeader = (Header*)fHitsChain->GetTree()->GetUserInfo()->FindObject("PixelMapBuilder::Header")->Clone();
121
122    fHitsTree = fHitsChain;
123    return kTRUE;
124}
125
126//______________________________________________________________________________
127Bool_t PixelMapBuilder::OpenRoot( const char* name, Bool_t write ) {
128    //
129    // Open rootfile with the correct name
130    //
131    CloseRoot();
132
133    TString fname(name);
134    TString ext = ".root";
135
136    if ( write) {
137        if ( fname.EndsWith(ext)) fname.Resize(fname.Length()-ext.Length());
138        Msg(EsafMsg::Info) << "Opening root file " << fname << ext << MsgDispatch;
139        fHitsRootFileName = fname+ext;
140        fHitsRootFile = new TFile(fHitsRootFileName,"CREATE");
141
142        if (fHitsRootFile->IsZombie()) {
143            Msg(EsafMsg::Warning) << MsgDispatch;
144            Msg(EsafMsg::Warning) << "Error opening file " << fname << MsgDispatch;
145            Msg(EsafMsg::Warning) << "Probably it already exists or there is no " << MsgDispatch;
146            Msg(EsafMsg::Warning) << "write access permission to this directory" << MsgDispatch;
147            Msg(EsafMsg::Warning) << MsgDispatch;
148
149            // try adding date and time to the name
150            TTimeStamp time;
151            TString now = time.AsString("s");
152            now[now.Last(' ')]=':';
153            TString newFileName = Form("%s-", fname.Data())+now+ext;
154            fHitsRootFile = new TFile(newFileName,"CREATE");
155            if (fHitsRootFile->IsZombie()) {
156                MsgForm(EsafMsg::Warning,"Cannot open root file");
157                return kFALSE;
158            }
159        }
160        Msg(EsafMsg::Info) << "Root file " << fHitsRootFile->GetName() << " successfully opened\n" << MsgDispatch;
161
162        fHitsTree = new TTree("oapxtree","Angle Pixel Map photons' database");
163        fHitsTree->Branch("Hit",&fHit,"fUID/I:fTheta/F:fPhi/F:fLambda/F");
164        fHeader = new Header;
165        fHitsTree->GetUserInfo()->Add(fHeader);
166
167        return kTRUE;
168    } else {
169        if ( fname.EndsWith(ext)) fname.Resize(fname.Length()-ext.Length());
170        fHitsRootFile = new TFile((fname+ext));
171
172        if ( fHitsRootFile->IsZombie() ) return kFALSE;
173
174        fHitsTree = (TTree*)fHitsRootFile->Get("oapxtree");
175        fHitsTree->SetBranchAddress("Hit", &fHit);
176
177        fHeader = (Header*)fHitsTree->GetUserInfo()->FindObject("PixelMapBuilder::Header");
178        if ( !fHeader ) {
179            Msg(EsafMsg::Warning) << "No header in the rootfile" << MsgDispatch;
180            return kFALSE;
181        }
182
183        return kTRUE;
184    }
185}
186
187//______________________________________________________________________________
188Bool_t PixelMapBuilder::CloseRoot() {
189    //
190    // close root file
191    //
192
193    if ( fHitsRootFile && fHitsChain ) {
194        Msg(EsafMsg::Warning) << "Both fHitsRootFile and fHitsChain are present. Something's wronmg. Anorting" << MsgDispatch;
195        return kFALSE;
196    }
197
198    if ( fHitsRootFile ) {
199        if ( fHitsRootFile->IsWritable() ) {
200            fHitsRootFile->Write();
201
202            SafeDelete(fHitsTree);
203            fHeader = 0;
204
205            Msg(EsafMsg::Info) << "Root file saved" << MsgDispatch;
206        } else {
207            SafeDelete(fHitsTree);
208            fHeader = 0;
209        }
210        fHitsRootFile->Close();
211
212        Msg(EsafMsg::Info) << "Root file closed" << MsgDispatch;
213        SafeDelete( fHitsRootFile )
214    } else if (fHitsChain) {
215        // chain, reading by default
216        SafeDelete(fHitsChain);
217        SafeDelete(fHeader);
218        fHitsTree = 0;
219    }
220
221    return kTRUE;
222}
223
224
225//______________________________________________________________________________
226void PixelMapBuilder::TrackPhotons() {
227    //
228    // build the map
229    //
230
231    // disable electronics simulation, just find the photon's pad
232    fElectronics->EnableSimulation(kFALSE);
233    fHit.fUID = 0;
234    fHit.fTheta = 0;
235    fHit.fPhi = 0;
236    fHit.fLambda = 0;
237    fPupil.SetThetaFOVMax(fFoV);
238    fHeader->fFoV = fFoV;
239    fHeader->fFullFoV = fDetector->GetGeometry()->GetFoV();
240
241    Msg(EsafMsg::Info) << "Photons to track: " << fNphotons << MsgDispatch;
242    Msg(EsafMsg::Info) << "PixelMapBuilder: tracing..." << MsgDispatch;
243    Long_t hits(0);
244    EVector initdir, initpos;
245    EVector calcdir,dirdiff;
246    for(Long_t i(0); i<fNphotons; i++) {
247        Photon *p = fPupil.GetPhoton();
248        initdir=p->dir;
249        initpos=p->pos;
250        fTransporter->Transport(p);
251        calcdir=(p->pos-initpos).Unit();
252        dirdiff=calcdir-initdir;
253//        double dirdiff_0=(p->dir-initdir).Mag();
254//        if ( p->hitIfs && dirdiff_0<1.e-9 ){
255//            double theta=initdir.Theta();
256//            cout<<"strang photon hits ifs"<<endl;
257//            cout<<"initpos="<<initpos<<" rho="<<initpos.Perp()<<endl;
258//            cout<<"initdir="<<initdir<<" theta="<<theta<<endl;
259//            cout<<"fspos="<<p->posOnIfs<<" rho="<<p->posOnIfs.Perp()<<endl;
260//            cout<<"latest pos="<<p->pos<<" rho="<<p->pos.Perp()<<endl;
261//            cout<<"latest dir="<<p->dir<<" theta="<<p->dir.Theta()<<endl;
262//        }
263
264        fHeader->fNtracked++;
265
266        if ( p->pixelUid > 0 ) { // FIXME: should be > 0
267            TVector3 dir = p->parent->Dir();
268            hits++;
269            fHeader->fNhits++;
270
271            fHit.fUID =  p->pixelUid;
272            fHit.fTheta = dir.Theta();
273            fHit.fPhi = dir.Phi();
274            fHit.fLambda = p->wl;
275
276            fHitsTree->Fill();
277            if ( i%1000000 == 0 ) fHitsTree->AutoSave();
278
279            if ( fHitsTree->GetCurrentFile() != fHitsRootFile ) {
280                fHitsRootFile = fHitsTree->GetCurrentFile();
281                fHeader->fNtracked = 0;
282                fHeader->fNhits = 0;
283            }
284
285        }
286
287        if ( (i+1)%1000 == 0) {
288            Msg(EsafMsg::Info) << '\r' <<  MsgFlush;
289            Msg(EsafMsg::Info) << Form("%ld photons traced, %ld hits on the pixels",
290                    (i+1), hits) << MsgFlush;
291        }
292    }
293    Msg(EsafMsg::Info) << MsgDispatch;
294    Msg(EsafMsg::Info) << "Completed" << MsgDispatch;
295
296    MsgForm(EsafMsg::Info, "%d photons traced.",fNphotons);
297    MsgForm(EsafMsg::Info, "%d photons hit the pixels.",hits);
298    // FIXME: should be false if it was before TrackPhotons()
299    fElectronics->EnableSimulation(kTRUE);
300}
301
302
303//______________________________________________________________________________
304Bool_t PixelMapBuilder::MakePhotonsFile( const char* name, Long_t n) {
305
306
307    fNphotons = n;
308
309
310    Msg(EsafMsg::Info) << "Building Detector" << MsgDispatch;
311
312    fDetector = new EusoDetector();
313    fTransporter = fDetector->GetDetectorTransportManager();
314    //fElectronics = fDetector->GetEusoElectronics();
315    fElectronics = gEusoElectronics;
316
317    Msg(EsafMsg::Debug) << "TransportManager Modules" << MsgDispatch;
318    Msg(EsafMsg::Debug) << "Optics: " << fTransporter->GetOptics()->IsA()->GetName() << MsgDispatch;
319    Msg(EsafMsg::Debug) << "FocalPlane: " << fTransporter->GetFocalPlane()->IsA()->GetName() << MsgDispatch;
320    Msg(EsafMsg::Info) <<  "Building complete" << MsgDispatch;
321
322    if ( !OpenRoot(name, kTRUE) ) {
323        Msg(EsafMsg::Warning) << "Unable to open the root file. Aborting" << MsgDispatch;
324        return kFALSE;
325    }
326
327    // setup the tree
328
329    fHeader->fNpixels = fElectronics->NumOfCh();
330
331    // loop here
332    TrackPhotons();
333
334    CloseRoot();
335    SafeDelete(fHeader);
336
337    return kTRUE;
338}
339
340//______________________________________________________________________________
341Bool_t PixelMapBuilder::MakeMap( const char* hits, const char* mapfile) {
342
343    if ( !OpenRoot(hits) ) return kFALSE;
344    InitFit();
345
346//    Int_t first = 1;
347//    Int_t last = fHeader->fNpixels;
348    ComputeMap( 1, 0, "WEQ0" );
349    ClearFit();
350
351    CloseRoot();
352
353    WriteMap(mapfile);
354
355    return kTRUE;
356}
357
358//______________________________________________________________________________
359void PixelMapBuilder::ComputeMap(Int_t first, Int_t last, const char* opt) {
360    //
361    // Fits the histograms
362    //
363    TString options(opt);
364
365    if ( !fMapTheta.GetSize() ) {
366        Msg(EsafMsg::Warning) << "Map not defined. Aborting" << MsgDispatch;
367        return;
368    }
369    if ( !first ) first = 1;
370    if ( !last ) last = fHeader->fNpixels;
371
372    Msg(EsafMsg::Info) << "Processing map, pixels " << first <<
373                          " - " << last << MsgDispatch;
374    TF1 *g;
375    Double_t parTh[3], parPhi[3];
376
377    //
378    // Tentative display
379    //
380    Bool_t display = gClient && fDisplay;
381    TCanvas *cPoly = NULL, *cHist = NULL;
382    TPolyMarker* m = NULL;
383    TGProgressBar* p = NULL;
384    TH2F* h = NULL;
385    Stat_t entries(0);
386
387    if ( display ) {
388        //
389        // Polymarker display
390        //
391        cPoly = fDisplay->FindAddTab("PixelMapDisplayPoly","Map projection");
392        cPoly->cd();
393        cPoly->DrawFrame(-Sin(fHeader->fFoV),-Sin(fHeader->fFoV),
394                      Sin(fHeader->fFoV), Sin(fHeader->fFoV));
395        m = new TPolyMarker(fHeader->fNpixels);
396//        m->SetMarkerStyle(7);
397        m->SetBit(kCanDelete);
398        m->Draw();
399
400        //
401        // Histogram display
402        //
403        cHist = fDisplay->FindAddTab("PixelMapDisplayHist","Histogram");
404        cHist->cd();
405        h = new TH2F("PixelMap","PixelMap",
406                           100, -Sin(fHeader->fFoV),Sin(fHeader->fFoV),
407                           100, -Sin(fHeader->fFoV),Sin(fHeader->fFoV));
408        h->SetDirectory(0);
409        h->SetBit(kCanDelete);
410        h->SetXTitle("Sin(#theta) #times Cos(#phi)");
411        h->SetYTitle("Sin(#theta) #times Sin(#phi)");
412        h->Draw("col");
413
414        p = fDisplay->GetProgressBar();
415        p->SetRange(1,fHeader->fNpixels);
416        fDisplay->ShowProgressBarPosition(kTRUE,kFALSE,"%.0f pixels");
417
418        // no fit display if the display is active.
419        if ( !options.Contains("0")) options.Append("0");
420    }
421
422
423    for(Int_t k(first); k<=last; k++){
424        //
425        // fit mut be called before UidToIndex because the buffer
426        // may be not be filled yet
427        //
428        Bool_t fit = Fit(k, options);
429        //cout << "fit " << fit << endl;
430        Int_t j=UidToIndex(k);
431   //cout << "index " << j << endl;
432
433        //
434        // fit pixel data
435        //
436        if (!fit) {
437            fMapTheta[k-1]    = -1; // -2*TwoPi()
438            fMapThetaRMS[k-1] = -1; // 0
439            fMapPhi[k-1]      = -1; // -2*TwoPi()
440            fMapPhiRMS[k-1]   = -1; // 0
441        } else {
442            // retrieve results from histos
443            g = fH1FitTheta->GetFunction("gaus");
444
445            if ( g ) g->GetParameters(parTh);
446            // if g doesn't exists put dummy values
447            else  parTh[0] = parTh[1] = parTh[2] = -1;
448
449            g = fH1FitPhi->GetFunction("gaus");
450            if ( g ) g->GetParameters(parPhi);
451            else  parPhi[0] = parPhi[1] = parPhi[2] = -1;
452
453            // fill the containers
454            fMapTheta[k-1]    = parTh[1];
455            fMapThetaRMS[k-1] = parTh[2];
456            fMapPhi[k-1]      = parPhi[1];
457            fMapPhiRMS[k-1]   = parPhi[2];
458        }
459        if ( fFitHistsTree ) fFitHistsTree->Fill();
460
461        if ( display && fit ) {
462
463            Float_t x = Sin(fMapTheta[k])*Cos(fMapPhi[k]);
464            Float_t y = Sin(fMapTheta[k])*Sin(fMapPhi[k]);
465            m->SetPoint(j,x,y);
466            h->Fill(x,y);
467
468        }
469
470        if ( !(k%100) || k == last ) {
471            cout  << '\r' << Form("Pixel %d of %d completed (%d%%) ",
472                  k, fHeader->fNpixels, k*100/fHeader->fNpixels) << flush;
473
474            //
475            // update the display
476            //
477            if (display && ( h->GetEntries() != entries || k == last ) ) {
478                cPoly->Modified();
479                cPoly->Update();
480                cHist->Modified();
481                cHist->Update();
482                p->SetPosition((Float_t)k);
483                if (!TThread::Self())
484                    gSystem->ProcessEvents();
485
486                entries = h->GetEntries();
487            }
488        }
489
490     }
491     cout << endl;
492
493
494
495}
496
497//______________________________________________________________________________
498void PixelMapBuilder::OptimizeBuffer() {
499    //
500    // Tries to evaluate the expectes average number of hits per pixel and the
501    // length of the buffers to match fBufferSize. Then reserves the
502    // corresponding amount of memory in the bugger
503    //
504
505    if (!fHitsTree) return;
506
507    Long64_t entries = fHitsTree->GetEntries();
508
509    //
510    // expected average number of hits per pixel
511    //
512    Int_t average = Nint( entries*(1-Cos(fFullFoV))/(fHeader->fNpixels*(1-Cos(fHeader->fFoV))) );
513    average=average?average:1;
514
515    Msg(EsafMsg::Info) << "Expected average number of entries per pixel: " << average << MsgDispatch;
516
517    //
518    // fBufferSize = sizeof(float)*2*average*length
519    //
520    double bufsize=sizeof(float)*2*average;
521    Int_t length = fBufferSize/int(bufsize);
522    length = Min(length,fHeader->fNpixels);
523    Msg(EsafMsg::Info) << "Buffer lenght set to " << length << MsgDispatch;
524    Msg(EsafMsg::Info) << "Expected buffer size in memory " <<
525                          (sizeof(float)*2*average*length)/(1024*1024) << " Mb"<< MsgDispatch;
526
527    // prepare the buffers
528    fThetaFOVBuffer.resize(length);
529    fPhiFOVBuffer.resize(length);
530
531    for(size_t k(0); k<fThetaFOVBuffer.size(); k++) {
532        fThetaFOVBuffer[k].reserve(average);
533        fPhiFOVBuffer[k].reserve(average);
534    }
535
536}
537
538//______________________________________________________________________________
539Bool_t PixelMapBuilder::InitFit() {
540
541    if ( !fHitsTree ) return kFALSE;
542
543    // clear the content of the buffers.
544    ClearPixelsBuffers();
545    // resize the map arrays
546    fMapTheta.Set(fHeader->fNpixels);
547    fMapThetaRMS.Set(fHeader->fNpixels);
548    fMapPhi.Set(fHeader->fNpixels);
549    fMapPhiRMS.Set(fHeader->fNpixels);
550
551    fMapIndex.Set(fHeader->fNpixels+1);
552    fMapIndex[0] = 0;
553
554    OptimizeBuffer();
555
556    // note: upper limits in theta is a little bit more than Pi/6
557    Bool_t status = TH1::AddDirectoryStatus();
558    TH1::AddDirectory(kFALSE);
559
560    fH1FitTheta  = new TH1F("oapx_thfit","Temp fit histogram for Theta",fNbins, 0, 1);
561    fH1FitTheta->SetXTitle("theta");
562    fH1FitTheta->SetYTitle("hits");
563    fH1FitPhi    = new TH1F("oapx_phifit","Temp fit histogram for Phi", fNbins, 0, 1);
564    fH1FitPhi->SetXTitle("phi");
565    fH1FitPhi->SetYTitle("hits");
566    fH2FitThetaPhi = new TH2F("oapx_thphifit","Temp fit histogram for Theta and Phi",
567                             fNbins, 0, 1, fNbins, 0, 1);
568    fH2FitThetaPhi->SetXTitle("theta");
569    fH2FitThetaPhi->SetYTitle("phi");
570
571    Float_t r = Sin(fHeader->fFoV*1.1);
572    fH2ThetaPhi =  new TH2F("oapx_tp","tp distribution projected on a plane",
573                                     fNbins, -r, r, fNbins, -r, r);
574    fH2FitThetaPhi->SetXTitle("Sin(#theta) #times Cos(#phi)");
575    fH2FitThetaPhi->SetYTitle("Sin(#theta) #times Sin(#phi)");
576
577    TH1::AddDirectory(status);
578
579    if ( !fFitHistsFileName.IsNull() ) {
580        TString ext = ".root";
581        TString fname = fFitHistsFileName;
582
583        if ( fname.EndsWith(ext)) fname.Resize(fname.Length()-ext.Length());
584        Msg(EsafMsg::Info) << "Saving histograms in " << fname << ext << MsgDispatch;
585        fFitHistsFileName = fname+ext;
586        fFitHistsFile = new TFile(fFitHistsFileName,"RECREATE");
587
588        fFitHistsTree = new TTree("pxmap_hist","Pixel Map histograms");
589        fFitHistsTree->Branch("fTheta",&fH1FitTheta, 100000, 0);
590        fFitHistsTree->Branch("fPhi",&fH1FitPhi, 100000, 0);
591        fFitHistsTree->Branch("fThetaPhi",&fH2FitThetaPhi, 100000, 0);
592        fFitHistsTree->Branch("fThetaPhiWholeFoV",&fH2ThetaPhi, 100000, 0);
593        fFitHistsTree->Branch("fUID",&fFitUID, "fUID/I");
594
595    }
596
597    fFitUID = 0;
598
599    return kTRUE;
600}
601
602//______________________________________________________________________________
603Bool_t PixelMapBuilder::ClearFit() {
604    //
605    // Clears the histograms usend in the fit procedure
606    //
607
608    if ( fFitHistsTree ) {
609        fFitHistsFile = fFitHistsTree->GetCurrentFile();
610        fFitHistsFile->cd();
611        fFitHistsTree->Write();
612        SafeDelete(fFitHistsTree);
613        SafeDelete(fFitHistsFile);
614    }
615
616    fMapIndex.Set(0);
617//    fHeader = 0;
618
619    SafeDelete(fH1FitTheta);
620    SafeDelete(fH1FitPhi);
621    SafeDelete(fH2FitThetaPhi);
622    SafeDelete(fH2ThetaPhi);
623
624    ClearPixelsBuffers();
625
626    return kTRUE;
627}
628
629//______________________________________________________________________________
630Bool_t PixelMapBuilder::FillPixelsBuffer(Int_t startuid ) {
631    //
632    // Fills the buffer starting from
633    //
634
635    //
636    // Check boundaries
637    //
638    if ( startuid >  fHeader->fNpixels ){
639        Msg(EsafMsg::Warning) << "First out of bounds" << endl;
640        return kFALSE;
641    }
642
643
644    Int_t length = BufferEntries();
645    Int_t lastuid = (startuid+length-1) <= fHeader->fNpixels ?
646                 (startuid+length-1) : (fHeader->fNpixels-1);
647
648    Long64_t nhits = fHitsTree->GetEntriesFast();
649    //
650    // clear the buffer and set the buffer index to this channel
651    //
652    ClearPixelsBuffers();
653    fBufIndex = startuid-1;
654
655    Bool_t display = gClient && fDisplay;
656    TGProgressBar* p = NULL;
657    Bool_t show = false, percent = false;
658    TString format;
659    Float_t min = 0, max = 0, pos = 0;
660    if ( display ) {
661        p = fDisplay->GetProgressBar();
662        format = p->GetFormat();
663        min = p->GetMin();
664        max = p->GetMax();
665        show = p->GetShowPos();
666        percent = p->UsePercent();
667        pos = p->GetPosition();
668
669        p->SetRange(0,100);
670        fDisplay->ShowProgressBarPosition(kTRUE,kFALSE,"Buffering... (%.1f%%)");
671    }
672
673    for( Long64_t i(0); i<nhits; i++) {
674        fHitsTree->GetEntry(i);
675        Int_t uid = fHit.fUID;
676
677        if ( uid >= startuid && uid <= lastuid ) {
678            fThetaFOVBuffer[uid-fBufIndex-1].push_back(fHit.fTheta);
679            fPhiFOVBuffer[uid-fBufIndex-1].push_back(fHit.fPhi);
680        }
681
682        //
683        // if the map index has not been filled, fMapIndex->fArray[0] == 0
684        // fill it
685        //
686        if ( !fMapIndex[0] )
687            fMapIndex[uid]++;
688
689        if ( !((i+1)%100000) || i == nhits-1) {
690            Float_t x = ((Double_t)(i+1)*100)/(Double_t)nhits;
691            // EsafMsg useless here
692            cout << '\r' << i+1 << " out of " << nhits << " hits processed ("
693                 << Nint(Floor(x)) << "%)" << flush;
694
695            if (display) {
696                p->SetPosition(x);
697                gClient->ProcessEventsFor(p);
698            }
699        }
700    }
701
702    fMapIndex[0] = 1;
703
704    cout << endl;
705    Msg(EsafMsg::Info) << "Buffer filled" << MsgDispatch;
706
707    if ( display ) {
708        p->SetRange(min,max);
709        p->SetPosition(pos);
710        fDisplay->ShowProgressBarPosition(show,percent,format);
711        gClient->ProcessEventsFor(p);
712    }
713
714    return kTRUE;
715
716}
717
718//______________________________________________________________________________
719Bool_t PixelMapBuilder::CheckPixelInBuffer(Int_t uid) {
720    // load id if it is out of current buffer
721
722    // check uid to be valid
723    if ( uid < 1 || uid > fHeader->fNpixels) {
724        MsgForm(EsafMsg::Warning,"PixelMapBuilder: uid out of range");
725        return kFALSE;
726    }
727
728    // transform in array's index
729    Int_t idx = uid-1;
730
731    Int_t length = BufferEntries();
732    // check if idx if in the current buffer
733    if ( fBufIndex == -1 || idx < fBufIndex || idx >= fBufIndex+length){
734        Msg(EsafMsg::Info) << "Channel out of buffer. Switching to the right buffer" << MsgDispatch;
735        Msg(EsafMsg::Info) << "Buffer " << idx/length << MsgDispatch;
736        // load the right buffer
737        return FillPixelsBuffer((idx/length)*length+1);
738    }
739
740    return kTRUE;
741}
742
743//______________________________________________________________________________
744Bool_t PixelMapBuilder::ClearPixelsBuffers() {
745    //
746    // clears buffers
747    //
748
749    for(size_t k(0); k<fThetaFOVBuffer.size(); k++) {
750        fThetaFOVBuffer[k].resize(0);
751        fPhiFOVBuffer[k].resize(0);
752    }
753    fBufIndex = kBufEmpty;
754
755    return kTRUE;
756}
757
758//______________________________________________________________________________
759Bool_t PixelMapBuilder::Fit(Int_t uid, Option_t *option) {
760
761    TString opt = option;
762    // check on the index if there are enough it for fitting
763    if ( fMapIndex.GetSize() &&              // index allocated
764         fMapIndex[0] &&                     // index filled
765         fMapIndex[uid] < fHitsThreshold ) { // px under thres
766        if ( !opt.Contains("Q") )
767            MsgForm(EsafMsg::Debug,"Pixel %d under threshold. Skipping", uid);
768        return kFALSE;
769    }
770
771    if ( !CheckPixelInBuffer(uid) ) {
772         MsgForm(EsafMsg::Debug,"Pixel %d not in buffer. Skipping", uid);
773        return kFALSE;
774    }
775
776    // clear histogramis
777    fH2ThetaPhi->Reset();
778    fH1FitTheta->Reset();
779    fH1FitPhi->Reset();
780    fH2FitThetaPhi->Reset();
781    // translate uid in buffer index and check number of hits
782    Int_t idx = UidToIndex(uid);
783
784    if ( !opt.Contains("Q") )
785        MsgForm(EsafMsg::Info,"Fit: %d entries for channel %d (idx = %d).",
786                                         fThetaFOVBuffer[idx].size(), uid, idx);
787
788    if ( (Int_t)fThetaFOVBuffer[idx].size() < fHitsThreshold ) {
789        if ( !opt.Contains("Q") )
790            MsgForm(EsafMsg::Debug,"Pixel %d under threshold. Skipping", uid);
791        return kFALSE;
792    }
793
794
795    // transform (th,phi)->(x,y)
796    Float_t x,y;
797
798    for(size_t j(0); j<fThetaFOVBuffer[idx].size(); j++) {
799        x = Sin(fThetaFOVBuffer[idx][j])
800            *Cos(fPhiFOVBuffer[idx][j]);
801        y = Sin(fThetaFOVBuffer[idx][j])
802            *Sin(fPhiFOVBuffer[idx][j]);
803
804        fH2ThetaPhi->Fill(x,y);
805    }
806
807    // look for the most populated bin
808    Int_t lx, ly, lz;
809    fH2ThetaPhi->GetMaximumBin(lx, ly, lz);
810
811    Float_t xCenter = fH2ThetaPhi->GetXaxis()->GetBinCenter(lx);
812    Float_t yCenter = fH2ThetaPhi->GetYaxis()->GetBinCenter(ly);
813
814    // calculate errors
815    Float_t dr = fH2ThetaPhi->GetXaxis()->GetBinWidth(1);
816    Float_t r = Sqrt(xCenter*xCenter+yCenter*yCenter);
817
818    // xmean,xxmean,ymean,yymean->th, therr, ph, phierr
819    Float_t thCenter = ASin(r);
820    Float_t phiCenter = ATan2(yCenter, xCenter);
821    Float_t thErr = dr/Cos(thCenter);
822    Float_t phiErr = dr/r;
823
824
825    // reset fit histos
826
827    // fill and fit theta
828    fH1FitTheta->GetXaxis()->SetLimits(thCenter-fNsigma*thErr, thCenter+fNsigma*thErr);
829    fH1FitTheta->GetXaxis()->SetRange();
830    // fill and fit phi
831    fH1FitPhi->GetXaxis()->SetLimits(phiCenter-fNsigma*phiErr, phiCenter+fNsigma*phiErr);
832    fH1FitPhi->GetXaxis()->SetRange();
833
834    // fill and fit theta
835    fH2FitThetaPhi->GetXaxis()->SetLimits(thCenter-fNsigma*thErr, thCenter+fNsigma*thErr);
836    fH2FitThetaPhi->GetYaxis()->SetLimits(phiCenter-fNsigma*phiErr, phiCenter+fNsigma*phiErr);
837    fH2FitThetaPhi->GetXaxis()->SetRange();
838    fH2FitThetaPhi->GetYaxis()->SetRange();
839
840    // save quadrant where phi is
841    Short_t quad = (Short_t)(phiCenter/TMath::PiOver2());
842
843    // fill
844    for(size_t j(0); j<fThetaFOVBuffer[idx].size(); j++) {
845        Float_t phi = fPhiFOVBuffer[idx][j];
846        Float_t theta = fThetaFOVBuffer[idx][j];
847        switch ( quad ) {
848            case 1:
849                if ( fPhiFOVBuffer[idx][j] > 3*TMath::PiOver2() )
850                    phi -= TwoPi();
851                break;
852            case 4:
853                if ( fPhiFOVBuffer[idx][j] < TMath::PiOver2() )
854                    phi += TwoPi();
855                break;
856            default:
857                break;
858        }
859
860        fH1FitTheta->Fill(theta);
861        fH1FitPhi->Fill(phi);
862        fH2FitThetaPhi->Fill(theta, phi);
863    }
864
865    //FIXME: there can be a ghost; add a TH2F histo and use
866    //       fH1FitTheta and fH1FitPhi as projections
867    //       (TH2::GetProjections())
868    //       this is a low priority fix
869
870    // fit
871    fH1FitTheta->Fit("gaus",opt);
872    fH1FitPhi->Fit("gaus",opt);
873
874    fFitUID = uid;
875    return kTRUE;
876
877}
878//______________________________________________________________________________
879TCanvas *PixelMapBuilder::DrawHistPad() const {
880    // Draws one of the buffer's histos
881
882    TCanvas *c = new TCanvas("PixelMapFit","Pixel Dist fit");
883
884    TText txt;
885    txt.SetTextFont(72);
886    txt.SetTextAlign(13);
887    txt.SetTextSize(5.5e-2);
888    txt.DrawTextNDC(0.01, 0.99,Form("Pixel %d",fFitUID));
889
890    TPad* pth = new TPad("PixelMapFit_Theta","PixelMapFit #Theta", 0., 0.475, 0.5, 0.95);
891    TPad* pphi = new TPad("PixelMapFit_Phi","PixelMapFit #Phi", 0, 0, 0.5, 0.475);
892    TPad* p2D = new TPad("PixelMapFit_2D","PixelMapFit #Theta-#Phi", 0.5, 0.475, 1., 0.95);
893    TPad* p2DFoV = new TPad("PixelMapFit_2D_FoV","PixelMapFit 2D Whole FoV", 0.5, 0., 1., 0.475);
894
895    pth->Draw();
896    pphi->Draw();
897    p2D->Draw();
898    p2DFoV->Draw();
899
900    pth->cd();
901    fH1FitTheta->Draw();
902
903    pphi->cd();
904    fH1FitPhi->Draw();
905
906    p2D->cd();
907    fH2FitThetaPhi->Draw("colz");
908
909    p2DFoV->cd();
910    fH2ThetaPhi->Draw("colz");
911
912    TEllipse *e = new TEllipse(0,0,Sin( fHeader->fFoV ));
913    e->SetBit(kCanDelete);
914    e->Draw();
915
916    return c;
917}
918
919//______________________________________________________________________________
920void PixelMapBuilder::DumpHit() const {
921    //
922    //
923    //
924    Msg(EsafMsg::Info) << "fHit.fUID    : " << fHit.fUID << MsgDispatch;
925    Msg(EsafMsg::Info) << "fHit.fTheta  : " << fHit.fTheta << MsgDispatch;
926    Msg(EsafMsg::Info) << "fHit.fPhi    : " << fHit.fPhi << MsgDispatch;
927    Msg(EsafMsg::Info) << "fHit.fLambda : " << fHit.fLambda << MsgDispatch;
928}
929
930//______________________________________________________________________________
931void PixelMapBuilder::WriteMap(const char* name) {
932    //
933    // Write the map to file in the ESAF format
934    //
935
936
937    ofstream fout;
938    if ( name ) {
939        fout.open(name);
940        if (!fout.is_open()) {
941            Msg(EsafMsg::Warning) << "Unable to open " <<  name << MsgDispatch;
942            return;
943        }
944        Msg(EsafMsg::Info) << "Saving map in " << name << MsgDispatch;
945    }
946
947    fout << "# " <<  endl;
948    fout << "# $Id: PixelMapBuilder.cc 2927 2011-06-16 18:02:41Z mabl $" << endl;
949    fout << "# PixelMapBuilder" << endl;
950    fout << "" << endl;
951    fout << "###############################################################################" << endl;
952    fout << "# ESAF: Euso Simulation and Analysis Framework                                #" << endl;
953    fout << "#                                                                             #" << endl;
954    fout << "#  Package: Optics                                                            #" << endl;
955    fout << "#  Coordinator: Alessandro.Thea                                               #" << endl;
956    fout << "#                                                                             #" << endl;
957    fout << "###############################################################################" << endl;
958    fout << "#" << endl;
959    fout << "# uid   theta   thetaRMS  phi   phiRMS" << endl;
960
961
962    for ( Int_t k(0); k<MapSize(); k++)
963        fout << Form("%d    %f    %f    %f    %f",
964                k+1, fMapTheta[k], fMapThetaRMS[k], fMapPhi[k], fMapPhiRMS[k]) << endl;
965
966    fout.close();
967}
Note: See TracBrowser for help on using the repository browser.