////////////////////////////////////////////////////////////////
//
//  Histogram set for tracker-efficiency study.
//
//  It contains the functions:
//
//  - CreateHistos()
//  - RetrieveHistos()
//  - FillHistos(PamLevel2* event)
//  - Save Histos()
//
//  Usage (example):
//  [] .L My-histo-set-macro.C+
//  [] .L My-selection.C+
//  [] .L Loop.C+
//  [] Loop("data-dir/","file-list.txt",100000,"+ALL fillTree fillHistos","output-file.root")
//  
////////////////////////////////////////////////////////////////

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

#include <TString.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TMath.h>
#include <TLine.h>
#include <TPolyMarker.h>
#include <TSelector.h>
#include <TFile.h>

#include <stdlib.h>
#include <iostream>
using namespace std;

#include <PamLevel2.h>
#include <CaloAxis.h>

#endif
//======================
// HITOGRAMS declaration
//======================

TH2F *hchi2def;
TH2F *hchi2rig;
TH2F *hdedxtofvsdedxtrk;
TH2F *hdedxtrkvsrig ;
TH2F *hdedxtofvsrig ;
TH2F *hbetavsrig ;
TH2F *hdedxtrkvsdef ;
TH2F *hbetavsdef ;
TH2F *hbetavsdef2 ;

TH1F *herrdef;
TH1F *herrx[6];
TH1F *herry[6];

TH1F *hdef ;
TH1F *hrig1 ;
TH1F *hrig2 ;
TH1F *hen1 ;
TH1F *hen2 ;


Int_t tot=0;
Int_t sel[]={0,0,0,0,0,0,0,0};

//===============================================================
// Get ToF dE/dx
//===============================================================
// Float_t GetdEdx( PamTrack *trk ){
//     Float_t dedx = 0.;
//     Int_t nval=0;
//     Bool_t ISNULL = false;
//     for (Int_t i=0; i<trk->GetToFTrack()->npmtadc; i++){
// 	if(trk->GetToFTrack()->dedx[i] == 0)ISNULL=true;	
// 	dedx += (trk->GetToFTrack()->dedx).At(i)*(Int_t)(!ISNULL); 
// 	nval += (Int_t)(!ISNULL);
//     };
//     //
//     if(nval)dedx=dedx/nval;
// //    cout <<" -- dedx  "<<dedx<<endl;
//     return(dedx);    
// }

// float GetZ(float beta, float dedx){

//     int nspec=2;
//     float p0[]={0.904,3.373};
//     float p1[]={0.295,2.756};
//     float zz[]={1.,4.};

//     if(beta>1)beta=1;
//     float f0 = p0[0]/beta+p1[0];
//     float f1 = p0[1]/beta+p1[1];
//     float z0 = zz[0];
//     float z1 = zz[1];
//     float aa = (z0-z1)/(f0-f1);
//     float bb = -aa*f0+z0;
//     if( (bb+aa*dedx)==0 )return 0;
//     return sqrt(bb+aa*dedx);

// }


//===============================================================
// Create histograms
//===============================================================
void CreateHistos( PamLevel2* event, TFile* outf ){

    gROOT->cd();//create histos in memory

    TString hid;
    TString title;
    //===========================================
    // rigidity binning
    Int_t   nbins = 200;
    Float_t rmin  = .1;
    Float_t rmax  = 1000.;
    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();
    //===========================================

    hchi2def = new TH2F("hchi2def","(chi**2)**0.25 vs deflection",100,0.,2.5,200,0.,10.);
    hchi2rig = new TH2F("hchi2rig","#chi^{2} vs R",nbins,rbinpt,1000,0.,500.);
    hdedxtofvsdedxtrk  = new TH2F("hdedxtofvsdedxtrk","dE/dx (ToF-average) vs dE/dx (trk-average)",600,0.,40.,600,0.,40.);
    hdedxtrkvsrig      = new TH2F("hdedxtrkvsrig","dE/dx (Trk-average) vs R",nbins,rbinpt,600,0.,40.);
    hdedxtofvsrig      = new TH2F("hdedxtofvsrig","dE/dx (ToF-average) vs R",nbins,rbinpt,600,0.,40.);
    hbetavsrig         = new TH2F("hbetavsrig","beta vs R",nbins,rbinpt,200,0.4,1.5);
    hrig1              = new TH1F("hrig1","Rigidity distribution Z=1",nbins,rbinpt);
    hrig2              = new TH1F("hrig2","Rigidity distribution Z=2",nbins,rbinpt);
    hen1               = new TH1F("hen1","Kin.Energy per nucleon Z=1",nbins,rbinpt);
    hen2               = new TH1F("hen2","Kin.Energy per nucleon Z=2",nbins,rbinpt);

    hdedxtrkvsdef      = new TH2F("hdedxtrkvsdef","dE/dx (Trk-average) vs deflection",500,-10.,10.,600,0.,40.);
    hbetavsdef         = new TH2F("hbetavsdef","beta vs deflection",500,-10.,10.,200,0.4,1.5);
    hbetavsdef2        = new TH2F("hbetavsdef2","beta vs deflection (Z=1)",500,-10.,10.,200,0.4,1.5);
    hdef               = new TH1F("hdef","Deflection distribution",5000,-2.5,2.5);

    herrdef            = new TH1F("herrdef","Deflection error",1000,0.,0.01);
    
    for(Int_t ip=0; ip<6; ip++){
	hid.Form("herrx%i",ip);
	title = ""; title.Form("Spatial resolution - plane %i (x)",ip);
	herrx[ip] = new TH1F(hid.Data(),title.Data(),500,0.,0.005);
	hid.Form("herry%i",ip);
	title = ""; title.Form("Spatial resolution - plane %i (y)",ip);
	herry[ip] = new TH1F(hid.Data(),title.Data(),500,0.,0.005);
    }

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

    TString hid;

    hchi2def           = dynamic_cast<TH2F*>(gDirectory->FindObject("hchi2def"));
    hchi2rig           = dynamic_cast<TH2F*>(gDirectory->FindObject("hchi2rig"));
    hdedxtofvsdedxtrk  = dynamic_cast<TH2F*>(gDirectory->FindObject("hdedxtofvsdedxtrk"));
    hdedxtrkvsrig      = dynamic_cast<TH2F*>(gDirectory->FindObject("hdedxtrkvsrig"));
    hdedxtofvsrig      = dynamic_cast<TH2F*>(gDirectory->FindObject("hdedxtofvsrig"));
    hbetavsrig         = dynamic_cast<TH2F*>(gDirectory->FindObject("hbetavsrig"));
    hdedxtrkvsdef      = dynamic_cast<TH2F*>(gDirectory->FindObject("hdedxtrkvsdef"));
    hbetavsdef         = dynamic_cast<TH2F*>(gDirectory->FindObject("hbetavsdef"));
    hbetavsdef2        = dynamic_cast<TH2F*>(gDirectory->FindObject("hbetavsdef2"));
    hdef               = dynamic_cast<TH1F*>(gDirectory->FindObject("hdef"));
    hrig1              = dynamic_cast<TH1F*>(gDirectory->FindObject("hrig1"));
    hrig2              = dynamic_cast<TH1F*>(gDirectory->FindObject("hrig2"));
    hen1               = dynamic_cast<TH1F*>(gDirectory->FindObject("hen1"));
    hen2               = dynamic_cast<TH1F*>(gDirectory->FindObject("hen2"));
    herrdef            = dynamic_cast<TH1F*>(gDirectory->FindObject("herrdef"));
    for(Int_t ip=0; ip<6; ip++){
	hid.Form("herrx%i",ip);
	herrx[ip] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));
	hid.Form("herry%i",ip);
	herry[ip] = dynamic_cast<TH1F*>(gDirectory->FindObject(hid.Data()));
    }
}
//===============================================================
// Fill histograms (called event-by-event)
//===============================================================
 
void FillHistos( PamLevel2* event ){
//  ----------------------------------
//  retrieve histos
//  ----------------------------------
    RetrieveHistos();

//  ----------------------------------
//  single-track events
//  ----------------------------------
    float beta = 0;
    float dedxtof = 0;
    float dedxtrk = 0;
    float rigidity = 0;
    float deflection = 0;
    float chi2 = 0;
    float errdef = 0.;
    float errx[] = {10.,10.,10.,10.,10.,10.};
    float erry[] = {10.,10.,10.,10.,10.,10.};
    
    tot++;
	
    if( event->GetTrkLevel2()->GetNTracks()==0 ) return; 
    PamTrack *track = event->GetTrack(0);
    //------------------------------
    // track selection
    //------------------------------
    if( 
	true
	){
	
	sel[0]++;
	//-------------------
	// Calo variables
	//-------------------
	
	
	//-------------------
	// ToF variables
	//-------------------
	beta    = track->GetToFTrack()->beta[12];
//	dedxtof = GetdEdx(track);
	
	//-------------------
	// tracker variables
	//-------------------
	dedxtrk    = track->GetTrkTrack()->GetDEDX()/1.3; //<<<ATTENZIONE TEMPORANEO!!!!
	rigidity   = track->GetTrkTrack()->GetRigidity();
	deflection = track->GetTrkTrack()->GetDeflection();
	chi2       = track->GetTrkTrack()->chi2;
	errdef     = sqrt(track->GetTrkTrack()->coval[4][4]); 
	for(Int_t ip=0; ip<6; ip++){
	    if(track->GetTrkTrack()->XGood(ip))errx[ip]=track->GetTrkTrack()->resx[ip];
	    if(track->GetTrkTrack()->YGood(ip))erry[ip]=track->GetTrkTrack()->resy[ip];
	}

	//-------------------
	// chi2 selection
	//-------------------
	hchi2def->Fill( fabs(deflection), pow((double)chi2,0.25) );
	hchi2rig->Fill( rigidity, chi2 );
	if(pow((double)chi2,0.25) > (2.5+1.85*fabs(deflection)) )return;
	sel[1]++;
	
//	hdedxtofvsdedxtrk->Fill(dedxtrk,dedxtof);
	hdedxtrkvsrig->Fill(rigidity,dedxtrk);
//	hdedxtofvsrig->Fill(rigidity,dedxtof);
	hbetavsrig->Fill(rigidity,beta);

	hdedxtrkvsdef->Fill(deflection,dedxtrk);
	hbetavsdef->Fill(deflection,beta);

	//------------------------
	// rough charge selection
	//------------------------
	if(dedxtrk > (6.+16./(rigidity*rigidity)))return;
	//--------
	// helium
	//--------
	if( dedxtrk > (2.8+3.3/(rigidity*rigidity)) && deflection>0 ){
	    sel[5]++;
	    hrig2->Fill(rigidity);
	    float m=3.725;
	    float a=4.;
	    float z=2.;
	    float en=(sqrt( (z*rigidity)*(z*rigidity)+m*m )-m)/a;
	    hen2->Fill(en);
	}

	//--------
	// Z=1
	//--------
	if(dedxtrk > (2.8+3.3/(rigidity*rigidity)))return;
	if(dedxtrk < (0.4+0.35/(rigidity*rigidity)))return;
	sel[2]++;

	if(deflection>0){
	    hrig1->Fill(rigidity);
	    float m=0.938;
	    float a=1.;
	    float z=1.;
	    float en=(sqrt( (z*rigidity)*(z*rigidity)+m*m )-m)/a;
	    hen1->Fill(en);
	}
	//------------------------
	// 6-planes lever-arm
	//------------------------
	if( !(track->GetTrkTrack()->XGood(0) && track->GetTrkTrack()->XGood(5)) )return;
	sel[3]++;

	for(Int_t ip=0; ip<6; ip++)herrx[ip]->Fill(errx[ip]);
	for(Int_t ip=0; ip<6; ip++)herry[ip]->Fill(erry[ip]);

	//---------------------------
	// cut on spatial resolution
	//---------------------------
	if( errx[0] > 0.0004 )return;
	if( errx[5] > 0.0004 )return;
	sel[4]++;

	herrdef->Fill(errdef);
	hbetavsdef2->Fill(deflection,beta);
	hdef->Fill(deflection);
    }
    
//    delete track;
   
    
}
//===============================================================
// Save histograms
//===============================================================
void SaveHistos( TFile * outf ){

    TString hid; 

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

//  ----------------------------------
//  save on file
//  ----------------------------------
    outf->cd();
    
    hchi2def->Write();
    hchi2rig->Write();
    hdedxtofvsdedxtrk->Write(); 
    hdedxtrkvsrig->Write();   
    hdedxtofvsrig->Write(); 
    hbetavsrig->Write();  
    hdedxtrkvsdef->Write(); 
    hbetavsdef->Write(); 
    hbetavsdef2->Write(); 
    hdef->Write(); 
    hrig1->Write(); 
    hrig2->Write(); 
    hen1->Write(); 
    hen2->Write(); 
    herrdef->Write(); 
    for(Int_t ip=0; ip<6; ip++){
	herrx[ip]->Write(); 
	herry[ip]->Write(); 
    }

    cout << "------------------------------------------------------------------" <<endl;
    cout << "Pre-selected events :"<<tot<<endl;
    if(!tot)return;
    cout << "chi2 selection                            :"<<sel[1]<< "  (" << ((float)sel[1]/(float)tot) <<")"<<endl;
    cout << "Z=2 selection                             :"<<sel[5]<< "  (" << ((float)sel[5]/(float)tot) <<")"<<endl;
    cout << "Z=1 selection                             :"<<sel[2]<< "  (" << ((float)sel[2]/(float)tot) <<")"<<endl;
    cout << "(6-plane lever-arm)                       :"<<sel[3]<< "  (" << ((float)sel[3]/(float)tot) <<")"<<endl;
    cout << "ext.spatial res. on planes 1-6 < 4 micron :"<<sel[4]<< "  (" << ((float)sel[4]/(float)tot) <<")"<<endl;
    cout << "------------------------------------------------------------------" <<endl;

}