R__LOAD_LIBRARY ($PLUTOLIBDIR/libPluto.so)

int dielectron_readtree ()
{
  TClonesArray* partArray = new TClonesArray ("PParticle", 10);
  PParticle* parts[10] ;

  TFile* fileIn = new TFile ("omegaphi_ee.root");
  TTree* tree = (TTree*) fileIn->Get ("data");
  tree->SetBranchAddress ("Particles", &partArray );

  TH1F* hMass = new TH1F ("hmass", "", 150, 0., 1.5);

  for (int iR = 0 ; iR < tree->GetEntries() ; iR++ )
  {
    tree->GetEntry (iR);
    PParticle* ep_om, * em_om , * ep_phi , * em_phi ;

    for (int iP = 0; iP < partArray->GetEntries() ; iP++ )
      parts[iP] = (PParticle*) partArray->At (iP);

    for (int iP = 0; iP < partArray->GetEntries() ; iP++ ) {
      int pid         = parts[iP] -> ID() ,
          parentIndex = parts[iP] -> GetParentIndex();
      if ( parentIndex < 0 ) continue;
      if ( parts[parentIndex]->ID() == 52 ) {            // parent = omega
        if ( pid == 2 ) ep_om  = parts[iP] ;
        if ( pid == 3 ) em_om  = parts[iP] ;
      }
      if ( parts[parentIndex]->ID() == 55 ) {            // parent = phi
        if ( pid == 2 ) ep_phi = parts[iP] ;
        if ( pid == 3 ) em_phi = parts[iP] ;
      }
    }
    PParticle omega = *ep_om  + *em_om  ,   // caution: W of parent gets W^2 of children
              phi   = *ep_phi + *em_phi ;

    hMass->Fill ( omega.M() , ep_om ->W() );
    hMass->Fill ( phi  .M() , ep_phi->W() );
  }
  return 0;
}