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