Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

PixelMapBuilder - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: PixelMapBuilder.cc,v 1.25 2005/10/15 10:12:59 thea Exp $
// A.Thea created Apr, 15 2004

#include "PixelMapBuilder.hh"
#include "EusoDetector.hh"
#include "DetectorTransportManager.hh"
#include "EusoElectronics.hh"
//#include "OAPxBunch.hh"
#include "Photon.hh"
#include "ParentPhoton.hh"
#include "DetectorGeometry.hh"
#include <TTree.h>
#include <TChain.h>
#include "TTimeStamp.h"
#include "TString.h"
#include <TCanvas.h>
#include <TF1.h>
#include <TText.h>
#include <TEllipse.h>

#include <TPolyMarker.h>
#include <TSystem.h>
#include <TGClient.h>
#include <TGProgressBar.h>
#include <TMutex.h>
#include "EGViewer.hh"

using namespace TMath;

//______________________________________________________________________________
//
// PixelMapBuilder
//
// Works in 2 modes
// 
// Create the photon's file for the anaysis
//
// fNphotons
// fDetector
// fTransporter
// fElectronics
// fHitsRootFile
// fHitsTree
// fFoV
// fHit
// fHeader
// 
//

ClassImp(PixelMapBuilder)

//______________________________________________________________________________
 PixelMapBuilder::PixelMapBuilder() : fHitsRootFile(0), fHitsTree(0), fHitsChain(0), 
    fFitHistsFile(0), fFitHistsTree(0), fFitUID(0), fH1FitTheta(0),
    fH1FitPhi(0), fH2FitThetaPhi(0), fH2ThetaPhi(0), fDisplay(0) {
    // 
    // ctor

    fMaxBufSize = 5000000/(sizeof(Float_t)*2);
    fBufferSize = 250 /*Mb*/ *1024*1024;
    fNsigma       = 1.5;

    fHitsThreshold = 50;
    fNbins        = 100;

    fBufIndex = kBufEmpty;

    fFullFoV = Pi()/6.;
    fFoV = Pi()/6.;
//    fBufferSize = 1;
}

//______________________________________________________________________________
 PixelMapBuilder::~PixelMapBuilder() {
    //
    // Destructor
    //
}


//______________________________________________________________________________
 Bool_t PixelMapBuilder::AddToChain( const char* name ) {

    if ( fHitsChain ) fHitsChain->Add(name,0);
    return kTRUE;
}

//______________________________________________________________________________
 Bool_t PixelMapBuilder::OpenChain( const char* name ) {

    fHitsChain = new TChain("oapxtree");

    AddToChain(name);
    fHitsChain->GetEntries();

    fHitsChain->SetBranchAddress("Hit", &fHit);
    fHitsChain->GetEntry();
    fHeader = (Header*)fHitsChain->GetTree()->GetUserInfo()->FindObject("PixelMapBuilder::Header")->Clone();

    fHitsTree = fHitsChain;
    return kTRUE;
}

//______________________________________________________________________________
 Bool_t PixelMapBuilder::OpenChain( const char** name, Int_t n ) {

    fHitsChain = new TChain("oapxtree");

    for(Int_t i(0); i<n; i++) {
        AddToChain(name[i]);
    }
    fHitsChain->GetEntries();

    fHitsChain->SetBranchAddress("Hit", &fHit);
    fHitsChain->GetEntry();
    fHeader = (Header*)fHitsChain->GetTree()->GetUserInfo()->FindObject("PixelMapBuilder::Header")->Clone();

    fHitsTree = fHitsChain;
    return kTRUE;
}

//______________________________________________________________________________
 Bool_t PixelMapBuilder::OpenRoot( const char* name, Bool_t write ) {
    //
    // Open rootfile with the correct name
    //
    CloseRoot();

    TString fname(name);
    TString ext = ".root";

    if ( write) {
        if ( fname.EndsWith(ext)) fname.Resize(fname.Length()-ext.Length());
        Msg(EsafMsg::Info) << "Opening root file " << fname << ext << MsgDispatch;
        fHitsRootFileName = fname+ext;
        fHitsRootFile = new TFile(fHitsRootFileName,"CREATE");

        if (fHitsRootFile->IsZombie()) {
            Msg(EsafMsg::Warning) << MsgDispatch;
            Msg(EsafMsg::Warning) << "Error opening file " << fname << MsgDispatch;
            Msg(EsafMsg::Warning) << "Probably it already exists or there is no " << MsgDispatch;
            Msg(EsafMsg::Warning) << "write access permission to this directory" << MsgDispatch;
            Msg(EsafMsg::Warning) << MsgDispatch;

            // try adding date and time to the name
            TTimeStamp time;
            TString now = time.AsString("s");
            now[now.Last(' ')]=':';
            TString newFileName = Form("%s-", fname.Data())+now+ext;
            fHitsRootFile = new TFile(newFileName,"CREATE");
            if (fHitsRootFile->IsZombie()) {
                MsgForm(EsafMsg::Warning,"Cannot open root file");
                return kFALSE;
            }
        }
        Msg(EsafMsg::Info) << "Root file " << fHitsRootFile->GetName() << " successfully openedn" << MsgDispatch;

        fHitsTree = new TTree("oapxtree","Angle Pixel Map photons' database");
        fHitsTree->Branch("Hit",&fHit,"fUID/I:fTheta/F:fPhi/F:fLambda/F");
        fHeader = new Header;
        fHitsTree->GetUserInfo()->Add(fHeader);

        return kTRUE;
    } else {
        if ( fname.EndsWith(ext)) fname.Resize(fname.Length()-ext.Length());
        fHitsRootFile = new TFile((fname+ext));

        if ( fHitsRootFile->IsZombie() ) return kFALSE;

        fHitsTree = (TTree*)fHitsRootFile->Get("oapxtree");
        fHitsTree->SetBranchAddress("Hit", &fHit);

        fHeader = (Header*)fHitsTree->GetUserInfo()->FindObject("PixelMapBuilder::Header");

        if ( !fHeader ) {
            Msg(EsafMsg::Warning) << "No header in the rootfile" << MsgDispatch;
            return kFALSE;
        }

        return kTRUE;
    }
}

//______________________________________________________________________________
 Bool_t PixelMapBuilder::CloseRoot() { 
    // 
    // close root file
    //

    if ( fHitsRootFile && fHitsChain ) {
        Msg(EsafMsg::Warning) << "Both fHitsRootFile and fHitsChain are present. Something's wronmg. Anorting" << MsgDispatch;
        return kFALSE;
    }
    
    if ( fHitsRootFile ) {
        if ( fHitsRootFile->IsWritable() ) {
            fHitsRootFile->Write();

            SafeDelete(fHitsTree);
            fHeader = 0;

            Msg(EsafMsg::Info) << "Root file saved" << MsgDispatch;
        } else {
            SafeDelete(fHitsTree);
            fHeader = 0;
        }
        fHitsRootFile->Close();
        
        Msg(EsafMsg::Info) << "Root file closed" << MsgDispatch;
        SafeDelete( fHitsRootFile )
    } else if (fHitsChain) {
        // chain, reading by default
        SafeDelete(fHitsChain);
        SafeDelete(fHeader);
        fHitsTree = 0;
    }

    return kTRUE;
}


//______________________________________________________________________________
 void PixelMapBuilder::TrackPhotons() {
    // 
    // build the map
    //

    // disable electronics simulation, just find the photon's pad
    fElectronics->EnableSimulation(kFALSE);
    fHit.fUID = 0;
    fHit.fTheta = 0;
    fHit.fPhi = 0;
    fHit.fLambda = 0;
    fPupil.SetThetaFOVMax(fFoV);
    fHeader->fFoV = fFoV;
    fHeader->fFullFoV = fDetector->GetGeometry()->GetFoV();

    Msg(EsafMsg::Info) << "Photons to track: " << fNphotons << MsgDispatch;
    Msg(EsafMsg::Info) << "PixelMapBuilder: tracing..." << MsgDispatch;
    Long_t hits(0);
    for(Long_t i(0); i<fNphotons; i++) {

        Photon *p = fPupil.GetPhoton();

        fTransporter->Transport(p);

        fHeader->fNtracked++;

        if ( p->pixelUid > 0 ) { // FIXME: should be > 0
            TVector3 dir = p->parent->Dir();

            hits++;
            fHeader->fNhits++;

            fHit.fUID =  p->pixelUid;
            fHit.fTheta = dir.Theta();
            fHit.fPhi = dir.Phi();
            fHit.fLambda = p->wl;

            fHitsTree->Fill();
            if ( i%1000000 == 0 ) fHitsTree->AutoSave();

            if ( fHitsTree->GetCurrentFile() != fHitsRootFile ) {
                fHitsRootFile = fHitsTree->GetCurrentFile();
                fHeader->fNtracked = 0;
                fHeader->fNhits = 0;
            }

        }

        if ( (i+1)%1000 == 0) {
            Msg(EsafMsg::Info) << 'r' <<  MsgFlush;
            Msg(EsafMsg::Info) << Form("%d photons traced, %d hits on the pixels",
                    (i+1), hits) << MsgFlush;
        }
    }
    Msg(EsafMsg::Info) << MsgDispatch;
    Msg(EsafMsg::Info) << "Completed" << MsgDispatch;

    MsgForm(EsafMsg::Info, "%d photons traced.",fNphotons);
    MsgForm(EsafMsg::Info, "%d photons hit the pixels.",hits);
    // FIXME: should be false if it was before TrackPhotons()
    fElectronics->EnableSimulation(kTRUE);
}
   

//______________________________________________________________________________
 Bool_t PixelMapBuilder::MakePhotonsFile( const char* name, Long_t n) {


    fNphotons = n;

    
    Msg(EsafMsg::Info) << "Building Detector" << MsgDispatch;

    fDetector = new EusoDetector(); 
    fTransporter = fDetector->GetDetectorTransportManager();
    fElectronics = fDetector->GetEusoElectronics();

    Msg(EsafMsg::Debug) << "TransportManager Modules" << MsgDispatch;
    Msg(EsafMsg::Debug) << "Optics: " << fTransporter->GetOptics()->IsA()->GetName() << MsgDispatch;
    Msg(EsafMsg::Debug) << "FocalPlane: " << fTransporter->GetFocalPlane()->IsA()->GetName() << MsgDispatch;
    Msg(EsafMsg::Info) <<  "Building complete" << MsgDispatch;

    if ( !OpenRoot(name, kTRUE) ) {
        Msg(EsafMsg::Warning) << "Unable to open the root file. Aborting" << MsgDispatch;
        return kFALSE;
    }

    // setup the tree

    fHeader->fNpixels = fElectronics->NumOfCh();

    // loop here
    TrackPhotons();

    CloseRoot();
    SafeDelete(fHeader);

    return kTRUE;
}

//______________________________________________________________________________
 Bool_t PixelMapBuilder::MakeMap( const char* hits, const char* mapfile) {

    if ( !OpenRoot(hits) ) return kFALSE; 

    InitFit();

//    Int_t first = 1;
//    Int_t last = fHeader->fNpixels;
   
    ComputeMap( 0, 0, "WEQ0" );

    ClearFit();

    CloseRoot();

    WriteMap(mapfile);

    return kTRUE;
}

//______________________________________________________________________________
 void PixelMapBuilder::ComputeMap(Int_t first, Int_t last, const char* opt) {
    //
    // Fits the histograms 
    //

    TString options(opt);
    
    if ( !fMapTheta.GetSize() ) {
        Msg(EsafMsg::Warning) << "Map not defined. Aborting" << MsgDispatch;
        return;
    }
        
    if ( !first ) first = 1;
    if ( !last ) last = fHeader->fNpixels;

    Msg(EsafMsg::Info) << "Processing map, pixels " << first << 
                          " - " << last << MsgDispatch;
    TF1 *g;
    Double_t parTh[3], parPhi[3];

    //
    // Tentative display
    //
    Bool_t display = gClient && fDisplay;
    TCanvas *cPoly, *cHist;
    TPolyMarker* m;
    TGProgressBar* p;
    TH2F* h; 
    Stat_t entries(0);

    if ( display ) {
        //
        // Polymarker display
        //
        cPoly = fDisplay->FindAddTab("PixelMapDisplayPoly","Map projection");
        cPoly->cd();
        cPoly->DrawFrame(-Sin(fHeader->fFoV),-Sin(fHeader->fFoV),
                      Sin(fHeader->fFoV), Sin(fHeader->fFoV));
        m = new TPolyMarker(fHeader->fNpixels); 
//        m->SetMarkerStyle(7);
        m->SetBit(kCanDelete);
        m->Draw();
        
        //
        // Histogram display
        //
        cHist = fDisplay->FindAddTab("PixelMapDisplayHist","Histogram");
        cHist->cd();
        h = new TH2F("PixelMap","PixelMap",
                           100, -Sin(fHeader->fFoV),Sin(fHeader->fFoV),
                           100, -Sin(fHeader->fFoV),Sin(fHeader->fFoV));
        h->SetDirectory(0);
        h->SetBit(kCanDelete);
        h->SetXTitle("Sin(#theta) #times Cos(#phi)");
        h->SetYTitle("Sin(#theta) #times Sin(#phi)");
        h->Draw("col");

        p = fDisplay->GetProgressBar();
        p->SetRange(1,fHeader->fNpixels);
        fDisplay->ShowProgressBarPosition(kTRUE,kFALSE,"%.0f pixels");

        // no fit display if the display is active.
        if ( !options.Contains("0")) options.Append("0");
    }
    

    for(Int_t k(first); k<=last; k++){
        // 
        // fit mut be called before UidToIndex because the buffer 
        // may be not be filled yet
        // 
        Bool_t fit = Fit(k, options);
        Int_t j=UidToIndex(k);

        // 
        // fit pixel data
        // 
        if (!fit) {
            fMapTheta[k-1]    = -1; // -2*TwoPi()
            fMapThetaRMS[k-1] = -1; // 0
            fMapPhi[k-1]      = -1; // -2*TwoPi()
            fMapPhiRMS[k-1]   = -1; // 0
        } else {

            // retrieve results from histos
            g = fH1FitTheta->GetFunction("gaus");

            if ( g ) g->GetParameters(parTh);
            // if g doesn't exists put dummy values
            else  parTh[0] = parTh[1] = parTh[2] = -1;

            g = fH1FitPhi->GetFunction("gaus");
            if ( g ) g->GetParameters(parPhi);
            else  parPhi[0] = parPhi[1] = parPhi[2] = -1;

            // fill the containers
            fMapTheta[k-1]    = parTh[2];
            fMapThetaRMS[k-1] = parTh[2];
            fMapPhi[k-1]      = parPhi[1];
            fMapPhiRMS[k-1]   = parPhi[2];
        }
        if ( fFitHistsTree ) fFitHistsTree->Fill();
        
        if ( display && fit ) {

            Float_t x = Sin(fMapTheta[k])*Cos(fMapPhi[k]);
            Float_t y = Sin(fMapTheta[k])*Sin(fMapPhi[k]);
            m->SetPoint(j,x,y);
            h->Fill(x,y);
            
        }

        if ( !(k%100) || k == last ) { 
            cout  << 'r' << Form("Pixel %d of %d completed (%d%%) ",
                  k, fHeader->fNpixels, k*100/fHeader->fNpixels) << flush;

            //
            // update the display
            // 
            if (display && ( h->GetEntries() != entries || k == last ) ) {
                cPoly->Modified();
                cPoly->Update();
                cHist->Modified();
                cHist->Update();
                p->SetPosition((Float_t)k);
                if (!TThread::Self())
                    gSystem->ProcessEvents();

                entries = h->GetEntries();
            }
        }
         
     }
     cout << endl;


    
}

//______________________________________________________________________________
 void PixelMapBuilder::OptimizeBuffer() {
    //
    // Tries to evaluate the expectes average number of hits per pixel and the
    // length of the buffers to match fBufferSize. Then reserves the
    // corresponding amount of memory in the bugger
    //

    if (!fHitsTree) return;

    Long64_t entries = fHitsTree->GetEntries();

    // 
    // expected average number of hits per pixel
    // 
    Int_t average = Nint(entries*(1-Cos(fFullFoV))/(fHeader->fNpixels*(1-Cos(fHeader->fFoV))));
    Msg(EsafMsg::Info) << "Expected average number of entries per pixel: " << average << MsgDispatch;
    
    // 
    // fBufferSize = sizeof(float)*2*average*length
    //
    Int_t length = fBufferSize/((Int_t)sizeof(float)*2*average);
    length = Min(length,fHeader->fNpixels);
    Msg(EsafMsg::Info) << "Buffer lenght set to " << length << MsgDispatch;
    Msg(EsafMsg::Info) << "Expected buffer size in memory " << 
                          (sizeof(float)*2*average*length)/(1024*1024) << " Mb"<< MsgDispatch;

    // prepare the buffers
    fThetaFOVBuffer.resize(length);
    fPhiFOVBuffer.resize(length);
    
    for(size_t k(0); k<fThetaFOVBuffer.size(); k++) {
        fThetaFOVBuffer[k].reserve(average);
        fPhiFOVBuffer[k].reserve(average);
    }

}

//______________________________________________________________________________
 Bool_t PixelMapBuilder::InitFit() {

    if ( !fHitsTree ) return kFALSE;
    
    // clear the content of the buffers.
    ClearPixelsBuffers();
    // resize the map arrays
    fMapTheta.Set(fHeader->fNpixels);
    fMapThetaRMS.Set(fHeader->fNpixels);
    fMapPhi.Set(fHeader->fNpixels);
    fMapPhiRMS.Set(fHeader->fNpixels);
    
    fMapIndex.Set(fHeader->fNpixels+1);
    fMapIndex[0] = 0;
    
    OptimizeBuffer();
    
    // note: upper limits in theta is a little bit more than Pi/6
    Bool_t status = TH1::AddDirectoryStatus();
    TH1::AddDirectory(kFALSE);

    fH1FitTheta  = new TH1F("oapx_thfit","Temp fit histogram for Theta",fNbins, 0, 1);
    fH1FitTheta->SetXTitle("theta");
    fH1FitTheta->SetYTitle("hits");
    fH1FitPhi    = new TH1F("oapx_phifit","Temp fit histogram for Phi", fNbins, 0, 1);
    fH1FitPhi->SetXTitle("phi");
    fH1FitPhi->SetYTitle("hits");
    fH2FitThetaPhi = new TH2F("oapx_thphifit","Temp fit histogram for Theta and Phi",
                             fNbins, 0, 1, fNbins, 0, 1);
    fH2FitThetaPhi->SetXTitle("theta");
    fH2FitThetaPhi->SetYTitle("phi");

    Float_t r = Sin(fHeader->fFoV*1.1);
    fH2ThetaPhi =  new TH2F("oapx_tp","tp distribution projected on a plane",
                                     fNbins, -r, r, fNbins, -r, r);
    fH2FitThetaPhi->SetXTitle("Sin(#theta) #times Cos(#phi)");
    fH2FitThetaPhi->SetYTitle("Sin(#theta) #times Sin(#phi)");
    
    TH1::AddDirectory(status);

    if ( !fFitHistsFileName.IsNull() ) {
        TString ext = ".root";
        TString fname = fFitHistsFileName;

        if ( fname.EndsWith(ext)) fname.Resize(fname.Length()-ext.Length());
        Msg(EsafMsg::Info) << "Saving histograms in " << fname << ext << MsgDispatch;
        fFitHistsFileName = fname+ext;
        fFitHistsFile = new TFile(fFitHistsFileName,"RECREATE");

        fFitHistsTree = new TTree("pxmap_hist","Pixel Map histograms");
        fFitHistsTree->Branch("fTheta",&fH1FitTheta, 100000, 0);
        fFitHistsTree->Branch("fPhi",&fH1FitPhi, 100000, 0);
        fFitHistsTree->Branch("fThetaPhi",&fH2FitThetaPhi, 100000, 0);
        fFitHistsTree->Branch("fThetaPhiWholeFoV",&fH2ThetaPhi, 100000, 0);
        fFitHistsTree->Branch("fUID",&fFitUID, "fUID/I");

    }

    fFitUID = 0;
    
    return kTRUE;
}

//______________________________________________________________________________
 Bool_t PixelMapBuilder::ClearFit() {
    //
    // Clears the histograms usend in the fit procedure
    //

    if ( fFitHistsTree ) {
        fFitHistsFile = fFitHistsTree->GetCurrentFile();
        fFitHistsFile->cd();
        fFitHistsTree->Write();
        SafeDelete(fFitHistsTree);
        SafeDelete(fFitHistsFile);
    }

    fMapIndex.Set(0);
//    fHeader = 0;

    SafeDelete(fH1FitTheta);
    SafeDelete(fH1FitPhi);
    SafeDelete(fH2FitThetaPhi);
    SafeDelete(fH2ThetaPhi);

    ClearPixelsBuffers();

    return kTRUE;
}

//______________________________________________________________________________
 Bool_t PixelMapBuilder::FillPixelsBuffer(Int_t startuid ) {
    //
    // Fills the buffer starting from 
    //

    // 
    // Check boundaries
    // 
    if ( startuid >  fHeader->fNpixels ){
        Msg(EsafMsg::Warning) << "First out of bounds" << endl;
        return kFALSE;
    }


    Int_t length = BufferEntries();
    Int_t lastuid = (startuid+length-1) <= fHeader->fNpixels ? 
                 (startuid+length-1) : (fHeader->fNpixels-1);

    Long64_t nhits = fHitsTree->GetEntriesFast();
    //
    // clear the buffer and set the buffer index to this channel
    //
    ClearPixelsBuffers();
    fBufIndex = startuid-1;
    
    Bool_t display = gClient && fDisplay;
    TGProgressBar* p;
    Bool_t show, percent;
    TString format;
    Float_t min, max, pos;
    if ( display ) {
        p = fDisplay->GetProgressBar();
        format = p->GetFormat();
        min = p->GetMin();
        max = p->GetMax();
        show = p->GetShowPos();
        percent = p->UsePercent();
        pos = p->GetPosition();
        
        p->SetRange(0,100);
        fDisplay->ShowProgressBarPosition(kTRUE,kFALSE,"Buffering... (%.1f%%)");
    }

    for( Long64_t i(0); i<nhits; i++) {
        fHitsTree->GetEntry(i);
        Int_t uid = fHit.fUID; 

        if ( uid >= startuid && uid <= lastuid ) {
            fThetaFOVBuffer[uid-fBufIndex-1].push_back(fHit.fTheta);
            fPhiFOVBuffer[uid-fBufIndex-1].push_back(fHit.fPhi);
        }

        //
        // if the map index has not been filled, fMapIndex->fArray[0] == 0
        // fill it
        //
        if ( !fMapIndex[0] )
            fMapIndex[uid]++;

        if ( !((i+1)%100000) || i == nhits-1) {
            Float_t x = ((Double_t)(i+1)*100)/(Double_t)nhits;
            // EsafMsg useless here
            cout << 'r' << i+1 << " out of " << nhits << " hits processed ("
                 << Nint(Floor(x)) << "%)" << flush;

            if (display) {
                p->SetPosition(x);
                gClient->ProcessEventsFor(p);
            }
        }
    }

    fMapIndex[0] = 1;

    cout << endl;
    Msg(EsafMsg::Info) << "Buffer filled" << MsgDispatch;

    if ( display ) {
        p->SetRange(min,max);
        p->SetPosition(pos);
        fDisplay->ShowProgressBarPosition(show,percent,format);
        gClient->ProcessEventsFor(p);
    }

    return kTRUE;

}

//______________________________________________________________________________
 Bool_t PixelMapBuilder::CheckPixelInBuffer(Int_t uid) {
    // load id if it is out of current buffer

    // check uid to be valid
    if ( uid < 1 || uid > fHeader->fNpixels) {
        MsgForm(EsafMsg::Warning,"PixelMapBuilder: uid out of range");
        return kFALSE;
    }

    // transform in array's index
    Int_t idx = uid-1;
        
    Int_t length = BufferEntries();
    // check if idx if in the current buffer
    if ( fBufIndex == -1 || idx < fBufIndex || idx >= fBufIndex+length){
        Msg(EsafMsg::Info) << "Channel out of buffer. Switching to the right buffer" << MsgDispatch;
        Msg(EsafMsg::Info) << "Buffer " << idx/length << MsgDispatch;
        // load the right buffer
        return FillPixelsBuffer((idx/length)*length+1);
    }

    return kTRUE;
}

//______________________________________________________________________________
 Bool_t PixelMapBuilder::ClearPixelsBuffers() {
    //
    // clears buffers
    //

    for(size_t k(0); k<fThetaFOVBuffer.size(); k++) {
        fThetaFOVBuffer[k].resize(0);
        fPhiFOVBuffer[k].resize(0);
    }
    fBufIndex = kBufEmpty;

    return kTRUE;
}

//______________________________________________________________________________
 Bool_t PixelMapBuilder::Fit(Int_t uid, Option_t *option) {

    TString opt = option;

    // check on the index if there are enough it for fitting
    if ( fMapIndex.GetSize() &&              // index allocated
         fMapIndex[0] &&                     // index filled
         fMapIndex[uid] < fHitsThreshold ) { // px under thres
        if ( !opt.Contains("Q") ) 
            MsgForm(EsafMsg::Debug,"Pixel %d under threshold. Skipping", uid);
        return kFALSE;
    }

    if ( !CheckPixelInBuffer(uid) ) 
        return kFALSE;
    

    // clear histogramis
    fH2ThetaPhi->Reset(); 
    fH1FitTheta->Reset();
    fH1FitPhi->Reset();
    fH2FitThetaPhi->Reset();
    // translate uid in buffer index and check number of hits
    Int_t idx = UidToIndex(uid);

    if ( !opt.Contains("Q") ) 
        MsgForm(EsafMsg::Info,"Fit: %d entries for channel %d (idx = %d).",
                                         fThetaFOVBuffer[idx].size(), uid, idx);
    
    if ( (Int_t)fThetaFOVBuffer[idx].size() < fHitsThreshold ) {
        if ( !opt.Contains("Q") ) 
            MsgForm(EsafMsg::Debug,"Pixel %d under threshold. Skipping", uid);
        return kFALSE;
    }


    // transform (th,phi)->(x,y) 
    Float_t x,y;

    for(size_t j(0); j<fThetaFOVBuffer[idx].size(); j++) {
        x = Sin(fThetaFOVBuffer[idx][j])
            *Cos(fPhiFOVBuffer[idx][j]); 
        y = Sin(fThetaFOVBuffer[idx][j])
            *Sin(fPhiFOVBuffer[idx][j]); 

        fH2ThetaPhi->Fill(x,y);
    }

    // look for the most populated bin
    Int_t lx, ly, lz;
    fH2ThetaPhi->GetMaximumBin(lx, ly, lz);

    Float_t xCenter = fH2ThetaPhi->GetXaxis()->GetBinCenter(lx);
    Float_t yCenter = fH2ThetaPhi->GetYaxis()->GetBinCenter(ly);

    // calculate errors
    Float_t dr = fH2ThetaPhi->GetXaxis()->GetBinWidth(1);
    Float_t r = Sqrt(xCenter*xCenter+yCenter*yCenter);

    // xmean,xxmean,ymean,yymean->th, therr, ph, phierr
    Float_t thCenter = ASin(r); 
    Float_t phiCenter = ATan2(yCenter, xCenter); 
    Float_t thErr = dr/Cos(thCenter);
    Float_t phiErr = dr/r;


    // reset fit histos

    // fill and fit theta
    fH1FitTheta->GetXaxis()->SetLimits(thCenter-fNsigma*thErr, thCenter+fNsigma*thErr);
    fH1FitTheta->GetXaxis()->SetRange();
    // fill and fit phi
    fH1FitPhi->GetXaxis()->SetLimits(phiCenter-fNsigma*phiErr, phiCenter+fNsigma*phiErr);
    fH1FitPhi->GetXaxis()->SetRange();

    // fill and fit theta
    fH2FitThetaPhi->GetXaxis()->SetLimits(thCenter-fNsigma*thErr, thCenter+fNsigma*thErr);
    fH2FitThetaPhi->GetYaxis()->SetLimits(phiCenter-fNsigma*phiErr, phiCenter+fNsigma*phiErr);
    fH2FitThetaPhi->GetXaxis()->SetRange();
    fH2FitThetaPhi->GetYaxis()->SetRange();

    // save quadrant where phi is 
    Short_t quad = (Short_t)(phiCenter/TMath::PiOver2());

    // fill 
    for(size_t j(0); j<fThetaFOVBuffer[idx].size(); j++) {
        Float_t phi = fPhiFOVBuffer[idx][j];
        Float_t theta = fThetaFOVBuffer[idx][j];
        switch ( quad ) {
            case 1:
                if ( fPhiFOVBuffer[idx][j] > 3*TMath::PiOver2() )
                    phi -= TwoPi();
                break;
            case 4:
                if ( fPhiFOVBuffer[idx][j] < TMath::PiOver2() )
                    phi += TwoPi();
                break;
            default:
                break;
        }
        
        fH1FitTheta->Fill(theta);
        fH1FitPhi->Fill(phi);
        fH2FitThetaPhi->Fill(theta, phi);
    }

    //FIXME: there can be a ghost; add a TH2F histo and use 
    //       fH1FitTheta and fH1FitPhi as projections
    //       (TH2::GetProjections())
    //       this is a low priority fix

    // fit
    fH1FitTheta->Fit("gaus",opt);
    fH1FitPhi->Fit("gaus",opt);

    fFitUID = uid;
    return kTRUE;

}
//______________________________________________________________________________
 TCanvas *PixelMapBuilder::DrawHistPad() const {
    // Draws one of the buffer's histos

    TCanvas *c = new TCanvas("PixelMapFit","Pixel Dist fit");

    TText txt;
    txt.SetTextFont(72);
    txt.SetTextAlign(13);
    txt.SetTextSize(5.5e-2);
    txt.DrawTextNDC(0.01, 0.99,Form("Pixel %d",fFitUID));

    TPad* pth = new TPad("PixelMapFit_Theta","PixelMapFit #Theta", 0., 0.475, 0.5, 0.95);
    TPad* pphi = new TPad("PixelMapFit_Phi","PixelMapFit #Phi", 0, 0, 0.5, 0.475);
    TPad* p2D = new TPad("PixelMapFit_2D","PixelMapFit #Theta-#Phi", 0.5, 0.475, 1., 0.95);
    TPad* p2DFoV = new TPad("PixelMapFit_2D_FoV","PixelMapFit 2D Whole FoV", 0.5, 0., 1., 0.475);

    pth->Draw();
    pphi->Draw();
    p2D->Draw();
    p2DFoV->Draw();

    pth->cd();
    fH1FitTheta->Draw();

    pphi->cd();
    fH1FitPhi->Draw();

    p2D->cd();
    fH2FitThetaPhi->Draw("colz");

    p2DFoV->cd();
    fH2ThetaPhi->Draw("colz");

    TEllipse *e = new TEllipse(0,0,Sin( fHeader->fFoV ));
    e->SetBit(kCanDelete);
    e->Draw();

    return c;
}

//______________________________________________________________________________
 void PixelMapBuilder::DumpHit() const {
    //
    //
    //
    Msg(EsafMsg::Info) << "fHit.fUID    : " << fHit.fUID << MsgDispatch;
    Msg(EsafMsg::Info) << "fHit.fTheta  : " << fHit.fTheta << MsgDispatch;
    Msg(EsafMsg::Info) << "fHit.fPhi    : " << fHit.fPhi << MsgDispatch;
    Msg(EsafMsg::Info) << "fHit.fLambda : " << fHit.fLambda << MsgDispatch;
}

//______________________________________________________________________________
 void PixelMapBuilder::WriteMap(const char* name) {
    //
    // Write the map to file in the ESAF format
    //


    ofstream fout;
    if ( name ) {
        fout.open(name);
        if (!fout.is_open()) {
            Msg(EsafMsg::Warning) << "Unable to open " <<  name << MsgDispatch;
            return;
        }
        Msg(EsafMsg::Info) << "Saving map in " << name << MsgDispatch;
    }

    fout << "# " <<  endl;
    fout << "# $Id: PixelMapBuilder.cc,v 1.25 2005/10/15 10:12:59 thea Exp $" << endl;
    fout << "# PixelMapBuilder" << endl;
    fout << "" << endl;
    fout << "###############################################################################" << endl;
    fout << "# ESAF: Euso Simulation and Analysis Framework                                #" << endl;
    fout << "#                                                                             #" << endl;
    fout << "#  Package: Optics                                                            #" << endl;
    fout << "#  Coordinator: Alessandro.Thea                                               #" << endl;
    fout << "#                                                                             #" << endl;
    fout << "###############################################################################" << endl;
    fout << "#" << endl;
    fout << "# uid   theta   thetaRMS  phi   phiRMS" << endl;


    for ( Int_t k(0); k<MapSize(); k++)
        fout << Form("%d    %f    %f    %f    %f",
                k+1, fMapTheta[k], fMapThetaRMS[k], fMapPhi[k], fMapPhiRMS[k]) << endl;

    fout.close();
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005