Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

OpticsResponseBuilder - source file

// $Id: OpticsResponseBuilder.cc,v 1.6 2005/10/31 13:16:31 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>

using namespace sou;

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( fPsfZoom );
    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();
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005