KaliVeda  1.13/01
Heavy-Ion Analysis Toolkit
KVDataQualityAudit.cpp
Go to the documentation of this file.
1 #include "KVDataQualityAudit.h"
2 #include <algorithm>
3 
4 
7 
9 {
10  // Add this reconstructed nucleus to the audit.
11 
12  auto tel = telescopes.get_object<idtelescope>(N.GetIdentifyingTelescope()->GetName());
13  if (!tel) {
14  tel = new idtelescope(N.GetIdentifyingTelescope()->GetName());
15  telescopes.Add(tel);
16  }
17  tel->add(N);
18 }
19 
20 
21 
23 
25 {
26  for (auto p : telescopes) {
27  p->Print(opt);
28  }
29 }
30 
31 
32 
38 
40 {
41  // Used to merge audit objects in different ROOT files (using hadd).
42  //
43  // This method is called for the object in the first file, the TCollection
44  // contains the objects with the same name in all the other files.
45 
46  TIter next(l);
47  TObject* o;
48  Long64_t nmerge = 0;
49  while ((o = next())) {
50  nmerge += merge(dynamic_cast<KVDataQualityAudit*>(o));
51  }
52  return nmerge;
53 }
54 
55 
56 
59 
61 {
62  // Add copy of idtelescope data to this audit
63  telescopes.Add(new idtelescope(*idt));
64 }
65 
66 
67 
71 
73 {
74  // Merge other audit object with this one. Return 1 if merge is successful,
75  // 0 in case of problems.
76  if (!other) return 0; // possibly called for wrong object type
77 
78  TIter next_tel(other->GetTelescopeList());
79  idtelescope* idt;
80  Long64_t nmerge = 0;
81  while ((idt = (idtelescope*)next_tel())) {
82  if (HasTelescope(idt))
83  GetTelescope(idt)->merge(idt);
84  else
85  add(idt);
86  ++nmerge;
87  }
88  return nmerge;
89 }
90 
91 
92 
97 
99 {
100  // add an isotopically identified nucleus
101  //
102  // for calibrated particles, we store the minimum energy (=> threshold)
103  A = N.GetA();
104  ++counts;
105  if (N.IsCalibrated() && ((emin < 0) || (N.GetEnergy() < emin))) emin = N.GetEnergy();
106 }
107 
108 
109 
111 
113 {
114  std::cout << "\t\tA=" << A << "\t\t[E>" << emin << "]" << std::endl;
115 }
116 
117 
118 
125 
127 {
128  // add an element (Z-identified nucleus)
129  //
130  // if isotopically identified, add to list of isotopes
131  //
132  // for calibrated particles, we store the minimum energy (=> threshold)
133  if (N.GetZ() > 0) Z = N.GetZ();
134  ++counts;
135  if (N.IsCalibrated() && ((emin < 0) || (N.GetEnergy() < emin))) emin = N.GetEnergy();
136  if (N.IsAMeasured()) {
137  isotopes[N.GetA()].add(N);
138  }
139  else if (N.GetA() > 0) A = N.GetA(); // default mass for Z-only identification
140 }
141 
142 
143 
145 
147 {
148  std::cout << "\tZ=" << (int)Z << "\t\t[E>" << emin << "]";
149  if (A > 0) std::cout << " default A=" << A;
150  if (!HasIsotopes()) {
151  std::cout << std::endl;
152  return;
153  }
154  std::cout << "\t\tIsotopes: " << GetIsotopeList().AsString() << std::endl;
155  for (auto i : isotopes) i.second.print();
156 }
157 
158 
159 
167 
169 {
170  // Merge the informations in element elem with those in this one:
171  //
172  // - we keep the lower of the two threshold values
173  // - we sum the counts for the element
174  // - we merge the lists of isotopes of the two elements
175  // - we keep the default mass, if it has been set
176 
177  if (emin > 0) {
178  if (elem.emin > 0) emin = std::min(emin, elem.emin);
179  }
180  else if (elem.emin > 0) emin = elem.emin;
181  A = std::max(A, elem.A);
182  counts += elem.counts;
183  // look at isotopes in other's list. if any are not in our list they are added.
184  // if any are in both, they are merged
185  std::for_each(std::begin(elem.isotopes), std::end(elem.isotopes),
186  [&](const std::pair<int, isotope>& isotop) {
187  if (HasIsotope(isotop.first)) {
188  isotopes[isotop.first].merge(isotop.second);
189  }
190  else {
191  isotopes[isotop.first] = isotop.second;
192  }
193  }
194  );
195 }
196 
197 
198 
201 
203 {
204  // calculate and return mean mass of isotopes measured for this element
205 
206  if (isotopes.empty()) return 0;
207  double sum{0}, weight{0};
208  std::for_each(std::begin(isotopes), std::end(isotopes),
209  [&](const std::pair<int, isotope>& count) {
210  sum += count.second.counts;
211  weight += count.first * count.second.counts;
212  }
213  );
214  if (sum > 0) return weight / sum;
215  return 0;
216 }
217 
218 
219 
222 
224 {
225  // return max A of isotopes measured for this element
226 
227  if (isotopes.empty()) return 0;
228  auto last = isotopes.rbegin();
229  return (*last).first;
230 }
231 
232 
233 
236 
238 {
239  // return min A of isotopes measured for this element
240 
241  if (isotopes.empty()) return 0;
242  auto first = isotopes.begin();
243  return (*first).first;
244 }
245 
246 
247 
250 
252 {
253  // \returns the smallest recorded threshold (in MeV/u) for isotopic identification for the element
254  if (isotopes.empty()) return -1;
255  double minE = 9999;
256  for (auto& p : isotopes) {
257  double myE = p.second.emin / p.first;
258  if ((myE > 0) && (myE < minE)) minE = myE;
259  }
260  if (minE < 9999) return minE;
261  // uncalibrated particles: all thresholds are -1 MeV
262  return -1;
263 }
264 
265 
266 
269 
271 {
272  // \returns probability distribution for each isotope as a map between isotope A and P(A)
273 
274  double sum{0};
275  std::for_each(std::begin(isotopes), std::end(isotopes),
276  [&](const std::pair<int, isotope>& count) {
277  sum += count.second.counts;
278  }
279  );
280  std::map<int, double> proba;
281  for (auto& p : isotopes) {
282  proba[p.first] = p.second.counts / sum;
283  }
284  return proba;
285 }
286 
287 
288 
290 
292 {
293  elements[N.GetZ()].add(N);
294 }
295 
296 
297 
299 
301 {
302  std::cout << GetName() << " :" << std::endl;
303  std::cout << "Elements: " << GetElementList().AsString() << std::endl;
304  for (auto e : elements) e.second.print();
305 }
306 
307 
308 
313 
315 {
316  // Merge data in idtelescope idt with this one
317 
318  // look at elements in other's list. if any are not in our list they are added.
319  // if any are in both, they are merged
320  std::for_each(std::begin(idt->elements), std::end(idt->elements),
321  [&](const std::pair<int, element>& elem) {
322  if (HasElement(elem.first)) {
323  elements[elem.first].merge(elem.second);
324  }
325  else {
326  elements[elem.first] = elem.second;
327  }
328  }
329  );
330 }
331 
332 
333 
336 
338 {
339  // \returns probability distribution for all elements identified by telescope as a map between Z and P(Z)
340 
341  double sum{0};
342  std::for_each(std::begin(elements), std::end(elements),
343  [&](const std::pair<int, element>& count) {
344  sum += count.second.counts;
345  }
346  );
347  std::map<int, double> proba;
348  for (auto& p : elements) {
349  proba[p.first] = p.second.counts / sum;
350  }
351  return proba;
352 }
353 
354 
355 
358 
360 {
361  // \returns mean Z measured in telescope
362  double sum{0}, weight{0};
363  std::for_each(std::begin(elements), std::end(elements),
364  [&](const std::pair<int, element>& count) {
365  sum += count.second.counts;
366  weight += count.second.counts * count.first;
367  }
368  );
369  return weight / sum;
370 }
371 
372 
373 
376 
378 {
379  // \returns maximum Z measured in telescope
380  auto last = elements.rbegin();
381  return (*last).first;
382 }
383 
384 
387 
389 {
390  // \returns maximum Z for which isotopic masses were measured
391  auto it = elements.rbegin();
392  while (it != elements.rend()) {
393  if ((*it).second.HasIsotopes()) return (*it).first;
394  ++it;
395  }
396  return 0;
397 }
398 
399 
402 
404 {
405  // \returns minimum Z measured in telescope
406  auto first = elements.begin();
407  return (*first).first;
408 }
409 
410 
411 
414 
416 {
417  // \returns mean A measured in telescope
418  double sum{0}, weight{0};
419  std::for_each(std::begin(elements), std::end(elements),
420  [&](const std::pair<int, element>& count) {
421  sum += count.second.counts;
422  weight += count.second.counts * count.second.get_mean_isotopic_mass();
423  }
424  );
425  return weight / sum;
426 }
427 
428 
429 
432 
434 {
435  // \returns maximum A measured in telescope
436  int max = 0;
437  for (auto& e : elements) {
438  max = std::max(max, e.second.get_max_isotopic_mass());
439  }
440  return max;
441 }
442 
443 
444 
447 
449 {
450  // \returns smallest A measured in telescope
451  int min = 999;
452  for (auto& e : elements) {
453  if (auto other_A = e.second.get_min_isotopic_mass())
454  min = std::min(min, other_A);
455  }
456  return min;
457 }
458 
459 
460 
461 
463 
464 
465 
ClassImp(KVPartitionList) void KVPartitionList
Initialisation.
#define e(i)
const char Option_t
#define N
void add(const KVReconstructedNucleus &N)
void Print(Option_t *opt="") const
std::map< int, element > elements
void merge(const idtelescope *)
std::map< int, double > get_element_distribution() const
Audit of experimental data identification and calibrations.
Long64_t merge(KVDataQualityAudit *)
void Print(Option_t *opt="") const
Long64_t Merge(TCollection *)
Bool_t HasTelescope(const TString &tel_name) const
idtelescope * GetTelescope(const TString &tel_name) const
KVUniqueNameList telescopes
const KVSeqCollection * GetTelescopeList() const
void Add(const KVReconstructedNucleus &N)
Add this reconstructed nucleus to the audit.
void add(const idtelescope *)
Add copy of idtelescope data to this audit.
Nuclei reconstructed from data measured by a detector array .
T * get_object(const TString &name) const
virtual void Add(TObject *obj)
virtual const char * GetName() const
long long Long64_t
double get_mean_isotopic_mass() const
calculate and return mean mass of isotopes measured for this element
double get_minimum_isotopic_threshold_mev_per_nuc() const
std::map< int, double > get_isotopic_distribution() const
std::map< int, isotope > isotopes
Double_t counts
watch the alignment !
int get_min_isotopic_mass() const
return min A of isotopes measured for this element
void add(const KVReconstructedNucleus &N)
int get_max_isotopic_mass() const
return max A of isotopes measured for this element
void add(const KVReconstructedNucleus &N)
Double_t counts
watch the alignment !
auto * l