00001 /*************************************************************************** 00002 * 00003 * $Id: StTpcDedxPidAlgorithm.cxx,v 2.29 2012/05/07 14:42:58 fisyak Exp $ 00004 * 00005 * Author: Thomas Ullrich, Sep 1999 00006 *************************************************************************** 00007 * 00008 * Description: 00009 * 00010 *************************************************************************** 00011 * 00012 * $Log: StTpcDedxPidAlgorithm.cxx,v $ 00013 * Revision 2.29 2012/05/07 14:42:58 fisyak 00014 * Add handilings for Track to Fast Detectors Matching 00015 * 00016 * Revision 2.28 2010/11/15 22:11:45 fisyak 00017 * Restore proper nSigma 00018 * 00019 * Revision 2.27 2010/08/31 20:15:11 fisyak 00020 * Clean up 00021 * 00022 * Revision 2.26 2008/03/13 16:56:39 ullrich 00023 * Add protection against missing or wrong mTraits->mean() 00024 * 00025 * Revision 2.25 2005/05/13 20:31:14 fisyak 00026 * Calculated nSigma wrt modified Bichsel I70M to account saturation 00027 * 00028 * Revision 2.24 2004/03/30 02:45:07 fisyak 00029 * return DBL_MAX instead of 0 if No mTraits 00030 * 00031 * Revision 2.23 2003/10/25 00:12:48 fisyak 00032 * Replace BetheBloch::Sirrf by m_Bichsel->GetI70 for nSigma calculations 00033 * 00034 * Revision 2.22 2003/10/21 14:23:23 fisyak 00035 * Switch from primary to global track momentum for Nsigma calculations 00036 * 00037 * Revision 2.21 2003/09/02 17:58:06 perev 00038 * gcc 3.2 updates + WarnOff 00039 * 00040 * Revision 2.20 2003/08/02 01:14:11 perev 00041 * warnOff 00042 * 00043 * Revision 2.19 2003/04/30 18:05:55 fisyak 00044 * Add P03ia flag, which fixes P03ia MuDst 00045 * 00046 * Revision 2.18 2002/04/04 01:42:34 jeromel 00047 * Extra fix for multiple instantiation (and multiple delete). 00048 * 00049 * Revision 2.17 2002/04/01 20:13:56 jeromel 00050 * Heuu !! if (theBetheBloch) delete ... 00051 * 00052 * Revision 2.16 2002/03/26 23:10:09 ullrich 00053 * Changed instantiaton of BetheBloch. 00054 * 00055 * Revision 2.15 2002/03/25 21:02:51 ullrich 00056 * Made BetheBloch a static global variable. 00057 * 00058 * Revision 2.14 2002/02/06 23:00:54 ullrich 00059 * Added float.h. 00060 * 00061 * Revision 2.13 2001/04/27 21:41:07 ullrich 00062 * Fixed bug. 00063 * 00064 * Revision 2.12 2001/04/24 15:38:08 fisyak 00065 * Add switch (use length) for calibrated and non calibrated data 00066 * 00067 * Revision 2.11 2001/04/05 04:00:56 ullrich 00068 * Replaced all (U)Long_t by (U)Int_t and all redundant ROOT typedefs. 00069 * 00070 * Revision 2.10 2000/10/26 18:33:20 calderon 00071 * Return DBL_MAX in the numberOfSigma() function when 00072 * the number of dE/dx points is equal to zero. 00073 * The function sigmaPidFunction() already checked this, but 00074 * when this happens, it probably means that mTraits->mean() 00075 * is undefined, so it's better to exit immediately from the 00076 * function. 00077 * 00078 * Revision 2.9 2000/07/28 16:55:35 calderon 00079 * Mike's version of StTpcDedxPidAlgorithm using BetheBloch 00080 * class lookup table from 00081 * StarClassLibrary and parameterization of resolution obtained 00082 * from the first two weeks of July 2000 Data. 00083 * 00084 * Revision 2.8 2000/05/19 18:33:45 ullrich 00085 * Minor changes (add const) to cope with modified StArray. 00086 * 00087 * Revision 2.7 2000/04/20 16:47:31 ullrich 00088 * Check for null pointer added. 00089 * 00090 * Revision 2.6 2000/03/02 12:43:49 ullrich 00091 * Method can be passed as argument to constructor. Default 00092 * method is truncated mean. 00093 * 00094 * Revision 2.5 1999/12/21 15:09:11 ullrich 00095 * Modified to cope with new compiler version on Sun (CC5.0). 00096 * 00097 * Revision 2.4 1999/12/03 13:18:18 ullrich 00098 * Fixed problem on Sun CC4.2 (dynamic_cast) and switched 00099 * off cut on number of points. 00100 * 00101 * Revision 2.3 1999/12/02 16:35:34 ullrich 00102 * Added method to return the stored dE/dx traits 00103 * 00104 * Revision 2.2 1999/10/28 22:27:01 ullrich 00105 * Adapted new StArray version. First version to compile on Linux and Sun. 00106 * 00107 * Revision 2.1 1999/10/13 19:45:24 ullrich 00108 * Initial Revision 00109 * 00110 **************************************************************************/ 00111 #include <typeinfo> 00112 #include <cmath> 00113 #include <float.h> 00114 #include "StTpcDedxPidAlgorithm.h" 00115 #include "StTrack.h" 00116 #include "StGlobalTrack.h" 00117 #include "StParticleTypes.hh" 00118 #include "StEnumerations.h" 00119 #include "StDedxPidTraits.h" 00120 #include "StTrackNode.h" 00121 #include "StTrackGeometry.h" 00122 #include "BetheBloch.h" 00123 #include "StBichsel/Bichsel.h" 00124 00125 static Bichsel *m_Bichsel = 0; 00126 static const char rcsid[] = "$Id: StTpcDedxPidAlgorithm.cxx,v 2.29 2012/05/07 14:42:58 fisyak Exp $"; 00127 00128 StTpcDedxPidAlgorithm::StTpcDedxPidAlgorithm(StDedxMethod dedxMethod) 00129 : mTraits(0), mTrack(0), mDedxMethod(dedxMethod) 00130 { 00131 // 00132 // Add all particles we want to get 00133 // checked in operator(). 00134 // 00135 mParticles.push_back(StPionMinus::instance()); 00136 mParticles.push_back(StPionPlus::instance()); 00137 mParticles.push_back(StKaonMinus::instance()); 00138 mParticles.push_back(StKaonPlus::instance()); 00139 mParticles.push_back(StProton::instance()); 00140 } 00141 00142 StParticleDefinition* 00143 StTpcDedxPidAlgorithm::operator() (const StTrack& track, const StSPtrVecTrackPidTraits& vec) 00144 { 00145 // 00146 // Select the info we need. 00147 // Here we ignore different kinds of dE/dx calculations in 00148 // the TPC and select the method 00149 // 00150 mTraits = 0; 00151 mTrack = &track; 00152 for (UInt_t i=0; i<vec.size(); i++) { 00153 const StDedxPidTraits *p = dynamic_cast<const StDedxPidTraits*>(vec[i]); 00154 if (p && p->detector() == kTpcId && p->method() == mDedxMethod) mTraits = p; 00155 } 00156 if (!mTraits) return 0; // no info available 00157 00158 // 00159 // Scan the list of particles we want to check and 00160 // return the most probable. 00161 // 00162 Double_t sigma, minSigma = 100000; 00163 UInt_t minIndex = mParticles.size(); 00164 for (UInt_t k=0; k<mParticles.size(); k++) { 00165 if (mParticles[k]->charge()*mTrack->geometry()->charge() > 0) { // require same charge sign 00166 if ((sigma = fabs(numberOfSigma(mParticles[k]))) < minSigma) { 00167 minIndex = k; 00168 minSigma = fabs(sigma); 00169 } 00170 } 00171 } 00172 return minIndex < mParticles.size() ? mParticles[minIndex] : 0; 00173 } 00174 00175 Double_t 00176 StTpcDedxPidAlgorithm::numberOfSigma(const StParticleDefinition* particle) const 00177 { 00178 if (!mTraits) return DBL_MAX; 00179 00180 if (mTraits->numberOfPoints()==0) return DBL_MAX; 00181 // sigmaPidFunction already checks this, but when number of dE/dx points is = 0, 00182 // mTraits->mean() is probably undefined too (have to check this) 00183 // so might as well exit here. 00184 00185 // returns the number of sigma a tracks dedx is away from 00186 // the expected mean for a track for a particle of this mass 00187 Double_t dedx_expected; 00188 Double_t dedx_resolution; 00189 Double_t momentum; 00190 Double_t z = -999; 00191 if (mTraits->mean() > 0) { 00192 const StGlobalTrack *gTrack = 00193 static_cast<const StGlobalTrack*>( mTrack->node()->track(global)); 00194 if (gTrack && mTraits->length() > 0 ) { 00195 momentum = abs(gTrack->geometry()->momentum()); 00196 if (! m_Bichsel) m_Bichsel = Bichsel::Instance(); 00197 Double_t log2dX = mTraits->log2dX(); 00198 if (log2dX <= 0) log2dX = 1; 00199 dedx_expected = 1.e-6*m_Bichsel->GetI70M(log10(momentum/particle->mass()),log2dX); 00200 dedx_resolution = mTraits->errorOnMean(); 00201 if (dedx_resolution > 0) 00202 z = ::log(mTraits->mean()/dedx_expected)/dedx_resolution; 00203 } 00204 } 00205 return z; 00206 }
1.5.9