#include "EventShape.h"
#include "SimpleNtp.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"

Typedefs
typedef std::vector< SimpleCand >	CandList

Functions
void	setStyle ()

void	boostList (CandList &list, CandList &boostlist, TVector3 boost)

bool	isFinalState (int n)

double	dponline (TLorentzVector &v)

void	makePidTable (double pideff, double misid)

void	init (double pideff, double misid)

void	printCand (SimpleCand c)

void	printList (CandList &l, TString tit="")

void	printAllLists ()

int	massCrit (CandList &l, double m, double w)

bool	eventAccepted (CandList &l, double m=0, double w=0)

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	combine (CandList &l1, CandList &l2, CandList &out, int matchPdg=0)

int	combine (CandList &l1, CandList &l2, CandList &l3, CandList &out, int matchPdg=0)

int	combine (CandList &l1, CandList &l2, CandList &l3, CandList &l4, CandList &out, int matchPdg=0)

void	smearMom (SimpleCand &c, double dpr=0.05, double dtht=0.000, double dphi=0.000)

bool	isDetected (double eff=0.90)

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

void	selectMass (CandList &in, CandList &out, double mmean=0, double mw=0)

void	select (CandList &in, CandList &out, int chrg=0, int pdg=0, double mass=0.139)

void	makePidSelection (CandList &l)

void	makeIni4Vector (TLorentzVector &l, double s)

void	configHisto (TH1 *h)

void	qaCand (TString pre, SimpleCand &c, SimpleNtp *n)

void	writeTuple (TString pre, CandList &l, EventShape &evsh, SimpleNtp *n, double mass, double win)

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)

Variables
TRandom3	fRand

TLorentzVector	fIni

const int	fMAX = 100

int	nbins = 300

double	low = 0.0

double	high = 6.0

const double	fNeutEff = 0.95

const double	fdE = 0.05

const double	fNdtht = 0.003

const double	fNdphi = 0.003

int	nsigpiks = 5

int	nsig = 8

int	fPi0Pdg = 111

double	fPi0Mass = 0.13497

double	fPi0win = 0

double	fPi0sig = 0.0055

int	fKsPdg = 310

double	fKsMass = 0.497614

double	fKswin = 0.025

double	fKssig = 0.012

int	fDPdg = 411

double	fDMass = 1.8693

double	fDsig = 0.045

double	fDwin =0.25

int	fD0Pdg = 421

double	fD0Mass = 1.8645

double	fD0sig = 0.055

double	fD0win =0.25

int	fJpsiPdg = 443

double	fJMass = 3.09687

double	fJsig = 0.11

double	fJwin =0.7

int	feePdg = 88880

double	feeMass = 0

double	feesigconst = 0.0024

double	feesigslope = 0.036

double	feewin = 0.

int	fLamPdg = 3122

double	fLamMass = 1.115684

double	fLamsig = 0.0042

double	fLamwin =0.1

int	fDsPdg = 431

double	fDsMass = 1.9686

double	fDssig = 0.04

double	fDswin =0.2

int	fPhiPdg = 333

double	fPhiMass = 1.0195

double	fPhisig = 0.0065

double	fPhiwin =0.2

int	fEtacPdg = 441

double	fEtacMass = 2.9798

double	fEtacsig = 0.08

double	fEtacwin =0.15

int	fLamcPdg = 4122

double	fLamcMass = 2.2849

double	fLamcsig = 0.035

double	fLamcwin =0.2

double	pidtab [25]

double	pidreal [25]

double	mE

double	mMu

double	mPi

double	mK

double	mP

double	mPi0

TF1 *	f_res =0

int	npi

int	nk

int	npr

int	ne

int	nmu

int	nlep

int	evtno

int	globmode

std::map< int, int >	pdgidx

TDatabasePDG *	pdg_db

CandList	mclist

CandList	mcFSlist

CandList	recolist

CandList	recocmslist

CandList	eplus

CandList	eminus

CandList	muplus

CandList	muminus

CandList	piplus

CandList	piminus

CandList	kplus

CandList	kminus

CandList	pplus

CandList	pminus

CandList	gam

Typedef Documentation

typedef std::vector<SimpleCand> CandList

Definition at line 132 of file toy_core.C.

Function Documentation

void boostList	(	CandList &	list,
		CandList &	boostlist,
		TVector3	boost
	)

Definition at line 164 of file toy_core.C.

References c, i, SimpleCand::P4(), and SimpleCand::SetP4().

Referenced by toy_core().

 {
         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);
     }
 }

int combine	(	CandList &	l1,
		CandList &	l2,
		CandList &	out,
		int	matchPdg = `0`
	)

Definition at line 363 of file toy_core.C.

References SimpleCand::AddDau(), c, i, m2(), mclist, SimpleCand::SetCharge(), SimpleCand::SetEvtId(), SimpleCand::SetId(), SimpleCand::SetMarker(), SimpleCand::SetMcPid(), SimpleCand::SetMct(), SimpleCand::SetMotherIdx(), SimpleCand::SetNFS(), SimpleCand::SetP4(), and SimpleCand::SetPdg().

Referenced by toy_core().

 {
         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`
	)

Definition at line 409 of file toy_core.C.

References SimpleCand::AddDau(), c, m2(), mclist, SimpleCand::SetCharge(), SimpleCand::SetEvtId(), SimpleCand::SetId(), SimpleCand::SetMarker(), SimpleCand::SetMcPid(), SimpleCand::SetMct(), SimpleCand::SetMotherIdx(), SimpleCand::SetNFS(), SimpleCand::SetP4(), and SimpleCand::SetPdg().

 {
         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`
	)

Definition at line 474 of file toy_core.C.

References SimpleCand::AddDau(), c, m2(), mclist, SimpleCand::SetCharge(), SimpleCand::SetEvtId(), SimpleCand::SetId(), SimpleCand::SetMarker(), SimpleCand::SetMcPid(), SimpleCand::SetMct(), SimpleCand::SetMotherIdx(), SimpleCand::SetNFS(), SimpleCand::SetP4(), and SimpleCand::SetPdg().

 {
         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 configHisto ( TH1 * h )

Definition at line 704 of file toy_core.C.

Referenced by toy_core().

 {
         h->SetFillStyle(3001);
         h->SetFillColor(602);
         h->SetLineStyle(3);
         h->SetLineColor(602);
 /*      TGaxis *ax = (TGaxis*)h->GetYaxis();
         ax->SetMaxDigits(2);*/
 }

double dponline ( TLorentzVector & v )

Definition at line 189 of file toy_core.C.

References f_res.

Referenced by smearMom().

 {       
         return f_res->Eval(v.Pt())/100.;
 }

bool eventAccepted	(	CandList &	l,
		double	m = `0`,
		double	w = `0`
	)

Definition at line 320 of file toy_core.C.

References i, and m.

Referenced by toy_core().

 {
         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;
 }

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`
	)

Definition at line 337 of file toy_core.C.

References cnt, fIni, i, and m.

Referenced by toy_core().

 {
         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;
 }

void init	(	double	pideff,
		double	misid
	)

Definition at line 235 of file toy_core.C.

References f_res, makePidTable(), mE, mK, mMu, mP, mPi, mPi0, pdg_db, and pdgidx.

Referenced by toy_core().

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

bool isDetected ( double eff = 0.90 )

Definition at line 587 of file toy_core.C.

References fRand.

Referenced by makeRecoCands().

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

bool isFinalState ( int n )

Definition at line 180 of file toy_core.C.

Referenced by makeRecoCands(), and toy_core().

 {
         n=abs(n);
         return (n==11 || n==13  || n==211 || n==321 || n==2212 || n==22);
 }

void makeIni4Vector	(	TLorentzVector &	l,
		double	s
	)

Definition at line 697 of file toy_core.C.

References mP, p, sqrt(), and X.

Referenced by toy_core().

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

void makePidSelection ( CandList & l )

Definition at line 676 of file toy_core.C.

References eminus, eplus, gam, kminus, kplus, mE, mK, mMu, mP, mPi, muminus, muplus, piminus, piplus, pminus, pplus, and select().

Referenced by toy_core().

 {
         // *** 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 makePidTable	(	double	pideff,
		double	misid
	)

Definition at line 195 of file toy_core.C.

References i, pidtab, and printf().

Referenced by init().

 {
         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;
         }
 }

void makeRecoCands	(	CandList &	mc,
		CandList &	reco,
		double	eff = `0.90`,
		double	dp = `0.05`,
		double	dtht = `0.001`,
		double	dphi = `0.001`
	)

Definition at line 594 of file toy_core.C.

References c, SimpleCand::Charge(), fdE, fNdphi, fNdtht, fNeutEff, i, isDetected(), isFinalState(), pid(), SimpleCand::Pid(), SimpleCand::SetMarker(), SimpleCand::SetMass(), SimpleCand::SetPdg(), and smearMom().

Referenced by toy_core().

 {
         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);
                 }
         }
 }

int massCrit	(	CandList &	l,
		double	m,
		double	w
	)

Definition at line 307 of file toy_core.C.

References cnt, and i.

 {
         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;
 }

void printAllLists ( )

Definition at line 283 of file toy_core.C.

References eminus, eplus, gam, kminus, kplus, mclist, muminus, muplus, piminus, piplus, pminus, pplus, printList(), and recolist.

 {
         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");
 }

void printCand ( SimpleCand c )

Definition at line 266 of file toy_core.C.

References SimpleCand::Charge(), SimpleCand::Id(), SimpleCand::McPid(), SimpleCand::Mct(), SimpleCand::MotherIdx(), SimpleCand::NFS(), SimpleCand::NSiblings(), SimpleCand::P4(), and SimpleCand::Pid().

Referenced by printList().

 {
         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;
 }

void printList	(	CandList &	l,
		TString	tit = `""`
	)

Definition at line 275 of file toy_core.C.

References printCand().

Referenced by printAllLists().

 {
         if (tit!="") cout << tit<<endl;
         for (unsigned int j=0;j<l.size();++j) printCand(l[j]);
         cout <<endl;
 }

void qaCand	(	TString	pre,
		SimpleCand &	c,
		SimpleNtp *	n
	)

Definition at line 714 of file toy_core.C.

References SimpleCand::Charge(), SimpleNtp::Column(), fIni, SimpleCand::McPid(), SimpleCand::Mct(), SimpleCand::P4(), SimpleCand::Pdg(), and SimpleCand::Pid().

Referenced by writeTuple().

 {
         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() );
 }

void select	(	CandList &	in,
		CandList &	out,
		int	chrg = `0`,
		int	pdg = `0`,
		double	mass = `0.139`
	)

Definition at line 638 of file toy_core.C.

References c, SimpleCand::Charge(), fRand, i, pdg_db, pdgidx, SimpleCand::Pid(), pidtab, SimpleCand::SetMass(), SimpleCand::SetMcPid(), and SimpleCand::SetPdg().

Referenced by makePidSelection().

 {
         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 selectMass	(	CandList &	in,
		CandList &	out,
		double	mmean = `0`,
		double	mw = `0`
	)

Definition at line 628 of file toy_core.C.

References fabs(), and i.

Referenced by toy_core().

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

void setStyle ( )

Definition at line 141 of file toy_core.C.

Referenced by toy_core().

 {
         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);
 }

void smearMom	(	SimpleCand &	c,
		double	dpr = `0.05`,
		double	dtht = `0.000`,
		double	dphi = `0.000`
	)

Definition at line 557 of file toy_core.C.

References dponline(), fRand, m, p, SimpleCand::P4(), and SimpleCand::SetP4().

Referenced by makeRecoCands().

 {
   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);
         
 }

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`
	)

Definition at line 841 of file toy_core.C.

References boostList(), c, c1, c2, EventShape::ChrgPSumCms(), EventShape::ChrgPtSumLab(), combine(), configHisto(), Double_t, eminus, eplus, eventAccepted(), evtno, f, fD0Mass, fD0Pdg, fD0sig, fD0win, fDMass, fDPdg, fDsig, fDsMass, fDsPdg, fDssig, fDswin, fDwin, feesigconst, feesigslope, feewin, fEtacMass, fEtacPdg, fEtacsig, fEtacwin, filename, fillHisto(), fIni, fJMass, fJpsiPdg, fJsig, fJwin, fKsMass, fKsPdg, fKssig, fKswin, fLamcMass, fLamcPdg, fLamcsig, fLamcwin, fLamMass, fLamPdg, fLamsig, fLamwin, fMAX, EventShape::FoxWolfMomR(), fPhiMass, fPhiPdg, fPhisig, fPhiwin, fPi0Mass, fPi0Pdg, fPi0sig, fPi0win, fRand, gam, globmode, high, i, init(), isFinalState(), kminus, kplus, m, makeIni4Vector(), makePidSelection(), makeRecoCands(), max(), mcFSlist, mclist, min(), mode, EventShape::MultPminCms(), muminus, muplus, nbins, EventShape::NCharged(), ne, EventShape::NeutESumCms(), nk, nlep, nmu, EventShape::NNeutral(), EventShape::NParticles(), npi, npr, nsig, nsigpiks, offset(), p, pdg_db, piminus, piplus, plotname, EventShape::PmaxCms(), EventShape::PmaxLab(), EventShape::PminCms(), pminus, pos, pplus, pt(), EventShape::Ptmax(), EventShape::PtSumLab(), pz, recocmslist, recolist, selectMass(), SimpleCand::SetDau(), SimpleCand::SetEvtId(), SimpleCand::SetMotherIdx(), SimpleCand::SetNSiblings(), setStyle(), EventShape::SumChrgPminCms(), t, TString, and writeTuple().

 {
         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;
 }

void writeTuple	(	TString	pre,
		CandList &	l,
		EventShape &	evsh,
		SimpleNtp *	n,
		double	mass,
		double	win
	)

Definition at line 781 of file toy_core.C.

References EventShape::Aplanarity(), EventShape::ChrgPSumCms(), EventShape::ChrgPtSumLab(), EventShape::Circularity(), SimpleNtp::Column(), d, SimpleNtp::DumpData(), evtno, fabs(), EventShape::FoxWolfMomR(), globmode, i, EventShape::NCharged(), ne, nk, nmu, EventShape::NNeutral(), EventShape::NParticles(), npi, npr, EventShape::Planarity(), EventShape::PmaxCms(), EventShape::PmaxLab(), EventShape::PminCms(), EventShape::Ptmax(), EventShape::PtSumLab(), qaCand(), EventShape::Sphericity(), tag, and EventShape::Thrust().

Referenced by toy_core().

 {
         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();
         }
 }

Variable Documentation

CandList eminus

Definition at line 138 of file toy_core.C.

Referenced by makePidSelection(), printAllLists(), and toy_core().

CandList eplus

Definition at line 138 of file toy_core.C.

Referenced by makePidSelection(), printAllLists(), and toy_core().

int evtno

Definition at line 124 of file toy_core.C.

Referenced by toy_core(), and writeTuple().

TF1* f_res =0

Definition at line 123 of file toy_core.C.

Referenced by dponline(), and init().

double fD0Mass = 1.8645

Definition at line 65 of file toy_core.C.

Referenced by toy_core().

int fD0Pdg = 421

Definition at line 64 of file toy_core.C.

Referenced by toy_core().

double fD0sig = 0.055

Definition at line 66 of file toy_core.C.

Referenced by toy_core().

double fD0win =0.25

Definition at line 67 of file toy_core.C.

Referenced by toy_core().

const double fdE = 0.05

Definition at line 39 of file toy_core.C.

Referenced by makeRecoCands().

double fDMass = 1.8693

Definition at line 60 of file toy_core.C.

Referenced by toy_core().

int fDPdg = 411

Definition at line 59 of file toy_core.C.

Referenced by toy_core().

double fDsig = 0.045

Definition at line 61 of file toy_core.C.

Referenced by toy_core().

double fDsMass = 1.9686

Definition at line 86 of file toy_core.C.

Referenced by toy_core().

int fDsPdg = 431

Definition at line 85 of file toy_core.C.

Referenced by toy_core().

double fDssig = 0.04

Definition at line 87 of file toy_core.C.

Referenced by toy_core().

double fDswin =0.2

Definition at line 88 of file toy_core.C.

Referenced by toy_core().

double fDwin =0.25

Definition at line 62 of file toy_core.C.

Referenced by toy_core().

double feeMass = 0

Definition at line 75 of file toy_core.C.

int feePdg = 88880

Definition at line 74 of file toy_core.C.

double feesigconst = 0.0024

Definition at line 76 of file toy_core.C.

Referenced by toy_core().

double feesigslope = 0.036

Definition at line 77 of file toy_core.C.

Referenced by toy_core().

double feewin = 0.

Definition at line 78 of file toy_core.C.

Referenced by toy_core().

double fEtacMass = 2.9798

Definition at line 96 of file toy_core.C.

Referenced by toy_core().

int fEtacPdg = 441

Definition at line 95 of file toy_core.C.

Referenced by toy_core().

double fEtacsig = 0.08

Definition at line 97 of file toy_core.C.

Referenced by toy_core().

double fEtacwin =0.15

Definition at line 98 of file toy_core.C.

Referenced by toy_core().

TLorentzVector fIni

Definition at line 29 of file toy_core.C.

Referenced by fillHisto(), qaCand(), and toy_core().

double fJMass = 3.09687

Definition at line 70 of file toy_core.C.

Referenced by toy_core().

int fJpsiPdg = 443

Definition at line 69 of file toy_core.C.

Referenced by toy_core().

double fJsig = 0.11

Definition at line 71 of file toy_core.C.

Referenced by toy_core().

double fJwin =0.7

Definition at line 72 of file toy_core.C.

Referenced by toy_core().

double fKsMass = 0.497614

Definition at line 55 of file toy_core.C.

Referenced by toy_core().

int fKsPdg = 310

Definition at line 54 of file toy_core.C.

Referenced by toy_core().

double fKssig = 0.012

Definition at line 57 of file toy_core.C.

Referenced by toy_core().

double fKswin = 0.025

Definition at line 56 of file toy_core.C.

Referenced by toy_core().

double fLamcMass = 2.2849

Definition at line 101 of file toy_core.C.

Referenced by toy_core().

int fLamcPdg = 4122

Definition at line 100 of file toy_core.C.

Referenced by toy_core().

double fLamcsig = 0.035

Definition at line 102 of file toy_core.C.

Referenced by toy_core().

double fLamcwin =0.2

Definition at line 103 of file toy_core.C.

Referenced by toy_core().

double fLamMass = 1.115684

Definition at line 81 of file toy_core.C.

Referenced by toy_core().

int fLamPdg = 3122

Definition at line 80 of file toy_core.C.

Referenced by toy_core().

double fLamsig = 0.0042

Definition at line 82 of file toy_core.C.

Referenced by toy_core().

double fLamwin =0.1

Definition at line 83 of file toy_core.C.

Referenced by toy_core().

const int fMAX = 100

Definition at line 32 of file toy_core.C.

Referenced by toy_core().

const double fNdphi = 0.003

Definition at line 41 of file toy_core.C.

Referenced by makeRecoCands().

const double fNdtht = 0.003

Definition at line 40 of file toy_core.C.

Referenced by makeRecoCands().

const double fNeutEff = 0.95

Definition at line 38 of file toy_core.C.

Referenced by makeRecoCands().

double fPhiMass = 1.0195

Definition at line 91 of file toy_core.C.

Referenced by toy_core().

int fPhiPdg = 333

Definition at line 90 of file toy_core.C.

Referenced by toy_core().

double fPhisig = 0.0065

Definition at line 92 of file toy_core.C.

Referenced by toy_core().

double fPhiwin =0.2

Definition at line 93 of file toy_core.C.

Referenced by toy_core().

double fPi0Mass = 0.13497

Definition at line 50 of file toy_core.C.

Referenced by toy_core().

int fPi0Pdg = 111

Definition at line 49 of file toy_core.C.

Referenced by toy_core().

double fPi0sig = 0.0055

Definition at line 52 of file toy_core.C.

Referenced by toy_core().

double fPi0win = 0

Definition at line 51 of file toy_core.C.

Referenced by toy_core().

TRandom3 fRand

Definition at line 28 of file toy_core.C.

Referenced by isDetected(), select(), smearMom(), and toy_core().

CandList gam

Definition at line 138 of file toy_core.C.

Referenced by makePidSelection(), printAllLists(), and toy_core().

int globmode

Definition at line 124 of file toy_core.C.

Referenced by toy_core(), and writeTuple().

double high = 6.0

Definition at line 35 of file toy_core.C.

Referenced by toy_core().

CandList kminus

Definition at line 138 of file toy_core.C.

Referenced by makePidSelection(), printAllLists(), and toy_core().

CandList kplus

Definition at line 138 of file toy_core.C.

Referenced by makePidSelection(), printAllLists(), and toy_core().

double low = 0.0

Definition at line 34 of file toy_core.C.

CandList mcFSlist

Definition at line 135 of file toy_core.C.

Referenced by toy_core().

CandList mclist

Definition at line 135 of file toy_core.C.

Referenced by combine(), printAllLists(), and toy_core().

double mE

Definition at line 121 of file toy_core.C.

Referenced by init(), and makePidSelection().

double mK

Definition at line 121 of file toy_core.C.

Referenced by init(), and makePidSelection().

double mMu

Definition at line 121 of file toy_core.C.

Referenced by init(), and makePidSelection().

double mP

Definition at line 121 of file toy_core.C.

Referenced by init(), makeIni4Vector(), and makePidSelection().

double mPi

Definition at line 121 of file toy_core.C.

Referenced by init(), and makePidSelection().

double mPi0

Definition at line 121 of file toy_core.C.

Referenced by init().

CandList muminus

Definition at line 138 of file toy_core.C.

Referenced by makePidSelection(), printAllLists(), and toy_core().

CandList muplus

Definition at line 138 of file toy_core.C.

Referenced by makePidSelection(), printAllLists(), and toy_core().

int nbins = 300

Definition at line 33 of file toy_core.C.

Referenced by toy_core().

int ne

Definition at line 124 of file toy_core.C.

Referenced by toy_core(), and writeTuple().

int nk

Definition at line 124 of file toy_core.C.

Referenced by PndLmdTrkQTask::Exec(), main(), softtrigger_toy12(), toy_core(), and writeTuple().

int nlep

Definition at line 124 of file toy_core.C.

Referenced by softtrigger_toy12(), and toy_core().

int nmu

Definition at line 124 of file toy_core.C.

Referenced by toy_core(), and writeTuple().

int npi

Definition at line 124 of file toy_core.C.

Referenced by ana_Lambda_fit(), PndHypFullAna::Exec(), PndHypSimpleAna::Exec(), PndHypFullIdealAna::Exec(), softtrigger_toy12(), toy_core(), and writeTuple().

int npr

Definition at line 124 of file toy_core.C.

Referenced by softtrigger_toy12(), toy_core(), and writeTuple().

int nsig = 8

int nsigpiks = 5

Definition at line 44 of file toy_core.C.

Referenced by toy_core().

TDatabasePDG* pdg_db

Definition at line 130 of file toy_core.C.

Referenced by init(), select(), and toy_core().

std::map<int,int> pdgidx

Definition at line 128 of file toy_core.C.

Referenced by init(), and select().

double pidreal[25]

Initial value:

= {  95,  6,  10,   7,  10,    
                         5,  86,  27,   6,   7,    
                         9,  29,  85,  18,  16,    
                         9,  21,  40,  67,  24,    
                         12, 12,  20,  16,  90 }

Definition at line 114 of file toy_core.C.

double pidtab[25]

Initial value:

= {  80,   5,   5,   5,   5,    
                         5,  80,   5,   5,   5,    
                         5,   5,  80,   5,   5,    
                         5,   5,   5,  80,   5,    
                         5,   5,   5,   5,  80 }

Definition at line 107 of file toy_core.C.

Referenced by makePidTable(), and select().

CandList piminus

Definition at line 138 of file toy_core.C.

Referenced by makePidSelection(), printAllLists(), and toy_core().

CandList piplus

Definition at line 138 of file toy_core.C.

Referenced by makePidSelection(), printAllLists(), and toy_core().

CandList pminus

Definition at line 138 of file toy_core.C.

Referenced by makePidSelection(), printAllLists(), and toy_core().

CandList pplus