1 | #include <stdlib.h> |
2 | #include <iostream> |
3 | #include <map> |
4 | |
5 | #include <JANA/JApplication.h> |
6 | #include <JANA/JEvent.h> |
7 | #include "DLumi.h" |
8 | |
9 | |
10 | |
11 | |
12 | DLumi::DLumi(JEventLoop *loop) |
13 | { |
14 | |
15 | compute_lumi = 1; |
16 | |
17 | |
18 | std::vector<std::map<string,double> > result; |
19 | |
20 | loop->GetCalib("/PHOTON_BEAM/lumi/PS_accept", result); |
21 | |
22 | if ((int)result.size() != DETECTORS) { |
23 | jerr << "Error in DLumi constructor: " |
24 | << "failed to read constants for PS/PSC acceptances " |
25 | << "from calibdb at /PHOTON_BEAM/lumi/PS_accept" << std::endl; |
26 | m_psc_accept[0] = 0.; m_psc_accept[1] = 0.; m_psc_accept[2] = 0.; |
27 | m_ps_accept[0] = 0.; m_ps_accept[1] = 0.; m_ps_accept[2] = 0.; |
28 | } |
29 | else { |
30 | |
31 | m_psc_accept[0] = (result[0])["Norm"]; |
32 | m_psc_accept[1] = (result[0])["Emin"]; |
33 | m_psc_accept[2] = (result[0])["Emax"]; |
34 | |
35 | m_ps_accept[0] = (result[1])["Norm"]; |
36 | m_ps_accept[1] = (result[1])["Emin"]; |
37 | m_ps_accept[2] = (result[1])["Emax"]; |
38 | } |
39 | |
40 | loop->GetCalib("/PHOTON_BEAM/lumi/tagm_tagged", result); |
41 | |
42 | if ((int)result.size() != TAGM_CH) { |
43 | jerr << "Error in DLumi constructor: " |
44 | << "failed to read constants number of TAGM hits " |
45 | << "from calibdb at /PHOTON_BEAM/lumi/tagm_tagged" << std::endl; |
46 | for(int ii = 0; ii < TAGM_CH; ii++) |
47 | tagm_tagged[ii] = 0.; |
48 | |
49 | compute_lumi = 0; |
50 | |
51 | } |
52 | else { |
53 | for (int ii = 0; ii < static_cast<int>(result.size()); ii++) { |
54 | int cnt = (result[ii])["id"]; |
55 | |
56 | if( ((ii + 1) == cnt) && (cnt > 0) && (cnt <= TAGM_CH)) |
57 | tagm_tagged[ii] = (result[ii])["hit"]; |
58 | |
59 | else { |
60 | jerr << "Error in DLumi constructor: " |
61 | << "Invalid counter in the /PHOTON_BEAM/lumi/tagm_tagged table " |
62 | << std::endl; |
63 | tagm_tagged[ii] = 0; |
64 | compute_lumi = 0; |
65 | } |
66 | |
67 | } |
68 | } |
69 | |
70 | |
71 | loop->GetCalib("/PHOTON_BEAM/lumi/tagh_tagged", result); |
72 | |
73 | if ((int)result.size() != TAGH_CH) { |
74 | jerr << "Error in DLumi constructor: " |
75 | << "failed to read constants number of TAGH hits " |
76 | << "from calibdb at /PHOTON_BEAM/lumi/tagh_tagged" << std::endl; |
77 | for(int ii = 0; ii < TAGH_CH; ii++) |
78 | tagh_tagged[ii] = 0.; |
79 | |
80 | compute_lumi = 0; |
81 | |
82 | } |
83 | else { |
84 | for (int ii = 0; ii < static_cast<int>(result.size()); ii++) { |
85 | int cnt = (result[ii])["id"]; |
86 | |
87 | if( ((ii + 1) == cnt) && (cnt > 0) && (cnt <= TAGH_CH)) |
88 | tagh_tagged[ii] = (result[ii])["hit"]; |
89 | else { |
90 | jerr << "Error in DLumi constructor: " |
91 | << "Invalid counter in the /PHOTON_BEAM/lumi/tagh_tagged table " |
92 | << std::endl; |
93 | tagh_tagged[ii] = 0; |
94 | compute_lumi = 0; |
95 | } |
96 | |
97 | } |
98 | } |
99 | |
100 | loop->Get( taghGeomVect ); |
101 | if (taghGeomVect.size() < 1){ |
102 | jerr << "Error in DLumi constructor: " |
103 | << "Cannot get TAGH geometry " |
104 | << endl; |
105 | compute_lumi = 0; |
106 | |
107 | } |
108 | |
109 | loop->Get( tagmGeomVect ); |
110 | if (tagmGeomVect.size() < 1){ |
111 | jerr << "Error in DLumi constructor: " |
112 | << "Cannot get TAGM geometry " |
113 | << endl; |
114 | compute_lumi = 0; |
115 | } |
116 | |
117 | std::map<string,double> result1; |
118 | loop->GetCalib("/PHOTON_BEAM/endpoint_energy", result1); |
119 | if (result1.find("PHOTON_BEAM_ENDPOINT_ENERGY") == result1.end()) { |
120 | std::cerr << "Error in DLumi constructor: " |
121 | << "failed to read photon beam endpoint energy " |
122 | << "from calibdb at /PHOTON_BEAM/endpoint_energy" << std::endl; |
123 | Ebeam = 0; |
124 | } |
125 | else { |
126 | Ebeam = result1["PHOTON_BEAM_ENDPOINT_ENERGY"]; |
127 | } |
128 | |
129 | |
130 | if(compute_lumi) |
131 | CalcLumi(); |
132 | else { |
133 | jerr << "Error in DLumi constructor: " |
134 | << "Cannot compute Luminosity " |
135 | << std::endl; |
136 | } |
137 | |
138 | } |
139 | |
140 | DLumi::~DLumi() { } |
141 | |
142 | |
143 | void DLumi::CalcLumi(){ |
144 | |
145 | double Norm = m_psc_accept[0]; |
146 | double Emin = m_psc_accept[1]; |
147 | double Emax = m_psc_accept[2]; |
148 | |
149 | double Epeak = Emin + Emax; |
150 | |
151 | double accept = 0.; |
152 | |
153 | |
154 | |
155 | for(int ii = 0; ii < TAGM_CH; ii++){ |
156 | |
157 | double tagm_emin = tagmGeomVect[0]->getElow(1); |
158 | double tagm_emax = tagmGeomVect[0]->getEhigh(1); |
159 | |
160 | double tagm_en = (tagm_emin + tagm_emax)/2.; |
161 | |
162 | accept = 0.; |
163 | |
164 | if( (tagm_en > 2*Emin) && (tagm_en < Epeak)){ |
165 | |
166 | accept = 1. - 2.*Emin / tagm_en; |
167 | |
168 | if(accept < 0.) accept = 0.; |
169 | |
170 | } else if( (tagm_en >= Epeak) && (tagm_en < Ebeam)){ |
171 | |
172 | accept = 2.*Emax / tagm_en - 1.; |
173 | |
174 | if(accept < 0.) accept = 0.; |
175 | } else |
176 | accept = 0.; |
177 | |
178 | tagm_lumi[ii] = tagm_tagged[ii]*Norm*accept; |
179 | |
180 | } |
181 | |
182 | |
183 | |
184 | for(int ii = 0; ii < TAGH_CH; ii++){ |
185 | |
186 | double tagh_emin = taghGeomVect[0]->getElow(1); |
187 | double tagh_emax = taghGeomVect[0]->getEhigh(1); |
188 | |
189 | double tagh_en = (tagh_emin + tagh_emax) / 2.; |
190 | |
191 | accept = 0.; |
| Value stored to 'accept' is never read |
192 | |
193 | if( (tagh_en > 2*Emin) && (tagh_en < Epeak)){ |
194 | |
195 | accept = 1. - 2.*Emin / tagh_en; |
196 | |
197 | } else if( (tagh_en >= Epeak) && (tagh_en < Ebeam)){ |
198 | |
199 | accept = 2.*Emax / tagh_en - 1.; |
200 | |
201 | if(accept < 0.) accept = 0.; |
202 | |
203 | } else |
204 | accept = 0.; |
205 | |
206 | tagh_lumi[ii] = tagh_tagged[ii]*Norm*accept; |
207 | |
208 | } |
209 | |
210 | |
211 | } |
212 | |
213 | void DLumi::PrintLumi(){ |
214 | |
215 | std::cout << std::endl; |
216 | std::cout << " Lumi for TAGM counters (nb) " << std::endl; |
217 | std::cout << std::endl; |
218 | |
219 | for(int ii = 0; ii < TAGM_CH; ii++) |
220 | std::cout << " CH = " << ii + 1 << " Lumi = " << tagm_lumi[ii] << std::endl; |
221 | |
222 | |
223 | std::cout << std::endl; |
224 | std::cout << " Lumi for TAGH counters (nb) " << std::endl; |
225 | std::cout << std::endl; |
226 | |
227 | for(int ii = 0; ii < TAGH_CH; ii++) |
228 | std::cout << " CH = " << ii + 1 << " Lumi = " << tagh_lumi[ii] << std::endl; |
229 | } |
230 | |
231 | void DLumi::SaveLumi(){ |
232 | |
233 | } |
234 | |