src/smear/Acceptance.cxx

0001 /**
0002  \file
0003  Implementation of class Smear::Acceptance.
0004
0005  \author    Michael Savastio
0006  \date      2011-08-19
0007  \copyright 2011 Brookhaven National Lab
0008  */
0009
0010 #include "eicsmear/smear/Acceptance.h"
0011
0012 #include <TDatabasePDG.h>
0013 #include <TLorentzVector.h>
0014 #include <TString.h>
0015
0016 namespace Smear {
0017
0018 Acceptance::~Acceptance() {
0019 }
0020
0021 Acceptance::Acceptance(int genre)
0022 : mGenre(genre)
0023 , mCharge(kAllCharges) {
0024 }
0025
0026 void Acceptance::AddZone(const Zone& z) {
0027   mZones.push_back(z);
0028 }
0029
0030 void Acceptance::SetGenre(int n) {
0031   if (n > 0 && n < 3) {
0032     mGenre = n;
0033   } else {
0034     mGenre = 0;
0035   }  // if
0036 }
0037
0038 void Acceptance::SetCharge(ECharge charge) {
0039   mCharge = charge;
0040 }
0041
0042 void Acceptance::AddParticle(int n) {
0043   mParticles.insert(n);
0044 }
0045
0046 bool Acceptance::Is(const erhic::VirtualParticle& prt) const {
0047   // Check for genre first (em, hadronic, any)
0048   // if (PGenre(prt) == 0 || (mGenre != 0 && PGenre(prt) != mGenre)) {
0049   //   return false;
0050   // }  // if
0051   if ( mGenre != 0 && PGenre(prt) != mGenre ) {
0052     return false;
0053   }  // if
0054
0055   // Check if the particle charge matches the required value.
0056   if (mCharge != kAllCharges) { // Don't need to check if accepting all
0057     // Try to find the particle's charge via its TParticlePDG object.
0058     TParticlePDG* pdg = prt.Id().Info();
0059     if (pdg) {
0060       bool charged = fabs(pdg->Charge()) > 0.;
0061       // Check the charge against the requested value and return false
0062       // if it is incorrect.
0063       if ((kNeutral == mCharge && charged) ||
0064          (kCharged == mCharge && !charged)) {
0065         return false;
0066       }  // if
0067     } else {
0068       // The particle is unknown. We can't guarantee it's charge matches
0069       // the requested value, so return false.
0070       return false;
0071     }  // if
0072   }  // if
0073   // Check against exclusive particle list
0074   if (!mParticles.empty() && mParticles.count(prt.Id()) == 0) {
0075     return false;
0076   }  // if
0077   // If there are no Zones, accept everything that passed genre check
0078   if (mZones.empty()) {
0079     return true;
0080   }  // if
0081   for (unsigned i(0); i < mZones.size(); i++) {
0082     if (mZones.at(i).Contains(prt)) {
0083       return true;
0084     }  // if
0085   }  // for
0086   return false;
0087 }
0088
0089 //
0090 // class Acceptance::CustomCut
0091 //
0092
0093 Acceptance::CustomCut::~CustomCut() {
0094 }
0095
0096 Acceptance::CustomCut::CustomCut()
0097 : mFormula("CustomCutFormula", "0")
0098 , dim(0)
0099 , Kin1(kP)
0100 , Kin2(kTheta)
0101 , Min(-TMath::Infinity())
0102 , Max(TMath::Infinity()) {
0103 }
0104
0105 Acceptance::CustomCutFunction::~CustomCutFunction() {
0106 }
0107
0108 Acceptance::CustomCutFunction::CustomCutFunction( cutfunc func )
0109  : func(func)
0110 { }
0111
0112 // bool Acceptance::CustomCutFunction::Contains(const erhic::VirtualParticle& p) const {
0113 //   return func ( p );
0114 // }
0115
0116 Acceptance::CustomCut::CustomCut(const TString& formula,
0117                                  double min, double max)
0118 : mFormula("CustomCutFormula", "0")
0119 , dim(0)
0120 , Kin1(kP)
0121 , Kin2(kTheta)
0122 , Min(min)
0123 , Max(max) {
0124   TString s(formula);
0125   dim = ParseInputFunction(s, Kin1, Kin2);
0126   if (!IsCoreType(Kin1) || !IsCoreType(Kin2)) {
0127     std::cerr <<
0128     "ERROR! Custom acceptance is not a function of E, p, theta, phi"
0129     << std::endl;
0130   }  // if
0131   if (1 == dim || 2 == dim) {
0132     mFormula = TFormula("CustomCutFormula", s);
0133   } else {
0134     std::cerr <<
0135     "ERROR! Provided custom acceptance is not of dimension 1 or 2."
0136     << std::endl;
0137     return;
0138   }  // if
0139   std::cout << "Added custom cut " << formula << std::endl;
0140 }
0141
0142 bool Acceptance::CustomCut::Contains(
0143                                      const erhic::VirtualParticle& prt) const {
0144   double x = GetVariable(prt, Kin1);
0145   double y(0.);
0146   if (2 == dim) {
0147     y = GetVariable(prt, Kin2);
0148   }  // if
0149   double z = mFormula.Eval(x, y);
0150   return z >= Min && z < Max;
0151 }
0152
0153 //
0154 // class Acceptance::Zone
0155 //
0156
0157 Acceptance::Zone::~Zone() {
0158 }
0159
0160 Acceptance::Zone::Zone(double thMin, double thMax,
0161                        double phMin, double phMax,
0162                        double eMin, double eMax,
0163                        double pMin, double pMax,
0164                        double ptmin, double ptmax,
0165                        double pzmin, double pzmax)
0166 : thetaMin(thMin)
0167 , thetaMax(thMax)
0168 , phiMin(phMin)
0169 , phiMax(phMax)
0170 , EMin(eMin)
0171 , EMax(eMax)
0172 , PMin(pMin)
0173 , PMax(pMax)
0174 , pTMin(ptmin)
0175 , pTMax(ptmax)
0176 , pZMin(pzmin)
0177 , pZMax(pzmax) {
0178 }
0179
0180 void Acceptance::Zone::Add(const CustomCut& cut) {
0181   CustomCuts.push_back(cut);
0182 }
0183
0184 Bool_t Acceptance::Zone::Contains(const erhic::VirtualParticle& prt) const {
0185   bool accept(true);
0186   const double theta = FixTheta(prt.GetTheta());
0187   const double phi = FixPhi(prt.GetPhi());
0188   if (theta < thetaMin || theta > thetaMax) {
0189     accept = false;
0190   } else if (phi < phiMin || phi > phiMax) {
0191     accept = false;
0192   } else if (prt.GetE() < EMin || prt.GetE() > EMax) {
0193     accept = false;
0194   } else if (prt.GetP() < PMin || prt.GetP() > PMax) {
0195     accept = false;
0196   } else if (prt.GetPz() < pZMin || prt.GetPz() > pZMax) {
0197     accept = false;
0198   } else if (prt.GetPt() < pTMin || prt.GetPt() > pTMax) {
0199     accept = false;
0200   }  // if
0201   // If it made it this far, test the custom cut(s)
0202   if (accept) {
0203     for (unsigned j(0); j < CustomCuts.size(); ++j) {
0204       if (!CustomCuts.at(j).Contains(prt)) {
0205         accept = false;
0206         break;
0207       }  // if
0208     }  // for
0209   }  // if
0210   return accept;
0211 }
0212
0213 }  // namespace Smear