Bug Summary

File:libraries/FCAL/DFCALShower_factory.cc
Location:line 258, column 2
Description:Function call argument is an uninitialized value

Annotated Source Code

1//
2// File: DFCALShower_factory.cc
3// Created: Tue May 17 11:57:50 EST 2005
4// Creator: remitche (on Linux mantrid00 2.4.20-18.8smp i686)
5// Edited: B. Schaefer 3/23/2012 removed radiation hard insert functionality
6
7#include <math.h>
8#include <DVector3.h>
9//#include <DLorentzVector.h>
10using namespace std;
11
12#include "DFCALShower_factory.h"
13#include "DFCALGeometry.h"
14//#include "DFCALCluster.h"
15#include "DFCALHit.h"
16#include <JANA/JEvent.h>
17#include <JANA/JApplication.h>
18using namespace jana;
19
20//----------------
21// Constructor
22//----------------
23DFCALShower_factory::DFCALShower_factory()
24{
25
26// Set of coefficients for non-linear energy corrections
27
28 SHOWER_ENERGY_THRESHOLD = 50*k_MeV;
29
30 NON_LIN_COEF_A = 0.439287;
31 NON_LIN_COEF_B = 0.503378;
32 NON_LIN_COEF_C = 1.80842;
33 NON_LIN_COEF_alfa = 1+0.0724789;
34
35 //Loose timing cut
36 SHOWER_TIMING_WINDOW = 5.0;
37
38// Parameters to make shower-depth correction taken from Radphi,
39// slightly modifed to match photon-polar angle
40 FCAL_RADIATION_LENGTH = 3.1;
41 FCAL_CRITICAL_ENERGY = 0.035;
42 FCAL_SHOWER_OFFSET = 1.0;
43
44 gPARMS->SetDefaultParameter("FCAL:SHOWER_ENERGY_THRESHOLD", SHOWER_ENERGY_THRESHOLD);
45
46 gPARMS->SetDefaultParameter("FCAL:NON_LIN_COEF_A", NON_LIN_COEF_A);
47 gPARMS->SetDefaultParameter("FCAL:NON_LIN_COEF_B", NON_LIN_COEF_B);
48 gPARMS->SetDefaultParameter("FCAL:NON_LIN_COEF_C", NON_LIN_COEF_C);
49 gPARMS->SetDefaultParameter("FCAL:NON_LIN_COEF_alfa", NON_LIN_COEF_alfa);
50
51 gPARMS->SetDefaultParameter("FCAL:SHOWER_TIMING_WINDOW", SHOWER_TIMING_WINDOW);
52
53 gPARMS->SetDefaultParameter("FCAL:FCAL_RADIATION_LENGTH", FCAL_RADIATION_LENGTH);
54 gPARMS->SetDefaultParameter("FCAL:FCAL_CRITICAL_ENERGY", FCAL_CRITICAL_ENERGY);
55 gPARMS->SetDefaultParameter("FCAL:FCAL_SHOWER_OFFSET", FCAL_SHOWER_OFFSET);
56
57}
58
59//------------------
60// brun
61//------------------
62// take merging out
63jerror_t DFCALShower_factory::brun(JEventLoop *loop, int runnumber)
64{
65 /*// Get calibration constants
66 map<string, double> cluster_merging;
67 loop->GetCalib("FCAL/cluster_merging", cluster_merging);
68 if(cluster_merging.find("MIN_CLUSTER_SEPARATION")!=cluster_merging.end()){
69 MIN_CLUSTER_SEPARATION = cluster_merging["MIN_CLUSTER_SEPARATION"];
70 if(debug_level>0)jout<<"MIN_CLUSTER_SEPARATION = "<<MIN_CLUSTER_SEPARATION<<endl;
71 }else{
72 jerr<<"Unable to get from MIN_CLUSTER_SEPARATION FCAL/cluster_merging in Calib database!"<<endl;
73 loop->GetJApplication()->Quit();
74 }*/
75
76 // Get calibration constants
77 FCAL_C_EFFECTIVE=15.;
78 map<string, double> fcal_parms;
79 loop->GetCalib("FCAL/fcal_parms", fcal_parms);
80 if (fcal_parms.find("FCAL_C_EFFECTIVE")!=fcal_parms.end()){
81 FCAL_C_EFFECTIVE = fcal_parms["FCAL_C_EFFECTIVE"];
82 if(debug_level>0)jout<<"FCAL_C_EFFECTIVE = "<<FCAL_C_EFFECTIVE<<endl;
83 } else {
84 jerr<<"Unable to get FCAL_C_EFFECTIVE from FCAL/fcal_parms in Calib database!"<<endl;
85 }
86
87 m_zTarget=65.0*k_cm;
88 DApplication *dapp = dynamic_cast<DApplication*>(loop->GetJApplication());
89 if (dapp) {
90 const DGeometry *geom = dapp->GetDGeometry(runnumber);
91 geom->GetTargetZ(m_zTarget);
92 }
93
94 return NOERROR;
95}
96
97
98//------------------
99// evnt
100//------------------
101jerror_t DFCALShower_factory::evnt(JEventLoop *eventLoop, int eventnumber)
102{
103 vector<const DFCALCluster*> fcalClusters;
104 eventLoop->Get(fcalClusters);
105
106 //------------- the following is a temporary kludge...
107// const DVertex *vertex=NULL;
108
109 //loop->GetSingle(vertex);
110// vector<const DVertex *>vertices;
111// eventLoop->Get(vertices);
112// if (vertices.size()) vertex=vertices[0];
113
114 // Return immediately if there isn't even one cluster
115 if(fcalClusters.size()<1)return NOERROR;
116
117// ------------------------------------------------
118/*
119 // Here we try and merge FCAL clusters simply by looking to see if they are
120 // within a certain distance of one another. We do this by first looping
121 // through the list of clusters as many times as is needed until we get
122 // a set of lists for which all clusters are at least MIN_CLUSTER_SEPARATION
123 // from all clusters in all other lists.
124
125 // Start by filling out out lists with one cluster each
126 vector<vector<const DFCALCluster*> > merge_lists;
127 for ( unsigned int i = 0; i < fcalClusters.size(); i++ ) {
128 vector<const DFCALCluster*> merge_list;
129 merge_list.push_back(fcalClusters[i]);
130 merge_lists.push_back(merge_list);
131 }
132
133 // Loop until we find no more lists to merge
134 bool clusters_were_merged;
135 do{
136 clusters_were_merged = false;
137
138 // Loop over all pairs of cluster lists
139 for(unsigned int i=0; i<merge_lists.size(); i++){
140 vector<const DFCALCluster*> &merge_list1 = merge_lists[i];
141 for(unsigned int j=i+1; j<merge_lists.size(); j++){
142 vector<const DFCALCluster*> &merge_list2 = merge_lists[j];
143
144 // Loop over all elements of both lists to see if any are
145 // within MIN_CLUSTER_SEPARATION.
146 for(unsigned int k=0; k<merge_list1.size(); k++){
147 DVector3 pos1 = merge_list1[k]->getCentroid();
148 for(unsigned int m=0; m<merge_list2.size(); m++){
149 DVector3 pos2 = merge_list2[m]->getCentroid();
150
151 double separation_xy = (pos1-pos2).Perp();
152 if(separation_xy<MIN_CLUSTER_SEPARATION){
153 // Phew! if we got here then we need to merge the 2 clusters
154 // The easiest way to do this is to just add all clusters
155 // from merge_list2 to merge_list1 and clear merge_list2.
156 // Then, we ignore empty lists below.
157 merge_list1.insert(merge_list1.end(), merge_list2.begin(), merge_list2.end());
158 merge_list2.clear();
159 clusters_were_merged = true;
160 }
161 }
162 if(clusters_were_merged)break; // we'll need to do the outer "do" loop again anyway so bail now
163 }
164 if(clusters_were_merged)break; // we'll need to do the outer "do" loop again anyway so bail now
165 }
166 if(clusters_were_merged)break; // we'll need to do the outer "do" loop again anyway so bail now
167 }
168
169 }while(clusters_were_merged);
170
171 // Now we loop over the lists of clusters to merge and make a
172 // DFCALShower from the list. Note that it may well be that each
173 // list is still only 1 element long!
174 for ( unsigned int i = 0; i < merge_lists.size(); i++ ) {
175 vector<const DFCALCluster*> &merge_list = merge_lists[i];
176 if(merge_list.size()<1)continue; // ignore empty lists (see comment above)
177
178 DFCALShower* fcalPhoton = makePhoton( merge_list, vertex );
179
180 if ( fcalPhoton->getEnergy() <= 0 ) {
181 cout << "Deleting fcalPhoton " << endl;
182 delete fcalPhoton;
183 continue;
184 }else {
185 _data.push_back(fcalPhoton);
186 }
187 }
188*/
189// --------------------------
190 for( vector< const DFCALCluster* >::const_iterator clItr = fcalClusters.begin();
191 clItr != fcalClusters.end(); ++clItr ){
192
193 DFCALShower* fcalShower = makeFcalShower( *clItr );
194
195 //Make a very loose timing cut here to eliminate out-of-time EM bkgd.
196 //Photon travels dist1 (from center of target to shower maximum
197 //position) at the speed of light and the remainder of the distance,
198 //dist2, to the back of the FCAL (where timing is measured) at c_eff.
199 //This travel time is the expected_time. Since this time is only used
200 //to make a very loose timing cut, we can ignore other factors such as
201 //the fact that not all particles originate exactly at the center of
202 //the target.
203 //Later analysis should use a tighter timing cut incoporating vertex
204 //position and time.
205 DVector3 pos = fcalShower->getPosition();
206 DVector3 vertex(0.0, 0.0, m_zTarget);
207 double dist1 = (pos-vertex).Mag();
208 double dist2 = DFCALGeometry::fcalFaceZ() + DFCALGeometry::blockLength() - pos.Z();
209 double expected_time = dist1/SPEED_OF_LIGHT29.9792 + dist2/FCAL_C_EFFECTIVE;
210
211 if ( fcalShower->getEnergy() <= 0 ) {
212 cout << "Deleting fcalShower " << endl;
213 delete fcalShower;
214 continue;
215 } else if ( fabs(fcalShower->getTime() - expected_time) > SHOWER_TIMING_WINDOW ){
216 delete fcalShower;
217 continue;
218 } else {
219 _data.push_back(fcalShower);
220 }
221
222 }
223
224 return NOERROR;
225}
226
227
228//--------------------------------
229// makePhoton
230//--------------------------------
231DFCALShower* DFCALShower_factory::makeFcalShower( const DFCALCluster* cluster )
232{
233
234 // Loop over list of DFCALCluster objects and calculate the "Non-linear" corrected
235 // energy and position for each. We'll use a logarithmic energy-weighting to
236 // find the final position and error.
237
238// Use target center if vertex does not exist
239 DVector3 target(0.0, 0.0, m_zTarget);
240
241 double cTime = cluster->getTime();
242
243 double errX = cluster->getRMS_x();
244 double errY = cluster->getRMS_y();
245 double errZ; // will be filled by call to GetCorrectedEnergyAndPosition()
246
247 // Get corrected energy, position, and errZ
248 double Ecorrected;
249 DVector3 pos_corrected;
250 GetCorrectedEnergyAndPosition( cluster , Ecorrected, pos_corrected, errZ, &target);
1
Calling 'GetCorrectedEnergyAndPosition'
→
17
←
Returning from 'GetCorrectedEnergyAndPosition'
→
251
252
253 // Make the DFCALShower object
254 DFCALShower* shower = new DFCALShower;
255
256 shower->setEnergy( Ecorrected );
257 shower->setPosition( pos_corrected );
258 shower->setPosError( errX, errY, errZ );
18
←
Function call argument is an uninitialized value
259 shower->setTime ( cTime );
260
261 shower->AddAssociatedObject(cluster);
262
263 return shower;
264}
265
266//--------------------------------
267// GetCorrectedEnergyAndPosition
268//
269// Non-linear and depth corrections should be fixed within DFCALShower member functions
270//--------------------------------
271void DFCALShower_factory::GetCorrectedEnergyAndPosition(const DFCALCluster* cluster, double &Ecorrected, DVector3 &pos_corrected, double &errZ, const DVector3 *vertex)
272{
273// Non-linar energy correction are done here
274 int MAXITER = 1000;
275
276 DVector3 posInCal = cluster->getCentroid();
2
←
Calling 'getCentroid'
→
3
←
Returning from 'getCentroid'
→
277 float x0 = posInCal.Px();
4
←
Calling 'Px'
→
5
←
Returning from 'Px'
→
278 float y0 = posInCal.Py();
6
←
Calling 'Py'
→
7
←
Returning from 'Py'
→
279
280 double Eclust = cluster->getEnergy();
8
←
Calling 'getEnergy'
→
9
←
Returning from 'getEnergy'
→
281
282 double A = NON_LIN_COEF_A;
283 double B = NON_LIN_COEF_B;
284 double C = NON_LIN_COEF_C;
285 double alfa = NON_LIN_COEF_alfa;
286
287
288 double Egamma = 0.;
289
290 if ( A > 0 ) {
10
←
Taking true branch
→
291
292 Egamma = Eclust/A;
293
294 for ( int niter=0; 1; niter++) {
11
←
Loop condition is true. Entering loop body
→
295
296 double energy = Egamma;
297 double non_lin_part = pow(Egamma,1+alfa)/(B+C*Egamma);
298 Egamma = Eclust/A - non_lin_part;
299 if ( fabs( (Egamma-energy)/energy ) < 0.001 ) {
12
←
Taking true branch
→
300
301 break;
13
←
 Execution continues on line 324
→
302
303 }
304 else if ( niter > MAXITER ) {
305
306 cout << " Iteration failed for cluster energy " << Eclust << endl;
307 Egamma = 0;
308
309 break;
310
311 }
312
313 }
314
315 }
316 else {
317
318 cout << "Warning: DFCALShower : parameter A" << A << " is not valid" << endl;
319 }
320
321
322// then depth corrections
323
324 if ( Egamma > 0 ) {
14
←
Taking false branch
→
325 float xV = vertex->X();
326 float yV = vertex->Y();
327 float zV = vertex->Z();
328
329
330 double z0 = DFCALGeometry::fcalFaceZ() - zV;
331 double zMax = (FCAL_RADIATION_LENGTH*(
332 FCAL_SHOWER_OFFSET + log(Egamma/FCAL_CRITICAL_ENERGY)));
333 double zed = z0;
334 double zed1 = z0 + zMax;
335
336 double r0 = sqrt( (x0-xV)*(x0-xV) + (y0-yV)*(y0-yV) );
337
338 int niter;
339 for ( niter=0; niter<100; niter++) {
340
341 double tt = r0/zed1;
342 zed = z0 + zMax/sqrt( 1 + tt*tt );
343 if ( fabs( (zed-zed1) ) < 0.001) {
344 break;
345 }
346 zed1 = zed;
347 }
348
349 posInCal.SetZ( zed + zV );
350 errZ = zed - zed1;
351 }
352
353 Ecorrected = Egamma;
354 pos_corrected = posInCal;
15
←
Calling 'operator='
→
16
←
Returning from 'operator='
→
355}
356
357
358