inuke_reweight Namespace Reference

Classes

class  delta_fate
class  delta_scale
class  parameter_set
class  parameter_limits

Typedefs

typedef enum
inuke_reweight::inuke_particle 
inuke_particle_t

Enumerations

enum  inuke_particle { kPion, kNucleon, kOther }

Functions

int fill_stdhep (const NtpMCTruth &t, const TClonesArray &stdheps)
int fill_stdhep (const TClonesArray &stdheps, int ilow, int ihigh)
int add_to_stdhep (const NtpMCStdHep &p, int offset, const TClonesArray &stdheps)
void clear_stdhep (bool hardclear=false)
void print_stdhep (int iopt)
void test_fill_stdhep (const char *filename, int entry=0)
void stdhep_to_nuclear_coordinates ()
void truncate_geant_daughters ()
void rotate_vhep ()
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 &)
double calc_weights ()
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)
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)
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)
int get_idx_had ()
double get_nucrad (int idx_had)
double unif_to_ws_fe (double r)

Typedef Documentation


Enumeration Type Documentation

Enumerator:
kPion 
kNucleon 
kOther 

Definition at line 12 of file inuke_rw.h.

00012                              {
00013     kPion,
00014     kNucleon,
00015     kOther
00016   } inuke_particle_t;


Function Documentation

int inuke_reweight::add_to_stdhep ( const NtpMCStdHep p,
int  offset,
const TClonesArray &  stdheps 
)

Definition at line 51 of file fill_stdhep.cxx.

References NtpMCStdHep::child, MuELoss::e, hepevt_, NtpMCStdHep::IdHEP, hepevt::idhep, hepevt::isthep, NtpMCStdHep::IstHEP, hepevt::jdahep, hepevt::jmohep, Munits::m, hepevt::nhep, NtpMCStdHep::p4, NtpMCStdHep::parent, hepevt::phep, NtpMCStdHep::Print(), print_stdhep(), hepevt::vhep, and NtpMCStdHep::vtx.

Referenced by fill_stdhep().

00051                                                                                         {
00052   
00053   // add p to stdhep, returning 
00054   // position in list 
00055   // -1 if the particle was ignored.
00056   // -2 in case of trouble
00057   
00058   if(p.IstHEP>200) return -1; // geant made these, or whacky 999 entry
00059   // see if this particle was created from one with IstHEP>200
00060   // this is why
00061   if( (p.parent[0]==p.parent[1]) && p.parent[0]>0 && p.parent[1]>0){
00062     bool dont_add=false;  
00063     for(int j=0; j<2; j++){          
00064       const int m=p.parent[j];
00065       if(m>=0 && m<v.GetEntriesFast()){
00066         const NtpMCStdHep& mommy = dynamic_cast<const NtpMCStdHep&>( *(v[m]));
00067         if(mommy.IstHEP>200) dont_add=true;
00068       }
00069     }
00070     if(dont_add) return -1;
00071   }
00072   const int i=hepevt_.nhep; // postion of new entry
00073   hepevt_.isthep[i]=p.IstHEP;
00074   hepevt_.idhep[i]=p.IdHEP;
00075   // remember, C convention for 2D array indexing is opposite of FORTRAN
00076 
00077   for(int j=0; j<2; j++){        
00078     hepevt_.jmohep[i][j]=p.parent[j];
00079     hepevt_.jdahep[i][j]=p.child[j];
00080     // make indices start at 1
00081     if(p.parent[j]!=-1) hepevt_.jmohep[i][j] -=(offset-1);
00082     if(p.child[j]!=-1) hepevt_.jdahep[i][j]-=(offset-1);
00083     // stdhep package checks for stable particles
00084     // by looking for jdahep = 0 (not -1 and not a test on isthep!!)
00085     // translate -1 to 0 here
00086     if(hepevt_.jmohep[i][j]== -1) hepevt_.jmohep[i][j]=0;
00087     if(hepevt_.jdahep[i][j]== -1) hepevt_.jdahep[i][j]=0;
00088   }
00089   for(int j=0; j<4; j++){
00090     hepevt_.phep[i][j] = p.p4[j];
00091     hepevt_.vhep[i][j] = p.vtx[j];
00092   }
00093   double psq=0.0;
00094   for(int j=0; j<3; j++) psq+=pow(p.p4[j],2);
00095   hepevt_.phep[i][4] =  pow(p.p4[3],2) - psq; 
00096   // precision errors a problem here
00097   // we care about > 1 MeV
00098   const double epsilon = pow(1e-3,2);
00099   if(fabs(hepevt_.phep[i][4])<epsilon) hepevt_.phep[i][4]=0;
00100   // geantinos not on shell
00101   if(hepevt_.idhep[i]==28) hepevt_.phep[i][4]=0;  
00102   // mass zero particles
00103   if(abs(hepevt_.idhep[i])==14 || abs(hepevt_.idhep[i])==12 
00104      || abs(hepevt_.idhep[i])==16 
00105      || abs(hepevt_.idhep[i])==22) hepevt_.phep[i][4]=0;
00106   // by hand for electrons (rounding problems)
00107   if(abs(hepevt_.idhep[i])==11)    hepevt_.phep[i][4] = 0.000502;
00108 
00109   if( hepevt_.phep[i][4] < 0.0) {
00110       std::cerr<<"phep["<<i<<"][4] = "<<hepevt_.phep[i][4]<<" < 0!!"<<std::endl;
00111       p.Print(std::cerr);
00112       inuke_reweight::print_stdhep(2);
00113       return -2;
00114       
00115   }
00116   else hepevt_.phep[i][4] = sqrt(hepevt_.phep[i][4]);
00117   hepevt_.nhep++;
00118   
00119   return i;
00120 
00121   
00122 }

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

Definition at line 231 of file inuke_rw.cxx.

References calc_weights(), change_fate_probs(), change_inter_scales(), clear_stdhep(), fill_stdhep(), get_idx_had(), get_nucrad(), kNucleon, kPion, inuke_reweight::parameter_set::pi_fate, inuke_reweight::parameter_set::pi_scale, inuke_reweight::parameter_set::pn_fate, inuke_reweight::parameter_set::pn_scale, NtpMCTruth::Print(), print_stdhep(), NtpStRecord::stdhep, unif_to_ws_fe(), and NtpMCTruth::z.

00234               {
00235   if(!rec) return;
00236   if(!mct) return;
00237   if(psets.size()==0) return;
00238 
00239   // assure w is the size of psets
00240   weights.resize(psets.size());
00241   // loop over all parameter sets
00242   for (unsigned int i=0; i<psets.size(); i++){
00243     const parameter_set& p=psets[i];
00244     inuke_reweight::change_fate_probs(p.pi_fate,inuke_reweight::kPion);
00245     inuke_reweight::change_inter_scales(p.pi_scale,inuke_reweight::kPion);
00246     inuke_reweight::change_fate_probs(p.pn_fate,inuke_reweight::kNucleon);
00247     inuke_reweight::change_inter_scales(p.pn_scale,inuke_reweight::kNucleon);
00248 
00249     const TClonesArray& stdheps = *(rec->stdhep);
00250 
00251     if(verbose) mct->Print();
00252     int nfill = inuke_reweight::fill_stdhep(*mct, stdheps);
00253     if(verbose){
00254       std::cout<<"MC Event: "<<i<<"  filled "<<nfill<<" entries"<<std::endl;
00255       inuke_reweight::print_stdhep(2);
00256     }
00257     double w=inuke_reweight::calc_weights(); // for current stdhep
00258     weights[i]=w;
00259     if(verbose) std::cout<<"weights["<<i<<"] = "<<w<<std::endl;
00260     nucrad=-1.0; wrad=1.0;
00261     if(mct->z == 26){ // iron only
00262       int idx_had=inuke_reweight::get_idx_had();
00263       nucrad=inuke_reweight::get_nucrad(idx_had);
00264       wrad=inuke_reweight::unif_to_ws_fe(nucrad);
00265     }
00266 
00267     
00268     inuke_reweight::clear_stdhep();
00269   }
00270   
00271 }

void inuke_reweight::calc_weights ( const NtpStRecord rec,
std::vector< double > &  v,
bool  verbose = false 
)

Definition at line 274 of file inuke_rw.cxx.

References calc_weights(), clear_stdhep(), fill_stdhep(), NtpStRecord::mc, NtpStRecord::mchdr, NtpMCSummary::nmc, NtpMCTruth::Print(), print_stdhep(), and NtpStRecord::stdhep.

00274                                                                                            {
00275   
00276   const TClonesArray& mc_array= *(rec->mc);
00277   const TClonesArray& stdheps = *(rec->stdhep);
00278   const int nmc = rec->mchdr.nmc;
00279   for(int i=0; i<nmc; i++){
00280     const NtpMCTruth* mct = dynamic_cast<const NtpMCTruth*>(mc_array[i]);
00281     if(verbose) mct->Print();
00282 
00283     int nfill = inuke_reweight::fill_stdhep(*mct, stdheps);
00284     if(verbose){
00285       std::cout<<"MC Event: "<<i<<"  filled "<<nfill<<" entries"<<std::endl;
00286       inuke_reweight::print_stdhep(2);
00287     }
00288     double w=inuke_reweight::calc_weights(); // for current stdhep
00289     v.push_back(w);
00290     std::cout<<"Weight = "<<w<<std::endl;
00291     
00292     inuke_reweight::clear_stdhep();
00293   }
00294   return; 
00295 
00296 }

double inuke_reweight::calc_weights (  ) 

Definition at line 211 of file inuke_rw.cxx.

References fill_shower_ntuple_(), identify_hadronic_system_(), shwr_, summarise_prenuke_(), shwr::tot_w_fate, and shwr::tot_w_inter.

Referenced by calc_weights(), MadMKAnalysis::CreatePAN(), NuIntranuke::Reweight(), and test_calc_weights().

00211                                    {
00212   int idx_had, icode_had;
00213   
00214   identify_hadronic_system_(&idx_had,&icode_had);
00215   if(icode_had==4){ return 1.0 ;} // nucleon target
00216   if(icode_had==0){ return 1.0 ;} // probably IMD
00217   //  inuke_reweight::print_stdhep(2);
00218   //  std::cout<<"idx_had icode_had: "<<idx_had<<" "<<icode_had<<std::endl;
00219   summarise_prenuke_(&idx_had,&icode_had);
00220   //  inuke_reweight::print_stdhep(2);
00221   fill_shower_ntuple_();
00222   //  inuke_reweight::print_stdhep(2);
00223   double W=shwr_.tot_w_fate * shwr_.tot_w_inter;
00224   return W;
00225  
00226 }

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

Definition at line 161 of file inuke_rw.cxx.

References inuke_reweight::delta_fate::abs, inuke_reweight::delta_fate::cex, change_fate_probs_(), inuke_reweight::delta_fate::elas, inuke_reweight::delta_fate::inel, kNucleon, kPion, inuke_reweight::delta_fate::nnp, inuke_reweight::delta_fate::npnp, inuke_reweight::delta_fate::npp, inuke_reweight::delta_fate::piprod, and inuke_reweight::delta_fate::pp.

Referenced by calc_weights(), and test_calc_weights().

00162                                                                                 {
00163 
00164   
00165   float v[10];
00166   v[0]=p.cex;  v[1]=p.elas;  v[2]=p.inel;
00167   v[3]=p.abs;  v[4]=p.pp;  v[5]=p.npp;
00168   v[6]=p.nnp;  v[7]=p.npnp;  v[8]=p.piprod;
00169   v[9]=0.0;
00170 
00171   int i=-1;
00172   if(t == inuke_reweight::kPion) i=0;
00173   else if( t == inuke_reweight::kNucleon ) i=1;
00174   //  for(int i=0; i<10; i++) std::cout<<"v["<<i<<"] = "<<v[i]<<std::endl;
00175   change_fate_probs_(v,&i);
00176   
00177 }

void inuke_reweight::change_inter_scales ( const inuke_reweight::delta_scale p,
const inuke_reweight::inuke_particle_t t 
)
void inuke_reweight::clear_stdhep ( bool  hardclear = false  ) 

Definition at line 124 of file fill_stdhep.cxx.

References hepevt_, hepevt::idhep, hepevt::isthep, hepevt::jdahep, hepevt::jmohep, hepevt::nhep, hepevt::phep, and hepevt::vhep.

Referenced by calc_weights(), and test_fill_stdhep().

00124                                                {
00125   if(!hardclear) { hepevt_.nhep=0; return; }
00126   // hardclear erases the existing block
00127   for(int i=0; i<hepevt_.nhep; i++){
00128     hepevt_.isthep[i]=0;
00129     hepevt_.idhep[i]=0;
00130     for(int j=0; j<2; j++){
00131       hepevt_.jmohep[i][j]= -1;
00132       hepevt_.jdahep[i][j]= -1;
00133     }
00134     for(int j=0; j<4; j++){
00135       hepevt_.phep[i][j] = 0;
00136       hepevt_.vhep[i][j] = 0;      
00137     }
00138       hepevt_.phep[i][4] = 0;
00139   }
00140   hepevt_.nhep=0;
00141 }

int inuke_reweight::fill_stdhep ( const TClonesArray &  stdheps,
int  ilow,
int  ihigh 
)

Definition at line 24 of file fill_stdhep.cxx.

References add_to_stdhep(), rotate_vhep(), stdhep_to_nuclear_coordinates(), and truncate_geant_daughters().

00024                                                                          {
00025   // return number of entries filled
00026   // -1 in case of error
00027   // ilow is the starting index
00028   // ihigh is the ending index (not 1 past the end)
00029   if(ilow>ihigh || ilow<0 || (ihigh-1)>v.GetEntriesFast()){
00030     std::cerr<<"Index problem!! ilow ihigh entries : "
00031              <<ilow<<" "<<ihigh<<" "<<v.GetEntriesFast()<<std::endl;
00032     return -1;
00033   }
00034   int nfill=0;
00035   for(int i=ilow; i<=ihigh; i++){  
00036     const NtpMCStdHep* pp = dynamic_cast<const NtpMCStdHep* >(v[i]);
00037     if(!pp) return -1;  
00038     const NtpMCStdHep& p = *pp;
00039     int r=add_to_stdhep(p, ilow, v);
00040     if(r==-2) return -1;
00041     
00042     else if(r>-1) nfill++;
00043   }
00044   stdhep_to_nuclear_coordinates();
00045   truncate_geant_daughters();
00046   rotate_vhep();
00047 
00048   return nfill;
00049 }

int inuke_reweight::fill_stdhep ( const NtpMCTruth t,
const TClonesArray &  stdheps 
)

Definition at line 18 of file fill_stdhep.cxx.

References NtpMCTruth::stdhep.

Referenced by calc_weights(), and test_fill_stdhep().

00018                                                                          {
00019   const int ilow=t.stdhep[0];
00020   const int ihigh=t.stdhep[1];
00021   return fill_stdhep(v,ilow,ihigh);
00022 }

void inuke_reweight::generate_1sigma_shifts ( const inuke_reweight::parameter_set sigmas,
std::vector< inuke_reweight::parameter_set > &  v 
)

Definition at line 410 of file inuke_rw.cxx.

References get_1sigma_shifts(), inuke_reweight::parameter_set::pi_fate, inuke_reweight::parameter_set::pi_scale, inuke_reweight::parameter_set::pn_fate, and inuke_reweight::parameter_set::pn_scale.

Referenced by MadMKAnalysis::CreatePAN(), and NuIntranuke::InitReweight().

00411                                            {
00412   
00413   std::vector<inuke_reweight::delta_fate> w;
00414   get_1sigma_shifts(sigmas.pi_fate,w);
00415   for (unsigned int i=0; i<w.size(); i++){
00416     inuke_reweight::parameter_set p;
00417     p.pi_fate=w[i];
00418     v.push_back(p);
00419   }
00420   w.clear();
00421   get_1sigma_shifts(sigmas.pn_fate,w);
00422   for (unsigned int i=0; i<w.size(); i++){
00423     inuke_reweight::parameter_set p;
00424     p.pn_fate=w[i];
00425     v.push_back(p);
00426   }
00427 
00428   std::vector<inuke_reweight::delta_scale> y;
00429   get_1sigma_shifts(sigmas.pi_scale,y);
00430   for (unsigned int i=0; i<y.size(); i++){
00431     inuke_reweight::parameter_set p;
00432     p.pi_scale=y[i];
00433     v.push_back(p);
00434   }
00435   y.clear();
00436   get_1sigma_shifts(sigmas.pn_scale,y);
00437   for (unsigned int i=0; i<y.size(); i++){
00438     inuke_reweight::parameter_set p;
00439     p.pn_scale=y[i];
00440     v.push_back(p);
00441   }
00442   
00443   
00444 
00445 }

void inuke_reweight::generate_uncor_shifts ( const inuke_reweight::parameter_set sigmas,
const inuke_reweight::parameter_limits limits,
inuke_reweight::parameter_set p 
)
void inuke_reweight::get_1sigma_shifts ( const inuke_reweight::delta_fate f,
std::vector< inuke_reweight::delta_fate > &  v 
)

Definition at line 357 of file inuke_rw.cxx.

References inuke_reweight::delta_fate::abs, inuke_reweight::delta_fate::cex, inuke_reweight::delta_fate::elas, inuke_reweight::delta_fate::inel, inuke_reweight::delta_fate::nnp, inuke_reweight::delta_fate::npnp, inuke_reweight::delta_fate::npp, inuke_reweight::delta_fate::piprod, and inuke_reweight::delta_fate::pp.

00358                                                                        {
00359 
00360   if(f.cex!=0) {
00361     inuke_reweight::delta_fate x;
00362     x.cex=f.cex; v.push_back(x);
00363     x.cex*= -1; v.push_back(x);
00364   }
00365   if(f.elas!=0) {
00366     inuke_reweight::delta_fate x;
00367     x.elas=f.elas; v.push_back(x);
00368     x.elas*= -1; v.push_back(x);
00369   }
00370   if(f.inel!=0) {
00371     inuke_reweight::delta_fate x;
00372     x.inel=f.inel; v.push_back(x);
00373     x.inel*= -1; v.push_back(x);
00374   }
00375   if(f.abs!=0) {
00376     inuke_reweight::delta_fate x;
00377     x.abs=f.abs; v.push_back(x);
00378     x.abs*= -1; v.push_back(x);
00379   }
00380   if(f.pp!=0) {
00381     inuke_reweight::delta_fate x;
00382     x.pp=f.pp; v.push_back(x);
00383     x.pp*= -1; v.push_back(x);
00384   }
00385   if(f.npp!=0) {
00386     inuke_reweight::delta_fate x;
00387     x.npp=f.npp; v.push_back(x);
00388     x.npp*= -1; v.push_back(x);
00389   }
00390   if(f.nnp!=0) {
00391     inuke_reweight::delta_fate x;
00392     x.nnp=f.nnp; v.push_back(x);
00393     x.nnp*= -1; v.push_back(x);
00394   }
00395   if(f.npnp!=0) {
00396     inuke_reweight::delta_fate x;
00397     x.npnp=f.npnp; v.push_back(x);
00398     x.npnp*= -1; v.push_back(x);
00399   }
00400   if(f.piprod!=0) {
00401     inuke_reweight::delta_fate x;
00402     x.piprod=f.piprod; v.push_back(x);
00403     x.piprod*= -1; v.push_back(x);
00404   }
00405 
00406 
00407 }

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

Definition at line 335 of file inuke_rw.cxx.

References inuke_reweight::delta_scale::ft, inuke_reweight::delta_scale::rhonuc, and inuke_reweight::delta_scale::xsec.

Referenced by generate_1sigma_shifts().

00336                                                                                {
00337   
00338   if(f.rhonuc!=0){
00339     inuke_reweight::delta_scale x;
00340     x.rhonuc=f.rhonuc; v.push_back(x);
00341     x.rhonuc*=-1; v.push_back(x);    
00342   }
00343   if(f.ft!=0){
00344     inuke_reweight::delta_scale x;
00345     x.ft=f.ft; v.push_back(x);
00346     x.ft*=-1; v.push_back(x);    
00347   }
00348   if(f.xsec!=0){
00349     inuke_reweight::delta_scale x;
00350     x.xsec=f.xsec; v.push_back(x);
00351     x.xsec*=-1; v.push_back(x);    
00352   }
00353 
00354 }

int inuke_reweight::get_idx_had (  ) 

Definition at line 585 of file inuke_rw.cxx.

References identify_hadronic_system_().

Referenced by calc_weights().

00585                                {
00586   int idx_had, icode_had;
00587   identify_hadronic_system_(&idx_had,&icode_had);
00588   if(icode_had==4){ return -1 ;} // nucleon target
00589   if(icode_had==0){ return -1 ;} // probably IMD
00590   return idx_had;
00591   
00592 }

double inuke_reweight::get_nucrad ( int  idx_had  ) 

Definition at line 596 of file inuke_rw.cxx.

References hepevt_, hepevt::nhep, and hepevt::vhep.

Referenced by calc_weights().

00596                                             {
00597   double r=-1.0;
00598   int i=idx_had-1;
00599   if(i<0 || i>hepevt_.nhep) return r;
00600   r=0.0;
00601   for(int j=0; j<3; j++){
00602     r+=pow(hepevt_.vhep[i][j],2.0);
00603   }
00604   return sqrt(r);
00605 }

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 v 
)

Definition at line 476 of file inuke_rw.cxx.

References inuke_reweight::delta_fate::abs, inuke_reweight::delta_fate::cex, inuke_reweight::delta_fate::elas, in_limits(), inuke_reweight::delta_fate::inel, inuke_reweight::delta_fate::nnp, inuke_reweight::delta_fate::npnp, inuke_reweight::delta_fate::npp, inuke_reweight::delta_fate::piprod, and inuke_reweight::delta_fate::pp.

00479                                                                   {
00480 
00481   while(!in_limits(x.cex=gRandom->Gaus()*f.cex,llim.cex,ulim.cex));
00482   while(!in_limits(x.elas=gRandom->Gaus()*f.elas,llim.elas,ulim.elas));
00483   while(!in_limits(x.inel=gRandom->Gaus()*f.inel,llim.inel,ulim.inel));
00484   while(!in_limits(x.abs=gRandom->Gaus()*f.abs,llim.abs,ulim.abs));
00485   while(!in_limits(x.pp=gRandom->Gaus()*f.pp,llim.pp,ulim.pp));
00486   while(!in_limits(x.npp=gRandom->Gaus()*f.npp,llim.npp,ulim.npp));
00487   while(!in_limits(x.nnp=gRandom->Gaus()*f.nnp,llim.nnp,ulim.nnp));
00488   while(!in_limits(x.npnp=gRandom->Gaus()*f.npnp,llim.npnp,ulim.npnp));
00489   while(!in_limits(x.piprod=gRandom->Gaus()*f.piprod,llim.piprod,ulim.piprod));
00490 
00491   /*
00492   x.cex=gRandom->Gaus()*f.cex;
00493   x.elas=gRandom->Gaus()*f.elas;
00494   x.inel=gRandom->Gaus()*f.inel;
00495   x.abs=gRandom->Gaus()*f.abs;
00496   x.pp=gRandom->Gaus()*f.pp;
00497   x.npp=gRandom->Gaus()*f.npp;
00498   x.nnp=gRandom->Gaus()*f.nnp;
00499   x.npnp=gRandom->Gaus()*f.npnp;
00500   x.piprod=gRandom->Gaus()*f.piprod;
00501   */
00502 
00503 }

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 v 
)

Definition at line 457 of file inuke_rw.cxx.

References inuke_reweight::delta_scale::ft, in_limits(), inuke_reweight::delta_scale::rhonuc, and inuke_reweight::delta_scale::xsec.

Referenced by generate_uncor_shifts().

00460                                                                     {
00461   // keep calling till in_limits
00462   while(!in_limits(x.rhonuc=gRandom->Gaus()*f.rhonuc,llim.rhonuc,ulim.rhonuc));
00463   //  while(!in_limits(x.ft=gRandom->Gaus()*f.ft,llim.ft,ulim.ft));
00464   while(!in_limits(x.ft=gRandom->Uniform(-1,1)*f.ft,llim.ft,ulim.ft));
00465   while(!in_limits(x.xsec=gRandom->Gaus()*f.xsec,llim.xsec,ulim.xsec));
00466 
00467   /*
00468   x.rhonuc=gRandom->Gaus()*f.rhonuc;
00469   x.ft=gRandom->Gaus()*f.ft;
00470   x.xsec=gRandom->Gaus()*f.xsec;
00471   */
00472 
00473 }

void inuke_reweight::print_stdhep ( int  iopt  ) 

Definition at line 147 of file fill_stdhep.cxx.

References heplst_().

Referenced by add_to_stdhep(), calc_weights(), and test_fill_stdhep().

00147                                           {
00148   heplst_(&iopt);
00149 }

void inuke_reweight::rotate_vhep (  ) 

Definition at line 218 of file fill_stdhep.cxx.

References hepevt_, hepevt::nhep, hepevt::phep, and hepevt::vhep.

Referenced by fill_stdhep().

00218                                 {
00219   // neugen generates neutrinos oriented along z
00220   // neugen_kine.F then rotates phep[][] along the direction
00221   // indicated by gnumi but neglects to rotate vhep[][]
00222   // this isn't a problem for isthep=1 since those are assigned
00223   // global coordinates
00224   //
00225   // so: rotate vhep to agree with neutrino direction
00226   // 
00227   // this is a copy of geant's gdrot
00228   //
00229   
00230   const double costh = hepevt_.phep[0][2]/hepevt_.phep[0][3];
00231   const double sinth = sqrt(1.0-costh*costh);
00232   const double phi=atan2(hepevt_.phep[0][1],hepevt_.phep[0][0]);
00233   const double cosph = cos(phi);
00234   const double sinph = sin(phi);
00235   
00236   for(int i=1; i<hepevt_.nhep; i++){
00237     double v[3]={hepevt_.vhep[i][0],hepevt_.vhep[i][1],hepevt_.vhep[i][2]};
00238     hepevt_.vhep[i][0]=v[0]*costh*cosph - v[1]*sinph + v[2]*sinth*cosph;
00239     hepevt_.vhep[i][1]=v[0]*costh*sinph + v[1]*cosph + v[2]*sinth*sinph;
00240     hepevt_.vhep[i][2]=-v[0]*sinth + v[2]*costh;
00241 
00242   }
00243   
00244   
00245 }

void inuke_reweight::stdhep_to_nuclear_coordinates (  ) 

Definition at line 180 of file fill_stdhep.cxx.

References hepevt_, hepevt::isthep, hepevt::nhep, and hepevt::vhep.

Referenced by fill_stdhep().

00180                                                   {
00181 
00182   for(int i=1; i<hepevt_.nhep; i++){
00183     if( hepevt_.isthep[i]!=14 
00184         && hepevt_.isthep[i]!=11 && hepevt_.isthep[i]!=3){
00185       // remove external vertex position
00186       // preserve any relative changes
00187       for(int j=0; j<4; j++){
00188         hepevt_.vhep[i][j] -= hepevt_.vhep[0][j];
00189       }
00190     }
00191     else{
00192       for(int j=0; j<3; j++){
00193         // convert to fm
00194         hepevt_.vhep[i][j] *= 1e15;
00195       }
00196       hepevt_.vhep[i][3]=0.0;
00197     }
00198   }
00199   // zero first entry
00200   for(int j=0; j<4; j++){
00201     hepevt_.vhep[0][j] = 0;
00202   }
00203 
00204 
00205 }

void inuke_reweight::test_calc_weights ( const char *  filename,
int  entry,
const inuke_reweight::parameter_set p 
)
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 
)

Definition at line 306 of file inuke_rw.cxx.

References calc_weights(), change_fate_probs(), change_inter_scales(), kNucleon, and kPion.

Referenced by test_calc_weights().

00310                                            {
00311 
00312 
00313   TFile f(filename);
00314   TTree* t = (TTree*) f.Get("NtpSt");
00315   NtpStRecord* rec=0;
00316   t->SetBranchAddress("NtpStRecord",&rec);
00317   t->GetEntry(entry);
00318   
00319   change_fate_probs(pi_fate,inuke_reweight::kPion);
00320   change_fate_probs(pn_fate,inuke_reweight::kNucleon);
00321   
00322   change_inter_scales(pi_scale,inuke_reweight::kPion);
00323   change_inter_scales(pn_scale,inuke_reweight::kNucleon);
00324 
00325   std::vector<double> v;
00326   bool verbose=false;
00327   calc_weights(rec,v,verbose);
00328   for (unsigned int i=0; i<v.size(); i++){
00329     std::cout<<i<<" "<<v[i]<<std::endl;
00330   }
00331   
00332 
00333 }

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

Definition at line 152 of file fill_stdhep.cxx.

References clear_stdhep(), fill_stdhep(), NtpStRecord::mc, NtpStRecord::mchdr, NtpMCSummary::nmc, NtpMCStdHep::Print(), NtpMCTruth::Print(), print_stdhep(), and NtpStRecord::stdhep.

00152                                                                     {
00153 
00154   TFile f(filename);
00155   TTree* t = (TTree*) f.Get("NtpSt");
00156   NtpStRecord* rec=0;
00157   t->SetBranchAddress("NtpStRecord",&rec);
00158   t->GetEntry(entry);
00159 
00160   const TClonesArray& mc_array= *(rec->mc);
00161   const TClonesArray& stdheps = *(rec->stdhep);
00162 
00163   for(int i=0; i<stdheps.GetEntries(); i++){
00164     const NtpMCStdHep* s = dynamic_cast<const NtpMCStdHep*>(stdheps[i]);
00165     s->Print();
00166   }
00167   
00168   int nmc = rec->mchdr.nmc;
00169   for(int i=0; i<nmc; i++){
00170     const NtpMCTruth* mct = dynamic_cast<const NtpMCTruth*>(mc_array[i]);
00171     mct->Print();
00172     int nfill = inuke_reweight::fill_stdhep(*mct, stdheps);
00173     std::cout<<"MC Event: "<<i<<"  filled "<<nfill<<" entries"<<std::endl;
00174     inuke_reweight::print_stdhep(2);
00175     inuke_reweight::clear_stdhep();
00176   }
00177   return;
00178 }

void inuke_reweight::truncate_geant_daughters (  ) 

Definition at line 208 of file fill_stdhep.cxx.

References hepevt_, hepevt::jdahep, and hepevt::nhep.

Referenced by fill_stdhep().

00208                                              {
00209 
00210   for(int i=0; i<hepevt_.nhep; i++){
00211 
00212     if(hepevt_.jdahep[i][0] > hepevt_.nhep) hepevt_.jdahep[i][0]=0;
00213     if(hepevt_.jdahep[i][1] > hepevt_.nhep) hepevt_.jdahep[i][1]=0;
00214   }
00215 
00216 }

double inuke_reweight::unif_to_ws_fe ( double  r  ) 

Definition at line 575 of file inuke_rw.cxx.

Referenced by calc_weights().

00575                                             {
00576   double c=1.0;
00577   if(r<5.36) {
00578     //    c = (pow(r,2.0)*(0.498473)/(1.+exp((r-4.1)/0.56)) )/(pow(r,2.0)*0.244814);
00579     c = ( (0.498473)/(1.+exp((r-4.1)/0.56)) )/(0.244814);
00580   }
00581   return c;
00582 
00583 }


Generated on 17 Dec 2018 for loon by  doxygen 1.6.1