#include "TClonesArray.h"
#include "tmb_tree/TMBMCvtx.hpp"
#include "tmb_tree/TMBMCpart.hpp"
#include "TStopwatch.h"
#include "TTree.h"

#include <iostream>
using namespace std;

bool a_parent_is (const TMBMCpart *p, int pid)
{
    while (p != 0) {
      if (p->pdgid() == pid) {
        return true;
      }
      const TMBMCvtx *v = p->getPMCvtx();
      if (v == 0) {
        return false;
      }
      if (v->nparents() == 0) {
        return false;
      }
      p = v->getParent(0);
    }
    return false;
}

void test_3_count_w_to_c_finally (void)
{
  TStopwatch watcher;
  watcher.Start();

  TFile *f = new TFile ("caf_zcc.root", "READ");

  TTree *tree = (TTree*) f->Get("TMBTree");
  TClonesArray *cpart = new TClonesArray ("TMBMCpart");
  TClonesArray *cvtx = new TClonesArray ("TMBMCvtx");

  tree->SetBranchAddress ("MCpart", &cpart);
  tree->SetBranchAddress ("MCvtx", &cvtx);

  TBranch *b_part = tree->GetBranch("MCpart");
  TBranch *b_vtx = tree->GetBranch("MCvtx");

  int nevents = tree->GetEntries();

  //tree->SetBranchStatus ("*", 0); // doesn't seem to actually work! Weird!

  int count_events = 0;
  int count_occurs = 0;
  for (int i = 0; i < nevents; i++) {
    b_part->GetEntry(i);
    b_vtx->GetEntry(i);

    //
    // Now, loop through all the particles lookign for a charm.
    //

    bool found_it = false;
    for (int i_part = 0; i_part < cpart->GetEntriesFast(); i_part++) {
      TMBMCpart *p = static_cast<TMBMCpart*> (cpart->At(i_part));
      if (p->pdgid() == 13) {
        if (a_parent_is (p, 4)) {
          found_it = true;
          count_occurs++;
        }
      }
    }

    if (found_it) {
      count_events++;
    }

  }

  cout << "Number of times: " << count_occurs << endl;
  cout << "Number of events: " << count_events << endl;

  watcher.Stop();
  watcher.Print ("m");
  
}