source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/detector/tools/src/PhotonsDatabaseBuilder.cc @ 117

// ESAF : Euso Simulation and Analysis Framework
// $Id: PhotonsDatabaseBuilder.cc 2789 2008-02-11 07:26:05Z naumov $
// A.Thea, J.Watts created Mar, 19 2004

#include "PhotonsDatabaseBuilder.hh"
#include "Photon.hh"
#include "OpticalSystem.hh"
#include "OpticsFactory.hh"
#include "TTimeStamp.h"
#include <TMath.h>

ClassImp(PhotonsDatabaseBuilder)

//______________________________________________________________________________
PhotonsDatabaseBuilder::PhotonsDatabaseBuilder() {
    // ctor
    
    fSamples = (Int_t)Conf()->GetNum("PhotonsDatabaseBuilder.fSamples");
    fOptics = OpticsFactory::Get()->MakeOpticalSystem();
    fRootFileName = Conf()->GetStr("PhotonsDatabaseBuilder.fRootFileName");
}

//______________________________________________________________________________
PhotonsDatabaseBuilder::~PhotonsDatabaseBuilder() {
    // dtor

    delete fOptics;
    
}

//______________________________________________________________________________
void PhotonsDatabaseBuilder::CloseRoot() { 
    // close root file
    if ( fRootFile ) {
        fRootFile->Write();
        fRootFile->Close();
        cout << "Root file closed\n";
        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 opened\n";
    }

}
//______________________________________________________________________________
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.GetPhoton() ) {
                // 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();
}