File indexing completed on 2024-09-28 07:02:52
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010 #include "eicsmear/smear/ParticleID.h"
0011
0012 #include <sstream>
0013 #include <string>
0014 #include <vector>
0015
0016 namespace Smear {
0017
0018 ParticleID::ParticleID()
0019 : Ran(0)
0020 , PMatPath("PIDMatrix.dat")
0021 , bUseMC(false) {
0022 ReadP(PMatPath);
0023 }
0024
0025 ParticleID::ParticleID(TString filename)
0026 : Ran(0)
0027 , PMatPath(filename)
0028 , bUseMC(false) {
0029 ReadP(PMatPath);
0030 }
0031
0032 ParticleID::~ParticleID() {
0033 }
0034
0035 void ParticleID::Speak() {
0036 std::cout.setf(std::ios::fixed);
0037 std::cout.precision(6);
0038 for (unsigned i(0); i < PMax.size(); i++) {
0039 for (unsigned k(0); k < FalseIdent.size(); k++) {
0040 std::cout << 1 << '\t' << i + 1 << '\t' << PMin.at(i) << '\t' <<
0041 PMax.at(i) << '\t' << k + 1;
0042 for (unsigned j(0); j < TrueIdent.size(); j++) {
0043 std::cout << '\t' << PMatrix.at(i).at(j).at(k);
0044 }
0045 std::cout << std::endl;
0046 }
0047 }
0048 }
0049
0050 void ParticleID::SetPMatrixSize() {
0051 PMatrix.resize(PMin.size());
0052 for (unsigned i(0); i < PMatrix.size(); i++) {
0053 PMatrix.at(i).resize(TrueIdent.size());
0054 for (unsigned j(0); j < PMatrix.at(i).size(); j++) {
0055 PMatrix.at(i).at(j).resize(FalseIdent.size());
0056 }
0057 }
0058 }
0059
0060 void ParticleID::SetupProbabilityArray() {
0061 double t = 0.;
0062 Range.resize(PMatrix.size());
0063 for (unsigned i(0); i < Range.size(); i++) {
0064 Range.at(i).resize(PMatrix.at(i).size());
0065 for (unsigned j(0); j < Range.at(i).size(); j++) {
0066 Range.at(i).at(j).resize(PMatrix.at(i).at(j).size());
0067 for (unsigned k = 0; k < Range.at(i).at(j).size(); k++) {
0068 t += PMatrix.at(i).at(j).at(k);
0069 Range.at(i).at(j).at(k) = t;
0070 }
0071 t = 0.;
0072 }
0073 }
0074 }
0075
0076 int ParticleID::Wild(int pbin, int trueID) {
0077 const double r = Ran.Rndm();
0078 const int k = InListOfTrue(trueID);
0079 int falseid(-999999999);
0080
0081
0082 const std::vector<double>& values = Range.at(pbin).at(k);
0083 for (unsigned i(0); i < values.size(); i++) {
0084 if (0 == i) {
0085 if (r < values.at(i)) {
0086 falseid = FalseIdent.at(i);
0087 }
0088 } else if (r > values.at(i-1) && r < values.at(i)) {
0089 falseid = FalseIdent.at(i);
0090 }
0091 }
0092 return falseid;
0093 }
0094
0095 int ParticleID::InListOfTrue(int ID) {
0096 for (unsigned i(0); i < TrueIdent.size(); i++) {
0097 if (TrueIdent.at(i) == abs(ID)) {
0098 return i;
0099 }
0100 }
0101 return -1;
0102 }
0103
0104 int ParticleID::InListOfFalse(int ID) {
0105 for (unsigned i(0); i < FalseIdent.size(); i++) {
0106 if (FalseIdent.at(i) == abs(ID)) {
0107 return i;
0108 }
0109 }
0110
0111 if (ID < 5 && ID > 0) {
0112 return ID - 1;
0113 }
0114 return -1;
0115 }
0116
0117 void ParticleID::Clear(Option_t* ) {
0118 TrueIdent.clear();
0119 FalseIdent.clear();
0120 PMin.clear();
0121 PMax.clear();
0122 PMatrix.clear();
0123 Range.clear();
0124 }
0125
0126 void ParticleID::ReadP(TString filename) {
0127 Clear("");
0128
0129 std::ifstream Qfile;
0130 Qfile.open(filename);
0131 if (!Qfile.good()) {
0132 std::cerr << "Error in ParticleID: unable to open file "
0133 << filename << std::endl;
0134 return;
0135 }
0136
0137 struct StartsWith {
0138 bool operator()(const std::string& s,
0139 const std::string& pattern) const {
0140 return s.find(pattern, 0) == 0;
0141 }
0142 }
0143 starts;
0144 std::stringstream ss;
0145 std::string line, dummy;
0146 bool gotTrue(false), gotFalse(false), gotBins(false);
0147 while (std::getline(Qfile, line).good()) {
0148
0149 line.erase(0, line.find_first_not_of(" \t"));
0150
0151 ss.str("");
0152 ss.clear();
0153 ss << line;
0154 int tmpint(0);
0155
0156 if (starts(line, "!T")) {
0157 ss >> dummy;
0158 while ((ss >> tmpint).good()) {
0159 TrueIdent.push_back(tmpint);
0160 }
0161 gotTrue = !TrueIdent.empty();
0162 } else if (starts(line, "!F")) {
0163
0164 ss >> dummy;
0165 while ((ss >> tmpint).good()) {
0166 FalseIdent.push_back(tmpint);
0167 }
0168 gotFalse = !FalseIdent.empty();
0169 } else if (starts(line, "!P")) {
0170
0171 int pbin;
0172 ss >> dummy >> pbin;
0173 PMin.resize(pbin);
0174 PMax.resize(pbin);
0175 gotBins = !PMin.empty();
0176 } else if (starts(line, "1") || starts(line, "2") || starts(line, "3")) {
0177
0178 if (!(gotTrue && gotFalse && gotBins)) {
0179 std::cerr << "Error in ParticleID: " <<
0180 "P matrix input file has bad or missing format lines.\n";
0181 return;
0182 }
0183 int table, pbin, pid;
0184 double pmin, pmax, p1, p2, p3;
0185 ss >> table >> pbin >> pmin >> pmax >> pid >> p1 >> p2 >> p3;
0186 if ((unsigned)pbin > PMin.size()) {
0187 std::cerr << "Error in ParticleID: " <<
0188 "Out of bounds momentum bin listing.\n";
0189 return;
0190 }
0191 pbin -= 1;
0192 PMin.at(pbin) = pmin;
0193 PMax.at(pbin) = pmax;
0194 pid = InListOfFalse(pid);
0195 if ((unsigned)pid > FalseIdent.size()) {
0196 std::cerr << "Error in ParticleID: " <<
0197 "P matrix has bad particle listing.\n";
0198 return;
0199 }
0200 if (PMatrix.empty()) {
0201 SetPMatrixSize();
0202 }
0203
0204 if (1 == table) {
0205 PMatrix.at(pbin).at(0).at(pid) = p1;
0206 PMatrix.at(pbin).at(1).at(pid) = p2;
0207 PMatrix.at(pbin).at(2).at(pid) = p3;
0208 }
0209 }
0210 }
0211 SetupProbabilityArray();
0212 }
0213
0214 void ParticleID::Smear(const erhic::VirtualParticle& prt,
0215 ParticleMCS& prtOut) {
0216 double momentum(0.);
0217 if (bUseMC) {
0218 momentum = prt.GetP();
0219 } else {
0220 momentum = prtOut.GetP();
0221 }
0222 const int pid = prt.Id();
0223 if (InListOfTrue(pid) != -1 && Accept.Is(prt)) {
0224 for (unsigned i(0); i < PMin.size(); i++) {
0225 if (momentum > PMin.at(i) && momentum < PMax.at(i)) {
0226
0227
0228 if (pid > 0) {
0229 prtOut.SetId (Wild(i, pid) );
0230 } else {
0231 prtOut.SetId( -Wild(i, pid) );
0232 }
0233 }
0234 }
0235 }
0236 }
0237
0238 void ParticleID::Print(Option_t* ) const {
0239 std::cout << "ParticleID using " << PMatPath << std::endl;
0240 }
0241
0242 }