Bug Summary

File:plugins/monitoring/CDC_drift/JEventProcessor_CDC_drift.cc
Location:line 349, column 3
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
17
18using namespace std;
19using namespace jana;
20
21
22#include "CDC/DCDCHit.h"
23#include "CDC/DCDCDigiHit.h"
24#include "DAQ/Df125PulsePedestal.h"
25#include "DAQ/Df125PulseTime.h"
26
27#include <TDirectory.h>
28#include <TH2.h>
29#include <TH1.h>
30#include <TF1.h>
31#include <TTree.h>
32#include <TBranch.h>
33
34// root hist pointers
35
36static TH1I *cdc_num_events = NULL__null;
37
38static TH1I *cdc_tfit = NULL__null;
39static TH1I *cdc_afit = NULL__null;
40
41static TTree *tfit = NULL__null;
42static TTree *afit = NULL__null;
43
44static bool DISABLE_FITTING = true;
45
46//----------------------------------------------------------------------------------
47
48
49// Routine used to create our JEventProcessor
50extern "C"{
51 void InitPlugin(JApplication *app){
52 InitJANAPlugin(app);
53 app->AddProcessor(new JEventProcessor_CDC_drift());
54 }
55}
56
57
58//----------------------------------------------------------------------------------
59
60
61JEventProcessor_CDC_drift::JEventProcessor_CDC_drift() {
62}
63
64
65//----------------------------------------------------------------------------------
66
67
68JEventProcessor_CDC_drift::~JEventProcessor_CDC_drift() {
69}
70
71
72//----------------------------------------------------------------------------------
73
74jerror_t JEventProcessor_CDC_drift::init(void) {
75
76 // lock all root operations
77 //app->RootWriteLock();
78
79
80 // max values for histogram scales, modified fa250-format readout
81
82 // const Int_t RTMAX = 32768; //max for raw time, fa250-format, 15 bits
83 const Int_t RTMAX = 12000; //max for raw time, less than full field width
84
85 const Char_t rtunits[8] = "0.125ns"; //raw time is in units of sample/64 = ns/8
86
87 const float TUNITS = 0.125; // fa250 - time is in units of 0.125ns
88
89 const Int_t AMAX = 4096; //max for amplitude, fa250-format, 12 bits
90 //const Int_t AMAX = 512; //max for amplitude, fa125-format, 9 bits
91
92
93 /*
94
95 // raw quantities for read out (125 format) are
96 // time field max 2047 scaled x 1, units 0.8ns
97 // time qf field max 1
98 // overflow count field max 7
99 // pedestal field max 255 scaled x 1/4 initially
100 // max amplitude 9 bits, field max 511 scaled x 1/8
101 // integral field max 16383 scaled x 1/14
102
103
104 // max values for histogram scales, fa125-format readout
105
106 const Int_t RTMAX = 2048; //max for raw time, fa125-format, 11 bits
107 const Char_t rtunits[6] = "0.8ns"; //raw time is in units of sample/10 = 0.8ns
108
109 */
110
111
112
113
114 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
115
116 // create root folder for cdc and cd to it, store main dir
117 //TDirectory *main = gDirectory;
118 gDirectory(TDirectory::CurrentDirectory())->mkdir("CDC_fits")->cd();
119
120
121 // book histograms
122
123 cdc_num_events = new TH1I("cdc_num_events","CDC number of events",1, 0.5, 1.5);
124
125 cdc_tfit = new TH1I("cdc_tfit",Form("CDC raw time (units of %s); raw time (%s)",rtunits,rtunits),(Int_t)RTMAX*TUNITS,0,RTMAX);
126
127 cdc_afit = new TH1I("cdc_afit","CDC amplitude (ADC units); ADC units",AMAX,0,AMAX);
128
129
130 tfit = new TTree("timefit","drift time fit params");
131 afit = new TTree("ampfit","max amplitude Landau fit params");
132
133 // TTree *tfit = new TTree("tfit","drift time fit params");
134 // TTree *afit = new TTree("afit","max amplitude fit params");
135
136 Long64_t tentries;
137 tfit->Branch("entries",&tentries,"entries/L");
138
139 Double_t t0;
140 tfit->Branch("t0",&t0,"t0/D");
141
142 Double_t tmax;
143 tfit->Branch("tmax",&tmax,"tmax/D");
144
145 Double_t tmax_slope;
146 tfit->Branch("tmax_slope",&tmax_slope,"tmax_slope/D");
147
148 Double_t tdiff;
149 tfit->Branch("tdiff_ns",&tdiff,"tdiff/D");
150
151 Long64_t aentries;
152 afit->Branch("entries",&aentries,"entries/L");
153
154 Double_t normconst;
155 afit->Branch("normconst",&normconst,"normconst/D");
156
157 Double_t mpv;
158 afit->Branch("MPV",&mpv,"mpv/D");
159
160 Double_t sigma;
161 afit->Branch("sigma",&sigma,"sigma/D");
162
163
164 japp->RootUnLock(); //RELEASE ROOT LOCK!!
165
166 return NOERROR;
167}
168
169
170//----------------------------------------------------------
171
172
173jerror_t JEventProcessor_CDC_drift::brun(JEventLoop *eventLoop, int runnumber) {
174 // This is called whenever the run number changes
175 return NOERROR;
176}
177
178
179//----------------------------------------------------------------------------------
180
181
182
183jerror_t JEventProcessor_CDC_drift::evnt(JEventLoop *eventLoop, int eventnumber) {
184 // This is called for every event. Use of common resources like writing
185 // to a file or filling a histogram should be mutex protected. Using
186 // loop-Get(...) to get reconstructed objects (and thereby activating the
187 // reconstruction algorithm) should be done outside of any mutex lock
188 // since multiple threads may call this method at the same time.
189
190 //cosmics, estimate 15 mins ~ 4.4e5 events ~ 4.4e5*82/372 ~ 1e5 useful hits
191
192
193
194 const uint32_t MIN_EVENTS = 10000; //min events to collect before fitting drift time
195 const uint32_t UPDATE_INTERVAL = 10000; //incremental events required to update fit
196
197 const Bool_t RESET = kFALSE; // if true, zero histos after fitting
198 const Bool_t VERBOSE = kFALSE; // if true, print fits to stdout
199
200 const uint16_t MAXAMP_FIT_MIN = 100; // lower end of amplitude histo fit range
201 const uint16_t MAXAMP_FIT_MAX = 3800; // upper end of amplitude histo fit range
202
203 const float TUNITS = 0.125; // fa250 - time is in units of 0.125ns
204 // const float TUNITS = 0.8; // fa125 - time is in units of 0.8ns
205
206
207 uint16_t ring,straw; // ring and straw numbers from either dcdchits or dcdcdigihits
208 uint16_t n; // straw number, 1 to 3522
209 uint16_t j;
210
211 Long64_t nentries; // current number of entries
212
213 int64_t previous; // number of entries when histo was last fitted
214
215 uint32_t pulse_number=-1;
216
217 uint32_t maxamp;
218
219 //array to make straw number n; add extra 0 at front to use offset[1] for ring 1
220 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};
221
222 const uint16_t nstraws = 77; //size of strawlist - list of n of straws to include in fit
223
224 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};
225
226
227 Bool_t fillhisto; // fill histo if true
228 Bool_t fithisto; // fit histo if true
229
230
231
232 // get raw data for cdc
233 vector<const DCDCDigiHit*> digihits;
234 eventLoop->Get(digihits);
235
236
237
238
239 japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
240
241 fithisto = kFALSE;
242
243 previous = (uint32_t)(cdc_tfit->GetEntries()/UPDATE_INTERVAL);
244
245
246 for (uint32_t i=0; i<digihits.size(); i++) {
247
248 const DCDCDigiHit *digihit = digihits[i];
249 const Df125PulseTime *pt = NULL__null;
250
251
252 ring = digihit->ring;
253 straw = digihit->straw;
254 n = straw_offset[ring] + straw;
255
256 fillhisto = kFALSE;
257
258 digihit->GetSingle(pt);
259 if (pt) pulse_number = pt->pulse_number;
260
261 if ((digihit->pulse_time) && (pulse_number==0) && (!digihit->QF)) {
262
263 j=0;
264
265 while ((!fillhisto) && (j<nstraws)) {
266
267 if (n == strawlist[j]) fillhisto = kTRUE;
268 j++;
269
270 }
271
272
273 if (fillhisto) {
274
275 cdc_tfit->Fill(digihit->pulse_time);
276 nentries = cdc_tfit->GetEntries();
277
278 if ( !DISABLE_FITTING
279 && (nentries > MIN_EVENTS)
280 && (uint32_t(nentries/UPDATE_INTERVAL) > previous))
281 fithisto = kTRUE;
282
283 // find max amplitude
284 // Get pointers to the underlying objects of interest
285 const Df125PulsePedestal *pp = NULL__null;
286
287 maxamp = 0;
288
289 //get amplitude from pulse peak in pulse pedestal
290 digihit->GetSingle(pp);
291
292 if (pp) pulse_number = pp->pulse_number;
293 if ((pp) && (pulse_number==0)) maxamp = pp->pulse_peak - pp->pedestal;
294
295 if (maxamp) cdc_afit->Fill(maxamp);
296
297 }
298
299 }
300
301 }
302
303
304
305 if (fithisto) {
306
307
308 //*** the quantities to save are fitstatus, fitparams 0 to 9 and tdiff ***
309
310
311 Double_t fitparams[9];
312 Float_t startpar[9];
313 int32_t fitstatus;
314 int32_t imax; // bin with most hits in, used to find startparam
315 Double_t tdiff=0; // max drift time in ns
316
317
318
319 TF1 *f = new TF1("f","[0] * (1 + [1]*exp(([3]-x)/[2]) + [7]*exp(([3]-x)/[8]) ) / ( (1+exp(([3]-x)/[5])) * (1+exp((x-[4])/[6])) )");
320
321 f->SetLineWidth(1);
322 f->SetLineColor(6);
323
324 // set start values and limits here for all fit params except 0,3,4
325
326 startpar[1] = 15; //amplitude of first exp contrib to peak
327 startpar[7] = 3; //amplitude of second exp contrib to peak
328
329 f->SetParLimits(1,0,startpar[1]*2); //prev *10
330 f->SetParLimits(7,0,startpar[7]*2); //prev *10
331
332 startpar[5] = 5*0.8/TUNITS; //slope up of t0 edge
333 startpar[6] = 25*0.8/TUNITS; //slope down of tmax edge
334
335 f->SetParLimits(5,0,startpar[5]*2.0); //prev *10
336 f->SetParLimits(6,0,startpar[6]*2.0); //prev *10
337
338 startpar[2] = 20*0.8/TUNITS; //first exp fall-off
339 startpar[8] = 200*0.8/TUNITS; //second exp fall-off
340
341 f->SetParLimits(2,0,startpar[2]*3);
342 f->SetParLimits(8,startpar[2]*3,startpar[8]*3);
343
344 for (j=1;j<3;j++) f->SetParameter(j,startpar[j]);
345 for (j=5;j<9;j++) f->SetParameter(j,startpar[j]);
346
347
348
349 previous = (uint32_t)(nentries/UPDATE_INTERVAL);
Value stored to 'previous' is never read
350
351
352
353 // start values & limits for fit params 0,3,4 depend on nentries
354
355 startpar[0] = 0.0005*nentries; //overall scaling factor
356
357 f->SetParLimits(0,0,startpar[0]*100);
358 f->SetParameter(0,startpar[0]);
359
360 imax = cdc_tfit->GetMaximumBin();
361 imax = imax * cdc_tfit->GetBinWidth(imax);
362
363 startpar[3] = imax; //t0
364 startpar[4] = imax+500/TUNITS; //tmax //prev 550
365
366 f->SetParLimits(3,startpar[3]-(50/TUNITS),startpar[3]);
367 f->SetParLimits(4,startpar[3]+(150/TUNITS),startpar[3]+(1500/TUNITS)); //max drift 1.5us? min 0.15us???
368
369 f->SetParameter(3,startpar[3]);
370 f->SetParameter(4,startpar[4]);
371
372 if (!VERBOSE) fitstatus = cdc_tfit->Fit("f","Q");
373 if (VERBOSE) fitstatus = cdc_tfit->Fit("f");
374
375
376 if (fitstatus == 0) {
377
378 f->GetParameters(fitparams);
379
380 tdiff = (fitparams[4] - fitparams[3])*TUNITS;
381
382 //cdc_tfit->SetTitle(Form("Estimated max drift time is %3.2f ns",tdiff));
383
384 } else {
385
386 tdiff = 0;
387
388 }
389
390
391 if (VERBOSE) printf("fitstatus:%1i nentries:%5lld [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 [tmax] %3.0f\n",fitstatus,nentries,fitparams[0],fitparams[1],fitparams[2],fitparams[3],fitparams[4],fitparams[5],fitparams[6],fitparams[7],fitparams[8],tdiff);
392
393
394 tfit->SetBranchAddress("entries",&nentries);
395 tfit->SetBranchAddress("t0",&fitparams[3]);
396 tfit->SetBranchAddress("tmax",&fitparams[4]);
397 tfit->SetBranchAddress("tmax_slope",&fitparams[6]);
398 tfit->SetBranchAddress("tdiff_ns",&tdiff);
399
400 tfit->Fill();
401
402 // fit maxamp histogram
403
404 Double_t ampfitparams[3];
405
406 TF1 *lan = new TF1("lan","landau",MAXAMP_FIT_MIN,MAXAMP_FIT_MAX);
407 lan->SetLineColor(6);
408 lan->SetLineWidth(1);
409
410 if (!VERBOSE) fitstatus = cdc_afit->Fit(lan,"QRLL");
411 if (VERBOSE) fitstatus = cdc_afit->Fit(lan,"RLL");
412
413 if (!fitstatus) lan->GetParameters(ampfitparams);
414
415 nentries = cdc_afit->GetEntries();
416
417 if (VERBOSE) printf("maxamp fitstatus:%1i nentries:%5lld [0] %2.0f [1] %2.0f [2] %3.0f \n",fitstatus,nentries,ampfitparams[0],ampfitparams[1],ampfitparams[2]);
418
419 afit->SetBranchAddress("entries",&nentries);
420 afit->SetBranchAddress("normconst",&ampfitparams[0]);
421 afit->SetBranchAddress("MPV",&ampfitparams[1]);
422 afit->SetBranchAddress("sigma",&ampfitparams[2]);
423
424 afit->Fill();
425
426
427
428 // **** reset histograms ****
429 if (RESET) cdc_tfit->Reset();
430 if (RESET) cdc_afit->Reset();
431
432
433 }
434
435
436
437 japp->RootUnLock(); //RELEASE ROOT LOCK!!
438
439 //app->RootUnLock();
440
441 return NOERROR;
442}
443
444
445//----------------------------------------------------------------------------------
446
447
448jerror_t JEventProcessor_CDC_drift::erun(void) {
449 // This is called whenever the run number changes, before it is
450 // changed to give you a chance to clean up before processing
451 // events from the next run number.
452 return NOERROR;
453}
454
455
456//----------------------------------------------------------------------------------
457
458
459jerror_t JEventProcessor_CDC_drift::fini(void) {
460 // Called before program exit after event processing is finished.
461
462
463 return NOERROR;
464}
465
466
467//----------------------------------------------------------------------------------
468//----------------------------------------------------------------------------------