src/smear/Device.cxx

0001 /**
0002  \file
0003  Implementation of class Smear::Device.
0004
0005  \author    Michael Savastio
0006  \date      2011-08-19
0007  \copyright 2011 Brookhaven National Lab
0008  */
0009
0010 #include "eicsmear/smear/Device.h"
0011
0012 #include <algorithm>
0013 #include <functional>
0014 #include <iterator>
0015 #include <stdexcept>
0016 #include <string>
0017 #include <vector>
0018
0019 #include <TUUID.h>
0020 #include <TDatabasePDG.h>
0021
0022 #include "eicsmear/smear/FormulaString.h"
0023 #include "eicsmear/smear/ParticleMCS.h"
0024
0025 using std::cout;
0026 using std::cerr;
0027 using std::endl;
0028
0029 namespace Smear {
0030
0031
0032   /// TODO: KK It seems the use of kinematicFunction is unnecessary
0033   /// All it does is
0034   // - create a string that identifies the smeared observable (say, phi)
0035   // - create a string ("phi"), translate it to a TF1 ("x")
0036   // - relate the smeared value to the observable using GetX()
0037   //But this TF1 is always identity, modulo range restrictions that cause obscure errors
0038 bool Device::Init(const TString& kinematicFunction,
0039                   const TString& resolutionFunction, int genre) {
0040   Accept.SetGenre(genre);
0041   // Use FormulaString to parse the kinematic function,
0042   // though we can't use a FormulaString for the kinematic
0043   // function as we need TF1 functionality.
0044   FormulaString f(kinematicFunction.Data());
0045   // The expression has to have exactly one variable.
0046   mSmeared = f.Variables().front();
0047   // Use UUID for ROOT function name to avoid instances clashing
0048   mKinematicFunction = new TF1(TUUID().AsString(),
0049                                f.GetString().c_str(), -1e15, 1.e16);
0050   // Set the resolution function.
0051   mFormula = new FormulaString(resolutionFunction.Data());
0052   // cout << mFormula->GetString() << endl;
0053   return mFormula && mKinematicFunction;
0054 }
0055
0056 Device::Device(KinType type, const TString& formula, EGenre genre)
0057 : mSmeared(type)
0058 , mKinematicFunction(NULL)
0059 , mFormula(NULL) {
0060   Accept.SetGenre(genre);
0061   Init(FormulaString::GetKinName(type), formula, genre);
0062 }
0063
0064 Device::Device(const TString& variable, const TString& resolution,
0065                EGenre genre)
0066 : mSmeared(kInvalidKinType)
0067 , mKinematicFunction(NULL)
0068 , mFormula(NULL) {
0069   Init(variable, resolution, genre);
0070 }
0071
0072 Device::Device(const Device& that)
0073 : Smearer(that)
0074 , mSmeared(that.mSmeared)
0075 , mKinematicFunction(NULL)
0076 , mFormula(NULL)
0077 , mDimensions(that.mDimensions) {
0078   if (that.mKinematicFunction) {
0079     mKinematicFunction = static_cast<TF1*>(
0080         that.mKinematicFunction->Clone(TUUID().AsString()));
0081   }  // if
0082   if (that.mFormula) {
0083     mFormula = static_cast<FormulaString*>(that.mFormula->Clone());
0084   }  // if
0085 }
0086
0087 Device::~Device() {
0088   if (mFormula) {
0089     delete mFormula;
0090     mFormula = NULL;
0091   }  // if
0092   if (mKinematicFunction) {
0093     delete mKinematicFunction;
0094     mKinematicFunction = NULL;
0095   }  // if
0096 }
0097
0098 void Device::Smear(const erhic::VirtualParticle &prt, ParticleMCS &out) {
0099   // Test for acceptance and do nothing if it fails.
0100   if (!Accept.Is(prt)) {
0101     return;
0102   }  // if
0103   // Get each argument for the resolution function from the particle.
0104   const std::vector<KinType> vars = mFormula->Variables();
0105   std::vector<double> args;
0106   for (unsigned i(0); i < vars.size(); ++i) {
0107     args.push_back(GetVariable(prt, vars.at(i)));
0108   }  // for
0109   // Evaluate the quantity to smear and the resolution, then throw
0110   // a random smeared value.
0111   double unsmeared = mKinematicFunction->Eval(GetVariable(prt, mSmeared));
0112   double resolution = mFormula->Eval(args);
0113   double smeared = mDistribution.Generate(unsmeared, resolution);
0114   // mDistribution.Print();
0115   if ( false && abs(prt.Id())==11){
0116     // if ( mSmeared != kE ){std::cerr << " ======================  " << mSmeared << std::endl; throw(-1);}
0117     std::cout << " ======================  mSmeared = " << mSmeared << std::endl;
0118     // prt.Print();
0119     std::cout << " orig eta = " << prt.GetEta()<< std::endl;
0120     std::cout << "unsmeared " << unsmeared << std::endl;
0121     std::cout << "resolution " << resolution << std::endl;
0122     std::cout << "smeared " << smeared << std::endl;
0123     std::cout << "mSmeared " << mSmeared << std::endl;
0124   }
0125   // mKinematicFunction->Print();
0126   if ( mKinematicFunction->GetFormula()->GetExpFormula() != "x" ){
0127     std::cerr << "Formula = " << mKinematicFunction->GetFormula()->GetExpFormula() << std::endl;
0128     throw -1;
0129   }
0130   // if ( smeared < 0 || fabs(mKinematicFunction->GetX(smeared)-smeared)>1e-5 ) cout << mKinematicFunction->GetX(smeared) << " " << smeared << endl;
0131   out.SetVariable(mKinematicFunction->GetX(smeared), mSmeared);
0132   // Fix angles to the correct ranges.
0133   if (kTheta == mSmeared) {
0134     out.SetTheta ( FixTheta(out.GetTheta() ), false);
0135   } else if (kPhi == mSmeared) {
0136     out.SetPhi ( FixPhi(out.GetPhi() ), false);
0137   }  // else if
0138   // Ensure E, p are positive definite
0139   out.HandleBogusValues(mSmeared);
0140   if ( std::isnan( GetVariable (out, mSmeared ) ) ) cout << "Problem. smeared is " << smeared << " but propagated nan in " << mSmeared<< endl;
0141   // std::cout << "Bye Smear" << std::endl << std::endl;
0142 }
0143
0144 Device* Device::Clone(const char* /** Unused */) const {
0145   return new Device(*this);
0146 }
0147
0148 void Device::Print(Option_t* /* option */) const {
0149   const std::string name = FormulaString::GetKinName(mSmeared);
0150   std::cout << "Device smearing " << name << " with sigma(" << name <<
0151   ") = " << mFormula->GetInputString() << std::endl;
0152 }
0153
0154 }  // namespace Smear