inuke_rw.h

Go to the documentation of this file.

#ifndef INUKE_RW_H
#define INUKE_RW_H
#include <vector>
#include <ostream>


class NtpStRecord;
class NtpMCTruth;


namespace inuke_reweight {
  typedef enum inuke_particle{
    kPion,
    kNucleon,
    kOther
  } inuke_particle_t;

  class delta_fate {
    // simple C analogue of inuke_reweight.fi : fate_delta_pn
  public:
    float cex,elas,inel,abs,pp,npp,nnp,npnp,piprod;//npnp called 4n4p in NEUGEN
    delta_fate();
    void set_abs_fates(float f);// set all absorption fates together
    void set_abs_fates(){set_abs_fates(abs);}
    std::ostream& print(std::ostream&);
        
  };
  class delta_scale {
    // simple C analogue of inuke_reweight.fi: rw_*scale
  public:
    float rhonuc,ft,xsec;
    delta_scale();
    std::ostream& print(std::ostream&);
    
  };
  
  class parameter_set {
    
  public:
       
    inuke_reweight::delta_fate pi_fate;
    inuke_reweight::delta_fate pn_fate;
    inuke_reweight::delta_scale pi_scale;
    inuke_reweight::delta_scale pn_scale;

    std::ostream& print(std::ostream&);
    double get_par(int ipar);
  };
  
  class parameter_limits {
    
  public:
    parameter_limits();
    parameter_limits(const inuke_reweight::parameter_set& sigmas, 
                     float nsigma);
    inuke_reweight::parameter_set upper;
    inuke_reweight::parameter_set lower;

    std::ostream& print(std::ostream&);
      
  };
  

  // call inuke_reweight.F : change_fate_probs
  // changes fate probabilities in fates_rw (variables rw_*)
  // note: this doesn't change the defaults in neugen
  void change_fate_probs(const inuke_reweight::delta_fate& , 
                         const inuke_reweight::inuke_particle_t&);

  void change_inter_scales(const inuke_reweight::delta_scale&,
                           const inuke_reweight::inuke_particle_t&);
  
  // calculate the weight for whatever event is in stdhep
  // we call the following:
  //     call identify_hadronic_system(idx_had,icode_had)
  //     call summarise_prenuke(idx_had,icode_had)
  //     call fill_shower_ntuple
  //
  // common /shwr/ ends up filled with weights for each hadron
  // function returns the event weight.
  double calc_weights();
  
  // calculate the weights for the entire record
  // push them back, in order, onto the vector
  // optional verbose output
  void calc_weights(const NtpStRecord*, std::vector<double>&, 
                    bool verbose=false);

  void calc_weights(const NtpStRecord*, const NtpMCTruth*, 
                    const std::vector<inuke_reweight::parameter_set>& v,
                    std::vector<double>&, double& nucrad, double& wrad,
                    bool verbose=false);

  void 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);

  void test_calc_weights(const char* filename, int entry,
                         const inuke_reweight::parameter_set& p);
  
  // generate uncorrelated 1 sigma shifts
  // one parameter shifted +/- at a time.
  // sigmas holds the magnitude of the shifts
  void generate_1sigma_shifts(const inuke_reweight::parameter_set& sigmas,
                              std::vector<inuke_reweight::parameter_set>& v);
  
  void generate_uncor_shifts(const inuke_reweight::parameter_set& sigmas,
                             const inuke_reweight::parameter_limits& limits,
                             inuke_reweight::parameter_set& p);
  
  // utility functions
  // generate a vector of scale or fate parameter sets
  // from the shifts recorded in the input variable f
  // the sets are made by shifting each parameter in f by +/- 1 sigma
  void get_1sigma_shifts(const inuke_reweight::delta_scale& f,
                         std::vector<inuke_reweight::delta_scale>& v);
  
  void get_1sigma_shifts(const inuke_reweight::delta_fate& f,
                         std::vector<inuke_reweight::delta_fate>& v);

  void 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& v);

  
  void 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& v);

  // get the stdhep index of the hadronic system
  int get_idx_had();
  // get the radial position inside the nucleus
  double get_nucrad(int idx_had);

  // correction factor for uniform density to WS in iron
  double unif_to_ws_fe(double r);
  
}

#endif

---