BmnRoot
Loading...
Searching...
No Matches
BmnZDCRaw2Digit.cxx
Go to the documentation of this file.
1#include "BmnZDCRaw2Digit.h"
2
3#include "TBox.h"
4#include "TCanvas.h"
5#include "TMath.h"
6#include "TRandom.h"
7#include "TSystem.h"
8#include "TText.h"
9
10#define PRINT_ZDC_CALIBRATION 0
11
12static int nevents;
13static float amp_array[MAX_EVENTS][MAX_CHANNELS];
14static float pedestals[MAX_EVENTS][MAX_CHANNELS];
15static float profile_amp[MAX_EVENTS][MAX_CHANNELS];
16static float profile_err[MAX_EVENTS][MAX_CHANNELS];
17
18static void fcn1(Int_t& npar, Double_t* gin, Double_t& f, Double_t* par, Int_t iflag);
19
20BmnZDCRaw2Digit::BmnZDCRaw2Digit()
21{
22 n_rec = 0;
23 n_test = 0;
24 runId = 0;
25 periodId = 0;
26 for (int i = 0; i < 6; i++)
27 digiPar[i] = 0;
28}
29
30BmnZDCRaw2Digit::BmnZDCRaw2Digit(Int_t period,
31 Int_t run,
32 TString mappingFile,
33 TString CalibrationFile,
34 TString MaxPositionFile)
35{
36 periodId = period;
37 runId = run;
38 for (int i = 0; i < MAX_CHANNELS; i++)
39 zdc_amp[i] = -1.;
40 for (int i = 0; i < MAX_LOG_CHANNELS; i++)
41 log_amp[i] = -1.;
42 for (int i = 0; i < MAX_LOG_CHANNELS; i++)
43 test_chan[i] = -1;
44 for (int i = 0; i < MAX_LOG_CHANNELS; i++)
45 test_id[i] = -1;
46 for (int i = 0; i < 256; i++)
47 is_test[i] = -1;
48 n_test = 0;
49 n_rec = 0;
50 TString dummy;
51 ifstream in;
52
53 TString dir = getenv("VMCWORKDIR");
54 TString path = dir + "/input/";
55 in.open((path + mappingFile).Data());
56 if (!in.is_open()) {
57 printf("Loading ZDC Map from file: %s - file open error!\n", mappingFile.Data());
58 return;
59 }
60 printf("Loading ZDC Map from file: %s\n", mappingFile.Data());
61 in >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy;
62 maxchan = 0;
63 int ixmin = -1, ixmax = -1, iymin = -1, iymax = -1;
64 int xmin = 10000., xmax = -10000., ymin = 10000., ymax = -10000.;
65 while (!in.eof()) {
66 int id, chan, front_chan, size, ix, iy, used;
67 float x, y;
68 in >> std::hex >> id >> std::dec >> chan >> front_chan >> size >> ix >> iy >> x >> y >> used;
69 if (!in.good())
70 break;
71 // printf("%0x %d %d %d %d %d %f %f\n",id,chan,front_chan,size,ix,iy,x,y);
72 if (size >= 200 && size <= 223) {
73 if (n_test < MAX_LOG_CHANNELS && chan > 0) {
74 if (is_test[chan - 1] < 0)
75 is_test[chan - 1] = n_test;
76 else
77 is_test[128 + chan - 1] = n_test;
78 test_chan[n_test] = chan - 1;
79 test_id[n_test++] = id;
80 };
81 };
82 if (size > 2 || size < 0)
83 continue;
84 if (chan <= 0)
85 continue;
86 if (front_chan <= 0)
87 continue;
88 zdc_map_element[n_rec].id = id;
89 if (chan > 64)
90 chan -= 64;
91 zdc_map_element[n_rec].adc_chan = chan - 1;
92 if (front_chan > maxchan)
93 maxchan = front_chan;
94 zdc_map_element[n_rec].chan = front_chan - 1;
95 zdc_map_element[n_rec].size = size - 1;
96 zdc_map_element[n_rec].ix = ix;
97 zdc_map_element[n_rec].iy = iy;
98 zdc_map_element[n_rec].x = x;
99 zdc_map_element[n_rec].y = y;
100 zdc_map_element[n_rec].used = used;
101 if (x < xmin) {
102 xmin = x;
103 ixmin = n_rec;
104 }
105 if (x > xmax) {
106 xmax = x;
107 ixmax = n_rec;
108 }
109 if (y < ymin) {
110 ymin = y;
111 iymin = n_rec;
112 }
113 if (y > ymax) {
114 ymax = y;
115 iymax = n_rec;
116 }
117 n_rec++;
118 }
119 in.close();
120 cell_size[0] = 75.;
121 cell_size[1] = 150.;
122 if (ixmin >= 0) {
123 x_min = zdc_map_element[ixmin].x - cell_size[zdc_map_element[ixmin].size];
124 } else {
125 x_min = -1000.;
126 }
127 if (ixmax >= 0) {
128 x_max = zdc_map_element[ixmax].x + cell_size[zdc_map_element[ixmin].size];
129 } else {
130 x_max = +1000.;
131 }
132 if (iymin >= 0) {
133 y_min = zdc_map_element[iymin].y - cell_size[zdc_map_element[ixmin].size];
134 } else {
135 y_min = -1000.;
136 }
137 if (iymax >= 0) {
138 y_max = zdc_map_element[iymax].y + cell_size[zdc_map_element[ixmin].size];
139 } else {
140 y_max = +1000.;
141 }
142 //-------------------------------
143 FILE* fin = 0;
144 char filn[128], tit1[32] = {"Channel"}, tit2[32] = {"Calibration"}, tit3[32] = {"Error"};
145 TString path1 = dir + "/parameters/zdc/";
146
147 UniRun* Run = UniRun::GetRun(periodId, runId);
148 if (strlen(CalibrationFile.Data()) == 0) {
149 if (Run == nullptr)
150 sprintf(filn, "%szdc_muon_calibration.txt", path1.Data()); // default
151 else {
152 sprintf(filn, "%szdc_muon_calibration_run%d_%s.txt", path1.Data(), periodId, Run->GetBeamParticle().Data());
153 }
154 } else
155 sprintf(filn, "%s%s", path1.Data(), CalibrationFile.Data());
156 fin = fopen(filn, "r");
157 for (int i = 0; i < maxchan; i++) {
158 cal[i] = 1.;
159 cale[i] = 0.;
160 cal_out[i] = 1.;
161 cale_out[i] = 0.;
162 }
163 if (!fin) {
164 printf(" ZDC: Can't open calibration file %s, use default calibration coefficients 1.\n\n", filn);
165 } else {
166 printf("Calibration file %s applied.\n", filn);
167 Int_t ch1;
168 Float_t ca = 1., cae = 0.;
169
170 char tit0[200];
171 char* readptr = fgets(tit0, 200, fin);
172 int32_t read_res1 = fscanf(fin, "%d %d %d %d %d %d\n", &digiPar[0], &digiPar[1], &digiPar[2], &digiPar[3],
173 &digiPar[4], &digiPar[5]);
174 // printf("\n ZDC digi parameters %d %d %d %d %d %d\n", digiPar[0], digiPar[1], digiPar[2], digiPar[3],
175 // digiPar[4], digiPar[5]);
176
177 int32_t read_res2 = fscanf(fin, "%s %s %s\n", tit1, tit2, tit3);
178 if ((!readptr) || (read_res1 != 0) || (read_res2 != 0)) {
179 fclose(fin);
180 throw std::ios_base::failure("Calibration file read error");
181 ;
182 }
183 while (fscanf(fin, "%d %f %f\n", &ch1, &ca, &cae) == 3) {
184 if (ch1 > 0 && ch1 <= maxchan) {
185 cal[ch1 - 1] = ca;
186 cale[ch1 - 1] = cae;
187 }
188 };
189 fclose(fin);
190 }
191 if (PRINT_ZDC_CALIBRATION)
192 printf("\n ZDC calibration coefficients\n\n%s\t%s\t%s\n\n", tit1, tit2, tit3);
193 for (int i = 0; i < maxchan; i++) {
194 if (PRINT_ZDC_CALIBRATION)
195 printf("%d\t%f\t%f\n", i, cal[i], cale[i]);
196 }
197 //----------------------------------
198 path1 = dir + "/input/";
199 MaxPos_min = 0;
200 MaxPos_max = 1000;
201 if (strlen(MaxPositionFile.Data()) == 0) {
202 printf("\n ZDC: No cuts on samples wave max position.\n\n");
203 } else {
204 sprintf(filn, "%s%s", path1.Data(), MaxPositionFile.Data());
205 fin = fopen(filn, "r");
206 if (!fin) {
207 printf("\n ZDC: Can't open samples wave max position file %s !\n", filn);
208 printf(" ZDC: No cuts on samples wave max position.\n\n");
209 } else {
210 Int_t runn;
211 Float_t mpos = 10., msig = 2.;
212 while (fscanf(fin, "%d %f %f\n", &runn, &mpos, &msig) == 3) {
213 if (runn == runId) {
214 MaxPos_min = (int)(mpos - 2. * msig);
215 MaxPos_max = (int)(mpos + 2. * msig + 0.5);
216 printf("\n ZDC: Run %d Samples wave Max position range %d - %d\n\n", runn, MaxPos_min, MaxPos_max);
217 break;
218 }
219 };
220 fclose(fin);
221 }
222 }
223 //----------------------------------
224 nevents = 0;
225 use_log_function = 0;
226 thres = 40.;
227 wave2amp_flag = 1;
228 min_samples = 5;
229 ped_samples = 5;
230 use_meanxy = 0;
231 sigma_amp = 10.;
232 shower_energy = 48.;
233 shower_norm = 6.461048; // 7.5 cm modules
234 // shower_norm = 13.923893; // 15 cm modules
235 x_beam = 0.;
236 y_beam = 0.;
237
238 for (int j = 0; j < MAX_CHANNELS; j++) {
239 number[j] = 0;
240 index[j] = 0;
241 channel0[j] = 0;
242 channel1[j] = 0;
243 SampleProf[j] = NULL;
244 for (int i = 0; i < MAX_EVENTS; i++) {
245 amp_array[i][j] = 0.;
246 pedestals[i][j] = 0.;
247 profile_amp[i][j] = 0.;
248 profile_err[i][j] = 0.;
249 }
250 }
251
252 char tit[128], nam[128];
253
254 for (int i = 0; i < MAX_LOG_CHANNELS; i++)
255 htest[i] = NULL;
256 for (int i = 0; i < MAX_LOG_CHANNELS; i++)
257 TestProf[i] = NULL;
258 for (int i = 0; i < n_test; i++) {
259 htest[i] = NULL;
260 TestProf[i] = NULL;
261 if (test_chan[i] >= 0) {
262 sprintf(nam, "htest%d", i);
263 sprintf(tit, "Amplitude for adc test channel %d", test_chan[i]);
264 htest[i] = new TH1F(nam, tit, 2000, 0., 16000.);
265 htest[i]->SetDirectory(0);
266 sprintf(nam, "ptest%d", i);
267 sprintf(tit, "Average sampling wave for adc test channel %d", test_chan[i]);
268 TestProf[i] = new TProfile(nam, tit, 200, 0., 200., -100000., +100000., "s");
269 TestProf[i]->SetDirectory(0);
270 }
271 }
272 hsum_sim = new TH1F("hsumsim", "Sum of theoretical amplitudes", 1000, 0., 100.);
273 hsum_sim->SetDirectory(0);
274 hsum_raw = new TH1F("hsumraw", "Sum of raw amplitudes", 2000, 0., 20000.);
275 hsum_raw->SetDirectory(0);
276 hsum = new TH1F("hsumcal", "Sum of calibrated amplitudes", 1000, 0., 200.);
277 hsum->SetDirectory(0);
278
279 hxmean = new TH1F("hxmean", "Shower mean X", 1000, -500., +500.);
280 hxmean->SetDirectory(0);
281 hymean = new TH1F("hymean", "Shower mean Y", 1000, -500., +500.);
282 hymean->SetDirectory(0);
283
284 for (int i = 0; i < n_rec; i++) {
285 sprintf(tit, "samprof%d", zdc_map_element[i].chan + 1);
286 sprintf(nam, "Average sampling wave, module %d", zdc_map_element[i].chan + 1);
287 SampleProf[i] = new TProfile(tit, nam, 200, 0., 200., -100000., +100000., "s");
288 SampleProf[i]->SetDirectory(0);
289 }
290}
291
292BmnZDCRaw2Digit::~BmnZDCRaw2Digit()
293{
294 for (int i = 0; i < n_test; i++) {
295 if (htest[i])
296 delete htest[i];
297 htest[i] = NULL;
298 if (TestProf[i])
299 delete TestProf[i];
300 TestProf[i] = NULL;
301 }
302 if (hsum_sim)
303 delete hsum_sim;
304 if (hsum_raw)
305 delete hsum_raw;
306 if (hsum)
307 delete hsum;
308
309 if (hxmean)
310 delete hxmean;
311 if (hymean)
312 delete hymean;
313
314 for (int i = 0; i < n_rec; i++) {
315 if (SampleProf[i])
316 delete SampleProf[i];
317 }
318}
319
320void BmnZDCRaw2Digit::print()
321{
322 printf("id#\tZDC chan\tADC_chan\tsize\tix\tiy\tx\ty\tused\n");
323 for (int i = 0; i < n_rec; i++)
324 printf("0x%06lX\t%d\t%d\t%d\t%d\t%d\t%d\t%g\t%g\n", zdc_map_element[i].id, zdc_map_element[i].chan + 1,
325 zdc_map_element[i].adc_chan + 1, zdc_map_element[i].size + 1, zdc_map_element[i].ix,
326 zdc_map_element[i].iy, zdc_map_element[i].used, zdc_map_element[i].x, zdc_map_element[i].y);
327}
328
329void BmnZDCRaw2Digit::fillSampleProfiles(TClonesArray* data, Float_t x, Float_t y, Float_t e, Int_t clsize)
330{
331 Float_t amp = 0;
332 double r = 0., dx = 0., dy = 0.;
333 for (int i = 0; i < MAX_CHANNELS; i++)
334 zdc_amp[i] = -1.;
335 for (int i = 0; i < MAX_LOG_CHANNELS; i++)
336 log_amp[i] = -1.;
337 if (nevents == 0) {
338 ncells = 0;
339 x_beam = x;
340 y_beam = y;
341 shower_energy = e;
342 for (int i = 0; i < clsize * clsize; i++) {
343 channel1[i] = -1;
344 }
345 float xnear[MAX_CHANNELS], ynear[MAX_CHANNELS];
346 int ixhalf = -1, iyhalf = -1, i0 = -1, imin = -1;
347 float rmin = 10000.;
348 for (int ind = 0; ind < n_rec; ind++) {
349 number[ind] = -1;
350 channel0[ind] = -1;
351 dx = x - zdc_map_element[ind].x;
352 dy = y - zdc_map_element[ind].y;
353 r = TMath::Sqrt(dx * dx + dy * dy);
354 if (r < rmin) {
355 rmin = r;
356 imin = ind;
357 }
358 if (x >= zdc_map_element[ind].x - cell_size[zdc_map_element[ind].size] / 2.
359 && x < zdc_map_element[ind].x + cell_size[zdc_map_element[ind].size] / 2.
360 && y >= zdc_map_element[ind].y - cell_size[zdc_map_element[ind].size] / 2.
361 && y < zdc_map_element[ind].y + cell_size[zdc_map_element[ind].size] / 2)
362 {
363 if (x >= zdc_map_element[ind].x)
364 ixhalf = 1;
365 if (y >= zdc_map_element[ind].y)
366 iyhalf = 1;
367 number[ind] = ncells;
368 channel0[ind] = zdc_map_element[ind].chan;
369 channel1[ncells] = zdc_map_element[ind].chan;
370 index[ncells] = ind;
371 xnear[0] = zdc_map_element[index[ncells]].x;
372 ynear[0] = zdc_map_element[index[ncells]].y;
373 ncells++;
374 i0 = ind;
375 // printf(" ncells %d cell %d\n", ncells, zdc_map_element[ind].chan+1);
376 // break;
377 }
378 }
379 if (i0 == -1) {
380 if (x < x_min || x > x_max || y < y_min || y > y_max) {
381 printf("Beam entry point (%f,%f) outside ZDC area!\n", x, y);
382 return;
383 } else if (imin >= 0) {
384 if (x >= zdc_map_element[imin].x)
385 ixhalf = 1;
386 if (y >= zdc_map_element[imin].y)
387 iyhalf = 1;
388 number[imin] = ncells;
389 channel0[imin] = zdc_map_element[imin].chan;
390 channel1[ncells] = zdc_map_element[imin].chan;
391 index[ncells] = imin;
392 xnear[0] = zdc_map_element[index[ncells]].x;
393 ynear[0] = zdc_map_element[index[ncells]].y;
394 ncells++;
395 } else {
396 printf("Invalid beam entry point (%f,%f)!\n", x, y);
397 return;
398 }
399 }
400 for (int i = 1; i < clsize; i++) {
401 xnear[i] = zdc_map_element[index[0]].x + ixhalf * cell_size[zdc_map_element[index[0]].size];
402 ixhalf = -ixhalf;
403 }
404 for (int i = 1; i < clsize; i++) {
405 ynear[i] = zdc_map_element[index[0]].y + iyhalf * cell_size[zdc_map_element[index[0]].size];
406 iyhalf = -iyhalf;
407 }
408
409 for (int i = 0; i < clsize; i++) {
410 for (int j = 0; j < clsize; j++) {
411 for (int ind = 0; ind < n_rec; ind++) {
412 if (ind == index[0])
413 continue;
414 if (xnear[i] >= zdc_map_element[ind].x - cell_size[zdc_map_element[ind].size] / 2.
415 && xnear[i] < zdc_map_element[ind].x + cell_size[zdc_map_element[ind].size] / 2.
416 && ynear[j] >= zdc_map_element[ind].y - cell_size[zdc_map_element[ind].size] / 2.
417 && ynear[j] < zdc_map_element[ind].y + cell_size[zdc_map_element[ind].size] / 2)
418 {
419 number[ind] = ncells;
420 channel0[ind] = zdc_map_element[ind].chan;
421 channel1[ncells] = zdc_map_element[ind].chan;
422 index[ncells] = ind;
423 ncells++;
424 // printf(" ncells %d cell %d\n", ncells, zdc_map_element[ind].chan+1);
425 }
426 }
427 }
428 }
429 }
430 float ped = 0., xm = 0., ym = 0., s = 0.;
431 int num_test = 0;
432 if (data != NULL) {
433 for (int i = 0; i < data->GetEntriesFast(); i++) {
434 BmnADCDigit* digit = (BmnADCDigit*)data->At(i);
435 num_test = is_test[digit->GetChannel()];
436 if (num_test < 0)
437 num_test = is_test[digit->GetChannel() + 128];
438 if (num_test >= 0 && (digit->GetSerial() & 0xFFFFFF) == test_id[num_test]) {
439 UShort_t* samt = digit->GetUShortValue();
440 int j2 = 0;
441 for (UInt_t j1 = 0; j1 < digit->GetNSamples(); j1++) {
442 j2 = 1 - j1 % 2 + (j1 / 2) * 2;
443 if (TestProf[num_test])
444 TestProf[num_test]->Fill(j1, samt[j2]);
445 }
446 if ((amp = testwave2amp(digit->GetNSamples(), digit->GetUShortValue(), &ped)) >= 0.) {
447 if (htest[num_test])
448 htest[num_test]->Fill(amp);
449 log_amp[num_test] = amp;
450 }
451 continue;
452 }
453 int ind, num;
454 for (ind = 0; ind < n_rec; ind++)
455 if ((digit->GetSerial() & 0xFFFFFF) == zdc_map_element[ind].id
456 && digit->GetChannel() == (zdc_map_element[ind].adc_chan))
457 break;
458 if (ind == n_rec)
459 continue;
460 if (zdc_map_element[ind].used == 0)
461 continue;
462 if ((num = number[ind]) < 0)
463 continue;
464 UShort_t* sam = digit->GetUShortValue();
465 int j2 = 0;
466 for (UInt_t j = 0; j < digit->GetNSamples(); j++) {
467 j2 = 1 - j % 2 + (j / 2) * 2;
468 if (SampleProf[num])
469 SampleProf[num]->Fill(j, sam[j2] >> 4);
470 }
471
472 float signalMin = 0.;
473 float signalMax = 0.;
474 float signalPed = 0.;
475 float signalInt = 0.;
476
477 if ((amp = wave2amp(digit->GetNSamples(), digit->GetUShortValue(), &ped, &signalMin, &signalMax, &signalPed,
478 &signalInt))
479 >= 0.)
480 {
481 // printf("chan %d amp %f coef %f\n", zdc_map_element[ind].chan, amp,
482 // cal[zdc_map_element[ind].chan]);
483 xm += amp * cal[zdc_map_element[ind].chan] * zdc_map_element[ind].x;
484 ym += amp * cal[zdc_map_element[ind].chan] * zdc_map_element[ind].y;
485 s += amp * cal[zdc_map_element[ind].chan];
486 zdc_amp[zdc_map_element[ind].chan] = amp;
487 }
488 }
489 }
490 if (s > 0.)
491 hxmean->Fill(xm / s);
492 if (s > 0.)
493 hymean->Fill(ym / s);
494 nevents++;
495}
496
497void BmnZDCRaw2Digit::fillSampleProfilesAll(TClonesArray* data, Float_t x, Float_t y, Float_t e)
498{
499 Float_t amp = 0;
500 Int_t j = 0;
501 for (int i = 0; i < MAX_CHANNELS; i++)
502 zdc_amp[i] = -1.;
503 for (int i = 0; i < MAX_LOG_CHANNELS; i++)
504 log_amp[i] = -1.;
505 if (nevents >= MAX_EVENTS)
506 return;
507 if (nevents == 0) {
508 x_beam = x;
509 y_beam = y;
510 int i0 = -1;
511 for (int ind = 0; ind < n_rec; ind++) {
512 if (x >= zdc_map_element[ind].x - cell_size[zdc_map_element[ind].size] / 2.
513 && x < zdc_map_element[ind].x + cell_size[zdc_map_element[ind].size] / 2.
514 && y >= zdc_map_element[ind].y - cell_size[zdc_map_element[ind].size] / 2.
515 && y < zdc_map_element[ind].y + cell_size[zdc_map_element[ind].size] / 2)
516 {
517 i0 = ind;
518 break;
519 }
520 }
521 if (i0 == -1) {
522 printf("Warninig: Beam entry point (%f,%f) outside any ZDC module!\n", x, y);
523 // return;
524 }
525 ncells = n_rec;
526 shower_energy = e;
527 for (int i = 0; i < n_rec; i++) {
528 channel1[i] = -1;
529 }
530 for (int ind = 0; ind < n_rec; ind++) {
531 number[ind] = j;
532 channel0[ind] = zdc_map_element[ind].chan;
533 channel1[j] = zdc_map_element[ind].chan;
534 j++;
535 }
536 }
537 float ped = 0., xm = 0., ym = 0., s = 0.;
538 int num_test = 0;
539 if (data != NULL) {
540 for (int i = 0; i < data->GetEntriesFast(); i++) {
541 BmnADCDigit* digit = (BmnADCDigit*)data->At(i);
542 num_test = is_test[digit->GetChannel()];
543 if (num_test < 0)
544 num_test = is_test[digit->GetChannel() + 128];
545 if (num_test >= 0 && (digit->GetSerial() & 0xFFFFFF) == test_id[num_test]) {
546 UShort_t* samt = digit->GetUShortValue();
547 int j2 = 0;
548 for (UInt_t j1 = 0; j1 < digit->GetNSamples(); j1++) {
549 j2 = 1 - j1 % 2 + (j1 / 2) * 2;
550 if (TestProf[num_test])
551 TestProf[num_test]->Fill(j1, samt[j2]);
552 }
553 if ((amp = testwave2amp(digit->GetNSamples(), digit->GetUShortValue(), &ped)) >= 0.) {
554 if (htest[num_test])
555 htest[num_test]->Fill(amp);
556 log_amp[num_test] = amp;
557 }
558 continue;
559 }
560 int ind, num;
561 for (ind = 0; ind < n_rec; ind++)
562 if ((digit->GetSerial() & 0xFFFFFF) == zdc_map_element[ind].id
563 && digit->GetChannel() == (zdc_map_element[ind].adc_chan))
564 break;
565 if (ind == n_rec)
566 continue;
567 if (zdc_map_element[ind].used == 0)
568 continue;
569 if ((num = number[ind]) < 0)
570 continue;
571 UShort_t* sam = digit->GetUShortValue();
572 int j2 = 0;
573 for (UInt_t j1 = 0; j1 < digit->GetNSamples(); j1++) {
574 j2 = 1 - j1 % 2 + (j1 / 2) * 2;
575 if (SampleProf[num])
576 SampleProf[num]->Fill(j1, sam[j2] >> 4);
577 }
578
579 float signalMin = 0.;
580 float signalMax = 0.;
581 float signalPed = 0.;
582 float signalInt = 0.;
583
584 if ((amp = wave2amp(digit->GetNSamples(), digit->GetUShortValue(), &ped, &signalMin, &signalMax, &signalPed,
585 &signalInt))
586 >= 0.)
587 {
588 // printf("chan %d amp %f coef %f\n", zdc_map_element[ind].chan, amp,
589 // cal[zdc_map_element[ind].chan]);
590 xm += amp * cal[zdc_map_element[ind].chan] * zdc_map_element[ind].x;
591 ym += amp * cal[zdc_map_element[ind].chan] * zdc_map_element[ind].y;
592 s += amp * cal[zdc_map_element[ind].chan];
593 zdc_amp[zdc_map_element[ind].chan] = amp;
594 }
595 }
596 }
597 if (s > 0.)
598 hxmean->Fill(xm / s);
599 if (s > 0.)
600 hymean->Fill(ym / s);
601 nevents++;
602}
603
604void BmnZDCRaw2Digit::fillEvent(TClonesArray* data, TClonesArray* zdcdigit)
605{
606 for (int i = 0; i < MAX_CHANNELS; i++)
607 zdc_amp[i] = -1.;
608 for (int i = 0; i < MAX_LOG_CHANNELS; i++)
609 log_amp[i] = -1.;
610 Float_t amp = 0., ped = 0.;
611 int num_test = 0;
612 for (int i = 0; i < data->GetEntriesFast(); i++) {
613 BmnADCDigit* digit = (BmnADCDigit*)data->At(i);
614 num_test = is_test[digit->GetChannel()];
615 if (num_test < 0)
616 num_test = is_test[digit->GetChannel() + 128];
617 if (num_test >= 0 && (digit->GetSerial() & 0xFFFFFF) == test_id[num_test]) {
618 UShort_t* samt = digit->GetUShortValue();
619
620 int j2 = 0;
621 for (UInt_t j1 = 0; j1 < digit->GetNSamples(); j1++) {
622 j2 = 1 - j1 % 2 + (j1 / 2) * 2;
623 if (TestProf[num_test])
624 TestProf[num_test]->Fill(j1, samt[j2]);
625 }
626 if ((amp = testwave2amp(digit->GetNSamples(), digit->GetUShortValue(), &ped)) >= 0.) {
627 if (htest[num_test])
628 htest[num_test]->Fill(amp);
629 log_amp[num_test] = amp;
630 }
631 continue;
632 }
633 int ind;
634 for (ind = 0; ind < n_rec; ind++) {
635 if ((digit->GetSerial() & 0xFFFFFF) == zdc_map_element[ind].id
636 && digit->GetChannel() == (zdc_map_element[ind].adc_chan))
637 break;
638 }
639 if (ind == n_rec)
640 continue;
641 if (zdc_map_element[ind].used == 0)
642 continue;
643 TClonesArray& ar_zdc = *zdcdigit;
644
645 float signalMin = 0.;
646 float signalMax = 0.;
647 float signalPed = 0.;
648 float signalInt = 0.;
649
650 if ((amp = wave2amp(digit->GetNSamples(), digit->GetUShortValue(), &ped, &signalMin, &signalMax, &signalPed,
651 &signalInt))
652 >= 0.)
653 {
654 amp *= cal[zdc_map_element[ind].chan];
655 new (ar_zdc[zdcdigit->GetEntriesFast()])
656 BmnZDCDigit(zdc_map_element[ind].ix, zdc_map_element[ind].iy, zdc_map_element[ind].x,
657 zdc_map_element[ind].y, zdc_map_element[ind].size + 1, zdc_map_element[ind].chan + 1, amp,
658 signalMin, signalMax, signalPed, signalInt);
659 zdc_amp[zdc_map_element[ind].chan] = amp;
660 }
661 }
662}
663
664void BmnZDCRaw2Digit::fillAmplitudes(TClonesArray* data)
665{
666 for (int i = 0; i < MAX_CHANNELS; i++)
667 zdc_amp[i] = -1.;
668 for (int i = 0; i < MAX_LOG_CHANNELS; i++)
669 log_amp[i] = -1.;
670 Float_t amp = 0., ped = 0.;
671 int num_test = 0;
672 for (int i = 0; i < data->GetEntriesFast(); i++) {
673 BmnADCDigit* digit = (BmnADCDigit*)data->At(i);
674 num_test = is_test[digit->GetChannel()];
675 if (num_test < 0)
676 num_test = is_test[digit->GetChannel() + 128];
677 if (num_test >= 0 && (digit->GetSerial() & 0xFFFFFF) == test_id[num_test]) {
678 if ((amp = testwave2amp(digit->GetNSamples(), digit->GetUShortValue(), &ped)) >= 0.) {
679 log_amp[num_test] = amp;
680 }
681 continue;
682 }
683 int ind;
684 for (ind = 0; ind < n_rec; ind++)
685 if ((digit->GetSerial() & 0xFFFFFF) == zdc_map_element[ind].id
686 && digit->GetChannel() == (zdc_map_element[ind].adc_chan))
687 break;
688 if (ind == n_rec)
689 continue;
690 if (zdc_map_element[ind].used == 0)
691 continue;
692
693 float signalMin = 0.;
694 float signalMax = 0.;
695 float signalPed = 0.;
696 float signalInt = 0.;
697
698 if ((amp = wave2amp(digit->GetNSamples(), digit->GetUShortValue(), &ped, &signalMin, &signalMax, &signalPed,
699 &signalInt))
700 >= 0.)
701 {
702 zdc_amp[zdc_map_element[ind].chan] = amp;
703 }
704 }
705}
706
707int BmnZDCRaw2Digit::fillCalibrateCluster(TClonesArray* data, Float_t x, Float_t y, Float_t e, Int_t clsize)
708{
709 for (int i = 0; i < MAX_CHANNELS; i++)
710 zdc_amp[i] = -1.;
711 for (int i = 0; i < MAX_LOG_CHANNELS; i++)
712 log_amp[i] = -1.;
713 Float_t amp = 0;
714 double r = 0., dx = 0., dy = 0.;
715 static double coef[MAX_CHANNELS] = {1.};
716 if (nevents >= MAX_EVENTS)
717 return 1;
718 if (nevents == 0) {
719 ncells = 0;
720 x_beam = x;
721 y_beam = y;
722 shower_energy = e;
723 for (int i = 0; i < clsize * clsize; i++) {
724 channel1[i] = -1;
725 }
726 float xnear[MAX_CHANNELS], ynear[MAX_CHANNELS];
727 int ixhalf = -1, iyhalf = -1, i0 = -1, imin = -1;
728 float rmin = 10000.;
729 for (int ind = 0; ind < n_rec; ind++) {
730 number[ind] = -1;
731 channel0[ind] = -1;
732 dx = x - zdc_map_element[ind].x;
733 dy = y - zdc_map_element[ind].y;
734 r = TMath::Sqrt(dx * dx + dy * dy);
735 if (r < rmin) {
736 rmin = r;
737 imin = ind;
738 }
739 if (x >= zdc_map_element[ind].x - cell_size[zdc_map_element[ind].size] / 2.
740 && x < zdc_map_element[ind].x + cell_size[zdc_map_element[ind].size] / 2.
741 && y >= zdc_map_element[ind].y - cell_size[zdc_map_element[ind].size] / 2.
742 && y < zdc_map_element[ind].y + cell_size[zdc_map_element[ind].size] / 2)
743 {
744 if (x >= zdc_map_element[ind].x)
745 ixhalf = 1;
746 if (y >= zdc_map_element[ind].y)
747 iyhalf = 1;
748 number[ind] = ncells;
749 channel0[ind] = zdc_map_element[ind].chan;
750 channel1[ncells] = zdc_map_element[ind].chan;
751 index[ncells] = ind;
752 xnear[0] = zdc_map_element[index[ncells]].x;
753 ynear[0] = zdc_map_element[index[ncells]].y;
754 ncells++;
755 i0 = ind;
756 // printf(" ncells %d cell %d\n", ncells, zdc_map_element[ind].chan+1);
757 // break;
758 }
759 }
760 if (i0 == -1) {
761 if (x < x_min || x > x_max || y < y_min || y > y_max) {
762 printf("Beam entry point (%f,%f) outside ZDC area!\n", x, y);
763 return 1;
764 } else if (imin >= 0) {
765 if (x >= zdc_map_element[imin].x)
766 ixhalf = 1;
767 if (y >= zdc_map_element[imin].y)
768 iyhalf = 1;
769 number[imin] = ncells;
770 channel0[imin] = zdc_map_element[imin].chan;
771 channel1[ncells] = zdc_map_element[imin].chan;
772 index[ncells] = imin;
773 xnear[0] = zdc_map_element[index[ncells]].x;
774 ynear[0] = zdc_map_element[index[ncells]].y;
775 ncells++;
776 } else {
777 printf("Invalid beam entry point (%f,%f)!\n", x, y);
778 return 1;
779 }
780 }
781 for (int i = 1; i < clsize; i++) {
782 xnear[i] = zdc_map_element[index[0]].x + ixhalf * cell_size[zdc_map_element[index[0]].size];
783 ixhalf = -ixhalf;
784 }
785 for (int i = 1; i < clsize; i++) {
786 ynear[i] = zdc_map_element[index[0]].y + iyhalf * cell_size[zdc_map_element[index[0]].size];
787 iyhalf = -iyhalf;
788 }
789
790 for (int i = 0; i < clsize; i++) {
791 for (int j = 0; j < clsize; j++) {
792 for (int ind = 0; ind < n_rec; ind++) {
793 if (ind == i0)
794 continue;
795 if (xnear[i] >= zdc_map_element[ind].x - cell_size[zdc_map_element[ind].size] / 2.
796 && xnear[i] < zdc_map_element[ind].x + cell_size[zdc_map_element[ind].size] / 2.
797 && ynear[j] >= zdc_map_element[ind].y - cell_size[zdc_map_element[ind].size] / 2.
798 && ynear[j] < zdc_map_element[ind].y + cell_size[zdc_map_element[ind].size] / 2)
799 {
800 number[ind] = ncells;
801 channel0[ind] = zdc_map_element[ind].chan;
802 channel1[ncells] = zdc_map_element[ind].chan;
803 index[ncells] = ind;
804 ncells++;
805 // printf(" ncells %d cell %d\n", ncells, zdc_map_element[ind].chan+1);
806 }
807 }
808 }
809 }
810
811 for (int j = 0; j < MAX_EVENTS; j++) {
812 for (int i = 0; i < ncells; i++) {
813 amp_array[j][i] = 0.;
814 profile_amp[j][i] = 0.;
815 profile_err[j][i] = 1.;
816 }
817 }
818 // drawzdc();
819 // hsum_raw = new TH1F("hsumraw","Sum of raw amplitudes", 1000, 0., 100000.);
820 // hsum = new TH1F("hsumcal","Sum of calibrated amplitudes", 1000, 0., 100000.);
821 // hsum->Draw();
822 if (data == NULL) {
823 for (int i = 0; i < ncells; i++) {
824 coef[i] = (0.5 * ((i + 1) % 4) + 2.4) / 10.;
825 }
826 }
827 }
828 float ped = 0.;
829 float sum_raw = 0., sum_sim = 0.;
830 double el = 0.;
831 if (data != NULL) {
832 for (int i = 0; i < data->GetEntriesFast(); i++) {
833 BmnADCDigit* digit = (BmnADCDigit*)data->At(i);
834 int ind, num;
835 for (ind = 0; ind < n_rec; ind++)
836 if ((digit->GetSerial() & 0xFFFFFF) == zdc_map_element[ind].id
837 && digit->GetChannel() == (zdc_map_element[ind].adc_chan))
838 break;
839 if (ind == n_rec)
840 continue;
841 if (zdc_map_element[ind].used == 0)
842 continue;
843 if ((num = number[ind]) < 0)
844 continue;
845
846 float signalMin = 0.;
847 float signalMax = 0.;
848 float signalPed = 0.;
849 float signalInt = 0.;
850
851 if ((amp = wave2amp(digit->GetNSamples(), digit->GetUShortValue(), &ped, &signalMin, &signalMax, &signalPed,
852 &signalInt))
853 >= 0.)
854 {
855 // printf("chan %d amp %f coef %f\n", zdc_map_element[ind].chan, amp,
856 // cal[zdc_map_element[ind].chan]);
857 amp_array[nevents][num] = amp * cal[zdc_map_element[ind].chan];
858 pedestals[nevents][num] = ped;
859 sum_raw += amp_array[nevents][num];
860 }
861 }
862 for (int i = 0; i < ncells; i++) {
863 dx = x - zdc_map_element[index[i]].x;
864 dy = y - zdc_map_element[index[i]].y;
865 r = TMath::Sqrt(dx * dx + dy * dy);
866 el = shower(r, cell_size[zdc_map_element[index[i]].size]);
867 sum_sim += el;
868 }
869 } else {
870 for (int i = 0; i < ncells; i++) {
871 dx = x - zdc_map_element[index[i]].x;
872 dy = y - zdc_map_element[index[i]].y;
873 r = TMath::Sqrt(dx * dx + dy * dy);
874 el = shower(r, cell_size[zdc_map_element[index[i]].size]);
875 sum_sim += el;
876 el /= coef[i];
877 amp = el + gRandom->Gaus() * sigma_amp;
878 amp_array[nevents][i] = amp * cal[zdc_map_element[index[i]].chan];
879 pedestals[nevents][i] = ped;
880 sum_raw += amp_array[nevents][i];
881 }
882 }
883 hsum_raw->Fill(sum_raw);
884 hsum_sim->Fill(sum_sim);
885 cellWeight(nevents);
886 nevents++;
887 return 0;
888}
889
890int BmnZDCRaw2Digit::fillCalibrateNumbers(TClonesArray* data,
891 Float_t x,
892 Float_t y,
893 Float_t e,
894 Int_t nchan,
895 Int_t* cells)
896{
897 for (int i = 0; i < MAX_CHANNELS; i++)
898 zdc_amp[i] = -1.;
899 for (int i = 0; i < MAX_LOG_CHANNELS; i++)
900 log_amp[i] = -1.;
901 Float_t amp = 0;
902 double r = 0., dx = 0., dy = 0.;
903 if (nevents >= MAX_EVENTS)
904 return 1;
905 if (nevents == 0) {
906 ncells = nchan;
907 x_beam = x;
908 y_beam = y;
909 int i0 = -1;
910 for (int ind = 0; ind < n_rec; ind++) {
911 if (x >= zdc_map_element[ind].x - cell_size[zdc_map_element[ind].size] / 2.
912 && x < zdc_map_element[ind].x + cell_size[zdc_map_element[ind].size] / 2.
913 && y >= zdc_map_element[ind].y - cell_size[zdc_map_element[ind].size] / 2.
914 && y < zdc_map_element[ind].y + cell_size[zdc_map_element[ind].size] / 2)
915 {
916 i0 = ind;
917 break;
918 }
919 }
920 if (i0 == -1) {
921 if (x < x_min || x > x_max || y < y_min || y > y_max) {
922 printf("Warning: Beam entry point (%f,%f) outside ZDC area!\n", x, y);
923 // return 1;
924 }
925 }
926 shower_energy = e;
927 for (int i = 0; i < ncells; i++) {
928 channel1[i] = -1;
929 }
930 for (int ind = 0; ind < n_rec; ind++) {
931 number[ind] = -1;
932 channel0[ind] = -1;
933 for (int j = 0; j < ncells; j++) {
934 if ((cells[j] - 1) == zdc_map_element[ind].chan) {
935 number[ind] = j;
936 channel0[ind] = zdc_map_element[ind].chan;
937 channel1[j] = zdc_map_element[ind].chan;
938 index[j] = ind;
939 break;
940 }
941 }
942 }
943 for (int j = 0; j < MAX_EVENTS; j++) {
944 for (int i = 0; i < ncells; i++) {
945 amp_array[j][i] = 0.;
946 profile_amp[j][i] = 0.;
947 profile_err[j][i] = 1.;
948 }
949 }
950 }
951 float ped = 0.;
952 float sum_raw = 0., sum_sim = 0.;
953 double el = 0.;
954 for (int i = 0; i < data->GetEntriesFast(); i++) {
955 BmnADCDigit* digit = (BmnADCDigit*)data->At(i);
956 int ind, num;
957 for (ind = 0; ind < n_rec; ind++)
958 if ((digit->GetSerial() & 0xFFFFFF) == zdc_map_element[ind].id
959 && digit->GetChannel() == (zdc_map_element[ind].adc_chan))
960 break;
961 if (ind == n_rec)
962 continue;
963 if (zdc_map_element[ind].used == 0)
964 continue;
965 if ((num = number[ind]) < 0)
966 continue;
967
968 float signalMin = 0.;
969 float signalMax = 0.;
970 float signalPed = 0.;
971 float signalInt = 0.;
972
973 if ((amp = wave2amp(digit->GetNSamples(), digit->GetUShortValue(), &ped, &signalMin, &signalMax, &signalPed,
974 &signalInt))
975 >= 0.)
976 {
977 amp_array[nevents][num] = amp * cal[zdc_map_element[ind].chan];
978 pedestals[nevents][num] = ped;
979 sum_raw += amp_array[nevents][num];
980 }
981 }
982 for (int i = 0; i < ncells; i++) {
983 dx = x - zdc_map_element[index[i]].x;
984 dy = y - zdc_map_element[index[i]].y;
985 r = TMath::Sqrt(dx * dx + dy * dy);
986 el = shower(r, cell_size[zdc_map_element[index[i]].size]);
987 sum_sim += el;
988 }
989 hsum_raw->Fill(sum_raw);
990 hsum_sim->Fill(sum_sim);
991 cellWeight(nevents);
992 nevents++;
993 return 0;
994}
995
996int BmnZDCRaw2Digit::fillCalibrateAll(TClonesArray* data, Float_t x, Float_t y, Float_t e)
997{
998 for (int i = 0; i < MAX_CHANNELS; i++)
999 zdc_amp[i] = -1.;
1000 for (int i = 0; i < MAX_LOG_CHANNELS; i++)
1001 log_amp[i] = -1.;
1002 Float_t amp = 0;
1003 Int_t j = 0;
1004 if (nevents >= MAX_EVENTS)
1005 return 1;
1006 if (nevents == 0) {
1007 x_beam = x;
1008 y_beam = y;
1009 int i0 = -1;
1010 for (int ind = 0; ind < n_rec; ind++) {
1011 if (x >= zdc_map_element[ind].x - cell_size[zdc_map_element[ind].size] / 2.
1012 && x < zdc_map_element[ind].x + cell_size[zdc_map_element[ind].size] / 2.
1013 && y >= zdc_map_element[ind].y - cell_size[zdc_map_element[ind].size] / 2.
1014 && y < zdc_map_element[ind].y + cell_size[zdc_map_element[ind].size] / 2)
1015 {
1016 i0 = ind;
1017 break;
1018 }
1019 }
1020 if (i0 == -1) {
1021 printf("Warning: Beam entry point (%f,%f) outside any ZDC module!\n", x, y);
1022 // return 1;
1023 }
1024 ncells = n_rec;
1025 shower_energy = e;
1026 for (int i = 0; i < n_rec; i++) {
1027 channel1[i] = -1;
1028 }
1029 for (int ind = 0; ind < n_rec; ind++) {
1030 number[ind] = j;
1031 channel0[ind] = zdc_map_element[ind].chan;
1032 channel1[j] = zdc_map_element[ind].chan;
1033 j++;
1034 }
1035 for (int j1 = 0; j1 < MAX_EVENTS; j1++) {
1036 for (int i = 0; i < n_rec; i++) {
1037 amp_array[j1][i] = 0.;
1038 profile_amp[j1][i] = 0.;
1039 profile_err[j1][i] = 1.;
1040 }
1041 }
1042 }
1043 float ped = 0.;
1044 for (int i = 0; i < data->GetEntriesFast(); i++) {
1045 BmnADCDigit* digit = (BmnADCDigit*)data->At(i);
1046 int ind, num;
1047 for (ind = 0; ind < n_rec; ind++)
1048 if ((digit->GetSerial() & 0xFFFFFF) == zdc_map_element[ind].id
1049 && digit->GetChannel() == (zdc_map_element[ind].adc_chan))
1050 break;
1051 if (ind == n_rec)
1052 continue;
1053 if (zdc_map_element[ind].used == 0)
1054 continue;
1055 if ((num = number[ind]) < 0)
1056 continue;
1057
1058 float signalMin = 0.;
1059 float signalMax = 0.;
1060 float signalPed = 0.;
1061 float signalInt = 0.;
1062
1063 if ((amp = wave2amp(digit->GetNSamples(), digit->GetUShortValue(), &ped, &signalMin, &signalMax, &signalPed,
1064 &signalInt))
1065 >= 0.)
1066 {
1067 amp_array[nevents][num] = amp * cal[zdc_map_element[ind].chan];
1068 pedestals[nevents][num] = ped;
1069 }
1070 }
1071 cellWeight(nevents);
1072 nevents++;
1073 return 0;
1074}
1075
1076void BmnZDCRaw2Digit::cellWeight(int ievent)
1077{
1078 // int n;
1079 double s = 0., e = 0., r = 0., x0 = 0., y0 = 0.;
1080 if (use_meanxy) {
1081 for (int i = 0; i < ncells; i++) {
1082 x0 += zdc_map_element[index[i]].x * amp_array[ievent][i];
1083 y0 += zdc_map_element[index[i]].y * amp_array[ievent][i];
1084 s += amp_array[ievent][i];
1085 }
1086 x0 /= s;
1087 y0 /= s;
1088 } else {
1089 x0 = x_beam;
1090 y0 = y_beam;
1091 }
1092 for (int i = 0; i < ncells; i++) {
1093 r = TMath::Sqrt((zdc_map_element[index[i]].x - x0) * (zdc_map_element[index[i]].x - x0)
1094 + (zdc_map_element[index[i]].y - y0) * (zdc_map_element[index[i]].y - y0));
1095 e = shower(r, cell_size[zdc_map_element[index[i]].size]);
1096 profile_amp[ievent][i] = e;
1097 profile_err[ievent][i] = sigma_amp * sigma_amp;
1098 }
1099}
1100
1101void BmnZDCRaw2Digit::calibrate()
1102{
1103 // initialize TMinuit with a maximum of 5 params
1104
1105 Double_t vstart[MAX_CHANNELS], step[MAX_CHANNELS], vmin = 0.00001, vmax = 100000.;
1106
1107 // printf(" ncells = %d\n", ncells);
1108 TMinuit* gMinuit1 = new TMinuit(ncells + 1);
1109 gMinuit1->SetFCN(fcn1);
1110
1111 Double_t arglist[MAX_CHANNELS];
1112 Int_t ierflg = 0;
1113
1114 arglist[0] = 1;
1115 gMinuit1->mnexcm("SET ERR", arglist, 1, ierflg);
1116
1117 // Set starting values and step sizes for parameters
1118 char name[32] = {""};
1119 // printf(" ncells1 = %d\n", ncells);
1120 for (int i = 0; i < ncells; i++) {
1121 vstart[i] = pstart;
1122 step[i] = pstep;
1123 sprintf(name, "Coeff_%03d", channel1[i] + 1);
1124 gMinuit1->mnparm(i, name, vstart[i], step[i], vmin, vmax, ierflg);
1125 // printf("%d name %s min %f max %f start %f step %f\n",i,name,vmin,vmax,vstart[i],step[i]);
1126 }
1127 // Now ready for minimization step
1128 arglist[0] = 50000;
1129 arglist[1] = 0.1;
1130 gMinuit1->mnexcm("MIGRAD", arglist, 2, ierflg);
1131
1132 // Print results
1133 Double_t amin, edm, errdef;
1134 Int_t nvpar, nparx, icstat;
1135 gMinuit1->mnstat(amin, edm, errdef, nvpar, nparx, icstat);
1136 gMinuit1->mnprin(3, amin);
1137
1138 //----- write results
1139
1140 for (int i = 0; i < maxchan; i++) {
1141 cal_out[i] = cal[i];
1142 cale_out[i] = cale[i];
1143 }
1144 Double_t par[MAX_CHANNELS], epar[MAX_CHANNELS];
1145 for (int i = 0; i < ncells; i++) {
1146 gMinuit1->GetParameter(i, par[i], epar[i]);
1147 cal_out[channel1[i]] = cal[channel1[i]] * par[i];
1148 cale_out[channel1[i]] = cal[channel1[i]] * epar[i];
1149 // printf(" %d par %f epar %f\n",i,par[i],epar[i]);
1150 }
1151 TString dir = getenv("VMCWORKDIR");
1152 TString path = dir + "/parameters/zdc/";
1153 FILE* fout = 0;
1154 char filn[128], tit1[32] = {"Channel"}, tit2[32] = {"Calibration"}, tit3[32] = {"Error"};
1155 sprintf(filn, "%szdc_period%d_run%d_calibration_out.txt", path.Data(), periodId, runId);
1156 fout = fopen(filn, "w");
1157 if (!fout) {
1158 printf("Can't open output calibration file %s\n", filn);
1159 } else {
1160 fprintf(fout, "%s\t%s\t%s\n", tit1, tit2, tit3);
1161 for (int i = 0; i < maxchan; i++) {
1162 fprintf(fout, "%3d\t%f\t%f\n", i + 1, cal_out[i], cale_out[i]);
1163 }
1164 fclose(fout);
1165 }
1166 printf("%s\t%s\t%s\n", tit1, tit2, tit3);
1167 for (int i = 0; i < maxchan; i++) {
1168 printf("%3d\t%f\t%f\n", i + 1, cal_out[i], cale_out[i]);
1169 }
1170 //-------
1171 float sum = 0.;
1172 for (int j = 0; j < nevents; j++) {
1173 sum = 0.;
1174 for (int i = 0; i < ncells; i++) {
1175 // sum += amp_array[j][i]*cal_out[channel1[i]];
1176 sum += amp_array[j][i] * par[i];
1177 }
1178 hsum->Fill(sum);
1179 }
1180 // hsum->Draw();
1181 // drawzdc();
1182}
1183
1184static void fcn1(Int_t& npar, Double_t* gin, Double_t& f, Double_t* par, Int_t iflag)
1185{
1186 double chi = 0., le = 0., laA = 0.;
1187 for (int j = 0; j < nevents; j++) {
1188 for (int i = 0; i < npar; i++) {
1189 // if (j == nevents-1) printf(" nevents %d par %f prof %f eprof %f amp %f\n", nevents, par[i],
1190 // profile_amp[j][i], profile_err[j][i], amp_array[j][i]);
1191 if (use_log_function) {
1192 le = TMath::Log10(profile_amp[j][i]);
1193 laA = TMath::Log10(par[i] * amp_array[j][i]);
1194 } else {
1195 le = profile_amp[j][i];
1196 laA = par[i] * amp_array[j][i];
1197 }
1198 chi += (le - laA) * (le - laA) / profile_err[j][i];
1199 }
1200 }
1201 f = chi;
1202 return;
1203}
1204
1205double BmnZDCRaw2Digit::ch(double x, double r)
1206{
1207 double amp = (exp(x / r) + exp(-x / r)) / 2.;
1208 return amp;
1209}
1210double BmnZDCRaw2Digit::sh(double x, double r)
1211{
1212 double amp = (exp(x / r) - exp(-x / r)) / 2.;
1213 return amp;
1214}
1215
1216//-----------------------------------------------PP1
1217double BmnZDCRaw2Digit::PP1(double x, double h)
1218{ // Shower profile in cell with size h when R <= h/2
1219
1220 double a1, a2, b1, b2, B;
1221 double ss1 = 1.;
1222 double bb2 = 6. * ss1;
1223 double bb1 = 1.4 * ss1;
1224
1225 double E = shower_energy; // total hadron energy , GeV
1226
1227 a2 = 0.105 - 0.014 * log(E); // R -dependences
1228 b1 = bb1 - 0.120 * log(E); // R -dependences
1229 b2 = bb2 - 0.260 * log(E); // R -dependences
1230
1231 a1 = 1. - a2;
1232 B = (a1 * b1 + a2 * b2);
1233 double fff = 0.;
1234 if (x <= h / 2)
1235 fff = (a1 * b1 * (1. - exp(-h / (2. * b1)) * ch(x, b1)) + a2 * b2 * (1. - exp(-h / (2. * b2)) * ch(x, b2)))
1236 / (h * B);
1237 return fff;
1238}
1239
1240//---------------------------------------- PP2
1241double BmnZDCRaw2Digit::PP2(double x, double h)
1242{ // Shower profile in cell with size h when R > h/2
1243
1244 double a1, a2, b1, b2, B;
1245 double ss1 = 1.;
1246 double bb2 = 6. * ss1;
1247 double bb1 = 1.4 * ss1;
1248
1249 double E = shower_energy; // total hadron energy , GeV
1250
1251 a2 = 0.105 - 0.014 * log(E); // R -dependences
1252 b1 = bb1 - 0.120 * log(E); // R -dependences
1253 b2 = bb2 - 0.260 * log(E); // R -dependences
1254
1255 a1 = 1. - a2;
1256 B = (a1 * b1 + a2 * b2);
1257 double fff = 0.;
1258 if (x > h / 2)
1259 fff = (a1 * b1 * exp(-x / b1) * sh(h, 2. * b1)) / (h * B) + (a2 * b2 * exp(-x / b2) * sh(h, 2. * b2)) / (h * B);
1260 return fff;
1261}
1262
1263double BmnZDCRaw2Digit::shower(double x, double h)
1264{
1265 double amp = PP1(x / 10., h / 10.) + PP2(x / 10., h / 10.);
1266 return shower_energy * shower_norm * amp;
1267}
1268
1269float BmnZDCRaw2Digit::wave2amp(UChar_t ns,
1270 UShort_t* s,
1271 Float_t* pedestal,
1272 Float_t* sigMin,
1273 Float_t* sigMax,
1274 Float_t* sigPed,
1275 Float_t* sigInt)
1276{
1277 float pedest = 0., ampl = 0., /*signal = 0., */ signal_max = 0.;
1278 // int nsignal = 0, nsignal_max = 0, ismax = -1, ismin = -1, m1 = 0;
1279
1280 float ampl_max = -32768.;
1281 float ampl_min = 32768.;
1282
1283 int pedBegin = digiPar[0];
1284 int pedEnd = digiPar[1];
1285 int gateBegin = digiPar[2];
1286 int gateEnd = digiPar[3];
1287
1288 int signalThresh = digiPar[4];
1289 int signalType = digiPar[5]; // 0 - max, 1 - integral
1290
1291 float scaleSignal = -1.0 / 16.; // invert + shift >>4
1292
1293 float integral = 0.;
1294
1295 if (ns > 0) {
1296 for (int m = 0; m < ns; m++) {
1297
1298 float value = scaleSignal * (Short_t)s[m];
1299
1300 if (m >= pedBegin && m <= pedEnd) {
1301 pedest += value;
1302 if (m == pedEnd)
1303 pedest /= (pedEnd - pedBegin + 1);
1304 continue;
1305 } else if (m >= gateBegin && m <= gateEnd) {
1306 ampl = value - pedest;
1307 integral += ampl;
1308 if (ampl > ampl_max) {
1309 ampl_max = ampl;
1310 // ismax = m;
1311 }
1312 if (ampl < ampl_min) {
1313 ampl_min = ampl;
1314 // ismin = m;
1315 }
1316 }
1317 } // loop over samples
1318
1319 // printf("Amplmax %f Ped %f imax %d\n", ampl_max, pedest, ismax);
1320
1321 if (signalType == 0)
1322 signal_max = ampl_max;
1323 else
1324 signal_max = integral;
1325
1326 // store pure signal parameters
1327 *sigMin = ampl_min;
1328 *sigMax = ampl_max;
1329 *sigPed = pedest;
1330 *sigInt = integral;
1331
1332 if (fabs(ampl_max - ampl_min) < signalThresh)
1333 signal_max = 0.;
1334
1335 } // if samples exist
1336 else
1337 {
1338 signal_max = -1.;
1339 }
1340
1341 *pedestal = pedest;
1342 return signal_max;
1343
1344 // old code
1345 /*
1346 if (ns > 0)
1347 {
1348 for (int m = 0; m < ns; m++)
1349 {
1350 m1 = 1 - (m%2) + (m/2)*2;
1351 if (m < ped_samples)
1352 {
1353 pedest += ((Short_t)s[m1]>>4);
1354 if (m == (ped_samples-1))
1355 {
1356 pedest /= ped_samples;
1357 }
1358 continue;
1359 }
1360 else
1361 {
1362 ampl = -(float)((Short_t)s[m1]>>4) + pedest;
1363 if (ampl > thres)
1364 {
1365 signal += ampl;
1366 nsignal++;
1367 if (ampl > ampl_max)
1368 {
1369 ampl_max = ampl;
1370 ismax = m;
1371 }
1372 }
1373 else
1374 {
1375 if (nsignal < min_samples)
1376 {
1377 signal = 0.;
1378 nsignal = 0;
1379 }
1380 else if (nsignal > nsignal_max)
1381 {
1382 if (MaxPos_min <= ismax && MaxPos_max >= ismax)
1383 {
1384 signal_max = signal;
1385 nsignal_max = nsignal;
1386 }
1387 else
1388 {
1389 signal = 0.;
1390 nsignal = 0;
1391 nsignal_max = 0;
1392 }
1393 }
1394 }
1395 }
1396 } // loop over samples
1397 // printf("Chan %d Amplmax %f Ped %f imax %d nsam %d\n",l,ampl_max,pedest[l],ismax,nsamples[l]);
1398 if (nsignal_max > 0 || ampl_max > 0.)
1399 {
1400 if (wave2amp_flag)
1401 {
1402 if (MaxPos_min <= ismax && MaxPos_max >= ismax) signal_max = ampl_max;
1403 else signal_max = 0.;
1404 }
1405 }
1406 else
1407 {
1408 signal_max = 0.;
1409 }
1410 } // if samples exist
1411 else
1412 {
1413 signal_max = -1.;
1414 }
1415
1416 *pedestal = pedest;
1417 return signal_max;
1418
1419 */
1420}
1421
1422float BmnZDCRaw2Digit::testwave2amp(UChar_t ns, UShort_t* s, Float_t* pedestal)
1423{
1424 float pedest = 0., ampl = 0., signal = 0., signal_max = 0.;
1425 int nsignal = 0, nsignal_max = 0, /*ismax = -1, */ m1 = 0;
1426 float ampl_max = 0.;
1427 if (ns > 0) {
1428 for (int m = 0; m < ns; m++) {
1429 m1 = 1 - (m % 2) + (m / 2) * 2;
1430 if (m < ped_samples) {
1431 pedest += (Short_t)s[m1] >> 4;
1432 if (m == (ped_samples - 1)) {
1433 pedest /= ped_samples;
1434 }
1435 continue;
1436 } else {
1437 ampl = -(float)((Short_t)s[m1] >> 4) + pedest;
1438 if (ampl > 0) {
1439 signal += ampl;
1440 nsignal++;
1441 if (ampl > ampl_max) {
1442 ampl_max = ampl;
1443 // ismax = m;
1444 }
1445 } else {
1446 if (nsignal < min_samples) {
1447 signal = 0.;
1448 nsignal = 0;
1449 } else if (nsignal > nsignal_max) {
1450 signal_max = signal;
1451 nsignal_max = nsignal;
1452 }
1453 }
1454 }
1455 } // loop over samples
1456 // printf("Chan %d Amplmax %f Ped %f imax %d nsam %d\n",l,ampl_max,pedest[l],ismax,nsamples[l]);
1457 if (nsignal_max > 0 || ampl_max > 0.) {
1458 if (wave2amp_flag)
1459 signal_max = ampl_max;
1460 } else {
1461 signal_max = 0.;
1462 }
1463 } // if samples exist
1464 else
1465 {
1466 signal_max = -1.;
1467 }
1468 *pedestal = pedest;
1469 return signal_max;
1470}
1471
1472void BmnZDCRaw2Digit::drawzdc(int nohist)
1473{
1474 TCanvas* callbe0 = new TCanvas("callbe0", "ZDC mapping", 800, 800);
1475 gPad->Range(x_min - 10., y_min - 10., x_max + 10., y_max + 10.);
1476 int i;
1477 char text[16];
1478 double x1, x2, y1, y2;
1479 TBox* b[MAX_CHANNELS] = {0};
1480 TText* t[MAX_CHANNELS] = {0};
1481 TText* tc[MAX_CHANNELS] = {0};
1482 callbe0->cd();
1483 for (i = 0; i < n_rec; i++) {
1484 x1 = zdc_map_element[i].x - cell_size[zdc_map_element[i].size] / 2.;
1485 x2 = zdc_map_element[i].x + cell_size[zdc_map_element[i].size] / 2.;
1486 y1 = zdc_map_element[i].y - cell_size[zdc_map_element[i].size] / 2.;
1487 y2 = zdc_map_element[i].y + cell_size[zdc_map_element[i].size] / 2.;
1488 // printf("%f %f %f %f\n",x1,y1,x2,y2);
1489 b[i] = new TBox(x1, y1, x2, y2);
1490 b[i]->Draw();
1491 if (ncells > 0 && number[i] >= 0)
1492 b[i]->SetLineColor(kRed);
1493 b[i]->SetFillStyle(0);
1494 b[i]->SetLineStyle(0);
1495 sprintf(text, "%d", zdc_map_element[i].chan + 1);
1496 t[i] = new TText(zdc_map_element[i].x, zdc_map_element[i].y, text);
1497 t[i]->Draw();
1498 if (ncells > 0 && number[i] >= 0) {
1499 t[i]->SetTextColor(kRed);
1500 t[i]->SetTextAlign(21);
1501 } else {
1502 t[i]->SetTextAlign(21);
1503 }
1504 t[i]->SetTextSize(0.02);
1505 // if (ncells > 0 && number[i] >= 0)
1506 // {
1507 // sprintf(text,"%6.3f",cal_out[channel0[i]]);
1508 sprintf(text, "%6.3f", cal_out[zdc_map_element[i].chan]);
1509 tc[i] = new TText(zdc_map_element[i].x, zdc_map_element[i].y, text);
1510 tc[i]->Draw();
1511 if (ncells > 0 && number[i] >= 0)
1512 tc[i]->SetTextColor(kRed);
1513 tc[i]->SetTextAlign(23);
1514 tc[i]->SetTextSize(0.02);
1515 // }
1516 }
1517
1518 if ((hsum_raw == NULL && hsum == NULL) || nohist)
1519 return;
1520
1521 TCanvas* calres = new TCanvas("calres", "ZDC calibration", 800, 800);
1522 calres->cd();
1523 calres->Divide(1, 2);
1524 calres->cd(1);
1525 if (hsum_raw)
1526 hsum_raw->Draw();
1527 gPad->AddExec("exselt", "select_hist()");
1528 calres->cd(2);
1529 if (hsum)
1530 hsum->Draw();
1531 gPad->AddExec("exselt", "select_hist()");
1532 // calres->cd(3);
1533 // if (hsum_sim) hsum_sim->Draw();
1534 // gPad->AddExec("exselt","select_hist()");
1535 return;
1536}
1537
1538void BmnZDCRaw2Digit::drawprof()
1539{
1540 TCanvas* csampro = new TCanvas("csampro", "ZDC sample profiles", 800, 800);
1541 csampro->cd();
1542 if (ncells <= 4)
1543 csampro->Divide(2, 2);
1544 else if (ncells <= 9)
1545 csampro->Divide(3, 3);
1546 else if (ncells <= 16)
1547 csampro->Divide(4, 4);
1548 else if (ncells <= 25)
1549 csampro->Divide(5, 5);
1550 else if (ncells <= 36)
1551 csampro->Divide(6, 6);
1552 else if (ncells <= 49)
1553 csampro->Divide(7, 7);
1554 else if (ncells <= 49)
1555 csampro->Divide(7, 7);
1556 else
1557 csampro->Divide(8, 7);
1558 for (int i = 0; i < ncells; i++) {
1559 csampro->cd(i + 1);
1560 if (SampleProf[i])
1561 SampleProf[i]->Draw();
1562 gPad->AddExec("exselt", "select_hist()");
1563 }
1564
1565 TCanvas* calres = new TCanvas("calres", "ZDC mean X and Y", 800, 800);
1566 calres->cd();
1567 calres->Divide(1, 2);
1568 calres->cd(1);
1569 if (hxmean)
1570 hxmean->Draw();
1571 gPad->AddExec("exselt", "select_hist()");
1572 calres->cd(2);
1573 if (hymean)
1574 hymean->Draw();
1575 gPad->AddExec("exselt", "select_hist()");
1576 // hsum_raw->Print();
1577 return;
1578}
1579
1580void BmnZDCRaw2Digit::drawtest()
1581{
1582 if (htest[0] == NULL)
1583 return;
1584 TCanvas* ctest = new TCanvas("ctest", "ZDC test channel", 800, 800);
1585 ctest->cd();
1586 ctest->Divide(2, n_test);
1587 for (int i = 0; i < n_test; i++) {
1588 ctest->cd(i * 2 + 1);
1589 htest[i]->Draw();
1590 ctest->cd(i * 2 + 2);
1591 TestProf[i]->Draw();
1592 gPad->AddExec("exselt", "select_hist()");
1593 }
1594 return;
1595}
PRINT_ZDC_CALIBRATION
#define PRINT_ZDC_CALIBRATION
Definition BmnZDCRaw2Digit.cxx:10
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
i
int i
Definition P4_F32vec4.h:22
m
__m128 m
Definition P4_F32vec4.h:27
log
friend F32vec4 log(const F32vec4 &a)
Definition P4_F32vec4.h:123
exp
friend F32vec4 exp(const F32vec4 &a)
Definition P4_F32vec4.h:122
f
float f
Definition P4_F32vec4.h:21
BmnADCDigit::GetNSamples
UInt_t GetNSamples() const
Definition BmnADCDigit.h:35
BmnADCDigit::GetChannel
UShort_t GetChannel() const
Definition BmnADCDigit.h:33
BmnADCDigit::GetSerial
UInt_t GetSerial() const
Definition BmnADCDigit.h:31
BmnADCDigit::GetUShortValue
UShort_t * GetUShortValue() const
Definition BmnADCDigit.h:37
BmnZDCRaw2Digit::fillCalibrateCluster
int fillCalibrateCluster(TClonesArray *data, Float_t x, Float_t y, Float_t e, Int_t clsize)
Definition BmnZDCRaw2Digit.cxx:707
BmnZDCRaw2Digit::fillEvent
void fillEvent(TClonesArray *data, TClonesArray *zdcdigit)
Definition BmnZDCRaw2Digit.cxx:604
BmnZDCRaw2Digit::drawzdc
void drawzdc(int nohist=0)
Definition BmnZDCRaw2Digit.cxx:1472
BmnZDCRaw2Digit::fillSampleProfilesAll
void fillSampleProfilesAll(TClonesArray *data, Float_t x, Float_t y, Float_t e)
Definition BmnZDCRaw2Digit.cxx:497
BmnZDCRaw2Digit::fillCalibrateAll
int fillCalibrateAll(TClonesArray *data, Float_t x, Float_t y, Float_t e)
Definition BmnZDCRaw2Digit.cxx:996
BmnZDCRaw2Digit::fillCalibrateNumbers
int fillCalibrateNumbers(TClonesArray *data, Float_t x, Float_t y, Float_t e, Int_t ncells, Int_t *numbers)
Definition BmnZDCRaw2Digit.cxx:890
BmnZDCRaw2Digit::fillSampleProfiles
void fillSampleProfiles(TClonesArray *data, Float_t x, Float_t y, Float_t e, Int_t clsize)
Definition BmnZDCRaw2Digit.cxx:329
BmnZDCRaw2Digit::fillAmplitudes
void fillAmplitudes(TClonesArray *data)
Definition BmnZDCRaw2Digit.cxx:664
Bmn_ZDC_map_element::id
unsigned long id
Definition BmnZDCRaw2Digit.h:30
Bmn_ZDC_map_element::ix
unsigned char ix
Definition BmnZDCRaw2Digit.h:34
Bmn_ZDC_map_element::y
float y
Definition BmnZDCRaw2Digit.h:38
Bmn_ZDC_map_element::iy
unsigned char iy
Definition BmnZDCRaw2Digit.h:35
Bmn_ZDC_map_element::x
float x
Definition BmnZDCRaw2Digit.h:37
Bmn_ZDC_map_element::chan
unsigned char chan
Definition BmnZDCRaw2Digit.h:31
Bmn_ZDC_map_element::adc_chan
unsigned char adc_chan
Definition BmnZDCRaw2Digit.h:32
Bmn_ZDC_map_element::size
unsigned char size
Definition BmnZDCRaw2Digit.h:33
Bmn_ZDC_map_element::used
unsigned char used
Definition BmnZDCRaw2Digit.h:36
UniRun::GetRun
static UniRun * GetRun(int period_number, int run_number)
get run from the database
Definition UniRun.cxx:177
UniRun::GetBeamParticle
TString GetBeamParticle()
get beam particle of the current run
Definition UniRun.h:113
MAX_CHANNELS
#define MAX_CHANNELS
Definition BmnZDCRaw2Digit.h:19
MAX_EVENTS
#define MAX_EVENTS
Definition BmnZDCRaw2Digit.h:18
MAX_LOG_CHANNELS
#define MAX_LOG_CHANNELS
Definition BmnZDCRaw2Digit.h:20
MAX_CHANNELS
#define MAX_CHANNELS
Definition BmnZDCRaw2Digit.h:19
MAX_EVENTS
#define MAX_EVENTS
Definition BmnZDCRaw2Digit.h:18
centrality_clusterization.float
float
Definition centrality_clusterization.py:212
check_success.s
s
Definition check_success.py:47
run
-clang-format