MuParentHelper.cxx

Go to the documentation of this file.

#ifndef muparenthelper_cxx
#define muparenthelper_cxx
#include "StandardNtuple/NtpStRecord.h"
#include "MCNtuple/NtpMCRecord.h"
#include "MCNtuple/NtpMCTruth.h"
#include "MCNtuple/NtpMCStdHep.h"
#include "MCNtuple/NtpMCFluxInfo.h"
#include "AnalysisNtuples/ANtpTruthInfoBeam.h"
#include "MCReweight/MuParentHelper.h"
#include "MCReweight/NuParent.h"
#include <iostream>
#include "TString.h"
#include "TClonesArray.h"

Int_t GEANTtoPDG(Int_t);

using namespace std;

//****************************
MuParentHelper::MuParentHelper()
{
  fok = false;
  fFile = NULL;
  fTree = NULL;
  fusemini = true;
  char *dir = getenv("GNUMIAUX");
  if(dir!=NULL) fdir = dir;
  else fdir = ".";
  if(OpenFile()) {
    SetAddresses();
    SetVectorEntries();
  }
}

//****************************
MuParentHelper::~MuParentHelper()
{
  CloseFile();
}

//****************************
Bool_t MuParentHelper::SetFileDir(std::string dir, bool usemini)
{
  fdir = dir;
  fusemini=usemini;
  if(OpenFile()) {
    SetAddresses();
    SetVectorEntries();
    return true;
  }
  return false;
} 

//****************************
Bool_t MuParentHelper::OpenFile()
{
  if(fFile) {
    fFile->Close();
    delete fFile;       //delete TFile
  }
  char name[256];
  char slash = '/';
  TString st(fdir.c_str());  
  st.Remove(0,st.Last(slash)+1);
  if(fusemini){
    sprintf(name,"%s/mupiTreeMini_%s.root",fdir.c_str(),st.Data());
  }
  else{
    sprintf(name,"%s/mupiTree_%s.root",fdir.c_str(),st.Data());
  }
  cout << name << endl;
  fFile = new TFile(name,"READ");
  if(fFile->IsOpen() && !fFile->IsZombie()){
    fTree = (TTree*) fFile->Get("mupitree");
    fok = true;
    return true;
  }
  return false;
}

//****************************
void MuParentHelper::SetAddresses()
{
  if(!fok) return;  
  fTree->SetBranchAddress("run",&fRunNum);
  fTree->SetBranchAddress("evtno",&fEvtNum);
  fTree->SetBranchAddress("tvx",&fX);
  fTree->SetBranchAddress("tvy",&fY);
  fTree->SetBranchAddress("tvz",&fZ);
  fTree->SetBranchAddress("tpx",&fPX);
  fTree->SetBranchAddress("tpy",&fPY);
  fTree->SetBranchAddress("tpz",&fPZ);
  fTree->SetBranchAddress("tptype",&fPID);
  fTree->SetBranchAddress("tgen",&fGen);
  if(!fusemini){
    fTree->SetBranchAddress("Ndxdz",&fNdxdz);
    fTree->SetBranchAddress("Ndydz",&fNdydz);
    fTree->SetBranchAddress("Npz",&fNpz);
    fTree->SetBranchAddress("Nenergy",&fNenergy);
    fTree->SetBranchAddress("Ndxdznea",&fNdxdznear);
    fTree->SetBranchAddress("Ndydznea",&fNdydznear);
    fTree->SetBranchAddress("Nenergyn",&fNenergynear);
    fTree->SetBranchAddress("Nwtnear",&fNwtnear);
    fTree->SetBranchAddress("Ndxdzfar",&fNdxdzfar);
    fTree->SetBranchAddress("Ndydzfar",&fNdydzfar);
    fTree->SetBranchAddress("Nenergyf",&fNenergyfar);
    fTree->SetBranchAddress("Nwtfar",&fNwtfar);
    fTree->SetBranchAddress("Norig",&fNorig);
    fTree->SetBranchAddress("Ndecay",&fNdecay);
    fTree->SetBranchAddress("Ntype",&fNtype);
    fTree->SetBranchAddress("Vx",&fVx);
    fTree->SetBranchAddress("Vy",&fVy);
    fTree->SetBranchAddress("Vz",&fVz);
    fTree->SetBranchAddress("pdpx",&fPdpx);
    fTree->SetBranchAddress("pdpy",&fPdpy);
    fTree->SetBranchAddress("pdpz",&fPdpz);  
    fTree->SetBranchAddress("ppdxdz",&fPpdxdz);
    fTree->SetBranchAddress("ppdydz",&fPpdydz);
    fTree->SetBranchAddress("pppz",&fPppz);  
    fTree->SetBranchAddress("ppenergy",&fPpenergy);
    fTree->SetBranchAddress("ppmedium",&fPpmedium);
    fTree->SetBranchAddress("ptype",&fPtype);
    fTree->SetBranchAddress("ppvx",&fPpvx);
    fTree->SetBranchAddress("ppvy",&fPpvy);
    fTree->SetBranchAddress("ppvz",&fPpvz);
    fTree->SetBranchAddress("muparpx",&fMuparpx);
    fTree->SetBranchAddress("muparpy",&fMuparpy);
    fTree->SetBranchAddress("muparpz",&fMuparpz);
    fTree->SetBranchAddress("mupare",&fMupare);
    fTree->SetBranchAddress("Necm",&fNecm);
    fTree->SetBranchAddress("Nimpwt",&fNimpwt);
    fTree->SetBranchAddress("xpoint",&fXpoint);
    fTree->SetBranchAddress("ypoint",&fYpoint);
    fTree->SetBranchAddress("zpoint",&fZpoint);
  }

  fTree->SetBranchAddress("mu_tvx",&fMuX);
  fTree->SetBranchAddress("mu_tvy",&fMuY);
  fTree->SetBranchAddress("mu_tvz",&fMuZ);
  fTree->SetBranchAddress("mu_tpx",&fMuPX);
  fTree->SetBranchAddress("mu_tpy",&fMuPY);
  fTree->SetBranchAddress("mu_tpz",&fMuPZ);
  fTree->SetBranchAddress("mu_tptype",&fMuPID);
  fTree->SetBranchAddress("mu_tgen",&fMuGen);
  if(!fusemini){
    fTree->SetBranchAddress("mu_Ndxdz",&fMuNdxdz);
    fTree->SetBranchAddress("mu_Ndydz",&fMuNdydz);
    fTree->SetBranchAddress("mu_Npz",&fMuNpz);
    fTree->SetBranchAddress("mu_Nenergy",&fMuNenergy);
    fTree->SetBranchAddress("mu_Ndxdznea",&fMuNdxdznear);
    fTree->SetBranchAddress("mu_Ndydznea",&fMuNdydznear);
    fTree->SetBranchAddress("mu_Nenergyn",&fMuNenergynear);
    fTree->SetBranchAddress("mu_Nwtnear",&fMuNwtnear);
    fTree->SetBranchAddress("mu_Ndxdzfar",&fMuNdxdzfar);
    fTree->SetBranchAddress("mu_Ndydzfar",&fMuNdydzfar);
    fTree->SetBranchAddress("mu_Nenergyf",&fMuNenergyfar);
    fTree->SetBranchAddress("mu_Nwtfar",&fMuNwtfar);
    fTree->SetBranchAddress("mu_Norig",&fMuNorig);
    fTree->SetBranchAddress("mu_Ndecay",&fMuNdecay);
    fTree->SetBranchAddress("mu_Ntype",&fMuNtype);
    fTree->SetBranchAddress("mu_Vx",&fMuVx);
    fTree->SetBranchAddress("mu_Vy",&fMuVy);
    fTree->SetBranchAddress("mu_Vz",&fMuVz);
    fTree->SetBranchAddress("mu_pdpx",&fMuPdpx);
    fTree->SetBranchAddress("mu_pdpy",&fMuPdpy);
    fTree->SetBranchAddress("mu_pdpz",&fMuPdpz);  
    fTree->SetBranchAddress("mu_ppdxdz",&fMuPpdxdz);
    fTree->SetBranchAddress("mu_ppdydz",&fMuPpdydz);
    fTree->SetBranchAddress("mu_pppz",&fMuPppz);  
    fTree->SetBranchAddress("mu_ppenergy",&fMuPpenergy);
    fTree->SetBranchAddress("mu_ppmedium",&fMuPpmedium);
    fTree->SetBranchAddress("mu_ptype",&fMuPtype);
    fTree->SetBranchAddress("mu_ppvx",&fMuPpvx);
    fTree->SetBranchAddress("mu_ppvy",&fMuPpvy);
    fTree->SetBranchAddress("mu_ppvz",&fMuPpvz);
    fTree->SetBranchAddress("mu_muparpx",&fMuMuparpx);
    fTree->SetBranchAddress("mu_muparpy",&fMuMuparpy);
    fTree->SetBranchAddress("mu_muparpz",&fMuMuparpz);
    fTree->SetBranchAddress("mu_mupare",&fMuMupare);
    fTree->SetBranchAddress("mu_Necm",&fMuNecm);
    fTree->SetBranchAddress("mu_Nimpwt",&fMuNimpwt);
    fTree->SetBranchAddress("mu_xpoint",&fMuXpoint);
    fTree->SetBranchAddress("mu_ypoint",&fMuYpoint);
    fTree->SetBranchAddress("mu_zpoint",&fMuZpoint);
  }
}

//****************************
void MuParentHelper::SetVectorEntries()
{
  if(!fok) return;
  map<Int_t,map<Int_t,Int_t> >::iterator beg = fVec.begin();
  map<Int_t,map<Int_t,Int_t> >::iterator end = fVec.end();
  while(beg!=end){
    beg->second.clear();
  }
  fVec.clear();
  Int_t nentries = fTree->GetEntries();
  for(int i=0;i<nentries;i++){
    fTree->GetEntry(i);
    fVec[fRunNum][fEvtNum] = i;
  }
}

//****************************
void MuParentHelper::CloseFile()
{
  if(fFile) {
    fFile->Close();
    delete fFile;
  }
  fok = false;
}

//****************************
void MuParentHelper::GetMuParent(Int_t runnum, Int_t evtnum,
                                 Float_t mu_tpx, Float_t mu_tpy, 
                                 Float_t mu_tpz, NuParent &par)
{
  const std::map<Int_t,std::map<Int_t,Int_t> >::const_iterator run_it = fVec.find(runnum);
  if(run_it == fVec.end())
  {
    //     std::cerr << "MuParentHelper: failed to find run number " << runnum << endl;
     return;
  }

  const std::map<Int_t,Int_t> &evt_map = run_it -> second;
  const std::map<Int_t,Int_t>::const_iterator evt_it = evt_map.find(evtnum);
  if(evt_it == evt_map.end())
  {
    //     std::cerr << "MuParentHelper: failed to find event number " << evtnum << endl;
     return;
  }
  
  const int entry = evt_it -> second;

  fTree->GetEntry(entry);
  if(runnum==fRunNum && evtnum==fEvtNum &&
     2.*TMath::Abs(mu_tpx-fMuPX)/(mu_tpx+fMuPX)<1e-4 && 
     2.*TMath::Abs(mu_tpy-fMuPY)/(mu_tpy+fMuPY)<1e-4 && 
     2.*TMath::Abs(mu_tpz-fMuPZ)/(mu_tpz+fMuPZ)<1e-4) {    
    //have an entry in the tree:
    if(fPID!=0 && fGen!=0) {
      //found parent!
      //return parent info
      par.SetX(fX);
      par.SetY(fY);
      par.SetZ(fZ);
      par.SetPx(fPX);
      par.SetPy(fPY);
      par.SetPz(fPZ);
      par.SetPID(GEANTtoPDG(fPID));
      par.SetGen(fGen);
    }
    else {
      //couldn't find the parent
      //return muon parameters
      par.SetX(fMuX);
      par.SetY(fMuY);
      par.SetZ(fMuZ);
      par.SetPx(fMuPX);
      par.SetPy(fMuPY);
      par.SetPz(fMuPZ);
      par.SetPID(GEANTtoPDG(fMuPID));
      par.SetGen(fMuGen);
    }
  }
  else {
    //have no info about this event
    //return zeros
    par.SetX(0);
    par.SetY(0);
    par.SetZ(0);
    par.SetPx(0);
    par.SetPy(0);
    par.SetPz(0);
    par.SetPID(0);
    par.SetGen(0);
  }
}

//****************************
void MuParentHelper::GetMuParent(NtpMCRecord *rec,int mcevent,
                                 NuParent &par)
{
  TClonesArray& mcarray = *(rec->mc);
  NtpMCTruth *mc = static_cast<NtpMCTruth*>(mcarray[mcevent]);
  Int_t runnum = mc->flux.fluxrun;
  Int_t evtnum = mc->flux.fluxevtno;
  Float_t mu_tpx = mc->flux.tpx;
  Float_t mu_tpy = mc->flux.tpy;
  Float_t mu_tpz = mc->flux.tpz;
  this->GetMuParent(runnum,evtnum,mu_tpx,mu_tpy,mu_tpz,par);
}

//****************************
void MuParentHelper::GetMuParent(NtpStRecord *rec,int mcevent,
                                 NuParent &par)
{
  TClonesArray& mcarray = *(rec->mc);
  NtpMCTruth *mc = static_cast<NtpMCTruth*>(mcarray[mcevent]);
  Int_t runnum = mc->flux.fluxrun;
  Int_t evtnum = mc->flux.fluxevtno;
  Float_t mu_tpx = mc->flux.tpx;
  Float_t mu_tpy = mc->flux.tpy;
  Float_t mu_tpz = mc->flux.tpz;
  this->GetMuParent(runnum,evtnum,mu_tpx,mu_tpy,mu_tpz,par);
}

//****************************
void MuParentHelper::FillANtpTruth(NtpMCRecord *rec,int mcevent,
                                   ANtpTruthInfoBeam &ntpBeam)
{
  NuParent parent;
  this->GetMuParent(rec,mcevent,parent);
  ntpBeam.parentX   = parent.GetX();
  ntpBeam.parentY   = parent.GetY();
  ntpBeam.parentZ   = parent.GetZ();
  ntpBeam.parentPX  = parent.GetPx();
  ntpBeam.parentPY  = parent.GetPy();
  ntpBeam.parentPZ  = parent.GetPz();
  ntpBeam.parentPID = parent.GetPID();
  ntpBeam.parentGen = parent.GetGen();
}

//****************************
void MuParentHelper::FillANtpTruth(NtpStRecord *rec,int mcevent,
                                   ANtpTruthInfoBeam &ntpBeam)
{
  NuParent parent;
  this->GetMuParent(rec,mcevent,parent);
  ntpBeam.parentX   = parent.GetX();
  ntpBeam.parentY   = parent.GetY();
  ntpBeam.parentZ   = parent.GetZ();
  ntpBeam.parentPX  = parent.GetPx();
  ntpBeam.parentPY  = parent.GetPy();
  ntpBeam.parentPZ  = parent.GetPz();
  ntpBeam.parentPID = parent.GetPID();
  ntpBeam.parentGen = parent.GetGen();
}

//****************************
Int_t GEANTtoPDG(Int_t gCode)
{
  if(gCode==1)  return  22;
  if(gCode==2)  return -11;
  if(gCode==3)  return  11;
  if(gCode==4)  return  14;
  if(gCode==5)  return -13;
  if(gCode==6)  return  13;
  if(gCode==7)  return  111;
  if(gCode==8)  return  211;
  if(gCode==9)  return -211;
  if(gCode==10) return  130;
  if(gCode==11) return  321;
  if(gCode==12) return -321;
  if(gCode==13) return  2112;
  if(gCode==14) return  2212;
  if(gCode==15) return -2212;
  if(gCode==16) return  310;
  if(gCode==17) return  221;
  if(gCode==18) return  3122;
  if(gCode==19) return  3222;
  if(gCode==20) return  3212;
  if(gCode==21) return  3112;
  if(gCode==22) return  3322;
  if(gCode==23) return  3312;
  if(gCode==24) return  3332;
  if(gCode==25) return -2112;
  if(gCode==26) return -3122;
  if(gCode==27) return -3112;
  if(gCode==28) return -3212;
  if(gCode==29) return -3112;
  if(gCode==30) return -3322;
  if(gCode==31) return -3312;
  if(gCode==32) return -3332;
  if(gCode==33) return -15;
  if(gCode==34) return  15;
  if(gCode==35) return  411;
  if(gCode==36) return -411;
  if(gCode==37) return  421;
  if(gCode==38) return -421;
  if(gCode==39) return  431;
  if(gCode==40) return -431;
  if(gCode==41) return  4122;
  if(gCode==42) return  24;
  if(gCode==43) return -24;
  if(gCode==44) return  23;
  return 0;
}

//****************************
void MuParentHelper::GetMuParent(NtpMCFluxInfo *mufi, NtpMCFluxInfo &fi)
{
  GetMuParent(mufi->fluxrun,mufi->fluxevtno,mufi->tpx,mufi->tpy,mufi->tpz,fi);

}

//****************************
void MuParentHelper::GetMuParent(Int_t runnum, Int_t evtnum,
                                 Float_t mu_tpx, Float_t mu_tpy, 
                                 Float_t mu_tpz, NtpMCFluxInfo &fi)
{
  const std::map<Int_t,std::map<Int_t,Int_t> >::const_iterator run_it = fVec.find(runnum);
  if(run_it == fVec.end())
  {
    //     std::cerr << "MuParentHelper: failed to find run number " << runnum << endl;
     return;
  }

  const std::map<Int_t,Int_t> &evt_map = run_it -> second;
  const std::map<Int_t,Int_t>::const_iterator evt_it = evt_map.find(evtnum);
  if(evt_it == evt_map.end())
  {
    //     std::cerr << "MuParentHelper: failed to find event number " << evtnum << endl;
     return;
  }
  
  const int entry = evt_it -> second;

  fTree->GetEntry(entry);
  if(runnum==fRunNum && evtnum==fEvtNum &&
     2.*TMath::Abs(mu_tpx-fMuPX)/(mu_tpx+fMuPX)<1e-4 && 
     2.*TMath::Abs(mu_tpy-fMuPY)/(mu_tpy+fMuPY)<1e-4 && 
     2.*TMath::Abs(mu_tpz-fMuPZ)/(mu_tpz+fMuPZ)<1e-4) {    
    //have an entry in the tree:
    if(fPID!=0 && fGen!=0) {
      //found parent!
      //return parent info
      fi.index=-1;
      fi.fluxrun=runnum;
      fi.fluxevtno=evtnum;
      fi.ndxdz=fNdxdz;
      fi.ndydz=fNdydz;
      fi.npz=fNpz;
      fi.nenergy=fNenergy;
      fi.ndxdznear=fNdxdznear;
      fi.ndydznear=fNdydznear;
      fi.nenergynear=fNenergynear;
      fi.nwtnear=fNwtnear;
      fi.ndxdzfar=fNdxdzfar;
      fi.ndydzfar=fNdydzfar;
      fi.nenergyfar=fNenergyfar;
      fi.nwtfar=fNwtfar;
      fi.norig=fNorig;
      fi.ndecay=fNdecay;
      fi.ntype=fNtype;
      fi.vx=fVx;
      fi.vy=fVy;
      fi.vz=fVz;
      fi.pdpx=fPdpx;
      fi.pdpy=fPdpy;
      fi.pdpz=fPdpz;
      fi.ppdxdz=fPpdxdz;
      fi.ppdydz=fPpdydz;
      fi.pppz=fPppz;
      fi.ppenergy=fPpenergy;
      fi.ppmedium=fPpmedium;
      fi.ptype=fPtype;
      fi.ppvx=fPpvx;
      fi.ppvy=fPpvy;
      fi.ppvz=fPpvz;
      fi.muparpx=fMuparpx;
      fi.muparpy=fMuparpy;
      fi.muparpz=fMuparpz;
      fi.mupare=fMupare;
      fi.necm=fNecm;
      fi.nimpwt=fNimpwt;
      fi.xpoint=fXpoint;
      fi.ypoint=fYpoint;
      fi.zpoint=fZpoint;
      fi.tvx=fX;
      fi.tvy=fY;
      fi.tvz=fZ;
      fi.tpx=fPX;
      fi.tpy=fPY;
      fi.tpz=fPZ;
      fi.tptype=GEANTtoPDG(fPID);
      fi.tgen=fGen;
    }
    else {
      //couldn't find the parent
      //return muon parameters
      fi.index=-1;
      fi.fluxrun=runnum;
      fi.fluxevtno=evtnum;
      fi.ndxdz=fMuNdxdz;
      fi.ndydz=fMuNdydz;
      fi.npz=fMuNpz;
      fi.nenergy=fMuNenergy;
      fi.ndxdznear=fMuNdxdznear;
      fi.ndydznear=fMuNdydznear;
      fi.nenergynear=fMuNenergynear;
      fi.nwtnear=fMuNwtnear;
      fi.ndxdzfar=fMuNdxdzfar;
      fi.ndydzfar=fMuNdydzfar;
      fi.nenergyfar=fMuNenergyfar;
      fi.nwtfar=fMuNwtfar;
      fi.norig=fMuNorig;
      fi.ndecay=fMuNdecay;
      fi.ntype=fMuNtype;
      fi.vx=fMuVx;
      fi.vy=fMuVy;
      fi.vz=fMuVz;
      fi.pdpx=fMuPdpx;
      fi.pdpy=fMuPdpy;
      fi.pdpz=fMuPdpz;
      fi.ppdxdz=fMuPpdxdz;
      fi.ppdydz=fMuPpdydz;
      fi.pppz=fMuPppz;
      fi.ppenergy=fMuPpenergy;
      fi.ppmedium=fMuPpmedium;
      fi.ptype=fMuPtype;
      fi.ppvx=fMuPpvx;
      fi.ppvy=fMuPpvy;
      fi.ppvz=fMuPpvz;
      fi.muparpx=fMuMuparpx;
      fi.muparpy=fMuMuparpy;
      fi.muparpz=fMuMuparpz;
      fi.mupare=fMuMupare;
      fi.necm=fMuNecm;
      fi.nimpwt=fMuNimpwt;
      fi.xpoint=fMuXpoint;
      fi.ypoint=fMuYpoint;
      fi.zpoint=fMuZpoint;
      fi.tvx=fMuX;
      fi.tvy=fMuY;
      fi.tvz=fMuZ;
      fi.tpx=fMuPX;
      fi.tpy=fMuPY;
      fi.tpz=fMuPZ;
      fi.tptype=GEANTtoPDG(fMuPID);
      fi.tgen=fMuGen;
    }
  }
  else {
    //have no info about this event
    //return zeros
      fi.index=-1;
      fi.fluxrun=runnum;
      fi.fluxevtno=evtnum;
      fi.ndxdz=0.;
      fi.ndydz=0.;
      fi.npz=0.;
      fi.nenergy=0.;
      fi.ndxdznear=0.;
      fi.ndydznear=0.;
      fi.nenergynear=0.;
      fi.nwtnear=0.;
      fi.ndxdzfar=0.;
      fi.ndydzfar=0.;
      fi.nenergyfar=0.;
      fi.nwtfar=0.;
      fi.norig=0;
      fi.ndecay=0;
      fi.ntype=0;
      fi.vx=0.;
      fi.vy=0.;
      fi.vz=0.;
      fi.pdpx=0.;
      fi.pdpy=0.;
      fi.pdpz=0.;
      fi.ppdxdz=0.;
      fi.ppdydz=0.;
      fi.pppz=0.;
      fi.ppenergy=0.;
      fi.ppmedium=0;
      fi.ptype=0;
      fi.ppvx=0.;
      fi.ppvy=0.;
      fi.ppvz=0.;
      fi.muparpx=0.;
      fi.muparpy=0.;
      fi.muparpz=0.;
      fi.mupare=0.;
      fi.necm=0.;
      fi.nimpwt=0.;
      fi.xpoint=0.;
      fi.ypoint=0.;
      fi.zpoint=0.;
      fi.tvx=0.;
      fi.tvy=0.;
      fi.tvz=0.;
      fi.tpx=0.;
      fi.tpy=0.;
      fi.tpz=0.;
      fi.tptype=0;
      fi.tgen=0;
  }
}

//****************************
void MuParentHelper::GetMuDaughter(NtpMCFluxInfo *oldfi,
                                   NtpMCFluxInfo &newfi, int desirednu)
{
  Int_t runnum = oldfi->fluxrun;
  Int_t evtnum = oldfi->fluxevtno;
  Float_t tpx = oldfi->tpx;
  Float_t tpy = oldfi->tpy;
  Float_t tpz = oldfi->tpz;
  this->GetMuDaughter(runnum,evtnum,tpx,tpy,tpz,newfi,desirednu);
}

//****************************
void MuParentHelper::GetMuDaughter(Int_t runnum, Int_t evtnum,
                                 Float_t tpx, Float_t tpy, 
                                 Float_t tpz, NtpMCFluxInfo &fi, int desirednu)
{
  const std::map<Int_t,std::map<Int_t,Int_t> >::const_iterator run_it = fVec.find(runnum);
  if(run_it == fVec.end())
  {
    //     std::cerr << "MuParentHelper: failed to find run number " << runnum << endl;
     return;
  }

  const std::map<Int_t,Int_t> &evt_map = run_it -> second;
  const std::map<Int_t,Int_t>::const_iterator evt_it = evt_map.find(evtnum);
  if(evt_it == evt_map.end())
  {
    //     std::cerr << "MuParentHelper: failed to find event number " << evtnum << endl;
     return;
  }
  
  const int entry = evt_it -> second;

  bool foundit=false;

  fTree->GetEntry(entry);
  if(fMuNtype==desirednu&&runnum==fRunNum&&evtnum==fEvtNum){ 
    foundit=true;
  }
  else{
    fTree->GetEntry(entry-1);
    if(fMuNtype==desirednu&&runnum==fRunNum&&evtnum==fEvtNum){
      foundit=true;
    }
    else{
      fTree->GetEntry(entry+1);    
      if(fMuNtype==desirednu&&runnum==fRunNum&&evtnum==fEvtNum){
        foundit=true;
      }
      else{
        foundit=false;
      }
    }
  }

  if(runnum==fRunNum && evtnum==fEvtNum &&
     2.*TMath::Abs(tpx-fPX)/(tpx+fPX)<1e-4 && 
     2.*TMath::Abs(tpy-fPY)/(tpy+fPY)<1e-4 && 
     2.*TMath::Abs(tpz-fPZ)/(tpz+fPZ)<1e-4) {    
    //have an entry in the tree:
    if(fPID!=0 && fGen!=0 && foundit) {
      //found parent!
      //return parent info
      fi.index=-1;
      fi.fluxrun=runnum;
      fi.fluxevtno=evtnum;
      fi.ndxdz=fMuNdxdz;
      fi.ndydz=fMuNdydz;
      fi.npz=fMuNpz;
      fi.nenergy=fMuNenergy;
      fi.ndxdznear=fMuNdxdznear;
      fi.ndydznear=fMuNdydznear;
      fi.nenergynear=fMuNenergynear;
      fi.nwtnear=fMuNwtnear;
      fi.ndxdzfar=fMuNdxdzfar;
      fi.ndydzfar=fMuNdydzfar;
      fi.nenergyfar=fMuNenergyfar;
      fi.nwtfar=fMuNwtfar;
      fi.norig=fMuNorig;
      fi.ndecay=fMuNdecay;
      fi.ntype=fMuNtype;
      fi.vx=fMuVx;
      fi.vy=fMuVy;
      fi.vz=fMuVz;
      fi.pdpx=fMuPdpx;
      fi.pdpy=fMuPdpy;
      fi.pdpz=fMuPdpz;
      fi.ppdxdz=fMuPpdxdz;
      fi.ppdydz=fMuPpdydz;
      fi.pppz=fMuPppz;
      fi.ppenergy=fMuPpenergy;
      fi.ppmedium=fMuPpmedium;
      fi.ptype=fMuPtype;
      fi.ppvx=fMuPpvx;
      fi.ppvy=fMuPpvy;
      fi.ppvz=fMuPpvz;
      fi.muparpx=fMuMuparpx;
      fi.muparpy=fMuMuparpy;
      fi.muparpz=fMuMuparpz;
      fi.mupare=fMuMupare;
      fi.necm=fMuNecm;
      fi.nimpwt=fMuNimpwt;
      fi.xpoint=fMuXpoint;
      fi.ypoint=fMuYpoint;
      fi.zpoint=fMuZpoint;
      fi.tvx=fMuX;
      fi.tvy=fMuY;
      fi.tvz=fMuZ;
      fi.tpx=fMuPX;
      fi.tpy=fMuPY;
      fi.tpz=fMuPZ;
      fi.tptype=GEANTtoPDG(fPID);
      fi.tgen=fMuGen;
    }
    else {
      //couldn't find the daughter
      //return zeros
      fi.index=-1;
      fi.fluxrun=runnum;
      fi.fluxevtno=evtnum;
      fi.ndxdz=0.;
      fi.ndydz=0.;
      fi.npz=0.;
      fi.nenergy=0.;
      fi.ndxdznear=0.;
      fi.ndydznear=0.;
      fi.nenergynear=0.;
      fi.nwtnear=0.;
      fi.ndxdzfar=0.;
      fi.ndydzfar=0.;
      fi.nenergyfar=0.;
      fi.nwtfar=0.;
      fi.norig=0;
      fi.ndecay=0;
      fi.ntype=0;
      fi.vx=0.;
      fi.vy=0.;
      fi.vz=0.;
      fi.pdpx=0.;
      fi.pdpy=0.;
      fi.pdpz=0.;
      fi.ppdxdz=0.;
      fi.ppdydz=0.;
      fi.pppz=0.;
      fi.ppenergy=0.;
      fi.ppmedium=0;
      fi.ptype=0;
      fi.ppvx=0.;
      fi.ppvy=0.;
      fi.ppvz=0.;
      fi.muparpx=0.;
      fi.muparpy=0.;
      fi.muparpz=0.;
      fi.mupare=0.;
      fi.necm=0.;
      fi.nimpwt=0.;
      fi.xpoint=0.;
      fi.ypoint=0.;
      fi.zpoint=0.;
      fi.tvx=0.;
      fi.tvy=0.;
      fi.tvz=0.;
      fi.tpx=0.;
      fi.tpy=0.;
      fi.tpz=0.;
      fi.tptype=0;
      fi.tgen=0;
    }
  }
  else {
    //couldn't find the daughter
    //return zeros
    fi.index=-1;
    fi.fluxrun=runnum;
    fi.fluxevtno=evtnum;
    fi.ndxdz=0.;
    fi.ndydz=0.;
    fi.npz=0.;
    fi.nenergy=0.;
    fi.ndxdznear=0.;
    fi.ndydznear=0.;
    fi.nenergynear=0.;
    fi.nwtnear=0.;
    fi.ndxdzfar=0.;
    fi.ndydzfar=0.;
    fi.nenergyfar=0.;
    fi.nwtfar=0.;
    fi.norig=0;
    fi.ndecay=0;
    fi.ntype=0;
    fi.vx=0.;
    fi.vy=0.;
    fi.vz=0.;
    fi.pdpx=0.;
    fi.pdpy=0.;
    fi.pdpz=0.;
    fi.ppdxdz=0.;
    fi.ppdydz=0.;
    fi.pppz=0.;
    fi.ppenergy=0.;
    fi.ppmedium=0;
    fi.ptype=0;
    fi.ppvx=0.;
    fi.ppvy=0.;
    fi.ppvz=0.;
    fi.muparpx=0.;
    fi.muparpy=0.;
    fi.muparpz=0.;
    fi.mupare=0.;
    fi.necm=0.;
    fi.nimpwt=0.;
    fi.xpoint=0.;
    fi.ypoint=0.;
    fi.zpoint=0.;
    fi.tvx=0.;
    fi.tvy=0.;
    fi.tvz=0.;
    fi.tpx=0.;
    fi.tpy=0.;
    fi.tpz=0.;
    fi.tptype=0;
    fi.tgen=0;
  }
}

//****************************



#endif

---