source: JEM-EUSO/esaf_lal/branches/camille/camille_work/KIT_phase1_study/macros/Analysis.C @ 168

//
//  Analysis.C
//
//
//  Created by moretto on 29/05/13.
//
//

Int_t nb_entries = 0;

const char* inputfile="/Users/moretto/Documents/thesis/ESAF/Code/esaf_lal/esaf/camille_work/KIT_phase1_study/30_ns_TimeResolution__100_event/simu.root";
const char* inputfile2="/Users/moretto/Documents/thesis/ESAF/Code/esaf_lal/esaf/camille_work/KIT_phase1_study/10_ns_TimeResolution__100_event/simu.root";
const char* inputfile3="/Users/moretto/Documents/thesis/ESAF/Code/esaf_lal/esaf/camille_work/KIT_phase1_study/5_ns_TimeResolution__100_event/simu.root";

EEvent*          ev = 0;
EDetector*       gDetector = NULL;
EDetPhoton*      gDetPhoton = NULL;
EBunchPhotons*   gBunchPhotons = NULL;

void gettree();
void init();
void process_profile();
Bool_t loadevent(Int_t i=0);

//______________________________________________________________________________
void Analysis()
{
    TProfile *hprof_multi[2];
    
    //Initialisation & process
    init(inputfile);
    process_profile();
    hprof_multi[1] = hprof;
    
    init(inputfile3);
    process_profile();
    hprof_multi[2] = hprof;
    
    TCanvas *c4 = new TCanvas("c4", "Efficiency vs GTU");
    TProfile *hprof_efficiency = new TProfile("hprof_efficiency", "Efficiency vs GTU", 128, 0, 128);
    
    for (Int_t i=1; i<128; i++)
    {
        Float_t bin1 = hprof_multi[1]->GetBinContent(i);
        Float_t bin2 = hprof_multi[2]->GetBinContent(i);
        
        if (bin2 == 0)
        {
            hprof_efficiency->Fill(i, 0); //Il faut entrer des entiers dedans sinon ça ne marchera pas!!!!
        }
        
        else
        {
            hprof_efficiency->Fill(i, bin1/bin2); //Il faut entrer des entiers dedans sinon ça ne marchera pas!!!!
        }
    }
    
    hprof_efficiency->Draw();
}

//______________________________________________________________________________
void process_profile ()
{
    //Déclaration des histogrammes et autres
    TH1F *histo_total = new TH1F("histo_total", "Histogram of photons read by the electronics for 100 showers", 128, 0, 128);
    TH1F *histo_multi[100];
    TProfile *hprof = new TProfile("hprof", "Average number of photons per GTU", 128, 0, 128);
    
    //Boucle sur l'ensemble des événements pour créer une distribution totale
    for (Int_t  i=0; i<nb_entries ; i++)
    {
        loadevent(i);
        
        Int_t nb_photon = gDetector->GetNumPhotons();
        
        for (Int_t j=0; j<nb_photon; j++)
        {
            Int_t macrocell_status = gDetector->GetPhoton(j)->GetMacroCell();
            Int_t photon_type = gDetector->GetPhoton(j)->GetType();
            
            if (macrocell_status == 1 & photon_type != 2) //on ne sélectionne que les photons étant détectés par l'électronique & autre que les Cerenkov
            {
                Float_t gtu = gDetector->GetPhoton(j)->GetGtu();
                histo_total->Fill(gtu);
            }
            
        }
        
    }

    TCanvas *c1 = new TCanvas("c1", "Distribution totale");
    histo_total->GetXaxis()->SetTitle("GTU");
    histo_total->Draw();
    
    //Boucle sur l'ensemble des événements pour tracer séparément les histogrammes
    TCanvas *c2 = new TCanvas("c2", "Collection of GTU histograms for 100 events");
    c2->Divide(10,10);
    
    for (Int_t i=0; i<nb_entries; i++)
        histo_multi[i] = new TH1F("histo_multi", "Single shower histograms", 128, 0, 128);
    
    for (Int_t  i=0; i<100 ; i++)
    {
        loadevent(i);
        
        nb_photon = gDetector->GetNumPhotons();
        
        for (Int_t j=0; j<nb_photon; j++)
        {
            macrocell_status = gDetector->GetPhoton(j)->GetMacroCell();
            photon_type = gDetector->GetPhoton(j)->GetType();
            
            if (macrocell_status == 1 & photon_type != 2) //on ne sélectionne que les photons étant détectés par l'électronique & autre que les Cerenkov
            {
                gtu = gDetector->GetPhoton(j)->GetGtu();
                histo_multi[i]->Fill(gtu);
            }
        }
        
        c2->cd(i+1);
        histo_multi[i]->Draw();
        
    }
    
    //Boucle sur le tableau d'histogrammes précédent afin de réaliser un TProfile
    for (Int_t j=0; j<nb_entries; j++)
    {
        for (Int_t i=1; i<=128; i++)
        {
            hprof->Fill(i, histo_multi[j]->GetBinContent(i));
        }
        
    }
    
    TCanvas *c3 = new TCanvas("c3", "TProfile of the GTU distribution for 100 events");
    hprof->Draw();

}

//______________________________________________________________________________
void init( const char* fname)
{
    TFile *f = new TFile(fname);
    gettree();
}

//______________________________________________________________________________
void gettree()
{
    TTree *etree = (TTree*)gDirectory->Get("etree");
    if ( etree )
    {
        nb_entries = etree->GetEntries();
        Printf("etree found with %d events stored.", nb_entries);
    }
}

//______________________________________________________________________________
Bool_t loadevent( Int_t i)
{
    TTree *etree = (TTree*)gDirectory->Get("etree");
    
    if ( !etree )
    {
        Printf("init(): unable to retrieve etree from TFile. Exiting");
        return kFALSE;
    }
    
    // associate the tree and the event
    if ( !ev )
    {
        ev = new EEvent();
        ev->SetBranches(etree);
    }
    
    etree->GetEntry(i);
    
    EDetector *detector = ev->GetDetector();
    
    gDetector = ev->GetDetector();
    
    return kTRUE;
}