//#include "SimpleCand.C"
#include <iostream>
#include <vector>
#include <map>

#include "TFile.h"
#include "TTree.h"
#include "TRandom3.h"
#include "TDatabasePDG.h"
#include "TParticlePDG.h"
#include "TCanvas.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TLegend.h"
#include "TParticle.h"
#include "TClonesArray.h"
#include "TF1.h"
#include "TStyle.h"
#include "TNtuple.h"
#include "TLatex.h"
#include "TVector.h"
#include "TMatrixD.h"
#include "TSystem.h"
#include "TMatrixDEigen.h"
#include "TMath.h"

#include "EventShape.h"


using std::cout;
using std::endl;

TRandom3      fRand;
TLorentzVector fIni;

// *** some program parameters
const int		fMAX	   = 100;  		// *** max tracks per event
int 			nbins 		= 300;      // *** bins of mass histograms
double 			low  		= 0.0;      // *** lower edge of histo
double 			high 		= 5.0;      // *** higher edge of histo (will be modified acc. to resonance mass)

// *** some contants for neutrals
const double    fNeutEff   = 0.9;      // *** efficiency for neutrals
const double    fdE        = 0.03;      // *** neutral energy resolution
const double    fNdtht     = 0.003;     // *** neutrals tht resolution
const double    fNdphi     = 0.003;     // *** neutrals phi resolution
const double    fNeutThrsh = 0.02;		// *** energy threshold of neutrals

// *** constants for particles
int				fPi0Pdg	   = 111;		// *** PDG code for pi0
double		 	fPi0Mass   = 0.13497;    // *** center of mass cut window
double			fPi0win	   = 0.020;

int				fKsPdg	   = 310;		// *** PDG code for pi0
double		 	fKsMass   = 0.497614;  // *** center of mass cut window
double			fKswin	   = 0.025;

int             fDPdg      = 411;         // *** PDG code to match D
double		 	fDMass     = 1.8693;      // *** center of mass cut window
double			fDsig      = 0.88;        // *** slope of linear width function
double			fDwin		=0.25;

int             fD0Pdg     = 421;         // *** PDG code to match D0
double		 	fD0Mass    = 1.8645;      // *** center of mass cut window
double			fD0sig     = 1.04;        // *** slope of linear width function
double			fD0win		=0.25;

int             fJpsiPdg   = 443;         // *** PDG code to match J/psi
double		 	fJMass     = 3.09687;        // *** center of mass cut window
double			fJsig      = 2.24;        // *** slope of linear width function
double			fJsig3pi   = 1.5;        // *** slope of linear width function for 3pi mode
double			fJ3pioff   = 0.027;        // *** offset of linear function for 3pi mode
double			fJwin		=0.7;

int             fLamPdg    = 3122;         // *** PDG code to match Lam
double		 	fLamMass   = 1.115684;        // *** center of mass cut window
double			fLamsig    = 0.08;        // *** slope of linear width function
double			fLamwin		=0.2;

int             fDsPdg     = 431;         // *** PDG code to match Ds
double		 	fDsMass    = 1.9686;        // *** center of mass cut window
double			fDssig     = 0.75;        // *** slope of linear width function
double			fDswin		=0.2;

int             fPhiPdg     = 333;         // *** PDG code to match Phi
double		 	fPhiMass    = 1.0195;        // *** center of mass cut window
double			fPhisig     = 0.15;        // *** slope of linear width function
double			fPhiwin		=0.2;

int             fEtacPdg     = 441;         // *** PDG code to match etac
double		 	fEtacMass    = 2.9798;        // *** center of mass cut window
double			fEtacsig     = 0.15;        // *** slope of linear width function
double			fEtacwin      =0.15;

int             fLamcPdg     = 4122;         // *** PDG code to match Lambda_c
double		 	fLamcMass    = 2.2849;        // *** center of mass cut window
double			fLamcsig     = 0.15;        // *** slope of linear width function
double			fLamcwin	  =0.2;


// *** for event shape variables (global)
double			_sph, _pla, _apl, _thr;
TVector3		_sphAx, _secAx, _thrAx;

// *** max order for Fox-Wolfram Moments
const int    	_nfwmom = 5;
// *** result container for computation (global)
double 			_foxwolfR[_nfwmom];

// *** pid probability table; will be modified according to purity parameter value
// ***  true pid:        e   mu    pi    K   p      classified as [%]
double pidtab [25] = {  80,   5,   5,   5,   5,    //  e
                        5,  80,   5,   5,   5,    // mu
                        5,   5,  80,   5,   5,    // pi
                        5,   5,   5,  80,   5,    // K
                        5,   5,   5,   5,  80
                     };  // p

// *** table determined by Donghee 11/2012 for P>0.1
double pidreal [25] = {  95,  6,  10,   7,  10,    //  e
                         5,  86,  27,   6,   7,    // mu
                         9,  29,  85,  18,  16,    // pi
                         9,  21,  40,  67,  24,    // K
                         12, 12,  20,  16,  90
                      };  // p

// *** global stable particle masses used in pid selection; set in init()
double mE, mMu, mPi, mK, mP, mPi0;

TF1 *f_res=0;

// *** this map contains mapping from PDG code to indices 0...4 (for pidtab)
// *** 11 -> 0, 13 -> 1, 211 -> 2, 321 -> 3, 2212 -> 4
std::map<int,int> pdgidx;

TDatabasePDG *pdg_db;

typedef std::vector<SimpleCand> CandList;

// *** all global particle lists we need to do combinatorics
CandList mclist, mcFSlist, recolist, recocmslist;

// *** the global pid lists
CandList eplus, eminus, muplus, muminus, piplus, piminus, kplus, kminus, pplus, pminus, gam;

// *** boost full list with -boost and store in boostlist
void boostList(CandList &list, CandList &boostlist, TVector3 boost)
{
	int N=list.size();
	boostlist.clear();

    for (unsigned int i=0;i<N;i++)
    {
        SimpleCand c=list[i];
		TLorentzVector lv(c.P4());
        lv.Boost(-boost);
		c.SetP4(lv);
		boostlist.push_back(c);
    }
}

// *** aux for thrust; eps(x)=sign(x)
double eps(TVector3 v1, TVector3 v2)
{
  double prod = v1.Dot(v2);
  return prod/fabs(prod);
}

// *** routine to compute thrust vector and thrust
double compThrust(CandList& list, TVector3 &axis)
{
  TVector3 n0;
  int N=list.size();
  
  int i,j;
  double pmax=0;
  
  //find starting vector as maximum momentum vector
  for (i=0;i<N;++i)
  {
	if (list[i].P4().Vect().Mag()>pmax)
	{
	  n0=list[i].P4().Vect();
	  pmax=list[i].P4().Vect().Mag();
	}
  }
  //n0=list[0].P4().Vect();
  
  TVector3 nnew(0,0,0);
  // find thrust axis (10 iterations)
  for (i=0;i<10;++i)
  {
	for (j=0;j<N;++j)
	{
	  //nnew+=eps(n0,list[j].P4().Vect())*list[j].P4().Vect();
	  nnew+=n0*list[j].P4().Vect()>0 ? list[j].P4().Vect() : -list[j].P4().Vect();
	}
	n0=nnew.Unit();
  }
  
  axis = n0;
  double thr=0, sum=0;
  for (i=0;i<N;++i)
  {
	thr+=fabs(axis.Dot(list[i].P4().Vect()));
	sum+=list[i].P4().Vect().Mag();
  }
  return thr/sum;
}

// *** Legendre Function; auxiliary function for Fox-Wolfram-Moments
double Legendre( int l, double x )
{
    assert(fabs(x) <= 1.);

	if (l==0) return 1.;
	
    double pmm = 1.;
	double pmmp1 = x;

	if (l>1)
	{
		for(int ll=2; ll<=l; ll++)
		{
		  double pll = (x * (2 * ll - 1) * pmmp1 - (ll - 1) * pmm) / (ll);
		  pmm = pmmp1;
		  pmmp1 = pll;
		}
	}
   return pmmp1;
}

// *** compute relative Fox-Wolfram-Moments R up to 4th order (see _nfwmom contant above)
void compFoxWolfMoms(CandList& list)
{
  int n=list.size();
  if( n==0 ) return;

  double s = 0.;
  double _sumarray[_nfwmom]={0,0,0,0};
  int i,j,l;

  for (i=0; i<n-1; ++i)
  {
     // this candidate's 3-momentum
	  TVector3 p1(list[i].P4().Vect());
      double pmag1 = p1.Mag();

      // loop over other candidates, starting at the next one in the list
 	  
	  for (j=i+1; j<n; ++j)
	  {
		// this candidate's 3-momentum
		TVector3 p2(list[j].P4().Vect());
		double pmag2 = p2.Mag();

		// the cosine of the angle between the two candidates
		double cosPhi =  cos ( p1.Angle(p2) );

		// the contribution of this pair of track 
		// (note the factor 2 : the pair enters the sum twice)
		for( l=0; l<_nfwmom; l++ )
		  _sumarray[l] += 2 * pmag1 * pmag2 * Legendre( l, cosPhi );
		  
	  }
      // contribution of this track
      for( l=0; l<_nfwmom; l++ )
		_sumarray[l] += pmag1 * pmag1 * Legendre( l, 1. );
      
      // total energy
      s += list[i].P4().Energy();
   }
  
  // well ...
  if( s<=0. ) return;
  double ssq=s*s;
  // normalize Fox Wolfram Moments
  for( i=0; i<_nfwmom; i++)
    _foxwolfR[i] = _sumarray[i]/ssq ;
}

// *** routine to compute sphericity, planarity and aplanarity
double compSphericity(CandList &list)
{
  int N=list.size();
  
  double stot=0, sxx=0, sxy=0, sxz=0, syy=0, syz=0, szz=0;	
  int i;
  
  for (i=0;i<N;++i)
  {
	TVector3 v(list[i].P4().Vect());
	
	sxx += v.X()*v.X(); sxy += v.X()*v.Y(); sxz += v.X()*v.Z();
	syy += v.Y()*v.Y(); syz += v.Y()*v.Z();
	szz += v.Z()*v.Z();
	
	stot += v.Mag2();
  }
  
  TMatrixD sm(3,3);
  sm(0,0) = sxx/stot; sm(0,1) = sxy/stot; sm(0,2) = sxz/stot; 
  sm(1,0) = sxy/stot; sm(1,1) = syy/stot; sm(1,2) = syz/stot; 
  sm(2,0) = sxz/stot; sm(2,1) = syz/stot; sm(2,2) = szz/stot; 
  
  TMatrixDEigen ei(sm);
  TMatrixD eiv=ei.GetEigenValues();
  
  _sph = 1.5 * (eiv(1,1) + eiv(2,2));
  _apl = 1.5 * eiv(2,2);
  _pla = eiv(1,1) - eiv(2,2);
  
  return _sph;
}

// *** check PDG code is final state particle
bool isFinalState(int n)
{
    n=abs(n);
    return (n==11 || n==13  || n==211 || n==321 || n==2212 || n==22);
}

// *** online momentum resoltion based on talk of Yutie Lang, PANDA CM 12/2011, TRK oder DAQ session
// *** values (pt[GeV/c], dp/p[%]): (1.0, 3.8) (2.0, 5.5) (5.0, 12.2)
// *** see init routine for function parameters
double dponline(TLorentzVector &v)
{
    return f_res->Eval(v.Pt())/100.;
}

// *** computes the pid table for given efficiency and purity
void makePidTable(double pideff, double misid)
{
    cout<< "Using PID table:"<<endl;

    if (pideff<0 && misid<0)
    {
        for (int i=0;i<5;++i)
        {
            for (int j=0;j<5;++j)
            {
                //pidtab[j+5*i] = pidreal[j+5*i];
                printf("%4.0f",pidtab[j+5*i]);
                //cout << pidtab[i+5*j]<<"  ";
            }
            cout<< endl;
        }
        return;
    }

    if (pideff<0)   pideff=0;
    if (pideff>100) pideff=100;
    if (misid>100) misid=100;

    // *** neff is the efficiency for particles of wrong type
    //double neff = pideff*(100-pidpur)/(4*pidpur);
    //if (neff>100) neff=100;

    for (int i=0;i<5;++i)
    {
        for (int j=0;j<5;++j)
        {
            if (i==j) pidtab[i+5*j] = pideff;
            else pidtab[i+5*j] = misid;
            cout << pidtab[i+5*j]<<"  ";
        }
        cout<< endl;
    }
}

// *** initialization of some globals
void init(double pideff, double misid)
{
    // *** setting the mapping of PDG code to tab index
    pdgidx.clear();
    pdgidx[11]=0;
    pdgidx[13]=1;
    pdgidx[211]=2;
    pdgidx[321]=3;
    pdgidx[2212]=4;

    pdg_db=TDatabasePDG::Instance();

    // *** create pid table
    makePidTable(pideff, misid);

    // *** gets particle masses from pdg database and set them globally
    mE  = pdg_db->GetParticle(11)->Mass();
    mMu = pdg_db->GetParticle(13)->Mass();
    mPi = pdg_db->GetParticle(211)->Mass();
    mK  = pdg_db->GetParticle(321)->Mass();
    mP  = pdg_db->GetParticle(2212)->Mass();

    mPi0= pdg_db->GetParticle(111)->Mass();

    // *** resolution fcn based on Yutie Lang results; see routine dponline above
    f_res=new TF1("f_res","pol2(0)");
    f_res->SetParameters(2.367,1.3,0.133);
    f_res->SetRange(0,10);
}

// *** prints single candidate of type SimpleCand
void printCand(SimpleCand c)
{
    cout <<c.Id()<<" pid:"<<c.Pid()<<" ("<<c.P4().X()<<","<<c.P4().Y()<<","<<c.P4().Z()<<";"<<c.P4().T()<<") ";
    cout <<c.Charge()<<" "<<c.P4().M();
    cout <<" moth:"<<c.MotherIdx()<<" nsib:"<<c.NSiblings()<<" nfs:"<<c.NFS();
    cout <<" mcpid:"<<c.McPid()<<" mct:"<<c.Mct()<<endl;
}

// *** prints list of candidates
void printList(CandList &l, TString tit="")
{
    if (tit!="") cout << tit<<endl;
    for (unsigned int j=0;j<l.size();++j) printCand(l[j]);
    cout <<endl;
}

// *** prints all global lists
void printAllLists()
{
    printList(mclist, "MC List");
    printList(recolist, "Reco List");

    printList(eplus, "e+ List");
    printList(eminus, "e- List");

    printList(muplus, "mu+ List");
    printList(muminus, "mu- List");

    printList(piplus, "pi+ List");
    printList(piminus, "pi- List");

    printList(kplus, "k+ List");
    printList(kminus, "k- List");

    printList(pplus, "p+ List");
    printList(pminus, "p- List");

    printList(gam, "gamma List");
}

// *** counts how many combinations (=composite particles in l) fullfill the mass criterion
int massCrit(CandList &l, double m, double w)
{
    int cnt = 0;
    for (unsigned int i=0;i<l.size();++i)
    {
        double mass=l[i].P4().M();
        if (mass>m-w && mass<m+w) cnt++;
    }

    return cnt;
}

// *** fill the mass histos
int fillHisto(CandList &l, TH1 *hall, TH1* hbg=0, double m =0, double w=0,  TH1* hallsel=0, TH1* hbgsel=0, TH1* hsig=0, TH1* hmm=0, TH2* hmvsmm=0)
{
    int cnt =0;
    for (unsigned int i=0;i<l.size();++i)
    {
        double mass = l[i].P4().M();
        double miss = (fIni-l[i].P4()).M();
        hall->Fill(mass);

        //missing mass
        if (hmm) hmm->Fill(miss);
        if (hmvsmm) hmvsmm->Fill(mass, miss);

        bool mct = l[i].Mct();
        if (hbg && !mct) hbg->Fill(mass);
        if (hsig && mct) hsig->Fill(mass);
        if (mass>m-w && mass<m+w)
        {
            cnt++;
            if (hallsel) hallsel->Fill(mass);
            if (hbgsel && !mct) hbgsel->Fill(mass);
        }
    }
    return cnt;
}

int combine(CandList &l1, CandList &l2, CandList &out, int matchPdg=0)
{
    out.clear();

    bool sameList = (&l1 == &l2);

    for (unsigned int i=0;i<l1.size();++i)
    {
        for (unsigned int j=0;j<l2.size();++j)
        {
            if (sameList && j<=i) continue;
            if (l1[i].Marker() & l2[j].Marker()) continue;  // *** overlap
            //if (l1[i].EvtId()!=l2[j].EvtId()) continue;     // *** tracks from different events
            SimpleCand c;
            c.SetP4(l1[i].P4()+l2[j].P4());
            c.SetCharge(l1[i].Charge()+l2[j].Charge());
            c.SetNFS(l1[i].NFS()+l2[j].NFS());
            c.SetMarker(l1[i].Marker() | l2[j].Marker());

            int m1 = l1[i].MotherIdx();
            int m2 = l2[j].MotherIdx();
            int ev1 = l1[i].EvtId();
            int ev2 = l2[j].EvtId();

            if ( ev1==ev2 &&							   // *** particles are from same event
                    m1>=0 && m1 == m2 &&                      // *** there is a common mother of all particles
                    (l1[i].McPid()||l1[i].Mct()) &&            // *** all particles match (FS: PID, Comp:mct)
                    (l2[j].McPid()||l2[j].Mct()) &&            // *** "
                    2==mclist[m1].NDau() &&            // *** mother has as many daughters as we combine particles
                    matchPdg == mclist[m1].Pid() )             // *** mother has correct PDG
            {
                c.SetId(m1);
                c.SetMcPid();
                c.SetMct();
                c.SetMotherIdx(mclist[m1].MotherIdx());
            }

            out.push_back(c);
        }
    }
    return out.size();
}

int combine(CandList &l1, CandList &l2, CandList &l3, CandList &out, int matchPdg=0)
{
    out.clear();

    bool same12 = (&l1 == &l2);
    bool same13 = (&l1 == &l3);
    bool same23 = (&l2 == &l3);

    for (unsigned int i1=0;i1<l1.size();++i1)
    {
        for (unsigned int i2=0;i2<l2.size();++i2)
        {
            if (same12 && i2<i1) continue;
            if (l1[i1].Marker() & l2[i2].Marker()) continue;

            for (unsigned int i3=0;i3<l3.size();++i3)
            {
                if ((same13 && i3<i1) || (same23 && i3<i2)) continue;

                if ( l1[i1].Marker() & l3[i3].Marker() ) continue;
                if ( l2[i2].Marker() & l3[i3].Marker() ) continue;

                //if (l1[i1].EvtId()!=l2[i2].EvtId() ||
                //    l1[i1].EvtId()!=l3[i3].EvtId() ||
                //    l2[i2].EvtId()!=l3[i3].EvtId()   ) continue;     // *** tracks from different events

                SimpleCand c;
                c.SetP4( l1[i1].P4() + l2[i2].P4() + l3[i3].P4());
                c.SetCharge( l1[i1].Charge() + l2[i2].Charge() + l3[i3].Charge() );
                c.SetNFS( l1[i1].NFS() + l2[i2].NFS() + l3[i3].NFS());
                c.SetMarker( l1[i1].Marker() | l2[i2].Marker() | l3[i3].Marker() );

                int m1 = l1[i1].MotherIdx();
                int m2 = l2[i2].MotherIdx();
                int m3 = l3[i3].MotherIdx();

                int ev1 = l1[i1].EvtId();
                int ev2 = l2[i2].EvtId();
                int ev3 = l3[i3].EvtId();

                if ( ev1==ev2 && ev2==ev3 &&				   // *** particles are from same event
                        m1>=0 && m1 == m2 && m2 == m3 &&          // *** there is a common mother of all particles
                        (l1[i1].McPid()||l1[i1].Mct()) &&          // *** all particles match (FS: PID, Comp:mct)
                        (l2[i2].McPid()||l2[i2].Mct()) &&          // *** "
                        (l3[i3].McPid()||l3[i3].Mct()) &&          // *** "
                        3==mclist[m1].NDau() &&            // *** mother has as many daughters as we combine particles
                        matchPdg == mclist[m1].Pid() )             // *** mother has correct PDG
                {
                    c.SetId(m1);
                    c.SetMcPid();
                    c.SetMct();
                    c.SetMotherIdx(mclist[m1].MotherIdx());
                }

                out.push_back(c);
            }
        }
    }

    return out.size();
}

void smearMom(SimpleCand &c, double dpr=0.05, double dtht=0.000, double dphi=0.000)
{
    TLorentzVector p4=c.P4();

    // *** use resolution function
    if (dpr<0) dpr = dponline(p4);

    double p    = p4.P();
    double m    = p4.M();

    double dp   = p*dpr;
    double newp = -1;

    if (dp==0) newp=p;
    while (newp<0) newp = p4.P() + fRand.Gaus(0.0,dp);

    p4.SetVectM( p4.Vect().Unit()*newp, m );

    double newtht = p4.Theta();
    if (dtht!=0) newtht += fRand.Gaus(0.,dtht);
    p4.SetTheta(newtht);

    double newphi = p4.Phi();
    if (dphi!=0) newphi += fRand.Gaus(0.,dphi);
    p4.SetPhi(newphi);

    c.SetP4(p4);

}

bool isDetected(double eff=0.90)
{
    // very simple efficiency cut; might be dependent on candidate itself later
    if (fRand.Rndm()<eff) return true;
    return false;
}

void makeRecoCands(CandList &mc, CandList &reco, double eff=0.90, double dp=0.05, double dtht=0.001, double dphi=0.001)
{
    reco.clear();

    for (unsigned int i=0;i<mc.size();i++)
    {
        SimpleCand c=mc[i];
        int pid = abs(c.Pid());
        if (!isFinalState(pid)) continue;

        if (pid!=22 && isDetected(eff))  // *** charged track
        {
            smearMom(c,dp,dtht,dphi);
			// set pion mass as default
			c.SetMass(0.13957);
            reco.push_back(c);
        }
        else if (pid==22 && isDetected(fNeutEff) && c.P4().E()>fNeutThrsh)  // *** neutral
        {
            smearMom(c,fdE, fNdtht, fNdphi);
            reco.push_back(c);
        }
    }
}

// **** A simple mass selector
void selectMass(CandList &in, CandList &out, double mmean=0, double mw=0)
{
    out.clear();

    for (unsigned int i=0;i<in.size();++i)
        if (fabs(in[i].P4().M()-mmean)<mw) out.push_back(in[i]);

}

// *** perform pid selection according to the pidtable
void select(CandList &in, CandList &out, int chrg=0, int pdg=0, double mass=0.139)
{
    out.clear();

    TParticlePDG *part = pdg_db->GetParticle(pdg);

    //pdg *= pdg<20 ? -chrg : chrg;

    if (part) mass = part->Mass();
    if (mass<0) mass = 0.13957;

    for (unsigned int i=0;i<in.size();++i)
    {
        SimpleCand c=in[i];

        // select a specific charge in any case
        if (c.Charge()==chrg)
        {
            c.SetMass(mass);
            if (c.Pid()==pdg) c.SetMcPid();
            if (pdg==0 || pdg==22) out.push_back(c); // do we want select a pdg code? 0=only charge
            else
            {
                // invalid pdg code?
                if (pdgidx.find(abs(pdg))==pdgidx.end()) continue;

                int selidx   = pdgidx[abs(pdg)];
                int partidx  = pdgidx[abs(c.Pid())];
                double prob  = pidtab[5*selidx + partidx]/100.;
                //cout <<pdg <<" "<<c.Pid()<<" : "<<selidx<<" "<<partidx<<" P:"<<prob<<endl;

                if (fRand.Rndm()<prob) out.push_back(c);
            }
        }
    }
}

void makePidSelection(CandList &l)
{
    // *** prepare globally defined PID lists
    select(l, eplus, 1, -11, mE);
    select(l, eminus, -1, 11, mE);

    select(l, muplus, 1, -13, mMu);
    select(l, muminus, -1, 13, mMu);

    select(l, piplus, 1, 211, mPi);
    select(l, piminus, -1, -211, mPi);

    select(l, kplus, 1, 321, mK);
    select(l, kminus, -1, -321, mK);

    select(l, pplus, 1, 2212, mP);
    select(l, pminus, -1, -2212, mP);

    select(l, gam, 0, 22, 0);
}

void makeIni4Vector(TLorentzVector &l, double s)
{
    double X = (s*s-2*mP*mP)/(2*mP);
    double p = sqrt(X*X - mP*mP);
    l.SetXYZM(0,0,p,s);
}


void softtrigger_kin5(TString fsig, int mode=0, double sqrts=3.77, int nev=1000, bool evcut=false, double dp=0.03,
                     double trkeff=100., double pideff=95., double misid =5., double P_mix=0.00 )
{
    fRand.SetSeed();

    init(pideff, misid);

    unsigned int i,j,k;

    makeIni4Vector(fIni, sqrts);

    //gStyle->SetOptStat(1111);

    // *** open file and get tree
    TFile *f=new TFile(fsig,"READ");
    if (f->IsZombie())
    {
        cout <<"File not found:"<<fsig<<endl;
        return;
    }

    bool dpm=fsig.Contains("DPM");            // *** DPM events
	bool bkg=dpm || fsig.Contains("BKG");     // *** Background events from EvtGen OR DPM events
    TString treename = dpm ? "data" : "ntp";  // *** set treename
    TTree *t=(TTree*)f->Get(treename);

    if (t==0)
    {
        cout <<"Tree '"<<treename<<"' not found in '"<<fsig<<"'."<<endl;
        return;
    }

    // *** how many events ?
    int Nevt = t->GetEntriesFast();
    if (nev==0) nev = Nevt;

    double mmean=0, mwsig=0, mwoff=0;
    int matchPDG=0;

    // *** choose histogram scale according to resonance
    high = mmean+2;

    double min=0, max=5.5;

    TString name = "ntp_";
    if (!bkg) name+="sig_mode_";
    else name+="bkg_mode_";
    name+=mode; name+=".root";
	
	TFile *ggg=new TFile(name, "RECREATE");

	//branches
    //                   0   1     2     3     4   5   6   7  8   9   10   11    12   13   14   15   16   17   18   19   20    21    22    23    24    25    26    27   28  29  30  31  32  33
	TString branches = "sig:mode:npart:nneut:nchrg:ne:nmu:npi:nk:npr:pmax:ptmax:np05:np10:np15:np20:np25:np30:np40:np50:npt05:npt10:npt15:npt20:npt25:npt30:npt40:npt50:sph:apl:pla:thr:sne:scp:";

    //                  34  35  36   37     38    39    40    41    42    43    44    45    46    47   48   49    50    51    52    53   54   55    56
	branches        += "fw1:fw2:fw3:pmaxl:ptmaxl:np05l:np10l:np15l:np20l:np25l:np30l:np40l:np50l:ne10:ne20:nmu10:nmu20:npi10:npi20:nk10:nk20:npr10:npr20:";
	
    //                   57     58     59     60     61      62      63    64    65     66    67  68  69  70   71    72     73      74      75      76      77       78      79
	branches        += "sumpt:sumptl:sumptc:sumptcl:sumetn:sumetnl:sumpc:sumpcl:sumen:sumenl:cir:fw4:fw5:pmin:ptmin:pminl:prapmax:sumpt05:sumpt10:sumpc05:sumpc10:sumpc05l:sumpc10l:";
    
	//                    80        81       82      83    84   85    86     87     88     89
	branches        += "sumen05:sumen10:sumen05l:sumen10l:npi0:nks0:npi005:nks005:npi010:nks010";
	
    TNtuple *ntpev = new TNtuple("ntpev","ntpev", branches);

    Float_t ntpval[100];

    // *** prepare the tree to loop through

    Int_t nTrk,Id[fMAX],DF[fMAX],DL[fMAX];
    Double_t px[fMAX],py[fMAX],pz[fMAX],E[fMAX],m[fMAX];
    TClonesArray* DpmPart=0;

    if (dpm)
    {
        DpmPart=new TClonesArray("TParticle",100);
        t->SetBranchAddress("Npart",&nTrk);
        t->SetBranchAddress("Particles", &DpmPart);
    }
    else
    {
        t->SetBranchAddress("nTrk",&nTrk);
        t->SetBranchAddress("Id",Id);
        t->SetBranchAddress("DF",DF);
        t->SetBranchAddress("DL",DL);
        t->SetBranchAddress("px",px);
        t->SetBranchAddress("py",py);
        t->SetBranchAddress("pz",pz);
        t->SetBranchAddress("E",E);
        t->SetBranchAddress("m",m);
    }


    // *** the event loop

    int mixevts=1;
    int currmix=0;

    // how many events are tagged by individual algo?
    bool channelselect[6];
    double channelcount[6];
    TString channelname[6]={"J/psi","D+-","Ds","Phi","D0","Lamc"};

    for (i=0;i<6;++i) channelcount[i]=0.;

    if (fRand.Rndm()<P_mix) mixevts=2;

    for (i=0;i<nev;++i)
    {
        t->GetEntry(i);

        if (i%5000==0) cout <<"ev "<<i<<endl;

        // *** read in from file the mc candidate
        if (0==currmix)
        {
            mclist.clear();   // *** contains all MC cands (also non-final-states)
            mcFSlist.clear(); // *** contains only final state MC cands
        }

        bool good = true;

        if (dpm)
        {
            for (j=0;j<nTrk;++j)
            {
                TParticle *pt=(TParticle*)DpmPart->At(j);
                Id[j] = pt->GetPdgCode();
                int pdg = abs(Id[j]);
                px[j] = pt->Px();
                py[j] = pt->Py();
                pz[j] = pt->Pz();
                E[j]  = pt->Energy();
                m[j]  = pt->GetMass();
                DF[j] = 0;
                DL[j] = -1;
            }
        }

		for (j=0;j<nTrk;++j)
        {
            int chrg = 0;
            int pdg = abs(Id[j]);
            TParticlePDG *p=pdg_db->GetParticle(Id[j]);
            if (p) chrg = p->Charge();
            if (abs(chrg)>1) chrg/=3;

            SimpleCand c(px[j], py[j], pz[j], E[j], Id[j], j, chrg);
            c.SetDau(DF[j], DL[j]);

            for (k=0;k<j;++k) if (j>=DF[k] && j<=DL[k])
                {
                    c.SetMotherIdx(k);
                    c.SetNSiblings(DL[k]-DF[k]+1);
                }

            c.SetEvtId(i);
			
            mclist.push_back(c);

            // final state particle?
            if (isFinalState(pdg)) mcFSlist.push_back(c);
        }

        if (++currmix<mixevts) continue;

        currmix=0;

		// *** the boost of the system
		TVector3 boost = fIni.BoostVector();

        // *** prepare the toy reco: smear momentum and apply track efficiency (both flat at the moment)
        makeRecoCands(mcFSlist, recolist, trkeff/100., dp );

/*		int nneut = 0;
		int nchrg = 0;

		for (j=0;j<recolist.size();++j)
		{
			int chrg = recolist[j].Charge();
			if (abs(chrg)==0) nneut++; else nchrg++;
		}
*/		
		boostList(recolist, recocmslist, boost);

        // *** prepare all pid list according to the pidtable probabilities
        makePidSelection(recolist);
		
		EventShape evs(recolist, fIni);
		
		// *** the lists we need for combinatorics
		CandList comb1;
		CandList LPi0, LKs;				// pi0 -> gg (mass cut), KS -> pi+ pi- (mass cut)
		
		// *** do combinatorics (pi0 -> gg) 
		combine(gam, gam, comb1, fPi0Pdg);
		//fillHisto(comb1, h_pi0, 0, fPi0Mass, fPi0win, h_pi0sel, 0, h_pi0sig ); 
		selectMass(comb1, LPi0, fPi0Mass, fPi0win);
		
		// *** do combinatorics (pi0 -> gg) 
		combine(piplus, piminus, comb1, fKsPdg);
		//fillHisto(comb1, h_ks, 0, fKsMass, fKswin, h_kssel, 0, h_kssig ); 
		selectMass(comb1, LKs, fKsMass, fKswin);
		
		// *** Make Event Selection
// 		if (evcut)
// 		{
// 			EventShape evshape(recolist, fIni);
// 			if (evshape.ChrgPtSumLab()<1.5) continue;
// 		}

        //        0   1     2     3     4   5   6   7  8  9   10   11    12   13   14   15   16   17   18   19   20    21    22    23    24    25    26    27   28  29  30  31  32  33
	    //      "sig:mode:npart:nneut:nchrg:ne:nmu:npi:nk:npr:pmax:ptmax:np05:np10:np15:np20:np25:np30:np40:np50:npt05:npt10:npt15:npt20:npt25:npt30:npt40:npt50:sph:apl:pla:thr:sne:scp";

        //       34  35  36   37     38    39    40    41    42    43    44    45    46    47   48   49    50    51    52    53   54   55    56
	    //      "fw1:fw2:fw3:pmaxl:ptmaxl:np05l:np10l:np15l:np20l:np25l:np30l:np40l:np50l:ne10:ne20:nmu10:nmu20:npi10:npi20:nk10:nk20:npr10:npr20");
        
        //        57     58     59     60     61      62      63    64    65     66    67  68  69  70   71    72     73      74      75      76      77       78      79
	    //      "sumpt:sumptl:sumptc:sumptcl:sumetn:sumetnl:sumpc:sumpcl:sumen:sumenl:cir:fw4:fw5:pmin:ptmin:pminl:prapmax:sumpt05:sumpt10:sumpc05:sumpc10:sumpc05l:sumpc10l:";
    
	    //       80        81       82      83    84   85    86     87     88     89
	    //   "sumen05:sumen10:sumen05l:sumen10l:npi0:nks0:npi005:nks005:npi010:nks010";
		
		ntpval[0]  = !bkg;                          // signal or background events
		ntpval[1]  = mode;                          // decay mode (parameter of script - set by hand)
		
		ntpval[2]  = evs.NParticles();              // #particles (multiplicity)
		ntpval[3]  = evs.NNeutral();                // #neutrals cands
		ntpval[4]  = evs.NCharged();                // #charged cands
		
		ntpval[5]  = eplus.size()+eminus.size();    // #electron candidates
		ntpval[6]  = muplus.size()+muminus.size();  // #muon candidates
		ntpval[7]  = piplus.size()+piminus.size();  // #pion candidates
		ntpval[8]  = kplus.size()+kminus.size();    // #kaon candidates
		ntpval[9]  = pplus.size()+pminus.size();    // #proton candidates
		ntpval[84] = LPi0.size();					// #pi0 candidates
		ntpval[85] = LKs.size();					// #KS0 candidates
		
		ntpval[37] = evs.PmaxLab();                 // maximum momentum (lab)
		ntpval[38] = evs.Ptmax();                   // maximum pt (lab)
		ntpval[72] = evs.PminLab();                 // minimum momentum (lab)
		ntpval[71] = evs.Ptmin();                   // minimum pt (lab)
		ntpval[73] = evs.PRapmax();                 // maximum pseudo rapidity (lab)
		
		ntpval[39] = evs.MultPminLab(0.5);          // #particles with p>0.5 GeV (lab)
		ntpval[40] = evs.MultPminLab(1.0);          // #particles with p>1.0 GeV (lab)
		ntpval[41] = evs.MultPminLab(1.5);          // #particles with p>1.5 GeV (lab)
		ntpval[42] = evs.MultPminLab(2.0);          // #particles with p>2.0 GeV (lab)
		ntpval[43] = evs.MultPminLab(2.5);          // #particles with p>2.5 GeV (lab)
		ntpval[44] = evs.MultPminLab(3.0);          // #particles with p>3.0 GeV (lab)
		ntpval[45] = evs.MultPminLab(4.0);          // #particles with p>4.0 GeV (lab)
		ntpval[46] = evs.MultPminLab(5.0);          // #particles with p>5.0 GeV (lab)
		
		ntpval[58] = evs.PtSumLab();				// sum pt of all particles (lab)
		ntpval[60] = evs.ChrgPtSumLab();			// sum pt of chrg particles (lab)
		ntpval[62] = evs.NeutEtSumLab();			// sum Et of neut particles (lab)
		ntpval[64] = evs.ChrgPSumLab();				// sum p of chrg particles (lab)
		ntpval[66] = evs.NeutESumLab();				// sum E of neut particles (lab)

		ntpval[78] = evs.SumChrgPminLab(0.5);		// sum p of charged particles with p>0.5 (lab)
		ntpval[79] = evs.SumChrgPminLab(1.0);		// sum p of charged particles with p>1.0 (lab)
		
		ntpval[82] = evs.SumNeutEminLab(0.5);		// sum E of neutral particles with E>0.5 (lab)
		ntpval[83] = evs.SumNeutEminLab(1.0);		// sum E of neutral particles with E>1.0 (lab)

		ntpval[58] = evs.PtSumLab();				// sum pt of all particles (lab)
		ntpval[74] = evs.SumPtminLab(0.5);			// sum pt of all particles with pt>0.5 GeV
		ntpval[75] = evs.SumPtminLab(1.0);			// sum pt of all particles with pt>1.0 GeV
					 
		// all following values are in center-of-momentum system
		ntpval[10] = evs.PmaxCms();                 // maximum momentum
		ntpval[11] = evs.Ptmax();                   // maximum pt
		ntpval[70] = evs.PminCms();                 // minimum momentum
		
		ntpval[12] = evs.MultPminCms(0.5);          // #particles with p>0.5 GeV
		ntpval[13] = evs.MultPminCms(1.0);          // #particles with p>1.0 GeV
		ntpval[14] = evs.MultPminCms(1.5);          // #particles with p>1.5 GeV
		ntpval[15] = evs.MultPminCms(2.0);          // #particles with p>2.0 GeV
		ntpval[16] = evs.MultPminCms(2.5);          // #particles with p>2.5 GeV
		ntpval[17] = evs.MultPminCms(3.0);          // #particles with p>3.0 GeV
		ntpval[18] = evs.MultPminCms(4.0);          // #particles with p>4.0 GeV
		ntpval[19] = evs.MultPminCms(5.0);          // #particles with p>5.0 GeV
				
		ntpval[57] = evs.PtSumCms();				// sum pt of all particles 
		ntpval[59] = evs.ChrgPtSumCms();			// sum pt of chrg particles 
		ntpval[61] = evs.NeutEtSumCms();			// sum Et of neut particles 
		ntpval[63] = evs.ChrgPSumCms();				// sum p of chrg particles
		ntpval[65] = evs.NeutESumCms();				// sum E of neut particles
		
		ntpval[76] = evs.SumChrgPminCms(0.5);		// sum p of charged particles with p>0.5 (cms)
		ntpval[77] = evs.SumChrgPminCms(1.0);		// sum p of charged particles with p>1.0 (cms)
		
		ntpval[80] = evs.SumNeutEminCms(0.5);		// sum E of neutral particles with E>0.5 (cms)
		ntpval[81] = evs.SumNeutEminCms(1.0);		// sum E of neutral particles with E>1.0 (cms)
		
		ntpval[20] = 0.0;                           // #particles with pt>0.5 GeV
		ntpval[21] = 0.0;                           // #particles with pt>1.0 GeV
		ntpval[22] = 0.0;                           // #particles with pt>1.5 GeV
		ntpval[23] = 0.0;                           // #particles with pt>2.0 GeV
		ntpval[24] = 0.0;                           // #particles with pt>2.5 GeV
		ntpval[25] = 0.0;                           // #particles with pt>3.0 GeV
		ntpval[26] = 0.0;                           // #particles with pt>4.0 GeV
		ntpval[27] = 0.0;                           // #particles with pt>5.0 GeV

		ntpval[47] = 0.0;							// #electrons with p>1.0
		ntpval[48] = 0.0;							// #electrons with p>2.0

		ntpval[49] = 0.0;							// #muons with p>1.0
		ntpval[50] = 0.0;							// #muons with p>2.0
		
		ntpval[51] = 0.0;							// #pions with p>1.0
		ntpval[52] = 0.0;							// #pions with p>2.0

		ntpval[53] = 0.0;							// #kaons with p>1.0
		ntpval[54] = 0.0;							// #kaons with p>2.0
		
		ntpval[55] = 0.0;							// #protons with p>1.0
		ntpval[56] = 0.0;							// #protons with p>2.0
		
		ntpval[86] = 0.0;							// #pi0s with p>0.5
		ntpval[88] = 0.0;							// #pi0s with p>1.0
		
		ntpval[87] = 0.0;							// #KSs with p>0.5
		ntpval[89] = 0.0;							// #KSs with p>1.0
		
		//compSphericity(recocmslist);
		
		//TVector3 thrAxis;
		//_thr = compThrust(recocmslist, thrAxis);
		
		// quantities based on momentum tensor
		ntpval[28] = evs.Sphericity();				// sphericity = 3/2 * (lam_2+lam_3)
		ntpval[29] = evs.Aplanarity();				// aplanarity = 3/2 * lam_3
		ntpval[30] = evs.Planarity();				// planarity  = lam_2 - lam_3
		ntpval[67] = evs.Circularity();				// circularity = 2 * min(lam_1, lam_2)/(lam_1 + lam_2)
		ntpval[31] = evs.Thrust();					// thrust
		
		//compFoxWolfMoms(recocmslist);
		
		ntpval[34] = evs.FoxWolfMomR(1);					// Fox Wolfram Moment R1 = H1/H0
		ntpval[35] = evs.FoxWolfMomR(2);					// Fox Wolfram Moment R2 = H2/H0
		ntpval[36] = evs.FoxWolfMomR(3);					// Fox Wolfram Moment R3 = H3/H0
		ntpval[68] = evs.FoxWolfMomR(4);					// Fox Wolfram Moment R4 = H4/H0
		ntpval[69] = evs.FoxWolfMomR(5);					// Fox Wolfram Moment R5 = H5/H0
		
		// count electrons
		for (j=0; j<eplus.size(); ++j) 
		{
			if (eplus[j].P4().P()>1.0) ntpval[47]+=1.0;
			if (eplus[j].P4().P()>2.0) ntpval[48]+=1.0;
		}
		for (j=0; j<eminus.size(); ++j) 
		{
			if (eminus[j].P4().P()>1.0) ntpval[47]+=1.0;
			if (eminus[j].P4().P()>2.0) ntpval[48]+=1.0;
		}
		
		// count muons
		for (j=0; j<muplus.size(); ++j) 
		{
			if (muplus[j].P4().P()>1.0) ntpval[49]+=1.0;
			if (muplus[j].P4().P()>2.0) ntpval[50]+=1.0;
		}
		for (j=0; j<muminus.size(); ++j) 
		{
			if (muminus[j].P4().P()>1.0) ntpval[49]+=1.0;
			if (muminus[j].P4().P()>2.0) ntpval[50]+=1.0;
		}
		
		// count pions
		for (j=0; j<piplus.size(); ++j) 
		{
			if (piplus[j].P4().P()>1.0) ntpval[51]+=1.0;
			if (piplus[j].P4().P()>2.0) ntpval[52]+=1.0;
		}
		for (j=0; j<piminus.size(); ++j) 
		{
			if (piminus[j].P4().P()>1.0) ntpval[51]+=1.0;
			if (piminus[j].P4().P()>2.0) ntpval[52]+=1.0;
		}
		
		// count kaons
		for (j=0; j<kplus.size(); ++j) 
		{
			if (kplus[j].P4().P()>1.0) ntpval[53]+=1.0;
			if (kplus[j].P4().P()>2.0) ntpval[54]+=1.0;
		}
		for (j=0; j<kminus.size(); ++j) 
		{
			if (kminus[j].P4().P()>1.0) ntpval[53]+=1.0;
			if (kminus[j].P4().P()>2.0) ntpval[54]+=1.0;
		}
		
		// count protons
		for (j=0; j<pplus.size(); ++j) 
		{
			if (pplus[j].P4().P()>1.0) ntpval[55]+=1.0;
			if (pplus[j].P4().P()>2.0) ntpval[56]+=1.0;
		}
		for (j=0; j<pminus.size(); ++j) 
		{
			if (pminus[j].P4().P()>1.0) ntpval[55]+=1.0;
			if (pminus[j].P4().P()>2.0) ntpval[56]+=1.0;
		}
		
		// count pi0
		for (j=0; j<LPi0.size(); ++j) 
		{
			if (LPi0[j].P4().P()>0.5) ntpval[86]+=1.0;
			if (LPi0[j].P4().P()>1.0) ntpval[88]+=1.0;
		}
		// count KS
		for (j=0; j<LKs.size(); ++j) 
		{
			if (LKs[j].P4().P()>0.5) ntpval[87]+=1.0;
			if (LKs[j].P4().P()>1.0) ntpval[89]+=1.0;
		}
		
		
		
		for (j=0;j<recolist.size();++j)
		{
			// *** candidate in lab system
			SimpleCand     cnd    = recolist[j];
			TLorentzVector lv     = cnd.P4();
			
			// *** candidate in cms system
			SimpleCand     cndcms = recocmslist[j];
			TLorentzVector lvcm   = cndcms.P4();
			
			// compute in LAB
			
			
			if (lvcm.Pt()>0.5) ntpval[20]+=1.0;
			if (lvcm.Pt()>1.0) ntpval[21]+=1.0;
			if (lvcm.Pt()>1.5) ntpval[22]+=1.0;
			if (lvcm.Pt()>2.0) ntpval[23]+=1.0;
			if (lvcm.Pt()>2.5) ntpval[24]+=1.0;
			if (lvcm.Pt()>3.0) ntpval[25]+=1.0;
			if (lvcm.Pt()>4.0) ntpval[26]+=1.0;
			if (lvcm.Pt()>5.0) ntpval[27]+=1.0;
			
        }

		ntpev->Fill(ntpval);

        mixevts=1;
        if (fRand.Rndm()<P_mix) mixevts=2;
    }

	
    ggg->cd();
    ntpev->Write();
    ggg->Close();

}