// *** TCandLists for the analysis
TCandList eplus, eminus;

...

while (theAnalysis->GetEvent() && i++<nevts){
    if ((i%1000)==0) cout << "Number of events analysed: " << i << endl;

    // *** Select with no PID info ('All'); type and mass are set
    theAnalysis->FillList(eplus, "ElectronTightPlus");
    //theAnalysis->FillList(eplus, "ElectronTightPlus","PidAlgoEmcBayes");
    theAnalysis->FillList(eminus, "ElectronTightPlus");
    //theAnalysis->FillList(eminus, "ElectronTightMinus","PidAlgoEmcBayes");

    // ELECTRONS ***
    for (Int_t j=0;j<eminus.GetLength();++j){
      // Fill momentum of reconstructed particles
      hP->Fill(eminus[j].GetMomentum().Mag());
      // Fill energy of reconstructed particles
      hE->Fill(eminus[j].Energy());
      // Fill momentum vs E/p
      hEP->Fill(eminus[j].GetMomentum().Mag(), eminus[j].Energy()/eminus[j].GetMomentum().Mag());
      // *** CM frame *** //
      // boost from LAB to CM
      lv_charged_cm.SetVect(eminus[j].GetMomentum());
      lv_charged_cm.SetE(TMath::Sqrt(cElectronMass*cElectronMass + eminus[j].GetMomentum().Mag()*eminus[j].GetMomentum().Mag()));
      lv_charged_cm.Boost(-(lv_bar_beam + lv_p_target).BoostVector()); 

      hCosTheta_N_CM->Fill(lv_charged_cm.CosTheta());
    }
    for (Int_t j=0;j<eplus.GetLength();++j){
      // Fill momentum of reconstructed particles
      hP->Fill(eplus[j].GetMomentum().Mag());
      // Fill energy of reconstructed particles
      hE->Fill(eplus[j].Energy());
      // Fill momentum vs E/p
      hEP->Fill(eplus[j].GetMomentum().Mag(), eplus[j].Energy()/eplus[j].GetMomentum().Mag());
      // *** CM frame *** //
      // boost from LAB to CM
      lv_charged_cm.SetVect(eplus[j].GetMomentum());
      lv_charged_cm.SetE(TMath::Sqrt(cElectronMass*cElectronMass + eplus[j].GetMomentum().Mag()*eplus[j].GetMomentum().Mag()));
      lv_charged_cm.Boost(-(lv_bar_beam + lv_p_target).BoostVector()); 

      hCosTheta_P_CM->Fill(lv_charged_cm.CosTheta());
    }
  }