2.4. Monte Carlo Truth Match
==============================
A Monte Carlo Truth Match is a check whether your reconstructed decay tree resembles the original one created by the generator concerning mass hypothesis assignments and topology. It can of course only be applied for simulated data, where the truth about a certain decay is known. This match is of particular
importance whenever the number of reconstructed signal particles in data has to be determined precisely, since it might be the only way to access the
reconstruction efficiency of the analysis procedure.(Example code in **tut_ana_mcmatch.C/tut_ana.C** in the tutorial directory.) 

To perform a proper and complete MC truth match sometimes can be tricky, in particular when you reconstruct complex decay trees. In our case it looks like this::

                      pi-           mu+
                    /             /
                  /             /
  p pbar  ---> psi(2S)  ---> J/psi
                  \             \
                    \             \
                      pi+            mu-

You can see, that this tree has four final state particles (*mu+, mu-, pi+_* and *_*pi-*) and two intermediate resonances (*psi(2S)* and *J/psi*), and for a complete tree match, we need
to identify the counterparts for all these particles in the MC list. We however have to take into account, that initially a ``pbarpSystem`` instead of a ``psi(2S)`` was generated.

2.4.1. Accessing the MC truth
-------------------------------

The starting point for doing any kind of MC truth match is to find the counterpart of the reconstructed final state particles in the list of Monte Carlo Truth objects. Therefore we can simply ask any final state reco particle for its Mc Truth counterpart, being
also a **RhoCandidate** object, with the code::

  RhoCandList muplus;
  ...
  // ... in event loop ...
    theAnalysis->FillList( muplus, "MuonLoosePlus");
    ...
    RhoCandidate *truth = muplus[i]->GetMcTruth();                 // *** here we ask for the MC truth object for track i
    if ( 0x0!=truth )  // *** check whether the truth object exists, in which case the match was successful
    { 
      Double_t p_diff = truth->P() - muplus[i]->P();               // *** do whatever kind of calculations ...
      ...

There might be even cases where such a MC truth particle does not really exist, since a reconstructed track might be based on noise hits only
(hopefully in rare cases). Thus it should always be checked, whether the pointer ``truth`` is meaningful (i.e. unequal to zero).

This kind of match is useful, if you e.g. want to compare the reconstructed and the generated particle to study mass or momentum resolution etc.

2.4.2. Full MC truth tree match
---------------------------------
As mentioned above a complete match needs a bit more effort than just matching the PID of the final state particles.

What needs to be checked in our example case is, that

   * the final state particles (*mu+, mu-, pi+, pi-*) have the correct PID
   * the muons have the same mother
   * this mother is of type *J/psi*
   * the *J/psi* doesn't have any additional daughters
   * the *J/psi* and the pions have the same mother
   * the type of this one is *pbarpSystem*
   * and the *pbarpSystem* doesn't have any additional daughters

If you don't know immediately how to perform such kind of match, don't worry! Fortunately there is a method available within **PndAnalysis** which does this job for you, called ``PndAnalysis::McTruthMatch``. This method needs a particle as input and does the tree match for arbitrary decay trees.

In order to use that full match by calling ``PndAnalysis::McTruthMatch`` some preparation has to be done during combinatorics. Since the matching routine also checks whether intermediate resonances
have the correct type (at least in default setting), we have to tell our composite candidates, of which type they are supposed to be, either by using the methods ``RhoCandidate::SetType`` or for a full list ``RhoCandList::SetType``, or even more comfortable, just together with the combinatorics. It is strongly recommended to initialize the singleton of **TDatabasePDG** with the EvtGen properties beforehand, since otherwise the names you specify might have different codes in the default TDatabasePDG setup than in EvtGen, leading to failing truth matches.
This is how it could look like for our example::

  ...
  // initialize TDatabasePDG with EvtGen names 
  // (otherwise there are inconsistencies between generator and analysis)
  TString wd = gSystem->Getenv("VMCWORKDIR");
  TDatabasePDG::Instance()->ReadPDGTable(wd+"/input/pdg_table_evtgen.txt");
  
  // create shortcut for TDatabasePDG::Instance()
  TDatabasePDG *pdb = TDatabasePDG::Instance();

  RhoCandList muplus, muminus, piplus, piminus, jpsi, psi2s;
  ...
    // ... in event loop ...
    // fill all the lists
    theAnalysis->FillList(muplus,  "MuonAllPlus");
    theAnalysis->FillList(muminus, "MuonAllMinus");
    theAnalysis->FillList(piplus,  "PionAllPlus");
    theAnalysis->FillList(piminus, "PionAllMinus");

    // do the combinatorics for J/psi and set type of composite to "J/psi"
    jpsi.Combine(muplus, muminus, pdb->GetParticle("J/psi")); 
    
    // do the combinatorics for psi(2S) and set type of composite to "pbarpSystem"
    psi2s.Combine(jpsi, piplus, piminus, pdb->GetParticle("pbarpSystem"));
    
    for (i=0; i<psi2s.GetLength(); ++i)
    if ( theAnalysis->McTruthMatch(psi2s[i]) ) // *** here the match is done
    {
        // ... do whatever
  ...

Another nice feature of the truth match routine is, that the match assigns an MC truth object also to the intermediate composite
states in addition. This means you are able to access e.g. the true *J/psi* from above with this code::

  ...
    if ( theAnalysis->McTruthMatch(psi2s[i]) ) // *** here the match is done
    {
      // get the true intermediate J/psi
      RhoCandidate *trueJpsi = psi2s[j]->Daughter(0)->GetMcTruth();

      // get the true initial psi(2S) 
      RhoCandidate *truePsi  = psi2s[i]->GetMcTruth();              

      // e.g. compute the mass difference of the true object to its reco counterpart
      Double mass_diff    = truePsi->M() - psi2s[i]->M();    
  ...

The mass spectra for truth matched J/psi and psi(2S) look like

.. image:: ./Images/hjpsim_ftm.png
   :width: 350
   :alt: J/psi with full truth match

.. image:: ./Images/hpsim_ftm.png
   :width: 350
   :alt: psi(2S) with full truth match

The wrong combinatorics, in particular formerly visible in the *J/psi* spectrum, are gone.

2.4.3. MC Truth List
------------------------
In case you want to study decay patterns directly from monte carlo truth without doing something with the reco particles, you also have the possibility to get the mc truth list containing all particles from generation and simulation with proper mother-daughter relations. (Example code in **tut_ana_mclist.C** in the tutorial directory.) 

This can be done with the following code::

  RhoCandList mclist;
  ...
    // ... in event loop ...
    theAnalysis->FillList(mclist,  "McTruth");

If you e.g. want to print out the mc truth list with mother-daughter relations, you can do it with something like this code::

  // *** print mc list
  for (int j=0;j<mclist.GetLength();++j)
  {
    RhoCandidate *mcmother = mclist[j]->TheMother(); // mother of mc particle         
    int muid = -1;
    if (0x0!=cmoth) muid = cmoth->GetTrackNumber();  // track ID of mother
        
    cout << "Track " << mclist[j]->GetTrackNumber()
         << " (PDG:" << mclist[j]->PdgCode() << ") has mother "<<muid;

    if (mclist[j]->NDaughters()>0) 
      cout <<" and daughter(s) ";

    for (k=0;k<mclist[j]->NDaughters();++k) 
      cout <<mclist[j]->Daughter(k)->GetTrackNumber()<<"  ";
      
    cout<<endl;        
  }