#if !defined(__CINT__) || defined(__MAKECINT__)


// ROOT includes //
#include <TString.h>
#include <TH1F.h>
#include <TH2F.h>

// standard libraries includes //
#include <stdlib.h>
#include <iostream>
//#include <fstream>
#include <iomanip>
using namespace std;

// PAMELA includes //
#include <PamLevel2.h>

#endif

//======================
// global variables
//======================
Int_t tot=0;
Int_t sel[]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};


//======================
// HISTOGRAMS declaration
//======================
TH1F* hp0gen;
TH1F* hp0trk;
TH1F* hpmtrk;



//===============================================================================
//
//
//
//
//===============================================================================
bool Select( PamLevel2* event ){
//    return false;
    TrkParams::SetVerboseMode();

    //-------------------------------------------
    // check if the event contains selection info
    //-------------------------------------------
    if( event->GetTrkLevel2()==0 ){
	cout << "Missing TrkLevel2 object "<<endl;
	return false;
    }
    if( event->GetCaloLevel2()==0 ){
	cout << "Missing CaloLevel2 object "<<endl;
	return false;
    }
    if( event->GetRunInfo()==0 ){
	cout << "Missing RunInfo object "<<endl;
	return false;
    }
    if( event->GetGPamela()==0 ){
	cout << "Missing GPamela object "<<endl;
	return false;
    }

    tot++;


    // -------------------------------------------------
    // no hits on CAT and CARD
    // -------------------------------------------------
    if( event->GetGPamela()->Nthcat > 0 ||  event->GetGPamela()->Nthcard > 0) return false;
    sel[0]++;

    // --------------------------------------------------
    // select simulated events where the primary particle
    // traverse S2 and S3
    // --------------------------------------------------
    bool  S2OK = false; // the primary particle hits S3
    bool  S3OK = false; // the primary particle hits S3

    for(int ih=0; ih<event->GetGPamela()->Nthtof; ih++){ //loop over tof hits

	// S2
	if( 
	    (Int_t) event->GetGPamela()->Ipartof[ih] == event->GetGPamela()->Ipa &&                                //primary type
	    ( (Int_t) event->GetGPamela()->Ipltof[ih] == 4 || (Int_t) event->GetGPamela()->Ipltof[ih] == 3) &&     //S21 or S22
	    fabs( (event->GetGPamela()->P0 - event->GetGPamela()->P0tof[ih]) / event->GetGPamela()->P0 ) < 0.05 && //energy consistent with the primary (1%)
	    true){
	    S2OK = true;
	}
	// S3
	if( 
	    (Int_t) event->GetGPamela()->Ipartof[ih] == event->GetGPamela()->Ipa &&                                //primary type
	    ( (Int_t) event->GetGPamela()->Ipltof[ih] == 5 || (Int_t) event->GetGPamela()->Ipltof[ih] == 6) &&     //S31 or S32
	    fabs( (event->GetGPamela()->P0 - event->GetGPamela()->P0tof[ih]) / event->GetGPamela()->P0 ) < 0.05 && //energy consistent with the primary (1%)
	    true){
	    S3OK = true;
	}
    }
    if(!S2OK || !S3OK)return false;
    sel[1]++;
    

    return true;
    

}
//===============================================================
// Create histograms
//
//
//
//
//
//===============================================================
void CreateHistos( PamLevel2* event , TFile* outf ){
    cout << "Creating histos..."<<endl;
    gROOT->cd();//create histos in memory
    TString hid;
    TString title;

    //===========================================
    // rigidity binning
    Int_t   nbins = 10;
    Float_t rmin  = 1.;
    Float_t rmax  = 100.;
    Float_t dr = (log10(rmax)-log10(rmin))/nbins;
    TArrayD* arr = new TArrayD(nbins+1);
    arr->AddAt((Double_t)rmin,0);
    for(Int_t i=1; i<=nbins; i++){
        Double_t previous = arr->At(i-1);
        Double_t value=0;
        value = (Double_t) pow(10,log10((Float_t)previous)+dr);
        arr->AddAt(value,i);
    };
    Double_t *rbinpt = arr->GetArray();
    //===========================================

    //--------------------------------
    // histos
    //--------------------------------
    hp0gen = new TH1F("hp0gen","Generated momentum (all events)",nbins,rbinpt);
    hp0trk = new TH1F("hp0trk","Generated momentum ",nbins,rbinpt);
    hpmtrk = new TH1F("hpmtrk","Measured  momentum ",nbins,rbinpt);



    //--------------------------------
    // initialization
    //--------------------------------
    cout << "Opening DB connection..."<<endl;
    event->SetDBConnection();

}
//===============================================================
// Retrieve histograms
//
//
//
//
//
//===============================================================
void RetrieveHistos(){

    gROOT->cd();
    TString hid;
    hp0gen=dynamic_cast<TH1F*>(gDirectory->FindObject("hp0gen"));
    hp0trk=dynamic_cast<TH1F*>(gDirectory->FindObject("hp0trk"));
    hpmtrk=dynamic_cast<TH1F*>(gDirectory->FindObject("hpmtrk"));

}
//===============================================================
// Get histograms (... non ho capito la differenza con retrieve)
//
//
//
//
//
//===============================================================
void GetHistos(){

    TString hid;
}
//===============================================================
// Fill histograms (called event-by-event)
//
//
//
//
//
//
//=============================================================== 
void FillHistos( PamLevel2* event ){
//  ----------------------------------
//  retrieve histos
//  ----------------------------------
    RetrieveHistos();




    cout << "================================================="<<endl;
    cout << "part        "<<event->GetGPamela()->Ipa <<endl;
    cout << "P0 (GeV/c)  "<<event->GetGPamela()->P0 <<endl;
    cout << "X0 (cm)     "<<event->GetGPamela()->X0 << endl;
    cout << "Y0 (cm)     "<<event->GetGPamela()->Y0    << endl; 
    cout << "Z0 (cm)     "<<event->GetGPamela()->Z0    << endl; 
    cout << "Theta (rad) "<<event->GetGPamela()->Theta  << endl;
    cout << "Phi (rad)   "<<event->GetGPamela()->Phi    << endl;
    cout <<"Spectrometer hits: "<< endl;
    cout << "hit --";
    cout << "part  ";
    cout << "plane ";
    cout << "X           ";
    cout << "Y           ";
    cout << "Z           ";
    cout << "En.release  ";
    cout << "P           ";
    cout << endl;
    for(int ih=0; ih<event->GetGPamela()->Nthspe; ih++){
	cout << setw(4) << ih << "--";
	cout << setw(6) <<   (Int_t) event->GetGPamela()->Iparspe[ih];//part
	cout << setw(6) <<   (Int_t) event->GetGPamela()->Itrpb[ih];//0-5 o 1-6 piano
	cout << setw(12)<< event->GetGPamela()->Xavspe[ih];
	cout << setw(12)<< event->GetGPamela()->Yavspe[ih];
	cout << setw(12)<< event->GetGPamela()->Zavspe[ih];
	cout << setw(12)<< event->GetGPamela()->Erelspe[ih];//en.release
	cout << setw(12)<< event->GetGPamela()->P0spe[ih];//momentum
	cout << endl;
    }
    




    hp0gen->Fill(event->GetGPamela()->P0);
    
    // one fitted track
    if( event->GetTrkLevel2()->GetNTracks()>0 ){
	sel[2]++;
	hp0trk->Fill(event->GetGPamela()->P0);

//  ---------------------------------
//  Retrieve the track
//  ---------------------------------
	event->SetSortingMethod("+CAL");        // use only calorimeter for image rejection
	PamTrack  *track = event->GetTrack(0);  // get the whole track 
	TrkTrack  *trk   = track->GetTrkTrack();// get tracker track variables
	
	hpmtrk->Fill(trk->GetRigidity());
	
    }
        
}
//===============================================================
// Save histograms
//
//
//
//
//
//===============================================================
void SaveHistos( TFile * outf ){

TString hid; 

gROOT->cd();
//  ----------------------------------
//  retrieve histos
//  ----------------------------------
    RetrieveHistos();

//  ----------------------------------
//  save on file
//  ----------------------------------
    outf->cd();

    hp0gen->Write();
    hp0trk->Write();
    hpmtrk->Write();

    cout << "####------------------------------------------------------------------" <<endl;
    cout << "#### total :"<<tot<<endl;
    if(!tot)return;
    cout << "#### no GP hits on CAT/CARD   :"<<sel[0]<< "  (" << ((float)sel[0]/(float)tot) <<")"<<endl;
    cout << "#### GP hits on S2/S3         :"<<sel[1]<< "  (" << ((float)sel[1]/(float)tot) <<")"<<endl;
    cout << "#### 1 fitted track           :"<<sel[2]<< "  (" << ((float)sel[2]/(float)tot) <<")"<<endl;
    cout << "####------------------------------------------------------------------" <<endl;

}