// $Id: OpticsResponseBuilder.cc,v 1.3 2005/03/17 12:42:46 naumov Exp $
// Author: ale 2005/01/19
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: OpticsResponseBuilder *
* Package: <packagename> *
* Coordinator: <coordinator> *
* *
*****************************************************************************/
//_____________________________________________________________________________
//
// OpticsResponseBuilder
//
// <extensive class description>
//
// Config file parameters
// ======================
//
// <parameter name>: <parameter description>
// -Valid options: <available options>
//
#include "OpticsResponseBuilder.hh"
#include "ParentPhoton.hh"
#include "Photon.hh"
#include "EORSample.hh"
#include "EORHeader.hh"
#include "EsafRandom.hh"
#include "OpticalSystem.hh"
#include "IdealFocalSurface.hh"
#include <TFile.h>
#include <TStyle.h>
#include <TH1.h>
#include <TH2F.h>
#include <TVector3.h>
#include <TGraph.h>
#include <TMultiGraph.h>
#include <TCanvas.h>
#include <TGaxis.h>
#include <TLegend.h>
#include <TRandom.h>
#include <TROOT.h>
#include <TMath.h>
#include <TTree.h>
#include <TBenchmark.h>
#include <TClonesArray.h>
ClassImp(OpticsResponseBuilder)
//_____________________________________________________________________________
OpticsResponseBuilder::OpticsResponseBuilder() : fDBName(""), fDBFile(0),
fDBTree(0), fOpticalSystem(0), fIdealFocalSurface(0), fPsf(0), fPsfZoom(0),
fPsfRProj(0), fCentroid(0.,0.,0.), fPoints(0), fParent(0), fSample(0) {
//
// Constructor
//
// default parameters
fDiscRadius = 1500.*mm;
fEusoRadius = 1250.*mm;
fRmax = 30.;
fTriggDiam = 5.;
// default settings
fThetaMin = 0*deg;
fThetaMax = 30*deg;
fLambdaMin = 300*nm;
fLambdaMax = 400*nm;
fBinsTheta = 30;
fBinsLambda = 100;
fPsfBins = 100;
fPsfZoomSide = 15.*mm;
fPoints = new TClonesArray("TVector3");
SetNrays(10000);
fParent = new ParentPhoton;
}
//_____________________________________________________________________________
OpticsResponseBuilder::~OpticsResponseBuilder() {
//
// Destructor
//
SafeDelete( fDBFile );
SafeDelete( fPsf );
SafeDelete( fPsfRProj );
SafeDelete( fPoints );
SafeDelete( fParent );
SafeDelete( fSample );
}
//_____________________________________________________________________________
void OpticsResponseBuilder::SetNrays(Int_t n){
//
// Reset clear fPoints and reset its size
//
fNrays = n;
fPoints->Clear();
fPoints->Expand(n);
fPoints->BypassStreamer(kFALSE);
}
//______________________________________________________________________________
void OpticsResponseBuilder::Clear() {
//
// Clears all temporary histograms
//
SafeDelete( fPsf );
SafeDelete( fPsfZoom );
SafeDelete( fPsfRProj );
fPoints->Clear();
}
//______________________________________________________________________________
void OpticsResponseBuilder::OpenDB() {
//
// Helper function. Opens and initializes the database
//
if ( fDBName.IsNull() ) {
MsgForm(EsafMsg::Warning, "Database file name is empty. File not opened.");
return;
}
fDBFile = new TFile(fDBName,"RECREATE");
if ( fDBFile->IsZombie() ) {
MsgForm(EsafMsg::Warning, "File %s not opened.",fDBName.Data());
return;
}
EORHeader *header = new EORHeader;
// instance database objects
fSample = new EORSample;
fDBTree = new TTree("ORDatabase","ORDatabase");
fDBTree->Branch("EORSample","EORSample",&fSample);
// fill the header
// simulation parameters
header->SetOpticsName(fOpticalSystem->IsA()->GetName());
header->SetIdealSurfaceName(fIdealFocalSurface->IsA()->GetName());
header->SetDiscRadius( fDiscRadius );
header->SetEusoRadius( fEusoRadius );
header->SetRmax( fRmax );
header->SetTriggDiam( fTriggDiam );
// run parameters
header->SetBins(fBinsTheta,fBinsLambda);
header->SetNrays( fNrays );
// header will be deleted by TTree destructor
fDBTree->GetUserInfo()->Add(header);
MsgForm(EsafMsg::Info,"Database file %s opened.",fDBName.Data());
}
//______________________________________________________________________________
void OpticsResponseBuilder::CloseDB() {
//
//
//
if ( !fDBFile ) {
MsgForm(EsafMsg::Warning,"No database file opened.");
return;
}
if ( fDBFile != fDBTree->GetCurrentFile() ) {
MsgForm(EsafMsg::Warning,"Mismatch between database file and current TTree file.");
Msg(EsafMsg::Warning) << "Maybe a file change has occurred. ";
Msg(EsafMsg::Warning) << "Only the latter will be saved." << MsgDispatch;
fDBFile = fDBTree->GetCurrentFile();
}
fDBTree->Write();
SafeDelete( fDBTree );
SafeDelete( fDBFile );
SafeDelete( fSample );
}
//_____________________________________________________________________________
Bool_t OpticsResponseBuilder::ComputeEfficacy(Float_t theta , Float_t phi, Float_t lambda) {
//
// Simulate and collect a uniform beam of photons for the chosen optical
// system.
//
if ( !fOpticalSystem ) {
Msg(EsafMsg::Warning) << "No OpticalSystem defined. Exiting."<< MsgDispatch;
return kFALSE;
}
if ( !fIdealFocalSurface ) {
Msg(EsafMsg::Warning) << "No IdealFocalSurface defined. Exiting."<< MsgDispatch;
return kFALSE;
}
Double_t pos[3] = {0,0,0};
Double_t spos[3] = {0,0,0};
TRandom* rndm = EsafRandom::Get();
Int_t ntrans(0);
Photon* ph;
for( Int_t i(0); i<fNrays; i++) {
// new photon position
Double_t radius = sqrt(rndm->Rndm())*fDiscRadius;
Double_t alpha = rndm->Rndm()*TMath::TwoPi();
pos[0] = radius*TMath::Cos(alpha);
pos[1] = radius*TMath::Sin(alpha);
// update the fParent
fParent->Build(0, spos, pos, theta, phi, lambda, 0, Fluorescence);
ph = &(fParent->GetPhoton());
ph = fOpticalSystem->Transport(ph);
// reject dead photons and going downward
if (!ph || ph->dir.Z() < 0.) continue;
fIdealFocalSurface->HitPosition(ph);
if (!ph->hitIfs) continue;
TClonesArray &ref = *fPoints;
new(ref[ntrans++]) TVector3(ph->posOnIfs);
}
Float_t pflux = fNrays/(TMath::Pi()*fDiscRadius*fDiscRadius*TMath::Cos(theta));
fTotalEfficacy = ntrans/pflux;
if ( fSample ) {
// fill the sample with the current results
fSample->SetTheta( theta );
fSample->SetLambda( lambda );
fSample->SetTotalEfficacy( fTotalEfficacy );
}
FillPsfHist();
if ( !fPsfRProj ) {
MsgForm(EsafMsg::Warning,"No photons in the point spread function");
return kTRUE;
}
Int_t bin = fPsfRProj->GetXaxis()->FindBin(fTriggDiam/2.);
fTriggEfficacy = fPsfRProj->Integral(1,bin)/pflux;
fPsf->Scale(1./pflux);
fPsfZoom->Scale(1./pflux);
if ( fSample ) {
fSample->SetTriggEfficacy( fTriggEfficacy );
fSample->SetCentroid( fCentroid );
fSample->SetPsfHist( fPsfZoom );
}
return kTRUE;
}
//______________________________________________________________________________
void OpticsResponseBuilder::FillPsfHist() {
//
// Build histograms
//
Float_t xmin, ymin;
Float_t xmax, ymax;
xmin = ymin =-fEusoRadius;
xmax = ymax = fEusoRadius;
Int_t npos = fPoints->GetEntriesFast();
if ( fPsf ) {
//TOCHECK: is it fast?
fPsf->SetBins(20*fPsfBins, xmin, xmax, 20*fPsfBins, ymin, ymax);
fPsf->Reset();
} else {
fPsf = new TH2F("psf", "Psf", 20*fPsfBins, xmin, xmax, 20*fPsfBins, ymin, ymax);
fPsf->SetDirectory(0);
}
// search for the peak
for (Int_t i=0; i< npos; i++) {
TVector3 *pos = (TVector3*)fPoints->At(i);
fPsf->Fill(pos->X(),pos->Y());
}
Int_t bx, by, bz;
fPsf->GetMaximumBin(bx, by, bz);
Double_t peakX = fPsf->GetXaxis()->GetBinCenter( bx );
Double_t peakY = fPsf->GetYaxis()->GetBinCenter( by );
TVector3 fAverage(0,0,0);
Int_t nnear(0);
Float_t spacethreshold = 20.*mm;
// find the center of mass around the peak
for (Int_t i=0; i< npos; i++ ){
TVector3 *pos = (TVector3*)fPoints->At(i);
Double_t dist = TMath::Sqrt(TMath::Power(pos->X()-peakX,2)
+TMath::Power(pos->Y()-peakY,2));
if ( dist > spacethreshold) continue;
fAverage += *pos;
nnear++;
}
if (nnear == 0)
return;
fAverage *= 1./nnear;
xmin = fAverage.X()-fPsfZoomSide;
xmax = fAverage.X()+fPsfZoomSide;
ymin = fAverage.Y()-fPsfZoomSide;
ymax = fAverage.Y()+fPsfZoomSide;
// build the projection
if ( fPsfRProj ) {
fPsfRProj->SetBins( fPsfBins, 0, fPsfZoomSide );
fPsfRProj->Reset();
} else {
fPsfRProj = new TH1F("PsfRProj", "Psf R-projection", fPsfBins, 0, fPsfZoomSide);
fPsfRProj->SetDirectory(0);
}
// build the zoom
if ( fPsfZoom ) {
fPsfZoom->SetBins( fPsfBins, xmin, xmax, fPsfBins, ymin, ymax);
fPsfZoom->Reset();
} else {
fPsfZoom = new TH2F("PsfZoom", "Psf Zoom", fPsfBins, xmin, xmax, fPsfBins, ymin, ymax);
fPsfZoom->SetDirectory(0);
}
fCentroid = fAverage;
Float_t r;
for (Int_t i=0; i< nnear; i++) {
TVector3 *pos = (TVector3*)fPoints->At(i);
fPsfZoom->Fill(pos->X(), pos->Y());
r = (*pos - fAverage).Mag();
if( r < fRmax)
fPsfRProj->Fill( r );
}
fPsf->GetXaxis()->SetTitle("mm");
fPsf->GetYaxis()->SetTitle("mm");
fPsf->GetYaxis()->SetTitleOffset(1.2);
fPsfRProj->GetXaxis()->SetTitle("mm");
}
//______________________________________________________________________________
void OpticsResponseBuilder::BuildEfficacyTables() {
//
//
//
OpenDB();
Float_t th, wl;
Float_t thstep = (fThetaMax-fThetaMin)/fBinsTheta;
Float_t wlstep = (fLambdaMax-fLambdaMin)/fBinsLambda;
if ( fDBTree ) {
EORHeader* header = (EORHeader*)fDBTree->GetUserInfo()->FindObject("EORHeader");
if ( header ) {
header->SetThetaMin(fThetaMin);
header->SetThetaMax(fThetaMax);
header->SetThetaIndexMin(0);
header->SetThetaIndexMax(fBinsTheta);
header->SetLambdaMin(fLambdaMin);
header->SetLambdaMax(fLambdaMax);
header->SetLambdaIndexMin(0);
header->SetLambdaIndexMax(fBinsLambda);
header->SetPsfBins(fPsfBins);
}
}
for (Int_t i(0); i <=fBinsTheta; i++) {
th = fThetaMin + thstep*i;
Msg(EsafMsg::Info) << Form("Theta: %.2f deg ", th/deg) << MsgFlush;
for (Int_t j(0); j <=fBinsLambda; j++) {
wl = fLambdaMin + wlstep*j;
if ( (j+1)*10%fBinsLambda == 0 )
Msg(EsafMsg::Info) << " ... "<<(j+1)*100/fBinsLambda << '%'<< MsgFlush;
ComputeEfficacy(th, 0. ,wl);
if ( fSample ) {
fSample->SetThetaIndex( i );
fSample->SetLambdaIndex( j );
}
if ( fDBTree ) fDBTree->Fill();
if ( fSample ) fSample->Clear();
Clear();
}
Msg(EsafMsg::Info) << " -done. " << MsgDispatch;
}
CloseDB();
}