9 : fDistance(distance), fFolding(folding), fGrouping(grouping), fAddback(addback)
11 SetName(
"GriffinAngles");
15 if(settings !=
nullptr) {
19 }
catch(std::out_of_range&) {}
22 }
catch(std::out_of_range&) {}
26 }
catch(std::out_of_range&) {}
28 std::cout <<
"Failed to find user settings in TGRSIOptions, can't get user settings for excluded detectors/crystals" << std::endl;
31 std::cout <<
"Failed to find TGRSIOptions, can't get user settings for excluded detectors/crystals" << std::endl;
36 std::cout <<
DRED <<
"Warning, grouping is only possible if folding is also active. Setting folding to true!" <<
RESET_COLOR << std::endl;
41 std::cout <<
"Creating angles for detectors " <<
fDistance <<
" mm from center of array, " << (
fAddback ?
"" :
"not ") <<
"using addback, " << (
fFolding ?
"" :
"not ") <<
"using folding, and " << (
fGrouping ?
"" :
"not ") <<
"using grouping. Any angles less than " <<
fRounding <<
" degrees apart will be considered the same." << std::endl;
45 for(
int firstDet = 1; firstDet <= 16; ++firstDet) {
47 for(
int firstCry = 0; firstCry < 4; ++firstCry) {
49 for(
int secondDet = 1; secondDet <= 16; ++secondDet) {
51 for(
int secondCry = 0; secondCry < 4; ++secondCry) {
54 if(firstDet == secondDet && (firstCry == secondCry ||
fAddback)) {
continue; }
57 std::cout <<
"det./cry. " << firstDet <<
"/" << firstCry <<
" with " << secondDet <<
"/" << secondCry <<
", at " <<
fDistance <<
" mm = " << angle << std::endl;
77 std::cout <<
"after folding and rounding: angle " << angle <<
", counts " <<
fAngleCount[
static_cast<int>(std::round(angle /
fRounding))] << std::endl;
98 std::ostringstream str;
100 throw std::runtime_error(str.str());
103 std::ostringstream str;
105 throw std::runtime_error(str.str());
108 std::ostringstream str;
110 throw std::runtime_error(str.str());
113 auto end = std::unique(tmp.begin(), tmp.end());
115 std::ostringstream str;
116 str <<
DRED <<
"Custom grouping with " <<
fCustomGrouping.size() <<
" entries supplied, but the last entry is " <<
fCustomGrouping[0] <<
" which does not match the number of unique groups minus one (" <<
end - tmp.begin() - 1 <<
")!" <<
RESET_COLOR;
117 throw std::runtime_error(str.str());
120 fCustomGrouping = {0, 1, 2, 2, 3, 3, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10};
122 fCustomGrouping = {0, 1, 2, 3, 3, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 9, 9, 10, 10, 10, 11, 11};
124 std::ostringstream str;
125 str <<
DRED <<
"No custom grouping supplied and distance is not 110 mm or 145 mm, but " <<
fDistance <<
"!" <<
RESET_COLOR;
126 throw std::runtime_error(str.str());
131 std::set<double> groupedAngles;
133 double previousAngle = *angle;
136 groupedAngles.insert(previousAngle);
139 previousAngle = *angle;
142 groupedAngles.insert(previousAngle);
184#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 20, 0)
186 std::cout <<
"-------------------- original --------------------" << std::endl;
187 for(
int i = 0; i < z0->GetN(); ++i) {
188 std::cout << std::setw(8) << z0->GetPointX(i) <<
": " << std::setw(8) << z0->GetPointY(i) <<
" " << std::setw(8) << z2->GetPointY(i) <<
" " << std::setw(8) << z4->GetPointY(i) << std::endl;
190 std::cout <<
"--------------------------------------------------" << std::endl;
197 auto* angle0 = z0->GetX();
198 auto* angle2 = z2->GetX();
199 auto* angle4 = z4->GetX();
200 for(
int i = 0; i < z0->GetN(); ++i) {
201 if(angle0[i] > 90.) {
202#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 20, 0)
204 std::cout << i <<
": Folding from " << std::setw(8) << z0->GetPointX(i) <<
" " << std::setw(8) << z2->GetPointX(i) <<
" " << std::setw(8) << z4->GetPointX(i);
207 angle0[i] = 180. - angle0[i];
208 angle2[i] = angle0[i];
209 angle4[i] = angle0[i];
210#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 20, 0)
212 std::cout <<
" to " << std::setw(8) << z0->GetPointX(i) <<
" " << std::setw(8) << z2->GetPointX(i) <<
" " << std::setw(8) << z4->GetPointX(i) << std::endl;
220#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 20, 0)
222 std::cout <<
"-------------------- sorted --------------------" << std::endl;
223 for(
int i = 0; i < z0->GetN(); ++i) {
224 std::cout << std::setw(8) << z0->GetPointX(i) <<
": " << std::setw(8) << z0->GetPointY(i) <<
" " << std::setw(8) << z2->GetPointY(i) <<
" " << std::setw(8) << z4->GetPointY(i) << std::endl;
226 std::cout <<
"------------------------------------------------" << std::endl;
230 auto* counts0 = z0->GetY();
231 auto* errors0 = z0->GetEY();
232 auto* counts2 = z2->GetY();
233 auto* errors2 = z2->GetEY();
234 auto* counts4 = z4->GetY();
235 auto* errors4 = z4->GetEY();
236 for(
int i = 1; i < z0->GetN(); ++i) {
237 if(std::abs(angle0[i] - angle0[i - 1]) <
fRounding) {
239 std::cout << i <<
": Adding angles " << angle0[i - 1] <<
" and " << angle0[i] << std::endl;
242 counts0[i - 1] += counts0[i];
243 counts0[i - 1] /= 2.;
244 errors0[i - 1] = TMath::Sqrt(TMath::Power(errors0[i - 1], 2) + TMath::Power(errors0[i], 2)) / 2.;
245 counts2[i - 1] += counts2[i];
246 counts2[i - 1] /= 2.;
247 errors2[i - 1] = TMath::Sqrt(TMath::Power(errors2[i - 1], 2) + TMath::Power(errors2[i], 2)) / 2.;
248 counts4[i - 1] += counts4[i];
249 counts4[i - 1] /= 2.;
250 errors4[i - 1] = TMath::Sqrt(TMath::Power(errors4[i - 1], 2) + TMath::Power(errors4[i], 2)) / 2.;
255 counts0 = z0->GetY();
256 errors0 = z0->GetEY();
257 counts2 = z2->GetY();
258 errors2 = z2->GetEY();
259 counts4 = z4->GetY();
260 errors4 = z4->GetEY();
264 std::cout << i <<
": Not adding angles " << angle0[i - 1] <<
" and " << angle0[i] << std::endl;
269#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 20, 0)
271 std::cout <<
"-------------------- after folding --------------------" << std::endl;
272 for(
int i = 0; i < z0->GetN(); ++i) {
273 std::cout << std::setw(8) << z0->GetPointX(i) <<
": " << std::setw(8) << z0->GetPointY(i) <<
" " << std::setw(8) << z2->GetPointY(i) <<
" " << std::setw(8) << z4->GetPointY(i) << std::endl;
275 std::cout <<
"-------------------------------------------------------" << std::endl;
283 auto* counts0 = z0->GetY();
284 auto* errors0 = z0->GetEY();
285 auto* counts2 = z2->GetY();
286 auto* errors2 = z2->GetEY();
287 auto* counts4 = z4->GetY();
288 auto* errors4 = z4->GetEY();
289 for(
int i = 0; i < z0->GetN(); ++i) {
293#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 20, 0)
295 std::cout << i <<
": Leaving as is " << z0->GetPointX(i) <<
": " << std::setw(8) << z0->GetPointY(i) <<
" " << std::setw(8) << z2->GetPointY(i) <<
" " << std::setw(8) << z4->GetPointY(i) << std::endl;
304#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 20, 0)
306 std::cout << i <<
": Grouping " << z0->GetPointX(i) <<
" and " << z0->GetPointX(i + 1) <<
" from " << std::setw(8) << z0->GetPointY(i) <<
" " << std::setw(8) << z2->GetPointY(i) <<
" " << std::setw(8) << z4->GetPointY(i);
309 counts0[i + 1] += counts0[i];
310 counts0[i + 1] /= 2.;
311 errors0[i + 1] = TMath::Sqrt(TMath::Power(errors0[i], 2) + TMath::Power(errors0[i + 1], 2)) / 2.;
312 counts2[i + 1] += counts2[i];
313 counts2[i + 1] /= 2.;
314 errors2[i + 1] = TMath::Sqrt(TMath::Power(errors2[i], 2) + TMath::Power(errors2[i + 1], 2)) / 2.;
315 counts4[i + 1] += counts4[i];
316 counts4[i + 1] /= 2.;
317 errors4[i + 1] = TMath::Sqrt(TMath::Power(errors4[i], 2) + TMath::Power(errors4[i + 1], 2)) / 2.;
321 counts0 = z0->GetY();
322 errors0 = z0->GetEY();
323 counts2 = z2->GetY();
324 errors2 = z2->GetEY();
325 counts4 = z4->GetY();
326 errors4 = z4->GetEY();
327#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 20, 0)
329 std::cout <<
" to " << z0->GetPointX(i) <<
" with " << std::setw(8) << z0->GetPointY(i) <<
" " << std::setw(8) << z2->GetPointY(i) <<
" " << std::setw(8) << z4->GetPointY(i) << std::endl;
336#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 20, 0)
338 std::cout << i <<
": Grouping " << z0->GetPointX(i) <<
", " << z0->GetPointX(i + 1) <<
", and " << z0->GetPointX(i + 2) <<
" from " << std::setw(8) << z0->GetPointY(i) <<
" " << std::setw(8) << z2->GetPointY(i) <<
" " << std::setw(8) << z4->GetPointY(i);
341 counts0[i + 2] += counts0[i] + counts0[i + 1];
342 counts0[i + 2] /= 3.;
343 errors0[i + 2] = TMath::Sqrt(TMath::Power(errors0[i], 2) + TMath::Power(errors0[i + 1], 2) + TMath::Power(errors0[i + 2], 2)) / 3.;
344 counts2[i + 2] += counts2[i] + counts2[i + 1];
345 counts2[i + 2] /= 3.;
346 errors2[i + 2] = TMath::Sqrt(TMath::Power(errors2[i], 2) + TMath::Power(errors2[i + 1], 2) + TMath::Power(errors0[i + 2], 2)) / 3.;
347 counts4[i + 2] += counts4[i] + counts4[i + 1];
348 counts4[i + 2] /= 3.;
349 errors4[i + 2] = TMath::Sqrt(TMath::Power(errors4[i], 2) + TMath::Power(errors4[i + 1], 2) + TMath::Power(errors0[i + 2], 2)) / 3.;
353 z0->RemovePoint(i + 1);
354 z2->RemovePoint(i + 1);
355 z4->RemovePoint(i + 1);
356 counts0 = z0->GetY();
357 errors0 = z0->GetEY();
358 counts2 = z2->GetY();
359 errors2 = z2->GetEY();
360 counts4 = z4->GetY();
361 errors4 = z4->GetEY();
362#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 20, 0)
364 std::cout <<
" to " << z0->GetPointX(i) <<
" with " << std::setw(8) << z0->GetPointY(i) <<
" " << std::setw(8) << z2->GetPointY(i) <<
" " << std::setw(8) << z4->GetPointY(i) << std::endl;
372#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 20, 0)
374 std::cout <<
"-------------------- after folding and grouping --------------------" << std::endl;
375 for(
int i = 0; i < z0->GetN(); ++i) {
376 std::cout << std::setw(8) << z0->GetPointX(i) <<
": " << std::setw(8) << z0->GetPointY(i) <<
" " << std::setw(8) << z2->GetPointY(i) <<
" " << std::setw(8) << z4->GetPointY(i) << std::endl;
378 std::cout <<
"--------------------------------------------------------------------" << std::endl;