1515
1616using namespace o2 ::zdc;
1717
18+ std::array<float , NChannels> RecEventScale::fe = {1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 };
19+ std::array<float , NTDCChannels> RecEventScale::fa = {1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 };
20+
21+ void RecEventScale::reset ()
22+ {
23+ for (int ich = 0 ; ich < NChannels; ich++) {
24+ fe[ich] = 1 ;
25+ }
26+ for (int itdc = 0 ; itdc < NTDCChannels; itdc++) {
27+ fa[itdc] = 1 ;
28+ }
29+ }
30+
31+ void RecEventScale::setGeV ()
32+ {
33+ for (int ich = 0 ; ich < NChannels; ich++) {
34+ fe[ich] = 1 ;
35+ }
36+ for (int itdc = 0 ; itdc < NTDCChannels; itdc++) {
37+ fa[itdc] = 1000 .;
38+ }
39+ }
40+
41+ void RecEventScale::setGeVZN ()
42+ {
43+ for (int ich = IdZNAC; ich <= IdZNASum; ich++) {
44+ fe[ich] = 1 ;
45+ }
46+ for (int ich = IdZNCC; ich <= IdZNCSum; ich++) {
47+ fe[ich] = 1 ;
48+ }
49+ fa[TDCZNAC] = 1000 .;
50+ fa[TDCZNAS] = 1000 .;
51+ fa[TDCZNCC] = 1000 .;
52+ fa[TDCZNCS] = 1000 .;
53+ }
54+
55+ void RecEventScale::setGeVZP ()
56+ {
57+ for (int ich = IdZPAC; ich <= IdZPASum; ich++) {
58+ fe[ich] = 1 ;
59+ }
60+ for (int ich = IdZPCC; ich <= IdZPCSum; ich++) {
61+ fe[ich] = 1 ;
62+ }
63+ fa[TDCZPAC] = 1000 .;
64+ fa[TDCZPAS] = 1000 .;
65+ fa[TDCZPCC] = 1000 .;
66+ fa[TDCZPCS] = 1000 .;
67+ }
68+
69+ void RecEventScale::setTeV ()
70+ {
71+ for (int ich = 0 ; ich < NChannels; ich++) {
72+ fe[ich] = 0.001 ;
73+ }
74+ for (int itdc = 0 ; itdc < NTDCChannels; itdc++) {
75+ fa[itdc] = 1 .;
76+ }
77+ };
78+
79+ void RecEventScale::setTeVZN ()
80+ {
81+ for (int ich = IdZNAC; ich <= IdZNASum; ich++) {
82+ fe[ich] = 0.001 ;
83+ }
84+ for (int ich = IdZNCC; ich <= IdZNCSum; ich++) {
85+ fe[ich] = 0.001 ;
86+ }
87+ fa[TDCZNAC] = 1 .;
88+ fa[TDCZNAS] = 1 .;
89+ fa[TDCZNCC] = 1 .;
90+ fa[TDCZNCS] = 1 .;
91+ };
92+
93+ void RecEventScale::setTeVZP ()
94+ {
95+ for (int ich = IdZPAC; ich <= IdZPASum; ich++) {
96+ fe[ich] = 0.001 ;
97+ }
98+ for (int ich = IdZPCC; ich <= IdZPCSum; ich++) {
99+ fe[ich] = 0.001 ;
100+ }
101+ fa[TDCZPAC] = 1 .;
102+ fa[TDCZPAS] = 1 .;
103+ fa[TDCZPCC] = 1 .;
104+ fa[TDCZPCS] = 1 .;
105+ };
106+
107+ void RecEventScale::setNucleonEnergy (float energy)
108+ { // In GeV
109+ for (int ich = 0 ; ich < NChannels; ich++) {
110+ fe[ich] = 1 . / energy;
111+ }
112+ for (int itdc = 0 ; itdc < NTDCChannels; itdc++) {
113+ fa[itdc] = 1000 . / energy;
114+ }
115+ };
116+
117+ void RecEventScale::setNucleonEnergyZN (float energy)
118+ { // In GeV
119+ for (int ich = IdZNAC; ich <= IdZNASum; ich++) {
120+ fe[ich] = 1 . / energy;
121+ }
122+ for (int ich = IdZNCC; ich <= IdZNCSum; ich++) {
123+ fe[ich] = 1 . / energy;
124+ }
125+ fa[TDCZNAC] = 1000 . / energy;
126+ fa[TDCZNAS] = 1000 . / energy;
127+ fa[TDCZNCC] = 1000 . / energy;
128+ fa[TDCZNCS] = 1000 . / energy;
129+ };
130+
131+ void RecEventScale::setNucleonEnergyZP (float energy)
132+ { // In GeV
133+ for (int ich = IdZPAC; ich <= IdZPASum; ich++) {
134+ fe[ich] = 1 . / energy;
135+ }
136+ for (int ich = IdZPCC; ich <= IdZPCSum; ich++) {
137+ fe[ich] = 1 . / energy;
138+ }
139+ fa[TDCZPAC] = 1000 . / energy;
140+ fa[TDCZPAS] = 1000 . / energy;
141+ fa[TDCZPCC] = 1000 . / energy;
142+ fa[TDCZPCS] = 1000 . / energy;
143+ };
144+
18145void RecEventFlat::init (const std::vector<o2::zdc::BCRecData>* RecBC, const std::vector<o2::zdc::ZDCEnergy>* Energy, const std::vector<o2::zdc::ZDCTDCData>* TDCData, const std::vector<uint16_t >* Info)
19146{
20147 mRecBC = *RecBC;
@@ -143,7 +270,7 @@ int RecEventFlat::at(int ientry)
143270 for (int i = mFirstE ; i < mStopE ; i++) {
144271 auto myenergy = mEnergy [i];
145272 auto ch = myenergy.ch ();
146- ezdc[ch] = myenergy.energy ();
273+ ezdc[ch] = myenergy.energy () * RecEventScale::fe[ch] ;
147274 // Assign implicit event info
148275 if (adcPedOr[ch] == false && adcPedQC[ch] == false && adcPedMissing[ch] == false ) {
149276 adcPedEv[ch] = true ;
@@ -163,7 +290,7 @@ int RecEventFlat::at(int ientry)
163290 isEnd[ch] = true ;
164291 }
165292 TDCVal[ch].push_back (mytdc.val );
166- TDCAmp[ch].push_back (mytdc.amp );
293+ TDCAmp[ch].push_back (mytdc.amp * RecEventScale::fa[ch] );
167294 // Assign implicit event info
168295 if (tdcPedQC[ch] == false && tdcPedMissing[ch] == false ) {
169296 tdcPedOr[ch] = true ;
@@ -386,13 +513,44 @@ float RecEventFlat::yZNC()
386513
387514float RecEventFlat::xZPA ()
388515{
389- // X coordinates of tower centers
390- // Positive because calorimeter is on the right of ZNA when looking
391- // from IP
392- const static float xc[4 ] = {2.8 , 8.4 , 14.0 , 19.6 };
516+ float x, rms;
517+ centroidZPA (x, rms);
518+ return x;
519+ }
520+
521+ float RecEventFlat::xZPC ()
522+ {
523+ float x, rms;
524+ centroidZPC (x, rms);
525+ return x;
526+ }
527+
528+ float RecEventFlat::rmsZPA ()
529+ {
530+ float x, rms;
531+ centroidZPA (x, rms);
532+ return rms;
533+ }
534+
535+ float RecEventFlat::rmsZPC ()
536+ {
537+ float x, rms;
538+ centroidZPC (x, rms);
539+ return rms;
540+ }
541+
542+ void RecEventFlat::centroidZPA (float & x, float & rms)
543+ {
544+ // x coordinates of tower centers
545+ // Negative because ZPA calorimeter is on the left of ZNA when looking from IP
546+ // rms can result from 0 to sqrt((2.8^2+19.6^2)/2-((2.6+19.4)/2)^2) = 8.66
547+ const static float xc[4 ] = {-19.6 , -14.0 , -8.4 , -2.8 };
548+ const static float xq[4 ] = {19.6 * 19.6 , 14.0 * 14.0 , 8.4 * 8.4 , 2.8 * 2.8 };
393549 static float c = 0 ;
394- if (mComputed [2 ]) {
395- return c;
550+ static float d = 0 ;
551+ if (mComputed [2 ] == true ) {
552+ x = c;
553+ rms = d;
396554 }
397555 mComputed [2 ] = true ;
398556 float e[4 ], w[4 ];
@@ -402,35 +560,56 @@ float RecEventFlat::xZPA()
402560 e[3 ] = EZDC (IdZPA4);
403561 if (e[0 ] < -1000000 || e[1 ] < -1000000 || e[2 ] < -1000000 || e[3 ] < -1000000 ) {
404562 c = -std::numeric_limits<float >::infinity ();
405- return c;
563+ d = -std::numeric_limits<float >::infinity ();
564+ x = c;
565+ rms = d;
566+ return ;
406567 }
407568 float sumw = 0 ;
408569 c = 0 ;
570+ d = 0 ;
409571 for (int i = 0 ; i < 4 ; i++) {
410572 if (e[i] < 0 ) {
411- e[i] = 0 ;
573+ w[i] = 0 ;
574+ } else {
575+ // w[i] = std::sqrt(e[i]);
576+ w[i] = e[i];
412577 }
413- w[i] = std::sqrt (e[i]);
414578 sumw += w[i];
415- c += w[i] * xc[i];
579+ c = c + w[i] * xc[i];
580+ d = d + w[i] * xq[i];
416581 }
417582 if (sumw <= 0 ) {
418583 c = -std::numeric_limits<float >::infinity ();
584+ d = -std::numeric_limits<float >::infinity ();
419585 } else {
420- c /= sumw;
586+ c = c / sumw;
587+ d = d / sumw;
588+ d = d - c * c; // variance
589+ if (d >= 0 ) {
590+ d = TMath::Sqrt (d);
591+ } else {
592+ LOGF (error, " %s FOP exception @ %d" , __FILE__, __LINE__);
593+ d = -std::numeric_limits<float >::infinity ();
594+ }
421595 }
422- return c;
596+ x = c;
597+ rms = d;
598+ return ;
423599}
424600
425- float RecEventFlat::xZPC ( )
601+ void RecEventFlat::centroidZPC ( float & x, float & rms )
426602{
427- // X coordinates of tower centers
428- // Negative because calorimeter is on the left of ZNC when looking
429- // from IP
430- const static float xc[4 ] = {-2.8 , -8.4 , -14.0 , -19.6 };
603+ // x coordinates of tower centers
604+ // Positive because ZPC calorimeter is on the right of ZNC when looking from IP
605+ // x can result in a value from 2.8 to 19.6
606+ const static float xc[4 ] = {2.8 , 8.4 , 14.0 , 19.6 };
607+ const static float xq[4 ] = {2.8 * 2.8 , 8.4 * 8.4 , 14.0 * 14.0 , 19.6 * 19.6 };
431608 static float c = 0 ;
609+ static float d = 0 ;
432610 if (mComputed [3 ]) {
433- return c;
611+ x = c;
612+ rms = d;
434613 }
435614 mComputed [3 ] = true ;
436615 float e[4 ], w[4 ];
@@ -440,24 +619,42 @@ float RecEventFlat::xZPC()
440619 e[3 ] = EZDC (IdZPC4);
441620 if (e[0 ] < -1000000 || e[1 ] < -1000000 || e[2 ] < -1000000 || e[3 ] < -1000000 ) {
442621 c = -std::numeric_limits<float >::infinity ();
443- return c;
622+ d = -std::numeric_limits<float >::infinity ();
623+ x = c;
624+ rms = d;
625+ return ;
444626 }
445627 float sumw = 0 ;
446628 c = 0 ;
629+ d = 0 ;
447630 for (int i = 0 ; i < 4 ; i++) {
448631 if (e[i] < 0 ) {
449- e[i] = 0 ;
632+ w[i] = 0 ;
633+ } else {
634+ // w[i] = std::sqrt(e[i]);
635+ w[i] = e[i];
450636 }
451- w[i] = std::sqrt (e[i]);
452637 sumw += w[i];
453638 c += w[i] * xc[i];
639+ d += w[i] * xq[i];
454640 }
455641 if (sumw <= 0 ) {
456642 c = -std::numeric_limits<float >::infinity ();
643+ d = -std::numeric_limits<float >::infinity ();
457644 } else {
458645 c /= sumw;
646+ d /= sumw;
647+ d = d - c * c; // variance
648+ if (d >= 0 ) {
649+ d = TMath::Sqrt (d);
650+ } else {
651+ LOGF (error, " %s FOP exception @ %d" , __FILE__, __LINE__);
652+ d = -std::numeric_limits<float >::infinity ();
653+ }
459654 }
460- return c;
655+ x = c;
656+ rms = d;
657+ return ;
461658}
462659
463660void RecEventFlat::print () const
0 commit comments