Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

EOpticsResponse - source file

// $Id: EOpticsResponse.cc,v 1.11 2005/10/03 17:53:59 naumov Exp $
// Author: ale   2005/01/19

/*****************************************************************************
 * ESAF: Euso Simulation and Analysis Framework                              *
 *                                                                           *
 *  Id: EOpticsResponse                                                      *
 *  Package: <packagename>                                                   *
 *  Coordinator: <coordinator>                                               *
 *                                                                           *
 *****************************************************************************/

//_____________________________________________________________________________
//
// OpticsResponse
//
// <extensive class description>
//
//   Config file parameters
//   ======================
//
//   <parameter name>: <parameter description>
//   -Valid options: <available options>
//

#include "EOpticsResponse.hh"
#include "EORSample.hh"
#include "EORHeader.hh"

#include <TFile.h>
#include <TTree.h>
#include <TH2F.h>
#include <TClonesArray.h>
#include <stdexcept>
#include <iostream>
using namespace std;

ClassImp(EOpticsResponse)

//_____________________________________________________________________________
 EOpticsResponse::EOpticsResponse( const char *dbname ) : 
    fFile(0), fTree(0), fSample(0), fHeader(0), fPsf(0), fPoints(0) {
    //
    // Constructor
    //

    fFile = new TFile(dbname);

    if( fFile->IsZombie() ) MakeZombie();
    
    if( !(fTree = (TTree*)fFile->Get("ORDatabase")) ) MakeZombie();

    if( IsZombie() ) {
       Error("EOpticsResponse","Can not open or find the Optics DataBase.");
    
       Error("EOpticsResponse","Please run macros/tools/opticsresponse.C click Build button");
       
       throw runtime_error(Form("Cannot open file %s",dbname));
    }
    fHeader = (EORHeader*)fTree->GetUserInfo()->FindObject("EORHeader");

    fTree->SetBranchAddress("EORSample",&fSample);
    
    fTree->BuildIndex("EORSample.fThetaIndex","EORSample.fLambdaIndex");
       
    fPoints = new TClonesArray("TVector3",1000000);
    
    fPoints->BypassStreamer(kFALSE);
    
    Info("EOpticsResponse","Enabled. Reading from DB file: %s",fFile->GetName());
  
}

//_____________________________________________________________________________
 EOpticsResponse::~EOpticsResponse() {
    //
    // Destructor
    //

    SafeDelete(fFile);
    SafeDelete(fTree);
    SafeDelete(fSample);
    SafeDelete(fPoints);
}

//_____________________________________________________________________________
 Bool_t EOpticsResponse::GetEntry(Float_t theta, Float_t lambda) {
    //
    // Reads the entry from the OR DB file for given theta (in rad) and lambda (in nm)
    //

    if ( IsZombie() ) return kFALSE;

    Float_t ThetaStep  = ( fHeader->GetThetaMax() - fHeader->GetThetaMin() ) / fHeader->GetBinsTheta();
     
    Float_t LambdaStep = ( fHeader->GetLambdaMax() - fHeader->GetLambdaMin() ) / fHeader->GetBinsLambda(); 
    
    
    Int_t IndexTheta  = Int_t(theta/ThetaStep); 
    
    Int_t IndexLambda = Int_t((lambda-fHeader->GetLambdaMin())/LambdaStep); 
   
    if ( theta  <= fHeader->GetThetaMin()  )  IndexTheta  = fHeader->GetThetaIndexMin();
    if ( theta  >= fHeader->GetThetaMax()  )  IndexTheta  = fHeader->GetThetaIndexMax();
    if ( lambda <= fHeader->GetLambdaMin() )  IndexLambda = fHeader->GetLambdaIndexMin();
    if ( lambda >= fHeader->GetLambdaMax() )  IndexLambda = fHeader->GetLambdaIndexMax();
    
    fTree->GetEntryWithIndex(IndexTheta,IndexLambda);
    
    return kTRUE; 
}

//_____________________________________________________________________________
 TH2F* EOpticsResponse::MakePsfHist(Float_t theta, Float_t phi, TH1F* spectrum){
    //
    // returnes a histo with weighted sum of Psf distributions for the entry spectrum
    //
    if ( IsZombie() ) return NULL;
    if( !spectrum )               return (TH2F*) NULL;
    if( !spectrum->GetEntries() ) return (TH2F*) NULL;

    Float_t lambda, weigth, weigthtot(0);
    Int_t entry(0);

    for (Int_t i = 1; i <= spectrum->GetNbinsX(); i++) {
        lambda = spectrum->GetBinCenter(i);
        weigth = spectrum->GetBinContent(i);
        if (!weigth) continue;
        weigthtot += weigth; 
        GetEntry(theta,lambda);

        TH2F* htemp = GetSample()->GetPsfHist(); // TOCHECK if should be deleted in the end!!!

	if (!htemp) {
	    Warning("MakePsfHist(Float_t theta, Float_t phi, TH1F* spectrum)","can not find the histo");
	    return (TH2F*) NULL;
	}
        for (Int_t j=1; j<=htemp->GetNbinsX(); j++)
            for (Int_t k=1; k<=htemp->GetNbinsY(); k++) {
                if (entry > fPoints->GetEntriesFast()) break;

                Float_t x = htemp->GetXaxis()->GetBinCenter(j);
                Float_t y = htemp->GetYaxis()->GetBinCenter(k);
                Float_t z = htemp->GetBinContent(j,k);

                x = x*TMath::Cos(phi) - y*TMath::Sin(phi);
                y = x*TMath::Sin(phi) + y*TMath::Cos(phi);
                z = z*weigth;

                TClonesArray &ref = *fPoints;

                new(ref[entry++]) TVector3(x,y,z);
            } 
    }

    // define corners of the histo

    Float_t xmin(100000), xmax(-xmin), ymin(xmin), ymax(-ymin);
    for (Int_t i=0; i<entry; i++) {
        TVector3 *pos = (TVector3*)fPoints->At(i);    
	if (xmin > pos->X()) xmin = pos->X();
	if (xmax < pos->X()) xmax = pos->X();
	if (ymin > pos->Y()) ymin = pos->Y();
	if (ymax < pos->Y()) ymax = pos->Y();
    }

    Int_t nx,ny;

    nx = ny = fHeader->GetPsfBins();

    if ( fPsf ) {
        //TOCHECK: is it fast?

        fPsf->SetBins(nx, xmin, xmax, ny, ymin, ymax);
        fPsf->Reset();
    } 

    else {	

        fPsf = new TH2F("psf", "Psf for the spectrum", nx, xmin, xmax, ny, ymin, ymax);
//         fPsf->SetDirectory(0);
    }

    for (Int_t i=0; i < entry; i++) {

        TVector3 *pos = (TVector3*)fPoints->At(i);             
        fPsf->Fill(pos->X(),pos->Y(),pos->Z()/weigthtot);
    }
    return fPsf;
}

//______________________________________________________________________________
 void EOpticsResponse::Clear( Option_t* ) {
    //
    // Clears all temporary histograms
    //

    SafeDelete( fPsf );
    fPoints->Clear();

}

//_____________________________________________________________________________
 TH2F* EOpticsResponse::MakePsfHist(Float_t theta, Float_t phi, Float_t lambda){
    //
    // returnes Psf distribution for the given theta and phi
    //
    if ( !GetEntry(theta,lambda) ) return NULL; 
    TH2F* fPsf = fSample->MakePsfHistRotated(phi);
//     fPsf->SetDirectory(0);
    
    return fPsf;
    
}
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