Bug Summary

File:libraries/PID/DNeutralParticleHypothesis_factory.cc
Location:line 129, column 27
Description:Dereference of null pointer

Annotated Source Code

1// $Id$
2//
3// File: DNeutralParticleHypothesis_factory.cc
4// Created: Thu Dec 3 17:27:55 EST 2009
5// Creator: pmatt (on Linux ifarml6 2.6.18-128.el5 x86_64)
6//
7
8
9#include <iostream>
10#include <iomanip>
11using namespace std;
12
13#include <TMath.h>
14
15#include "DNeutralParticleHypothesis_factory.h"
16using namespace jana;
17
18//------------------
19// init
20//------------------
21jerror_t DNeutralParticleHypothesis_factory::init(void)
22{
23 return NOERROR;
24}
25
26//------------------
27// brun
28//------------------
29jerror_t DNeutralParticleHypothesis_factory::brun(jana::JEventLoop *locEventLoop, int runnumber)
30{
31 // Get Target parameters from XML
32 DApplication *locApplication = dynamic_cast<DApplication*> (locEventLoop->GetJApplication());
33 DGeometry *locGeometry = locApplication ? locApplication->GetDGeometry(runnumber):NULL__null;
34 dTargetLength = 30.0;
35 double locTargetCenterZ = 65.0;
36 dTargetRadius = 1.5; //FIX: grab from database!!!
37 if(locGeometry)
38 {
39 locGeometry->GetTargetZ(locTargetCenterZ);
40 locGeometry->GetTargetLength(dTargetLength);
41 }
42 dTargetCenter.SetXYZ(0.0, 0.0, locTargetCenterZ);
43
44
45 // Get the particle ID algorithms
46 vector<const DParticleID *> locPIDAlgorithms;
47 locEventLoop->Get(locPIDAlgorithms);
48 if(locPIDAlgorithms.size() < 1){
49 _DBG_std::cerr<<"DNeutralParticleHypothesis_factory.cc"<<
":"<<49<<" "<<"Unable to get a DParticleID object! NO PID will be done!"<<endl;
50 return RESOURCE_UNAVAILABLE;
51 }
52 // Drop the const qualifier from the DParticleID pointer (I'm surely going to hell for this!)
53 dPIDAlgorithm = const_cast<DParticleID*>(locPIDAlgorithms[0]);
54
55 // Warn user if something happened that caused us NOT to get a dPIDAlgorithm object pointer
56 if(!dPIDAlgorithm){
57 _DBG_std::cerr<<"DNeutralParticleHypothesis_factory.cc"<<
":"<<57<<" "<<"Unable to get a DParticleID object! NO PID will be done!"<<endl;
58 return RESOURCE_UNAVAILABLE;
59 }
60
61 return NOERROR;
62}
63
64//------------------
65// evnt
66//------------------
67jerror_t DNeutralParticleHypothesis_factory::evnt(jana::JEventLoop *locEventLoop, int eventnumber)
68{
69 unsigned int loc_i, loc_j, loc_k, locNDF = 1;
70 bool locShowerMatchFlag;
71 float locMass, locMomentum, locShowerEnergy, locParticleEnergy, locPathLength, locFlightTime, locProjectedTime, locTimeDifference;
72 float locParticleEnergyUncertainty, locShowerEnergyUncertainty, locTimeDifferenceVariance, locChiSq, locFOM;
73 DVector3 locPathVector;
74
75 const DNeutralShower *locNeutralShower;
76 const DChargedTrackHypothesis *locChargedTrackHypothesis;
77 DNeutralParticleHypothesis *locNeutralParticleHypothesis;
78 DMatrixDSym locVariances, locErrorMatrix;
79 vector<const DBCALShower*> locAssociatedBCALShowers_NeutralShower;
80 vector<const DFCALShower*> locAssociatedFCALShowers_NeutralShower;
81 vector<const DBCALShower*> locAssociatedBCALShowers_ChargedTrack;
82 vector<const DFCALShower*> locAssociatedFCALShowers_ChargedTrack;
83
84 vector<const DChargedTrack*> locChargedTracks;
85 vector<const DNeutralShower*> locNeutralShowers;
86 locEventLoop->Get(locChargedTracks);
87 locEventLoop->Get(locNeutralShowers);
88
89 vector<Particle_t> locPIDHypotheses;
90 locPIDHypotheses.push_back(Gamma);
91 locPIDHypotheses.push_back(Neutron);
92
93 double locStartTime = 0.0; //fix when rf info available!!
94 DLorentzVector locSpacetimeVertex(dTargetCenter, locStartTime);
95 double locVertexTimeUncertainty = 0.0;
96
97 // Loop over DNeutralShowers
98 for (loc_i = 0; loc_i < locNeutralShowers.size(); loc_i++){
1
Loop condition is true. Entering loop body
→
5
←
Loop condition is true. Entering loop body
→
9
←
Loop condition is true. Entering loop body
→
13
←
Loop condition is true. Entering loop body
→
99 locNeutralShower = locNeutralShowers[loc_i];
100 locNeutralShower->GetT(locAssociatedBCALShowers_NeutralShower);
101 locNeutralShower->GetT(locAssociatedFCALShowers_NeutralShower);
102
103 // If the DNeutralShower is matched to the DChargedTrackHypothesis with the highest FOM of ANY DChargedTrack, skip it
104 locShowerMatchFlag = false;
105 for (loc_j = 0; loc_j < locChargedTracks.size(); loc_j++){
2
←
Loop condition is false. Execution continues on line 122
→
6
←
Loop condition is false. Execution continues on line 122
→
10
←
Loop condition is false. Execution continues on line 122
→
14
←
Loop condition is false. Execution continues on line 122
→
106 locChargedTrackHypothesis = locChargedTracks[loc_j]->dChargedTrackHypotheses[0];
107 locChargedTrackHypothesis->GetT(locAssociatedBCALShowers_ChargedTrack);
108 locChargedTrackHypothesis->GetT(locAssociatedFCALShowers_ChargedTrack);
109 if ((locAssociatedBCALShowers_ChargedTrack.size() > 0) && (locAssociatedBCALShowers_NeutralShower.size() > 0)){
110 if (locAssociatedBCALShowers_ChargedTrack[0]->id == locAssociatedBCALShowers_NeutralShower[0]->id){
111 locShowerMatchFlag = true;
112 break;
113 }
114 }
115 if ((locShowerMatchFlag == false) && (locAssociatedFCALShowers_ChargedTrack.size() > 0) && (locAssociatedFCALShowers_NeutralShower.size() > 0)){
116 if (locAssociatedFCALShowers_ChargedTrack[0]->id == locAssociatedFCALShowers_NeutralShower[0]->id){
117 locShowerMatchFlag = true;
118 break;
119 }
120 }
121 }
122 if (locShowerMatchFlag == true)
3
←
Taking false branch
→
7
←
Taking false branch
→
11
←
Taking false branch
→
15
←
Taking false branch
→
123 continue; //shower matched to a DChargedTrackHypothesis with the highest FOM, not a neutral
124
125 // Loop over vertices and PID hypotheses & create DNeutralParticleHypotheses for each combination
126 for (loc_k = 0; loc_k < locPIDHypotheses.size(); loc_k++){
4
←
Loop condition is false. Execution continues on line 98
→
8
←
Loop condition is false. Execution continues on line 98
→
12
←
Loop condition is false. Execution continues on line 98
→
16
←
Loop condition is true. Entering loop body
→
127
128 // Calculate DNeutralParticleHypothesis Quantities (projected time at vertex for given id, etc.)
129 locMass = ParticleMass(locPIDHypotheses[loc_k]);
17
←
Dereference of null pointer
130 locShowerEnergy = locNeutralShower->dEnergy;
131 locParticleEnergy = locShowerEnergy; //need to correct this for neutrons!
132 if (locParticleEnergy < locMass)
133 continue; //not enough energy for PID hypothesis
134
135 locShowerEnergyUncertainty = locNeutralShower->dEnergyUncertainty;
136 locParticleEnergyUncertainty = locShowerEnergyUncertainty; //need to correct this for neutrons!
137
138 locPathVector = locNeutralShower->dSpacetimeVertex.Vect() - locSpacetimeVertex.Vect();
139 locPathLength = locPathVector.Mag();
140 if(!(locPathLength > 0.0))
141 continue; //invalid, will divide by zero when creating error matrix, so skip!
142 locMomentum = sqrt(locParticleEnergy*locParticleEnergy - locMass*locMass);
143 locFlightTime = locPathLength*locParticleEnergy/(locMomentum*SPEED_OF_LIGHT29.9792);
144 locProjectedTime = locSpacetimeVertex.T() - locFlightTime;
145
146 // Calculate DNeutralParticleHypothesis FOM
147 locTimeDifference = locSpacetimeVertex.T() - locProjectedTime;
148 locTimeDifferenceVariance = 1.0; //completely random, ok because ID disabled for neutrons anyway
149 locChiSq = locTimeDifference*locTimeDifference/locTimeDifferenceVariance;
150 locFOM = TMath::Prob(locChiSq, locNDF);
151 if(locPIDHypotheses[loc_k] == Neutron)
152 locFOM = -1.0; //disables neutron ID until the neutron energy is calculated correctly from the deposited energy in the shower
153
154 // Build DNeutralParticleHypothesis // dEdx not set
155 locNeutralParticleHypothesis = new DNeutralParticleHypothesis;
156 locNeutralParticleHypothesis->AddAssociatedObject(locNeutralShower);
157
158 locNeutralParticleHypothesis->setMass(locMass);
159 locNeutralParticleHypothesis->setCharge(0.0);
160
161 Calc_Variances(locNeutralShower, locParticleEnergyUncertainty, locVariances);
162 Build_ErrorMatrix(locPathVector, locParticleEnergy, locVariances, locErrorMatrix);
163
164 locNeutralParticleHypothesis->setErrorMatrix(locErrorMatrix);
165 locNeutralParticleHypothesis->clearTrackingErrorMatrix();
166 locPathVector.SetMag(locMomentum);
167 locNeutralParticleHypothesis->setMomentum(locPathVector);
168 locNeutralParticleHypothesis->setPosition(locSpacetimeVertex.Vect());
169 locNeutralParticleHypothesis->setT0(locSpacetimeVertex.T(), locVertexTimeUncertainty, SYS_NULL);
170 locNeutralParticleHypothesis->setT1(locNeutralShower->dSpacetimeVertex.T(), locNeutralShower->dSpacetimeVertexUncertainties.T(), locNeutralShower->dDetectorSystem);
171 locNeutralParticleHypothesis->setPathLength(locPathLength, 0.0); //zero uncertainty (for now)
172
173 locNeutralParticleHypothesis->dPID = locPIDHypotheses[loc_k];
174 locNeutralParticleHypothesis->dChiSq = locChiSq;
175 locNeutralParticleHypothesis->dNDF = locNDF;
176 locNeutralParticleHypothesis->dFOM = locFOM;
177
178 _data.push_back(locNeutralParticleHypothesis);
179 } //end PID loop
180 } //end DNeutralShower loop
181
182 return NOERROR;
183}
184
185//------------------
186// erun
187//------------------
188jerror_t DNeutralParticleHypothesis_factory::erun(void)
189{
190 return NOERROR;
191}
192
193//------------------
194// fini
195//------------------
196jerror_t DNeutralParticleHypothesis_factory::fini(void)
197{
198 return NOERROR;
199}
200
201#define DELTA(i,j)((i==j) ? 1 : 0) ((i==j) ? 1 : 0)
202void DNeutralParticleHypothesis_factory::Calc_Variances(const DNeutralShower *locNeutralShower, double locParticleEnergyUncertainty, DMatrixDSym &locVariances){
203 DLorentzVector locShowerPositionUncertainties = locNeutralShower->dSpacetimeVertexUncertainties;
204
205 // create the simplest error matrix:
206 // At this point, it is assumed that error matrix of measured quantities is diagonal,
207 // with elements like: sigma_Z_t = L/sqrt(12) sigma_X_t = sigma_Y_t = r0/2
208 // L=target length, r0 = target radius...
209 // This means that energy-depth-polar angle relation is neglected.
210 // the order of sigmas is: x_c, y_c, z_c, E, x_t, y_t, z_t
211
212 locVariances.Clear();
213 locVariances.ResizeTo(7, 7);
214
215 locVariances(0,0) = pow(locShowerPositionUncertainties.X(), 2.0);
216 locVariances(1,1) = pow(locShowerPositionUncertainties.Y(), 2.0);
217 locVariances(2,2) = pow(locShowerPositionUncertainties.Z(), 2.0);
218
219 locVariances[3][3] = pow(locParticleEnergyUncertainty, 2.0);
220
221 locVariances[4][4] = pow(0.5*dTargetRadius, 2.0) ; // x_t, y_t
222 locVariances[5][5] = pow(0.5*dTargetRadius, 2.0) ; // x_t, y_t
223 locVariances[6][6] = pow(dTargetLength/sqrt(12.0), 2.0) ; // z_t
224}
225
226void DNeutralParticleHypothesis_factory::Build_ErrorMatrix(const DVector3 &locPathVector, double locEnergy, const DMatrixDSym& locVariances, DMatrixDSym& locErrorMatrix)
227{
228 unsigned int loc_i, loc_j, loc_ik, loc_jk;
229 double R = locPathVector.Mag(); //path length
230 double R2= locPathVector.Mag2(); //path length ^ 2
231 double R3 = R*R2; //path length ^ 3
232 double E_R3 = locEnergy/R3;
233
234 // init and fill rotation matrix, first with momentum derivatives
235 DMatrix A(7, 7);
236 for (loc_i = 0; loc_i < 3; loc_i++) {
237 for (loc_j = 0; loc_j <3; loc_j++) {
238 A[loc_i][loc_j] = E_R3*( R2*DELTA(loc_i, loc_j)((loc_i==loc_j) ? 1 : 0) - locPathVector(loc_i) * locPathVector(loc_j) );
239 A[loc_i][loc_j + 4] = - A[loc_i][loc_j];
240 }
241 }
242
243 // fill energy part and remember: relation between energy and photon position in calorimeter is neglected!
244 A[3][3] = 1.;
245 for (loc_j = 0; loc_j < 3; loc_j++)
246 A[loc_j][3] = locPathVector(loc_j)/R;
247
248 // fill spatial part where: dp_r_x/dp_x_c = - dp_r_x/dp_x_v ....
249 for (loc_i = 0; loc_i < 3; loc_i++) {
250 for ( loc_j = 0; loc_j < 3; loc_j++) {
251 loc_ik = loc_i + 4;
252 loc_jk = loc_j + 4;
253 A[loc_ik][loc_j] = DELTA(loc_i, loc_j)((loc_i==loc_j) ? 1 : 0);
254 A[loc_ik][loc_jk] = - A[loc_ik][loc_j];
255 }
256 }
257
258 locErrorMatrix.ResizeTo(7, 7);
259 locErrorMatrix = locVariances;
260 locErrorMatrix = locErrorMatrix.Similarity(A);
261}
262