src/smear/ParticleID.cxx

0001 /**
0002  \file
0003  Implementation of class Smear::ParticleID.
0004
0005  \author    Michael Savastio
0006  \date      2011-08-12
0007  \copyright 2011 BNL. All rights reserved.
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       }  // for
0045       std::cout << std::endl;
0046     }  // for
0047   }  // for
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     }  // for
0057   }  // for
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       }  // for
0071       t = 0.;
0072     }  // for
0073   }  // for
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   // Get the cumulative probability values for this momentum bin
0081   // and true ID
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       }  // if
0088     } else if (r > values.at(i-1) && r < values.at(i)) {
0089       falseid = FalseIdent.at(i);
0090     }  // if
0091   }  // for
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     }  // if
0100   }  // for
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     }  // if
0109   }  // for
0110      // This is to accomodate the old HERMES format
0111   if (ID < 5 && ID > 0) {
0112     return ID - 1;
0113   }  // for
0114   return -1;
0115 }
0116
0117 void ParticleID::Clear(Option_t* /* option */) {
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   // Open the file and check for errors
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   }  // if
0136      // Returns true if s begins with pattern.
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     // Strip leading whitespace
0149     line.erase(0, line.find_first_not_of(" \t"));
0150     // Feed the line to the stringstream, clearing existing contents first
0151     ss.str("");  // Remove contents
0152     ss.clear();  // Clear flags
0153     ss << line;
0154     int tmpint(0);
0155     // Read true particle IDs
0156     if (starts(line, "!T")) {
0157       ss >> dummy;
0158       while ((ss >> tmpint).good()) {
0159         TrueIdent.push_back(tmpint);
0160       }  // while
0161       gotTrue = !TrueIdent.empty();
0162     } else if (starts(line, "!F")) {
0163       // Read misidentified particle IDs
0164       ss >> dummy;
0165       while ((ss >> tmpint).good()) {
0166         FalseIdent.push_back(tmpint);
0167       }  // while
0168       gotFalse = !FalseIdent.empty();
0169     } else if (starts(line, "!P")) {
0170       // Read number of momentum bins
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       // Read the probabilities for this momentum bin/true ID/false ID
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       }  // if
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       }  // if
0191       pbin -= 1;  // Go from [1, N] to [0, N) for vector index range
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       }  // if
0200       if (PMatrix.empty()) {
0201         SetPMatrixSize();
0202       }  // if
0203       // Set identification probabilities
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       }  // if
0209     }  // if
0210   }  // while
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   }  // if
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         // Generated ID is always positive.
0227         // Keep same sign as input PID i.e. no error in charge sign
0228         if (pid > 0) {
0229           prtOut.SetId (Wild(i, pid) );
0230         } else {
0231           prtOut.SetId( -Wild(i, pid) );
0232         }  // if
0233       }  // if
0234     }  // for
0235   }  // if
0236 }
0237
0238 void ParticleID::Print(Option_t* /* option */) const {
0239   std::cout << "ParticleID using " << PMatPath << std::endl;
0240 }
0241
0242 }  // namespace Smear