Bug Summary

File:libraries/PID/DNeutralParticleHypothesis_factory.cc
Location:line 34, column 2
Description:Called C++ object pointer is null

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, int32_t 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;
1
Assuming 'locApplication' is null
→
2
←
'?' condition is false
→
3
←
'locGeometry' initialized to a null pointer value
→
34 locGeometry->GetTargetZ(dTargetCenterZ);
4
←
Called C++ object pointer is null
35 locEventLoop->GetSingle(dParticleID);
36
37 return NOERROR;
38}
39
40//------------------
41// evnt
42//------------------
43jerror_t DNeutralParticleHypothesis_factory::evnt(jana::JEventLoop *locEventLoop, uint64_t eventnumber)
44{
45 vector<const DNeutralShower*> locNeutralShowers;
46 locEventLoop->Get(locNeutralShowers);
47
48 vector<Particle_t> locPIDHypotheses;
49 locPIDHypotheses.push_back(Gamma);
50 locPIDHypotheses.push_back(Neutron);
51
52 const DEventRFBunch* locEventRFBunch = NULL__null;
53 locEventLoop->GetSingle(locEventRFBunch);
54
55 const DVertex* locVertex = NULL__null;
56 locEventLoop->GetSingle(locVertex);
57
58 // Loop over DNeutralShowers
59 for(size_t loc_i = 0; loc_i < locNeutralShowers.size(); ++loc_i)
60 {
61 const DNeutralShower *locNeutralShower = locNeutralShowers[loc_i];
62 // Loop over vertices and PID hypotheses & create DNeutralParticleHypotheses for each combination
63 for(size_t loc_j = 0; loc_j < locPIDHypotheses.size(); ++loc_j)
64 {
65 DNeutralParticleHypothesis* locNeutralParticleHypothesis = Create_DNeutralParticleHypothesis(locEventLoop, locNeutralShower, locPIDHypotheses[loc_j], locEventRFBunch, locVertex);
66 if(locNeutralParticleHypothesis != NULL__null)
67 _data.push_back(locNeutralParticleHypothesis);
68 }
69 }
70
71 return NOERROR;
72}
73
74DNeutralParticleHypothesis* DNeutralParticleHypothesis_factory::Create_DNeutralParticleHypothesis(JEventLoop *locEventLoop, const DNeutralShower* locNeutralShower, Particle_t locPID, const DEventRFBunch* locEventRFBunch, const DVertex* locVertex) const
75{
76 DVector3 locVertexGuess = locVertex->dSpacetimeVertex.Vect();
77 double locStartTime = locVertex->dSpacetimeVertex.T();
78
79 double locHitTime = locNeutralShower->dSpacetimeVertex.T();
80 DVector3 locHitPoint = locNeutralShower->dSpacetimeVertex.Vect();
81
82 // Calculate DNeutralParticleHypothesis Quantities (projected time at vertex for given id, etc.)
83 double locMass = ParticleMass(locPID);
84
85 DVector3 locPath = locHitPoint - locVertexGuess;
86 double locPathLength = locPath.Mag();
87 if(!(locPathLength > 0.0))
88 return NULL__null; //invalid, will divide by zero when creating error matrix, so skip!
89
90 DVector3 locMomentum(locPath);
91
92 uint64_t locEventNumber = locEventLoop->GetJEvent().GetEventNumber();
93 TMatrixFSym* locParticleCovariance = (dynamic_cast<DApplication*>(japp))->Get_CovarianceMatrixResource(7, locEventNumber);
94
95 double locProjectedTime = 0.0, locPMag = 0.0;
96 if(locPID != Gamma)
97 {
98 double locDeltaT = locHitTime - locStartTime;
99 double locBeta = locPathLength/(locDeltaT*29.9792458);
100 if(locBeta >= 1.0)
101 locBeta = 0.9999;
102 if(locBeta < 0.0)
103 locBeta = 0.0;
104 double locGamma = 1.0/sqrt(1.0 - locBeta*locBeta);
105 locPMag = locGamma*locBeta*locMass;
106 locMomentum.SetMag(locPMag);
107 locProjectedTime = locStartTime + (locVertexGuess.Z() - dTargetCenterZ)/29.9792458;
108 Calc_ParticleCovariance_Massive(locNeutralShower, locVertex, locMass, locDeltaT, locMomentum, locPath, locParticleCovariance);
109 }
110 else
111 {
112 locPMag = locNeutralShower->dEnergy;
113 double locFlightTime = locPathLength/29.9792458;
114 locProjectedTime = locHitTime - locFlightTime;
115 locMomentum.SetMag(locPMag);
116 Calc_ParticleCovariance_Photon(locNeutralShower, locVertex, locMomentum, locPath, locParticleCovariance);
117 }
118
119 // Build DNeutralParticleHypothesis // dEdx not set
120 DNeutralParticleHypothesis* locNeutralParticleHypothesis = new DNeutralParticleHypothesis;
121 locNeutralParticleHypothesis->AddAssociatedObject(locNeutralShower);
122
123 locNeutralParticleHypothesis->dNeutralShowerID = locNeutralShower->dShowerID;
124 locNeutralParticleHypothesis->setPID(locPID);
125 locNeutralParticleHypothesis->setMass(locMass);
126 locNeutralParticleHypothesis->setCharge(0.0);
127 locNeutralParticleHypothesis->setMomentum(locMomentum);
128 locNeutralParticleHypothesis->setPosition(locVertexGuess);
129 locNeutralParticleHypothesis->setT0(locStartTime, locVertex->dCovarianceMatrix(3, 3), locEventRFBunch->dTimeSource);
130 locNeutralParticleHypothesis->setTime(locProjectedTime);
131 locNeutralParticleHypothesis->setT1(locNeutralShower->dSpacetimeVertex.T(), sqrt(locNeutralShower->dCovarianceMatrix(4, 4)), locNeutralShower->dDetectorSystem);
132 locNeutralParticleHypothesis->setPathLength(locPathLength, 0.0); //zero uncertainty (for now)
133 locNeutralParticleHypothesis->setErrorMatrix(locParticleCovariance);
134
135 // Calculate DNeutralParticleHypothesis FOM
136 unsigned int locNDF = 0;
137 double locChiSq = 0.0;
138 double locFOM = -1.0; //undefined for non-photons
139 if(locPID == Gamma)
140 {
141 double locTimePull = 0.0;
142 locChiSq = dParticleID->Calc_TimingChiSq(locNeutralParticleHypothesis, locNDF, locTimePull);
143 locFOM = TMath::Prob(locChiSq, locNDF);
144 }
145 locNeutralParticleHypothesis->dChiSq = locChiSq;
146 locNeutralParticleHypothesis->dNDF = locNDF;
147 locNeutralParticleHypothesis->dFOM = locFOM;
148 return locNeutralParticleHypothesis;
149}
150
151void DNeutralParticleHypothesis_factory::Calc_ParticleCovariance_Photon(const DNeutralShower* locNeutralShower, const DVertex* locVertex, const DVector3& locMomentum, const DVector3& locPathVector, TMatrixFSym* locParticleCovariance) const
152{
153 //build 8x8 matrix: 5x5 shower, 3x3 vertex position
154 TMatrixFSym locShowerPlusVertCovariance(8);
155 for(unsigned int loc_l = 0; loc_l < 5; ++loc_l) //shower: e, x, y, z, t
156 {
157 for(unsigned int loc_m = 0; loc_m < 5; ++loc_m)
158 locShowerPlusVertCovariance(loc_l, loc_m) = locNeutralShower->dCovarianceMatrix(loc_l, loc_m);
159 }
160 for(unsigned int loc_l = 0; loc_l < 3; ++loc_l) //vertex xyz
161 {
162 for(unsigned int loc_m = 0; loc_m < 3; ++loc_m)
163 locShowerPlusVertCovariance(5 + loc_l, 5 + loc_m) = locVertex->dCovarianceMatrix(loc_l, loc_m);
164 }
165
166 //the input delta X is defined as "hit - start"
167 //however, the documentation and derivations define delta x as "start - hit"
168 //so, reverse the sign of the inputs to match the documentation
169 //and then the rest will follow the documentation
170 DVector3 locDeltaX = -1.0*locPathVector;
171 DVector3 locDeltaXOverDeltaXSq = (1.0/locDeltaX.Mag2())*locDeltaX;
172 DVector3 locUnitP = (1.0/locNeutralShower->dEnergy)*locMomentum;
173 DVector3 locUnitDeltaXOverC = (1.0/(29.9792458*locDeltaX.Mag()))*locDeltaX;
174
175 //build transform matrix
176 TMatrix locTransformMatrix(7, 8);
177
178 locTransformMatrix(0, 0) = locUnitP.X(); //partial deriv of px wrst shower-e
179 locTransformMatrix(0, 1) = locMomentum.Px()*(locDeltaXOverDeltaXSq.X() - 1.0/locDeltaX.X()); //partial deriv of px wrst shower-x
180 locTransformMatrix(0, 2) = locMomentum.Px()*locDeltaX.Y()/locDeltaX.Mag2(); //partial deriv of px wrst shower-y
181 locTransformMatrix(0, 3) = locMomentum.Px()*locDeltaX.Z()/locDeltaX.Mag2(); //partial deriv of px wrst shower-z
182 locTransformMatrix(0, 5) = -1.0*locTransformMatrix(0, 1); //partial deriv of px wrst vert-x
183 locTransformMatrix(0, 6) = -1.0*locTransformMatrix(0, 2); //partial deriv of px wrst vert-y
184 locTransformMatrix(0, 7) = -1.0*locTransformMatrix(0, 3); //partial deriv of px wrst vert-z
185
186 locTransformMatrix(1, 0) = locUnitP.Y(); //partial deriv of py wrst shower-e
187 locTransformMatrix(1, 1) = locMomentum.Py()*locDeltaX.X()/locDeltaX.Mag2(); //partial deriv of py wrst shower-x
188 locTransformMatrix(1, 2) = locMomentum.Py()*(locDeltaXOverDeltaXSq.Y() - 1.0/locDeltaX.Y()); //partial deriv of py wrst shower-y
189 locTransformMatrix(1, 3) = locMomentum.Py()*locDeltaX.Z()/locDeltaX.Mag2(); //partial deriv of py wrst shower-z
190 locTransformMatrix(1, 5) = -1.0*locTransformMatrix(1, 1); //partial deriv of py wrst vert-x
191 locTransformMatrix(1, 6) = -1.0*locTransformMatrix(1, 2); //partial deriv of py wrst vert-y
192 locTransformMatrix(1, 7) = -1.0*locTransformMatrix(1, 3); //partial deriv of py wrst vert-z
193
194 locTransformMatrix(2, 0) = locUnitP.Z(); //partial deriv of pz wrst shower-e
195 locTransformMatrix(2, 1) = locMomentum.Pz()*locDeltaX.X()/locDeltaX.Mag2(); //partial deriv of pz wrst shower-x
196 locTransformMatrix(2, 2) = locMomentum.Pz()*locDeltaX.Y()/locDeltaX.Mag2(); //partial deriv of pz wrst shower-y
197 locTransformMatrix(2, 3) = locMomentum.Pz()*(locDeltaXOverDeltaXSq.Z() - 1.0/locDeltaX.Z()); //partial deriv of pz wrst shower-z
198 locTransformMatrix(2, 5) = -1.0*locTransformMatrix(2, 1); //partial deriv of pz wrst vert-x
199 locTransformMatrix(2, 6) = -1.0*locTransformMatrix(2, 2); //partial deriv of pz wrst vert-y
200 locTransformMatrix(2, 7) = -1.0*locTransformMatrix(2, 3); //partial deriv of pz wrst vert-z
201
202 locTransformMatrix(3, 5) = 1.0; //partial deriv of x wrst vertex-x
203 locTransformMatrix(4, 6) = 1.0; //partial deriv of y wrst vertex-y
204 locTransformMatrix(5, 7) = 1.0; //partial deriv of z wrst vertex-z
205
206 locTransformMatrix(6, 0) = 0.0; //partial deriv of t wrst shower-e //beta defined = 1
207 locTransformMatrix(6, 1) = locUnitDeltaXOverC.X(); //partial deriv of t wrst shower-x
208 locTransformMatrix(6, 2) = locUnitDeltaXOverC.Y(); //partial deriv of t wrst shower-y
209 locTransformMatrix(6, 3) = locUnitDeltaXOverC.Z(); //partial deriv of t wrst shower-z
210 locTransformMatrix(6, 4) = 1.0; //partial deriv of t wrst shower-t
211 locTransformMatrix(6, 5) = -1.0*locTransformMatrix(6, 1); //partial deriv of t wrst vert-x
212 locTransformMatrix(6, 6) = -1.0*locTransformMatrix(6, 2); //partial deriv of t wrst vert-y
213 locTransformMatrix(6, 7) = -1.0*locTransformMatrix(6, 3); //partial deriv of t wrst vert-z
214
215 //convert
216 *locParticleCovariance = locShowerPlusVertCovariance.Similarity(locTransformMatrix);
217}
218
219void DNeutralParticleHypothesis_factory::Calc_ParticleCovariance_Massive(const DNeutralShower* locNeutralShower, const DVertex* locVertex, double locMass, double locDeltaT, const DVector3& locMomentum, const DVector3& locPathVector, TMatrixFSym* locParticleCovariance) const
220{
221 //build 9x9 matrix: 5x5 shower, 4x4 vertex position & time
222 TMatrixFSym locShowerPlusVertCovariance(9);
223 for(unsigned int loc_l = 0; loc_l < 5; ++loc_l) //shower: e, x, y, z, t
224 {
225 for(unsigned int loc_m = 0; loc_m < 5; ++loc_m)
226 locShowerPlusVertCovariance(loc_l, loc_m) = locNeutralShower->dCovarianceMatrix(loc_l, loc_m);
227 }
228 for(unsigned int loc_l = 0; loc_l < 4; ++loc_l) //vertex xyzt
229 {
230 for(unsigned int loc_m = 0; loc_m < 4; ++loc_m)
231 locShowerPlusVertCovariance(5 + loc_l, 5 + loc_m) = locVertex->dCovarianceMatrix(loc_l, loc_m);
232 }
233
234 //the input delta X/t are defined as "hit - start"
235 //however, the documentation and derivations define delta x/t as "start - hit"
236 //so, reverse the sign of the inputs to match the documentation
237 //and then the rest will follow the documentation
238 DVector3 locDeltaX = -1.0*locPathVector;
239 locDeltaT *= -1.0;
240 double locCSq = 29.9792458*29.9792458;
241 double locCDeltaTSq = locCSq*locDeltaT*locDeltaT;
242 double locDeltaX4Sq = locCDeltaTSq - locDeltaX.Mag2();
243 DVector3 locDeltaXOverDeltaX4Sq = (1.0/locDeltaX4Sq)*locDeltaX;
244 DVector3 locEPVecOverPSq = (locNeutralShower->dEnergy/locMomentum.Mag2())*locMomentum;
245 DVector3 locEPVecOverCPMagDeltaXMag = (locNeutralShower->dEnergy/(29.9792458*locDeltaX.Mag()*locMomentum.Mag()))*locDeltaX;
246
247 //build transform matrix
248 TMatrix locTransformMatrix(7, 9);
249
250 locTransformMatrix(0, 0) = 0.0; //partial deriv of px wrst shower-e
251 locTransformMatrix(0, 1) = locMomentum.Px()*(locDeltaXOverDeltaX4Sq.X() + 1.0/locDeltaX.X()); //partial deriv of px wrst shower-x
252 locTransformMatrix(0, 2) = locMomentum.Px()*locDeltaX.Y()/locDeltaX4Sq; //partial deriv of px wrst shower-y
253 locTransformMatrix(0, 3) = locMomentum.Px()*locDeltaX.Z()/locDeltaX4Sq; //partial deriv of px wrst shower-z
254 locTransformMatrix(0, 4) = locDeltaT*locCSq*locMomentum.Px()/locDeltaX4Sq; //partial deriv of px wrst shower-t
255 locTransformMatrix(0, 5) = -1.0*locTransformMatrix(0, 1); //partial deriv of px wrst vert-x
256 locTransformMatrix(0, 6) = -1.0*locTransformMatrix(0, 2); //partial deriv of px wrst vert-y
257 locTransformMatrix(0, 7) = -1.0*locTransformMatrix(0, 3); //partial deriv of px wrst vert-z
258 locTransformMatrix(0, 8) = -1.0*locTransformMatrix(0, 4); //partial deriv of px wrst vert-t
259
260 locTransformMatrix(1, 0) = 0.0; //partial deriv of py wrst shower-e
261 locTransformMatrix(1, 1) = locMomentum.Py()*locDeltaX.X()/locDeltaX4Sq; //partial deriv of py wrst shower-x
262 locTransformMatrix(1, 2) = locMomentum.Py()*(locDeltaXOverDeltaX4Sq.Y() + 1.0/locDeltaX.Y()); //partial deriv of py wrst shower-y
263 locTransformMatrix(1, 3) = locMomentum.Py()*locDeltaX.Z()/locDeltaX4Sq; //partial deriv of py wrst shower-z
264 locTransformMatrix(1, 4) = locDeltaT*locCSq*locMomentum.Py()/locDeltaX4Sq; //partial deriv of py wrst shower-t
265 locTransformMatrix(1, 5) = -1.0*locTransformMatrix(1, 1); //partial deriv of py wrst vert-x
266 locTransformMatrix(1, 6) = -1.0*locTransformMatrix(1, 2); //partial deriv of py wrst vert-y
267 locTransformMatrix(1, 7) = -1.0*locTransformMatrix(1, 3); //partial deriv of py wrst vert-z
268 locTransformMatrix(1, 8) = -1.0*locTransformMatrix(1, 4); //partial deriv of py wrst vert-t
269
270 locTransformMatrix(2, 0) = 0.0; //partial deriv of pz wrst shower-e
271 locTransformMatrix(2, 1) = locMomentum.Pz()*locDeltaX.X()/locDeltaX4Sq; //partial deriv of pz wrst shower-x
272 locTransformMatrix(2, 2) = locMomentum.Pz()*locDeltaX.Y()/locDeltaX4Sq; //partial deriv of pz wrst shower-y
273 locTransformMatrix(2, 3) = locMomentum.Pz()*(locDeltaXOverDeltaX4Sq.Z() + 1.0/locDeltaX.Z()); //partial deriv of pz wrst shower-z
274 locTransformMatrix(2, 4) = locDeltaT*locCSq*locMomentum.Pz()/locDeltaX4Sq; //partial deriv of pz wrst shower-t
275 locTransformMatrix(2, 5) = -1.0*locTransformMatrix(2, 1); //partial deriv of pz wrst vert-x
276 locTransformMatrix(2, 6) = -1.0*locTransformMatrix(2, 2); //partial deriv of pz wrst vert-y
277 locTransformMatrix(2, 7) = -1.0*locTransformMatrix(2, 3); //partial deriv of pz wrst vert-z
278 locTransformMatrix(2, 8) = -1.0*locTransformMatrix(2, 4); //partial deriv of pz wrst vert-t
279
280 locTransformMatrix(3, 5) = 1.0; //partial deriv of x wrst vertex-x
281 locTransformMatrix(4, 6) = 1.0; //partial deriv of y wrst vertex-y
282 locTransformMatrix(5, 7) = 1.0; //partial deriv of z wrst vertex-z
283 locTransformMatrix(6, 8) = 1.0; //partial deriv of t wrst vertex-t
284
285 //convert
286 *locParticleCovariance = locShowerPlusVertCovariance.Similarity(locTransformMatrix);
287}