10#define PRINT_ZDC_CALIBRATION 0
26 for (
int i = 0;
i < 6;
i++)
46 for (
int i = 0;
i < 256;
i++)
53 TString dir = getenv(
"VMCWORKDIR");
55 in.open((path + mappingFile).Data());
57 printf(
"Loading ZDC Map from file: %s - file open error!\n", mappingFile.
Data());
60 printf(
"Loading ZDC Map from file: %s\n", mappingFile.
Data());
61 in >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy >> dummy;
63 int ixmin = -1, ixmax = -1, iymin = -1, iymax = -1;
64 int xmin = 10000., xmax = -10000., ymin = 10000., ymax = -10000.;
66 int id, chan, front_chan, size, ix, iy, used;
68 in >> std::hex >>
id >> std::dec >> chan >> front_chan >> size >> ix >> iy >> x >> y >> used;
72 if (size >= 200 && size <= 223) {
82 if (size > 2 || size < 0)
144 char filn[128], tit1[32] = {
"Channel"}, tit2[32] = {
"Calibration"}, tit3[32] = {
"Error"};
145 TString path1 = dir +
"/parameters/zdc/";
148 if (strlen(CalibrationFile.
Data()) == 0) {
150 sprintf(filn,
"%szdc_muon_calibration.txt", path1.
Data());
155 sprintf(filn,
"%s%s", path1.
Data(), CalibrationFile.
Data());
156 fin = fopen(filn,
"r");
164 printf(
" ZDC: Can't open calibration file %s, use default calibration coefficients 1.\n\n", filn);
166 printf(
"Calibration file %s applied.\n", filn);
171 char* readptr = fgets(tit0, 200, fin);
177 int32_t read_res2 = fscanf(fin,
"%s %s %s\n", tit1, tit2, tit3);
178 if ((!readptr) || (read_res1 != 0) || (read_res2 != 0)) {
180 throw std::ios_base::failure(
"Calibration file read error");
183 while (fscanf(fin,
"%d %f %f\n", &ch1, &ca, &cae) == 3) {
184 if (ch1 > 0 && ch1 <=
maxchan) {
192 printf(
"\n ZDC calibration coefficients\n\n%s\t%s\t%s\n\n", tit1, tit2, tit3);
198 path1 = dir +
"/input/";
201 if (strlen(MaxPositionFile.
Data()) == 0) {
202 printf(
"\n ZDC: No cuts on samples wave max position.\n\n");
204 sprintf(filn,
"%s%s", path1.
Data(), MaxPositionFile.
Data());
205 fin = fopen(filn,
"r");
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");
212 while (fscanf(fin,
"%d %f %f\n", &runn, &mpos, &msig) == 3) {
216 printf(
"\n ZDC: Run %d Samples wave Max position range %d - %d\n\n", runn,
MaxPos_min,
MaxPos_max);
252 char tit[128], nam[128];
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]);
272 hsum_sim =
new TH1F(
"hsumsim",
"Sum of theoretical amplitudes", 1000, 0., 100.);
274 hsum_raw =
new TH1F(
"hsumraw",
"Sum of raw amplitudes", 2000, 0., 20000.);
276 hsum =
new TH1F(
"hsumcal",
"Sum of calibrated amplitudes", 1000, 0., 200.);
277 hsum->SetDirectory(0);
279 hxmean =
new TH1F(
"hxmean",
"Shower mean X", 1000, -500., +500.);
281 hymean =
new TH1F(
"hymean",
"Shower mean Y", 1000, -500., +500.);
286 sprintf(nam,
"Average sampling wave, module %d",
zdc_map_element[
i].chan + 1);
322 printf(
"id#\tZDC chan\tADC_chan\tsize\tix\tiy\tx\ty\tused\n");
332 double r = 0., dx = 0., dy = 0.;
342 for (
int i = 0;
i < clsize * clsize;
i++) {
346 int ixhalf = -1, iyhalf = -1, i0 = -1, imin = -1;
348 for (
int ind = 0; ind <
n_rec; ind++) {
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);
383 }
else if (imin >= 0) {
396 printf(
"Invalid beam entry point (%f,%f)!\n", x, y);
400 for (
int i = 1;
i < clsize;
i++) {
404 for (
int i = 1;
i < clsize;
i++) {
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++) {
430 float ped = 0., xm = 0., ym = 0., s = 0.;
438 if (num_test >= 0 && (digit->
GetSerial() & 0xFFFFFF) ==
test_id[num_test]) {
442 j2 = 1 - j1 % 2 + (j1 / 2) * 2;
444 TestProf[num_test]->Fill(j1, samt[j2]);
448 htest[num_test]->Fill(amp);
454 for (ind = 0; ind <
n_rec; ind++)
462 if ((num =
number[ind]) < 0)
467 j2 = 1 - j % 2 + (j / 2) * 2;
472 float signalMin = 0.;
473 float signalMax = 0.;
474 float signalPed = 0.;
475 float signalInt = 0.;
511 for (
int ind = 0; ind <
n_rec; ind++) {
522 printf(
"Warninig: Beam entry point (%f,%f) outside any ZDC module!\n", x, y);
530 for (
int ind = 0; ind <
n_rec; ind++) {
537 float ped = 0., xm = 0., ym = 0., s = 0.;
545 if (num_test >= 0 && (digit->
GetSerial() & 0xFFFFFF) ==
test_id[num_test]) {
549 j2 = 1 - j1 % 2 + (j1 / 2) * 2;
551 TestProf[num_test]->Fill(j1, samt[j2]);
555 htest[num_test]->Fill(amp);
561 for (ind = 0; ind <
n_rec; ind++)
569 if ((num =
number[ind]) < 0)
574 j2 = 1 - j1 % 2 + (j1 / 2) * 2;
579 float signalMin = 0.;
580 float signalMax = 0.;
581 float signalPed = 0.;
582 float signalInt = 0.;
617 if (num_test >= 0 && (digit->
GetSerial() & 0xFFFFFF) ==
test_id[num_test]) {
622 j2 = 1 - j1 % 2 + (j1 / 2) * 2;
624 TestProf[num_test]->Fill(j1, samt[j2]);
628 htest[num_test]->Fill(amp);
634 for (ind = 0; ind <
n_rec; ind++) {
645 float signalMin = 0.;
646 float signalMax = 0.;
647 float signalPed = 0.;
648 float signalInt = 0.;
658 signalMin, signalMax, signalPed, signalInt);
677 if (num_test >= 0 && (digit->
GetSerial() & 0xFFFFFF) ==
test_id[num_test]) {
684 for (ind = 0; ind <
n_rec; ind++)
693 float signalMin = 0.;
694 float signalMax = 0.;
695 float signalPed = 0.;
696 float signalInt = 0.;
714 double r = 0., dx = 0., dy = 0.;
723 for (
int i = 0;
i < clsize * clsize;
i++) {
727 int ixhalf = -1, iyhalf = -1, i0 = -1, imin = -1;
729 for (
int ind = 0; ind <
n_rec; ind++) {
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);
764 }
else if (imin >= 0) {
777 printf(
"Invalid beam entry point (%f,%f)!\n", x, y);
781 for (
int i = 1;
i < clsize;
i++) {
785 for (
int i = 1;
i < clsize;
i++) {
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++) {
824 coef[
i] = (0.5 * ((
i + 1) % 4) + 2.4) / 10.;
829 float sum_raw = 0., sum_sim = 0.;
835 for (ind = 0; ind <
n_rec; ind++)
843 if ((num =
number[ind]) < 0)
846 float signalMin = 0.;
847 float signalMax = 0.;
848 float signalPed = 0.;
849 float signalInt = 0.;
902 double r = 0., dx = 0., dy = 0.;
910 for (
int ind = 0; ind <
n_rec; ind++) {
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);
930 for (
int ind = 0; ind <
n_rec; ind++) {
933 for (
int j = 0; j <
ncells; j++) {
952 float sum_raw = 0., sum_sim = 0.;
957 for (ind = 0; ind <
n_rec; ind++)
965 if ((num =
number[ind]) < 0)
968 float signalMin = 0.;
969 float signalMax = 0.;
970 float signalPed = 0.;
971 float signalInt = 0.;
1010 for (
int ind = 0; ind <
n_rec; ind++) {
1021 printf(
"Warning: Beam entry point (%f,%f) outside any ZDC module!\n", x, y);
1029 for (
int ind = 0; ind <
n_rec; ind++) {
1047 for (ind = 0; ind <
n_rec; ind++)
1055 if ((num =
number[ind]) < 0)
1058 float signalMin = 0.;
1059 float signalMax = 0.;
1060 float signalPed = 0.;
1061 float signalInt = 0.;
1079 double s = 0., e = 0., r = 0., x0 = 0., y0 = 0.;
1115 gMinuit1->
mnexcm(
"SET ERR", arglist, 1, ierflg);
1118 char name[32] = {
""};
1123 sprintf(name,
"Coeff_%03d",
channel1[
i] + 1);
1124 gMinuit1->
mnparm(
i, name, vstart[
i], step[
i], vmin, vmax, ierflg);
1130 gMinuit1->
mnexcm(
"MIGRAD", arglist, 2, ierflg);
1134 Int_t nvpar, nparx, icstat;
1135 gMinuit1->
mnstat(amin, edm, errdef, nvpar, nparx, icstat);
1136 gMinuit1->
mnprin(3, amin);
1151 TString dir = getenv(
"VMCWORKDIR");
1152 TString path = dir +
"/parameters/zdc/";
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");
1158 printf(
"Can't open output calibration file %s\n", filn);
1160 fprintf(fout,
"%s\t%s\t%s\n", tit1, tit2, tit3);
1166 printf(
"%s\t%s\t%s\n", tit1, tit2, tit3);
1172 for (
int j = 0; j <
nevents; j++) {
1186 double chi = 0., le = 0., laA = 0.;
1187 for (
int j = 0; j <
nevents; j++) {
1188 for (
int i = 0;
i < npar;
i++) {
1207 double amp = (
exp(x / r) +
exp(-x / r)) / 2.;
1212 double amp = (
exp(x / r) -
exp(-x / r)) / 2.;
1220 double a1, a2, b1, b2, B;
1222 double bb2 = 6. * ss1;
1223 double bb1 = 1.4 * ss1;
1227 a2 = 0.105 - 0.014 *
log(E);
1228 b1 = bb1 - 0.120 *
log(E);
1229 b2 = bb2 - 0.260 *
log(E);
1232 B = (a1 * b1 + a2 * b2);
1235 fff = (a1 * b1 * (1. -
exp(-h / (2. * b1)) *
ch(x, b1)) + a2 * b2 * (1. -
exp(-h / (2. * b2)) *
ch(x, b2)))
1244 double a1, a2, b1, b2, B;
1246 double bb2 = 6. * ss1;
1247 double bb1 = 1.4 * ss1;
1251 a2 = 0.105 - 0.014 *
log(E);
1252 b1 = bb1 - 0.120 *
log(E);
1253 b2 = bb2 - 0.260 *
log(E);
1256 B = (a1 * b1 + a2 * b2);
1259 fff = (a1 * b1 *
exp(-x / b1) *
sh(h, 2. * b1)) / (h * B) + (a2 * b2 *
exp(-x / b2) *
sh(h, 2. * b2)) / (h * B);
1265 double amp =
PP1(x / 10., h / 10.) +
PP2(x / 10., h / 10.);
1270 vector<Short_t>& scale,
1277 float pedest = 0., ampl = 0., signal_max = 0.;
1280 float ampl_max = -32768.;
1281 float ampl_min = 32768.;
1288 int signalThresh =
digiPar[4];
1291 float scaleSignal = -1.0 / 16.;
1293 float integral = 0.;
1296 for (
int m = 0;
m < ns;
m++) {
1298 float value = scaleSignal * scale[
m];
1300 if (
m >= pedBegin &&
m <= pedEnd) {
1303 pedest /= (pedEnd - pedBegin + 1);
1305 }
else if (
m >= gateBegin &&
m <= gateEnd) {
1306 ampl = value - pedest;
1308 if (ampl > ampl_max) {
1312 if (ampl < ampl_min) {
1321 if (signalType == 0)
1322 signal_max = ampl_max;
1324 signal_max = integral;
1332 if (
fabs(ampl_max - ampl_min) < signalThresh)
1424 float pedest = 0., ampl = 0., signal = 0., signal_max = 0.;
1425 int nsignal = 0, nsignal_max = 0, m1 = 0;
1426 float ampl_max = 0.;
1428 for (
int m = 0;
m < ns;
m++) {
1429 m1 = 1 - (
m % 2) + (
m / 2) * 2;
1431 pedest += (
Short_t)scale[m1] >> 4;
1437 ampl = -(float)((
Short_t)scale[m1] >> 4) + pedest;
1441 if (ampl > ampl_max) {
1449 }
else if (nsignal > nsignal_max) {
1450 signal_max = signal;
1451 nsignal_max = nsignal;
1457 if (nsignal_max > 0 || ampl_max > 0.) {
1459 signal_max = ampl_max;
1474 TCanvas* callbe0 =
new TCanvas(
"callbe0",
"ZDC mapping", 800, 800);
1478 double x1, x2, y1, y2;
1489 b[
i] =
new TBox(x1, y1, x2, y2);
1492 b[
i]->SetLineColor(
kRed);
1493 b[
i]->SetFillStyle(0);
1494 b[
i]->SetLineStyle(0);
1521 TCanvas* calres =
new TCanvas(
"calres",
"ZDC calibration", 800, 800);
1527 gPad->AddExec(
"exselt",
"select_hist()");
1531 gPad->AddExec(
"exselt",
"select_hist()");
1540 TCanvas* csampro =
new TCanvas(
"csampro",
"ZDC sample profiles", 800, 800);
1562 gPad->AddExec(
"exselt",
"select_hist()");
1565 TCanvas* calres =
new TCanvas(
"calres",
"ZDC mean X and Y", 800, 800);
1571 gPad->AddExec(
"exselt",
"select_hist()");
1575 gPad->AddExec(
"exselt",
"select_hist()");
1582 if (
htest[0] == NULL)
1584 TCanvas* ctest =
new TCanvas(
"ctest",
"ZDC test channel", 800, 800);
1588 ctest->
cd(
i * 2 + 1);
1590 ctest->
cd(
i * 2 + 2);
1592 gPad->AddExec(
"exselt",
"select_hist()");
static float profile_err[MAX_EVENTS][MAX_CHANNELS]
static void fcn1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
static float pedestals[MAX_EVENTS][MAX_CHANNELS]
static float profile_amp[MAX_EVENTS][MAX_CHANNELS]
#define PRINT_ZDC_CALIBRATION
static float amp_array[MAX_EVENTS][MAX_CHANNELS]
friend F32vec4 fabs(const F32vec4 &a)
R__EXTERN TRandom * gRandom
uint32_t GetSerial() const
vector< IntType > & GetValue()
uint16_t GetChannel() const
int channel1[MAX_CHANNELS]
int fillCalibrateCluster(TClonesArray *data, Float_t x, Float_t y, Float_t e, Int_t clsize)
float wave2amp(UChar_t ns, vector< Short_t > &scale, Float_t *p, Float_t *sigMin, Float_t *sigMax, Float_t *sigPed, Float_t *sigInt)
TH1F * htest[MAX_LOG_CHANNELS]
UInt_t test_id[MAX_LOG_CHANNELS]
double ch(double x, double r)
TProfile * TestProf[MAX_LOG_CHANNELS]
void FillEvent(TClonesArray *data, TClonesArray *zdcdigit)
void cellWeight(int ievent)
void drawzdc(int nohist=0)
float cal_out[MAX_CHANNELS]
double shower(double x, double h)
Bmn_ZDC_map_element zdc_map_element[256]
void fillSampleProfilesAll(TClonesArray *data, Float_t x, Float_t y, Float_t e)
int fillCalibrateAll(TClonesArray *data, Float_t x, Float_t y, Float_t e)
TProfile * SampleProf[MAX_CHANNELS]
double PP2(double x, double h)
float testwave2amp(UChar_t ns, vector< Short_t > &scale, Float_t *p)
double sh(double x, double r)
float log_amp[MAX_LOG_CHANNELS]
int test_chan[MAX_LOG_CHANNELS]
int fillCalibrateNumbers(TClonesArray *data, Float_t x, Float_t y, Float_t e, Int_t ncells, Int_t *numbers)
float cale_out[MAX_CHANNELS]
virtual ~BmnZDCRaw2Digit()
int channel0[MAX_CHANNELS]
void fillSampleProfiles(TClonesArray *data, Float_t x, Float_t y, Float_t e, Int_t clsize)
float zdc_amp[MAX_CHANNELS]
void fillAmplitudes(TClonesArray *data)
double PP1(double x, double h)
virtual void SetTextAlign(Short_t align=11)
virtual void SetTextColor(Color_t tcolor=1)
virtual void SetTextSize(Float_t tsize=1)
TVirtualPad * cd(Int_t subpadnumber=0) override
virtual Int_t GetParameter(Int_t parNo, Double_t ¤tValue, Double_t ¤tError) const
virtual void SetFCN(void(*fcn)(Int_t &, Double_t *, Double_t &f, Double_t *, Int_t))
virtual void mnexcm(const char *comand, Double_t *plist, Int_t llist, Int_t &ierflg)
virtual void mnstat(Double_t &fmin, Double_t &fedm, Double_t &errdef, Int_t &npari, Int_t &nparx, Int_t &istat)
virtual void mnprin(Int_t inkode, Double_t fval)
virtual void mnparm(Int_t k, TString cnamj, Double_t uk, Double_t wk, Double_t a, Double_t b, Int_t &ierflg)
Int_t GetEntriesFast() const
TObject * At(Int_t idx) const override
virtual void Draw(Option_t *option="")
void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0) override
const char * Data() const
UniRun.h Class for the table: run_.
static UniRun * GetRun(int period_number, int run_number)
get run from the database
TString GetBeamParticle()
get beam particle of the current run
static int use_log_function
RVec< PromoteType< T > > log(const RVec< T > &v)
RVec< PromoteType< T > > exp(const RVec< T > &v)
static int use_log_function
Double_t Sqrt(Double_t x)
Double_t Log10(Double_t x)