Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

OAPxPlayer - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: OAPxPlayer.cc,v 1.11 2005/01/12 18:45:52 thea Exp $
// A.Thea created Apr, 21 2004

#include "OAPxPlayer.hh"
#include "ERunParameters.hh"
#include "Riostream.h"
#include "TBranch.h"
#include "TF1.h"

ClassImp(OAPxPlayer)

//______________________________________________________________________________
OAPxPlayer::OAPxPlayer( TTree* t ) : fTree(t), fBunch(0),// fPixelMap(0),
    // ctor
    fH1DistTheta(0), fH1DistPhi(0), fH1FitTheta(0), fH1FitPhi(0), fH2DistXY(0) {

    fNSigma       = 1.5;
    fFitThreshold = 50;
    fNBins        = 500;
    fNBinsFit     = 100;
        
    // -1 means the buffer is empty
    fFirst = -1;
    fBufferSize = 1;
    // max size = 500 Mb
    fMaxSize = 50000000/sizeof(Float_t);
    Printf("MaxSize = %d", fMaxSize);

    TBranch *bBunch = fTree->GetBranch("hitbunches");
    if (!bBunch) {
        Printf("hitsbunch not found. Init aborted");
        return;
    }
    bBunch->SetAddress(&fBunch);
        
    Init();
}

//______________________________________________________________________________
OAPxPlayer::~OAPxPlayer() {
    // dtor
    Printf("OAPxPlayer Destructor called");
//    if ( fPixelMap ) delete fPixelMap;

    if ( fH1DistTheta ) delete fH1DistTheta;
    if ( fH1DistPhi ) delete fH1DistPhi;
    if ( fH1FitTheta ) delete fH1FitTheta;
    if ( fH1FitPhi ) delete fH1FitPhi;
    if ( fH2DistXY ) delete fH2DistXY;
}

//______________________________________________________________________________
 void OAPxPlayer::Init() {
// initialize buffers
    Printf("Init()");

    // check all the bunches to have the same fPixels
    Int_t nPixels;
    TBranch *b_fNumChannels = fTree->GetBranch("fNumChannels");
    fTree->SetMakeClass(1);
    b_fNumChannels->SetAddress(&nPixels);
    b_fNumChannels->GetEntry();
    fNumPixels = nPixels;
    
    for(Int_t i(0); i<fTree->GetEntriesFast(); i++){
        b_fNumChannels->GetEntry();
        if (fNumPixels != nPixels) {
            Printf("Wrong number of pixel in bunch %d.nAborting...",i);
            return;
        }
            
    }
    fTree->SetMakeClass(0);

    Printf("Number of pixels: %d", fNumPixels);

    // create the map containers
    
    fTheta.Set(fNumPixels);
    fThetaRMS.Set(fNumPixels);
    fPhi.Set(fNumPixels);
    fPhiRMS.Set(fNumPixels);
    
    
    // evaluate optimal buffer size
    fTree->GetEntry();
    Long_t nHits = (Long_t)fTree->GetEntriesFast()*fBunch->GetNumHits();
    Int_t nHitPPx = nHits/fNumPixels+1;
    fBufferSize = fMaxSize/(2*nHitPPx);
    Printf("Evaluated BufferSize = %d", fBufferSize);
    // buffer larger than number of pixels isn't useful
    if (fBufferSize > fNumPixels ) fBufferSize = fNumPixels;
    Printf("Number of Hits : %d", nHits);
    Printf("Average Hit per Pixel: %d", nHitPPx);
    Printf("BufferSize set to %d", fBufferSize); 
    
    // buffers
    fThetaFOVBuffer.resize(fBufferSize);
    fPhiFOVBuffer.resize(fBufferSize);


    // histograms
    // note: upper limits in theta is a little bit more than Pi/6
    fH1DistTheta = new TH1F("oapx_theta","Temp theta distribution",
                         fNBins, 0, TMath::Pi()/5);
    fH1DistPhi   = new TH1F("oapx_phi","Temp phi histogram",
                         fNBins, -TMath::Pi(), TMath::Pi());
    fH1FitTheta  = new TH1F("oapx_thfit","Temp fit histogram for Theta",
                         fNBinsFit, 0, 1);
    fH1FitPhi    = new TH1F("oapx_phifit","Temp fit histogram for Phi",
                         fNBinsFit, 0, 1);
    // sligthly more than EUSO FOV FIXME: should it be read from file?
    Float_t r = TMath::Sin(TMath::Pi()/5);
    fH2DistXY    = new TH2F("oapx_xy","XY distribution on cartesian plane",
                         100, -r, r, 100, -r, r);
}

//______________________________________________________________________________
 Bool_t OAPxPlayer::LoadBuffer(Int_t firstid, Long_t entries) {
    // load channel data from first to first+fBufferSize
    if (firstid < 1 || firstid > fNumPixels) {
        Printf("LoadBuffer: firstid out of range");
        return kFALSE;
    }
        
    Printf("Loading buffer from channel %d", firstid);
    // if entries == 0 load all events
    if (entries == 0 ) entries = (Long_t)fTree->GetEntriesFast();
    
    // clear buffers
    ClearBuffers();
    fFirst = firstid-1;

    OAPxHit *h;
    
    for(Int_t i(0); i<entries; i++) {
        fTree->GetEntry(i);

        Printf("Event %d, hits %d",i,fBunch->GetNumHits());
        for (Int_t j(0); j<fBunch->GetNumHits(); j++){

            h = fBunch->GetHit(j);
            Int_t id = h->ChUID()-1;
            if ( id >= fFirst && id < fFirst+fBufferSize ) {
                fThetaFOVBuffer[id-fFirst].push_back(h->ThetaFOV());
                fPhiFOVBuffer[id-fFirst].push_back(h->PhiFOV());
            } 
        } 
    }
     return kTRUE;
}

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

    // check uid to be valid
    if ( uid < 1 || uid > fNumPixels) {
        Printf("LoadChannel: firstid out of range");
        return kFALSE;
    }
    // transform in array's index
    uid--;
    // check if uid if in the current buffer
    if ( fFirst == -1 || uid < fFirst || uid >= fFirst+fBufferSize){
        Printf("Channel out of buffer. Loading right one.");
        Printf("buffer number %d", uid/fBufferSize);
        // load the right buffer
        return LoadBuffer((uid/fBufferSize)*fBufferSize+1);
    }

    return kTRUE;
}

//______________________________________________________________________________
 void OAPxPlayer::ClearBuffers() {
    // clears buffers

    for(Int_t k(0); k<fBufferSize; k++) {
        fThetaFOVBuffer[k].resize(0);
        fPhiFOVBuffer[k].resize(0);
    }
    fFirst = -1;
}

//______________________________________________________________________________
 Bool_t OAPxPlayer::OldFit(Int_t uid, Option_t *option) {

    TString opt = option;
    if ( !LoadChannel(uid) ) return kFALSE;
    
    // translate uid in buffer index
    Int_t idx = UidToIndex(uid);

    fH1DistTheta->Reset(); 
    fH1DistPhi->Reset(); 

    for(size_t j(0); j<fThetaFOVBuffer[idx].size(); j++) {
        fH1DistTheta->Fill(fThetaFOVBuffer[idx][j]);
        fH1DistPhi->Fill(fPhiFOVBuffer[idx][j]);
    }
    
    Bool_t fitTh = fH1DistTheta->GetEntries() > fFitThreshold;
    Bool_t fitPhi = fH1DistPhi->GetEntries() > fFitThreshold;
    if ( !fitTh && !opt.Contains("Q") ) 
        Printf("Fit: too few hits in theta. Skipping");

    if ( !fitPhi && !opt.Contains("Q") )
        Printf("Fit: too few hits in phi. Skipping");


    // get stats
    Float_t meanTh = fH1DistTheta->GetMean();
    Float_t rmsTh = fH1DistTheta->GetRMS();
    Float_t meanPhi = fH1DistPhi->GetMean();
    Float_t rmsPhi = fH1DistPhi->GetRMS();


    if ( fitTh ) {
        // fill and fit theta
        fH1FitTheta->GetXaxis()->SetLimits(meanTh-fNSigma*rmsTh, meanTh+fNSigma*rmsTh);
        fH1FitTheta->GetXaxis()->SetRange();
        fH1FitTheta->Reset();

        for(size_t j(0); j<fThetaFOVBuffer[idx].size(); j++)
            fH1FitTheta->Fill(fThetaFOVBuffer[idx][j]);

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

    if ( fitPhi ) {
        // fill and fit phi
        fH1FitPhi->GetXaxis()->SetLimits(meanPhi-fNSigma*rmsPhi, meanPhi+fNSigma*rmsPhi);
        fH1FitPhi->GetXaxis()->SetRange();
        fH1FitPhi->Reset();

        for(size_t j(0); j<fPhiFOVBuffer[idx].size(); j++)
            fH1FitPhi->Fill(fPhiFOVBuffer[idx][j]);

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

    return kTRUE;
}

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

    TString opt = option;
    if ( !LoadChannel(uid) ) return kFALSE;

    // clear histogramis
    fH2DistXY->Reset(); 
    fH1FitTheta->Reset();
    fH1FitPhi->Reset();
    
    // translate uid in buffer index and check number of hits
    Int_t idx = UidToIndex(uid);
    if ( (Int_t)fThetaFOVBuffer[idx].size() < fFitThreshold ) {
        if ( !opt.Contains("Q") ) Printf("Fit: %d entries for channel %d (idx = %d). Skipping",
                                         fThetaFOVBuffer[idx].size(), uid, idx);
        return kFALSE;
    }


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

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

        fH2DistXY->Fill(x,y);
    }

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

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

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

    // xmean,xxmean,ymean,yymean->th, therr, ph, phierr
    Float_t thCenter = TMath::ASin(r); 
    Float_t phiCenter = TMath::ATan2(yCenter, xCenter); 
    Float_t thErr = dr/TMath::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();

    // 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];
        switch ( quad ) {
            case 1:
                if ( fPhiFOVBuffer[idx][j] > 3*TMath::PiOver2() )
                    phi -= TMath::TwoPi();
                break;
            case 4:
                if ( fPhiFOVBuffer[idx][j] < TMath::PiOver2() )
                    phi += TMath::TwoPi();
                break;
            default:
                break;
        }
        
        fH1FitTheta->Fill(fThetaFOVBuffer[idx][j]);
        fH1FitPhi->Fill(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);

    return kTRUE;

}

//______________________________________________________________________________
 Bool_t OAPxPlayer::DrawChannel(Int_t uid, Option_t *option) {
    if ( !LoadChannel(uid) ) return kFALSE;
    
    // translate uid in buffer index
    Int_t idx = UidToIndex(uid);

    TString opt = option;
    TH1F *h, *hf;
    vector< vector<Float_t> >* pBuffer;
    if ( opt == "theta" ) {
        pBuffer = &fThetaFOVBuffer;
        h = fH1DistTheta;
        hf = fH1FitTheta;
    } else if ( opt == "phi" ) { 
        pBuffer = &fPhiFOVBuffer;
        h = fH1DistPhi;
        hf = fH1FitPhi;
    } else {
        Printf(opt+" unknown option");
        return kFALSE;
    }

    vector< vector<Float_t> >& buffer = *pBuffer;
    h->Reset();
    
    // fill histograms
    for(size_t k(0); k< buffer[idx].size(); k++)
        h->Fill(buffer[idx][k]);
    // get stats 
    Float_t mean = h->GetMean();
    Float_t rms = h->GetRMS();

    // resize histogram
    hf->GetXaxis()->SetLimits(mean-fNSigma*rms, mean+fNSigma*rms);
    hf->Reset();

    for(size_t k(0); k<buffer[idx].size(); k++)
        hf->Fill(buffer[idx][k]);

    // draw
    hf->Draw();

    return kTRUE;
}


//______________________________________________________________________________
 void OAPxPlayer::ExportMap(const char* filename) {
    ofstream fout;
    fout.open(filename);
    // loop on the pixelmap
//    EPixel *px;
    for(Int_t j(0); j<fNumPixels; j++){
//        px = (EPixel*)fPixelMap->At(j);
//        fout << Form("%d    %f    %f    %f    %f",
//                j+1, px->GetTheta(), px->GetThetaSigma(),
//                px->GetPhi(), px->GetPhiSigma()) << endl;
        
        fout << Form("%d    %f    %f    %f    %f",
                j+1, fTheta[j], fThetaRMS[j], fPhi[j], fPhiRMS[j]) << endl;
    }
    fout.close();
}


//______________________________________________________________________________
 void OAPxPlayer::MakeMap(Int_t first, Int_t last, const char* filename) {
    // loop on pixels and calls NewFit for each one.
    // if filemane is not null it
    Printf("MakeMap() called.");
    
    ofstream fout;
    if ( filename ) {
        Printf("Saving map in %s",filename);
        fout.open(filename);
        if (!fout.is_open())
            Printf("MakeMap: Unable to open %s", filename);
    }

    if ( first == 0 ) first = 1;
    if ( last == 0 ) last = fNumPixels;

    if ( first < 1 || last > fNumPixels || first > last) {
        Printf("Wrong indexes, exiting.");
        return;
    }
        
    TF1 *g;
    Double_t parTh[3], parPhi[3];

    // loop over pixel
    for(Int_t k(first); k<=last; k++){
        Int_t j=UidToIndex(k);
        // fit pixel data
        if (!Fit(k, "WEQ0")) {
            fTheta[j] = -1;
            fThetaRMS[j] = -1;
            fPhi[j] = -1;
            fPhiRMS[j] = -1;
        } 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
            fTheta[j] = parTh[1];
            fThetaRMS[j] = parTh[2];
            fPhi[j] = parPhi[1];
            fPhiRMS[j] = parPhi[2];
        }

        if (k % 100 == 0 ) Printf("ch %d done",k);
        
        if (fout.is_open())
            fout << Form("%d    %f    %f    %f    %f",
                k, fTheta[j], fThetaRMS[j], fPhi[j], fPhiRMS[j]) << endl;
    }
    
    if (fout.is_open()) fout.close();
}

//______________________________________________________________________________
 void OAPxPlayer::DrawHisto(Int_t hid) {
    // Draws one of the buffer's histos
    // hid = 1  oapx_theta
    //       2  oapx_thfit
    //       3  oapx_phi
    //       4  oapx_phifit 
    //       5  oapx_xy 

    switch (hid) {
        case 1:
            fH1DistTheta->Draw();;
            break;
        case 2:
            fH1FitTheta->Draw();
            break;
        case 3:
            fH1DistPhi->Draw();
            break;
        case 4:
            fH1FitPhi->Draw();
            break;
        case 5:
            fH2DistXY->Draw("colz");
            break;
        default:
            return;

    }

}

About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update 14-May-2005 21:31