fill_stdhep.cxx

Go to the documentation of this file.

#include "NeugenInterface/fill_stdhep.h" 
#include <iostream>
// include stdhep.h from $CERN_DIR/include (possibly)
// defines the stdhep common as an extern struct
#include "stdhep.h"

#include "TClonesArray.h"
#include "TFile.h"
#include "TTree.h"

#include "MCNtuple/NtpMCStdHep.h"
#include "MCNtuple/NtpMCTruth.h"
#include "StandardNtuple/NtpStRecord.h"
#include "StandardNtuple/NtpStRecord.h"

#include <cmath>

int inuke_reweight::fill_stdhep(const NtpMCTruth& t, const TClonesArray& v){
  const int ilow=t.stdhep[0];
  const int ihigh=t.stdhep[1];
  return fill_stdhep(v,ilow,ihigh);
}

int inuke_reweight::fill_stdhep(const TClonesArray& v, int ilow, int ihigh){
  // return number of entries filled
  // -1 in case of error
  // ilow is the starting index
  // ihigh is the ending index (not 1 past the end)
  if(ilow>ihigh || ilow<0 || (ihigh-1)>v.GetEntriesFast()){
    std::cerr<<"Index problem!! ilow ihigh entries : "
             <<ilow<<" "<<ihigh<<" "<<v.GetEntriesFast()<<std::endl;
    return -1;
  }
  int nfill=0;
  for(int i=ilow; i<=ihigh; i++){  
    const NtpMCStdHep* pp = dynamic_cast<const NtpMCStdHep* >(v[i]);
    if(!pp) return -1;  
    const NtpMCStdHep& p = *pp;
    int r=add_to_stdhep(p, ilow, v);
    if(r==-2) return -1;
    
    else if(r>-1) nfill++;
  }
  stdhep_to_nuclear_coordinates();
  truncate_geant_daughters();
  rotate_vhep();

  return nfill;
}

int inuke_reweight::add_to_stdhep(const NtpMCStdHep& p, int offset, const TClonesArray& v){
  
  // add p to stdhep, returning 
  // position in list 
  // -1 if the particle was ignored.
  // -2 in case of trouble
  
  if(p.IstHEP>200) return -1; // geant made these, or whacky 999 entry
  // see if this particle was created from one with IstHEP>200
  // this is why
  if( (p.parent[0]==p.parent[1]) && p.parent[0]>0 && p.parent[1]>0){
    bool dont_add=false;  
    for(int j=0; j<2; j++){          
      const int m=p.parent[j];
      if(m>=0 && m<v.GetEntriesFast()){
        const NtpMCStdHep& mommy = dynamic_cast<const NtpMCStdHep&>( *(v[m]));
        if(mommy.IstHEP>200) dont_add=true;
      }
    }
    if(dont_add) return -1;
  }
  const int i=hepevt_.nhep; // postion of new entry
  hepevt_.isthep[i]=p.IstHEP;
  hepevt_.idhep[i]=p.IdHEP;
  // remember, C convention for 2D array indexing is opposite of FORTRAN

  for(int j=0; j<2; j++){        
    hepevt_.jmohep[i][j]=p.parent[j];
    hepevt_.jdahep[i][j]=p.child[j];
    // make indices start at 1
    if(p.parent[j]!=-1) hepevt_.jmohep[i][j] -=(offset-1);
    if(p.child[j]!=-1) hepevt_.jdahep[i][j]-=(offset-1);
    // stdhep package checks for stable particles
    // by looking for jdahep = 0 (not -1 and not a test on isthep!!)
    // translate -1 to 0 here
    if(hepevt_.jmohep[i][j]== -1) hepevt_.jmohep[i][j]=0;
    if(hepevt_.jdahep[i][j]== -1) hepevt_.jdahep[i][j]=0;
  }
  for(int j=0; j<4; j++){
    hepevt_.phep[i][j] = p.p4[j];
    hepevt_.vhep[i][j] = p.vtx[j];
  }
  double psq=0.0;
  for(int j=0; j<3; j++) psq+=pow(p.p4[j],2);
  hepevt_.phep[i][4] =  pow(p.p4[3],2) - psq; 
  // precision errors a problem here
  // we care about > 1 MeV
  const double epsilon = pow(1e-3,2);
  if(fabs(hepevt_.phep[i][4])<epsilon) hepevt_.phep[i][4]=0;
  // geantinos not on shell
  if(hepevt_.idhep[i]==28) hepevt_.phep[i][4]=0;  
  // mass zero particles
  if(abs(hepevt_.idhep[i])==14 || abs(hepevt_.idhep[i])==12 
     || abs(hepevt_.idhep[i])==16 
     || abs(hepevt_.idhep[i])==22) hepevt_.phep[i][4]=0;
  // by hand for electrons (rounding problems)
  if(abs(hepevt_.idhep[i])==11)    hepevt_.phep[i][4] = 0.000502;

  if( hepevt_.phep[i][4] < 0.0) {
      std::cerr<<"phep["<<i<<"][4] = "<<hepevt_.phep[i][4]<<" < 0!!"<<std::endl;
      p.Print(std::cerr);
      inuke_reweight::print_stdhep(2);
      return -2;
      
  }
  else hepevt_.phep[i][4] = sqrt(hepevt_.phep[i][4]);
  hepevt_.nhep++;
  
  return i;

  
}

void inuke_reweight::clear_stdhep(bool hardclear){
  if(!hardclear) { hepevt_.nhep=0; return; }
  // hardclear erases the existing block
  for(int i=0; i<hepevt_.nhep; i++){
    hepevt_.isthep[i]=0;
    hepevt_.idhep[i]=0;
    for(int j=0; j<2; j++){
      hepevt_.jmohep[i][j]= -1;
      hepevt_.jdahep[i][j]= -1;
    }
    for(int j=0; j<4; j++){
      hepevt_.phep[i][j] = 0;
      hepevt_.vhep[i][j] = 0;      
    }
      hepevt_.phep[i][4] = 0;
  }
  hepevt_.nhep=0;
}


//
extern "C" void heplst_(int* iopt);

void inuke_reweight::print_stdhep(int iopt) {
  heplst_(&iopt);
}


void inuke_reweight::test_fill_stdhep(const char* filename, int entry){

  TFile f(filename);
  TTree* t = (TTree*) f.Get("NtpSt");
  NtpStRecord* rec=0;
  t->SetBranchAddress("NtpStRecord",&rec);
  t->GetEntry(entry);

  const TClonesArray& mc_array= *(rec->mc);
  const TClonesArray& stdheps = *(rec->stdhep);

  for(int i=0; i<stdheps.GetEntries(); i++){
    const NtpMCStdHep* s = dynamic_cast<const NtpMCStdHep*>(stdheps[i]);
    s->Print();
  }
  
  int nmc = rec->mchdr.nmc;
  for(int i=0; i<nmc; i++){
    const NtpMCTruth* mct = dynamic_cast<const NtpMCTruth*>(mc_array[i]);
    mct->Print();
    int nfill = inuke_reweight::fill_stdhep(*mct, stdheps);
    std::cout<<"MC Event: "<<i<<"  filled "<<nfill<<" entries"<<std::endl;
    inuke_reweight::print_stdhep(2);
    inuke_reweight::clear_stdhep();
  }
  return;
}

void inuke_reweight::stdhep_to_nuclear_coordinates(){

  for(int i=1; i<hepevt_.nhep; i++){
    if( hepevt_.isthep[i]!=14 
        && hepevt_.isthep[i]!=11 && hepevt_.isthep[i]!=3){
      // remove external vertex position
      // preserve any relative changes
      for(int j=0; j<4; j++){
        hepevt_.vhep[i][j] -= hepevt_.vhep[0][j];
      }
    }
    else{
      for(int j=0; j<3; j++){
        // convert to fm
        hepevt_.vhep[i][j] *= 1e15;
      }
      hepevt_.vhep[i][3]=0.0;
    }
  }
  // zero first entry
  for(int j=0; j<4; j++){
    hepevt_.vhep[0][j] = 0;
  }


}


void inuke_reweight::truncate_geant_daughters(){

  for(int i=0; i<hepevt_.nhep; i++){

    if(hepevt_.jdahep[i][0] > hepevt_.nhep) hepevt_.jdahep[i][0]=0;
    if(hepevt_.jdahep[i][1] > hepevt_.nhep) hepevt_.jdahep[i][1]=0;
  }

}

void inuke_reweight::rotate_vhep(){
  // neugen generates neutrinos oriented along z
  // neugen_kine.F then rotates phep[][] along the direction
  // indicated by gnumi but neglects to rotate vhep[][]
  // this isn't a problem for isthep=1 since those are assigned
  // global coordinates
  //
  // so: rotate vhep to agree with neutrino direction
  // 
  // this is a copy of geant's gdrot
  //
  
  const double costh = hepevt_.phep[0][2]/hepevt_.phep[0][3];
  const double sinth = sqrt(1.0-costh*costh);
  const double phi=atan2(hepevt_.phep[0][1],hepevt_.phep[0][0]);
  const double cosph = cos(phi);
  const double sinph = sin(phi);
  
  for(int i=1; i<hepevt_.nhep; i++){
    double v[3]={hepevt_.vhep[i][0],hepevt_.vhep[i][1],hepevt_.vhep[i][2]};
    hepevt_.vhep[i][0]=v[0]*costh*cosph - v[1]*sinph + v[2]*sinth*cosph;
    hepevt_.vhep[i][1]=v[0]*costh*sinph + v[1]*cosph + v[2]*sinth*sinph;
    hepevt_.vhep[i][2]=-v[0]*sinth + v[2]*costh;

  }
  
  
}

---