// ESAF : Euso Simulation and Analysis Framework
// $Id: PhotonsDatabaseBuilder.cc,v 1.9 2004/06/16 13:35:18 thea Exp $
// A.Thea, J.Watts created Mar, 19 2004
#include "PhotonsDatabaseBuilder.hh"
#include "Photon.hh"
#include "OpticalSystem.hh"
#include "OpticsFactory.hh"
#include "TTimeStamp.h"
ClassImp(PhotonsDatabaseBuilder)
//______________________________________________________________________________
PhotonsDatabaseBuilder::PhotonsDatabaseBuilder() {
// ctor
fSamples = (Int_t)Conf()->GetNum("PhotonsDatabaseBuilder.fSamples");
fOptics = OpticsFactory::Get()->GetOpticalSystem();
fRootFileName = Conf()->GetStr("PhotonsDatabaseBuilder.fRootFileName");
}
//______________________________________________________________________________
PhotonsDatabaseBuilder::~PhotonsDatabaseBuilder() {
// dtor
}
//______________________________________________________________________________
void PhotonsDatabaseBuilder::CloseRoot() {
// close root file
if ( fRootFile ) {
fRootFile->Write();
fRootFile->Close();
cout << "Root file closedn";
delete fRootFile; fRootFile = 0;
}
}
//______________________________________________________________________________
void PhotonsDatabaseBuilder::OpenRoot() {
// open rootfile with the correct name
CloseRoot();
cout << "Opening root file " << fRootFileName << endl;
fRootFile = new TFile((fRootFileName+".root").c_str(),"CREATE");
if (fRootFile->IsZombie()) {
cout << endl;
cout << "Error opening file " << fRootFileName << endl;
cout << "Probably it already exists or there is no " << endl;
cout << "write access permission to this directory" << endl;
cout << endl;
// try adding date and time to the name
TTimeStamp time;
TString now = time.AsString("s");
now[now.Last(' ')]=':';
TString newFileName = Form("%s-", fRootFileName.c_str())+now+".root";
fRootFile = new TFile(newFileName,"CREATE");
if (fRootFile->IsZombie()) {
throw invalid_argument("Cannot open root file");
}
cout << "Root file " << newFileName << " successfully openedn";
}
}
//______________________________________________________________________________
void PhotonsDatabaseBuilder::InitDBTree() {
// create the tree
fTree = new OADBTree;
fPhotons = new OADBPhotons();
cout << fPupil.GetNumPos() << " spots to scan" << endl;
fPhotons->SetSize(fPupil.GetNumPos());
fTree->Branch("photons","OADBPhotons",&fPhotons);
// fill the header
OADBHeader &header = fTree->GetHeader();
header.SetXYStep( fPupil.GetXYStep() );
header.SetOpticsRadius( fOptics->Radius() );
header.SetEntranceDiscRadius( fPupil.GetEntranceDiscRadius() );
header.SetFirstLensBottom( fOptics->FirstLensBottom() );
header.SetFirstLensTop( fOptics->FirstLensTop() );
header.SetSecondLensBottom( fOptics->SecondLensBottom() );
header.SetSecondLensTop( fOptics->SecondLensTop() );
header.SetEntranceZ(fOptics->Bottom());
header.SetExitZ(fOptics->Top());
header.SetPosZ(fOptics->Position().Z());
vector<Double_t> dummyD;
dummyD = fPupil.GetAngles();
header.SetNumAngles(dummyD.size());
for(size_t i(0); i<dummyD.size(); i++ )
header.SetAngle(i, dummyD[i]);
cout << header.GetNumAngles() << " angles saved in the header" << endl;
dummyD = fPupil.GetWavelengths();
header.SetNumWavelengths(dummyD.size());
for(size_t i(0); i<dummyD.size(); i++ )
header.SetWavelength(i, dummyD[i]);
vector<Int_t> dummyI = fPupil.GetRows();
header.SetNumRows(dummyI.size());
for(size_t i(0); i<dummyI.size(); i++ )
header.SetNumCols(i, dummyI[i]);
}
//______________________________________________________________________________
void PhotonsDatabaseBuilder::FillDBTree() {
// calls fTree->Fill() and then clears the photon container
fTree->Fill();
fPhotons->Clean();
// check rootfile
if ( fRootFile != fTree->GetCurrentFile() ) {
fRootFile = fTree->GetCurrentFile();
cout << Form("RootFile capacity limit reached. Moving to %s",
fRootFile->GetName()) << endl;
}
}
//______________________________________________________________________________
void PhotonsDatabaseBuilder::CloseDBTree() {
// write and close the tree
cout << "Saving rootfile" << endl;
fTree->Write();
delete fTree; fTree = 0;
delete fPhotons; fPhotons = 0;
}
//______________________________________________________________________________
Bool_t PhotonsDatabaseBuilder::Sample( Photon *ph ) {
// resample the same photon and save the most rated path in fLastPath
vector<pair<Photon,Int_t> > buffer;
// fire fSamples Photons
for(Int_t k(0); k< fSamples; k++){
// reset the photon
ph->parent->CreatePhoton();
// transport the ph
Photon* newph = fOptics->Transport(ph);
// if survives and it is out of the second lens add to the list
// FIXME it does not take in account the phs that go out from the sides
if ( newph && TMath::Abs(newph->pos[Z]-fOptics->SecondLensTop()) < TOLERANCE) {
vector< pair<Photon, Int_t> >::iterator it;
Bool_t isNew(kTRUE);
// look for a similar ph in the buffer
for(it=buffer.begin(); it!=buffer.end(); it++){
Float_t dirDiff = (it->first.dir.Unit()-newph->dir.Unit()).Mag();
Float_t posDiff = (it->first.pos-newph->pos).Mag();
if (dirDiff < TOLERANCE && posDiff < TOLERANCE ){
isNew = kFALSE;
it->second++;
break;
}
}
// this ph is new, add an entry
if (isNew) {
pair<Photon,Int_t> dummy(*newph, 1);
buffer.push_back(dummy);
}
}
}
vector< pair<Photon, Int_t> >::iterator it;
fLastPath.second = 0;
for (it = buffer.begin(); it != buffer.end(); it++) {
if (it->second > fLastPath.second) {
// calculate transmission and put it in the path
fLastPath.first = it->first;
fLastPath.second = it->second/(Float_t)fSamples;
}
}
if ( fLastPath.second == 0 )
return kFALSE;
else
return kTRUE;
}
//______________________________________________________________________________
void PhotonsDatabaseBuilder::Go() {
// build the database
OpenRoot();
InitDBTree();
for(; !fPupil.IsWlEnd() ;fPupil.NextWl() ) {
cout << "Switching to wavelength " << fPupil.GetCurWl() << endl;
for(; !fPupil.IsAngleEnd() ;fPupil.NextAngle() ) {
cout << "Switching to angle " << TMath::RadToDeg()*fPupil.GetCurAngle() << endl;
fPhotons->SetAngle(fPupil.GetCurAngle(), fPupil.GetCurAngleId());
fPhotons->SetWavelength(fPupil.GetCurWl(), fPupil.GetCurWlId());
while ( Photon *ph = fPupil.Get() ) {
// find the "top rated photon"
if ( Sample( ph ) ) {
// get the "top rated photon" and save it
fPhotons->SavePhoton(fPupil.GetLastPosId(), &(fLastPath.first), fLastPath.second );
} else {
//if there are no ph in the database, save a ghost
fPhotons->SaveGhost(fPupil.GetLastPosId());
}
}
FillDBTree();
}
}
CloseDBTree();
CloseRoot();
}