Bug Summary

File:plugins/monitoring/CDC_drift/JEventProcessor_CDC_drift.cc
Location:line 400, column 5
Description:Value stored to 'previous' is never read

Annotated Source Code

1// $Id$
2//
3// File: JEventProcessor_CDC_drift.cc
4// Created: Wed Oct 22 2014
5// Creator: Naomi Jarvis
6
7
8#include <stdint.h>
9#include <vector>
10
11#include <TMath.h>
12
13
14#include "JEventProcessor_CDC_drift.h"
15#include <JANA/JApplication.h>
16#include <JANA/JFactory.h>
17
18
19using namespace std;
20using namespace jana;
21
22
23#include "CDC/DCDCHit.h"
24#include "CDC/DCDCDigiHit.h"
25#include "DAQ/Df125CDCPulse.h"
26
27#include "TRIGGER/DTrigger.h"
28
29#include <TDirectory.h>
30#include <TH2.h>
31#include <TH1.h>
32#include <TF1.h>
33#include <TTree.h>
34#include <TBranch.h>
35
36// root hist pointers
37
38static TH1F *cdc_time = NULL__null;
39static TH1I *cdc_rawtime = NULL__null;
40
41static TTree *tfit = NULL__null;
42static TTree *rtfit = NULL__null;
43
44static bool FIT_TIME = true;
45
46
47//----------------------------------------------------------------------------------
48
49
50// Routine used to create our JEventProcessor
51extern "C"{
52 void InitPlugin(JApplication *app){
53 InitJANAPlugin(app);
54 app->AddProcessor(new JEventProcessor_CDC_drift());
55 }
56}
57
58
59//----------------------------------------------------------------------------------
60
61
62JEventProcessor_CDC_drift::JEventProcessor_CDC_drift() {
63}
64
65
66//----------------------------------------------------------------------------------
67
68
69JEventProcessor_CDC_drift::~JEventProcessor_CDC_drift() {
70}
71
72
73//----------------------------------------------------------------------------------
74
75jerror_t JEventProcessor_CDC_drift::init(void) {
76
77 /*
78 // max values for histogram scales, modified fa250-format readout
79 const Int_t RTMAX = 12000; //max for raw time, less than full field width
80 const Char_t rtunits[8] = "0.125ns"; //raw time is in units of sample/64 = ns/8
81 */
82
83
84 // raw quantities for read out (125 format) are
85 // time field max 2047 scaled x 1, units 0.8ns
86 // time qf field max 1
87 // overflow count field max 7
88 // pedestal field max 255 scaled x 1/4 initially
89 // max amplitude 9 bits, field max 511 scaled x 1/8
90 // integral field max 16383 scaled x 1/14
91
92
93 // max values for histogram scales, fa125-format readout
94
95
96 const Char_t rtunits[6] = "0.8ns"; //raw time is in units of sample/10 = 0.8n
97
98 const Int_t RTMAX = 2048; //max for raw time histo, fa125-format, 11 bits
99 const Int_t TMAX = 2000; //max for time
100
101
102 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
103
104 // create root folder for cdc and cd to it, store main dir
105 TDirectory *main = gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()));
106 gDirectory(ROOT::Internal::TDirectoryAtomicAdapter(TDirectory::CurrentDirectory
()))->mkdir("CDC_drift")->cd();
107
108
109 // book histograms
110
111 cdc_time = new TH1F("cdc_time","CDC time (units of ns); time (ns)",TMAX,-500,TMAX-500);
112 cdc_rawtime = new TH1I("cdc_rawtime",Form("CDC raw time (units of %s); raw time (%s)",rtunits,rtunits),RTMAX,0,RTMAX);
113
114
115 if (FIT_TIME) {
116
117 rtfit = new TTree("rawtimefit","raw drift time fit params");
118
119 Long64_t tentries;
120 rtfit->Branch("entries",&tentries,"entries/L");
121
122 Double_t t0;
123 rtfit->Branch("t0",&t0,"t0/D");
124
125 Double_t tmax;
126 rtfit->Branch("tmax",&tmax,"tmax/D");
127
128 Double_t tmax_slope;
129 rtfit->Branch("tmax_slope",&tmax_slope,"tmax_slope/D");
130
131 Double_t tdiff;
132 rtfit->Branch("tdiff_ns",&tdiff,"tdiff/D");
133
134
135 tfit = new TTree("timefit","drift time fit params");
136
137 tfit->Branch("entries",&tentries,"entries/L");
138
139 tfit->Branch("t0",&t0,"t0/D");
140
141 tfit->Branch("tmax",&tmax,"tmax/D");
142
143 tfit->Branch("tmax_slope",&tmax_slope,"tmax_slope/D");
144
145 tfit->Branch("tdiff_ns",&tdiff,"tdiff/D");
146
147 }
148
149 japp->RootUnLock(); //RELEASE ROOT LOCK!!
150
151 main->cd();
152
153 return NOERROR;
154}
155
156
157//----------------------------------------------------------------------------------
158
159
160jerror_t JEventProcessor_CDC_drift::brun(JEventLoop *eventLoop, int32_t runnumber) {
161 // This is called whenever the run number changes
162 return NOERROR;
163}
164
165
166//----------------------------------------------------------------------------------
167
168
169jerror_t JEventProcessor_CDC_drift::evnt(JEventLoop *eventLoop, uint64_t eventnumber) {
170 // This is called for every event. Use of common resources like writing
171 // to a file or filling a histogram should be mutex protected. Using
172 // loop-Get(...) to get reconstructed objects (and thereby activating the
173 // reconstruction algorithm) should be done outside of any mutex lock
174 // since multiple threads may call this method at the same time.
175
176 // cosmics, estimate 15 mins ~ 4.4e5 events ~ 4.4e5*82/372 ~ 1e5 useful hits
177
178
179
180 const uint32_t MIN_EVENTS = 50000; //min events to collect before fitting drift time
181 const uint32_t UPDATE_INTERVAL = 25000; //incremental events required to update fit
182
183 const Bool_t RESET = kFALSE; // if true, zero histos after fitting
184 const Bool_t VERBOSE = kFALSE; // if true, print fits to stdout
185
186 uint16_t ring,straw; // ring and straw numbers from either dcdchits or dcdcdigihits
187 uint16_t n; // straw number, 1 to 3522
188 uint16_t j;
189
190 Long64_t nentries; // current number of entries
191
192 int64_t previous,tprevious; // number of entries when histo was last fitted
193
194
195 //array to make straw number n; add extra 0 at front to use offset[1] for ring 1
196 int straw_offset[29] = {0,0,42,84,138,192,258,324,404,484,577,670,776,882,1005,1128,1263,1398,1544,1690,1848,2006,2176,2346,2528,2710,2907,3104,3313};
197
198 const uint16_t nstraws = 77; //size of strawlist - list of n of straws to include in fit
199
200 const uint16_t strawlist[] = {176, 237, 496, 497, 775, 776, 777, 782, 879, 881, 882, 895, 900, 1021, 1026, 1047, 1052, 1056, 1057, 1130, 1241, 1252, 1266, 1318, 1340, 1376, 1567, 1568, 1679, 1682, 1701, 1849, 1853, 1864, 1918, 1998, 2088, 2242, 2244, 2248, 2255, 2256, 2430, 2445, 2556, 2585, 2748, 2767, 2770, 2772, 2774, 2782, 2788, 2789, 2793, 2796, 2943, 2951, 2952, 2962, 2963, 2965, 2969, 2973, 2985, 3159, 3160, 3176, 3177, 3184, 3214, 3361, 3363, 3365, 3369, 3428, 3429};
201
202
203 Bool_t fillhisto; // fill histo if true
204 Bool_t fithisto; // fit histo if true
205
206 const DTrigger* locTrigger = NULL__null;
207 eventLoop->GetSingle(locTrigger);
208 if(locTrigger->Get_L1FrontPanelTriggerBits() != 0)
209 return NOERROR;
210 if (!locTrigger->Get_IsPhysicsEvent()){ // do not look at PS triggers
211 return NOERROR;
212 }
213
214
215 // get raw data for cdc
216 vector<const DCDCDigiHit*> digihits;
217 eventLoop->Get(digihits);
218
219
220 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
221
222 fithisto = kFALSE;
223
224 previous = (uint32_t)(cdc_rawtime->GetEntries()/UPDATE_INTERVAL);
225
226
227 for (uint32_t i=0; i<digihits.size(); i++) {
228
229 const DCDCDigiHit *digihit = digihits[i];
230 const Df125CDCPulse *cp = NULL__null;
231
232 digihit->GetSingle(cp);
233
234 if (!cp) continue;
235
236 ring = digihit->ring;
237 straw = digihit->straw;
238 n = straw_offset[ring] + straw;
239
240 fillhisto = kFALSE;
241
242 if ((digihit->pulse_time) && (!cp->time_quality_bit)) {
243
244 j=0;
245
246 while ((!fillhisto) && (j<nstraws)) {
247
248 if (n == strawlist[j]) fillhisto = kTRUE;
249 j++;
250
251 }
252
253
254 if (fillhisto) {
255
256 cdc_rawtime->Fill(digihit->pulse_time);
257 nentries = cdc_rawtime->GetEntries();
258 if ((nentries > MIN_EVENTS) && (uint32_t(nentries/UPDATE_INTERVAL) > previous)) fithisto = kTRUE;
259 if (!FIT_TIME) fithisto = kFALSE;
260
261 }
262
263 }
264
265 }
266
267 Bool_t fitthisto; // fit histo if true
268
269 // get raw data for cdc
270 vector<const DCDCHit*> hits;
271 eventLoop->Get(hits);
272
273 fitthisto = kFALSE;
274
275 tprevious = (uint32_t)(cdc_time->GetEntries()/UPDATE_INTERVAL);
276
277
278 for (uint32_t i=0; i<hits.size(); i++) {
279
280 const DCDCHit *hit = hits[i];
281
282 ring = hit->ring;
283 straw = hit->straw;
284 n = straw_offset[ring] + straw;
285
286 fillhisto = kFALSE;
287
288 j=0;
289
290 while ((!fillhisto) && (j<nstraws)) {
291
292 if (n == strawlist[j]) fillhisto = kTRUE;
293 j++;
294
295 }
296
297
298 if (fillhisto) {
299
300 cdc_time->Fill(hit->t);
301 nentries = cdc_time->GetEntries();
302 if ((nentries > MIN_EVENTS) && (uint32_t(nentries/UPDATE_INTERVAL) > tprevious)) fitthisto = kTRUE;
303 if (!FIT_TIME) fitthisto = kFALSE;
304
305 }
306
307 }
308
309
310
311 if (fithisto && FIT_TIME) {
312
313 //*** the quantities to save are fitstatus, fitparams 0 to 9 and tdiff ***
314
315 if (VERBOSE) printf("\n\nFitting cdc_rawtime\n");
316
317 const float TUNITS = 0.8; // fa125 - time is in units of 0.8ns
318 // const float TUNITS = 0.125; // fa250 - time is in units of 0.125ns
319
320 Double_t fitparams[10];
321 Float_t startpar[10];
322 Int_t fitstatus;
323 Double_t tdiff=0; // max drift time in ns
324
325 Int_t imin = cdc_rawtime->FindFirstBinAbove(); // first histogram bin with counts
326 Int_t imax; // last bin with counts
327 Int_t ipeak=0; // bin with t0 peak in, used to find startparam
328
329 Double_t xmin; // x value of imin - used for fit range
330 Double_t xmax; // x value of imax
331
332 Int_t bgpeakwidth = 10; // width in bins of initial background peak, if there is one
333 Int_t bgrange = 40; //scan 4 full samples - include this many bins in background estimate
334
335 Int_t i;
336
337 Int_t chunk1 = 0;
338 for (i=0; i<bgpeakwidth; i++) chunk1 += cdc_rawtime->GetBinContent(imin+i);
339
340 Int_t chunk2 = 0;
341 for (i=0; i<bgpeakwidth; i++) chunk2 += cdc_rawtime->GetBinContent(imin+bgpeakwidth+i);
342
343 // skip over noise peak
344 if (chunk1>1.5*chunk2) imin += bgpeakwidth;
345 xmin = cdc_rawtime->GetXaxis()->GetBinLowEdge(imin);
346
347 // find max content bin
348 // search on from peak to find counts=0, end of fit range
349
350 Int_t maxcontent=0;
351 for (i=imin; i<cdc_rawtime->GetNbinsX(); i++) {
352 if (cdc_rawtime->GetBinContent(i) > maxcontent) {
353 maxcontent = cdc_rawtime->GetBinContent(i);
354 ipeak = i;
355 }
356 }
357
358 Double_t xpeak = cdc_rawtime->GetXaxis()->GetBinLowEdge(ipeak);
359
360 i=ipeak;
361 while (cdc_rawtime->GetBinContent(i) > 0 && i<cdc_rawtime->GetNbinsX() ) i++;
362
363 imax = i;
364 xmax = imax*cdc_rawtime->GetBinWidth(imax);
365
366 //starting point for background height
367 Double_t bg = 0;
368 for (i=imin; i<imin+bgrange; i++) bg += cdc_rawtime->GetBinContent(i);
369 bg = bg/(Double_t)bgrange;
370
371
372 TF1 *f = new TF1("fdrift","[9] + [0] * (1 + [1]*exp(([3]-x)/[2]) + [7]*exp(([3]-x)/[8]) ) / ( (1+exp(([3]-x)/[5])) * (1+exp((x-[4])/[6])) )",xmin,xmax);
373
374 f->SetLineWidth(1);
375 f->SetLineColor(6);
376
377 // set start values and limits here for all fit params except 0,3,4
378
379 startpar[1] = 15; //amplitude of first exp contrib to peak
380 startpar[7] = 3; //amplitude of second exp contrib to peak
381
382 f->SetParLimits(1,0,startpar[1]*2); //prev *10
383 f->SetParLimits(7,0,startpar[7]*2); //prev *10
384
385 startpar[5] = 5*0.8/TUNITS; //slope up of t0 edge
386 startpar[6] = 25*0.8/TUNITS; //slope down of tmax edge
387
388 f->SetParLimits(5,0,startpar[5]*2.5); //prev *2
389 f->SetParLimits(6,0,startpar[6]*2.5); //prev *2
390
391 startpar[2] = 20*0.8/TUNITS; //first exp fall-off
392 startpar[8] = 200*0.8/TUNITS; //second exp fall-off
393
394 f->SetParLimits(2,0,startpar[2]*3);
395 f->SetParLimits(8,startpar[2]*3,startpar[8]*3);
396
397 for (j=1;j<3;j++) f->SetParameter(j,startpar[j]);
398 for (j=5;j<9;j++) f->SetParameter(j,startpar[j]);
399
400 previous = (uint32_t)(nentries/UPDATE_INTERVAL);
Value stored to 'previous' is never read
401
402 // start values & limits for fit params 0,3,4 depend on nentries
403
404 startpar[0] = 0.0005*nentries; //overall scaling factor
405 startpar[9] = bg; //noise background
406
407 f->SetParLimits(0,0,startpar[0]*100);
408 f->SetParameter(0,startpar[0]);
409
410 f->SetParLimits(9,0,bg*2);
411 f->SetParameter(9,startpar[9]);
412
413 startpar[3] = xpeak; //t0
414 startpar[4] = xmax; //xpeak+500*0.8/TUNITS; //tmax //prev 550
415
416 f->SetParLimits(3,startpar[3]-(50*0.8/TUNITS),startpar[3]);
417 f->SetParLimits(4,startpar[3]+(500*0.8/TUNITS),xmax); //min 0.5us
418
419 f->SetParameter(3,startpar[3]);
420 f->SetParameter(4,startpar[3] + (700*0.8/TUNITS));
421
422 if (!VERBOSE) fitstatus = cdc_rawtime->Fit("f","QRLL");
423 if (VERBOSE) fitstatus = cdc_rawtime->Fit("f","RLL");
424
425 if (fitstatus == 0 || fitstatus == 2) { //fitstatus 0=good, 2=error matrix not posdef, fit params are correlated
426
427 f->GetParameters(fitparams);
428
429 tdiff = (fitparams[4] - fitparams[3])*TUNITS;
430
431 //cdc_rawtime->SetTitle(Form("Estimated max drift time is %3.2f ns",tdiff));
432
433 } else {
434
435 tdiff = 0;
436
437 }
438
439 if (VERBOSE) printf("fitstatus:%1i nentries:%5lli [0] %2.0f [1] %2.0f [2] %3.0f [3] %4.0f [4] %4.0f [5] %4.1f [6] %3.0f [7] %3.1f [8] %4.0f [9] %4.0f [tmax] %3.0f\n",fitstatus,nentries,fitparams[0],fitparams[1],fitparams[2],fitparams[3],fitparams[4],fitparams[5],fitparams[6],fitparams[7],fitparams[8],fitparams[9],tdiff);
440
441
442 rtfit->SetBranchAddress("entries",&nentries);
443 rtfit->SetBranchAddress("t0",&fitparams[3]);
444 rtfit->SetBranchAddress("tmax",&fitparams[4]);
445 rtfit->SetBranchAddress("tmax_slope",&fitparams[6]);
446 rtfit->SetBranchAddress("tdiff_ns",&tdiff);
447
448 rtfit->Fill();
449
450 delete f;
451
452 // **** reset histogram ****
453 if (RESET) cdc_rawtime->Reset();
454
455 }
456
457
458 if (fitthisto && FIT_TIME) {
459
460 //*** the quantities to save are fitstatus, fitparams 0 to 9 and tdiff ***
461
462 if (VERBOSE) printf("\n\nFitting cdc_time\n");
463
464 const float TUNITS = 1.0; // time in ns
465
466 Double_t fitparams[10];
467 Float_t startpar[10];
468 Int_t fitstatus;
469 Double_t tdiff=0; // max drift time in ns
470
471 Int_t imin = cdc_time->FindFirstBinAbove(); // first histogram bin with counts
472 Int_t imax; // last bin with counts
473 Int_t ipeak=0; // bin with t0 peak in, used to find startparam
474
475 Double_t xmin; // x value of imin - used for fit range
476 Double_t xmax; // x value of imax
477
478 Int_t bgpeakwidth = 8; // width in bins of initial background peak, if there is one
479 Int_t bgrange = 32; //scan 4 x 8ns samples - include this many bins in background estimate
480
481 Int_t i;
482
483 Double_t chunk1 = 0;
484 for (i=0; i<bgpeakwidth; i++) chunk1 += cdc_time->GetBinContent(imin+i);
485
486 Double_t chunk2 = 0;
487 for (i=0; i<bgpeakwidth; i++) chunk2 += cdc_time->GetBinContent(imin+bgpeakwidth+i);
488 if (chunk1>1.5*chunk2) imin += bgpeakwidth;
489 xmin = cdc_time->GetXaxis()->GetBinLowEdge(imin);
490
491 // find max content bin
492 // search on from peak to find counts=0, end of fit range
493
494 Int_t maxcontent=0;
495 for (i=imin; i<cdc_time->GetNbinsX(); i++) {
496 if (cdc_time->GetBinContent(i) > maxcontent) {
497 maxcontent = cdc_time->GetBinContent(i);
498 ipeak = i;
499 }
500 }
501
502 Double_t xpeak = cdc_time->GetXaxis()->GetBinLowEdge(ipeak);
503
504 i=ipeak;
505 while (cdc_time->GetBinContent(i) > 0 && i<cdc_time->GetNbinsX() ) i++;
506
507 imax = i;
508 xmax = cdc_time->GetBinLowEdge(imax);
509
510 //starting point for background height
511 Double_t bg = 0;
512 for (i=imin; i<imin+bgrange; i++) bg += cdc_time->GetBinContent(i);
513 bg = bg/(Double_t)bgrange;
514
515
516 TF1 *f = new TF1("f","[9] + [0] * (1 + [1]*exp(([3]-x)/[2]) + [7]*exp(([3]-x)/[8]) ) / ( (1+exp(([3]-x)/[5])) * (1+exp((x-[4])/[6])) )",xmin,xmax);
517
518 f->SetLineWidth(1);
519 f->SetLineColor(6);
520
521 // set start values and limits here for all fit params except 0,3,4
522
523 startpar[1] = 15; //amplitude of first exp contrib to peak
524 startpar[7] = 3; //amplitude of second exp contrib to peak
525
526 f->SetParLimits(1,0,startpar[1]*2); //prev *10
527 f->SetParLimits(7,0,startpar[7]*2); //prev *10
528
529 startpar[5] = 5*0.8/TUNITS; //slope up of t0 edge
530 startpar[6] = 25*0.8/TUNITS; //slope down of tmax edge
531
532 f->SetParLimits(5,0,startpar[5]*2.5); //prev *2
533 f->SetParLimits(6,0,startpar[6]*2.5); //prev *2
534
535 startpar[2] = 20*0.8/TUNITS; //first exp fall-off
536 startpar[8] = 200*0.8/TUNITS; //second exp fall-off
537
538 f->SetParLimits(2,0,startpar[2]*3);
539 f->SetParLimits(8,startpar[2]*3,startpar[8]*3);
540
541 for (j=1;j<3;j++) f->SetParameter(j,startpar[j]);
542 for (j=5;j<9;j++) f->SetParameter(j,startpar[j]);
543
544 previous = (uint32_t)(nentries/UPDATE_INTERVAL);
545
546 // start values & limits for fit params 0,3,4 depend on nentries
547
548 startpar[0] = 0.0005*nentries; //overall scaling factor
549 startpar[9] = bg; //noise background
550
551 f->SetParLimits(0,0,startpar[0]*100);
552 f->SetParameter(0,startpar[0]);
553
554 f->SetParLimits(9,0,bg*2);
555 f->SetParameter(9,startpar[9]);
556
557 startpar[3] = xpeak; //t0
558 startpar[4] = xmax; //xpeak + 500*0.8/TUNITS; //tmax //prev 550
559
560 f->SetParLimits(3,startpar[3]-(50*0.8/TUNITS),startpar[3]);
561 f->SetParLimits(4,startpar[3]+(500*0.8/TUNITS),xmax); // min 0.5us
562
563 f->SetParameter(3,startpar[3]);
564 f->SetParameter(4,startpar[3] + (700*0.8/TUNITS));
565
566 if (!VERBOSE) fitstatus = cdc_time->Fit("f","QRLL");
567 if (VERBOSE) fitstatus = cdc_time->Fit("f","RLL");
568
569 if (fitstatus == 0 || fitstatus == 2) { //fitstatus 0=good, 2=error matrix not posdef, fit params are correlated
570
571 f->GetParameters(fitparams);
572
573 tdiff = (fitparams[4] - fitparams[3])*TUNITS;
574
575 //cdc_time->SetTitle(Form("Estimated max drift time is %3.2f ns",tdiff));
576
577 } else {
578
579 tdiff = 0;
580
581 }
582
583 if (VERBOSE) printf("fitstatus:%1i nentries:%5lli [0] %2.0f [1] %2.0f [2] %3.0f [3] %4.0f [4] %4.0f [5] %4.1f [6] %3.0f [7] %3.1f [8] %4.0f [9] %4.0f [tmax] %3.0f\n",fitstatus,nentries,fitparams[0],fitparams[1],fitparams[2],fitparams[3],fitparams[4],fitparams[5],fitparams[6],fitparams[7],fitparams[8],fitparams[9],tdiff);
584
585
586 tfit->SetBranchAddress("entries",&nentries);
587 tfit->SetBranchAddress("t0",&fitparams[3]);
588 tfit->SetBranchAddress("tmax",&fitparams[4]);
589 tfit->SetBranchAddress("tmax_slope",&fitparams[6]);
590 tfit->SetBranchAddress("tdiff_ns",&tdiff);
591
592 tfit->Fill();
593
594 // **** reset histogram ****
595 if (RESET) cdc_time->Reset();
596
597 }
598
599
600
601
602
603 japp->RootUnLock(); //RELEASE ROOT LOCK!!
604
605 return NOERROR;
606}
607
608
609//----------------------------------------------------------------------------------
610
611
612jerror_t JEventProcessor_CDC_drift::erun(void) {
613 // This is called whenever the run number changes, before it is
614 // changed to give you a chance to clean up before processing
615 // events from the next run number.
616 return NOERROR;
617}
618
619
620//----------------------------------------------------------------------------------
621
622
623jerror_t JEventProcessor_CDC_drift::fini(void) {
624 // Called before program exit after event processing is finished.
625
626
627 return NOERROR;
628}
629
630
631//----------------------------------------------------------------------------------
632//----------------------------------------------------------------------------------