AncillaryFilter.C

Go to the documentation of this file.

#define AncillaryFilter_cxx
// The class definition in AncillaryFilter.h has been generated automatically
// by the ROOT utility TTree::MakeSelector(). This class is derived
// from the ROOT class TSelector. For more information on the TSelector
// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.

// The following methods are defined in this file:
//    Begin():        called every time a loop on the tree starts,
//                    a convenient place to create your histograms.
//    SlaveBegin():   called after Begin(), when on PROOF called only on the
//                    slave servers.
//    Process():      called for each event, in this function you decide what
//                    to read and fill your histograms.
//    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
//                    called only on the slave servers.
//    Terminate():    called at the end of the loop on the tree,
//                    a convenient place to draw/fit your histograms.
//
// To use this file, try the following session on your Tree T:
//
// Root > T->Process("AncillaryFilter.C")
// Root > T->Process("AncillaryFilter.C","some options")
// Root > T->Process("AncillaryFilter.C+")
//

#include "AncillaryFilter.h"
#include <TH2.h>
#include <TStyle.h>
#include "TSystem.h"
#include "TEventList.h"
#include "TEntryList.h"
#include "TH1.h"
#include "TH2.h"
#include <TDatabasePDG.h>

#include <iostream>
#include <iomanip>
#include <fstream>


namespace {
Int_t gPDG;
Int_t nPDG;
Int_t pPDG;
Int_t aPDG;

// event or entry list
TEventList *evlist; TEntryList *enlist;
}

void AncillaryFilter::Begin(TTree * tree)
{
    // The Begin() function is called at the start of the query.
    // When running with PROOF Begin() is only called on the client.
    // The tree argument is deprecated (on PROOF 0 is passed).

    Init(tree);

    ReadConfFile();

    AGATACountsGlobNorm = new TH2F("AGATACountsGlobNorm","AGATACountsGlobNorm",50,0,50,50,0,50);
    AGATACountsGlobNorm->GetXaxis()->SetTitle("Number of emited gammas");
    AGATACountsGlobNorm->GetYaxis()->SetTitle("Number of hits in AGATA");
    AGATACountsGlobNorm->GetXaxis()->SetTitleOffset(0.9);
    AGATACountsGlobNorm->GetXaxis()->SetTitleSize(0.05);
    AGATACountsGlobNorm->GetYaxis()->SetTitleOffset(0.8);
    AGATACountsGlobNorm->GetYaxis()->SetTitleSize(0.05);

    DiamantCountsGlobNorm = new TH2F("DiamantCountsGlobNorm","DiamantCountsGlobNorm",10,0,10,20,0,20);
    DiamantCountsGlobNorm->GetXaxis()->SetTitle("Number of emited charged particles");
    DiamantCountsGlobNorm->GetYaxis()->SetTitle("Number of hits in DIAMANT");
    DiamantCountsGlobNorm->GetXaxis()->SetTitleOffset(0.9);
    DiamantCountsGlobNorm->GetXaxis()->SetTitleSize(0.05);
    DiamantCountsGlobNorm->GetYaxis()->SetTitleOffset(0.8);
    DiamantCountsGlobNorm->GetYaxis()->SetTitleSize(0.05);

    NWCountsGlobNorm = new TH2F("NWCountsGlobNorm","NWCountsGlobNorm",10,0,10,20,0,20);
    NWCountsGlobNorm->GetXaxis()->SetTitle("Number of emited neutrons");
    NWCountsGlobNorm->GetYaxis()->SetTitle("Number of hits in NW");
    NWCountsGlobNorm->GetXaxis()->SetTitleOffset(0.9);
    NWCountsGlobNorm->GetXaxis()->SetTitleSize(0.05);
    NWCountsGlobNorm->GetYaxis()->SetTitleOffset(0.8);
    NWCountsGlobNorm->GetYaxis()->SetTitleSize(0.05);

    AGATACountsPrNorm = new TH2F("AGATACountsPrNorm","AGATACountsPrNorm",50,0,50,50,0,50);
    AGATACountsPrNorm->GetXaxis()->SetTitle("Number of emited gammas");
    AGATACountsPrNorm->GetYaxis()->SetTitle("Number of hits in AGATA");
    AGATACountsPrNorm->GetXaxis()->SetTitleOffset(0.9);
    AGATACountsPrNorm->GetXaxis()->SetTitleSize(0.05);
    AGATACountsPrNorm->GetYaxis()->SetTitleOffset(0.8);
    AGATACountsPrNorm->GetYaxis()->SetTitleSize(0.05);

    DiamantCountsPrNorm = new TH2F("DiamantCountsPrNorm","DiamantCountsPrNorm",10,0,10,20,0,20);
    DiamantCountsPrNorm->GetXaxis()->SetTitle("Number of emited charged particles");
    DiamantCountsPrNorm->GetYaxis()->SetTitle("Number of hits in DIAMANT");
    DiamantCountsPrNorm->GetXaxis()->SetTitleOffset(0.9);
    DiamantCountsPrNorm->GetXaxis()->SetTitleSize(0.05);
    DiamantCountsPrNorm->GetYaxis()->SetTitleOffset(0.8);
    DiamantCountsPrNorm->GetYaxis()->SetTitleSize(0.05);

    NWCountsPrNorm = new TH2F("NWCountsPrNorm","NWCountsPrNorm",10,0,10,20,0,20);
    NWCountsPrNorm->GetXaxis()->SetTitle("Number of emited neutrons");
    NWCountsPrNorm->GetYaxis()->SetTitle("Number of hits in NW");
    NWCountsPrNorm->GetXaxis()->SetTitleOffset(0.9);
    NWCountsPrNorm->GetXaxis()->SetTitleSize(0.05);
    NWCountsPrNorm->GetYaxis()->SetTitleOffset(0.8);
    NWCountsPrNorm->GetYaxis()->SetTitleSize(0.05);

    AGATACountsEvNorm = new TH2F("AGATACountsEvNorm","AGATACountsEvNorm",50,0,50,50,0,50);
    AGATACountsEvNorm->GetXaxis()->SetTitle("Number of emited gammas");
    AGATACountsEvNorm->GetYaxis()->SetTitle("Number of hits in AGATA");
    AGATACountsEvNorm->GetXaxis()->SetTitleOffset(0.9);
    AGATACountsEvNorm->GetXaxis()->SetTitleSize(0.05);
    AGATACountsEvNorm->GetYaxis()->SetTitleOffset(0.8);
    AGATACountsEvNorm->GetYaxis()->SetTitleSize(0.05);

    DiamantCountsEvNorm = new TH2F("DiamantCountsEvNorm","DiamantCountsEvNorm",10,0,10,20,0,20);
    DiamantCountsEvNorm->GetXaxis()->SetTitle("Number of emited charged particles");
    DiamantCountsEvNorm->GetYaxis()->SetTitle("Number of hits in DIAMANT");
    DiamantCountsEvNorm->GetXaxis()->SetTitleOffset(0.9);
    DiamantCountsEvNorm->GetXaxis()->SetTitleSize(0.05);
    DiamantCountsEvNorm->GetYaxis()->SetTitleOffset(0.8);
    DiamantCountsEvNorm->GetYaxis()->SetTitleSize(0.05);

    NWCountsEvNorm = new TH2F("NWCountsEvNorm","NWCountsEvNorm",10,0,10,20,0,20);
    NWCountsEvNorm->GetXaxis()->SetTitle("Number of emited neutrons");
    NWCountsEvNorm->GetYaxis()->SetTitle("Number of hits in NW");
    NWCountsEvNorm->GetXaxis()->SetTitleOffset(0.9);
    NWCountsEvNorm->GetXaxis()->SetTitleSize(0.05);
    NWCountsEvNorm->GetYaxis()->SetTitleOffset(0.8);
    NWCountsEvNorm->GetYaxis()->SetTitleSize(0.05);

    GlobalExitChannelDistrib = new TH1F("ExitChannelDistrib","ExitChannelDistrib",50,0,50);

    for(int i=0 ; i<=5 ; i++)
    {
        for(int j=0 ; j<=5 ; j++)
        {
            for(int k=0 ; k<=1 ; k++)
            {
                TString Name = Form("ECDist_%dn%dp%da",i,j,k);
                ExitChannelDistrib[i][j][k] = new TH1F(Name,Name,50,0,50);
            }
        }
    }

    NEntries = tree->GetEntries();
    Fact = NEntries/100;

    gPDG = TDatabasePDG::Instance()->GetParticle("gamma")->PdgCode();
    nPDG = TDatabasePDG::Instance()->GetParticle("neutron")->PdgCode();
    pPDG = TDatabasePDG::Instance()->GetParticle("proton")->PdgCode();
    aPDG = 1000020040;

    //ELIST
    evlist = 0x0;
    enlist = 0x0;

    if ( gDirectory->GetList()->FindObject("SToGS_entrylist") ) {
        delete gDirectory->GetList()->FindObject("SToGS_entrylist");
    }
    fChain->SetEntryList(0);

    // case when one creates/fills the event list
    enlist = new TEntryList("SToGS_entrylist","selection from TSelector");
    enlist->SetDirectory(0);

    enlist->SetTreeName(fChain->GetName());
    enlist->SetFileName(gSystem->BaseName(fChain->GetDirectory()->GetName()));

    cout<<endl;
    cout<<"Starting to process events : "<<endl;
    cout<<endl;
}

void AncillaryFilter::SlaveBegin(TTree * /*tree*/)
{
    // The SlaveBegin() function is called after the Begin() function.
    // When running with PROOF SlaveBegin() is called on each slave server.
    // The tree argument is deprecated (on PROOF 0 is passed).

    TString option = GetOption();

}

Bool_t AncillaryFilter::Process(Long64_t entry)
{
    // The Process() function is called for each entry in the tree (or possibly
    // keyed object in the case of PROOF) to be processed. The entry argument
    // specifies which entry in the currently loaded tree is to be processed.
    // It can be passed to either AncillaryFilter::GetEntry() or TBranch::GetEntry()
    // to read either all or the required parts of the data. When processing
    // keyed objects with PROOF, the object is already loaded and is available
    // via the fObject pointer.
    //
    // This function should contain the "body" of the analysis. It can contain
    // simple or elaborate selection criteria, run algorithms on the data
    // of the event and typically fill histograms.
    //
    // The processing can be stopped by calling Abort().
    //
    // Use fStatus to set the return value of TTree::Process().
    //
    // The return value is currently not used.

    GetEntry(entry);

    if(fChain->GetChainEntryNumber(entry+1)%Fact==0)
    {
        cout<<"\r"<<"Analysis progress : "<<setw(3)<<(double)fChain->GetChainEntryNumber(entry+1)/(double)NEntries*100<<" %"<<" ("<<fChain->GetChainEntryNumber(entry+1)<<"/"<<NEntries<<")"<<flush;
    }

    Int_t g_pr_hits, n_pr_hits, p_pr_hits, a_pr_hits, nb_AGATA, nb_NW, nb_DIAM, nb_DIAM_proton, nb_DIAM_alpha;

    g_pr_hits=0;n_pr_hits=0;p_pr_hits=0;a_pr_hits=0;nb_AGATA=0;nb_NW=0;nb_DIAM=0;nb_DIAM_proton=0;nb_DIAM_alpha=0;


    for (int i=0; i< Pr_fHits_; i++)
    {
        if (Pr_fHits_fPDG[i]  == gPDG)
            g_pr_hits++;
        else if (Pr_fHits_fPDG[i]  == nPDG)
        {
            n_pr_hits++;
        }
        else if (Pr_fHits_fPDG[i]  == pPDG)
            p_pr_hits++;

        else if (Pr_fHits_fPDG[i]  == aPDG)
            a_pr_hits++;
    }
    for (int i=0; i< Ev_fHits_; i++)
    {
        if ( Ev_fHits_fUID[i] >=0 && Ev_fHits_fUID[i] <=180 ) // AGATA, compute number of hits in the array [tracking mode]
            nb_AGATA++;
        else if (Ev_fHits_fUID[i]>=1000 && Ev_fHits_fUID[i] <2000)
        {
            nb_DIAM++;
            if(Pr_fHits_fPDG[Ev_fHits_fFlag[i]] == pPDG) nb_DIAM_proton++;
            else if(Pr_fHits_fPDG[Ev_fHits_fFlag[i]] == aPDG) nb_DIAM_alpha++;
        }
        else if (Ev_fHits_fUID[i]>=2000 && Ev_fHits_fUID[i] <3000)
            nb_NW++;
    }

    if(PlotEfficiencies == 2 && (nb_AGATA >= fN_Agata[0] && nb_AGATA <= fN_Agata[1]))
    {
        AGATACountsPrNorm->Fill(g_pr_hits,nb_AGATA);
        DiamantCountsPrNorm->Fill(p_pr_hits+a_pr_hits,nb_DIAM);
        NWCountsPrNorm->Fill(n_pr_hits,nb_NW);

        AGATACountsEvNorm->Fill(g_pr_hits,nb_AGATA);
        DiamantCountsEvNorm->Fill(p_pr_hits+a_pr_hits,nb_DIAM);
        NWCountsEvNorm->Fill(n_pr_hits,nb_NW);

        AGATACountsGlobNorm->Fill(g_pr_hits,nb_AGATA);
        DiamantCountsGlobNorm->Fill(p_pr_hits+a_pr_hits,nb_DIAM);
        NWCountsGlobNorm->Fill(n_pr_hits,nb_NW);

        TString ExitChannel = Form("%dn%dp%da",n_pr_hits,p_pr_hits,a_pr_hits);
        ExitChannel.ReplaceAll("0n","");
        ExitChannel.ReplaceAll("0p","");
        ExitChannel.ReplaceAll("0a","");

        GlobalExitChannelDistrib->GetXaxis()->FindBin(ExitChannel);
        GlobalExitChannelDistrib->Fill(ExitChannel,1);


        TString HistName = Form("ECDist_%dn%dp%da",nb_NW,nb_DIAM_proton,nb_DIAM_alpha);
        TH1 *h = (TH1*)gROOT->FindObject(HistName);
        if(h != 0x0)
        {
            h->GetXaxis()->FindBin(ExitChannel);
            h->Fill(ExitChannel,1);
        }
    }

    //    if(entry>10000)
    //    {
    //        Terminate();
    //        gROOT->ProcessLine(".QQQQQQQQQQQQQQQQQQQQ");
    //    }

    if(  (g_pr_hits >= fN_Gamma[0] && g_pr_hits <= fN_Gamma[1]) &&
         (n_pr_hits >= fN_Neutron[0] && n_pr_hits <= fN_Neutron[1]) &&
         (p_pr_hits >= fN_Proton[0] && p_pr_hits <= fN_Proton[1]) &&
         (a_pr_hits >= fN_Alpha[0] && a_pr_hits <= fN_Alpha[1]) &&
         (nb_AGATA >= fN_Agata[0] && nb_AGATA <= fN_Agata[1]) &&
         (nb_DIAM >= fN_Diamant[0] && nb_DIAM <= fN_Diamant[1]) &&
         (nb_DIAM_proton >= fN_Diamant_proton[0] && nb_DIAM_proton <= fN_Diamant_proton[1]) &&
         (nb_DIAM_alpha >= fN_Diamant_alpha[0] && nb_DIAM_alpha <= fN_Diamant_alpha[1]) &&
         (nb_NW >= fN_NW[0] && nb_NW <= fN_NW[1]))
    {
        enlist->Enter(fChain->GetChainEntryNumber(entry));

        if(PlotEfficiencies == 1)
        {
            AGATACountsPrNorm->Fill(g_pr_hits,nb_AGATA);
            DiamantCountsPrNorm->Fill(p_pr_hits+a_pr_hits,nb_DIAM);
            NWCountsPrNorm->Fill(n_pr_hits,nb_NW);

            AGATACountsEvNorm->Fill(g_pr_hits,nb_AGATA);
            DiamantCountsEvNorm->Fill(p_pr_hits+a_pr_hits,nb_DIAM);
            NWCountsEvNorm->Fill(n_pr_hits,nb_NW);

            AGATACountsGlobNorm->Fill(g_pr_hits,nb_AGATA);
            DiamantCountsGlobNorm->Fill(p_pr_hits+a_pr_hits,nb_DIAM);
            NWCountsGlobNorm->Fill(n_pr_hits,nb_NW);

            TString ExitChannel = Form("%dn%dp%da",n_pr_hits,p_pr_hits,a_pr_hits);
            ExitChannel.ReplaceAll("0n","");
            ExitChannel.ReplaceAll("0p","");
            ExitChannel.ReplaceAll("0a","");

            GlobalExitChannelDistrib->GetXaxis()->FindBin(ExitChannel);
            GlobalExitChannelDistrib->Fill(ExitChannel,1);
        }

        return kTRUE;
    }
    else
    {
        return kFALSE;
    }
}

void AncillaryFilter::SlaveTerminate()
{
    // The SlaveTerminate() function is called after all entries or objects
    // have been processed. When running with PROOF SlaveTerminate() is called
    // on each slave server.

}

void AncillaryFilter::Terminate()
{
    // The Terminate() function is the last function to be called during
    // a query. It always runs on the client, it can be used to present
    // the results graphically or save the results to file.

    cout<<endl;
    cout<<"Entry list completed ! "<<endl;
    cout<<"Number of events in the list : "<<enlist->GetN()<<" ("<<(double)enlist->GetN()/(double)NEntries*100<<" % of input events)"<<endl;

    TFile *FileOut = new TFile(fFileName,"UPDATE");

    enlist->Write();
    FileOut->Close();
    delete FileOut;

    if(PlotEfficiencies != 0)
    {
        for(int i=1 ; i<=50 ; i++)
        {
            TH1 *h = AGATACountsPrNorm->ProjectionX();
            Int_t EntriesAGATA = h->GetBinContent(i);
            delete h;

            h = DiamantCountsPrNorm->ProjectionX();
            Int_t EntriesDiam = h->GetBinContent(i);
            delete h;

            h = NWCountsPrNorm->ProjectionX();
            Int_t EntriesNW = h->GetBinContent(i);
            delete h;

            for(int j=1 ; j<=50 ; j++)
            {
                if(EntriesAGATA!=0) AGATACountsPrNorm->SetBinContent(i,j,(double)AGATACountsPrNorm->GetBinContent(i,j)/((double)EntriesAGATA));
                if(EntriesDiam!=0) DiamantCountsPrNorm->SetBinContent(i,j,(double)DiamantCountsPrNorm->GetBinContent(i,j)/((double)EntriesDiam));
                if(EntriesNW!=0) NWCountsPrNorm->SetBinContent(i,j,(double)NWCountsPrNorm->GetBinContent(i,j)/((double)EntriesNW));
            }
        }

        for(int j=1 ; j<=50 ; j++)
        {
            TH1 *h = AGATACountsEvNorm->ProjectionY();
            Int_t EntriesAGATA = h->GetBinContent(j);
            delete h;

            h = DiamantCountsEvNorm->ProjectionY();
            Int_t EntriesDiam = h->GetBinContent(j);
            delete h;

            h = NWCountsEvNorm->ProjectionY();
            Int_t EntriesNW = h->GetBinContent(j);
            delete h;

            for(int i=1 ; i<=50 ; i++)
            {
                if(EntriesAGATA!=0) AGATACountsEvNorm->SetBinContent(i,j,(double)AGATACountsEvNorm->GetBinContent(i,j)/((double)EntriesAGATA));
                if(EntriesDiam!=0) DiamantCountsEvNorm->SetBinContent(i,j,(double)DiamantCountsEvNorm->GetBinContent(i,j)/((double)EntriesDiam));
                if(EntriesNW!=0) NWCountsEvNorm->SetBinContent(i,j,(double)NWCountsEvNorm->GetBinContent(i,j)/((double)EntriesNW));
            }
        }

        Int_t Entries = AGATACountsGlobNorm->GetEntries();

        for(int i=1 ; i<=50 ; i++)
        {
            for(int j=1 ; j<=50 ; j++)
            {
                AGATACountsGlobNorm->SetBinContent(i,j,(double)AGATACountsGlobNorm->GetBinContent(i,j)/((double)Entries));
                DiamantCountsGlobNorm->SetBinContent(i,j,(double)DiamantCountsGlobNorm->GetBinContent(i,j)/((double)Entries));
                NWCountsGlobNorm->SetBinContent(i,j,(double)NWCountsGlobNorm->GetBinContent(i,j)/((double)Entries));
            }
        }

        TFile *fHists = new TFile(EfficiencyPlotFileName,"recreate");

        TString Name = GlobalExitChannelDistrib->GetName();
        TString Title = GlobalExitChannelDistrib->GetName();
        Int_t NBins = GlobalExitChannelDistrib->FindLastBinAbove(0);

        TH1F *htemp = (TH1F*)GlobalExitChannelDistrib->Clone("htemp");
        delete GlobalExitChannelDistrib;
        GlobalExitChannelDistrib = new TH1F(Name,Title,NBins,0,NBins);
        GlobalExitChannelDistrib->GetXaxis()->SetLabelSize(0.1);
        GlobalExitChannelDistrib->GetYaxis()->SetTitle("Rates");
        GlobalExitChannelDistrib->GetYaxis()->SetTitleOffset(0.8);
        GlobalExitChannelDistrib->GetYaxis()->SetTitleSize(0.05);

        for(int i=0 ; i<NBins ; i++)
        {
            int bin = htemp->GetMaximumBin();
            Double_t Value = (double)htemp->GetBinContent(bin)/((double)Entries);
            TString Label = htemp->GetXaxis()->GetBinLabel(bin);
            htemp->SetBinContent(bin,0);
            GlobalExitChannelDistrib->SetBinContent(i+1,Value);
            GlobalExitChannelDistrib->GetXaxis()->SetBinLabel(i+1,Label);
        }
        delete htemp;

        if(PlotEfficiencies == 2)
        {
            fHists->mkdir("ExitChannelsDist");
            fHists->cd("ExitChannelsDist");

            for(int i=0 ; i<=5 ; i++)
            {
                for(int j=0 ; j<=5 ; j++)
                {
                    for(int k=0 ; k<=1 ; k++)
                    {
                        TH1 *h = ExitChannelDistrib[i][j][k];
                        h->Scale(1./((double)h->GetEntries()));
                        TString Name = h->GetName();
                        TString Title = h->GetName();

                        TH1F *htemp = (TH1F*)h->Clone("htemp");
                        delete h;
                        ExitChannelDistrib[i][j][k] = new TH1F(Name,Title,NBins,0,NBins);
                        h = ExitChannelDistrib[i][j][k];
                        h->GetXaxis()->SetLabelSize(0.1);
                        h->GetYaxis()->SetTitle("Rates");
                        h->GetYaxis()->SetTitleOffset(0.8);
                        h->GetYaxis()->SetTitleSize(0.05);

                        for(int i=0 ; i<NBins ; i++)
                        {
                            int bin = htemp->GetXaxis()->FindBin(GlobalExitChannelDistrib->GetXaxis()->GetBinLabel(i+1));
                            h->SetBinContent(i+1,htemp->GetBinContent(bin));
                            h->GetXaxis()->SetBinLabel(i+1,GlobalExitChannelDistrib->GetXaxis()->GetBinLabel(i+1));
                        }
                        delete htemp;
                        h->Write();
                    }
                }
            }
        }

        GlobalExitChannelDistrib->Write();

        fHists->mkdir("GlobNorm");
        fHists->cd("GlobNorm");
        AGATACountsGlobNorm->Write();
        DiamantCountsGlobNorm->Write();
        NWCountsGlobNorm->Write();

        fHists->mkdir("PrNorm");
        fHists->cd("PrNorm");
        AGATACountsPrNorm->Write();
        DiamantCountsPrNorm->Write();
        NWCountsPrNorm->Write();

        fHists->mkdir("EvNorm");
        fHists->cd("EvNorm");
        AGATACountsEvNorm->Write();
        DiamantCountsEvNorm->Write();
        NWCountsEvNorm->Write();

        fHists->Close();
    }
}

void AncillaryFilter::ReadConfFile()
{
    cout<<endl;
    cout<<"******************************"<<endl;
    cout<<"* Reading Configuration File *"<<endl;
    cout<<"******************************"<<endl;
    cout<<endl;

    std::ifstream FileIn;
    FileIn.open("AncillaryFilter.conf", std::ifstream::in);

    fN_Gamma[0]   = 0;
    fN_Neutron[0] = 0;
    fN_Proton[0]  = 0;
    fN_Alpha[0]   = 0;

    fN_Agata[0]   = 0;
    fN_Diamant[0] = 0;
    fN_Diamant_proton[0] = 0;
    fN_Diamant_alpha[0] = 0;
    fN_NW[0]      = 0;

    fN_Gamma[1]   = kMaxPr_fHits;
    fN_Neutron[1] = kMaxPr_fHits;
    fN_Proton[1]  = kMaxPr_fHits;
    fN_Alpha[1]   = kMaxPr_fHits;

    fN_Agata[1]   = kMaxEv_fHits;
    fN_Diamant[1] = kMaxEv_fHits;
    fN_Diamant_proton[1] = kMaxEv_fHits;
    fN_Diamant_alpha[1] = kMaxEv_fHits;
    fN_NW[1]      = kMaxEv_fHits;

    PlotEfficiencies = 0;
    EfficiencyPlotFileName = "";

    fFileName     = Form("%s_%sEntryList.root",fChain->GetName(),GetOption());

    if(!FileIn)
    {
        cout<<"No Configuration file ==> Exit"<<endl;
        cout<<"An example of configuration file is available in the Conf folder"<<endl;
        gROOT->ProcessLine(".qqqqqqqqqqqqqqqqq");
    }

    string line;
    TString Buffer;

    Int_t *CurrentParameter;

    while(FileIn)
    {
        getline(FileIn,line);
        Buffer = line;

        if(Buffer.BeginsWith("#")) continue;
        if(Buffer.BeginsWith("N_Gamma")) CurrentParameter = fN_Gamma;
        else if(Buffer.BeginsWith("N_Neutron")) CurrentParameter = fN_Neutron;
        else if(Buffer.BeginsWith("N_Proton")) CurrentParameter = fN_Proton;
        else if(Buffer.BeginsWith("N_Alpha")) CurrentParameter = fN_Alpha;
        else if(Buffer.BeginsWith("Hits_Agata")) CurrentParameter = fN_Agata;
        else if(Buffer.BeginsWith("Hits_Diamant")) CurrentParameter = fN_Diamant;
        else if(Buffer.BeginsWith("Hits_Proton")) CurrentParameter = fN_Diamant_proton;
        else if(Buffer.BeginsWith("Hits_Alpha")) CurrentParameter = fN_Diamant_alpha;
        else if(Buffer.BeginsWith("Hits_NW")) CurrentParameter = fN_NW;
        else if(Buffer.BeginsWith("FilterName")) CurrentParameter = 0x0;
        else if(Buffer.BeginsWith("EffPlot")) CurrentParameter = 0x0;
        else continue;

        Buffer.ReplaceAll("\t"," ");
        TObjArray *loa=Buffer.Tokenize(" ");
        Int_t NValues = loa->GetEntries();

        if(CurrentParameter == 0x0)
        {
            if(Buffer.BeginsWith("FilterName"))
            {
                if(NValues != 2)
                {
                    cout<<"Wrong Filter name definition, default one is taken"<<endl;
                }
                else
                {
                    fFileName = (TString)fChain->GetName() + "_" + (TString)loa->Last()->GetName() + "_EntryList.root";
                }
                cout<<"EntryList will be stored in "<<fFileName<<endl;
                continue;
            }
            else if(Buffer.BeginsWith("EffPlot"))
            {
                if(NValues != 2)
                {
                    cout<<"Wrong Efficiencies plot definition, no plot will be stored"<<endl;
                }
                else
                {
                    PlotEfficiencies = ((TString)loa->Last()->GetName()).Atoi();
                    if(PlotEfficiencies<0) PlotEfficiencies = 0;
                    if(PlotEfficiencies>2) PlotEfficiencies = 2;
                }
                if(PlotEfficiencies == 0)
                {
                    cout<<"Efficiencies will not be plotted"<<endl;
                    continue;
                }
                if(PlotEfficiencies == 1)
                {
                    EfficiencyPlotFileName = fFileName.Copy().ReplaceAll("EntryList","Efficiencies");
                    cout<<"Efficiencies will be plotted for filtered events in "<<EfficiencyPlotFileName<<endl;
                }
                else
                {
                    EfficiencyPlotFileName = (TString)fChain->GetName() + "_Global_Efficiencies.root";
                    cout<<"Global Efficiencies will be plotted in "<<EfficiencyPlotFileName<<endl;
                }

                continue;
            }

        }
        else
        {
            if(NValues == 2)
            {
                CurrentParameter[0] = ((TString)loa->Last()->GetName()).Atoi();
                CurrentParameter[1] = CurrentParameter[0];
            }
            else if(NValues == 3)
            {
                CurrentParameter[0] = ((TString)loa->At(1)->GetName()).Atoi();
                if((TString)loa->Last()->GetName() != "Max") CurrentParameter[1] = ((TString)loa->Last()->GetName()).Atoi();
            }
            else
            {
                cout<<setw(13)<<loa->First()->GetName()<<" : Wrong number of entries in the configuration file"<<endl;
            }
            cout<<left<<setw(13)<<loa->First()->GetName()<<" restricted to [ "<<setw(4)<<CurrentParameter[0]<<" ; "<<setw(4)<<CurrentParameter[1]<<"]"<<endl;

        }

        delete loa;
    }
    FileIn.close();
}