MC Event Generation Information in AATrack package
           --------------------------------------------------

MC Event Generation Information is contained in the following classes:

class AA::PtlMC;
class AA::VrtMC;

Auxilary classes:

  typedef std::list<PtlMC*> PtlMCLst;
  typedef std::list<PtlMC*>::iterator PtlMCLstIt;
  typedef std::list<PtlMC*>::const_iterator PtlMCLstCIt;
  typedef std::list<PtlMC*>::reverse_iterator PtlMCLstRIt;

  typedef std::list<VrtMC*> VrtMCLst;
  typedef std::list<VrtMC*>::iterator VrtMCLstIt;
  typedef std::list<VrtMC*>::const_iterator VrtMCLstCIt;
  typedef std::list<VrtMC*>::reverse_iterator VrtMCLstRIt;

The objects of these types are stored in:

class AA::PtlMCBox;

The complete description of these classes can be found in:

#include "PtlMC.hpp"

The most usefule methods of PtlMCBox class are:

namespace AA{
  class AA::PtlMCBox{
    int sizePtlMC(); //Returns number of PtlMC objects stored
    int sizeVrtMC(); //Returns number of VrtMC objects stored
    
// Get the list of pointers to all stored PtlMC objects
// Optional: giving idPDG != 0 will select only particles of given
// type. idPDG - Id number of particle type, following PDG recepie
// also followed by Dzero.

    void particles (PtlMCLst& lst, int idPDG = 0) const;

// Fast method of getting the list of all stored objects

    const PtlMCLst* particles() const {return &_boxPtl;}
  }
}
     
The most useful methods of PtlMC class are:

namespace AA{
  class PtlMC {
  public:

//Get 4-momentum of particle (px,py,pz,e)
    const HepVector& mom4() const;

//Get i-th component of 4-momentum
    double mom4(int i) const;

//Get PDG index (also Dzero convention)
    int idPDG() const;

//Get pointer to the origin vertex of particle
//returns 0 if no origin vertex exist
    VrtMC* vrtStart() const;

//Get pointer to the end vertex of particle
//returns 0 if no origin vertex exist
    VrtMC* vrtEnd() const {return _pen;}

// returns true if given particle is from Geant output
    bool geantOutput() const {return _geant;}

// Total momentum of particle
    double ptot() const;

// Transverse momentum
    double pt()   const;

// pseudorapidity
    double eta()  const;

// Axial angle of particle
    double phi()  const;

// Produce State5 object for given particle
    void state5(State5* pst);

// Associate given MC particle to the recosntructed particle
// From the list TrkLst lst. This list can be obtained like this:
//
// TrkLst lst;
// AA::ptlBox.particles(lst);
//
// Returns pointer to reconstructed track in Trk* ptr.
// Trk* ptr can be converted to pointer to particle:
//
// Ptl* ptl = AA::particle(ptr);
//
// Returns chi2 of match between generated and reconstructed 
// particles. If match is not found (chi2 > 100) ptr is set to 0;

    void associate(const TrkLst& lst, Trk*& ptr, double& chi2);

// Printout of the PtlMC content. The complete decay chain
// of given particle is printed out.
    void debug() const;

  }
}

The most useful methods of VrtMC class are:

namespace AA{
  class VrtMC {
  public:

//List of pointers to all parents of given vertex;
    const std::list<PtlMC*>& parents() const;

//List of pointers to all children of given vertex;
    const std::list<PtlMC*>& children() const{return _vc;}

//Position of vertex
    const HepVector& x() const {return _v;}

// i-th component of position
    double x(int i) const {return _v(i);}

// Printout of information in VrtMC. The complete decay
// chain is printed out.
    void debug() const;
  }
}

An example of access to MC information can look like this:

// In this method we find PtlMC corresponding to B-hadron decaying
// B -> D* mu nu X. D* should decay to D0 pi, D0 -> Kpi. 
// The pointer to this particle is returned in PtlMC* ptl;
// returns false if given particle is not found.

bool checkMC(PtlMC*& ptl);
void checkBDst(VrtMC* pv, bool& okmu, bool& okD, bool& okDst);

bool checkMC(PtlMC*& ptl){

//Loop over all PtlMC objects of given event
  for(AA::PtlMCLstCIt p = AA::ptlMCBox.particles()->begin();
      p != AA::ptlMCBox.particles()->end(); ++p){
    ptl = *p;

//Select B-hadron (Bu, Bd, Bs, Bc)
    if(abs(ptl->idPDG()) == 511 || abs(ptl->idPDG()) == 521 ||
       abs(ptl->idPDG()) == 531 || abs(ptl->idPDG()) == 541 ) {

      bool okmu  = false;
      bool okD   = false;
      bool okDst = false;

//Decay vertex of given B-hadron
      VrtMC* ven = ptl->vrtEnd();

//Check the decay type of given B-hadron
      checkBDst(ven,okmu,okD,okDst);

// If B hadron decays to muon (okmu==true), D* -> D0 pi (okDst == true)
// and D0->Kpi (okD = true), the required B-hadron is found
      if(okmu && okD && okDst) return true;
    }
  }

// B-hadron is not found
  return false;
}

void checkBDst(VrtMC* pv, bool& okmu, bool& okD, bool& okDst){

//List of all children of given vertex pv
  const PtlMCLst* pch = &pv->children();

//Loop over all children
  for(PtlMCLstCIt p = pch->begin(); p != pch->end(); ++p){
    PtlMC* ptl = *p;

    //Check Muon in the decay products
    if(abs(ptl->idPDG()) == 13 && okmu == false) {
      okmu = true;
    }

    //Check D* -> D0 pi in the decay products
    if(abs(ptl->idPDG()) == 413 && okDst == false) {

      //Decay vertex of D*
      VrtMC* ven = ptl->vrtEnd();

      if(ven != 0) {
	okDst = true;

        //children of D* decay
	const PtlMCLst* pchD = &ven->children();

        //number of decay products of D*
	if(pchD->size() != 2) okDst = false;

        Loop over all decay products of D*
	for(PtlMCLstCIt p1 = pchD->begin(); p1 != pchD->end(); ++p1){
	  PtlMC* ptl1 = *p1;

          //Require D* -> D0 pi decay
	  if(abs(ptl1->idPDG()) != 211 && 
	     abs(ptl1->idPDG()) != 421) okDst = false;
	}
      }
    }

    //Check D0 -> K pi decay
    if(abs(ptl->idPDG()) == 421 && okD == false) {

      // Decay vertex of D0
      VrtMC* ven = ptl->vrtEnd();

      if(ven != 0) {
	okD = true;

        //Children of D0 decay
	const PtlMCLst* pchD = &ven->children();

        //Number of decay products
	if(pchD->size() != 2) okD = false;

        //Loop over all children
	for(PtlMCLstCIt p1 = pchD->begin(); p1 != pchD->end(); ++p1){
	  PtlMC* ptl1 = *p1;

          //Require D0 -> Kpi decay
	  if(abs(ptl1->idPDG()) != 211 && 
	     abs(ptl1->idPDG()) != 321) okD = false;
	}
      }
    }
    
    //Recursive loop over complete decay chain
    if(ptl->vrtEnd() != 0) checkBDst(ptl->vrtEnd(),okmu,okD,okDst);
  }
}