File indexing completed on 2024-09-28 07:02:51
0001
0002
0003
0004
0005
0006
0007
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 }
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
0048
0049
0050
0051 if ( mGenre != 0 && PGenre(prt) != mGenre ) {
0052 return false;
0053 }
0054
0055
0056 if (mCharge != kAllCharges) {
0057
0058 TParticlePDG* pdg = prt.Id().Info();
0059 if (pdg) {
0060 bool charged = fabs(pdg->Charge()) > 0.;
0061
0062
0063 if ((kNeutral == mCharge && charged) ||
0064 (kCharged == mCharge && !charged)) {
0065 return false;
0066 }
0067 } else {
0068
0069
0070 return false;
0071 }
0072 }
0073
0074 if (!mParticles.empty() && mParticles.count(prt.Id()) == 0) {
0075 return false;
0076 }
0077
0078 if (mZones.empty()) {
0079 return true;
0080 }
0081 for (unsigned i(0); i < mZones.size(); i++) {
0082 if (mZones.at(i).Contains(prt)) {
0083 return true;
0084 }
0085 }
0086 return false;
0087 }
0088
0089
0090
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
0113
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 }
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 }
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 }
0149 double z = mFormula.Eval(x, y);
0150 return z >= Min && z < Max;
0151 }
0152
0153
0154
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 }
0201
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 }
0208 }
0209 }
0210 return accept;
0211 }
0212
0213 }