Bug Summary

File:libraries/PID/DVertex_factory.cc
Location:line 285, column 62
Description:Dereference of null pointer

Annotated Source Code

1// $Id$
2//
3// File: DVertex_factory.cc
4// Created: Tue Apr 6 17:01:54 EDT 2010
5// Creator: davidl (on Darwin Amelia.local 9.8.0 i386)
6//
7
8
9#include <iostream>
10#include <iomanip>
11using namespace std;
12
13#include <TROOT.h>
14#include <TMath.h>
15#include "DVertex_factory.h"
16using namespace jana;
17
18//------------------
19// init
20//------------------
21jerror_t DVertex_factory::init(void)
22{
23 return NOERROR;
24}
25
26//------------------
27// brun
28//------------------
29jerror_t DVertex_factory::brun(jana::JEventLoop *loop, int runnumber)
30{
31
32 // Get Target parameters from XML
33 DApplication *locApplication = dynamic_cast<DApplication*> (loop->GetJApplication());
34 DGeometry *locGeometry = locApplication ? locApplication->GetDGeometry(runnumber):NULL__null;
35 dTargetCenter_Z = 0.0;
36 if(locGeometry)
37 locGeometry->GetTargetZ(dTargetCenter_Z);
38
39 // Configuration parameters
40 GROUP_NUM_SIGMAS_TIME = 100.0; // originally was 3, but changed as temporary measure
41 GROUP_NUM_SIGMAS_Z = 100.0; // originally was 3, but changed as temporary measure
42 DEBUG_HISTS = false;
43 gPARMS->SetDefaultParameter("PID:GROUP_NUM_SIGMAS_TIME", GROUP_NUM_SIGMAS_TIME, "Number of sigmas particles (tracks and or photons) can be apart in time and still be associated with same vertex");
44 gPARMS->SetDefaultParameter("PID:GROUP_NUM_SIGMAS_Z" , GROUP_NUM_SIGMAS_Z , "Number of sigmas particles (tracks and or photons) can be apart in z and still be associated with same vertex");
45 gPARMS->SetDefaultParameter("PID:DEBUG_HISTS" , DEBUG_HISTS , "Enable generation and filling if PID related debugging histos");
46
47 // Specify the size of the vertexInfo pool.
48 // This is not a strict limit. Rather, it is the size the pool will be reduced
49 // to if it has grown larger on the previous event and the current event does
50 // not require this many. This prevents memory-leak-like behavior when running
51 // many threads where each thread keeps allocating bigger and bigger pools as
52 // it comes across slightly busier and busier events.
53 MAX_VERTEXINFOS = 10;
54
55 Nbinst = 600;
56 tmin = -100.0;
57 tmax = 500.0;
58 Nbinsz = 50;
59 zmin = 0.0;
60 zmax = 100.0;
61
62 if (DEBUG_HISTS){
63 locApplication->Lock();
64
65 Nsigmas_t_particles=(TH1F *)gROOT->FindObject("Nsigmas_t_particles");
66 if (!Nsigmas_t_particles) Nsigmas_t_particles=new TH1F("Nsigmas_t_particles","Number of sigmas in t (AssignParticlesToGroups)",5000, 0.0 , 500.0);
67
68 Nsigmas_z_particles=(TH1F *)gROOT->FindObject("Nsigmas_z_particles");
69 if (!Nsigmas_z_particles) Nsigmas_z_particles=new TH1F("Nsigmas_z_particles","Number of sigmas in z (AssignParticlesToGroups)",5000, 0.0 , 500.0);
70
71 locApplication->Unlock();
72 }
73
74 return NOERROR;
75}
76
77//------------------
78// evnt
79//------------------
80jerror_t DVertex_factory::evnt(JEventLoop *loop, int eventnumber)
81{
82 unsigned int loc_i;
83 vector<const DChargedTrack *> locChargedTracks;
84 loop->Get(locChargedTracks);
85
86 // To minimize memory usage and time in allocation, we maintain a
87 // pool of vertexInfo_t objects. Make sure the pool is large enough to hold
88 // all of the particles we have in this event.
89 for(loc_i = dVertexInfoPool.size(); loc_i < locChargedTracks.size(); loc_i++){
90 vertexInfo_t *locVertexInfo = new vertexInfo_t();
91 locVertexInfo->SetLimits(tmin, tmax, zmin, zmax, Nbinst, Nbinsz);
92 dVertexInfoPool.push_back(locVertexInfo);
93 }
94
95 // Periodically delete some vertexInfo_t objects if the pool gets too large.
96 // This prevents memory-leakage-like behavor.
97 if((locChargedTracks.size() < MAX_VERTEXINFOS) && (dVertexInfoPool.size() > MAX_VERTEXINFOS)){
98 for(loc_i = MAX_VERTEXINFOS; loc_i < dVertexInfoPool.size(); loc_i++) delete dVertexInfoPool[loc_i];
99 dVertexInfoPool.resize(MAX_VERTEXINFOS);
100 }
101
102 // Make list of vertices
103 MakeVertices(locChargedTracks);
104
105 // If no vertices but neutral showers present, create vertex at center of target
106 if(_data.size() == 0){
107 vector<const DNeutralShower *> locNeutralShowers;
108 loop->Get(locNeutralShowers);
109 if(locNeutralShowers.size() > 0){ //use center of target
110 DVertex *locVertex = new DVertex;
111 locVertex->dTimeUncertainty = 0.0;
112 DLorentzVector locSpacetimeVertex(0.0, 0.0, dTargetCenter_Z, 0.0);
113 locVertex->dSpacetimeVertex = locSpacetimeVertex;
114 _data.push_back(locVertex);
115 }
116 }
117
118 return NOERROR;
119}
120
121//------------------
122// erun
123//------------------
124jerror_t DVertex_factory::erun(void)
125{
126 return NOERROR;
127}
128
129//------------------
130// fini
131//------------------
132jerror_t DVertex_factory::fini(void)
133{
134 return NOERROR;
135}
136
137// Form vertices from grouping charged particle tracks together according to
138// proximity in time and position of the closest approach to the beam line
139jerror_t DVertex_factory::MakeVertices(vector<const DChargedTrack*> &locChargedTracks){
140 unsigned int loc_i, loc_j;
141 // Vector to hold list of vertexInfo_t objects for all charged tracks
142 vector<vertexInfo_t*> locVertexInfos;
143
144 // Assign and fill vertexInfo_t objects for each charged track
145 for(loc_i = 0; loc_i < locChargedTracks.size(); loc_i++){
146 vertexInfo_t *locVertexInfo = dVertexInfoPool[locVertexInfos.size()];
147 // Use the fit result with the highest figure of merit
148 FillVertexInfoChargedTrack(locVertexInfo, locChargedTracks[loc_i]);
149 locVertexInfos.push_back(locVertexInfo);
150 }
151 // Group tracks together
152 vector< vector<vertexInfo_t *> > locVertexInfoGroups;
153 AssignParticlesToGroups(locVertexInfos, locVertexInfoGroups);
154
155 // OK, we've now grouped the particles together into groups. Create a new
156 // DVertex for each group
157 for (loc_i = 0; loc_i < locVertexInfoGroups.size(); loc_i++){
158 vector<vertexInfo_t *> &locVertexInfoGroup = locVertexInfoGroups[loc_i];
159
160 DVertex *locVertex = new DVertex;
161 locVertex->locCovarianceMatrix.ResizeTo(3,3);
162
163 // Simply average POCA to beamline for all tracks
164 DVector3 temp;
165 double t0 = 0.;
166 double sum_invar = 0.0, invar = 0.0;
167 for(loc_j = 0; loc_j < locVertexInfoGroup.size(); loc_j++){
168 invar = 1.0/(locVertexInfoGroup[loc_j]->sigmat*locVertexInfoGroup[loc_j]->sigmat);
169 sum_invar += invar;
170 t0 += locVertexInfoGroup[loc_j]->t*invar;
171 temp += locVertexInfoGroup[loc_j]->dChargedTrack->dChargedTrackHypotheses[0]->position();
172 }
173 t0 *= 1.0/sum_invar;
174 temp *= 1.0/double(locVertexInfoGroup.size());
175 locVertex->dSpacetimeVertex.SetVect(temp);
176 locVertex->dSpacetimeVertex.SetT(t0);
177 locVertex->dTimeUncertainty = 0.; // <------ this needs to be fixed
178
179 // Add list of tracks used to create this vertex
180 for(loc_j = 0; loc_j < locVertexInfoGroup.size(); loc_j++)
181 locVertex->dChargedTracks.push_back(locVertexInfoGroup[loc_j]->dChargedTrack);
182 _data.push_back(locVertex);
183 }
184
185 return NOERROR;
186}
187
188//------------------
189// FillVertexInfoChargedTrack
190//------------------
191void DVertex_factory::FillVertexInfoChargedTrack(DVertex_factory::vertexInfo_t *locVertexInfo, const DChargedTrack *locChargedTrack)
192{
193 locVertexInfo->Reset();
194 locVertexInfo->dChargedTrack = locChargedTrack;
195 const DChargedTrackHypothesis* locChargedTrackHypothesis = locChargedTrack->dChargedTrackHypotheses[0];
196
197 locVertexInfo->t = locChargedTrackHypothesis->t0();
198 locVertexInfo->sigmat = locChargedTrackHypothesis->t0_err();
199 locVertexInfo->z = locChargedTrackHypothesis->z();
200 locVertexInfo->sigmaz = 0.8/sin(locChargedTrackHypothesis->momentum().Theta()); // in cm. For now, use 3mm wide angle track resolution scaled by sin(theta)
201
202 locVertexInfo->Fill(locVertexInfo->t, locVertexInfo->sigmat, locVertexInfo->z, locVertexInfo->sigmaz);
203}
204
205// Group particles together by time and z position
206void DVertex_factory::AssignParticlesToGroups(vector<vertexInfo_t*> &locVertexInfos, vector< vector<vertexInfo_t *> > &locVertexInfoGroups){
207 unsigned int loc_i;
208
209 // Each particle has a histogram of t vs.z
210 // values filled using approriate uncertainties (no covariance). We
211 // can now use this list to identify resonances in the t/z plane which
212 // indicate a vertex location. Particles within 3sigma in both t and
213 // z of the resonance will be grouped together as belonging to the
214 // same vertex. We loop until all particles have been assigned
215 // to a group, even if that means assigning particles to their own
216 // "group of one".
217
218 // Loop until all particles have been assigned to a group.
219 while(!AllInGroups(locVertexInfos)){
1
Loop condition is true. Entering loop body
→
7
←
Loop condition is true. Entering loop body
→
13
←
Loop condition is true. Entering loop body
→
19
←
Loop condition is true. Entering loop body
→
220 // Make a list of all particles that have not been assigned to a group
221 vector<const DHoughFind*> locUnassignedHoughs;
222 for(loc_i = 0; loc_i < locVertexInfos.size(); loc_i++){
2
←
Loop condition is false. Execution continues on line 227
→
8
←
Loop condition is false. Execution continues on line 227
→
14
←
Loop condition is false. Execution continues on line 227
→
20
←
Loop condition is false. Execution continues on line 227
→
223 if(!locVertexInfos[loc_i]->is_in_group) locUnassignedHoughs.push_back(locVertexInfos[loc_i]);
224 }
225 // Find the maximum t,z coordinate by adding all unassigned
226 // particle's histos together
227 DVector2 maxloc = DHoughFind::GetMaxBinLocation(locUnassignedHoughs);
228
229 if(debug_level>0)_DBG_std::cerr<<"DVertex_factory.cc"<<":"<<229<<
" "<<"Location of maximum: t="<<maxloc.X()<<" z="<<maxloc.Y()<<endl;
3
←
Taking false branch
→
9
←
Taking false branch
→
15
←
Taking false branch
→
21
←
Taking false branch
→
230
231 // Loop over all unassigned particles, assigning any within
232 // 3 sigma in both t and z to the new group. We loop over
233 // the locVertexInfos vector just because it saves a dynamic_cast
234 // if we were to use the unassigned vector.
235 vector<vertexInfo_t *> locVertexInfoGroup;
236 for(loc_i = 0; loc_i < locVertexInfos.size(); loc_i++){
4
←
Loop condition is false. Execution continues on line 265
→
10
←
Loop condition is false. Execution continues on line 265
→
16
←
Loop condition is false. Execution continues on line 265
→
22
←
Loop condition is false. Execution continues on line 265
→
237 vertexInfo_t *locVertexInfo = locVertexInfos[loc_i];
238 if (locVertexInfo->is_in_group) continue;
239 if (locVertexInfo->is_matched_to_vertex) continue;
240
241 double delta_t = fabs(maxloc.X() - locVertexInfo->t);
242 double Nsigmas_t = delta_t/locVertexInfo->sigmat;
243 if(DEBUG_HISTS) Nsigmas_t_particles->Fill(Nsigmas_t);
244 if(Nsigmas_t > GROUP_NUM_SIGMAS_TIME) continue;
245
246 double delta_z = fabs(maxloc.Y() - locVertexInfo->z);
247 double Nsigmas_z = delta_z/locVertexInfo->sigmaz;
248 if(DEBUG_HISTS) Nsigmas_z_particles->Fill(Nsigmas_z);
249 if(Nsigmas_z > GROUP_NUM_SIGMAS_Z) continue;
250
251 // Assign this particle to the group
252 locVertexInfo->is_in_group = true;
253 locVertexInfoGroup.push_back(locVertexInfo);
254 }
255
256 // At this point it's possible (but hopefully unlikely) that the
257 // maximum in the t,z sum histo was generated at an in-between place
258 // with no single particle nearby. In that case, the new_group is
259 // empty, even though there are unassigned particles. The best we
260 // can do here is to assign one particle to the new_group and hope
261 // that the next iteration groups the remaining ones appropriately.
262 // To try and minimize the chances of placing a particle from the
263 // L1 trigger event in its own group, we choose the particle with a time
264 // the furthest away from t=0.
265 if(locVertexInfoGroup.size() == 0){
5
←
Taking false branch
→
11
←
Taking false branch
→
17
←
Taking false branch
→
23
←
Taking true branch
→
266 vertexInfo_t *vi_with_max_t = NULL__null;
267 double delta_t_max = 0.0;
268 for(loc_i = 0; loc_i < locVertexInfos.size(); loc_i++){
24
←
Loop condition is true. Entering loop body
→
27
←
Loop condition is false. Execution continues on line 277
→
269 vertexInfo_t *locVertexInfo = locVertexInfos[loc_i];
270 if(locVertexInfo->is_in_group) continue;
25
←
Taking false branch
→
271 double delta_t = fabs(maxloc.X() - locVertexInfo->t);
272 if(delta_t>delta_t_max || vi_with_max_t==NULL__null){
26
←
Taking true branch
→
273 delta_t_max = delta_t;
274 vi_with_max_t = locVertexInfo;
275 }
276 }
277 if(vi_with_max_t == NULL__null){
28
←
Taking false branch
→
278 _DBG_std::cerr<<"DVertex_factory.cc"<<":"<<278<<
" "<<"vi_with_max_t==NULL. This should never happen! Complain to davidl@jlab.org"<<endl;
279 break;
280 }
281 locVertexInfoGroup.push_back(vi_with_max_t);
282 }
283
284 // Set the is_in_group flags for all of the members of the new group
285 for(loc_i = 0; loc_i < locVertexInfoGroup.size(); loc_i++) locVertexInfoGroup[loc_i]->is_in_group = true;
6
←
Loop condition is false. Execution continues on line 288
→
12
←
Loop condition is false. Execution continues on line 288
→
18
←
Loop condition is false. Execution continues on line 288
→
29
←
Loop condition is true. Entering loop body
→
30
←
Dereference of null pointer
286
287 // Add the new group to the list of groups
288 locVertexInfoGroups.push_back(locVertexInfoGroup);
289 }
290}
291
292//------------------
293// AllInGroups
294//------------------
295bool DVertex_factory::AllInGroups(vector<vertexInfo_t*> &locVertexInfos)
296{
297 for(unsigned int loc_i = 0; loc_i < locVertexInfos.size(); loc_i++){
298 if(!locVertexInfos[loc_i]->is_in_group) return false;
299 }
300 return true;
301}
302