#define CPmix_cxx
#include "CPmix.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TLorentzVector.h>
#include <TMath.h>
#include <iostream>
#include <vector>
#include "HSum.h"
#include <stdio.h> 
#include "mypart.h"

using namespace std ;

int nTrue(0),nTrueSel(0),nTrueBest(0) ;

int iMaxTag(8), iMaxDcy(4) ;
double dMaxTag = (double) iMaxTag ;
double dMaxDcy = (double) iMaxDcy ;

mypart theLep, thePst ;
vector <mypart> TrackContainer ;
vector <int> TrackSource ;

vector <TH1F*> ListTH1 ;
vector <TH2F*> ListTH2 ;

// =============================================================================
inline int CPmix::Hfill( char* hName , 
			 Double_t xvar, Int_t nx, Double_t xmin, Double_t xmax, 
			 Double_t yvar, Int_t ny, Double_t ymin, Double_t ymax, 
			 Double_t w ) 
{
  Int_t retval(1) ;
  TH2F* h = (TH2F*) gDirectory->Get( hName ) ;
  if( yvar < ymin ) yvar = ymin*1.001 ;
  if( yvar > ymax ) yvar = ymax*0.999 ;
  if( h == 0 ) { 
    cout << "create non existing histogram " << hName << endl ;
    TH2F *h = new TH2F( hName , hName, nx, xmin, xmax, ny, ymin,ymax ) ; 
    h->Sumw2() ;
    h->Fill( xvar, yvar, w ) ;
    ListTH2.push_back( h ) ;
    return 0 ;
  }
  h->Fill( xvar, yvar, w ) ;
  return retval ;
} 
// =============================================================================
inline int CPmix::Hfill( char* hName , Double_t var, 
			 Int_t n, Double_t xmin, Double_t xmax, Double_t w ) 
{
  Int_t retval(1) ;
  TH1F* h = (TH1F*) gDirectory->Get( hName ) ;
  if( var < xmin ) var = xmin*1.001 ;
  if( var > xmax ) var = xmax*0.999 ;
  if( h == 0 ) { 
    cout << "create non existing histogram " << hName << endl ;
    TH1F *h = new TH1F( hName , hName, n, xmin, xmax ) ; 
    ListTH1.push_back( h ) ;
    h->Sumw2() ;
    h->Fill( var, w ) ;
    return 0 ;
  }
  h->Fill( var, w ) ;
  return retval ;
} 
// =============================================================================
inline int CPmix::Hfill( char* hNam , int ijk, Double_t var, 
			 Int_t n, Double_t xmin, Double_t xmax, Double_t w ) 
{
  Int_t retval(1) ;

  char hName[100] ;
  int k = sprintf( hName,"%s_%d",hNam,ijk) ;
  
  TH1F* h = (TH1F*) gDirectory->Get( hName ) ;
  if( var < xmin ) var = xmin*1.001 ;
  if( var > xmax ) var = xmax*0.999 ;
  if( h == 0 ) { 
    cout << "create non existing histogram " << hName << " " << ijk << endl ;
    TH1F *h = new TH1F( hName , hName, n, xmin, xmax ) ; 
    ListTH1.push_back( h ) ;
    h->Sumw2() ;
    h->Fill( var, w ) ;
    return 0 ;
  }
  h->Fill( var, w ) ;
  return retval ;
} 
// =============================================================================
void CPmix::Loop( int nmax, int step, char* OutputFileName )
{
   if (fChain == 0) return;

   Long64_t nentries = fChain->GetEntriesFast();
   if( nentries > nmax ) nentries = nmax ;

   int processed (0) ;
   for (Long64_t jentry=0; jentry<nentries;jentry+=step)
     {
       if( jentry >= nentries ) break ; 
       Long64_t ientry = LoadTree(jentry);
       if (ientry < 0) break;
       Long64_t nbytes = 0, nb = 0;
       nb = fChain->GetEntry(jentry);   nbytes += nb;
       if( (processed++ )%1000 == 0 ) 
	 { cout << processed << " " << jentry << " " << ientry << endl ;}

       Hfill("Cut",(Double_t)Cut(ientry),10,-8.,2.,1);
       if (Cut(ientry) < 0) continue;
       analyze() ;
     } 
   cout << " end of Loop ; processed " << processed << " events " << endl ;
   EndJob(OutputFileName ) ;

}
// =============================================================================
int CPmix::GetTrackSource( int i ) {  return abs(TkTag[i]) / 1000 ; }
// =============================================================================
int CPmix::GetCharge( int i )
{
  if( TkTag[i]>0 ) return 1;
  return -1 ;
}
// =============================================================================
double CPmix::GetMass( int i )
{ 
 double mass(0.1395) ;
  int tkpid = abs(TkTag[i])%1000 ;

  if( tkpid == 666 )     { return mass ; }
  else if( tkpid>=100 )  { mass = 0.495 ;}
  else if( tkpid >=10  ) { mass = 0.0005;  }
  else if( tkpid == 1 )  { mass = 0.105 ; }

  return mass ;
}
// =============================================================================
bool CPmix::GoodTrack( int itk ) 
{
  if( abs(TkTag[itk]) == 666 ) return false ;   // skip soft pion
  if( TkeVt[itk] <= 0. ) return false ;         // vertex quality
  if( TkeVt[itk] > 0.5 ) return false ;         // vertex quality
  if( sqrt( TkPx[itk]*TkPx[itk]+TkPy[itk]*TkPy[itk]+TkPz[itk]*TkPz[itk] ) < 0.2 ) return false ;
  return true ;
}
// =============================================================================
double CPmix::CosAng( TLorentzVector *a , TLorentzVector *b )
{
  double PaPb =a->P()*b->P() ;
  Hfill("testa",a->P(),100,0.,5.,1.) ;
  Hfill("testb",b->P(),100,0.,5.,1.) ;
  if( PaPb == 0 ) return 0 ;
  return ( a->Px()*b->Px()+a->Py()*b->Py()+a->Pz()*b->Pz() )/PaPb ;
}
// =============================================================================
double CPmix::DeltaM( TLorentzVector *a , TLorentzVector *b )
{
  TLorentzVector c = (*a)+(*b) ;
  return c.M()-a->M() ;
}
// =============================================================================
double CPmix::MissingMass( TLorentzVector *a , TLorentzVector *b )
{
  TLorentzVector B( 0,0,0,5.27 ) ; // B meson at rest
  TLorentzVector Y = (*a)+(*b) ;
  TLorentzVector N = B-Y ;
  return N.M() ;
}
// =============================================================================
double CPmix::SquareMissingMass( TLorentzVector *a , TLorentzVector *b )
{
  TLorentzVector B( 0,0,0,5.27 ) ; // B meson at rest
  TLorentzVector Y = (*a)+(*b) ;
  TLorentzVector N = B-Y ;
  return N.M2() ;
}
// =============================================================================
void CPmix::analyze() 
{
  if( IsSgn != 1 ) return ;
  if( PrbPil < 0.001 ) return ;
  if( eVDec <= 0 ) return ;
  if( eVDec > 0.05 ) return ;

  double LepMass(0.0005 ) ;
  if( abs(LPId) == 2 ) LepMass = 0.105 ;
  int QLep = LPId / abs(LPId) ;
  mypart aLep( QLep, PxLep, PyLep, PzLep, LepMass, ZVDec, eVDec ) ;
  theLep = aLep ;

  Hfill("Plep",theLep.GetFourVector().P(),100,0.,3.,1.) ;

  TrackContainer.erase( TrackContainer.begin(), TrackContainer.end() ) ;
  TrackSource.erase( TrackSource.begin(), TrackSource.end() ) ;

  for( int itk = 0 ; itk< nTk ; itk++ ) 
    {
      mypart aTrack( GetCharge(itk), TkPx[itk], TkPy[itk], TkPz[itk], 
		     GetMass(itk), TkZVt[itk], TkeVt[itk] ) ;      

      if( abs(TkTag[itk]) == 666 ) { thePst = aTrack ; continue ;}      // soft pion
      if( !GoodTrack( itk )  ) continue ;

      TrackSource.push_back( GetTrackSource(itk) ) ;
      TrackContainer.push_back( aTrack ) ;
    }

  if( TrackContainer.size() <1 ) return ;
  CheckTruth() ;
}

// =============================================================================
void CPmix::CheckTruth( )
{
  mypart TagSide, DcySide ;
  int nTag (0) ;
  int nDcy (0) ;
  for( int itk = 0 ; itk < TrackContainer.size() ; itk++ )
    {								
      if( TrackSource.at(itk) > 2 ) { TagSide += TrackContainer.at(itk) ; nTag++ ; }
      else if( TrackSource.at(itk) > 0 ) { DcySide += TrackContainer.at(itk) ; nDcy++ ; }
    }

  MakeTrueHisto( TagSide, DcySide ) ;
}
// =============================================================================
void CPmix::MakeTrueHisto( mypart & TagSide, mypart & DcySide ) 
{
  mypart TagDst, DcyDst ;

  Hfill("PTag",TagSide.GetFourVector().P(),100,0.,5.,1.) ;
  Hfill("PDcy",DcySide.GetFourVector().P(),100,0.,5.,1.) ;

  int nTag = TagSide.GetMultiplicity() ;
  int nDcy = DcySide.GetMultiplicity() ;
  int nTot = min(12,nTag+nDcy) ;

  Hfill("nTag",nTag,20,0.,20.,1.) ;
  Hfill("nDcy",nDcy,20,0.,20.,1.) ;

  Hfill("NumDcyVsTot",(double) nDcy+nTag, 20,0.,20., (double) nDcy ,8,0.,8.,1. ) ;
  Hfill("NumDcy",nTot,(double) nDcy ,10,0.,10.,1. ) ;

  int xTag = min( nTag, iMaxTag ) ;
  int xDcy = min( nDcy, iMaxDcy ) ;
  Hfill("DcyVsTot",xDcy+xTag, iMaxDcy+iMaxTag,0.,dMaxDcy+dMaxTag, xDcy,iMaxDcy,0.,iMaxDcy,1. ) ;

  double ZVTag(999) ;
  double eVTag = TagSide.GetVertexError() ;

  if( eVTag > 0 ) 
    {
      ZVTag = TagSide.GetVertexPosition() ;
      double PulTag = ((zDec-zTag) - (ZVDec-ZVTag))/sqrt(eVTag*eVTag+eVDec*eVDec) ;   
      Hfill("TagZV",xTag,ZVDec-ZVTag,100,-0.2,0.2,1.) ;
      Hfill("TagPl",xTag,(ZVDec-ZVTag)/eVTag,100,-20.,20.,1.) ;
      Hfill("TageV",xTag,eVTag,50,0.,0.10,1.) ;
      Hfill("TagMs",xTag,TagSide.GetFourVector().M(),100,0.,6.,1.) ;
      Hfill("TagEn",xTag,TagSide.GetFourVector().E(),100,0.,6.,1.) ;
      Hfill("TagQ", xTag,(double) TagSide.GetCharge() , 10,-5.,5.,1. ) ;
      Hfill("TagCos",xTag,CosAng( &TagSide.GetFourVector(), &theLep.GetFourVector() ) , 100,-1.,1.,1. ) ;
      Hfill("TagDm" ,xTag,DeltaM( &TagSide.GetFourVector(), &thePst.GetFourVector() ), 100,0.140,0.290,1. ) ;
      Hfill("TagNuM",xTag,MissingMass( &TagSide.GetFourVector(), &theLep.GetFourVector()),100,-4.,4.,1.) ;
      Hfill("TagNuM2",xTag,SquareMissingMass( &TagSide.GetFourVector(), &theLep.GetFourVector()),100,-15.,10.,1.) ;

      double TagX2(0),TagPrb(-1),LTagPrb(-1)  ;
      if( nTag>1 ) 
	{ TagX2 =TagSide.GetChiSquare() ;
	  TagPrb=TMath::Prob( TagX2, nTag-1 ) ;
	  LTagPrb = -TMath::Log10(TagPrb) ;
	  if( LTagPrb>9.999 ) LTagPrb = 9.999 ;
	}
      Hfill("TagX2",xTag,TagX2,50,0.,10.,1.) ;
      Hfill("TagPrb",xTag,TagPrb,100,0.,1.,1.) ;
      Hfill("TagLPrb",xTag,LTagPrb,100,0.,10.,1.) ;
      Hfill("TagPull",xTag,PulTag,100,-10.,10.,1.) ;
      Hfill("TagCorr",zDec-zTag,50,-0.25,0.25,ZVDec-ZVTag,50,-0.25,0.25,1.) ;
      Hfill("TagRes",(ZVDec-ZVTag)-(zDec-zTag),100,-0.15,0.15,1.) ;

      TagDst += TagSide ;
      TagDst += thePst ;
      Hfill("TagDstNuM",xTag,MissingMass( &TagDst.GetFourVector(), &theLep.GetFourVector() ),100,-4.,4.,1.) ;
      Hfill("TagDstNuM2",xTag,SquareMissingMass( &TagDst.GetFourVector(), &theLep.GetFourVector() ),100,-5.,20.,1.) ;
    }
    
    double ZVDcy(999) ;
    double eVDcy = DcySide.GetVertexError() ;
    if( eVDcy > 0 )
      {
	ZVDcy = DcySide.GetVertexPosition() ;
	double PulDcy = ((zDec-zTag) - (ZVDec-ZVDcy))/sqrt(eVDcy*eVDcy+eVDec*eVDec) ;   
	Hfill("DcyZV",xDcy,ZVDec-ZVDcy,100,-0.2,0.2,1.) ;
	Hfill("DcyPl",xDcy,(ZVDec-ZVDcy)/eVDcy,100,-20.,20.,1.) ;
	Hfill("DcyeV",xDcy,eVDcy,50,0.,0.10,1.) ;
	Hfill("DcyMs",xDcy,DcySide.GetFourVector().M(),100,0.,2.5,1.) ;
	Hfill("DcyEn",xDcy,DcySide.GetFourVector().E(),100,0.,4.0,1.) ;
	Hfill("DcyQ", xDcy,(double) DcySide.GetCharge() , 10,-5.,5.,1. ) ;
	Hfill("DcyCos",xDcy,CosAng( &DcySide.GetFourVector(), &theLep.GetFourVector() ) , 100,-1.,1.,1. ) ;
	Hfill("DcyDm" ,xDcy,DeltaM( &DcySide.GetFourVector(), &thePst.GetFourVector() ), 100,0.140,0.290,1. ) ;
	Hfill("DcyNuM",xDcy,MissingMass( &DcySide.GetFourVector(), &theLep.GetFourVector()),100,-2.,4.,1.) ;
	Hfill("DcyNuM2",xDcy,SquareMissingMass( &DcySide.GetFourVector(), &theLep.GetFourVector()),100,-2.,12.,1.) ;
	
	double DcyX2(0),DcyPrb(-1),LDcyPrb(-1)  ;
	if( nDcy>1 ) 
	  { DcyX2 =DcySide.GetChiSquare() ;
	    DcyPrb=TMath::Prob( DcyX2, nDcy-1 ) ;
	    LDcyPrb = -TMath::Log10(DcyPrb) ;
	    if( LDcyPrb>9.999 ) LDcyPrb = 9.999 ;
	  }
	Hfill("DcyX2",xDcy,DcyX2,50,0.,10.,1.) ;
	Hfill("DcyPrb",xDcy,DcyPrb,100,0.,1.,1.) ;
	Hfill("DcyLPrb",xDcy,LDcyPrb,100,0.,10.,1.) ;
	Hfill("DcyPull",xDcy,PulDcy,100,-10.,10.,1.) ;
	Hfill("DcyCorr",zDec-zTag,50,-0.25,0.25,ZVDec-ZVDcy,50,-0.25,0.25,1.) ;
	Hfill("DcyRes",(ZVDec-ZVDcy)-(zDec-zTag),100,-0.15,0.15,1.) ;
	
	DcyDst += thePst ;	
	DcyDst += DcySide ;

	Hfill("DcyDstMs",xDcy,DcyDst.GetFourVector().M(),100,0.,2.5,1.) ;
	Hfill("DcyDstEn",xDcy,DcyDst.GetFourVector().E(),100,0.,4.0,1.) ;
	Hfill("DcyDstNuM",xDcy,MissingMass( &DcyDst.GetFourVector(), &theLep.GetFourVector() ),100,-4.,4.,1.) ;
	Hfill("DcyDstNuM2",xDcy,SquareMissingMass( &DcyDst.GetFourVector(), &theLep.GetFourVector() ),100,-5.,20.,1.) ;
      } 
    
    Hfill("Ppst",thePst.GetFourVector().P()*thePst.GetCharge(),100,-0.2,0.2,1.) ;
    Hfill("Mpst",thePst.GetFourVector().Mag(),100,0.1,0.2,1.) ;
    if( eVDcy > 0 && eVTag > 0 ) 
      {
	mypart TagB = TagSide ;
	mypart DcyB = DcySide ;
	DcyB += theLep ;
	DcyB += thePst ;
	Hfill( "DeltaE",TagB.GetFourVector().E()-DcyB.GetFourVector().E(),200,-5.,5.,1.) ; 
	Hfill( "CosB",CosAng( &TagB.GetFourVector(),&DcyB.GetFourVector() ),200,-1.,1.,1.) ;
	TLorentzVector B1B2 = TagB.GetFourVector() - DcyB.GetFourVector() ;
	Hfill("CosTheta",cos(B1B2.Theta()),200,-1.,1.,1.);
      }
}
inline void CPmix::SaveHist( char* FileName ) 
{
  cout << " Output to file " << FileName <<"\t" ;
  // TList* hList = gDirectory->GetList() ;
  TFile* f = new TFile(FileName,"RECREATE") ;
  for( unsigned int i = 0 ; i <ListTH2.size() ; i++ ) 
    {
      TH2F* h = ListTH2.at(i) ;
      cout << "(TH2) " <<  h->GetTitle() << "\t" ;
      NormalizeBySlice( h ) ;
      h->Write() ;
    }
  for( unsigned int i = 0 ; i <ListTH1.size() ; i++ ) 
    {
      TH1F* h = ListTH1.at(i) ;
      cout << "(TH1) " << h->GetTitle() << "\t" ;
      NormalizeHist( h ) ;
      h->Write() ;
    }
  f->Close() ;
  cout << endl ;
}
void CPmix::EndJob( char* OutputFileName) 
{
  // HSum a( "CosB" ) ;
  // a.SumFirstToLast("CosB_Int") ;
  // HSum b( "DeltaE" ) ;
  // b.SumToMaximum("DeltaE_Int") ;
  SaveHist( OutputFileName ) ;
}
void CPmix::NormalizeBySlice( TH2F* h )
{
  // cout << h->GetSum() << " " << h->GetEntries() << endl ;
  for( int i = 0 ; i<= h->GetNbinsX() ; i++ )
    {
      double sum(0) ;
      for( int j=0 ; j<= h->GetNbinsY() ; j++ ) 
	{ sum += h->GetBinContent(i,j) ; }
      if( sum == 0 ) continue ;
      for( int j=0 ; j<= h->GetNbinsY() ; j++ ) 
	{ 
	  double val = h->GetBinContent(i,j)/sum ;
	  double err = h->GetBinError(i,j)/sum ;
	  h->SetBinContent(i,j,val) ;
	  h->SetBinError(i,j,err) ;
	}
    }

}
void CPmix::NormalizeHist( TH1F* h )
{
  double sum(0) ;
  for( int j=0 ; j<= h->GetNbinsX() ; j++ ) 
    { sum += h->GetBinContent(j) ; }
  if( sum == 0 ) return ;
  for( int j=0 ; j<= h->GetNbinsX() ; j++ ) 
    { 
      double val = h->GetBinContent(j)/sum ;
      double err = h->GetBinError(j)/sum ;
      h->SetBinContent(j,val) ;
      h->SetBinError(j,err) ;
    }
}