#define macro_PhiAnalysis_cxx
#include "macro_PhiAnalysis.h"
#include <TH2.h>
#include <TH1.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TLorentzVector.h>
#include <TMath.h>
#include <vector>
#include <math.h>
#include "myHist.h"

using namespace std;
TH2D *looper_m_low ,*looper_m_high ,*looper_p_low ,*looper_p_high ;

static const Int_t ndof_cut(5) ;
static const Int_t nchi2_cut(2) ;
static const Double_t dxy_cut(1.) ;
static const Double_t pt_cut(0.1) ;
static const Double_t eta_cut(2.5) ;

static const Double_t pimass(0.1396), kamass(0.4937) ; 

static const Double_t leftinf(0.99178) ; 
static const Double_t leftsup=1.00369;
static const Double_t rightinf=1.03545;
static const Double_t rightsup=1.04736;
static const Double_t centerinf=1.00766;
static const Double_t centersup=1.03148;

static const Double_t cst_looper_cut(0.9995) ;
static const Double_t dpt_looper_cut(0.04 ) ;

inline bool macro_PhiAnalysis::Looper(Int_t i, Int_t j, TLorentzVector *pi, TLorentzVector *pj) 
{
  if( charge->at(i)*charge->at(j) < 0 ) return false ;
  Double_t dpt = pi->Pt()-pj->Pt() ;
  Double_t cst = cos( pi->Angle( pj->Vect() ) ) ;
  if( cst > cst_looper_cut && fabs(dpt) < dpt_looper_cut ) return true ;
  return false ;
}
inline void macro_PhiAnalysis::FillSequentialHisto( Int_t i ) 
{
}

inline bool macro_PhiAnalysis::SideBand( Double_t mass ) 
{
  if(( mass>leftinf && mass<leftsup ) || 
     ( mass>rightinf && mass<rightsup) ) return true ;
  return false ;
}
inline bool macro_PhiAnalysis::SignalRegion( Double_t mass ) 
{ return ( mass>centerinf && mass<centersup ); }

inline Int_t macro_PhiAnalysis::RejectTrack( Int_t i ) 
{
  if( !highPurity_track->at(i) ) return 7 ;
  if( pt->at(i) < pt_cut ) return 2 ; 
  if( fabs( dxy_PV->at(i)) > dxy_cut ) return 3 ;
  if( normalizedChi2->at(i) > nchi2_cut) return 4 ;
  if( ndof->at(i) < ndof_cut ) return 5;
  if( fabs( eta->at(i) ) > eta_cut ) return 6 ;
  return 10 ;
}
inline bool macro_PhiAnalysis::GoodTrack( Int_t i ) 
{  
  return ( RejectTrack(i) == 10 ) ;
}
void macro_PhiAnalysis::Loop()
{
   if (fChain == 0) return;
 
   //VARIABLE DECLARATIONS
   Int_t mass_bins = 4000;
   Double_t mass_min = 0;
   Double_t mass_max = 4;

   //HISTO DECLARATIONS
   TH1D *phimass_2_05= new TH1D("phimass_2_05","mass of phi candidates chi2 2 pt 0.5", mass_bins, mass_min, mass_max);
   TH1D *phimass_2_05_nodedx= new TH1D("phimass_2_05_nodedx","mass of phi candidates chi2 2 pt 0.5", mass_bins, mass_min, mass_max);

   TH2D *seldedxvsp= new TH2D("seldedxvsp","sel dedx 1 vs logPt", 200, -1, 2, 200, 0, 10);
   
   TH1D * cands = new TH1D ("cands","Phi->kk candidates", 200, 0.98, 1.18);

   TH1D * Eta_mw = new TH1D("Eta_mw","Eta Phi cands mw", 10, -2.5, 2.5);
   TH1D * Pt_mw  = new TH1D("Pt_mw", "Pt Phi cands mw", 10, 0., 10.);
   TH1D * Eta_sb = new TH1D("Eta_sb","Eta Phi cands sb", 10, -2.5, 2.5);
   TH1D * Pt_sb  = new TH1D("Pt_sb", "Pt Phi cands sb", 10, 0., 10.);

   looper_p_low = new TH2D("looper_p_low","cos(theta) vs drho +",100,-.1,.1,100,0.9991,1.0001);
   looper_p_high= new TH2D("looper_p_high","cos(theta) vs drho +",100,-.1,.1,100,0.9991,1.0001);
   looper_m_low = new TH2D("looper_m_low","cos(theta) vs drho -",100,-.1,.1,100,0.9991,1.0001);   
   looper_m_high = new TH2D("looper_m_high","cos(theta) vs drho -",100,-.1,.1,100,0.9991,1.0001);

   TH2D * hvsM = new TH2D("hvsM","hEly vs Mass",100, 0.98, 1.08,50,-1.,1.);
   TH1D * h_sb = new TH1D("h_sb","Elycity Side Band", 100,-1.,1.);
   TH1D * h_mw = new TH1D("h_mw","Elycity Signal Region", 100,-1.,1.);

   // BE
   TH1D* q_pp = new TH1D("q_pp","qsquare ++",100,0.,2.) ;
   TH1D* q_pm = new TH1D("q_pm","qsquare +-",100,0.,2.) ;
   TH1D* q_mp = new TH1D("q_mp","qsquare -+",100,0.,2.) ;
   TH1D* q_mm = new TH1D("q_mm","qsquare --",100,0.,2.) ;
   TH1D* looper_q = new TH1D("looper_q","looper qsquare",100,0.,2.) ;

   // debug

   TH1D*    eta_ptp = new TH1D("eta_ptp","eta +",200,-3.,3.) ;
   TH1D*    eta_dxyp = new TH1D("eta_dxyp","eta +",200,-3.,3.) ;
   TH1D*    eta_x2np = new TH1D("eta_x2np","eta +",200,-3.,3.) ;
   TH1D*    eta_HPTp = new TH1D("eta_HPTp","eta +",200,-3.,3.) ;
   TH1D*    eta_ndofp = new TH1D("eta_ndofp","eta +",200,-3.,3.) ;

   TH1D*    eta_ptm = new TH1D("eta_ptm","eta +",200,-3.,3.) ;
   TH1D*    eta_dxym = new TH1D("eta_dxym","eta +",200,-3.,3.) ;
   TH1D*    eta_x2nm = new TH1D("eta_x2nm","eta +",200,-3.,3.) ;
   TH1D*    eta_HPTm = new TH1D("eta_HPTm","eta +",200,-3.,3.) ;
   TH1D*    eta_ndofm = new TH1D("eta_ndofm","eta +",200,-3.,3.) ;
 
   TH1D* ptp = new TH1D("ptp","pt +",200,0.,5.) ;
   TH1D* ptm = new TH1D("ptm","pt -",200,0.,5.) ;
   TH1D* ptp_sel = new TH1D("ptp_sel","pt + sel",200,0.,5.) ;
   TH1D* ptm_sel = new TH1D("ptm_sel","pt - sel",200,0.,5.) ;

   TH1D* etap = new TH1D("etap","eta +",200,-3.,3.) ;
   TH1D* etam = new TH1D("etam","eta -",200,-3.,3.) ;
   TH1D* etap_sel = new TH1D("etap_sel","eta + sel",200,-3.,3.) ;
   TH1D* etam_sel = new TH1D("etam_sel","eta - sel",200,-3.,3.) ;

   TH1D* phip = new TH1D("phip","phi +",126,-3.15,3.15) ;
   TH1D* phim = new TH1D("phim","phi -",126,-3.15,3.15) ;
   TH1D* phip_sel = new TH1D("phip_sel","phi + sel",126,-3.15,3.15) ;
   TH1D* phim_sel = new TH1D("phim_sel","phi - sel",126,-3.15,3.15) ;

   TH1D* dxyp = new TH1D("dxyp","dxy +",100,-1.5,1.5) ;
   TH1D* dxym = new TH1D("dxym","dxy -",100,-1.5,1.5) ;
   TH1D* dxyp_sel = new TH1D("dxyp_sel","dxy + sel",100,-1.5,1.5) ;
   TH1D* dxym_sel = new TH1D("dxym_sel","dxy - sel",100,-1.5,1.5) ;

   TH1D* dzp = new TH1D("dzp","dz +",100,-10.,10.) ;
   TH1D* dzm = new TH1D("dzm","dz -",100,-10.,10.) ;
   TH1D* dzp_sel = new TH1D("dzp_sel","dz + sel",100,-10.,10.) ;
   TH1D* dzm_sel = new TH1D("dzm_sel","dz - sel",100,-10.,10.) ;

   TH1D* x2np = new TH1D("x2np","x2n +",100,0.,5.) ;
   TH1D* x2nm = new TH1D("x2nm","x2n -",100,0.,5.) ;
   TH1D* x2np_sel = new TH1D("x2np_sel","x2n + sel",100,0.,5.) ;
   TH1D* x2nm_sel = new TH1D("x2nm_sel","x2n - sel",100,0.,5.) ;

   TH1D* ndofp = new TH1D("ndofp","ndof +",70,0.,70.) ;
   TH1D* ndofm = new TH1D("ndofm","ndof -",70,0.,70.) ;
   TH1D* ndofp_sel = new TH1D("ndofp_sel","ndof + sel",70,0.,70.) ;
   TH1D* ndofm_sel = new TH1D("ndofm_sel","ndof - sel",70,0.,70.) ;

   TH1D* looper_dz = new TH1D("looper_dz","dz looper",100,-0.1,0.1) ;
   TH1D* no_looper_dz = new TH1D("no_looper_dz","dz looper",100,-0.1,0.1) ;
 
   cout << " End Booking " << endl ;

   Int_t cout_n = 100;
   Bool_t debugg = true;
   
   //Long64_t nentries = fChain->GetEntriesFast();
   Long64_t nentries = fChain->GetEntries();
   if(debugg){
      nentries = 100;
      cout << "debug mode: nentries= " << nentries << endl;
   }
   std::cout << "total entries: " << nentries << std::endl;
   std::cout << "std::cout every " << cout_n << " lines" << std::endl;
   std::cout << "==========================" << std::endl;
   std::cout << "" << std::endl;

   Long64_t nbytes = 0, nb = 0;
   
   for (Long64_t jentry=0; jentry<nentries;jentry++) { // event loop

      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      cout << "cfr " << ientry << " " << jentry << endl ;
      nb = fChain->GetEntry(jentry);   nbytes += nb;

      if( (jentry%cout_n)==0 ) 	std::cout << "processing entry " << jentry+1 
					  << " ; Run " << runNumber << " event " << evNumber << std::endl;
      
      if( debugg ) std::cout << "l'evento ha " << numberOfTrack <<" tracce" << std::endl;
      if( numberOfTrack>150 || numberOfTrack<2 ) continue;
      
      TLorentzVector p_pi_i, p_pi_j, p_ka_i, p_ka_j ;
      TLorentzVector p_pi_sum, p_ka_sum ;

      // first loop on charged tracks
      for( int i = 0; i < numberOfTrack; i++ )
	{
	  Int_t k = RejectTrack(i) ;
	  if( charge->at(i) > 0 ) 
	    { 
	      if( k>1 ) eta_HPTp->Fill( eta->at(i) ) ;
	      if( k>3 ) eta_dxyp->Fill( eta->at(i) ) ;
	      if( k>2 ) eta_ptp->Fill( eta->at(i) ) ; 
	      if( k>5 ) eta_ndofp->Fill( eta->at(i) ) ;
	      if( k>4 ) eta_x2np->Fill( eta->at(i) ) ;
	    }
	  else
	    { 
	      if( k>1 ) eta_HPTm->Fill( eta->at(i) ) ;
	      if( k>3 ) eta_dxym->Fill( eta->at(i) ) ;
	      if( k>2 ) eta_ptm->Fill( eta->at(i) ) ; 
	      if( k>5 ) eta_ndofm->Fill( eta->at(i) ) ;
	      if( k>4 ) eta_x2nm->Fill( eta->at(i) ) ;
	    }


	  if( charge->at(i) > 0 ) 
	    {
	      ptp->Fill( pt->at(i) ) ;
	      etap->Fill( eta->at(i) ) ;
	      phip->Fill( phi->at(i) ) ;
	      dxyp->Fill( dxy_PV->at(i) ) ;
	      dzp->Fill( dz->at(i) ) ;
	      x2np->Fill( normalizedChi2->at(i) ) ;
	      ndofp->Fill( ndof->at(i) ) ;
	    }
	  else
	    {
	      ptm->Fill( pt->at(i) ) ;
	      etam->Fill( eta->at(i) ) ;
	      phim->Fill( phi->at(i) ) ;
	      dxym->Fill( dxy_PV->at(i) ) ;
	      dzm->Fill( dz->at(i) ) ;
	      ndofm->Fill( ndof->at(i) ) ;
	      x2nm->Fill( normalizedChi2->at(i) ) ;
	    }
	  if( !GoodTrack(i) ) continue ;

	  if( charge->at(i) > 0 ) 
	    {
	      ptp_sel->Fill( pt->at(i) ) ;
	      etap_sel->Fill( eta->at(i) ) ;
	      phip_sel->Fill( phi->at(i) ) ;
	      dxyp_sel->Fill( dxy_PV->at(i) ) ;
	      dzp_sel->Fill( dz->at(i) ) ;
	      ndofp_sel->Fill( ndof->at(i) ) ;
	      x2np_sel->Fill( normalizedChi2->at(i) ) ;
	    }
	  else
	    {
	      ptm_sel->Fill( pt->at(i) ) ;
	      etam_sel->Fill( eta->at(i) ) ;
	      phim_sel->Fill( phi->at(i) ) ;
	      dxym_sel->Fill( dxy_PV->at(i) ) ;
	      dzm_sel->Fill( dz->at(i) ) ;
	      x2nm_sel->Fill( normalizedChi2->at(i) ) ;
	      ndofm_sel->Fill( ndof->at(i) ) ;
	    }

	  Bool_t passed_i = true;
	  if( p->at(i)<1 ) passed_i = K_comp_dEdx_->at(i);
	  
	  //	  if(debugg) std::cout << "dopo caricamento variabili 1" << std::endl;
	  // if(debugg) std::cout << "\t nel primo loop" << std::endl;
	  
	  p_pi_i.SetXYZM( px->at(i), py->at(i), pz->at(i), pimass ) ;
	  p_ka_i.SetXYZM( px->at(i), py->at(i), pz->at(i), kamass ) ;

	  for( int j = i+1; j < numberOfTrack; j++)
	    {	     
	      if( !GoodTrack(j) ) continue ;
	      // if(debugg) std::cout << "\t dopo caricamento variabili 2" << std::endl;

	      Bool_t passed_j = true;
	      if( p->at(j)<1 ) passed_j = K_comp_dEdx_->at(j);	  
	      
	      p_pi_j.SetXYZM( px->at(j), py->at(j), pz->at(j), pimass ) ;
	      p_ka_j.SetXYZM( px->at(j), py->at(j), pz->at(j), kamass ) ;
	      p_pi_sum = p_pi_i + p_pi_j ;
	      p_ka_sum = p_ka_i + p_ka_j ;
	      
	      Double_t Q = p_pi_sum.Mag2() - 4*pimass*pimass ;
	      Q = ( Q>=0 ? sqrt(Q) : -sqrt(-Q) ) ;

	      // remove loopers
	      if( charge->at(i) == charge->at(j) && Looper(i,j, &p_pi_i, &p_pi_j ) )
		{ 
		  looper_q->Fill( Q ) ;
		  looper_dz->Fill( dz->at(i)-dz->at(j) ) ;
		  continue ;
		}
	      no_looper_dz->Fill( dz->at(i)-dz->at(j) ) ;
	      if( charge->at(i)>0 && charge->at(j) >0 ) q_pp->Fill(Q) ; 
	      if( charge->at(i)<0 && charge->at(j) >0 ) q_mp->Fill(Q) ; 
	      if( charge->at(i)>0 && charge->at(j) <0 ) q_pm->Fill(Q) ; 
	      if( charge->at(i)<0 && charge->at(j) <0 ) q_mm->Fill(Q) ; 

	      Double_t CosHely(999) ;
	      if( charge->at(i)>0 ) 
		{ CosHely = GetHelicity( &p_ka_i, &p_ka_j) ;}
	      else
		{ CosHely = GetHelicity( &p_ka_j, &p_ka_i) ;}
	      
	      // if( debugg ) std::cout << "\t\t prima del filling!" << std::endl;

	      if( charge->at(i) == charge->at(j) ) continue ;

	      Double_t mass = p_ka_sum.Mag() ;
	      if( mass >1.01 && mass<1.03 ) 
		{
		  seldedxvsp->Fill(log10(p->at(i)),dEdx1->at(i));
		  seldedxvsp->Fill(log10(p->at(j)),dEdx1->at(j));
		}
	      if( passed_i && passed_j ) 
		{
		  phimass_2_05->Fill(mass);
		  hvsM->Fill(mass,CosHely) ;
		  cands->Fill(mass);
		}
	      else
		{
		  phimass_2_05_nodedx->Fill(mass) ;
		}
	      if( SideBand( mass ) ) 
		{ 
		  Eta_sb->Fill(p_ka_sum.Rapidity());
		  Pt_sb->Fill( p_ka_sum.Pt());
		  h_sb->Fill(CosHely) ;
		}
	      else if( SignalRegion( mass ) ) 
		{
		  Eta_mw->Fill(p_ka_sum.Rapidity());
		  Pt_mw->Fill( p_ka_sum.Pt());
		  h_mw->Fill(CosHely) ;
		}
	     
	      //   if(debugg) std::cout << "\t\t dopo il filling!" << std::endl;
	    }//second track (j) loop end	  
	}//first track loop end
      
     if(debugg) std::cout << "fine loop evento" << std::endl;
     if(debugg) std::cout << "====================" << std::endl;
          
   } // event loop end
   
   outputRootFile->Write();
   
}// Loop method end