Bug Summary

File:libraries/BCAL/DBCALClump_factory.cc
Location:line 310, column 7
Description:Value stored to 'Mt' is never read

Annotated Source Code

1/*
2 * DBCALClump.h
3 *
4 * Created by Beni Zihlmann Tue Mar 12 2013
5 *
6 */
7
8
9#include <JANA/JApplication.h>
10using namespace jana;
11
12#include <DBCALClump_factory.h>
13
14
15//----------------
16// init
17//----------------
18jerror_t DBCALClump_factory::init(void)
19{
20 return NOERROR;
21}
22
23//----------------
24// init
25//----------------
26jerror_t DBCALClump_factory::brun(JEventLoop *loop, int32_t runnumber)
27{
28 map<string, double> bcalparms;
29
30 if ( !loop->GetCalib("BCAL/mc_parms", bcalparms)){
31 cout<<"DBCALClump_factory: loading values from TOF data base"<<endl;
32 } else {
33 cout << "DBCALClumpo_factory: Error loading values from BCAL MC data base" <<endl;
34
35 VELOCITY = 15.; // set to some reasonable value
36
37 return NOERROR;
38 }
39
40 VELOCITY = bcalparms["C_EFFECTIVE"];
41
42
43 return NOERROR;
44}
45
46
47
48//----------------
49// evnt
50//----------------
51jerror_t DBCALClump_factory::evnt(JEventLoop *loop, int eventnumber) {
52
53 const DBCALHit *BcalMatrixU[48*4][4]; // Up stream
54 const DBCALHit *BcalMatrixD[48*4][4]; // Down stream
55
56 int MatrixUD[48*4][4];
57
58 float EU[48*4];
59 float ED[48*4];
60 int EUi[48*4];
61 int EDi[48*4];
62
63 vector<const DBCALHit*> AllBcalHits;
64 loop->Get(AllBcalHits);
65
66 for (int i=0;i<48*4;i++){
67 EU[i] = 0.;
68 ED[i] = 0.;
69 EUi[i] = 0;
70 EDi[i] = 0;
71 for (int m=0;m<4;m++){
72 BcalMatrixU[i][m] = NULL__null;
73 BcalMatrixD[i][m] = NULL__null;
74 MatrixUD[i][m] = 0;
75 }
76 }
77
78 // setup BCAL Matrix
79
80 for (int i=0; i< (int)AllBcalHits.size();i++){
81 const DBCALHit* hit = AllBcalHits[i];
82 //cout<<"Module: "<<hit->module<<" Sector: "<<hit->sector<<" Layer: "<<hit->layer<<endl;
83 int chan = hit->sector-1 + (hit->module-1)*4;
84 if ((hit->t<200.) && (hit->t>0.)) { // HARD CODED
85 if (hit->end == DBCALGeometry::kUpstream){
86 BcalMatrixU[chan][hit->layer-1] = hit;
87 } else{
88 BcalMatrixD[chan][hit->layer-1] = hit;
89 }
90 }
91 }
92
93 // Build the up/down stream coincidence Matrix
94 for (int i=0;i<48*4;i++){
95 for (int m=0;m<4;m++){
96 if ((BcalMatrixU[i][m]) && (BcalMatrixD[i][m])){
97 float meant = (BcalMatrixU[i][m]->t + BcalMatrixD[i][m]->t - 390./16.75)/2.; // HARD CODED VALUES
98 // this mean time has to be positive!!
99 // and the cut could be made stronger.
100 if (meant>0){
101 MatrixUD[i][m] = 1;
102 } else {
103 BcalMatrixU[i][m] = NULL__null;
104 BcalMatrixD[i][m] = NULL__null;
105 }
106 }
107 }
108 }
109
110
111 // first get rid of single-ended isolated hits
112 for (int i=0;i<48*4;i++){
113 int bef = i-1;
114 int aft = i+1;
115 if (bef<0){
116 bef+=48*4;
117 }
118 if (aft>=48*4){
119 aft -=48*4;
120 }
121 // layer 1;
122 if (!MatrixUD[i][0]){
123 if (BcalMatrixU[i][0]){
124 if ((!BcalMatrixU[bef][0]) &&
125 (!BcalMatrixU[aft][0]) &&
126 (!BcalMatrixU[i][1]) &&
127 (!BcalMatrixU[aft][1]) &&
128 (!BcalMatrixU[aft][1])){
129 BcalMatrixU[i][0] = NULL__null;
130 }
131 }
132 if (BcalMatrixD[i][0]){
133 if ((!BcalMatrixD[bef][0]) &&
134 (!BcalMatrixD[aft][0]) &&
135 (!BcalMatrixD[i][1]) &&
136 (!BcalMatrixD[aft][1]) &&
137 (!BcalMatrixD[aft][1])){
138 BcalMatrixD[i][0] = NULL__null;
139 }
140 }
141 }
142 // layer 2;
143 if (!MatrixUD[i][1]){
144 if (BcalMatrixU[i][1]){
145 if ((!BcalMatrixU[bef][1]) &&
146 (!BcalMatrixU[aft][1]) &&
147 (!BcalMatrixU[i][2]) &&
148 (!BcalMatrixU[aft][2]) &&
149 (!BcalMatrixU[aft][2])){
150 BcalMatrixU[i][1] = NULL__null;
151 }
152 }
153 if (BcalMatrixD[i][1]){
154 if ((!BcalMatrixD[bef][1]) &&
155 (!BcalMatrixD[aft][1]) &&
156 (!BcalMatrixD[i][2]) &&
157 (!BcalMatrixD[aft][2]) &&
158 (!BcalMatrixD[aft][2])){
159 BcalMatrixD[i][1] = NULL__null;
160 }
161 }
162 }
163 // layer 3;
164 if (!MatrixUD[i][2]){
165 if (BcalMatrixU[i][2]){
166 if ((!BcalMatrixU[bef][2]) &&
167 (!BcalMatrixU[aft][2]) &&
168 (!BcalMatrixU[i][3]) &&
169 (!BcalMatrixU[aft][3]) &&
170 (!BcalMatrixU[aft][3])){
171 BcalMatrixU[i][2] = NULL__null;
172 }
173 }
174 if (BcalMatrixD[i][2]){
175 if ((!BcalMatrixD[bef][2]) &&
176 (!BcalMatrixD[aft][2]) &&
177 (!BcalMatrixD[i][3]) &&
178 (!BcalMatrixD[aft][3]) &&
179 (!BcalMatrixD[aft][3])){
180 BcalMatrixD[i][2] = NULL__null;
181 }
182 }
183 }
184 // layer 4;
185 if (!MatrixUD[i][3]){
186 if (BcalMatrixU[i][3]){
187 if ((!BcalMatrixU[bef][3]) &&
188 (!BcalMatrixU[aft][3]) &&
189 (!BcalMatrixU[i][2]) &&
190 (!BcalMatrixU[aft][2]) &&
191 (!BcalMatrixU[aft][2])){
192 BcalMatrixU[i][3] = NULL__null;
193 }
194 }
195 if (BcalMatrixD[i][3]){
196 if ((!BcalMatrixD[bef][3]) &&
197 (!BcalMatrixD[aft][3]) &&
198 (!BcalMatrixD[i][2]) &&
199 (!BcalMatrixD[aft][2]) &&
200 (!BcalMatrixD[aft][2])){
201 BcalMatrixD[i][3] = NULL__null;
202 }
203 }
204 }
205 }
206
207 for (int i=0;i<48*4;i++){
208 for (int m=0;m<4;m++){
209
210 if (BcalMatrixU[i][m]){
211 EU[i] += BcalMatrixU[i][m]->E;
212 EUi[i] += 1;
213 }
214 if (BcalMatrixD[i][m]){
215 ED[i] += BcalMatrixD[i][m]->E;
216 EDi[i] += 1;
217 }
218
219 }
220 }
221
222 // clean up Matrix from obvious stray hits
223 for (int i=1;i<48*4-1;i++) {
224 int a = EDi[i-1]+EUi[i-1];
225 int b = EDi[i] + EUi[i];
226 int c = EDi[i+1]+EUi[i+1];
227 if ((b<2)&&(a==0)&&(c==0)){
228 EDi[i] = 0;
229 EUi[i] = 0;
230 EU[i] = 0;
231 ED[i] = 0;
232 }
233 }
234
235 int CNT=1;
236 while (CNT) {
237
238 //cout<<"Find Next Clump"<<endl;
239 int MaxU=999;
240 int MaxD=999;
241
242 float mU=0.;
243 float mD=0.;
244 // find largest energy deposition in the matrix
245 for (int i=0;i<48*4;i++){
246 if ((EU[i]>0.0)&&(ED[i]>0.0)){
247 if (EU[i]>mU){
248 mU = EU[i];
249 MaxU = i;
250 }
251 if (ED[i]>mD){
252 mD = ED[i];
253 MaxD = i;
254 }
255 }
256 }
257
258 if (MaxU==999 || MaxD==999){
259 return NOERROR;
260 }
261
262 if (mU>mD){
263 MaxD = MaxU;
264 } else {
265 MaxU = MaxD;
266 }
267
268 //cout<<MaxD<<" / "<<MaxU<<endl;
269
270 int mod = MaxD;
271 float Mt = 0;
272 float cnt = 0;
273 float EmaxU = 0.;
274 float EmaxD = 0.;
275 int idxMaxU[2] = {999,999};
276 int idxMaxD[2] = {999,999};
277
278 // loop over all layers in the seed sector
279 // and its neighouring sectors and use all
280 // instances with hits on both ends to determine
281 // the mean time
282 for (int k=mod-1; k<mod+2; k++) {
283 int idx = k;
284 if (k<0){
285 idx += 48*4;
286 } else if (k>=4*48){
287 idx -= 48*4;
288 }
289 for (int n=0;n<3;n++){
290 const DBCALHit* hitD = BcalMatrixD[idx][n];
291 const DBCALHit* hitU = BcalMatrixU[idx][n];
292 if (hitD && hitU){
293 if (hitD->E > EmaxD){
294 EmaxD = hitD->E;
295 idxMaxD[0] = idx; // index for seed
296 idxMaxD[1] = n;
297 }
298 if (hitU->E > EmaxU){
299 EmaxU = hitU->E;
300 idxMaxU[0] = idx; // index for seed
301 idxMaxU[1] = n;
302 }
303 float mt = hitD->t + hitU->t;
304 Mt += mt;
305 cnt +=1.;
306 }
307 }
308 }
309 if (cnt>0){
310 Mt /= cnt; // need improve to use Energy weight
Value stored to 'Mt' is never read
311 }
312
313 //cout<<Mt<<" "<<cnt<<endl;
314
315 if (EmaxU>EmaxD){
316 idxMaxD[0] = idxMaxU[0];
317 idxMaxD[1] = idxMaxU[1];
318 } else {
319 idxMaxU[0] = idxMaxD[0];
320 idxMaxU[1] = idxMaxD[1];
321 }
322
323 // Note: Mt is the flight time of the patricle creating this shower!
324 // This means Mt is the same for up and down stream!.
325 //
326 // Now collect all hits related to the seed hit for up and down stream
327 // separately. Only hits that appear in time (+/- 2ns) with the seed hit
328 // are added to the lists for layers 1,2 and 3. Layer 4 is always added.
329 // This may be changed if the timing information from the FADC is good enough.
330 // All hits in the list regardless if they contribut or not are then removed
331 // from the list of hits so they can not be used twice.
332
333 CNT = (int)cnt;
334
335 if (cnt>0) {
336 vector <const DBCALHit*> HitListU;
337 vector <const DBCALHit*> HitListD;
338 vector <float> MeanTime;
339 vector <float> DeltaTime;
340 vector <int> sector;
341 vector <int> layer;
342
343// float seedTime = BcalMatrixU[idxMaxU[0]][idxMaxU[1]]->t;
344 int Idx = MaxU;
345 // collect hits for up stream shower
346 // first move upwards in sectors from the seed
347 while(EUi[Idx]) {
348 //cout<<"EUi> "<<Idx<<" "<<EUi[Idx]<<endl;
349 for (int i=0;i<4;i++){
350 const DBCALHit* hitU = BcalMatrixU[Idx][i];
351 const DBCALHit* hitD = BcalMatrixD[Idx][i];
352 if ((hitU)&&(hitD)){
353 float mt = hitD->t + hitU->t;
354 //mt -= (390./VELOCITY);
355 mt -= (390./16.75); // subract detector internal path HARD CODED VALUES!!!!!!
356 mt /= 2.;
357 MeanTime.push_back(mt);
358
359 // Note below D-U means positive dt is upstream
360 // from bcal center
361 float dt = hitD->t - hitU->t;
362 //dt = dt/2./VELOCITY;
363 dt = dt/2.*16.75; // convert to cm HARD CODED VALUE!!!!!
364 DeltaTime.push_back(dt);
365 sector.push_back(Idx);
366 layer.push_back(i);
367 }
368
369 if (hitU){
370 //if ((hitU->t - seedTime)<2.){
371 HitListU.push_back(hitU);
372 //}
373 BcalMatrixU[Idx][i] = NULL__null;
374 }
375 }
376 EUi[Idx] = 0;
377 EU[Idx] = 0.;
378 Idx++;
379 if (Idx>=4*48){
380 Idx -= 4*48;
381 }
382 }
383 // now move downwards from the seed center
384 Idx = MaxU-1;
385 if (Idx<0){
386 Idx += 4*48;
387 }
388 while(EUi[Idx]) {
389 //cout<<"EUi< "<<Idx<<" "<<EUi[Idx]<<endl;
390 for (int i=0;i<4;i++){
391 const DBCALHit* hitU = BcalMatrixU[Idx][i];
392 const DBCALHit* hitD = BcalMatrixD[Idx][i];
393 if ((hitU)&&(hitD)){
394 float mt = hitD->t + hitU->t;
395 //mt -= (390./VELOCITY);
396 mt -= (390./16.5); // subract detector internal drift
397 mt /= 2.;
398 MeanTime.push_back(mt);
399 float dt = hitD->t - hitU->t;
400 //dt = dt/2./VELOCITY;
401 dt = dt/2.*16.5; // convert to cm
402 DeltaTime.push_back(dt);
403 sector.push_back(Idx);
404 layer.push_back(i);
405 }
406
407 if (hitU){
408 //if ((hitU->t - seedTime)<2.){
409 HitListU.push_back(hitU);
410 //}
411 BcalMatrixU[Idx][i] = NULL__null;
412 }
413 }
414
415 EUi[Idx] = 0;
416 EU[Idx] = 0.;
417 Idx--;
418 if (Idx<0){
419 Idx += 4*48;
420 }
421 }
422
423// seedTime = BcalMatrixD[idxMaxD[0]][idxMaxD[1]]->t;
424 // collect hits for down stream shower
425 // first move upwards in sectors
426 Idx = MaxD;
427 while(EDi[Idx]) {
428 //cout<<"EDi> "<<Idx<<" "<<EUi[Idx]<<endl;
429 for (int i=0;i<4;i++){
430 const DBCALHit* hit = BcalMatrixD[Idx][i];
431 if (hit){
432 //if ((hit->t - seedTime)<2.){
433 HitListD.push_back(hit);
434 //}
435 BcalMatrixD[Idx][i] = NULL__null;
436 }
437 }
438 EDi[Idx] = 0;
439 ED[Idx] = 0.;
440 Idx++;
441 if (Idx>=4*48){
442 Idx -= 4*48;
443 }
444 }
445 // now move downwrads from center
446 Idx = MaxD-1;
447 if (Idx<0){
448 Idx += 4*48;
449 }
450 while(EDi[Idx]) {
451 //cout<<"EDi< "<<Idx<<" "<<EUi[Idx]<<endl;
452 for (int i=0;i<4;i++){
453 const DBCALHit* hit = BcalMatrixD[Idx][i];
454 if (hit){
455 //if ((hit->t - seedTime)<2.){
456 HitListD.push_back(hit);
457 //}
458 BcalMatrixD[Idx][i] = NULL__null;
459 }
460 }
461 EDi[Idx] = 0;
462 ED[Idx] = 0.;
463 Idx--;
464 if (Idx<0){
465 Idx += 4*48;
466 }
467 }
468
469 if ( (HitListU.size()>0) && (HitListD.size()>0) ) {
470 //cout<<"Make New Clump:"<<cnt<<endl;
471 DBCALClump* myClump = new DBCALClump(HitListU, HitListD);
472 myClump->MeanTime = MeanTime;
473 myClump->DeltaTime = DeltaTime;
474 myClump->Sector = sector;
475 myClump->Layer = layer;
476 myClump->AnalyzeClump();
477
478 int oK = 1;
479 if ((HitListU.size()==1) || (HitListD.size()==1)){
480 if (myClump->ClumpE[0]<50.) { // HARD CODED
481 oK = 0;
482 }
483 }
484 if (oK){
485 _data.push_back(myClump);
486 } else {
487 delete myClump;
488 }
489 } else {
490 cout<<"Error no hits in this Clump!!!! Event: "<<eventnumber<<" cnt= "<<cnt <<endl;
491 }
492 }
493
494 }
495
496 return NOERROR;
497}
498