#include "EventMixer.h"
#include <TH2D.h>
#include <TH2.h>
#include <TVector3.h>
#include <TLorentzVector.h>
#include "TMatrixD.h" 
#include "TMatrixDSym.h" 
#include "TMatrixDSymEigen.h" 
#include <iostream>
#include <vector>
#include <math.h>
#include <string>
#include "myHist.C"
using namespace std;

// track quality cuts
static const int nTrackCut(6) ;

static const Int_t ndof_cut(0) ;
static const Int_t nchi2_cut(5) ;
//static const Double_t dxy_cut(1.) ;
static const Double_t dxy_cut(.15) ;
static const Double_t rmin_cut(20.) ;
static const Double_t pt_cut(0.2) ;
static const Double_t eta_cut(2.4) ;
// split track cuts
static const Double_t cst_split_cut(0.99996) ;
static const Double_t dpt_split_cut(0.04) ;

static const Double_t pimass(0.1396), kamass(0.4937) ; 
static const int MaxEvt( 1500000 ) ;

static const Long64_t offset = (Long64_t) pow(10,7) ;

bool Used[MaxEvt] = { MaxEvt * 0 } ;

vector <bool> GoodTrack ; // list of good tracks
vector <bool> SplitTrack ; // list of split tracks
vector <TLorentzVector> Tk4V ; // track four vectors (pi mass hyp.)
vector <TLorentzVector> TkGood ; // track four vectors (pi mass hyp.)
vector <int> ChGood ;
vector <Long64_t> EventSorter ;
vector <int> SortedList ;
vector <int> SortedCont ;
EventMixer::EventMixer() 
{}

EventMixer::~EventMixer() 
{}
// =======================================================================
void EventMixer::SelectPair( int entry, int & entry1, int & entry2 ) 
{
  if( (uint) entry >= SortedList.size() )    
    { entry1 = -1 ; entry2 = -1 ; return ;}

  entry1 = SortedList.at( entry )  ;
  entry2 = SortedList.at( entry+1 );

  if( SortedList.size()-entry < 10 ) 
  { cout << " sorted list " << entry << " " << entry1 << " " << entry2 
	 << " " << SortedCont.at(entry) << " " << SortedCont.at(entry+1) << endl ; }
}
// =======================================================================
void EventMixer::MixEvents( vector <int> ch1, vector <TLorentzVector> tkv1, 
			    vector <int> ch2, vector <TLorentzVector> tkv2 ) 
{
  Hfill( "Vsize", tkv1.size(), 100, 0.,100., tkv2.size(), 100, 0.,100., 1. ) ;

  // signal (1st)
  for( uint i=0 ; i< tkv1.size() ; i++ ) {
    for( uint j=i+1 ; j< tkv1.size() ; j++ ) 
      {
	double q = GetQ( tkv1.at(i) , tkv1.at(j) ) ;
	if( ch1.at(i) == ch1.at(j) ) 
	  { 
	    Hfill("QsigSameCharge",q,200,0.,2.,1.) ; 
	    Hfill("QsigSameChargeCC",q,200,0.,2.,CoulombW(q) ) ;
	  }
	else
	  { 
	    Hfill("QsigOppoCharge",q,200,0.,2.,1.) ;
	    Hfill("QsigOppoChargeCC",q,200,0.,2.,CoulombWpm(q) ) ; 
	  }
      }
  }      

  // signal (2nd)
  for( uint i=0 ; i< tkv2.size() ; i++ ) {
    for( uint j=i+1 ; j< tkv2.size() ; j++ ) 
      {
	double q = GetQ( tkv2.at(i) , tkv2.at(j) ) ;
	if( ch2.at(i) == ch2.at(j) ) 
	  { 
	    Hfill("QsigSameCharge",q,200,0.,2.,1.) ; 
	    Hfill("QsigSameChargeCC",q,200,0.,2.,CoulombW(q) ) ;
	  }
	else
	  { 
	    Hfill("QsigOppoCharge",q,200,0.,2.,1.) ;
	    Hfill("QsigOppoChargeCC",q,200,0.,2.,CoulombWpm(q) ) ; 
	  }
      }
  }      

  // mix (1st+2nd)
  for( uint i=0 ; i< tkv1.size() ; i++ ) {
    for( uint j=0 ; j< tkv2.size() ; j++ ) 
      {
	double q = GetQ( tkv1.at(i) , tkv2.at(j) ) ;
	Hfill("Qmix",q,200,0.,2.,1. ) ; 
	if( ch1.at(i) == ch2.at(j) ) 
	  { Hfill("QmixSameCharge",q,200,0.,2.,1. ) ; }
	else
	  { Hfill("QmixOppoCharge",q,200,0.,2.,1. ) ; }
      }
  }
}

// =======================================================================
void EventMixer::analyze( int entry , vector<int>& chv ,vector <TLorentzVector>& tkv )
{
  analyze( entry ) ;
  tkv = TkGood ;
  chv = ChGood ;
  return; 
}
// =======================================================================
Long64_t EventMixer::Encode( int w )
{
  int nbins(10) ;
  int Range[]={3,5,8,10,13,16,20,25,30,200} ;
  int i(0), j(0) ;
  while( w >= Range[i++] ) j++ ;
  Hfill("enc",(double) w,50,0.,50., (double) j , nbins, 0., (double) nbins, 1. ) ;
  return (Long64_t) j ;
}
// =======================================================================
Long64_t EventMixer::GetMultEtaRange() 
{
  int nFwd(0), nBkw(0), nCtr(0) ;
  for( uint i=0 ; i< TkGood.size() ; i++ ) 
    {
      double eta = (TkGood.at(i)).Eta() ;
      if( eta < -0.8 ) 
	{ nBkw++ ;}
      else if( eta > 0.8 ) 
	{ nFwd++ ;}
      else
	{ nCtr++ ;}
    }
  Hfill("nFwd",(double) nFwd, 150, 0.,150.,1. ) ;
  Hfill("nBkw",(double) nBkw, 150, 0.,150.,1. ) ;
  Hfill("nCtr",(double) nCtr, 150, 0.,150.,1. ) ;
  Hfill("nTotGood",(double) TkGood.size(), 150, 0., 150.,1.) ;
 
  Long64_t ReturnValue = (100*Encode( nFwd ) + 10 * Encode( nCtr ) + Encode( nBkw )) ;
  return ReturnValue ;
}
// =======================================================================
vector <int> EventMixer::EvtSort( )
{
  cout << " Enter Event Sorter " << endl ;
  sort( EventSorter.begin(), EventSorter.end() ) ;  

  vector <int> ListSorted ;
  for( uint i = 0 ; i< EventSorter.size() ; i++ ) 
    {
      int k = EventSorter.at(i)%offset ;
      int l = EventSorter.at(i)/offset ;
      SortedList.push_back(k) ;
      SortedCont.push_back(l) ;
      ListSorted.push_back(k) ;
      // cout << strg ;
    }
  cout << " Events are sorted ! " << endl ;
  ofstream OutputSortedList("SortList.txt",ios::out) ;
  for( uint i = 0 ; i< SortedList.size() ; i++ ) 
    {
      char strg[100] ;
      int s = sprintf(strg,"%10d %10d %10d \n",i,SortedList.at(i),SortedCont.at(i)  ) ;
      OutputSortedList << strg ;
    }
  OutputSortedList.close() ;
  return ListSorted ;

}
// =======================================================================
void EventMixer::analyze( int entry ) {

  MC=true;
  if (runNumber>124000 && runNumber<125000) MC=false;
  if( numberOfTrack > 150 ) return ;   // sanity cut 

  SelectGoodTracks() ;  // as defined in GoodTrackNoDofCut
  RemoveSplitTracks() ; 
  FillGoodTracks() ;

  if( entry<0 ) return ;

  // sort by sphericity
  ComputeSphericity() ;
  return ;

  // sort by eta range multiplicity
  Long64_t MultEtaRange = GetMultEtaRange() ; 
  MultEtaRange = MultEtaRange* offset + entry ;
  EventSorter.push_back( MultEtaRange ) ;
  return ;

  // sort by event mass 
  Long64_t EvtMas = min( (int) ComputeEventMass() ,100 ) ;
  EvtMas = EvtMas*offset + entry ;
  EventSorter.push_back( EvtMas ) ;

  //AnalyzeTracks() ; 
  //ComputeSphericity() ;
  //ComputeBoseEinsteinMC() ; 
  //ComputeBoseEinstein() ; 

  return;
}
// =======================================================================
double EventMixer::ComputeEventMass()
{
  TLorentzVector Evt4V(0,0,0,0) ;
  for( uint i=0 ; i< TkGood.size() ; i++ ) 
    {
      Evt4V += TkGood.at(i) ;
    }
  Hfill("MassVsN",Evt4V.Mag(),100,0.,100.,(double)TkGood.size(),100.,0.,100.,1.) ;
  return Evt4V.Mag() ;
}
// =======================================================================
void EventMixer::ComputeSphericity() 
{
  vector <TLorentzVector> TkV ;
  for( int i=0 ; i < numberOfTrack ; i++ )
    {
      if( !GoodTrack.at(i) ) continue ;
      TkV.push_back( Tk4V.at(i) ) ;
    }
  TMatrixD SphM(3,3) ;
  double_t sph = GetSphericity( TkV , &SphM ) ;
  Hfill( "SphVsNtk",TkV.size(),100,0.,100., sph,100,0.,1.,1. ) ;

  const TVector3 SphV3( SphM[0][0],SphM[1][0],SphM[2][0] ) ;

  // look at vectors in the same and in the opposite hemisphere wrt sph axis
  vector <TLorentzVector> TkVSame ;
  vector <TLorentzVector> TkVOppo ;
  for( unsigned int i=0 ; i < TkV.size()  ; i++ )
    {
      TVector3 TkV3 = (TkV.at(i)).Vect() ;
      double cst = cos( TkV3.Angle( SphV3 ) ) ;
      if( cst > 0 ) 
	{ TkVSame.push_back( TkV.at(i) ) ; }
      else
	{ TkVOppo.push_back( TkV.at(i) ) ; }
      //      if( TkV.size() == 1 ) cout << " single " << cst << "\n" ;
      Hfill("Cst",cst,202,-1.01,1.01,1.) ;
    }

  TMatrixD SphMS(3,3) ;
  double_t sphs = GetSphericity( TkVSame , &SphMS ) ;
  TMatrixD SphMO(3,3) ;
  double_t spho = GetSphericity( TkVOppo , &SphMO ) ;
  //  cout << sphs << " " << TkVSame.size() << " " << spho << " " << TkVOppo.size() << endl ;
  Hfill( "SphSVsSphO",spho,100,0.,1., sphs,100,0.,1.,1. ) ;
  double Phi = atan2( SphM[1][0], SphM[0][0] ) ;
  double Pxy = sqrt(  SphM[0][0]*SphM[0][0] + SphM[1][0]*SphM[1][0] ) ;
  double Theta = atan2( SphM[2][0], Pxy ) ;
  Hfill( "STP",Phi,100,-3.2,3.2,Theta,100,0.,3.2,1.) ;
  return ;
}
// =======================================================================
double EventMixer::GetSphericity( vector <TLorentzVector> TrackV, TMatrixD *SphEigenMatrix )
{
  double SumP2(0) ;
  double SphT[]={0,0,0,0,0,0,0,0,0} ;
  for( unsigned int i=0 ; i< TrackV.size() ; i++ )
    {
      TLorentzVector V = TrackV.at(i) ;
      SphT[0] += V.Px()*V.Px() ;
      SphT[1] += V.Px()*V.Py() ;
      SphT[2] += V.Px()*V.Pz() ;

      SphT[3] += V.Py()*V.Px() ;
      SphT[4] += V.Py()*V.Py() ; 
      SphT[5] += V.Py()*V.Pz() ; 

      SphT[6] += V.Pz()*V.Px() ;
      SphT[7] += V.Pz()*V.Py() ; 
      SphT[8] += V.Pz()*V.Pz() ;
 
      SumP2 += V.P() * V.P() ;
    }

  if( SumP2 == 0 ) return -1 ;

  TMatrixDSym SphM( 3, SphT ) ;

  SphM *= 1./SumP2 ;
  TMatrixDSymEigen SphE( SphM ) ;
  *SphEigenMatrix = SphE.GetEigenVectors() ;
  TVectorD SphEigenValues = SphE.GetEigenValues() ;
  double sph = 1.5 * (SphEigenValues[1]+SphEigenValues[2]);

  return sph ;

}

// =======================================================================
void EventMixer::AnalyzeTracks() { 
  int goodTracks=0;
  for( int i = 0 ; i < numberOfTrack ; i++ ) {
    //cout << " Genp_idx_clst " << Genp_idx_clst->at(i) << endl;
    Hfill("trk_charge",charge->at(i),4,-1.5,1.5,1. ) ;
    Hfill("trk_pt",pt->at(i),200,0.,50.,1. ) ;
    Hfill("trk_p",p->at(i),200,0.,50.,1. ) ;
    Hfill("trk_eta",eta->at(i),200,-3.,3.,1. ) ;
    Hfill("trk_phi",phi->at(i),200,-4.,4.,1. ) ;
    Hfill("trk_dxy",dxy_PV->at(i),200,0.,0.5,1. );
    Hfill("trk_dz",dz->at(i),200,-20.,20.,1. ) ;
    Hfill("trk_d0",d0->at(i),200,0.,0.5,1. );
    Hfill("trk_chi2",normalizedChi2->at(i),200,0.,10.,1. ) ;
    //Hfill("trk_dedx",olddEdx1->at(i),200,0.,10.,1. ) ;
    Hfill("trk_validHits",numberOfvalidHits->at(i),50,0.,50.,1. ) ;
    Hfill("trk_LostHits",highPurity_track->at(i),50,0.,50.,1. ) ;
    if( !GoodTrack.at(i) ) continue ;
    goodTracks++;
    Hfill("goodtrk_charge",charge->at(i),4,-1.5,1.5,1. ) ;
    Hfill("goodtrk_pt",pt->at(i),200,0.,50.,1. ) ;
    Hfill("goodtrk_p",p->at(i),200,0.,50.,1. ) ;
    Hfill("goodtrk_eta",eta->at(i),200,-3.,3.,1. ) ;
    Hfill("goodtrk_phi",phi->at(i),200,-4.,4.,1. ) ;
    Hfill("goodtrk_dxy",dxy_PV->at(i),200,0.,0.5,1. );
    Hfill("goodtrk_dz",dz->at(i),200,-20.,20.,1. ) ;
    Hfill("goodtrk_d0",d0->at(i),200,0.,0.5,1. );
    Hfill("goodtrk_chi2",normalizedChi2->at(i),200,0.,10.,1. ) ;
    //Hfill("trk_dedx",olddEdx1->at(i),200,0.,10.,1. ) ;
    Hfill("goodtrk_validHits",numberOfvalidHits->at(i),50,0.,50.,1. ) ;
    Hfill("goodtrk_LostHits",highPurity_track->at(i),50,0.,50.,1. ) ;
    
  }
  Hfill("GoodTracks",goodTracks,150,0.,150.,1. ) ;
}

void EventMixer::ComputeBoseEinsteinMC() { 
  for( int i = 0 ; i < numberOfTrack ; i++ ) {
    double eta_t=eta->at(i);
    double phi_t=phi->at(i);
    double Rmin=9999;
    int igen=0;
    for( unsigned int j = 0 ; j < Genp_status->size() ; j++ ) {
      if (Genp_charge->at(j)!=charge->at(i)) continue;
      double eta_g=Genp_eta->at(j);
      double phi_g=Genp_phi->at(j);
      double R=sqrt((eta_t-eta_g)*(eta_t-eta_g)+(phi_t-phi_g)*(phi_t-phi_g));
      if (R<Rmin) {
        Rmin=R;
        igen=j;
      }
    }
    //cout << " Genp_idx_clst " << igen << " id " << Genp_Id->at(igen) << " R " << Rmin << endl;
    Hfill("Rmin_gen_reco",Rmin,200,0.,1.,1. ) ;
    if (Rmin<0.03) (*Genp_idx_clst)[i]=igen;
    else (*Genp_idx_clst)[i]=0;
  }
}

// =======================================================================
// Compute B.E. correlations for the siganl and some control samples
// =======================================================================
void EventMixer::ComputeBoseEinstein() { 
  // for( int i = 0 ; i < Genp_idx_clst->size() ; i++ ) {
  //   cout << " Genp_idx_clst " << Genp_idx_clst->at(i) << endl;
  // }
  for( int i = 0 ; i < numberOfTrack ; i++ ) {
    //cout << " Genp_idx_clst " << Genp_idx_clst->at(i) << endl;
    if( !GoodTrack.at(i) ) continue ;
    TLorentzVector pv_i = Tk4V.at(i) ;

    int imom=0;
    int idmom=0;
    if (MC) {
      int ig = Genp_idx_clst->at(i) ;
      //cout << "ig " << ig << endl;
      if (ig>0) Hfill("Reco_id",Genp_Id->at(ig),3000,0.,3000.,1. ) ;
      if ((uint) ig<Genp_vMom->size()) {
        unsigned int mom=Genp_vMom->at(ig)[0];
        //cout << " genp " << Genp_Id->at(ig) ;
        if ( (uint) mom<Genp_Id->size()) {
          imom=mom;
          idmom=Genp_Id->at(mom);
          //cout << " mom " << Genp_Id->at(mom);
          Hfill("Reco_mom",Genp_Id->at(mom),3000,0.,3000.,1. ) ;
        }
        //cout << endl;
      }
    }

    

    // count the number of tracks in a cone of 0.1, 0.2, 0.3, 0.4, 0.5 around track i
    int ntrkR01=0;
    int ntrkR02=0;
    int ntrkR03=0;
    int ntrkR04=0;
    int ntrkR05=0;
    // count the number of tracks in a cone of 0.1, 0.2, 0.3, 0.4, 0.5 around eta = 0
    int ntrk_eta25=0;
    int ntrk_eta50=0;
    int ntrk_eta75=0;
    int ntrk_eta100=0;
    int ntrk_eta150=0;

    for( int j = 0 ; j< numberOfTrack ; j++ ) {
      if (i==j) continue;
      if( !GoodTrack.at(j) ) continue ;
      TLorentzVector pv_j = Tk4V.at(j) ;
      double deltaR=sqrt(pow(pv_i.Eta()-pv_j.Eta(),2)+
                         pow(pv_i.Phi()-pv_j.Phi(),2));
      if (deltaR<0.1) ntrkR01++;
      if (deltaR<0.2) ntrkR02++;
      if (deltaR<0.3) ntrkR03++;
      if (deltaR<0.4) ntrkR04++;
      if (deltaR<0.5) ntrkR05++;

      if (fabs(pv_j.Eta())<0.25) ntrk_eta25++;
      if (fabs(pv_j.Eta())<0.5)  ntrk_eta50++;
      if (fabs(pv_j.Eta())<0.75) ntrk_eta75++;
      if (fabs(pv_j.Eta())<1.)   ntrk_eta100++;
      if (fabs(pv_j.Eta())<1.5)  ntrk_eta150++;
    }

    for( int j = i+1 ; j< numberOfTrack ; j++ ) {
      if( !GoodTrack.at(j) ) continue ;
      TLorentzVector pv_j = Tk4V.at(j) ;

      int idmompair=0;
      if (MC) {
        int jg = Genp_idx_clst->at(j) ;
        if ((uint) jg<Genp_vMom->size()) {
          unsigned int mom=Genp_vMom->at(jg)[0];
          if ( (uint) mom<Genp_Id->size()) {
            if ( mom==imom) {
              idmompair = idmom;
              //cout << " pair mom " << idmompair;
            }
          }
        }
        Hfill("Reco_pair",idmompair,3000,0.,3000.,1. ) ;
      }

      double q = GetQ( pv_i, pv_j ) ;
      double q0 = fabs(pv_i.Energy()-pv_j.Energy()) ;
      if( q0> 0.99 ) q0= 0.99 ;

      // invert II 4 vector (control sample) 
      TLorentzVector j_vp( -pv_j.Px(), -pv_j.Py(), -pv_j.Pz() , pv_j.Energy() ) ;
      double qinv = GetQ( pv_i, j_vp ) ;
      double q0inv = fabs(pv_i.Energy()-j_vp.Energy()) ;

      // control sample, same cuts as signal
      if( fabs( pv_i.Pt()-j_vp.Pt() ) < dpt_split_cut &&
          cos( pv_i.Angle( j_vp.Vect() ) )> cst_split_cut ) qinv = 10. ;
      if( q0inv >0.99 ) q0inv=0.99  ;

      if( charge->at(i) != charge->at(j) ) 
      { 
        Hfill("qsigpm_vs_idmom",q,200,0.,2.,idmompair,3000,0.,3000.,1. ) ;
        Hfill("qinvpm_vs_idmom",qinv,200,0.,2.,idmompair,3000,0.,3000.,1. ) ;

        // plot q vs q0
        Hfill("qsigpm_vs_q0",q,200,0.,2.,q0,50,0.,1.,1. ) ;     // signal ++
        Hfill("qsigpmc_vs_q0",q,200,0.,2.,q0,50,0.,1.,CoulombW(q) ) ;     // signal ++
        Hfill("qinvpm_vs_q0",qinv,200,0.,2.,q0,50,0.,1.,1. ) ;  // control sample ++

        // plot q _vs_ eta
        Hfill("qsigpm_vs_eta",q,200,0.,2.,fabs(eta->at(i)),50,0.,2.5,1. ) ;     // signal ++
        Hfill("qinvpm_vs_eta",qinv,200,0.,2.,fabs(eta->at(i)),50,0.,2.5,1. ) ;  // control sample ++

        // plot q _vs_ pt min
        double ptmin=min(pv_i.Pt(),j_vp.Pt());
        Hfill("qsigpm_vs_ptmin",q,200,0.,2.,ptmin,100,0.,5.,1. ) ;     // signal ++
        Hfill("qinvpm_vs_ptmin",qinv,200,0.,2.,ptmin,100,0.,5.,1. ) ;  // control sample ++

        // plot q _vs_ pt avg
        double ptmean=(pv_i.Pt()+j_vp.Pt())/2.;
        Hfill("qsigpm_vs_ptmean",q,200,0.,2.,ptmean,100,0.,5.,1. ) ;     // signal ++
        Hfill("qinvpm_vs_ptmean",qinv,200,0.,2.,ptmean,100,0.,5.,1. ) ;  // control sample ++

        // plot q _vs_ n tracks 
        Hfill("qsigpm_vs_ntrk",q,200,0.,2.,numberOfTrack,100,0.,150.,1. ) ;     // signal ++
        Hfill("qinvpm_vs_ntrk",qinv,200,0.,2.,numberOfTrack,100,0.,150.,1. ) ;  // control sample ++

        // plot q _vs_ n tracks in DeltaR cones
        Hfill("qsigpm_vs_ntrk01",q,200,0.,2.,ntrkR01,25,0.,50.,1. ) ;     // signal ++
        Hfill("qinvpm_vs_ntrk01",qinv,200,0.,2.,ntrkR01,25,0.,50.,1. ) ;  // control sample ++

        Hfill("qsigpm_vs_ntrk02",q,200,0.,2.,ntrkR02,25,0.,50.,1. ) ;     // signal ++
        Hfill("qinvpm_vs_ntrk02",qinv,200,0.,2.,ntrkR02,25,0.,50.,1. ) ;  // control sample ++

        Hfill("qsigpm_vs_ntrk03",q,200,0.,2.,ntrkR03,25,0.,50.,1. ) ;     // signal ++
        Hfill("qinvpm_vs_ntrk03",qinv,200,0.,2.,ntrkR03,25,0.,50.,1. ) ;  // control sample ++

        Hfill("qsigpm_vs_ntrk04",q,200,0.,2.,ntrkR04,25,0.,50.,1. ) ;     // signal ++
        Hfill("qinvpm_vs_ntrk04",qinv,200,0.,2.,ntrkR04,25,0.,50.,1. ) ;  // control sample ++

        Hfill("qsigpm_vs_ntrk05",q,200,0.,2.,ntrkR05,25,0.,50.,1. ) ;     // signal ++
        Hfill("qinvpm_vs_ntrk05",qinv,200,0.,2.,ntrkR05,25,0.,50.,1. ) ;  // control sample ++

        // plot q _vs_ n tracks in DeltaEta cones
        Hfill("qsigpm_vs_ntrketa025",q,200,0.,2.,ntrk_eta25 ,25,0.,25.,1. ) ;     // signal ++
        Hfill("qinvpm_vs_ntrketa025",qinv,200,0.,2.,ntrk_eta25 ,25,0.,25.,1. ) ;  // control sample ++

        Hfill("qsigpm_vs_ntrketa050",q,200,0.,2.,ntrk_eta50 ,25,0.,25.,1. ) ;     // signal ++
        Hfill("qinvpm_vs_ntrketa050",qinv,200,0.,2.,ntrk_eta50 ,25,0.,25.,1. ) ;  // control sample ++

        Hfill("qsigpm_vs_ntrketa075",q,200,0.,2.,ntrk_eta75 ,25,0.,25.,1. ) ;     // signal ++
        Hfill("qinvpm_vs_ntrketa075",qinv,200,0.,2.,ntrk_eta75 ,25,0.,25.,1. ) ;  // control sample ++

        Hfill("qsigpm_vs_ntrketa100",q,200,0.,2.,ntrk_eta100 ,25,0.,25.,1. ) ;     // signal ++
        Hfill("qinvpm_vs_ntrketa100",qinv,200,0.,2.,ntrk_eta100 ,25,0.,25.,1. ) ;  // control sample ++

        Hfill("qsigpm_vs_ntrketa150",q,200,0.,2.,ntrk_eta150 ,25,0.,25.,1. ) ;     // signal ++
        Hfill("qinvpm_vs_ntrketa150",qinv,200,0.,2.,ntrk_eta150 ,25,0.,25.,1. ) ;  // control sample ++

        // plot q _vs_ n jet energy
        int jet_i=trk_jetRef->at(i);
        if (jet_i>0) {
          Hfill("qsigpm_vs_ntrkjet",q,200,0.,2.,jet_energy->at(jet_i),100,0.,50.,1. ) ;     // signal ++
          Hfill("qinvpm_vs_ntrkjet",qinv,200,0.,2.,jet_energy->at(jet_i),100,0.,50.,1. ) ;  // control sample ++
        }

        // plto q _vs_ dxy
        Hfill("qsigpm_vs_dxy",q,200,0.,2.,dxy_PV->at(i),100,0.,5.,1. ) ;     // signal ++
        Hfill("qinvpm_vs_dxy",qinv,200,0.,2.,dxy_PV->at(i),100,0.,5.,1. ) ;  // control sample ++
        continue ;
      }   // control sample +-


      if (charge->at(i)>0 && charge->at(j)>0) { //++
        // plot q vs q0
        Hfill("qsigpp_vs_q0",q,200,0.,2.,q0,50,0.,1.,1. ) ;     // signal ++
        Hfill("qsigppc_vs_q0",q,200,0.,2.,q0,50,0.,1.,CoulombW(q) ) ;     // signal ++
        Hfill("qinvpp_vs_q0",qinv,200,0.,2.,q0,50,0.,1.,1. ) ;  // control sample ++

        // plot q _vs_ eta
        Hfill("qsigpp_vs_eta",q,200,0.,2.,fabs(eta->at(i)),50,0.,2.5,1. ) ;     // signal ++
        Hfill("qinvpp_vs_eta",qinv,200,0.,2.,fabs(eta->at(i)),50,0.,2.5,1. ) ;  // control sample ++

        // plot q _vs_ pt min
        double ptmin=min(pv_i.Pt(),j_vp.Pt());
        Hfill("qsigpp_vs_ptmin",q,200,0.,2.,ptmin,100,0.,5.,1. ) ;     // signal ++
        Hfill("qinvpp_vs_ptmin",qinv,200,0.,2.,ptmin,100,0.,5.,1. ) ;  // control sample ++

        // plot q _vs_ pt avg
        double ptmean=(pv_i.Pt()+j_vp.Pt())/2.;
        Hfill("qsigpp_vs_ptmean",q,200,0.,2.,ptmean,100,0.,5.,1. ) ;     // signal ++
        Hfill("qinvpp_vs_ptmean",qinv,200,0.,2.,ptmean,100,0.,5.,1. ) ;  // control sample ++

        // plot q _vs_ n tracks 
        Hfill("qsigpp_vs_ntrk",q,200,0.,2.,numberOfTrack,100,0.,150.,1. ) ;     // signal ++
        Hfill("qinvpp_vs_ntrk",qinv,200,0.,2.,numberOfTrack,100,0.,150.,1. ) ;  // control sample ++

        // plot q _vs_ n tracks in DeltaR cones
        Hfill("qsigpp_vs_ntrk01",q,200,0.,2.,ntrkR01,25,0.,50.,1. ) ;     // signal ++
        Hfill("qinvpp_vs_ntrk01",qinv,200,0.,2.,ntrkR01,25,0.,50.,1. ) ;  // control sample ++

        Hfill("qsigpp_vs_ntrk02",q,200,0.,2.,ntrkR02,25,0.,50.,1. ) ;     // signal ++
        Hfill("qinvpp_vs_ntrk02",qinv,200,0.,2.,ntrkR02,25,0.,50.,1. ) ;  // control sample ++

        Hfill("qsigpp_vs_ntrk03",q,200,0.,2.,ntrkR03,25,0.,50.,1. ) ;     // signal ++
        Hfill("qinvpp_vs_ntrk03",qinv,200,0.,2.,ntrkR03,25,0.,50.,1. ) ;  // control sample ++

        Hfill("qsigpp_vs_ntrk04",q,200,0.,2.,ntrkR04,25,0.,50.,1. ) ;     // signal ++
        Hfill("qinvpp_vs_ntrk04",qinv,200,0.,2.,ntrkR04,25,0.,50.,1. ) ;  // control sample ++

        Hfill("qsigpp_vs_ntrk05",q,200,0.,2.,ntrkR05,25,0.,50.,1. ) ;     // signal ++
        Hfill("qinvpp_vs_ntrk05",qinv,200,0.,2.,ntrkR05,25,0.,50.,1. ) ;  // control sample ++

        // plot q _vs_ n tracks in DeltaEta cones
        Hfill("qsigpp_vs_ntrketa025",q,200,0.,2.,ntrk_eta25 ,25,0.,25.,1. ) ;     // signal ++
        Hfill("qinvpp_vs_ntrketa025",qinv,200,0.,2.,ntrk_eta25 ,25,0.,25.,1. ) ;  // control sample ++

        Hfill("qsigpp_vs_ntrketa050",q,200,0.,2.,ntrk_eta50 ,25,0.,25.,1. ) ;     // signal ++
        Hfill("qinvpp_vs_ntrketa050",qinv,200,0.,2.,ntrk_eta50 ,25,0.,25.,1. ) ;  // control sample ++

        Hfill("qsigpp_vs_ntrketa075",q,200,0.,2.,ntrk_eta75 ,25,0.,25.,1. ) ;     // signal ++
        Hfill("qinvpp_vs_ntrketa075",qinv,200,0.,2.,ntrk_eta75 ,25,0.,25.,1. ) ;  // control sample ++

        Hfill("qsigpp_vs_ntrketa100",q,200,0.,2.,ntrk_eta100 ,25,0.,25.,1. ) ;     // signal ++
        Hfill("qinvpp_vs_ntrketa100",qinv,200,0.,2.,ntrk_eta100 ,25,0.,25.,1. ) ;  // control sample ++

        Hfill("qsigpp_vs_ntrketa150",q,200,0.,2.,ntrk_eta150 ,25,0.,25.,1. ) ;     // signal ++
        Hfill("qinvpp_vs_ntrketa150",qinv,200,0.,2.,ntrk_eta150 ,25,0.,25.,1. ) ;  // control sample ++

        // plot q _vs_ n jet energy
        int jet_i=trk_jetRef->at(i);
        if (jet_i>0) {
          Hfill("qsigpp_vs_ntrkjet",q,200,0.,2.,jet_energy->at(jet_i),100,0.,50.,1. ) ;     // signal ++
          Hfill("qinvpp_vs_ntrkjet",qinv,200,0.,2.,jet_energy->at(jet_i),100,0.,50.,1. ) ;  // control sample ++
        }

        // plto q _vs_ dxy
        Hfill("qsigpp_vs_dxy",q,200,0.,2.,dxy_PV->at(i),100,0.,5.,1. ) ;     // signal ++
        Hfill("qinvpp_vs_dxy",qinv,200,0.,2.,dxy_PV->at(i),100,0.,5.,1. ) ;  // control sample ++
      }
      else { // --
        Hfill("qsigmm_vs_q0",q,200,0.,2.,q0,50,0.,1.,1. ) ;     // signal ++
        Hfill("qsigmmc_vs_q0",q,200,0.,2.,q0,50,0.,1.,CoulombW(q) ) ;     // signal ++
        Hfill("qinvmm_vs_q0",qinv,200,0.,2.,q0,50,0.,1.,1. ) ;  // control sample ++

        // plot q _vs_ eta
        Hfill("qsigmm_vs_eta",q,200,0.,2.,fabs(eta->at(i)),50,0.,2.5,1. ) ;     // signal ++
        Hfill("qinvmm_vs_eta",qinv,200,0.,2.,fabs(eta->at(i)),50,0.,2.5,1. ) ;  // control sample ++

        // plot q _vs_ pt min
        double ptmin=min(pv_i.Pt(),j_vp.Pt());
        Hfill("qsigmm_vs_ptmin",q,200,0.,2.,ptmin,100,0.,5.,1. ) ;     // signal ++
        Hfill("qinvmm_vs_ptmin",qinv,200,0.,2.,ptmin,100,0.,5.,1. ) ;  // control sample ++

        // plot q _vs_ pt avg
        double ptmean=(pv_i.Pt()+j_vp.Pt())/2.;
        Hfill("qsigmm_vs_ptmean",q,200,0.,2.,ptmean,100,0.,5.,1. ) ;     // signal ++
        Hfill("qinvmm_vs_ptmean",qinv,200,0.,2.,ptmean,100,0.,5.,1. ) ;  // control sample ++

        // plot q _vs_ n tracks 
        Hfill("qsigmm_vs_ntrk",q,200,0.,2.,numberOfTrack,100,0.,150.,1. ) ;     // signal ++
        Hfill("qinvmm_vs_ntrk",qinv,200,0.,2.,numberOfTrack,100,0.,150.,1. ) ;  // control sample ++

        // plot q _vs_ n tracks in DeltaR cones
        Hfill("qsigmm_vs_ntrk01",q,200,0.,2.,ntrkR01,25,0.,50.,1. ) ;     // signal ++
        Hfill("qinvmm_vs_ntrk01",qinv,200,0.,2.,ntrkR01,25,0.,50.,1. ) ;  // control sample ++

        Hfill("qsigmm_vs_ntrk02",q,200,0.,2.,ntrkR02,25,0.,50.,1. ) ;     // signal ++
        Hfill("qinvmm_vs_ntrk02",qinv,200,0.,2.,ntrkR02,25,0.,50.,1. ) ;  // control sample ++

        Hfill("qsigmm_vs_ntrk03",q,200,0.,2.,ntrkR03,25,0.,50.,1. ) ;     // signal ++
        Hfill("qinvmm_vs_ntrk03",qinv,200,0.,2.,ntrkR03,25,0.,50.,1. ) ;  // control sample ++

        Hfill("qsigmm_vs_ntrk04",q,200,0.,2.,ntrkR04,25,0.,50.,1. ) ;     // signal ++
        Hfill("qinvmm_vs_ntrk04",qinv,200,0.,2.,ntrkR04,25,0.,50.,1. ) ;  // control sample ++

        Hfill("qsigmm_vs_ntrk05",q,200,0.,2.,ntrkR05,25,0.,50.,1. ) ;     // signal ++
        Hfill("qinvmm_vs_ntrk05",qinv,200,0.,2.,ntrkR05,25,0.,50.,1. ) ;  // control sample ++

        // plot q _vs_ n tracks in DeltaEta cones
        Hfill("qsigmm_vs_ntrketa025",q,200,0.,2.,ntrk_eta25 ,25,0.,25.,1. ) ;     // signal ++
        Hfill("qinvmm_vs_ntrketa025",qinv,200,0.,2.,ntrk_eta25 ,25,0.,25.,1. ) ;  // control sample ++

        Hfill("qsigmm_vs_ntrketa050",q,200,0.,2.,ntrk_eta50 ,25,0.,25.,1. ) ;     // signal ++
        Hfill("qinvmm_vs_ntrketa050",qinv,200,0.,2.,ntrk_eta50 ,25,0.,25.,1. ) ;  // control sample ++

        Hfill("qsigmm_vs_ntrketa075",q,200,0.,2.,ntrk_eta75 ,25,0.,25.,1. ) ;     // signal ++
        Hfill("qinvmm_vs_ntrketa075",qinv,200,0.,2.,ntrk_eta75 ,25,0.,25.,1. ) ;  // control sample ++

        Hfill("qsigmm_vs_ntrketa100",q,200,0.,2.,ntrk_eta100 ,25,0.,25.,1. ) ;     // signal ++
        Hfill("qinvmm_vs_ntrketa100",qinv,200,0.,2.,ntrk_eta100 ,25,0.,25.,1. ) ;  // control sample ++

        Hfill("qsigmm_vs_ntrketa150",q,200,0.,2.,ntrk_eta150 ,25,0.,25.,1. ) ;     // signal ++
        Hfill("qinvmm_vs_ntrketa150",qinv,200,0.,2.,ntrk_eta150 ,25,0.,25.,1. ) ;  // control sample ++

        // plot q _vs_ n jet energy
        int jet_i=trk_jetRef->at(i);
        if (jet_i>0) {
          Hfill("qsigmm_vs_ntrkjet",q,200,0.,2.,jet_energy->at(jet_i),100,0.,50.,1. ) ;     // signal ++
          Hfill("qinvmm_vs_ntrkjet",qinv,200,0.,2.,jet_energy->at(jet_i),100,0.,50.,1. ) ;  // control sample ++
        }

        // plto q _vs_ dxy
        Hfill("qsigmm_vs_dxy",q,200,0.,2.,dxy_PV->at(i),100,0.,5.,1. ) ;     // signal ++
        Hfill("qinvmm_vs_dxy",qinv,200,0.,2.,dxy_PV->at(i),100,0.,5.,1. ) ;  // control sample ++
      }
    }
  }
  return ;
}

// =======================================================================
// Compute correlated variable Q (= sqrt(qsq))
// =======================================================================
double EventMixer::GetQ( const TLorentzVector &p1, const TLorentzVector &p2 ) {
  TLorentzVector Sum4V = p1+p2 ;
  double q = Sum4V.Mag2() - 4*pimass*pimass  ;
  return (  q>0 ?  sqrt(q) : -sqrt(-q)  ) ; 
}

// =======================================================================
// return the wight factor due to Coloumb repulsion [Gamow] same charge
// =======================================================================
double EventMixer::CoulombW( const double& q) const {
  const double alpha=1./137.;
  double x=2*TMath::Pi()*(alpha*pimass/q);
  return (exp(x)-1.)/x;
}
// =======================================================================
// return the wight factor due to Coloumb attraction [Gamow] opposite charge
// =======================================================================
double EventMixer::CoulombWpm( const double& q) const {
  const double alpha=1./137.;
  double x=2*TMath::Pi()*(alpha*pimass/q);
  return (1.-exp(-x))/x;
}
// =======================================================================
// Select good track, all cuts but ndof
// =======================================================================
inline void EventMixer::SelectGoodTracks() {
  Tk4V.erase( Tk4V.begin(), Tk4V.end() ) ;
  GoodTrack.erase( GoodTrack.begin() , GoodTrack.end() ) ;

  TLorentzVector pTk ;
  for ( int i = 0; i < numberOfTrack ; i++ ) {
    pTk.SetXYZM( px->at(i), py->at(i), pz->at(i), pimass ) ;
    Tk4V.push_back( pTk ) ;

    GoodTrack.push_back(false) ;          
    if( !GoodTrackNoDofCut(i) ) continue ; // apply all cuts but that on NDof
    GoodTrack[i] = true ;

    if( charge->at(i) > 0 ) { 
      Hfill( "EtaVsPtGoodP", pt->at(i), 100,0.,2.,eta->at(i), 60,-3.,3.,1. ) ; 
      Hfill( "EtaVsNdofGoodP", (double) ndof->at(i), 50,0.,50.,eta->at(i), 60,-3.,3.,1. ) ; 
      if( ndof->at(i) > ndof_cut ) Hfill( "EtaVsPtGoodP_ndof", pt->at(i), 100,0.,2.,eta->at(i), 60,-3.,3.,1. ) ;
      // if( fabs(eta->at(i)) < 1.5 && pt->at(i)>0.25 && pt->at(i)<0.45 ) 
      // { Hfill( "DdxVsPtP",p->at(i),50,0.25,1.25,newdEdx1->at(i),100,0.,15.,1.) ;}
    }
    else { 
      Hfill( "EtaVsPtGoodM", pt->at(i), 100,0.,2.,eta->at(i), 60,-3.,3.,1. ) ;
      Hfill( "EtaVsNdofGoodM", (double) ndof->at(i), 50,0.,50.,eta->at(i), 60,-3.,3.,1. ) ; 
      if( ndof->at(i) > ndof_cut ) Hfill( "EtaVsPtGoodM_ndof", pt->at(i), 100,0.,2.,eta->at(i), 60,-3.,3.,1. ) ; 
      // if( fabs(eta->at(i)) < 1.5 && pt->at(i)>0.25 && pt->at(i)<0.45 ) 
      // { Hfill( "DdxVsPtM",p->at(i),50,0.25,1.25,newdEdx1->at(i),100,0.,15.,1.) ;}
    }
  } 
  return ;
}

inline void  EventMixer::FillGoodTracks() 
{  
  TkGood.erase( TkGood.begin(), TkGood.end() ) ;
  ChGood.erase( ChGood.begin(), ChGood.end() ) ;
  for( uint i = 0 ; i<Tk4V.size() ; i++ )
    {
      if( !GoodTrack.at(i) ) continue ;
      ChGood.push_back( charge->at(i) ) ;
      TkGood.push_back( Tk4V.at(i) ) ;
    }
  return ;
}
// =======================================================================
// remove tracks having cos(theta)>cst_split_cut and 
// abs(delta p) < dpt_splt_cut  with another equal charge track 
// fill also some control histo
// =======================================================================
inline void  EventMixer::RemoveSplitTracks() {

  SplitTrack.erase( SplitTrack.begin(), SplitTrack.end() ) ;
  for( int i = 0 ; i<numberOfTrack ; i++ ) { SplitTrack.push_back( false ) ;}
  
  for( int i = 0 ; i<numberOfTrack ; i++ ) {
    if( !GoodTrack.at(i) ) continue ;

    TLorentzVector pv_i = Tk4V.at(i) ;
    
    for( int j = i+1 ; j<numberOfTrack ; j++ ) {
      if( !GoodTrack.at(j) ) continue ;

      TLorentzVector pv_j = Tk4V.at(j) ;
      double dpt = pv_i.Pt()-pv_j.Pt() ;
      double cst = cos( pv_i.Angle( pv_j.Vect() ) ) ;	 
      double q = GetQ( pv_i, pv_j ) ;

      if( charge->at(i) != charge->at(j) ) {
        Hfill("DptVsCst_pm",cst,50,0.9991,1.0001, dpt,50,-.5,.5,1.) ;
        continue ;
      }

      // cos theta vs delta pt (/pt) for all charged tracks
      string htit = "split" ;
      if( charge->at(i) < 0 ) { htit.append("_m") ; }
      Hfill((char*) htit.c_str(),cst,50,0.9991,1.0001, dpt,50,-.5,.5,1.) ;
      htit.append("R") ;
      Hfill((char*) htit.c_str(),cst,50,0.9991,1.0001, 2*dpt/(pv_i.Pt()+pv_j.Pt()),50,-1.,1.,1.) ;

      if( cst >cst_split_cut && fabs(dpt)< dpt_split_cut ) { 
        // remove both tracks!
        GoodTrack[j] = false ; 
        SplitTrack[j] = true ;
        GoodTrack[i] = false ; 
        SplitTrack[i] = true ;
        Hfill("q_split",q,200,0.,.2,1. ) ;
        Hfill("ndof_split",ndof->at(j),50,0.,50.,1.) ;
        Hfill("eta_vs_pt_split",pt->at(j),50,0.,1.,eta->at(j),60,-3.,3.,1.) ;
        Hfill("eta_vs_phi_split",phi->at(j),65,-3.25,3.25,eta->at(j),60,-3.,3.,1.) ;
        Hfill("eta_vs_dz_split",dz->at(j),50,-5.,5.,eta->at(j),60,-3.,3.,1.) ;
        Hfill("eta_vs_ddz_split",dz->at(j)-dz->at(i),100,-1.,1.,eta->at(j),60,-3.,3.,1.) ;
        Hfill("eta_vs_ddxy_split",dxy_PV->at(j)-dxy_PV->at(i),100,-1.,1.,eta->at(j),60,-3.,3.,1.) ;
        Hfill("NdfjVsNdfi",(double) ndof->at(j), 50,0.,50., ndof->at(i), 50, 0., 50., 1.) ;
        Hfill("NdfDiff",(double) ndof->at(j)-ndof->at(i), 41, -20.5, 20.5, 1.) ;
      }
    }
  }
  return ;
}

// =======================================================================
inline void EventMixer::SaveToFile( char* OutputFileName ) { SaveHist( OutputFileName ) ;}

// =======================================================================
inline void EventMixer::end_analysis() {}

// =======================================================================
inline void EventMixer::begin_analysis() { return ; }

// =======================================================================
inline int EventMixer::GoodTrackNoDofCut( int i ) {
  if( RejectTrack(i) == 0 ) return 1 ;
  //if( RejectTrack(i) == 0 || RejectTrack(i) == (int) pow(2,4) ) return 1 ;
  return 0 ;
}

// =======================================================================
inline int EventMixer::RejectTrack( int i ) {
  int ReturnValue(0) ;

  if( !highPurity_track->at(i) )         ReturnValue += (int) pow(2,0) ;
  if( pt->at(i) < pt_cut )               ReturnValue += (int) pow(2,1) ;
  double inner_r=sqrt(innerPosition_x->at(i)*innerPosition_x->at(i)+
                      innerPosition_y->at(i)*innerPosition_y->at(i));
  if( fabs( dxy_PV->at(i)) > dxy_cut || inner_r>rmin_cut ) 
                                         ReturnValue += (int) pow(2,2) ;
  if( normalizedChi2->at(i) > nchi2_cut) ReturnValue += (int) pow(2,3) ;
  if( ndof->at(i) < ndof_cut )           ReturnValue += (int) pow(2,4) ;
  if( fabs( eta->at(i) ) > eta_cut )     ReturnValue += (int) pow(2,5) ;

  FillCutControlHisto( i, ReturnValue ) ;

  return ReturnValue ;
}

// =======================================================================
inline void EventMixer::FillCutControlHisto( int itrack, int rcut ) {
  char htitle[100] ;
  for( int jcut=0 ; jcut<nTrackCut ; jcut++ ) { 
    if( (rcut&(int) pow(2,jcut) ) ) continue ;
    sprintf( htitle, "etacut%d", jcut ) ;
    if( charge->at(itrack) > 0 ) 
    { strcat( htitle,"_plus" ) ; }
    else
    { strcat( htitle,"_minus" ) ; }
    Hfill( htitle, eta->at(itrack), 60,-3.,3.,1. ) ;      
  }
}