toy_core.C

Go to the documentation of this file.

#include "EventShape.h"
#include "SimpleNtp.C"
//#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 "TGaxis.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            = 6.0;      // *** higher edge of histo (will be modified acc. to resonance mass)

// *** some contants for neutrals
const double    fNeutEff   = 0.95;      // *** efficiency for neutrals
const double    fdE        = 0.05;      // *** neutral energy resolution
const double    fNdtht     = 0.003;     // *** neutrals tht resolution
const double    fNdphi     = 0.003;     // *** neutrals phi resolution

// number of sigmas selection for pi0, KS
int nsigpiks = 5;
// number of sigmas selection for channels
int nsig     = 8;

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

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

int             fDPdg      = 411;         // *** PDG code to match D
double                  fDMass     = 1.8693;      // *** center of mass cut window
double                  fDsig      = 0.045;        // *** sigma of peak
double                  fDwin           =0.25;

int             fD0Pdg     = 421;         // *** PDG code to match D0
double                  fD0Mass    = 1.8645;      // *** center of mass cut window
double                  fD0sig     = 0.055;        // *** sigma of peak
double                  fD0win          =0.25;

int             fJpsiPdg   = 443;         // *** PDG code to match J/psi
double                  fJMass     = 3.09687;        // *** center of mass cut window
double                  fJsig      = 0.11;        // *** sigma of peak
double                  fJwin           =0.7;
        
int             feePdg    = 88880;         // *** PDG code to match pbarpSystem
double                  feeMass   = 0;             // *** center of mass cut window, here = sqrt(s)
double                  feesigconst = 0.0024;        // *** const of linear width function
double                  feesigslope = 0.036;        // *** slope of linear width function
double                  feewin    = 0.;             // *** to be determined by function 

int             fLamPdg    = 3122;         // *** PDG code to match Lam
double                  fLamMass   = 1.115684;        // *** center of mass cut window
double                  fLamsig    = 0.0042;        // *** sigma of peak
double                  fLamwin         =0.1;

int             fDsPdg     = 431;         // *** PDG code to match Ds
double                  fDsMass    = 1.9686;        // *** center of mass cut window
double                  fDssig     = 0.04;        // *** sigma of peak
double                  fDswin          =0.2;
        
int             fPhiPdg     = 333;         // *** PDG code to match Phi
double                  fPhiMass    = 1.0195;        // *** center of mass cut window
double                  fPhisig     = 0.0065;        // *** sigma of peak
double                  fPhiwin         =0.2;
        
int             fEtacPdg     = 441;         // *** PDG code to match etac
double                  fEtacMass    = 2.9798;        // *** center of mass cut window
double                  fEtacsig     = 0.08;        // ***  sigma of peak
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.035;        // *** sigma of peak
double                  fLamcwin          =0.2;

// *** 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;
int npi, nk, npr, ne, nmu, nlep, evtno, globmode;

// *** 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;

// *** setup  style
void setStyle()
{
        gStyle->SetPadTopMargin(0.13);
        gStyle->SetPadBottomMargin(0.18);
        gStyle->SetPadLeftMargin(0.2);
        gStyle->SetPadRightMargin(0.015);

        gStyle->SetLabelSize(0.075,"X");
        gStyle->SetLabelSize(0.07,"Y");

        gStyle->SetStatY(0.87);
        gStyle->SetStatX(0.985);
        gStyle->SetStatW(0.4);
        gStyle->SetStatH(0.12);
        gStyle->SetOptStat("e"); // only entries
        
        gStyle->SetTitleH(0.1);
        gStyle->SetTitleX(0.08);
        gStyle->SetTitleOffset(1.1,"x");
        gStyle->SetTitleXSize(0.075);
}

// *** 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);
    }
}

// *** 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;
}

// *** is there one mass combination in the window?
bool eventAccepted(CandList &l,  double m =0, double w=0)
{
        bool accepted = false;
        unsigned int i=0;
        
        while (!accepted && i<l.size())
        //for (unsigned int i=0;i<l.size();++i) 
        {
                double mass = l[i++].P4().M();
                if (mass>m-w && mass<m+w) accepted =true;
        }
        
        return accepted;
}


// *** 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());
                        c.AddDau(&(l1[i]));c.AddDau(&(l2[j]));
                        c.SetPdg(matchPdg);
                        
                        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());
                                c.SetEvtId(ev1);
                        }
                        
                        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() );
                                c.AddDau(&(l1[i1]));c.AddDau(&(l2[i2]));c.AddDau(&(l3[i3]));
                                c.SetPdg(matchPdg);
                                
                                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());
                                        c.SetEvtId(ev1);
                                }
                                
                                out.push_back(c);
                        }
                }
        }
        
        return out.size();
}

int combine(CandList &l1, CandList &l2, CandList &l3, CandList &l4, CandList &out, int matchPdg=0)
{
        out.clear();
        
        bool same12 = (&l1 == &l2);
        bool same13 = (&l1 == &l3);
        bool same14 = (&l1 == &l4);
        
        bool same23 = (&l2 == &l3);
        bool same24 = (&l2 == &l4);
        
        bool same34 = (&l3 == &l4);
        
        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;
                                
                                for (unsigned int i4=0;i4<l4.size();++i4)
                                {
                                        if ((same14 && i4<i1) || (same24 && i4<i2) || (same34 && i4<i3)) continue;
                                        
                                        if ( l1[i1].Marker() & l4[i4].Marker() ) continue;
                                        if ( l2[i2].Marker() & l4[i4].Marker() ) continue;
                                        if ( l3[i3].Marker() & l4[i4].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() + l4[i4].P4() );
                                        c.SetCharge( l1[i1].Charge() + l2[i2].Charge() + l3[i3].Charge() + l4[i4].Charge() );
                                        c.SetNFS( l1[i1].NFS() + l2[i2].NFS() + l3[i3].NFS() + l4[i4].NFS() );
                                        c.SetMarker( l1[i1].Marker() | l2[i2].Marker() | l3[i3].Marker() | l4[i4].Marker() );
                                        c.AddDau(&(l1[i1]));c.AddDau(&(l2[i2]));c.AddDau(&(l3[i3]));c.AddDau(&(l4[i4]));
                                        c.SetPdg(matchPdg);
                                        
                                        int m1 = l1[i1].MotherIdx();
                                        int m2 = l2[i2].MotherIdx();
                                        int m3 = l3[i3].MotherIdx();
                                        int m4 = l4[i4].MotherIdx();
                                        
                                        int ev1 = l1[i1].EvtId();
                                        int ev2 = l2[i2].EvtId();
                                        int ev3 = l3[i3].EvtId();
                                        int ev4 = l4[i4].EvtId();
                                        
                                        if ( ev1==ev2 && ev2==ev3 && ev3==ev4 &&                        // *** particles are from same event
                                                m1>=0 && m1 == m2 && m2 == m3 && m3 == m4 &&    // *** 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()) &&                               // *** "
                                                (l4[i4].McPid()||l4[i4].Mct()) &&                               // *** "
                                                                4==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());
                                                c.SetEvtId(ev1);
                                        }
                                        
                                        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();
        
        unsigned int marker = 1;
        
        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);
                        c.SetMarker(marker);
                        c.SetPdg(211*c.Charge());
                        marker *=2;
                        reco.push_back(c);
                }
                else if (pid==22 && isDetected(fNeutEff))  // *** neutral
                {
                        smearMom(c,fdE, fNdtht, fNdphi);
                        c.SetMarker(marker);
                        c.SetPdg(22);
                        marker *=2;
                        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.139;

        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();
                        c.SetPdg(pdg);
                        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 configHisto(TH1* h)
{
        h->SetFillStyle(3001);
        h->SetFillColor(602);
        h->SetLineStyle(3);
        h->SetLineColor(602);
/*      TGaxis *ax = (TGaxis*)h->GetYaxis();
        ax->SetMaxDigits(2);*/
}

void qaCand(TString pre, SimpleCand &c, SimpleNtp *n)
{
        TLorentzVector l = c.P4();
        TVector3 boost = fIni.BoostVector();
        TLorentzVector lb = l;
        lb.Boost(-boost);
        
        n->Column(pre+"px",  l.Px());
        n->Column(pre+"py",  l.Py());
        n->Column(pre+"pz",  l.Pz());
        n->Column(pre+"p",   l.P());
        n->Column(pre+"e",   l.E());
        n->Column(pre+"m",   l.M());
        n->Column(pre+"pdg", c.Pdg());
        n->Column(pre+"tpdg", c.Pid());
        n->Column(pre+"ch",  c.Charge());
        
        n->Column(pre+"pt",  l.Pt());
        n->Column(pre+"tht", l.Theta());
        n->Column(pre+"phi", l.Phi());
        n->Column(pre+"pcm", lb.P());
        
        n->Column(pre+"mct", c.Mct()||c.McPid() );
}

// int getDau( int idx, SimpleCand &c)
// {
//      int num=0;
//      int curridx=-1;
//      unsigned int marker = c.Marker();
//      
//      while (curridx<idx)
//      {
//              while (!(marker&1) && num<32) {num++;marker>>=1;}
//              curridx++;
//      }
//      return num;
// }
// 
// int countDau(SimpleCand &c)
// {
//      int num=0;
//      unsigned int marker = c.Marker();
//      while (marker) {marker &= (marker-1); num++;}
//      
//      return num;
// }
// 
// int getDauIdx(SimpleCand &c, int* dauidx)
// {
//      int num=0;
//      int idx=0;
//      unsigned int marker = c.Marker();
// 
//      while (marker)
//      {
//              if (marker & 1) 
//              {
//                      dauidx[num]=idx;
//                      num++;
//              }
//              marker >>= 1;
//              idx++;
//      }
//      return num;
// }

void writeTuple(TString pre, CandList &l, EventShape &evsh, SimpleNtp *n, double mass, double win)
{
        int dauidx[32];
        for (int i=0;i<l.size();++i)
        {
                qaCand(pre,l[i], n);
                float tag = 0;
                float mct = l[i].Mct();
                
                if (fabs(l[i].P4().M()-mass)<win) tag=1.0;
                
                n->Column("tag",tag);
                n->Column("mct",mct);
        
                //int nDau = getDauIdx(l[i], dauidx);
                
                //cout <<"c("<<l[i].Marker()<<") "<<i<<" made from ";
                for (int j=0;j<l[i].NDau();++j)
                {
                        //cout <<dauidx[j]<<"  ";
                        //SimpleCand d = recolist[dauidx[j]];
                        SimpleCand d=*(l[i].Dau(j));
                        qaCand(TString::Format("%sd%d",pre.Data(),j),d, n);
                }
                //cout <<endl;
                
                n->Column("evt", evtno);
                n->Column("mode", globmode);
                
                // store eventshape
                n->Column("npart", evsh.NParticles());
                n->Column("nneut", evsh.NNeutral());
                n->Column("nchrg", evsh.NCharged());
                n->Column("sumpc", evsh.ChrgPSumCms());
                n->Column("sumpt", evsh.PtSumLab());
                n->Column("sumptc", evsh.ChrgPtSumLab());
                n->Column("pmaxl" ,evsh.PmaxLab());
                n->Column("pmax", evsh.PmaxCms());
                n->Column("pmin", evsh.PminCms());
                n->Column("ptmax",evsh.Ptmax());
                n->Column("thr",  evsh.Thrust());
                n->Column("pla",  evsh.Planarity());
                n->Column("apl",  evsh.Aplanarity());
                n->Column("sph",  evsh.Sphericity());
                n->Column("cir",  evsh.Circularity());
                n->Column("ne",   ne);
                n->Column("nmu",  nmu);
                n->Column("npi",  npi);
                n->Column("nk",   nk);
                n->Column("npr",  npr);
                n->Column("fw1", evsh.FoxWolfMomR(1));
                n->Column("fw2", evsh.FoxWolfMomR(2));
                n->Column("fw3", evsh.FoxWolfMomR(3));
                n->Column("fw4", evsh.FoxWolfMomR(4));
                n->Column("fw5", evsh.FoxWolfMomR(5));
                
                n->DumpData();
        }
}

int toy_core(TString fsig, int nev=2, double sqrts=3.77, bool simcut=true, bool evcut=false, double dp=0.05, 
                        double trkeff=95., double pideff=95., double misid =5., double P_mix=0.00, bool writentp=false, int mode =9999)
{
        setStyle();
        
        globmode= mode;
        
        fRand.SetSeed(314);
        
        init(pideff, misid);
        
        int i,j,k;
        
        makeIni4Vector(fIni, sqrts);
        
        fIni.Print();
        
        fPi0win = nsigpiks*fPi0sig;
        fKswin  = nsigpiks*fKssig;
        
        fPhiwin  = nsig*fPhisig;        
        fDwin    = nsig*fDsig;
        fD0win   = nsig*fD0sig;
        fDswin   = nsig*fDssig;
        fJwin    = nsig*fJsig;  
        fEtacwin = nsig*fEtacsig;       
        feewin   = nsig*(feesigconst+sqrts*feesigslope);        
        fLamcwin = nsig*fLamcsig;       
        fLamwin  = nsig*fLamsig;        
        
        //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
        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;
        
        TCanvas *c1=new TCanvas("c1","c1",10,10,600,300);
        c1->Divide(2,1);
        TCanvas *c2=new TCanvas("c2","c2",20,20,1800,800);
        c2->Divide(5,2,0.005,0.01);
        
/*      TCanvas *c3=new TCanvas("c3","c3",1010,10,600,600);
        c3->Divide(2,2);*/
        
//      TCanvas *c4=new TCanvas("c4","c4",1220,20,500,500);
        
        double min=0, max=6.0, min2=0.5;
        
        // *** create some histograms
        TH1F *h_mom= new TH1F("h_mom","Momenta",200,0,5);
        TH1F *h_momc= new TH1F("h_momc","Momenta charged",200,0,5);
        TH1F *h_momn= new TH1F("h_momn","Momenta neutral",200,0,5);
        TH1F *h_mult= new TH1F("h_mult","Multiplicity",50,-0.5,49.5);
        TH1F *h_multrec= new TH1F("h_multrec","Multiplicity reco",50,-0.5,49.5);
        TH1F *h_multc= new TH1F("h_multc","Multiplicity charged",30,-0.5,29.5);
        TH1F *h_multn= new TH1F("h_multn","Multiplicity neutral",50,-0.5,49.5);
        TH1F *h_fw1= new TH1F("h_fw1","Fox Wolfram R1",300,-1.1,1.1);
        TH1F *h_pmax= new TH1F("h_pmax","pmax cms",300,0,5);
        
        TH1F *h_mult_k = new TH1F("h_mult_k","Multiplicity Kaons",10,-0.5,9.5);
        TH1F *h_mult_lep = new TH1F("h_mult_lep","Multiplicity Leptons",10,-0.5,9.5);
        TH1F *h_sum_klep = new TH1F("h_sum_lep","Sum Kaons + Leptons",15,-0.5,14.5);
        TH2F *h_mult_klep=new TH2F("h_mult_klep","#Kaons vs. #Leptons",10,-0.5,9.5, 10, -0.5, 9.5);
        h_mult_klep->SetYTitle("#Leptons");
        h_mult_klep->SetXTitle("#Kaons");
        
        // ****** pi0
        TH1F *h_pi0=new TH1F("h_pi0","#pi^{0} #rightarrow #gamma#gamma",200,min,0.4);
        TH1F *h_pi0sel=new TH1F("h_pi0sel","Pi0",200,min,0.4);
        TH1F *h_pi0sig=new TH1F("h_pi0sig","Pi0",200,min,0.4);
        h_pi0sig->SetLineColor(2);
        configHisto(h_pi0sel);
        
        // ****** Ks
        TH1F *h_ks=new TH1F("h_ks","K_{S}^{0} #rightarrow #pi^{+}#pi^{-}",200,0.2,0.8);
        TH1F *h_kssel=new TH1F("h_kssel","Pi0",200,0.2,0.8);
        TH1F *h_kssig=new TH1F("h_kssig","Pi0",200,0.2,0.8);
        h_kssig->SetLineColor(2);
        configHisto(h_kssel);
        
        // ****** Jpsi -> ee
        TH1F *h_jpsi=new TH1F("h_jpsi","J/#psi #rightarrow e^{+}e^{-}",nbins,min,max);
        TH1F *h_jpsisel=new TH1F("h_jpsisel","J/psi",nbins,min,max);
        TH1F *h_jpsisig=new TH1F("h_jpsisig","J/psi",nbins,min,max);
        h_jpsisig->SetLineColor(2);
        configHisto(h_jpsisel);
        
        // ****** Jpsi -> mu mu
        TH1F *h_jpsi2=new TH1F("h_jpsi2","J/#psi #rightarrow #mu^{+}#mu^{-}",nbins,min,max);
        TH1F *h_jpsi2sel=new TH1F("h_jpsi2sel","J/psi",nbins,min,max);
        TH1F *h_jpsi2sig=new TH1F("h_jpsi2sig","J/psi",nbins,min,max);
        h_jpsi2sig->SetLineColor(2);
        configHisto(h_jpsi2sel);
        
        // ****** eta_c -> Ks K pi
        TH1F *h_etac=new TH1F("h_etac","#eta_{c} #rightarrow K_{S} K^{+} #pi^{-}",nbins,min,max);
        TH1F *h_etacsel=new TH1F("h_etacsel","etac",nbins,min,max);
        TH1F *h_etacsig=new TH1F("h_etacsig","etac",nbins,min,max);
        h_etacsig->SetLineColor(2);
        configHisto(h_etacsel);
        
        // ****** D0 -> K- pi+
        TH1F *h_d0=new TH1F("h_d0","D^{0} #rightarrow K^{-}#pi^{+}",nbins,min,max);
        TH1F *h_d0sel=new TH1F("h_d0sel","D0",nbins,min,max);
        TH1F *h_d0sig=new TH1F("h_d0sig","D0",nbins,min,max);
        h_d0sig->SetLineColor(2);
        configHisto(h_d0sel);
                
        // ****** D+ -> K- pi+ pi+
        TH1F *h_dpm=new TH1F("h_dpm","D^{+} #rightarrow K^{-}#pi^{+}#pi^{+}",nbins,min,max);
        TH1F *h_dpmsel=new TH1F("h_dpmsel","D#pm",nbins,min,max);
        TH1F *h_dpmsig=new TH1F("h_dpmsig","D#pm",nbins,min,max);
        h_dpmsig->SetLineColor(2);
        configHisto(h_dpmsel);
                
        // ****** Ds+ -> K+ K- pi+
        TH1F *h_ds=new TH1F("h_ds","D_{s}^{+} #rightarrow K^{+}K^{-}#pi^{+}",nbins,min,max);
        TH1F *h_dssel=new TH1F("h_dssel","Ds",nbins,min,max);
        TH1F *h_dssig=new TH1F("h_dssig","Ds",nbins,min,max);
        h_dssig->SetLineColor(2);
        configHisto(h_dssel);
        
        // ****** phi -> K+ K-
        TH1F *h_phi=new TH1F("h_phi","#phi #rightarrow K^{+}K^{-}",nbins,min2,max/2);
        TH1F *h_phisel=new TH1F("h_phisel","Phi",nbins,min2,max/2);
        TH1F *h_phisig=new TH1F("h_phisig","Phi",nbins,min2,max/2);
        h_phisig->SetLineColor(2);
        configHisto(h_phisel);
        
        // ****** Lambda_c -> p K- pi+
        TH1F *h_lamc=new TH1F("h_lamc","#Lambda_{c} #rightarrow pK^{-}#pi^{+}",nbins,min,max);
        TH1F *h_lamcsel=new TH1F("h_lamcsel","Lam_c",nbins,min,max);
        TH1F *h_lamcsig=new TH1F("h_lamcsig","Lam_c",nbins,min,max);
        h_lamcsig->SetLineColor(2);
        configHisto(h_lamcsel);
        
        // ****** Lambda -> p pi-
        TH1F *h_lam=new TH1F("h_lam","#Lambda #rightarrow p#pi^{-}",nbins,min2,max/2);
        TH1F *h_lamsel=new TH1F("h_lamsel","Lam",nbins,min2,max/2);
        TH1F *h_lamsig=new TH1F("h_lamsig","Lam",nbins,min2,max/2);
        h_lamsig->SetLineColor(2);
        configHisto(h_lamsel);
        
        // ****** pp -> e+ e-
        TH1F *h_ee=new TH1F("h_ee","p#bar{p} #rightarrow e^{+}e^{-}",nbins,min,max*1.2);
        TH1F *h_eesel=new TH1F("h_eesel","ee",nbins,min,max*1.2);
        TH1F *h_eesig=new TH1F("h_eesig","ee",nbins,min,max*1.2);
        h_eesig->SetLineColor(2);
        configHisto(h_eesel);
        
        
        h_jpsi->SetXTitle("m(e^{+}e^{-}) [GeV/c^{2}]");
        h_jpsi2->SetXTitle("m(#mu^{+}#mu^{-}) [GeV/c^{2}]");
        h_etac->SetXTitle("m(K_{S} K^{+}#pi^{-}) [GeV/c^{2}]");

        h_d0->SetXTitle("m(K^{-}#pi^{+}) [GeV/c^{2}]");

        h_dpm->SetXTitle("m(K^{-}#pi^{+}#pi^{+}) [GeV/c^{2}]");

        h_ds->SetXTitle("m(K^{+}K^{-}#pi^{+}) [GeV/c^{2}]");

        h_phi->SetXTitle("m(K^{+}K^{-}) [GeV/c^{2}]");
        h_lamc->SetXTitle("m(pK^{-}#pi^{+}) [GeV/c^{2}]");
        h_lam->SetXTitle("m(p#pi^{-}) [GeV/c^{2}]");
        h_ee->SetXTitle("m(e^{+}e^{-}) [GeV/c^{2}]");
        
        // *** 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);     
        }
        
        TFile *fn=0;
        
        TString filename(fsig(fsig.Index("/")+1,fsig.Length())+"_ntp.root");
        
        if (writentp) fn=new TFile(filename, "RECREATE");
        
        SimpleNtp *nphi         = new SimpleNtp("nphi","phi -> K+K-");
        SimpleNtp *nd0          = new SimpleNtp("nd0","D0 -> K+ pi+");
        SimpleNtp *ndpm         = new SimpleNtp("ndpm","D+ -> K- pi+ pi+");
        SimpleNtp *nds          = new SimpleNtp("nds","Ds -> K+ K- pi+");
        SimpleNtp *njpsi1       = new SimpleNtp("njpsi1","J/psi -> e+ e-");
        SimpleNtp *njpsi2       = new SimpleNtp("njpsi2","J/psi -> mu+ mu-");
        SimpleNtp *netac        = new SimpleNtp("netac","eta_c -> KS (pi+pi-) K- pi+");
        SimpleNtp *nee          = new SimpleNtp("nee","ppbar -> e+ e-");
        SimpleNtp *nlam         = new SimpleNtp("nlam","Lambda -> p pi-");
        SimpleNtp *nlamc        = new SimpleNtp("nlamc","Lambda_c -> p K- pi+");                
                
        int combCnt=0;  // how many combinatione are within our mass criterion
        int evCnt=0;    // how many _events_ (not combinations!) are within our mass criterion
        int evMultCut=0;
        
        // *** the event loop
        
        int mixevts=1;
        int currmix=0;
        
        int cntInLam=0;
        
        // how many events are tagged by individual algo?
        bool channelselect[10];
        double channelcount[10];
        TString channelname[10]={"Phi(KK)","pp->ee" ,"Lam(ppi)","etac(KS K pi)","J/psi(ee)","J/psi(mumu)","D+-(K2pi)","D0(Kpi)","Ds(KKpi)","Lamc(pKpi)"};
        
        for (i=0;i<10;++i) channelcount[i]=0.;
        
        if (fRand.Rndm()<P_mix) mixevts=2;
        
        for (i=0;i<nev;++i)
        {       
                evtno = 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)
                {
                        int npi = 0;
                        int nstab = 0;
                        
                        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;
                                if (pdg==211) npi++;
                                if (isFinalState(pdg)) nstab++; 
                        }
                        //if (npi!=4 || nstab!=npi ) good = false; 
                }
                
                if (!good) continue;
                
                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;
                
                h_mult->Fill(mcFSlist.size());
                
                // *** if mass cut is negative (= number of sigmas), compute absolute window size
                //if (mwin<0) mwidth=-mwin*(dp*mwsig+mwoff);
                
                // *** prepare the toy reco: smear momentum and apply track efficiency (both flat at the moment)
                makeRecoCands(mcFSlist, recolist, trkeff/100., dp );
                
                // *** the boost of the system
                TVector3 boost = fIni.BoostVector();
                boostList(recolist, recocmslist, boost);
                
                // *** Make Event Selection
                EventShape evsh(recolist, fIni);
/*              if (evcut)
                {
                        if (evshape.ChrgPtSumLab()<1.5) continue;
                }*/
                
                h_multrec->Fill(recolist.size());
                //cout <<mcFSlist.size()<<"  "<<recolist.size()<<endl;
                
                // *** prepare all pid list according to the pidtable probabilities
                makePidSelection(recolist);
                
                npi = piplus.size()+piminus.size();
                nk = kplus.size()+kminus.size();
                npr = pplus.size()+pminus.size();
                ne = eplus.size()+eminus.size();
                nmu = muplus.size()+muminus.size();
                nlep = ne+nmu;
                
/*              int nk10 = 0;  // # kaons with p>1.0
                for (j=0;j<kplus.size();++j) if (kplus[j].P4().P()>1.0) nk10++;
                for (j=0;j<kminus.size();++j) if (kminus[j].P4().P()>1.0) nk10++;
                
                int npr10 = 0;  // # protons with p>1.0
                for (j=0;j<pplus.size();++j) if (pplus[j].P4().P()>1.0) npr10++;
                for (j=0;j<pminus.size();++j) if (pminus[j].P4().P()>1.0) npr10++;*/
                
                h_mult_k->Fill(nk);
                h_mult_lep->Fill(nlep);
                h_mult_klep->Fill(nk, nlep);
                h_sum_klep->Fill(nk + nlep);
                
                                                                        
                // *** the lists we need for combinatorics
                CandList comb1, comb2, comb3;

                CandList LPhi, LPi0, LKs;                               // phi -> KK , pi0 -> gg (mass cut), KS -> pi+ pi- (mass cut)
                CandList LJpsiE, LJpsiMu;                       // J/psi -> ee/mu mu
                CandList LEtac, LEtac2;                 // etac
                CandList LDp, LDm, LDp2, LDm2, LDp3, LDm3, LDp4, LDm4;  // D+- -> K2pi, K2pipi0, KS pi pi0, KS pi pi pi
                CandList LDsp, LDsm, LDsp2, LDsm2;              // Ds -> KKpi, KKpipi0                          
                CandList LD0, LD02, LD0b, LD02b, LD03, LD03b;   // D0 -> Kpi, Kpipi0, K3pi
                CandList LLamcp, LLamcm;                                // Lc -> pKpi
                CandList LLam, LLamb;                                   // Lam -> ppi
                CandList Lee;                                   // Lam -> ppi
                
                int ncomb1=0, ncomb2=0;
                bool selected=false;
                
                for (j=0;j<13;++j) channelselect[j]=false;
                                
                // Apply cut on kaon/lepton multiplicity
/*              if ( multcut>0 && (N_K+N_lep<multcut) )
                {
                        mixevts=1;
                        if (fRand.Rndm()<P_mix) mixevts=2;
                        continue;
                }*/
                
                evMultCut++;

                bool accepted = false;
                
                //
                //  ********* pi0 / KS
                //
                // *** 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);
                

                // **** cache some eventshape vars
                int npart      = evsh.NParticles();
                int nchrg      = evsh.NCharged();
                int nneut      = evsh.NNeutral();
                int npi0       = LPi0.size();                                   // #pi0 candidates
                int nks0       = LKs.size();                                    // #KS0 candidates
                int npi005     = 0.0;
                int     npi010     = 0.0;
                int nks005     = 0.0;
                int nks010     = 0.0;
                int np05       = evsh.MultPminCms(0.5);
                int np10       = evsh.MultPminCms(1.0);

                double pmaxl   = evsh.PmaxLab();
                double pmax    = evsh.PmaxCms();
                double pmin    = evsh.PminCms();
                double ptmax   = evsh.Ptmax();
                
                double sumpc   = evsh.ChrgPSumCms();
                double sumpc05 = evsh.SumChrgPminCms(0.5);
                double sumpt   = evsh.PtSumLab();
                double sumptc  = evsh.ChrgPtSumLab();
                double sumptl  = sumpt;
                double sumen   = evsh.NeutESumCms();
                
                double fw1     = evsh.FoxWolfMomR(1);
                double fw2     = evsh.FoxWolfMomR(2);
                double fw3     = evsh.FoxWolfMomR(3);
                double fw4     = evsh.FoxWolfMomR(4);
                double fw5     = evsh.FoxWolfMomR(5);
                
                // count pi0 with threshold
                for (j=0; j<LPi0.size(); ++j) 
                {
                        if (LPi0[j].P4().P()>0.5) npi005++;
                        if (LPi0[j].P4().P()>1.0) npi010++;
                }
                // count KS with threshold
                for (j=0; j<LKs.size(); ++j) 
                {
                        if (LKs[j].P4().P()>0.5) nks005++;
                        if (LKs[j].P4().P()>1.0) nks010++;
                }
                
                h_fw1->Fill(fw2);
                h_pmax->Fill(pmax);

                
                // **** find out, whether event was accepted
                
                // ********* PHI -> K+ K-
                if ( !evcut 
                     || (sqrts==2.4  &&   nk>1 && npart<7 && pmax<0.648   )
                         || (sqrts==3.77 &&   nk>1 && fw2>0.5456 )
                         || (sqrts==4.28 &&   nk>1 && 1 )
                         || (sqrts==5.0  &&   nk>1 && 1 )
                         || (sqrts==5.5  &&   nk>1 && 1 )
                        )
                {
                  // *** do combinatorics (phi -> K K) 
                  combine (kplus, kminus, LPhi, fPhiPdg); 
                  
                  accepted |= eventAccepted(LPhi,fPhiMass, fPhiwin);
                  
                  writeTuple("phi",LPhi,evsh, nphi,fPhiMass, fPhiwin);
                }


                // ********* LAMBDA -> p pi-
                if ( !evcut 
                     || (sqrts==2.4  &&   npart>3 && npr>0 && npi>0 && fw1>0.132 && fw4>0.2024 )
                         || (sqrts==3.77 &&   npart>3 && pmax>1.056 &&     fw1>0.0   && fw2>0.6248 && fw5>0.22  )
                         || (sqrts==4.28 &&   npart>3 && 1 )
                         || (sqrts==5.0  &&   npart>3 && 1 )
                         || (sqrts==5.5  &&   npart>3 && 1 )
                        )
                {
                  // *** do combinatorics (Lam -> p pi) 
                  combine (pplus, piminus, LLam, fLamPdg); 
                  combine (pminus, piplus, LLamb, -fLamPdg);            
                  
                  accepted |= eventAccepted(LLam, fLamMass, fLamwin);
                  accepted |= eventAccepted(LLamb, fLamMass, fLamwin);
                  
                  writeTuple("lam",LLam,evsh, nlam, fLamMass, fLamwin);
                  writeTuple("lam",LLamb,evsh, nlam, fLamMass, fLamwin);
                }

                // ********* pp -> e+ e-
                if ( !evcut 
                        )
                {
                  // *** do combinatorics (e+ e-) 
                  combine(eminus, eplus, Lee, 88880);
                  
                  accepted |= eventAccepted(Lee,  fIni.M(), feewin);
                  
                  writeTuple("ee",Lee,evsh, nee,  fIni.M(), feewin);
                }

                // ********* ETAC
                if ( !evcut 
/*                       || (sqrts==3.77 &&  np10>1 && npi<4 && pmin<0.3 && sumpc>3.296 )
                         || (sqrts==4.28 &&  1 )
                         || (sqrts==5.0  &&  1 )
                         || (sqrts==5.5  &&  1 )*/
                        )
                {
                  // *** do combinatorics (ETAC) 
                  combine(LKs, kplus, piminus, LEtac, fEtacPdg);
                  combine(LKs, kminus, piplus, LEtac2, fEtacPdg);
                  
                  accepted |= eventAccepted(LEtac,  fEtacMass, fEtacwin);
                  accepted |= eventAccepted(LEtac2,  fEtacMass, fEtacwin);
                  
                  writeTuple("etac",LEtac,evsh, netac,  fEtacMass, fEtacwin);
                  writeTuple("etac",LEtac2,evsh, netac,  fEtacMass, fEtacwin);
                }

                // ********* JPSI -> e+e-
                if ( !evcut 
                         || (sqrts==3.77 &&  np10>1 && npi<4 && pmin<0.3 && sumpc>3.296 )
                         || (sqrts==4.28 &&  1 )
                         || (sqrts==5.0  &&  1 )
                         || (sqrts==5.5  &&  1 )
                        )
                {
                  // *** do combinatorics (J/psi -> ee) 
                  combine(eminus, eplus, LJpsiE, fJpsiPdg);
                  
                  accepted |= eventAccepted(LJpsiE,  fJMass, fJwin);

                  writeTuple("jpsi1",LJpsiE,evsh, njpsi1,  fJMass, fJwin);
                }
                
                // ********* JPSI -> mu+mu-
                if ( !evcut 
                         || (sqrts==3.77 &&  np10>1 && npi<4 && pmin<0.3 && sumpc>3.296 )
                         || (sqrts==4.28 &&  1 )
                         || (sqrts==5.0  &&  1 )
                         || (sqrts==5.5  &&  1 )
                        )
                {
                  // *** do combinatorics (J/psi -> mu mu) 
                  combine(muminus, muplus, LJpsiMu, fJpsiPdg);
                  
                  accepted |= eventAccepted(LJpsiMu,  fJMass, fJwin);

                  writeTuple("jpsi2",LJpsiMu,evsh, njpsi2,  fJMass, fJwin);
                }
                
                
                // ********* D+ -> K- pi+ pi+
                if ( !evcut 
                         || (sqrts==3.77 &&   nk>0 && npart<12 && pmax<0.936 && sumpc>1.76 )
                         || (sqrts==4.28 &&   nk>0 && 1 )
                         || (sqrts==5.0  &&   nk>0 && 1 )
                         || (sqrts==5.5  &&   nk>0 && 1 )
                        )
                {
                  // *** do combinatorics (D+- -> K pi pi) 
                  combine (piplus,  piplus,  kminus, LDp, fDPdg); 
                  combine (piminus, piminus, kplus,  LDm, -fDPdg); 
                  
                  accepted |= eventAccepted(LDp,fDMass, fDwin);
                  accepted |= eventAccepted(LDm,fDMass, fDwin);
                  
                  writeTuple("d",LDp,evsh, ndpm,fDMass, fDwin);
                  writeTuple("d",LDm,evsh, ndpm,fDMass, fDwin);
                }
                
                
                // ********* D0 -> K- pi+
                if ( !evcut 
                         || (sqrts==3.77 &&   nk>0 && sumpc05>3.18 && pmax>0.768 && pmax<1.056 && fw3<0.1848 )
                         || (sqrts==4.28 &&   nk>0 && 1 )
                         || (sqrts==5.0  &&   nk>0 && 1 )
                         || (sqrts==5.5  &&   nk>0 && 1 )
                        )
                {
                  // *** do combinatorics (D0 -> K pi) 
                  combine (kplus, piminus, LD0, -fD0Pdg); 
                  combine (kminus, piplus, LD0b, fD0Pdg); 

                  accepted |= eventAccepted(LD0, fD0Mass, fD0win);
                  accepted |= eventAccepted(LD0b, fD0Mass, fD0win);
                  
                  writeTuple("d0",LD0,evsh, nd0, fD0Mass, fD0win);
                  writeTuple("d0",LD0b,evsh, nd0, fD0Mass, fD0win);
                }       


                // ********* DS -> K+ K- pi+
                if ( !evcut 
                         || (sqrts==4.28 &&   nk>1 && 1 )
                         || (sqrts==5.0  &&   nk>1 && 1 )
                         || (sqrts==5.5  &&   nk>1 && 1 )
                        )
                {
                  // *** do combinatorics (Ds+- -> K K pi) 
                  combine (kplus, kminus, piplus, LDsp, fDsPdg); 
                  combine (kplus, kminus, piminus, LDsm, -fDsPdg); 

                  accepted |= eventAccepted(LDsp, fDsMass, fDswin);
                  accepted |= eventAccepted(LDsm, fDsMass, fDswin);
                  
                  writeTuple("ds",LDsp,evsh, nds, fDsMass, fDswin);
                  writeTuple("ds",LDsm,evsh, nds, fDsMass, fDswin);
                }


                // ********* LAMBDA_C -> p K- pi+
                if ( !evcut 
                         || (sqrts==5.0  &&  1 )
                         || (sqrts==5.5  &&  1 )
                        )
                {
                  // *** do combinatorics (Lamc) 
                  combine (pplus, kminus, piplus, LLamcp, fLamcPdg); 
                  combine (pminus, kplus, piminus, LLamcm, -fLamcPdg);          
                
                  accepted |= eventAccepted(LLamcp, fLamcMass, fLamcwin);
                  accepted |= eventAccepted(LLamcm, fLamcMass, fLamcwin);
                  
                  writeTuple("lamc",LLamcp,evsh, nlamc, fLamcMass, fLamcwin);
                  writeTuple("lamc",LLamcm,evsh, nlamc, fLamcMass, fLamcwin);
                }

                
                // ********* NOW do counting and drawing, ONLY for events which are accepted by any selector                    
            // {"Phi(KK)", "ee", "Lam(ppi)", "etac", "J/psi(ee)","J/psi(mumu)","D+-(K2pi)", "D0(Kpi)","Ds(KKpi)","Lamc(pKpi)"};
        //     0        1        2         3         4             5           6           7           8          9        
                if (!simcut || accepted)
                {
                        // **** Phi
                        if (sqrts>2.2)
                        {
                          ncomb1 = fillHisto(LPhi, h_phi, 0, fPhiMass, fPhiwin, h_phisel, 0, h_phisig  ); 
                          combCnt += ncomb1;
                          if (ncomb1>0) {selected=true; channelselect[0]=true;}
                        }
                        
                        // **** ee
                        if (sqrts>2.2)
                        {
                          ncomb1 = fillHisto(Lee, h_ee, 0, fIni.M(), feewin, h_eesel, 0, h_eesig  ); 
                          combCnt += ncomb1;
                          if (ncomb1>0) {selected=true; channelselect[1]=true;}
                        }

                        // **** Lam
                        if (sqrts>2.2)
                        {
                          ncomb1 = fillHisto(LLam, h_lam, 0, fLamMass, fLamwin, h_lamsel, 0, h_lamsig ); 
                          ncomb2 = fillHisto(LLamb, h_lam, 0, fLamMass, fLamwin, h_lamsel, 0, h_lamsig ); 
                          combCnt += ncomb1 + ncomb2;
                          if (ncomb1>0 || ncomb2>0) {selected=true; channelselect[2]=true;}
                        }
                        
                        // **** etac
                        if (sqrts>3.1)
                        {
                          ncomb1 = fillHisto(LEtac, h_etac, 0, fEtacMass, fEtacwin, h_etacsel, 0, h_etacsig ); 
                          ncomb2 = fillHisto(LEtac2, h_etac, 0, fEtacMass, fEtacwin, h_etacsel, 0, h_etacsig ); 
                          combCnt += ncomb1 + ncomb2;
                          if (ncomb1>0 || ncomb2>0 ) {selected=true; channelselect[3]=true;}
                        }
                        
                        // **** J/psi ee
                        if (sqrts>3.1)
                        {
                          ncomb1 = fillHisto(LJpsiE, h_jpsi, 0, fJMass, fJwin, h_jpsisel, 0, h_jpsisig ); 
                          combCnt += ncomb1;
                          if (ncomb1>0 ) {selected=true; channelselect[4]=true;}
                        }
                        
                        // **** J/psi mumu
                        if (sqrts>3.1)
                        {
                          ncomb1 = fillHisto(LJpsiMu, h_jpsi2, 0, fJMass, fJwin, h_jpsi2sel, 0, h_jpsi2sig ); 
                          combCnt += ncomb1 ;
                          if (ncomb1>0 ) {selected=true; channelselect[5]=true;}
                        }
                        
                        // **** Dp
                        if (sqrts>3.75)
                        {
                          ncomb1 = fillHisto(LDp, h_dpm, 0, fDMass, fDwin, h_dpmsel, 0, h_dpmsig ); 
                          ncomb2 = fillHisto(LDm, h_dpm, 0, fDMass, fDwin, h_dpmsel, 0, h_dpmsig );     
                          combCnt += ncomb1 +ncomb2;
                          if (ncomb1>0 || ncomb2>0) {selected=true; channelselect[6]=true;}               
                        }

                        // **** D0 
                        if (sqrts>3.75)
                        {
                          ncomb1 = fillHisto(LD0, h_d0, 0, fD0Mass, fD0win, h_d0sel, 0, h_d0sig ); 
                          ncomb2 = fillHisto(LD0b, h_d0, 0, fD0Mass, fD0win, h_d0sel, 0, h_d0sig );     
                          combCnt += ncomb1 +ncomb2;
                          if (ncomb1>0 || ncomb2>0) {selected=true; channelselect[7]=true;}       
                        }

                        // **** Ds
                        if (sqrts>4.)
                        {
                          ncomb1 = fillHisto(LDsp, h_ds, 0, fDsMass, fDswin, h_dssel, 0, h_dssig ); 
                          ncomb2 = fillHisto(LDsm, h_ds, 0, fDsMass, fDswin, h_dssel, 0, h_dssig ); 
                          combCnt += ncomb1 +ncomb2;
                          if (ncomb1>0 || ncomb2>0) {selected=true; channelselect[8]=true;}
                        }
                
                        // **** Lamc
                        if (sqrts>4.57)
                        {
                          ncomb1 = fillHisto(LLamcp, h_lamc, 0, fLamcMass, fLamcwin, h_lamcsel, 0, h_lamcsig ); 
                          ncomb2 = fillHisto(LLamcm, h_lamc, 0, fLamcMass, fLamcwin, h_lamcsel, 0, h_lamcsig ); 
                          combCnt += ncomb1 + ncomb2;
                          if (ncomb1>0 || ncomb2>0) {selected=true; channelselect[9]=true;}
                        }
                }
            // ****************
            
                if (selected) evCnt+=mixevts; // *** if any combination is in window count event
                
                for (j=0;j<10;++j) if (channelselect[j]) channelcount[j]+=mixevts;
                
                //hmix->Fill(mixevts);
                // new round of mixing
                mixevts=1;
                if (fRand.Rndm()<P_mix) mixevts=2;
                
                //if (i<3) printAllLists();
        }
        //cout <<"cntInLam "<<cntInLam<<endl;
        cout<< "#Combinations:"<<combCnt<<"  #Events:"<<evCnt<<"  ";
        cout <<"Efficiency   (ev):"<<(double)evCnt/(double)nev<<endl;
        cout <<"Mult Cut Eff (ev):"<<(double)evMultCut/(double)nev<<endl;
        
        for (i=0;i<10;++i) 
        {
                channelcount[i]/=nev;
                cout <<channelname[i]<<":"<<channelcount[i]<<endl;
        }
        cout <<endl;
        
        TLatex latex;
        latex.SetTextSize(0.08);
        char tmp[200];
        double offset = 0.1;
        double pos=3.3, pos2=1.85;
        double hfac = 0.75;
        // *** plot all stuff and write out info
        c1->cd(1); h_pi0->Draw(); h_pi0sel->Draw("same"); h_pi0sig->Draw("same");
        c1->cd(2); h_ks->Draw(); h_kssel->Draw("same"); h_kssig->Draw("same");

        c2->cd(1); h_phi->Draw(); h_phisel->Draw("same");   h_phisig->Draw("same");
        sprintf(tmp,"#epsilon = %4.1f%%",channelcount[0]*100.);
        latex.DrawLatex(pos2,h_phi->GetMaximum()*hfac+offset,tmp);
        
        c2->cd(2); h_ee->Draw(); h_eesel->Draw("same"); h_eesig->Draw("same");
        sprintf(tmp,"#epsilon = %4.1f%%",channelcount[1]*100.);
        latex.DrawLatex(pos*1.2,h_ee->GetMaximum()*hfac+offset,tmp);
        
        c2->cd(3); h_lam->Draw(); h_lamsel->Draw("same"); h_lamsig->Draw("same");
        sprintf(tmp,"#epsilon = %4.1f%%",channelcount[2]*100.);
        latex.DrawLatex(pos2,h_lam->GetMaximum()*hfac+offset,tmp);
        
        c2->cd(4); h_etac->Draw(); h_etacsel->Draw("same"); h_etacsig->Draw("same");
        sprintf(tmp,"#epsilon = %4.1f%%",channelcount[3]*100.);
        latex.DrawLatex(pos,h_etac->GetMaximum()*hfac+offset,tmp);
        
        c2->cd(5); h_jpsi->Draw(); h_jpsisel->Draw("same"); h_jpsisig->Draw("same");
        sprintf(tmp,"#epsilon = %4.1f%%",channelcount[4]*100.);
        latex.DrawLatex(pos,h_jpsi->GetMaximum()*hfac+offset,tmp);
        
        c2->cd(6); h_jpsi2->Draw(); h_jpsi2sel->Draw("same"); h_jpsi2sig->Draw("same");
        sprintf(tmp,"#epsilon = %4.1f%%",channelcount[5]*100.);
        latex.DrawLatex(pos,h_jpsi2->GetMaximum()*hfac+offset,tmp);
        
        c2->cd(7); h_dpm->Draw();  h_dpmsel->Draw("same");  h_dpmsig->Draw("same");
        sprintf(tmp,"#epsilon = %4.1f%%",channelcount[6]*100.);
        latex.DrawLatex(pos,h_dpm->GetMaximum()*hfac+offset,tmp);
        
        c2->cd(8); h_d0->Draw();   h_d0sel->Draw("same");   h_d0sig->Draw("same");
        sprintf(tmp,"#epsilon = %4.1f%%",channelcount[7]*100.);
        latex.DrawLatex(pos,h_d0->GetMaximum()*hfac+offset,tmp);

        c2->cd(9); h_ds->Draw();   h_dssel->Draw("same");   h_dssig->Draw("same");
        sprintf(tmp,"#epsilon = %4.1f%%",channelcount[8]*100.);
        latex.DrawLatex(pos,h_ds->GetMaximum()*hfac+offset,tmp);
        
        c2->cd(10); h_lamc->Draw(); h_lamcsel->Draw("same"); h_lamcsig->Draw("same");
        sprintf(tmp,"#epsilon = %4.1f%%",channelcount[9]*100.);
        latex.DrawLatex(pos,h_lamc->GetMaximum()*hfac+offset,tmp);
        
        // total efficiency in first plot (phi)
        c2->cd(5);
        double posy = h_jpsi->GetMaximum()*1.12;
        if (posy=0) posy=1.12;
        sprintf(tmp,"(#epsilon_{t} = %4.1f%%)",double(evCnt)/double(nev)*100.);
        latex.SetTextColor(kMagenta+2);
        latex.DrawLatex(3.0,h_jpsi->GetMaximum()*1.12,tmp);

        
        //sprintf(tmp,"#epsilon = %4.1f%%",channelcount[6]*100.);
        //latex.DrawLatex(2.5,h_pi0->GetMaximum()*0.95,tmp);
        
/*      c2->cd();
        TPad *peff=new TPad("peff","peff",0.45,0.47,0.55,0.53);
        peff->Draw();
        peff->cd();
        sprintf(tmp,"#epsilon_{tot} = %4.1f%%",double(evCnt)/double(nev)*100.);
        latex.SetTextSize(0.5);
        latex.SetTextColor(2);
        latex.DrawLatex(0.1,0.3,tmp);
        */
        c2->cd();
        
        h_mult_klep->Scale(1.0/nev);
        
        TString plotname(fsig(fsig.Index("/")+1,fsig.Length()-5));
        
        c2->SaveAs("fig/"+plotname+"_core.gif");
        c2->SaveAs("fig/"+plotname+"_core.pdf");
        
        if (writentp) 
        {
                fn->Write();
                fn->Close();
        }
        
  return 0;
}