eicsmear/smear/ParticleID.h

0001 /**
0002  \file
0003  Declaration of class Smear::ParticleID.
0004
0005  \author    Michael Savastio
0006  \date      2011-08-12
0007  \copyright 2011 BNL. All rights reserved.
0008  */
0009
0010 #ifndef INCLUDE_EICSMEAR_SMEAR_PARTICLEID_H_
0011 #define INCLUDE_EICSMEAR_SMEAR_PARTICLEID_H_
0012
0013 #include <cmath>
0014 #include <fstream>
0015 #include <iostream>
0016 #include <sstream>
0017 #include <vector>
0018
0019 #include <TF1.h>
0020 #include <TF2.h>
0021 #include <TLorentzVector.h>
0022 #include <TRandom3.h>
0023 #include <TROOT.h>
0024 #include <TSystem.h>
0025 #include <TString.h>
0026
0027 #include "eicsmear/smear/Device.h"
0028 #include "eicsmear/erhic/Kinematics.h"
0029 #include "eicsmear/erhic/Particle.h"
0030 #include "eicsmear/erhic/ParticleIdentifier.h"
0031 #include "eicsmear/smear/Smear.h"
0032 #include "eicsmear/smear/Smearer.h"
0033 #include "eicsmear/erhic/VirtualParticle.h"
0034
0035 namespace Smear {
0036
0037 /**
0038  This structure is used to generate particle ID.
0039
0040  The input file containing the P matrix must begin with format
0041  lines beginning with "!T", "!F" and "!P".  For example
0042
0043  !T 211 321 2212
0044  !F 211 321 2212 0
0045  !P 15
0046
0047  The first line tells the PID generator which particles false ID
0048  will be generated for. The second line tells it which particles the
0049  particles in the first line can be identified as. The third line tells
0050  it how many momentum bins there are.
0051
0052  Lines of data should appear as:
0053
0054  1 pbinNumber  pmin pmax  FalseID  P1  P2  P3
0055
0056  The line must begin with a 1 to be read.
0057  pbinNumber is the number of the momentum bin.
0058  pmin and pmax are the bounds of the momentum bin.
0059  FalseID is the ID that the particle will be mis-identified as,
0060  this should be a PDG particle code (unless it is an old Hermes file).
0061  P1 is the probability that the first particle appearing in the "!T"
0062  line will be misidentified as FalseID. Likewise for P2 and P3.
0063
0064  \todo Implement data hiding
0065  \remark Why is the bitwise AND assignment operator&= overloaded?!
0066  */
0067 struct ParticleID : public Smearer {
0068   /**
0069    Default constructor.
0070    Default path for probability matrix is "PIDMatrix.dat".
0071    Also, by default, the particle ID uses smeared kinematic
0072    variables as input, rather than Monte Carlo values.
0073    */
0074   ParticleID();
0075
0076   /**
0077    Constructor.
0078    Initialise the particle misidentification matrix from the named file.
0079    */
0080   explicit ParticleID(TString filename);
0081
0082   /**
0083    Destructor.
0084    */
0085   virtual ~ParticleID();
0086
0087   /**
0088    Set the path to the file containing the particle ID probability matrix,
0089    and read it in.
0090    */
0091   void SetPMatrixPath(TString);
0092
0093   /**
0094    If true, the ParticleID will use the original Monte Carlo values
0095    to generate ID's, rather than smeared values.
0096    Default is false.
0097    */
0098   void SetPIDUseMC(bool useMc);
0099
0100   /**
0101    Return the TRandom3 instance used by the ParticleID.
0102    */
0103   TRandom3& GetRandomGenerator();
0104
0105   /**
0106    Set the seed of the random number generator instance used by the
0107    ParticleID.
0108    See ROOT TRandom3 documentation for more information.
0109    */
0110   void SetRanSeed(int seed);
0111
0112   /**
0113    Set the ParticleID instance to have the same acceptance as the device.
0114    This includes all acceptance zones as well as the list of specific
0115    particles to be detected (if there is one).
0116    */
0117   void GetAcceptanceFromDevice(const Device&);
0118
0119   /**
0120    Resize the PMatrix, (# of pbins, # of true, # of false).
0121    */
0122   void SetPMatrixSize();
0123
0124   /**
0125    Setup the range array used by wild.
0126
0127    Range is used to set up zones in [0,1],
0128    if a random number falls in one of the zones, the associated ID is used.
0129    */
0130   void SetupProbabilityArray();
0131
0132   /**
0133    Returns a copy of this ParticleID.
0134    Inherited from TObject.
0135    The const char* argument has no effect.
0136    */
0137   virtual ParticleID* Clone(const char* = "") const;
0138
0139   /**
0140    Randomly generates a false ID based on the current, a momentum bin
0141    and a true particle ID.
0142    */
0143   int Wild(int pbin, int trueID);
0144
0145   int InListOfTrue(int ID);
0146
0147   int InListOfFalse(int ID);
0148
0149   /**
0150    Read in a P matrix and set up the ParticleID to be ready to generate.
0151    See the documentation for the detector function ReadPIDMatrix.
0152    */
0153   void ReadP(TString filename);
0154
0155   /**
0156    Generates particle ID if the particle is in the list of particles
0157    to be identified, falls within the momentum range of the PMatrix
0158    and falls within acceptence.
0159    New id is stored in prtOut.id.
0160    By default, this will use the momentum prtOut.p to make its
0161    determination, but you can set it to use the values stored in prt
0162    instead using SetPIDUseMC(true).
0163    */
0164   void Smear(const erhic::VirtualParticle&, ParticleMCS&);
0165
0166   /**
0167    Dump the contents of the table to the screen.
0168    */
0169   void Speak();
0170
0171   /**
0172    Clears existing table contents.
0173    */
0174   virtual void Clear(Option_t* = "");
0175
0176   /**
0177    Print information about this device to standard output.
0178    */
0179   virtual void Print(Option_t* = "") const;
0180
0181   TRandom3 Ran;
0182   TString PMatPath;
0183   std::vector<int> TrueIdent;
0184   std::vector<int> FalseIdent;
0185   std::vector<double> PMin;
0186   std::vector<double> PMax;
0187   // Indices are for momentum bin, true ID, and false ID
0188   // PMatrix[i][j][k] stores the probability for a particle of type
0189   // TrueIdent[j] in the momentum bin (PMin[i], PMax[i]) to be
0190   // identified as a particle of type FalseIdent[k].
0191   std::vector< std::vector<std::vector<double> > > PMatrix;
0192   // Range is the cumulative probability distribution of PMatrix.
0193   std::vector< std::vector<std::vector<double> > > Range;
0194   bool bUseMC;
0195
0196   ClassDef(Smear::ParticleID, 1)
0197 };
0198
0199 inline void ParticleID::SetPMatrixPath(TString str) {
0200   PMatPath = str;
0201   ReadP(PMatPath);
0202 }
0203
0204 inline void ParticleID::SetPIDUseMC(bool d) {
0205   bUseMC = d;
0206 }
0207
0208 inline TRandom3& ParticleID::GetRandomGenerator() {
0209   return Ran;
0210 }
0211
0212 inline void ParticleID::SetRanSeed(int n) {
0213   Ran.SetSeed(n);
0214 }
0215
0216 inline void ParticleID::GetAcceptanceFromDevice(const Device& dev) {
0217   Accept = dev.Accept;
0218 }
0219
0220 inline ParticleID* ParticleID::Clone(const char*) const {
0221   // const char* argument comes from TObject::Clone(), usused here.
0222   return new ParticleID(*this);
0223 }
0224
0225 }  // namespace Smear
0226
0227 #endif  // INCLUDE_EICSMEAR_SMEAR_PARTICLEID_H_