1 | |
2 | |
3 | |
4 | |
5 | |
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 | |
19 | using namespace std; |
20 | using 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 | |
37 | |
38 | static TH1F *cdc_time = NULL__null; |
39 | static TH1I *cdc_rawtime = NULL__null; |
40 | |
41 | static TTree *tfit = NULL__null; |
42 | static TTree *rtfit = NULL__null; |
43 | |
44 | static bool FIT_TIME = true; |
45 | |
46 | |
47 | |
48 | |
49 | |
50 | |
51 | extern "C"{ |
52 | void InitPlugin(JApplication *app){ |
53 | InitJANAPlugin(app); |
54 | app->AddProcessor(new JEventProcessor_CDC_drift()); |
55 | } |
56 | } |
57 | |
58 | |
59 | |
60 | |
61 | |
62 | JEventProcessor_CDC_drift::JEventProcessor_CDC_drift() { |
63 | } |
64 | |
65 | |
66 | |
67 | |
68 | |
69 | JEventProcessor_CDC_drift::~JEventProcessor_CDC_drift() { |
70 | } |
71 | |
72 | |
73 | |
74 | |
75 | jerror_t JEventProcessor_CDC_drift::init(void) { |
76 | |
77 | |
78 | |
79 | |
80 | |
81 | |
82 | |
83 | |
84 | |
85 | |
86 | |
87 | |
88 | |
89 | |
90 | |
91 | |
92 | |
93 | |
94 | |
95 | |
96 | const Char_t rtunits[6] = "0.8ns"; |
97 | |
98 | const Int_t RTMAX = 2048; |
99 | const Int_t TMAX = 2000; |
100 | |
101 | |
102 | japp->RootWriteLock(); |
103 | |
104 | |
105 | TDirectory *main = gDirectory(TDirectory::CurrentDirectory()); |
106 | gDirectory(TDirectory::CurrentDirectory())->mkdir("CDC_drift")->cd(); |
107 | |
108 | |
109 | |
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(); |
150 | |
151 | main->cd(); |
152 | |
153 | return NOERROR; |
154 | } |
155 | |
156 | |
157 | |
158 | |
159 | |
160 | jerror_t JEventProcessor_CDC_drift::brun(JEventLoop *eventLoop, int32_t runnumber) { |
161 | |
162 | return NOERROR; |
163 | } |
164 | |
165 | |
166 | |
167 | |
168 | |
169 | jerror_t JEventProcessor_CDC_drift::evnt(JEventLoop *eventLoop, uint64_t eventnumber) { |
170 | |
171 | |
172 | |
173 | |
174 | |
175 | |
176 | |
177 | |
178 | |
179 | |
180 | const uint32_t MIN_EVENTS = 50000; |
181 | const uint32_t UPDATE_INTERVAL = 25000; |
182 | |
183 | const Bool_t RESET = kFALSE; |
184 | const Bool_t VERBOSE = kFALSE; |
185 | |
186 | uint16_t ring,straw; |
187 | uint16_t n; |
188 | uint16_t j; |
189 | |
190 | Long64_t nentries; |
191 | |
192 | int64_t previous,tprevious; |
193 | |
194 | |
195 | |
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; |
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; |
204 | Bool_t fithisto; |
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()){ |
211 | return NOERROR; |
212 | } |
213 | |
214 | |
215 | |
216 | vector<const DCDCDigiHit*> digihits; |
217 | eventLoop->Get(digihits); |
218 | |
219 | |
220 | japp->RootWriteLock(); |
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; |
268 | |
269 | |
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 | |
314 | |
315 | if (VERBOSE) printf("\n\nFitting cdc_rawtime\n"); |
316 | |
317 | const float TUNITS = 0.8; |
318 | |
319 | |
320 | Double_t fitparams[10]; |
321 | Float_t startpar[10]; |
322 | Int_t fitstatus; |
323 | Double_t tdiff=0; |
324 | |
325 | Int_t imin = cdc_rawtime->FindFirstBinAbove(); |
326 | Int_t imax; |
327 | Int_t ipeak=0; |
328 | |
329 | Double_t xmin; |
330 | Double_t xmax; |
331 | |
332 | Int_t bgpeakwidth = 10; |
333 | Int_t bgrange = 40; |
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 | |
344 | if (chunk1>1.5*chunk2) imin += bgpeakwidth; |
345 | xmin = cdc_rawtime->GetXaxis()->GetBinLowEdge(imin); |
346 | |
347 | |
348 | |
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 | |
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 | |
378 | |
379 | startpar[1] = 15; |
380 | startpar[7] = 3; |
381 | |
382 | f->SetParLimits(1,0,startpar[1]*2); |
383 | f->SetParLimits(7,0,startpar[7]*2); |
384 | |
385 | startpar[5] = 5*0.8/TUNITS; |
386 | startpar[6] = 25*0.8/TUNITS; |
387 | |
388 | f->SetParLimits(5,0,startpar[5]*2.5); |
389 | f->SetParLimits(6,0,startpar[6]*2.5); |
390 | |
391 | startpar[2] = 20*0.8/TUNITS; |
392 | startpar[8] = 200*0.8/TUNITS; |
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 | |
403 | |
404 | startpar[0] = 0.0005*nentries; |
405 | startpar[9] = bg; |
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; |
414 | startpar[4] = xmax; |
415 | |
416 | f->SetParLimits(3,startpar[3]-(50*0.8/TUNITS),startpar[3]); |
417 | f->SetParLimits(4,startpar[3]+(500*0.8/TUNITS),xmax); |
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) { |
426 | |
427 | f->GetParameters(fitparams); |
428 | |
429 | tdiff = (fitparams[4] - fitparams[3])*TUNITS; |
430 | |
431 | |
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 | |
453 | if (RESET) cdc_rawtime->Reset(); |
454 | |
455 | } |
456 | |
457 | |
458 | if (fitthisto && FIT_TIME) { |
459 | |
460 | |
461 | |
462 | if (VERBOSE) printf("\n\nFitting cdc_time\n"); |
463 | |
464 | const float TUNITS = 1.0; |
465 | |
466 | Double_t fitparams[10]; |
467 | Float_t startpar[10]; |
468 | Int_t fitstatus; |
469 | Double_t tdiff=0; |
470 | |
471 | Int_t imin = cdc_time->FindFirstBinAbove(); |
472 | Int_t imax; |
473 | Int_t ipeak=0; |
474 | |
475 | Double_t xmin; |
476 | Double_t xmax; |
477 | |
478 | Int_t bgpeakwidth = 8; |
479 | Int_t bgrange = 32; |
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 | |
492 | |
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 | |
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 | |
522 | |
523 | startpar[1] = 15; |
524 | startpar[7] = 3; |
525 | |
526 | f->SetParLimits(1,0,startpar[1]*2); |
527 | f->SetParLimits(7,0,startpar[7]*2); |
528 | |
529 | startpar[5] = 5*0.8/TUNITS; |
530 | startpar[6] = 25*0.8/TUNITS; |
531 | |
532 | f->SetParLimits(5,0,startpar[5]*2.5); |
533 | f->SetParLimits(6,0,startpar[6]*2.5); |
534 | |
535 | startpar[2] = 20*0.8/TUNITS; |
536 | startpar[8] = 200*0.8/TUNITS; |
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 | |
547 | |
548 | startpar[0] = 0.0005*nentries; |
549 | startpar[9] = bg; |
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; |
558 | startpar[4] = xmax; |
559 | |
560 | f->SetParLimits(3,startpar[3]-(50*0.8/TUNITS),startpar[3]); |
561 | f->SetParLimits(4,startpar[3]+(500*0.8/TUNITS),xmax); |
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) { |
570 | |
571 | f->GetParameters(fitparams); |
572 | |
573 | tdiff = (fitparams[4] - fitparams[3])*TUNITS; |
574 | |
575 | |
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 | |
595 | if (RESET) cdc_time->Reset(); |
596 | |
597 | } |
598 | |
599 | |
600 | |
601 | |
602 | |
603 | japp->RootUnLock(); |
604 | |
605 | return NOERROR; |
606 | } |
607 | |
608 | |
609 | |
610 | |
611 | |
612 | jerror_t JEventProcessor_CDC_drift::erun(void) { |
613 | |
614 | |
615 | |
616 | return NOERROR; |
617 | } |
618 | |
619 | |
620 | |
621 | |
622 | |
623 | jerror_t JEventProcessor_CDC_drift::fini(void) { |
624 | |
625 | |
626 | |
627 | return NOERROR; |
628 | } |
629 | |
630 | |
631 | |
632 | |