Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

BeamEnergyCalculator Class Reference

#include <BeamEnergyCalculator.h>

Inheritance diagram for BeamEnergyCalculator:

WeightCalculator List of all members.

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

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

const double low_energy_cut = 1e-6

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.

Referenced by GetHist().

Constructor & Destructor Documentation

BeamEnergyCalculator::BeamEnergyCalculator (  Registry *  stdconfig = 0  ) 
 

Definition at line 43 of file BeamEnergyCalculator.cxx.

References Config().

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 (   )  [virtual]
 

Reimplemented from WeightCalculator.

Definition at line 58 of file BeamEnergyCalculator.cxx.

References fBeam, fDetector, fFluxFileChanged, fFluxFileName, fHighRange, fLowRange, Registry::Get(), Registry::GetCharString(), 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 Detector::AsString(), det, and s().

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  ,
double  ,
double 
)  [private]
 

Definition at line 231 of file BeamEnergyCalculator.cxx.

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 *  ,
double  ,
double  ,
double 
)  [private]
 

Definition at line 189 of file BeamEnergyCalculator.cxx.

References 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(), CacheMapIter, ConstructName(), det, fCache, fFluxDir, fFluxFileChanged, fFluxFileName, gSystem(), MAXMSG, MSG, and s().

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(), det, GetHist(), 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.

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, Registry::Get(), 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]
 

Reimplemented from WeightCalculator.

Definition at line 101 of file BeamEnergyCalculator.cxx.

References WeightCalculator::SetReweightConfig().

00101                                               {
00102   if(fStandardConfig) SetReweightConfig(fStandardConfig);
00103 }

Member Data Documentation

BeamType::BeamType_t BeamEnergyCalculator::fBeam [private]
 

Definition at line 68 of file BeamEnergyCalculator.h.

Referenced by Config().

CacheMap BeamEnergyCalculator::fCache [private]
 

Definition at line 87 of file BeamEnergyCalculator.h.

Referenced by GetHist().

Detector::Detector_t BeamEnergyCalculator::fDetector [private]
 

Definition at line 67 of file BeamEnergyCalculator.h.

Referenced by Config().

TDirectory* BeamEnergyCalculator::fFluxDir [private]
 

Definition at line 63 of file BeamEnergyCalculator.h.

Referenced by GetHist().

bool BeamEnergyCalculator::fFluxFileChanged [private]
 

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().

double BeamEnergyCalculator::fHighRange [private]
 

Definition at line 70 of file BeamEnergyCalculator.h.

Referenced by Config().

double BeamEnergyCalculator::fLowRange [private]
 

Definition at line 69 of file BeamEnergyCalculator.h.

Referenced by Config().

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

Definition at line 41 of file BeamEnergyCalculator.cxx.

The documentation for this class was generated from the following files:
Generated on Mon Feb 15 11:08:46 2010 for loon by  doxygen 1.3.9.1