inuke_rw.cxx

Go to the documentation of this file.

#include "NeugenInterface/inuke_rw.h"
#include "NeugenInterface/fill_stdhep.h"
#include <iostream>
#include <iomanip>

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

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


inuke_reweight::delta_fate::delta_fate():cex(0.0),elas(0.0),inel(0.0),abs(0.0),pp(0.0),npp(0.0),nnp(0.0),npnp(0.0),piprod(0.0){
  /*
  std::cout<<"delta_fate()"<<std::endl;
  std::cout<<"cex : "<<cex<<std::endl;
  std::cout<<"elase : "<<elas<<std::endl;
  std::cout<<"inel : "<<inel<<std::endl;
  std::cout<<"abs : "<<abs<<std::endl;
  std::cout<<"pp : "<<pp<<std::endl;
  std::cout<<"npp : "<<npp<<std::endl;
  std::cout<<"nnp : "<<nnp<<std::endl;
  std::cout<<"npnp : "<<npnp<<std::endl;
  std::cout<<"piprod : "<<piprod<<std::endl;
  */

}

std::ostream& inuke_reweight::delta_fate::print(std::ostream& os){
  //  os.width(4);
  //  os.precision(3);
  os.fill(' ');
  os<<"delta_fate:\t"<<std::setw(4)<<std::setprecision(3)<< cex<<" "<<std::setw(4)<<std::setprecision(3)<<elas<<" "<<std::setw(4)<<std::setprecision(3)<<inel<<" "<<std::setw(4)<<std::setprecision(3)<<abs<<" "<<std::setw(4)<<std::setprecision(3)<<pp<<" "<<std::setw(4)<<std::setprecision(3)<<npp<<" "<<std::setw(4)<<std::setprecision(3)<<std::setw(4)<<std::setprecision(3)<<nnp<<" "<<std::setw(4)<<std::setprecision(3)<<npnp<<" "<<std::setw(4)<<std::setprecision(3)<<piprod;
  return os;
}

void inuke_reweight::delta_fate::set_abs_fates(float f){ 
  abs=f; pp=f; npp=f; nnp=f; npnp=f;
}

inuke_reweight::delta_scale::delta_scale():rhonuc(0.0),ft(0.0),xsec(0.0){}

std::ostream& inuke_reweight::delta_scale::print(std::ostream& os){
  //  os.width(4);
  //  os.precision(3);
  //  setw(4)<<setprecision(3)<<
  os<<"delta_scale:\t"<<  std::setw(4)<<std::setprecision(3)<<rhonuc<<" "<<  std::setw(4)<<std::setprecision(3)<<ft<<" "<<std::setw(4)<<std::setprecision(3)<<xsec;
  return os;
}

std::ostream& inuke_reweight::parameter_set::print(std::ostream& os){
  os<<"pi "; pi_fate.print(os); os<<std::endl;
  os<<"pi "; pi_scale.print(os); os<<std::endl;
  os<<"pn "; pn_fate.print(os); os<<std::endl;
  os<<"pn "; pn_scale.print(os); os<<std::endl;
  os<<std::endl;

  return os;
}

double inuke_reweight::parameter_set::get_par(int ipar){
  double x=0;
  switch(ipar){
  case 0:
    x=pi_fate.cex;
    break;
  case 1:
    x=pi_fate.elas;
    break;
  case 2:
    x=pi_fate.inel;
    break;
  case 3:
    x=pi_fate.abs;
    break;
  case 4:
    x=pi_fate.pp;
    break;
  case 5:
    x=pi_fate.npp;
    break;
  case 6:
    x=pi_fate.nnp;
    break;
  case 7:
    x=pi_fate.npnp;
    break;
  case 8:
    x=pi_fate.piprod;
    break;
    // pn fates
  case 9:
    x=pn_fate.cex;
    break;
  case 10:
    x=pn_fate.elas;
    break;
  case 11:
    x=pn_fate.inel;
    break;
  case 12:
    x=pn_fate.abs;
    break;
  case 13:
    x=pn_fate.pp;
    break;
  case 14:
    x=pn_fate.npp;
    break;
  case 15:
    x=pn_fate.nnp;
    break;
  case 16:
    x=pn_fate.npnp;
    break;
  case 17:
    x=pn_fate.piprod;
    break;
    // pi scales
  case 18:
    x=pi_scale.rhonuc;
    break;
  case 19:
    x=pi_scale.ft;
    break;
  case 20:
    x=pi_scale.xsec;
    break;

    // pi scales
  case 21:
    x=pn_scale.rhonuc;
    break;
  case 22:
    x=pn_scale.ft;
    break;
  case 23:
    x=pn_scale.xsec;
    break;
  default:
    x=0;
    break;
  }
  return x;
}


std::ostream& inuke_reweight::parameter_limits::print(std::ostream& os){
  os<<"Upper limits: \n";upper.print(os); os<<std::endl;
  os<<"Lower limits: \n";lower.print(os); os<<std::endl;

  return os;
}

extern "C" int change_fate_probs_(float* vmod, int* itype);

void inuke_reweight::change_fate_probs(const inuke_reweight::delta_fate& p, 
                                       const inuke_reweight::inuke_particle_t& t) {

  
  float v[10];
  v[0]=p.cex;  v[1]=p.elas;  v[2]=p.inel;
  v[3]=p.abs;  v[4]=p.pp;  v[5]=p.npp;
  v[6]=p.nnp;  v[7]=p.npnp;  v[8]=p.piprod;
  v[9]=0.0;

  int i=-1;
  if(t == inuke_reweight::kPion) i=0;
  else if( t == inuke_reweight::kNucleon ) i=1;
  //  for(int i=0; i<10; i++) std::cout<<"v["<<i<<"] = "<<v[i]<<std::endl;
  change_fate_probs_(v,&i);
  
}


extern struct interaction_rw {
  float rw_pi_ft_scale; float rw_pn_ft_scale;
  float rw_pi_xsec_scale; float rw_pn_xsec_scale;
  float rw_rhonuc_scale;

} interaction_rw_;


void inuke_reweight::change_inter_scales(const inuke_reweight::delta_scale& p,
                                         const inuke_reweight::inuke_particle_t&t){

  // p.xsec = 0.2 means raise cross section by 20%
  
  interaction_rw_.rw_rhonuc_scale = (1+p.rhonuc);

  if(t==inuke_reweight::kPion){
    interaction_rw_.rw_pi_ft_scale=(1+p.ft);
    interaction_rw_.rw_pi_xsec_scale=(1+p.xsec);
  }
  else if (t==inuke_reweight::kNucleon){
    interaction_rw_.rw_pn_ft_scale=(1+p.ft);
    interaction_rw_.rw_pn_xsec_scale=(1+p.xsec);
  }

}

extern "C" void identify_hadronic_system_(int*, int*);
extern "C" void summarise_prenuke_(int*, int*);
extern "C" void fill_shower_ntuple_();
#include "shower_ntuple.h"

double inuke_reweight::calc_weights(){
  int idx_had, icode_had;
  
  identify_hadronic_system_(&idx_had,&icode_had);
  if(icode_had==4){ return 1.0 ;} // nucleon target
  if(icode_had==0){ return 1.0 ;} // probably IMD
  //  inuke_reweight::print_stdhep(2);
  //  std::cout<<"idx_had icode_had: "<<idx_had<<" "<<icode_had<<std::endl;
  summarise_prenuke_(&idx_had,&icode_had);
  //  inuke_reweight::print_stdhep(2);
  fill_shower_ntuple_();
  //  inuke_reweight::print_stdhep(2);
  double W=shwr_.tot_w_fate * shwr_.tot_w_inter;
  return W;
 
}



void inuke_reweight::calc_weights
(const NtpStRecord* rec, const NtpMCTruth* mct,
 const std::vector<inuke_reweight::parameter_set>& psets,
 std::vector<double>& weights, double& nucrad, double& wrad, 
 bool verbose){
  if(!rec) return;
  if(!mct) return;
  if(psets.size()==0) return;

  // assure w is the size of psets
  weights.resize(psets.size());
  // loop over all parameter sets
  for (unsigned int i=0; i<psets.size(); i++){
    const parameter_set& p=psets[i];
    inuke_reweight::change_fate_probs(p.pi_fate,inuke_reweight::kPion);
    inuke_reweight::change_inter_scales(p.pi_scale,inuke_reweight::kPion);
    inuke_reweight::change_fate_probs(p.pn_fate,inuke_reweight::kNucleon);
    inuke_reweight::change_inter_scales(p.pn_scale,inuke_reweight::kNucleon);

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

    if(verbose) mct->Print();
    int nfill = inuke_reweight::fill_stdhep(*mct, stdheps);
    if(verbose){
      std::cout<<"MC Event: "<<i<<"  filled "<<nfill<<" entries"<<std::endl;
      inuke_reweight::print_stdhep(2);
    }
    double w=inuke_reweight::calc_weights(); // for current stdhep
    weights[i]=w;
    if(verbose) std::cout<<"weights["<<i<<"] = "<<w<<std::endl;
    nucrad=-1.0; wrad=1.0;
    if(mct->z == 26){ // iron only
      int idx_had=inuke_reweight::get_idx_had();
      nucrad=inuke_reweight::get_nucrad(idx_had);
      wrad=inuke_reweight::unif_to_ws_fe(nucrad);
    }

    
    inuke_reweight::clear_stdhep();
  }
  
}


void inuke_reweight::calc_weights(const NtpStRecord* rec, std::vector<double>& v, bool verbose){
  
  const TClonesArray& mc_array= *(rec->mc);
  const TClonesArray& stdheps = *(rec->stdhep);
  const int nmc = rec->mchdr.nmc;
  for(int i=0; i<nmc; i++){
    const NtpMCTruth* mct = dynamic_cast<const NtpMCTruth*>(mc_array[i]);
    if(verbose) mct->Print();

    int nfill = inuke_reweight::fill_stdhep(*mct, stdheps);
    if(verbose){
      std::cout<<"MC Event: "<<i<<"  filled "<<nfill<<" entries"<<std::endl;
      inuke_reweight::print_stdhep(2);
    }
    double w=inuke_reweight::calc_weights(); // for current stdhep
    v.push_back(w);
    std::cout<<"Weight = "<<w<<std::endl;
    
    inuke_reweight::clear_stdhep();
  }
  return; 

}

void inuke_reweight::test_calc_weights(const char* filename, int entry,
                                       const inuke_reweight::parameter_set& p){
  inuke_reweight::test_calc_weights(filename, entry,
                                    p.pi_fate,p.pi_scale,p.pn_fate,p.pn_scale);
}


void inuke_reweight::test_calc_weights
(const char* filename, int entry,
 const inuke_reweight::delta_fate& pi_fate,
 const inuke_reweight::delta_scale& pi_scale,
 const inuke_reweight::delta_fate& pn_fate,
 const inuke_reweight::delta_scale& pn_scale){


  TFile f(filename);
  TTree* t = (TTree*) f.Get("NtpSt");
  NtpStRecord* rec=0;
  t->SetBranchAddress("NtpStRecord",&rec);
  t->GetEntry(entry);
  
  change_fate_probs(pi_fate,inuke_reweight::kPion);
  change_fate_probs(pn_fate,inuke_reweight::kNucleon);
  
  change_inter_scales(pi_scale,inuke_reweight::kPion);
  change_inter_scales(pn_scale,inuke_reweight::kNucleon);

  std::vector<double> v;
  bool verbose=false;
  calc_weights(rec,v,verbose);
  for (unsigned int i=0; i<v.size(); i++){
    std::cout<<i<<" "<<v[i]<<std::endl;
  }
  

}

void inuke_reweight::get_1sigma_shifts(const inuke_reweight::delta_scale& f,
                                       std::vector<inuke_reweight::delta_scale>& v){
  
  if(f.rhonuc!=0){
    inuke_reweight::delta_scale x;
    x.rhonuc=f.rhonuc; v.push_back(x);
    x.rhonuc*=-1; v.push_back(x);    
  }
  if(f.ft!=0){
    inuke_reweight::delta_scale x;
    x.ft=f.ft; v.push_back(x);
    x.ft*=-1; v.push_back(x);    
  }
  if(f.xsec!=0){
    inuke_reweight::delta_scale x;
    x.xsec=f.xsec; v.push_back(x);
    x.xsec*=-1; v.push_back(x);    
  }

}


void inuke_reweight::get_1sigma_shifts(const inuke_reweight::delta_fate& f,
                                std::vector<inuke_reweight::delta_fate>& v){

  if(f.cex!=0) {
    inuke_reweight::delta_fate x;
    x.cex=f.cex; v.push_back(x);
    x.cex*= -1; v.push_back(x);
  }
  if(f.elas!=0) {
    inuke_reweight::delta_fate x;
    x.elas=f.elas; v.push_back(x);
    x.elas*= -1; v.push_back(x);
  }
  if(f.inel!=0) {
    inuke_reweight::delta_fate x;
    x.inel=f.inel; v.push_back(x);
    x.inel*= -1; v.push_back(x);
  }
  if(f.abs!=0) {
    inuke_reweight::delta_fate x;
    x.abs=f.abs; v.push_back(x);
    x.abs*= -1; v.push_back(x);
  }
  if(f.pp!=0) {
    inuke_reweight::delta_fate x;
    x.pp=f.pp; v.push_back(x);
    x.pp*= -1; v.push_back(x);
  }
  if(f.npp!=0) {
    inuke_reweight::delta_fate x;
    x.npp=f.npp; v.push_back(x);
    x.npp*= -1; v.push_back(x);
  }
  if(f.nnp!=0) {
    inuke_reweight::delta_fate x;
    x.nnp=f.nnp; v.push_back(x);
    x.nnp*= -1; v.push_back(x);
  }
  if(f.npnp!=0) {
    inuke_reweight::delta_fate x;
    x.npnp=f.npnp; v.push_back(x);
    x.npnp*= -1; v.push_back(x);
  }
  if(f.piprod!=0) {
    inuke_reweight::delta_fate x;
    x.piprod=f.piprod; v.push_back(x);
    x.piprod*= -1; v.push_back(x);
  }


}

void inuke_reweight::generate_1sigma_shifts
(const inuke_reweight::parameter_set& sigmas,
 std::vector<inuke_reweight::parameter_set>& v){
  
  std::vector<inuke_reweight::delta_fate> w;
  get_1sigma_shifts(sigmas.pi_fate,w);
  for (unsigned int i=0; i<w.size(); i++){
    inuke_reweight::parameter_set p;
    p.pi_fate=w[i];
    v.push_back(p);
  }
  w.clear();
  get_1sigma_shifts(sigmas.pn_fate,w);
  for (unsigned int i=0; i<w.size(); i++){
    inuke_reweight::parameter_set p;
    p.pn_fate=w[i];
    v.push_back(p);
  }

  std::vector<inuke_reweight::delta_scale> y;
  get_1sigma_shifts(sigmas.pi_scale,y);
  for (unsigned int i=0; i<y.size(); i++){
    inuke_reweight::parameter_set p;
    p.pi_scale=y[i];
    v.push_back(p);
  }
  y.clear();
  get_1sigma_shifts(sigmas.pn_scale,y);
  for (unsigned int i=0; i<y.size(); i++){
    inuke_reweight::parameter_set p;
    p.pn_scale=y[i];
    v.push_back(p);
  }
  
  

}


bool in_limits(float& x, float l, float h);
bool in_limits(float& x, float l, float h){
  if(l==0.0 && h==0.0) {x=0.0; return true;}
  else return (x<h && x>l);
}




void inuke_reweight::get_uncor_shifts(const inuke_reweight::delta_scale& f,
                                      const inuke_reweight::delta_scale& llim,
                                      const inuke_reweight::delta_scale& ulim,
                                       inuke_reweight::delta_scale& x){
  // keep calling till in_limits
  while(!in_limits(x.rhonuc=gRandom->Gaus()*f.rhonuc,llim.rhonuc,ulim.rhonuc));
  //  while(!in_limits(x.ft=gRandom->Gaus()*f.ft,llim.ft,ulim.ft));
  while(!in_limits(x.ft=gRandom->Uniform(-1,1)*f.ft,llim.ft,ulim.ft));
  while(!in_limits(x.xsec=gRandom->Gaus()*f.xsec,llim.xsec,ulim.xsec));

  /*
  x.rhonuc=gRandom->Gaus()*f.rhonuc;
  x.ft=gRandom->Gaus()*f.ft;
  x.xsec=gRandom->Gaus()*f.xsec;
  */

}


void inuke_reweight::get_uncor_shifts(const inuke_reweight::delta_fate& f,
                                      const inuke_reweight::delta_fate& llim,
                                      const inuke_reweight::delta_fate& ulim,
                                      inuke_reweight::delta_fate& x){

  while(!in_limits(x.cex=gRandom->Gaus()*f.cex,llim.cex,ulim.cex));
  while(!in_limits(x.elas=gRandom->Gaus()*f.elas,llim.elas,ulim.elas));
  while(!in_limits(x.inel=gRandom->Gaus()*f.inel,llim.inel,ulim.inel));
  while(!in_limits(x.abs=gRandom->Gaus()*f.abs,llim.abs,ulim.abs));
  while(!in_limits(x.pp=gRandom->Gaus()*f.pp,llim.pp,ulim.pp));
  while(!in_limits(x.npp=gRandom->Gaus()*f.npp,llim.npp,ulim.npp));
  while(!in_limits(x.nnp=gRandom->Gaus()*f.nnp,llim.nnp,ulim.nnp));
  while(!in_limits(x.npnp=gRandom->Gaus()*f.npnp,llim.npnp,ulim.npnp));
  while(!in_limits(x.piprod=gRandom->Gaus()*f.piprod,llim.piprod,ulim.piprod));

  /*
  x.cex=gRandom->Gaus()*f.cex;
  x.elas=gRandom->Gaus()*f.elas;
  x.inel=gRandom->Gaus()*f.inel;
  x.abs=gRandom->Gaus()*f.abs;
  x.pp=gRandom->Gaus()*f.pp;
  x.npp=gRandom->Gaus()*f.npp;
  x.nnp=gRandom->Gaus()*f.nnp;
  x.npnp=gRandom->Gaus()*f.npnp;
  x.piprod=gRandom->Gaus()*f.piprod;
  */

}


void inuke_reweight::generate_uncor_shifts
(const inuke_reweight::parameter_set& sigmas,
 const inuke_reweight::parameter_limits& limits,
 inuke_reweight::parameter_set& p){
  
  inuke_reweight::get_uncor_shifts(sigmas.pi_fate,
                                   limits.lower.pi_fate, 
                                   limits.upper.pi_fate,
                                   p.pi_fate);
  inuke_reweight::get_uncor_shifts(sigmas.pn_fate,
                                   limits.lower.pn_fate, 
                                   limits.upper.pn_fate,
                                   p.pn_fate);

  inuke_reweight::get_uncor_shifts(sigmas.pi_scale,
                                   limits.lower.pi_scale, 
                                   limits.upper.pi_scale,
                                   p.pi_scale);
  inuke_reweight::get_uncor_shifts(sigmas.pn_scale,
                                   limits.lower.pn_scale, 
                                   limits.upper.pn_scale,
                                   p.pn_scale);

}

inuke_reweight::parameter_limits::parameter_limits() : upper(),lower(){}

inuke_reweight::parameter_limits::parameter_limits
(const inuke_reweight::parameter_set& sigmas, float nsigma){
  
  upper.pi_fate=sigmas.pi_fate;
  upper.pn_fate=sigmas.pn_fate;
  lower.pi_fate=sigmas.pi_fate;
  lower.pn_fate=sigmas.pn_fate;

  inuke_reweight::delta_fate* v[4]={&(upper.pi_fate),&(lower.pi_fate),
                                    &(upper.pn_fate),&(lower.pn_fate)};
  for(int j=0; j<4; j++){
    inuke_reweight::delta_fate& f = *v[j];
    //pow(-1,j) --> +,-,+,-
    f.cex=pow(-1,j)*(f.cex*nsigma);
    f.elas=pow(-1,j)*(f.elas*nsigma);
    f.inel=pow(-1,j)*(f.inel*nsigma);
    f.abs=pow(-1,j)*(f.abs*nsigma);
    f.pp=pow(-1,j)*(f.pp*nsigma);
    f.npp=pow(-1,j)*(f.npp*nsigma);
    f.nnp=pow(-1,j)*(f.nnp*nsigma);
    f.npnp=pow(-1,j)*(f.npnp*nsigma);
    f.piprod=pow(-1,j)*(f.piprod*nsigma);
  }


  upper.pi_scale=sigmas.pi_scale;
  upper.pn_scale=sigmas.pn_scale;
  lower.pi_scale=sigmas.pi_scale;
  lower.pn_scale=sigmas.pn_scale;

  inuke_reweight::delta_scale* w[4]={&(upper.pi_scale),&(lower.pi_scale),
                                    &(upper.pn_scale),&(lower.pn_scale)};

  for(int j=0; j<4; j++){
    inuke_reweight::delta_scale& f = *w[j];
    f.rhonuc=pow(-1,j)*(f.rhonuc*nsigma);
    f.ft=pow(-1,j)*(f.ft*nsigma);
    f.xsec=pow(-1,j)*(f.xsec*nsigma);
    
  }
}

double inuke_reweight::unif_to_ws_fe(double r){
  double c=1.0;
  if(r<5.36) {
    //    c = (pow(r,2.0)*(0.498473)/(1.+exp((r-4.1)/0.56)) )/(pow(r,2.0)*0.244814);
    c = ( (0.498473)/(1.+exp((r-4.1)/0.56)) )/(0.244814);
  }
  return c;

}

int inuke_reweight::get_idx_had(){
  int idx_had, icode_had;
  identify_hadronic_system_(&idx_had,&icode_had);
  if(icode_had==4){ return -1 ;} // nucleon target
  if(icode_had==0){ return -1 ;} // probably IMD
  return idx_had;
  
}

#include "stdhep.h"

double inuke_reweight::get_nucrad(int idx_had){
  double r=-1.0;
  int i=idx_had-1;
  if(i<0 || i>hepevt_.nhep) return r;
  r=0.0;
  for(int j=0; j<3; j++){
    r+=pow(hepevt_.vhep[i][j],2.0);
  }
  return sqrt(r);
}

---