ParticleID/RICHAnalysis/PIDProbabilities.cc

0001 // Some documentation //
0002
0003 #include "PIDProbabilities.h"
0004 #include "Poisson.h"
0005
0006 #include <TDatabasePDG.h>
0007
0008 using namespace std;
0009
0010
0011 PIDProbabilities::PIDProbabilities()
0012 {
0013   _pdg = new TDatabasePDG();
0014   _poisson = new Poisson();
0015 }
0016
0017 bool
0018 PIDProbabilities::particle_probs( vector<float> angles, double momentum, double index, long double probs[4] ){
0019
0020   bool use_reconstructed_momentum = false;
0021   bool use_truth_momentum = false;
0022   bool use_emission_momentum = true;
0023
0024
0025   /* Initialize particle masses, probably want to use pdg database */
0026   int pid[4] = { 11, 211, 321, 2212 };
0027   double m[4] = {
0028     _pdg->GetParticle(pid[0])->Mass(),
0029     _pdg->GetParticle(pid[1])->Mass(),
0030     _pdg->GetParticle(pid[2])->Mass(),
0031     _pdg->GetParticle(pid[3])->Mass()
0032   };
0033
0034   /* Set angle window by taking smallest difference, highest momentum in question (70 GeV) */
0035   double test_p = 70;
0036   double windowx2 = acos( sqrt( m[1]*m[1] + test_p*test_p ) / index / test_p ) - acos( sqrt( m[2]*m[2] + test_p*test_p ) / index / test_p );
0037   double window = windowx2/2;
0038
0039   /* Determine expectation value for each particle */
0040   double beta[4] = {
0041     momentum/sqrt( m[0]*m[0] + momentum*momentum ),
0042     momentum/sqrt( m[1]*m[1] + momentum*momentum ),
0043     momentum/sqrt( m[2]*m[2] + momentum*momentum ),
0044     momentum/sqrt( m[3]*m[3] + momentum*momentum )
0045   };
0046
0047   double theta_expect[4] = {
0048     acos( 1/( index * beta[0] ) ),
0049     acos( 1/( index * beta[1] ) ),
0050     acos( 1/( index * beta[2] ) ),
0051     acos( 1/( index * beta[3] ) )
0052   };
0053
0054
0055   /* Calculate average counts based on sims */
0056   double counts_cal[4] = {0,0,0,0};
0057   if ( use_reconstructed_momentum )
0058     {
0059       cout << "Reconstructed momentum function undefined now!" << endl;
0060       counts_cal[0] = 0;
0061       counts_cal[1] = 0;
0062       counts_cal[2] = 0;
0063       counts_cal[3] = 0;
0064     }
0065   if ( use_truth_momentum )
0066     {
0067       cout << "Truth momentum function undefined now!" << endl;
0068       counts_cal[0] = 0;
0069       counts_cal[1] = 0;
0070       counts_cal[2] = 0;
0071       counts_cal[3] = 0;
0072     }
0073   if ( use_emission_momentum )
0074     {
0075       // Proton function still needs to be implemented after sims are finished //
0076       counts_cal[0] = 29.5;
0077       counts_cal[1] = 25.8 * acos( sqrt(178.4 + momentum*momentum) / (2.62*momentum) );
0078       counts_cal[2] = 22.08 * acos( sqrt(982600 + momentum*momentum) / (52.93*momentum) );
0079       counts_cal[3] = 0;
0080     }
0081
0082
0083   /* Count hits within window */
0084   int counts[4] = {0,0,0,0};
0085   for (int i=0; i<300; i++){
0086     for (int j=0; j<4; j++){
0087       if ( angles[i] > (theta_expect[j] - window) && angles[i] < (theta_expect[j] + window) )
0088     counts[j]++;
0089     }
0090   }
0091
0092   /* Determine particle probabilities, Poisson probability mass function */
0093   long double probs_norm = 0;
0094   for (int k=1; k<3; k++)
0095     {
0096       probs[k] = _poisson->poisson_prob( counts_cal[k], counts[k] );
0097       probs_norm += probs[k];
0098     }
0099   for (int k=1; k<3; k++)
0100     probs[k] /= probs_norm;
0101
0102
0103   /* Alert that this may have been a glitched event if very few photons */
0104   double allcounts = counts[0] + counts[1] + counts[2] + counts[3];
0105   double allcal = counts_cal[0] + counts_cal[1] + counts_cal[2] + counts_cal[3];
0106   if ( allcounts / allcal < 0.2 )
0107     cout << "WARNING: Very few photons in event. May be an unintended track!" << endl;
0108
0109
0110   return true;
0111 }