Histo_Perfil_Area.C

// root script file
/* Programa para extraer los datos y realizar histogramas
    Elser Lopez, 23/05/2020
    Realiza el histograma del Area del pulso y energia depositada
    Asi como sus perfiles y el fit
*/

#include "TFile.h"
#include "TTree.h"
#include "TGraph.h"
#include "TCanvas.h"

#include <cstdlib>
#include <iostream>

void Histo_Perfil_Area(const char* FileName)
{
    TFile *RFile = new TFile(FileName);
    if( !RFile )
    {
        std::cout<<"Can't open " << FileName << std::endl;
        return;
    }

    TTree *Config;
    TTree *Data;

    Float_t ADC_Bits;
    Float_t ADC_Vref;
    Float_t ADC_Sample_Rate;
    Float_t ADC_VSignal_Offset;
    Float_t ADC_Samples_per_Pulse;
    Float_t ADC_Trigger_Voltage;
    Float_t ADC_Pre_Trigger_Samples;

    Config = (TTree*)RFile->Get("ElecSim_info");
    Config->SetBranchAddress("ADC_Bits",            &ADC_Bits);
    Config->SetBranchAddress("ADC_Vref",            &ADC_Vref);
    Config->SetBranchAddress("ADC_Vin_Offset",      &ADC_VSignal_Offset);
    Config->SetBranchAddress("Sample_Rate",         &ADC_Sample_Rate);
    Config->SetBranchAddress("Trigger_Voltage",     &ADC_Trigger_Voltage);
    Config->SetBranchAddress("Samples_per_Pulse",   &ADC_Samples_per_Pulse);
    Config->SetBranchAddress("Pre_Trigger_Samples", &ADC_Pre_Trigger_Samples);

    Config->GetEntry(0);
	
	Int_t 		PDG_Code;
    Double_t 	Energy;
	Double_t 	Zenith_angle;
	Int_t 		Direction;
	Double_t 	Deposited_Energy;
	Double_t 	Track_Length;
	Int_t	 	Cherenkov_Photon_Count;
	Int_t	 	PMT_Photon_Count;

    std::vector< Double_t > *Event_Data = 0;

    Data = (TTree*)RFile->Get("ElecSim_Output");
    Data->SetBranchAddress("PDG_Code",                  &PDG_Code);
    Data->SetBranchAddress("Energy",                    &Energy);
	Data->SetBranchAddress("Zenith_angle",              &Zenith_angle);
	Data->SetBranchAddress("Direction",                 &Direction);
	Data->SetBranchAddress("Deposited_Energy",          &Deposited_Energy);
	Data->SetBranchAddress("Track_Length",              &Track_Length);
	Data->SetBranchAddress("Cherenkov_Photon_Count",    &Cherenkov_Photon_Count);
	Data->SetBranchAddress("PMT_Photon_Count",          &PMT_Photon_Count);
    Data->SetBranchAddress("Digitalized_Data",          &Event_Data);

    Long_t NEntries = Data->GetEntries();

    //if(EventNumber >= NEntries)
    //{
      //  std::cout<< "Invalid event number" << std::endl;
      //  return;
    //}



Double_t Area_por_pulso[(Int_t)NEntries]; // Vector para guardar todas la areas de los pulsos
//Double_t Minimo_por_pulso[(Int_t)NEntries]; //Vector para guardos todos los minimos de los pulso
// Con esto saca los datos de todos los pulso
Int_t EventNumber;
for(Int_t i=0; i< NEntries; i++)
 {
    EventNumber=i;
    Data->GetEntry(EventNumber);

    Double_t t_array[(Int_t)ADC_Samples_per_Pulse];
    Double_t *d_array = Event_Data->data();

    //std::cout<< "yes" << *d_array << std::endl;

    //Genera los puntos en x
    for(Int_t i=0; i<ADC_Samples_per_Pulse; i++ )
    {
        t_array[i] = i;
       // std::cout<< "i=" << i << 'y' << *d_array++ << std::endl;
    }




/* Calculando el area del pulso */

    /* calculo la media para asi obtener el Base-Line*/
    Double_t Dmean[(Int_t)ADC_Pre_Trigger_Samples];
    Double_t Mean=0; //Media de los primeros puntos antes del cambio puntos analizados del eje y, para cada pulso
    for(Int_t i=0; i<(Int_t)ADC_Pre_Trigger_Samples; i++)
    {
        Dmean[i]=d_array[i];
    }

    //Calculo la media
    for(Int_t i=0; i<(Int_t)ADC_Pre_Trigger_Samples; i++)
    {
        Mean=Mean+Dmean[i]/(Int_t)ADC_Pre_Trigger_Samples;
    }



    //Calculo el area del pulso
    Double_t h; //Altura en la formula del trapecio
    Double_t Base; 
    Double_t AreaB=0; //Area bajo la curva

    Double_t AreaP=0; //Area del pulso
    for(Int_t i=0; i<ADC_Samples_per_Pulse-20; i++ )
        {//Con esto calculo el area bajo la curva
            h=(t_array[i+1]-t_array[i]);
            Base=d_array[i]+d_array[i+1];
            AreaB=AreaB+h*(Base)/2.0;
            //std::cout << "Area bajo la curva=" << AreaB << " valores en y=" << d_array[i] << std::endl;
        }
   // AreaS=Mean*( (Int_t)ADC_Samples_per_Pulse-20 );

    AreaP=Mean*((Int_t)ADC_Samples_per_Pulse-20 )-AreaB;
    

    Area_por_pulso[EventNumber]=AreaP; //Agrego al vector 

}

    //Creo un histograma
    TH2D *EDepvsArea= new TH2D("AreavsEd", " ", 100,0,0,100,0,0);
    //Creo un perfil
    auto *EDepvsAreaP= new TProfile("AreavsEd", " ", 100,0,1000,0,350000);
    //TH2D *EDepvsMax= new TH2D("EdvsMax", " ", 100,0,0,100,0,0);
    EDepvsArea->GetYaxis()->SetTitle("#acute{A}rea del pulso");
    EDepvsArea->GetXaxis()->SetTitle("Energ#acute{i}a Depositada [MeV]");

    EDepvsAreaP->GetYaxis()->SetTitle("#acute{A}rea del pulso");
    EDepvsAreaP->GetXaxis()->SetTitle("Energ#acute{i}a Depositada [MeV]");




    //Obtengo los valores de la energia depositada y el area de los puntos
    for(Long_t i=0; i<NEntries; i++)
    {
        Data->GetEntry(i);
        EDepvsArea->Fill(Deposited_Energy,Area_por_pulso[i]);
        EDepvsAreaP->Fill(Deposited_Energy,Area_por_pulso[i]);
    }


    //Creo la grafica 

    TF1 * f1 = new TF1("f1","pol1");
    std::cout<< "Parametros energia Depositada vs Area del Pulso" << std::endl;
    EDepvsAreaP->Fit(f1,"S");
    TFitResultPtr r = EDepvsAreaP->Fit(f1,"S");
    //TMatrixDSym C = r->GetCorrelationMatrix();
    r->Print("V");

    TCanvas *a0Canvas = new TCanvas();
    EDepvsArea->Draw();
    a0Canvas->SaveAs("His_EDepvsArea.pdf");
    EDepvsAreaP->Draw();
    a0Canvas->SaveAs("His_EDepvsArea_P.pdf");

}