// 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;
}
}