BeamEnergyCalculator Class Reference

#include <BeamEnergyCalculator.h>

Inheritance diagram for BeamEnergyCalculator:
WeightCalculator

List of all members.

Classes

struct  myhash

Public Member Functions

 BeamEnergyCalculator (Registry *stdconfig=0)
virtual void Config ()
virtual void ReweightConfigReset ()
virtual double GetWeight (Registry *eventinfo)
virtual double GetWeight (MCEventInfo *event, NuParent *parent)
virtual double GetWeight (BeamType::BeamType_t bt, Detector::Detector_t det, int nu_idhep, double nu_energy, double high_energy_limit=120.0, double low_energy_limit=low_energy_cut)

Private Types

typedef __gnu_cxx::hash_map
< std::string, TH1 *, myhash
CacheMap
typedef __gnu_cxx::hash_map
< std::string, TH1 *, myhash >
::iterator 
CacheMapIter

Private Member Functions

double DeriveWeight (TH1 *, double, double, double)
double DeriveInverseEWeight (double, double, double)
TH1 * GetHist (BeamType::BeamType_t bt, Detector::Detector_t det, int nu_idhep)

Static Private Member Functions

static void ConstructName (Detector::Detector_t det, int nu_idhep, std::string &s)

Private Attributes

TDirectory * fFluxDir
std::string fFluxFileName
bool fFluxFileChanged
Detector::Detector_t fDetector
BeamType::BeamType_t fBeam
double fLowRange
double fHighRange
CacheMap fCache

Static Private Attributes

static const double low_energy_cut = 1e-6

Detailed Description

Definition at line 40 of file BeamEnergyCalculator.h.


Member Typedef Documentation

typedef __gnu_cxx::hash_map<std::string, TH1*,myhash> BeamEnergyCalculator::CacheMap [private]

Definition at line 85 of file BeamEnergyCalculator.h.

typedef __gnu_cxx::hash_map<std::string, TH1*,myhash>::iterator BeamEnergyCalculator::CacheMapIter [private]

Definition at line 86 of file BeamEnergyCalculator.h.


Constructor & Destructor Documentation

BeamEnergyCalculator::BeamEnergyCalculator ( Registry stdconfig = 0  ) 

Definition at line 43 of file BeamEnergyCalculator.cxx.

References Config(), WeightCalculator::fStandardConfigChanged, and WeightCalculator::fWCname.

00044    : 
00045   WeightCalculator(stdconfig),
00046   fFluxDir(0),fFluxFileName(""),
00047   fFluxFileChanged(true),
00048   fDetector(Detector::kUnknown),
00049   fBeam(BeamType::kUnknown),
00050   fLowRange(low_energy_cut), fHighRange(120.0),fCache()
00051 {
00052   fWCname="BeamEnergyCalculator";
00053   fStandardConfigChanged=true;
00054   Config();
00055 }


Member Function Documentation

void BeamEnergyCalculator::Config ( void   )  [virtual]

Reimplemented from WeightCalculator.

Definition at line 58 of file BeamEnergyCalculator.cxx.

References fBeam, fDetector, fFluxFileChanged, fFluxFileName, fHighRange, fLowRange, WeightCalculator::fStandardConfig, WeightCalculator::fStandardConfigChanged, Registry::Get(), Registry::GetCharString(), Msg::kInfo, and MSG.

Referenced by BeamEnergyCalculator().

00058                                   {
00059   MSG("BEC", Msg::kInfo)<<"configuring..."<<endl;
00060 
00061   int tempi=-1;
00062   double tempd=0.0;
00063   std::string ff;
00064   if(fStandardConfigChanged){
00065     // this is a little clumsy
00066     const char* c=fStandardConfig->GetCharString("beam:flux_file");    
00067     if(c){
00068       ff=c;
00069       if(ff!=fFluxFileName){
00070         fFluxFileName=ff;
00071         fFluxFileChanged=true;
00072       }
00073     }
00074     /*    
00075     if(fStandardConfig->Get("beam:nu_id",tempi)){
00076       fNuID=tempi;
00077     }
00078     if(fStandardConfig->Get("beam:nu_energy",tempd)){
00079       fNuEnergy=tempd;
00080     }
00081     */
00082 
00083     if(fStandardConfig->Get("beam:beam_type",tempi)){
00084       fBeam=static_cast<BeamType::BeamType_t>(tempi);
00085     }
00086     if(fStandardConfig->Get("beam:detector",tempi)){
00087       fDetector=static_cast<Detector::Detector_t>(tempi);
00088     }    
00089     if(fStandardConfig->Get("beam:low_energy_limit",tempd)){
00090       fLowRange=tempd;
00091     }
00092     if(fStandardConfig->Get("beam:high_energy_limit",tempd)){
00093       fHighRange=tempd;
00094     }
00095     
00096     fStandardConfigChanged=0;
00097   }
00098 
00099 }

void BeamEnergyCalculator::ConstructName ( Detector::Detector_t  det,
int  nu_idhep,
std::string &  s 
) [static, private]

Definition at line 241 of file BeamEnergyCalculator.cxx.

References bfld::AsString().

Referenced by GetHist().

00242                                                                     {
00243   // given the detector and type of neutrino
00244   // construct (part of) a histogram name
00245   // note, adds to 's'
00246   
00247   using std::string;
00248   string nu;
00249   switch(nu_idhep){
00250   case 12: nu="nue"; break;
00251   case -12: nu="nuebar"; break;
00252   case 14: nu="numu"; break;
00253   case -14: nu="numubar"; break;
00254   case 16: nu="nutau"; break;
00255   case -16: nu="nutaubar"; break;
00256   default: nu="???"; break;
00257   }
00258   string dets=Detector::AsString(det);
00259   for(string::size_type i=0; i<dets.size(); i++){
00260     dets[i]=std::tolower(dets[i]);
00261   }
00262 
00263   s+=dets+'_'+nu;
00264 
00265 }

double BeamEnergyCalculator::DeriveInverseEWeight ( double  energy,
double  high,
double  low 
) [private]

Definition at line 231 of file BeamEnergyCalculator.cxx.

References low_energy_cut.

Referenced by GetWeight().

00232                                                                          {
00233   // 1/E flux over range [low,high]
00234   // range is necessary for proper normalization of weights!
00235   if(low<low_energy_cut) low = low_energy_cut;
00236   if(energy<low || energy>high) return 0;
00237   return (1/energy)/std::log(high/low);
00238 
00239 }

double BeamEnergyCalculator::DeriveWeight ( TH1 *  h,
double  energy,
double  high,
double  low 
) [private]

Definition at line 189 of file BeamEnergyCalculator.cxx.

References Msg::kFatal, low_energy_cut, and MSG.

Referenced by GetWeight().

00190                                                                   {
00191   // simply uses the binned histogram to estimate the weights
00192   // certainly, it's reliant on the quality of the histogram...
00193 
00194   if(energy>high || energy<low) return 0;
00195   if(!h) return 0;
00196   if(low<low_energy_cut) low=low_energy_cut;
00197   Double_t* integral=h->GetIntegral();
00198   if(!integral){ h->ComputeIntegral(); integral=h->GetIntegral();}
00199   Int_t binlow=h->FindBin(low);
00200   Int_t binhigh=h->FindBin(high);
00201   if(binlow<1 || binhigh>h->GetNbinsX()) return 0;
00202   Int_t bin = h->FindBin(energy);
00203   double weight = h->GetBinContent(bin);
00204   double norm = integral[binhigh]-integral[binlow-1];
00205   Stat_t stats[4];
00206   h->GetStats(stats);
00207   Stat_t sumw=stats[0];
00208   /*
00209   static bool first=true;
00210   if(first){
00211     std::cout<<"Histogram: "<<std::endl;
00212     h->Print("all");
00213     for(int i=1; i<=h->GetNbinsX(); i++){
00214       std::cout<<"integral["<<i<<"] = "<<integral[i]<<std::endl;
00215     }
00216     std::cout<<"sumw: "<<sumw<<std::endl;
00217     first=false;
00218   }
00219   */
00220   if(norm<=0) {
00221     MSG("BEC",Msg::kFatal)<<"normalization=0"<<endl;
00222     assert(0);
00223   }
00224 
00225   weight/=(norm*sumw);
00226 
00227   return weight;
00228 }

TH1 * BeamEnergyCalculator::GetHist ( BeamType::BeamType_t  bt,
Detector::Detector_t  det,
int  nu_idhep 
) [private]

Definition at line 267 of file BeamEnergyCalculator.cxx.

References BeamType::AsTag(), ConstructName(), fCache, fFluxDir, fFluxFileChanged, fFluxFileName, gSystem(), it, Msg::kFatal, Msg::kInfo, Msg::kWarning, MAXMSG, and MSG.

Referenced by GetWeight().

00269                                                 {
00270   // get a histogram from fFluxDir
00271   // histograms are taken from the file according to a naming convention
00272   // that encompasses the type of neutrino, Detector, and BeamType
00273 
00274   // open new file if needed
00275   if(fFluxFileChanged){
00276     
00277     if(fFluxDir){
00278       delete fFluxDir; fFluxDir=0;
00279     }
00280     fCache.clear();
00281     MSG("BEC",Msg::kInfo)<<"opening new flux file: "
00282                          <<fFluxFileName<<endl;
00283     // check that file exists
00284     if(gSystem->AccessPathName(fFluxFileName.c_str()) == kTRUE){
00285       // file *doesn't* exist... who established that convention!!
00286       MSG("BEC",Msg::kFatal)<<"cannot find flux file: "
00287                             <<fFluxFileName<<endl;
00288       assert(0); //die,die,die
00289     }
00290     else{
00291       fFluxDir = new TFile(fFluxFileName.c_str());
00292       MSG("BEC",Msg::kInfo)<<"successfully opened: "
00293                            <<fFluxFileName<<endl;
00294     }
00295     fFluxFileChanged=false;
00296   }
00297 
00298   // look for histogram in our cache
00299   // initially I just wanted to use the TFile for this purpose
00300   // however, if you do TH1* h = (TH1*) file->Get("some_dir/some_hist")
00301   // multiple times, the TFile code looks in the list of objects
00302   // for "some_dir/some_hist" after the second time,
00303   // but, of course, the actual thing in the list is called "some_hist"
00304   // so, it rereads the histogram from the file, which is real slow
00305   //  I couldn't find a workaround, so, hash_map..
00306   TH1* h = 0;
00307 
00308   // make the histogram name
00309   // beamtype tag is the directory name
00310   std::string s=BeamType::AsTag(bt);
00311   // pick a flux histogram
00312   // there are also 'evt' histograms but
00313   // they don't play well with 1/E spectrum.. must know xsection
00314   s+="/flux_";
00315   ConstructName(det,nu_idhep,s);
00316   // so, for LE beam, muon-neutrinos in the near detector
00317   // s=z_000/flux_near_numu
00318 
00319 
00320   CacheMapIter it = fCache.find(s);
00321   if(it!=fCache.end()){
00322     h=it->second;
00323   }
00324   else{
00325     // if not found, look in the file    
00326     //TH1* h = FindHist(fFluxDir,bt,det,nu_idhep);
00327     h = dynamic_cast<TH1*>(fFluxDir->Get(s.c_str()));
00328     if(h){
00329       //      fCache[s]=h; // add it!
00330       fCache.insert(std::pair<std::string,TH1*>(s,h));
00331     }
00332     else{
00333       MAXMSG("BEC",Msg::kWarning,10)
00334         <<"Could not find a flux histogram with the tag: "
00335         <<s<<"\n It's possible to continue but"
00336         <<" these events will not be weighted properly."<<endl;      
00337     }
00338   }
00339   return h;
00340 }

double BeamEnergyCalculator::GetWeight ( BeamType::BeamType_t  bt,
Detector::Detector_t  det,
int  nu_idhep,
double  nu_energy,
double  high_energy_limit = 120.0,
double  low_energy_limit = low_energy_cut 
) [virtual]

Definition at line 164 of file BeamEnergyCalculator.cxx.

References DeriveInverseEWeight(), DeriveWeight(), GetHist(), Msg::kFatal, BeamType::kInverseE, BeamType::kUnknown, and MSG.

00168                                                                {
00169   double w=0.0;
00170   if(bt==BeamType::kInverseE){
00171     // special case of 1/E flux
00172     w=DeriveInverseEWeight(nu_energy,high_energy_limit, low_energy_limit);
00173   }
00174   else if(bt==BeamType::kUnknown){
00175     MSG("BEC", Msg::kFatal)<<"Requested weight for unknown beam type"<<endl;
00176     assert(0);
00177   }
00178   else{
00179     TH1* h = GetHist(bt,det,nu_idhep);
00180     
00181     if(!h) {w=1.0; return w;}
00182     
00183     w=DeriveWeight(h,nu_energy,high_energy_limit,low_energy_limit);
00184   }
00185 
00186   return w;
00187 }

double BeamEnergyCalculator::GetWeight ( MCEventInfo event,
NuParent parent 
) [virtual]

Reimplemented from WeightCalculator.

Definition at line 156 of file BeamEnergyCalculator.cxx.

References Msg::kFatal, and MSG.

00157                                                               {
00158   MSG("BEC", Msg::kFatal)<<"Function Not Implemented Yet!"<<endl;
00159   assert(0);
00160   return 1;
00161 }

double BeamEnergyCalculator::GetWeight ( Registry eventinfo  )  [virtual]

Reimplemented from WeightCalculator.

Definition at line 106 of file BeamEnergyCalculator.cxx.

References det, fBeam, fDetector, fHighRange, fLowRange, Registry::Get(), Msg::kFatal, Detector::kUnknown, BeamType::kUnknown, and MSG.

00106                                                  {
00107   
00108   BeamType::BeamType_t bt=BeamType::kUnknown;
00109   Detector::Detector_t det = Detector::kUnknown;
00110   int nu_idhep=0;
00111   double nu_energy=0;
00112   double high=0;
00113   double low=0;
00114   double w=1.0;
00115 
00116   // these must be passed in
00117   bool bad=false;
00118   if(!r->Get("beam:nu_id",nu_idhep)) bad=true;
00119   if(!r->Get("beam:nu_energy",nu_energy)) bad=true;
00120 
00121   if(bad){
00122     MSG("BEC", Msg::kFatal)<<"neutrino id or energy not set"<<endl;
00123     assert(0);
00124   }
00125   // look for certain values set in the registry
00126   // if not found, use defaults
00127   int tmpi=0;
00128   if(!r->Get("beam:beam_type",tmpi)){
00129     //    if(!fStandardConfig->Get("beam:beam_type",bt)) return w;
00130     bt=fBeam;
00131   }
00132   else{
00133     bt=static_cast<BeamType::BeamType_t>(tmpi);
00134   }
00135   
00136   if(!r->Get("beam:detector",tmpi)){
00137     //    if(!fStandardConfig->Get("beam:detector",det)) return w;
00138     det=fDetector;
00139   }
00140   else{
00141     det=static_cast<Detector::Detector_t>(tmpi);
00142   }
00143   if(!r->Get("beam:high_energy_limit",high)){
00144     //    if(!fStandardConfig->Get("beam:high_energy_limit",high)) return w;
00145     high=fHighRange;
00146   }
00147   if(!r->Get("beam:low_energy_limit",low)){
00148     //    if(!fStandardConfig->Get("beam:low_energy_limit",low)) return w;
00149     low=fLowRange;
00150   }
00151 
00152   w=GetWeight(bt,det,nu_idhep,nu_energy,high,low);
00153   
00154   return w;
00155 }

void BeamEnergyCalculator::ReweightConfigReset (  )  [virtual]

Member Data Documentation

Definition at line 68 of file BeamEnergyCalculator.h.

Referenced by Config(), and GetWeight().

Definition at line 87 of file BeamEnergyCalculator.h.

Referenced by GetHist().

Definition at line 67 of file BeamEnergyCalculator.h.

Referenced by Config(), and GetWeight().

TDirectory* BeamEnergyCalculator::fFluxDir [private]

Definition at line 63 of file BeamEnergyCalculator.h.

Referenced by GetHist().

Definition at line 65 of file BeamEnergyCalculator.h.

Referenced by Config(), and GetHist().

std::string BeamEnergyCalculator::fFluxFileName [private]

Definition at line 64 of file BeamEnergyCalculator.h.

Referenced by Config(), and GetHist().

Definition at line 70 of file BeamEnergyCalculator.h.

Referenced by Config(), and GetWeight().

Definition at line 69 of file BeamEnergyCalculator.h.

Referenced by Config(), and GetWeight().

const double BeamEnergyCalculator::low_energy_cut = 1e-6 [static, private]

Definition at line 55 of file BeamEnergyCalculator.h.

Referenced by DeriveInverseEWeight(), and DeriveWeight().


The documentation for this class was generated from the following files:

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1