GammaFilter.cpp

Go to the documentation of this file.

#ifndef ROOT_TFile
#include <TFile.h>
#endif

#ifndef GW_GAMMAFILTER_H
#include "GammaFilter.h"
#endif

using namespace Gw;

ClassImp(GammaFilter);

// GammaFilter members /////////////////////////////////////////////////////
GammaFilter::GammaFilter() : BaseFilter()
{
        fCompton = fEfficiency = fDoppler = NULL;

        // TMP a profite for the compton scattering
        // linear from 0 to 0.1 and then constant to 1
        fComptonShape = new TH1F("GammaFilter_ComptonShape","ComptonShape",1000,0,1);
        for (Int_t i = 0; i < 1000; i++ ) {
                if ( i < 200 ) fComptonShape->SetBinContent(i+1,0.005*i);
                else fComptonShape->SetBinContent(i+1,1);
        }
        
        SetNameTitle("GammaFilter","Standard Experimental Filter for Raw Gamma-rays");  
}

GammaFilter::GammaFilter(Axis_t xmin, Axis_t xmax, Float_t eff, Float_t PT, Float_t dop) 
{
        Axis_t lxmin, lxmax;
        fCompton = fEfficiency = fDoppler = NULL;
        
        // check without saying anything
        lxmin = xmin; lxmax = xmax;
        if ( lxmin > lxmax ) { lxmin = xmax; lxmax = xmin; }

        Int_t bin = Int_t ( lxmax - lxmin );

        fEfficiency     = new TH1F("GammaFilter_Efficiency","Efficiency",bin,lxmin,lxmax);
        fCompton        = new TH1F("GammaFilter_Compton","Compton",bin,lxmin,lxmax);
        fDoppler        = new TH1F("GammaFilter_Doppler","Doppler",bin,lxmin,lxmax);

        for (Int_t i = 1; i < bin ; i++ ) {
                fEfficiency->SetBinContent(i,TMath::Abs(eff));
                fCompton->SetBinContent(i,TMath::Abs(PT));
                fDoppler->SetBinContent(i,TMath::Abs(dop));
        }
        // IMPORTANT to set the underflow and overflow to 0 !
        fEfficiency->SetBinContent(0,0); fEfficiency->SetBinContent(bin,0); 
        
        SetNameTitle("GammaFilter","Standard Experimental Filter for Raw Gamma-rays");
}

GammaFilter::~GammaFilter()
{
        if ( fCompton != NULL )    { delete fCompton; fCompton = NULL; }
        if ( fEfficiency != NULL ) { delete fEfficiency; fEfficiency = NULL; }
        if ( fDoppler != NULL )    { delete fDoppler;  fDoppler = NULL; }
        if ( fComptonShape != NULL ) { delete fComptonShape; fComptonShape =  NULL; }
}

Bool_t GammaFilter::InitFilter(const char *rootfile)
{
        TH1F *h;
        
        // delete the previous definition
        if ( fCompton != NULL )    { delete fCompton; fCompton = NULL; }
        if ( fDoppler != NULL )    { delete fDoppler;  fDoppler = NULL; }
        if ( fEfficiency != NULL ) { delete fEfficiency; fEfficiency = NULL; }
                
        TFile *rfile = TFile::Open(rootfile);
        if ( rfile ) { // the file is open
                h = (TH1F *)rfile->Get("GammaFilter_Efficiency");
                if ( h ) { // the correct histogram has been found
                        fEfficiency = h; fEfficiency->SetDirectory(0);
                        fEfficiency->SetBinContent(0,0); 
                        fEfficiency->SetBinContent(fEfficiency->GetNbinsX()+1,0);       
                }
                h = (TH1F *)rfile->Get("GammaFilter_Compton");
                if ( h ) { // the correct histogram has been found
                        fCompton = h; fCompton->SetDirectory(0);
                }       
                h = (TH1F *)rfile->Get("GammaFilter_Doppler");
                if ( h ) { // the correct histogram has been found
                        fDoppler = h; fDoppler->SetDirectory(0);
                }               
        }
        delete rfile;
        
        return IsEffective();
}

Bool_t GammaFilter::IsEffective()
{
        return (fEfficiency != NULL) && (fCompton != NULL) && (fDoppler != NULL) ;
}
// void GammaFilter::SetEfficiency(TH1F &)
//{
// The first bin (0) and the highest bin (dim+1) must be equal to 0 to automatically reject 
// (see the IsRejected member) any gamma-ray that are out of the range defined by the histogram.   
//}

Bool_t GammaFilter::Compton (GammaLink *gam, Float_t &e)
{
        Bool_t ok; Float_t etmp;

        etmp = gam->GetEnergy().Get(); // to get the gamma-ray energy
        
        ok = fCompton->GetBinContent(fCompton->FindBin(etmp)) < fRand->Current()->Rndm();
        if ( ok ) e = fComptonShape->GetRandom() * etmp;

        return ok;
}

Int_t GammaFilter::ApplyE(const TSeqCollection &li_gamma, TBuffer &b)
{
        TObject *obj; GammaLink *gam;
        Float_t e;
        Int_t nb;
        
        nb = 0; TIter next(&li_gamma);
        while( (obj = next()) != NULL ) {
                gam = (GammaLink *)obj;
                
                // efficiency
                if ( IsRejected(gam) ) continue;

                // compton
                if ( Compton(gam,e) ) { nb++; b << e; continue; }

                // doppler broadring
                nb++; b << Doppler(gam);        
        }
        return nb;
}

Int_t GammaFilter::ApplyE(const TSeqCollection &li_gamma, Float_t *energies)
{
        TObject *obj; GammaLink *gam;
        Float_t e;
        Int_t nb;
        
        nb = 0; TIter next(&li_gamma);
        while( (obj = next()) != NULL ) {
                gam = (GammaLink *)obj;
                
                // efficiency
                if ( IsRejected(gam) ) continue;

                // compton
                if ( Compton(gam,e) ) { energies[nb++] = e; continue; }

                // doppler broadring
                energies[nb++] = Doppler(gam);  
        }
        return nb;
}
// End GammaFilter members /////////////////////////////////////////////////