#if !defined(__CINT__) || defined(__MAKECINT__)
#include <TString.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TF1.h>
#include <TFile.h>
#include <TEventList.h>
#include <TMultiGraph.h>
#include <TCanvas.h>
#include <TCut.h>

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

#include <PamLevel2.h>
#include <TrkNuclei.h>

#endif

//-------------------------------------------------
// functions prototypes
//-------------------------------------------------
float GetDEDX_ToFPlane(int plane, ToFTrkVar *tof);
int   GetHitPMTsOutsideTrack(ToFLevel2* l2,ToFTrkVar* tof,int plane);
bool  ToFOK_New(ToFLevel2* tofl2, ToFTrkVar *tof);
int   GetClustersOusideTrack(int iv, int ip, TrkLevel2 *trkl2, TrkTrack *trk, float lr);
int   GetPlanesWithClustersOutsideTrack(TrkLevel2 *trkl2, TrkTrack *trk, float lr);
bool  MultipleTracksInsideTracker(TrkLevel2* trkl2, TrkTrack *trk, ToFLevel2 *tofl2, ToFTrkVar *tof);
//--------------------------------------------------------------------------------
//
//
//
//
//--------------------------------------------------------------------------------
TH2F* test(TString ddir, TString list){

    PamLevel2 *pam = new PamLevel2(ddir,list,"+AUTO");
    TChain *chain = pam->GetPamTree();

    TrkNuclei *trknuc = new TrkNuclei();
    TH2F *h = new TH2F("h","Z (trk) vs beta",100,0.3,1.1,100,0.,9.);
    TH2F *hh = new TH2F("hh","Z (trk) vs rigidity",100,0.1,1000,100,0.,9.);
    TH2F *hhh = new TH2F("hhh","Z vs Z",100,0.,9.,100,0.,9.);


    ULong64_t sel=0;
    ULong64_t tot=0;
// --------------------
// loop over the events
// --------------------
    for(int iev=0; iev<chain->GetEntries(); iev++){
//    for(int iev=0; iev<50000; iev++){

	// read the event
	pam->GetEntry(iev);
	tot++;

	// retrieve pamela objects
	ToFLevel2 *tofl2 = pam->GetToFLevel2();
	TrkLevel2 *trkl2 = pam->GetTrkLevel2();

	ToFTrkVar *tof = pam->GetTrack(0)->GetToFTrack();
	TrkTrack  *trk = pam->GetTrack(0)->GetTrkTrack();

	// -------------------------------------------------------
	// these are part of the cuts to reject interacting events
	// developed for antiprotons (slightly modified... search XXX)
	// -------------------------------------------------------
	bool INTERACTION = true;
	if( 
	    ToFOK_New(tofl2,tof) &&
	    !MultipleTracksInsideTracker(trkl2,trk,tofl2,tof) &&
	    true)INTERACTION=false;
	if(INTERACTION)continue; 
	

	// -------------------------------------------------------
	// retrieve the beta (Wolfgang)
	// -------------------------------------------------------
	// NB: this loop over the stored tracks IS NOT in the intention of 
	// the tracker classes implementation.
	// CalcBeta should be a method of the ToFTrkVar class and should
	// not need "notrack" as input. 
	int notrack;
	for(notrack=0; notrack<tofl2->ntrk(); notrack++){
	    if( ((ToFTrkVar*)(tofl2->ToFTrk->At(notrack)))->trkseqno == trk->seqno )break;
	};
	float b2 = pam->GetToFLevel2()->CalcBeta(notrack,5.,15.,4.);  // medium quality
	float b3 = pam->GetToFLevel2()->CalcBeta(notrack,3.,20.,3.);  // high quality

	if( b3<=0. || b3>2. )continue;

	// -------------------------------------------------------
	// now fill the TrkNuclei class
	// -------------------------------------------------------
	trknuc->Reset();
	if( !trknuc->Set(trk) )continue;

	sel++;

	// -------------------------------------------------------
	// the class contains a copy of the TrkTrack object
	// so that you can retrieve any TrkTrack variable/methods
	// as trknuc->GetTrkTrack()->....
	// -------------------------------------------------------
	float dedx = trknuc->GetTrkTrack()->GetDEDX();	
	// -------------------------------------------------------
	// the methods to get the charge (as a function of beta) are
	//
	//  TrkNuclei::GetZ_Beta(ip,iv,beta) // single view (x or y)
	//  TrkNuclei::GetZ_Beta(ip,beta)    // single plane (average between x and y)
	//  TrkNuclei::GetZ_Beta(beta)       // average over all planes
	//
	// -------------------------------------------------------
	float z1 = trknuc->GetZ_Beta( b3 );
	// -------------------------------------------------------
	// the calibration mip-to-charge is implemented as a 
	// TF2 static object, function of beta and dedx. 
	// You can access directly this function as follows:
	// -------------------------------------------------------
	float z2 = TrkNuclei_parameters::Get()->GetZ_Beta()->Eval(b3,dedx);
	// -------------------------------------------------------
	// or you can change the calibration function as you like, 
	// with the following methods (must be called once):
	//
	// TrkNuclei_parameters::Get()->SetZ_Beta(TF2 *myfunc);
	//
	// -------------------------------------------------------


	// -------------------------------------------------------
	// Also implemented is the mip-to-charge as a function 
	// of rigidity. 
	//
	//  TrkNuclei::GetZ_Rigidity(ip,iv,rig) // single view (x or y)
	//  ecc...
	//
	// At present however the TF2 object is NULL. 
	// You can set your own function as explained above 
	//
	// TrkNuclei_parameters::Get()->SetZ_Rigidity(TF2 *myfunc);
	//
	// -------------------------------------------------------
	float rig = trknuc->GetTrkTrack()->GetRigidity();

	float z3  = trknuc->GetZ_Rigidity( rig );
	float z4  = TrkNuclei_parameters::Get()->GetZ_Rigidity()->Eval(rig,dedx);

	h->Fill(b3,z1);
	hh->Fill(rig,z3);
	hhh->Fill(z1,z3);

	cout << setw(10) << b3 <<setw(10) << z1 <<setw(10)<<1./rig<<setw(10)<< z3 << endl;

       
    }

    cout << endl ;
    cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-"<<endl;
    cout << sel<<" selected events over "<<tot<<endl;
    cout << "+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-"<<endl;

//    return h;
//    return hh;
    return hhh;

}
/////////////////////////////////////////////////////////////////////////////
//
//
// Selection function from antiproton analysis (to be reviewed!!!!)
//
//
/////////////////////////////////////////////////////////////////////////////
//----------------------------------------------------------------------------
//
//----------------------------------------------------------------------------
float GetDEDX_ToFPlane(int plane, ToFTrkVar *tof){

    if(!tof)return 0.;
    float dedx = 0.;
    float ndedx = 0;
    ToFLevel2 l2 = ToFLevel2();
    for(int ip=0; ip<tof->npmtadc; ip++){

        Int_t iplane  = -1;
        Int_t ipaddle = -1;
        Int_t ipmt    = -1;
        Int_t id      = tof->pmtadc[ip];
	
	l2.GetPMTName(id,iplane,ipaddle,ipmt);

	if(iplane!=plane)continue;
        dedx  += ( tof->adcflag[ip]==0 )*tof->dedx[ip];
	ndedx += ( tof->adcflag[ip]==0 );
    }    
    if(ndedx>0)return dedx/ndedx;
    return 0.;
}
//--------------------------------------------------------------------------------
//
// conta il numero di PMT con segnale (tdc) fuori dalla traccia
//
//--------------------------------------------------------------------------------
int GetHitPMTsOutsideTrack(ToFLevel2* l2,ToFTrkVar* tof,int plane){
    int nn=0;
//    cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%"<<endl;
    for(int ip=0; ip<l2->npmt(); ip++){
        Int_t iplane  = -1;
        Int_t ipaddle = -1;
        Int_t ipmt    = -1;
        Int_t id      = l2->GetToFPMT(ip)->pmt_id;
	l2->GetPMTName(id,iplane,ipaddle,ipmt);

	if(iplane!=plane)continue;
//	cout << endl<<plane << " -- "<<id;
        float tdc = l2->GetToFPMT(ip)->tdc; if(tdc >= 4095.)continue;///to avoid adc pile-up

	int iht = -1;
	int iha = -1;
	for(iht=0; iht<tof->npmttdc; iht++)if(id==tof->pmttdc[iht])break;
	if( iht>=0 && iht <tof->npmttdc )continue;
	for(iha=0; iha<tof->npmtadc; iha++)if(id==tof->pmtadc[iha])break;
	if( iha>=0 && iha <tof->npmtadc )continue;
	nn++;
//	cout << " OUTSIDE ";
    }    
    return nn;
}
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//
//
//
//
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//
// it requires:
// - not more than 1 hit paddle on S11/S12/S21/S22
// - at least 1 hit paddle on S1/S2
// - if there is a hit paddle (S1/S2), this must be associated to the track
// 
//--------------------------------------------------------------------------------
bool ToFOK_New(ToFLevel2* tofl2, ToFTrkVar *tof){
    
    if( !tofl2)return false;
    
//    cout << "((( TOFOK )))"<<endl;

    bool TOF__OK=false;

    // -----------------------------------------
    // at least one hit paddle on both S1 and S2
    // -----------------------------------------
    if(	
	//-----------------------------------
	// not more than 1 paddle per layer 
	//-----------------------------------
	tofl2 &&	
	tofl2->GetNHitPaddles(0) <= 1 && //S11
	tofl2->GetNHitPaddles(1) <= 1 && //S12
	tofl2->GetNHitPaddles(2) <= 1 && //S21
	tofl2->GetNHitPaddles(3) <= 1 && //S22
//	tofl2->GetNHitPaddles(4) <= 1 && //S31
//	tofl2->GetNHitPaddles(5) <= 1 && //S32
	//-----------------------------------
	// 1 paddle per layer 
	//-----------------------------------
// 	tofl2->GetNHitPaddles(0) == 1 && //S11
// 	tofl2->GetNHitPaddles(1) == 1 && //S12
// 	tofl2->GetNHitPaddles(2) == 1 && //S21
// 	tofl2->GetNHitPaddles(3) == 1 && //S22
//	tofl2->GetNHitPaddles(4) == 1 && //S31
//	tofl2->GetNHitPaddles(5) == 1 && //S32
	//-----------------------------------
	// at least 1 paddle per plane
	//-----------------------------------
	(tofl2->GetNHitPaddles(0)==1 || tofl2->GetNHitPaddles(1)==1) && //S1
	(tofl2->GetNHitPaddles(2)==1 || tofl2->GetNHitPaddles(3)==1) && //S2
//	(tofl2->GetNHitPaddles(4)==1 || tofl2->GetNHitPaddles(5)==1) && //S3
	//-----------------------------------
	// n.hit pmts outside track
	//-----------------------------------
//  	GetHitPMTsOutsideTrack(tofl2,tof,0)<2 && //S11	 //troppo (puo` succedere )
//  	GetHitPMTsOutsideTrack(tofl2,tof,1)<2 && //S12	
//  	GetHitPMTsOutsideTrack(tofl2,tof,2)<2 && //S21	
//  	GetHitPMTsOutsideTrack(tofl2,tof,3)<2 && //S22	
 	GetHitPMTsOutsideTrack(tofl2,tof,0)<3 && //S11	 
 	GetHitPMTsOutsideTrack(tofl2,tof,1)<3 && //S12	
// 	GetHitPMTsOutsideTrack(tofl2,tof,2)<3 && //S21	
// 	GetHitPMTsOutsideTrack(tofl2,tof,3)<3 && //S22	
	true ) TOF__OK = true;
//      cout << "////"<<tofl2->GetNHitPaddles(0)<<endl;
//      cout << "////"<<tofl2->GetNHitPaddles(1)<<endl;
//      cout << "////"<<tofl2->GetNHitPaddles(2)<<endl;
//      cout << "////"<<tofl2->GetNHitPaddles(3)<<endl;
//      cout << "////"<<tofl2->GetNHitPaddles(4)<<endl;
//      cout << "////"<<tofl2->GetNHitPaddles(5)<<endl;

    if( !TOF__OK )return false;
        
// ----------------------------------------
// dopo che wolfgang mi ha illuminato....
// ----------------------------------------    
    int hitplane[]={0,0,0,0,0,0};
    int hitpaddle[]={-1,-1,-1,-1,-1,-1};
    for(Int_t ip=0; ip< tof->npmttdc; ip++){
        Int_t iplane  = -1;
        Int_t ipaddle = -1;
        Int_t ipmt    = -1;
        Int_t id      = tof->pmttdc[ip];
	tofl2->GetPMTName(id,iplane,ipaddle,ipmt);
	if(tof->tdcflag[ip]==0){
	    hitplane[iplane]=1; //there is a true tdc signal associated to the track
	    if( hitpaddle[iplane]>=0 && hitpaddle[iplane]!=ipaddle )cout<< "ORRORE!!!!"<<endl;
	    hitpaddle[iplane]=ipaddle;//store the id of the paddle associated to the track
	}
    };
    //
    for(int iplane=0; iplane<4; iplane++){//loop over S11/S12/S21/S22
	//retrieve the id of the paddle traversed by the track
	// (...there is some tolerance...)
	int ipaddle = tofl2->GetPaddleIdOfTrack(tof->xtr_tof[iplane],tof->ytr_tof[iplane],iplane);
	//check if the traversed paddle is hit
	bool OK = true;
	if(
	    tofl2->GetNHitPaddles(iplane)>0 &&//if there is a hit paddle in this plane...
	    !tofl2->HitPaddle(iplane,ipaddle) && //...and the paddle traversed by the track is not hit	    
	    true)OK = false;//..discard the event
	/// hence try to recover some events...
	if( 
	    tofl2->GetNHitPaddles(iplane)>0 &&//if there is a hit paddle...
	    ipaddle==-1 &&//...and the track does not traverse any paddle...
	    tofl2->HitPaddle(iplane,hitpaddle[iplane]) &&//... BUT there are tdc signals belonging to a hit paddle
	    true)OK=true;//recover
	TOF__OK = TOF__OK && OK;
    }    

    if( !TOF__OK )return false;


    return true;

}
///////////////////////////////////////////////////////////
// conta i cluster fuori traccia sulla singola vista
// se lr=+1(-1) conta solo quelli a destra(sinistra) della traccia
// vengono esclusi dal conteggio: 
// - i cluster bad
// - i cluster saturati (...misembrava ragionevole)
// - i cluster entro toll dalla traccia
// - i cluster con signal < dedxcut
//
int GetClustersOusideTrack(int iv, int ip, TrkLevel2 *trkl2, TrkTrack *trk, float lr){

    if(!trkl2)return -1;
    if(ip<-1||ip>5)return -1;
    if(lr!=0 && trk==0)return -1;


    float dedxcut = 0.7;//mip
    float toll    = 0.01;//cm (100 micron)

    int nhits=0;
    int nsing=0;
    
    bool LEFT = false;
    bool RIGHT = false;
    if(lr>0)RIGHT=true;
    if(lr<0)LEFT=true;

    //
    if(iv==0)nsing=trkl2->nclsx();
    if(iv==1)nsing=trkl2->nclsy();

    float trkcoord = 0.;
    if(trk && iv==0)trkcoord=trk->xv[ip];
    if(trk && iv==1)trkcoord=trk->yv[ip];

    //
    for(int ii=0; ii<nsing; ii++){
	TrkSinglet *s = 0;
	if(iv==0)s = trkl2->GetSingletX(ii);
	if(iv==1)s = trkl2->GetSingletY(ii);
	if(!s)continue;
	if(s->plane != ip+1)continue;
//	cout << ii << s->plane << s->IsBad() << endl;
//	if(!s->IsBad())cout << s->sgnl;
//	cout << s->plane << " - "<< (s->coord[0] - trkcoord) << "  "<<(s->coord[1] - trkcoord)<<" "<<s->sgnl<<endl;
	bool SIDEOK=true;
	if( 
	    LEFT && 
	    ((s->coord[0] - trkcoord)> -1*toll 
	     || 
	     (s->coord[1] - trkcoord)> -1*toll) &&
	    true)SIDEOK=false; ;
	if( 
	    RIGHT && 
	    ((s->coord[0] - trkcoord)< +1*toll 
	     ||
	     (s->coord[1] - trkcoord)< +1*toll) &&
	    true)SIDEOK=false; ;

	if( 
//	    (ip==-1 || s->plane == ip+1) &&
	    !s->IsBad()      &&
	    s->GetSignal() > dedxcut &&
	    !s->IsSaturated() &&
	    SIDEOK &&
	    true){
//	    cout << "*";	    
	    nhits++;	    
	} 
///	cout << endl;
    }
//    cout << "-----"<<endl;
    return nhits;

}
//----------------------------------------------------------------------------
//
//----------------------------------------------------------------------------
int GetPlanesWithClustersOutsideTrack(TrkLevel2 *trkl2, TrkTrack *trk, float lr){
    int np=0;
    for(int ip=0; ip<6; ip++)
 	np+=(GetClustersOusideTrack(0,ip,trkl2,trk,lr)>0||GetClustersOusideTrack(1,ip,trkl2,trk,lr)>0);
    return np;
}
///////////////////////////////////////////////////////////
bool MultipleTracksInsideTracker(TrkLevel2* trkl2, TrkTrack *trk, ToFLevel2 *tofl2, ToFTrkVar *tof){
    

    if(!trkl2)return true;
    if(!tofl2)return true;

    if(!trk)return true;
    if(!tof)return true;

//    cout << "MultipleTracksInsideTracker ?"<<endl;
    float rig = trk->GetRigidity();
    
    // n.PMTs outside the track
    int npmt = 0;
    if(tofl2&&tof)
	for(int ip=0; ip<4; ip++)npmt+=GetHitPMTsOutsideTrack(tofl2,tof,ip);
    
    // n.hits on the left/right of the track (x-view)
    int npxl=0;
    int npxr=0;
    for(int ip=0; ip<5; ip++)
	npxl+=(GetClustersOusideTrack(0,ip,trkl2,trk,-1)>0) ;
    for(int ip=0; ip<5; ip++)
	npxr+=(GetClustersOusideTrack(0,ip,trkl2,trk,+1)>0) ;
    
    // n.hits outside the track (y-view)
    int npy=0;
    for(int ip=0; ip<5; ip++)
	npy+=(GetClustersOusideTrack(1,ip,trkl2,trk,0)>0) ;

    // upper limit on tof dedx
    bool troppo1=false;
    bool troppo2=false;
    //
    //
//  ------------------------------------------------------
//   XXX
//  commented to adapt the routine to nuclei analysis !!!
//  ------------------------------------------------------
//
//     float cuts11up = F_ftofdedxcut(0,rig);
//     float cuts12up = F_ftofdedxcut(1,rig);
//     if( tofl2 && GetDEDX_ToFPlane(0,tof) > cuts11up)troppo1=true;
//     if( tofl2 && GetDEDX_ToFPlane(1,tof) > cuts12up)troppo1=true;
//     //
//     float cuts21up = F_ftofdedxcut(2,rig);
//     float cuts22up = F_ftofdedxcut(3,rig);
//     if( tofl2 && GetDEDX_ToFPlane(2,tof) > cuts21up)troppo2=true;
//     if( tofl2 && GetDEDX_ToFPlane(3,tof) > cuts22up)troppo2=true;

    // multiple tracks
    bool INTERACTING = false;
    if( (npxl>2 || npxr>2) && npy>2)INTERACTING=true;
    if( (npxl>1 || npxr>1) && npy>1 && (npmt!=0||troppo2) )INTERACTING=true;

//     cout << npxl << " "<<npxr<< " "<<npy<<endl;
//     cout << npmt << endl;
//     if( (npxl>2 || npxr>2) && npy>2)cout << " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 33"<<endl;
//     if( (npxl>1 || npxr>1) && npy>1 && npmt!=0 )cout << " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% npmt !!"<<endl;
//     if( (npxl>1 || npxr>1) && npy>1 && troppo2 )cout << " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% troppo2 !!"<<endl;
    
    return INTERACTING;
    
}