KaliVeda  1.13/01
Heavy-Ion Analysis Toolkit
KVIDTelescope.cpp
Go to the documentation of this file.
1 /***************************************************************************
2 $Id: KVIDTelescope.cpp,v 1.52 2009/05/05 15:54:04 franklan Exp $
3 Author : $Author: franklan $
4  KVIDTelescope.cpp - description
5  -------------------
6  begin : Wed Jun 18 2003
7  copyright : (C) 2003 by J.D Frankland
8  email : frankland@ganil.fr
9  ***************************************************************************/
10 
11 /***************************************************************************
12  * *
13  * This program is free software; you can redistribute it and/or modify *
14  * it under the terms of the GNU General Public License as published by *
15  * the Free Software Foundation; either version 2 of the License, or *
16  * (at your option) any later version. *
17  * *
18  ***************************************************************************/
19 #include "TROOT.h"
20 #include "KVIDTelescope.h"
21 #include "KVTelescope.h"
22 #include "KVGroup.h"
23 #include "KVNucleus.h"
24 #include "KVReconstructedNucleus.h"
25 #include "KVIDGraph.h"
26 #include "KVIDGrid.h"
27 #include "Riostream.h"
28 #include "TPluginManager.h"
29 #include "KVMultiDetArray.h"
30 #include "KVDataSet.h"
31 #include "KVIDGridManager.h"
32 #include "KVIDZALine.h"
33 #include "KVIDCutLine.h"
34 #include "KVIdentificationResult.h"
35 #include "TMath.h"
36 #include "TClass.h"
37 #include "TH2.h"
38 #include "KVParticleCondition.h"
39 
41 #include <KVIDZAGrid.h>
42 
43 using namespace std;
44 
47 
48 
50 
52  : fGroup(nullptr)
53 {
54  init();
55 }
56 
57 
58 
61 
63 {
64  //default init
67 }
68 
69 
70 
102 
104 {
105  // Default initialisation for ID telescopes.
106  // If telescope has at least 1 grid then it is ready to identify
107  // particles after initialising the grid(s) (kReadyForID=true);
108  // otherwise kReadyForID is set to kFALSE, unless the current dataset (if defined)
109  // has been declared to have no associated identification/calibration parameters,
110  // in which case kReadyForID is by default set to kTRUE (for filtering simulations).
111  //
112  // In order to enable mass identification for certain telescopes without a dedicated
113  // implementation (e.g. for simulating array response), put the following lines
114  // in your .kvrootrc:
115  //
116  // [dataset].[telescope label].MassID: yes
117  //
118  // If you want to limit mass identification to certain values of Z and/or A,
119  // add the following line:
120  //
121  // [dataset].[telescope label].MassID.Validity: [expression]
122  //
123  // where [expression] is some valid C++ boolean expression involving Z and/or A,
124  // for example
125  //
126  // [dataset].[telescope label].MassID.Validity: (Z>3)&&(A<20)
127  //
128  //For identifications using more than one grid, the default behaviour is to try identification
129  //with each grid in turn until a successful identification is obtained. The order in which
130  //the grids should be tried should be specified by a variable with the following format:
131  //
132  //~~~~~~~~~~~~~~~~
133  //[Dataset].[telescope label].GridOrder: [Grid1],[Grid2],...
134  //~~~~~~~~~~~~~~~~
135 
137 
138  // for datasets with no calib/ident infos, all id telescopes work
139  if (gDataSet && !gDataSet->HasCalibIdentInfos()) {
141  }
142  else { // for datasets with calib/ident infos, we need a grid & all detectors working
143  // looping over detectors to check they are working
144  // if one of them is not -> set kReadyForID to false
145  TIter it(GetDetectors());
146  KVDetector* det = 0;
147  while ((det = (KVDetector*)it())) if (!det->IsOK()) {
149  return;
150  }
151 
152  if (GetIDGrid()) {
153  KVIDGraph* gr;
154  TIter it(GetListOfIDGrids());
155  bool ok = kTRUE;
156  KVUniqueNameList tmp_list;// for re-ordering grids
157  bool mass_id = false;
158  while ((gr = (KVIDGraph*)it())) {
159  tmp_list.Add(gr);
160  if (gr->HasMassIDCapability()) mass_id = true;
161  gr->Initialize();
162  // make sure both x & y axes' signals are well set up
163  if (!fGraphCoords[gr].fVarX || !fGraphCoords[gr].fVarY) {
164  ok = kFALSE;
165  Warning("Initialize",
166  "ID tel. %s: grid %s has undefined VarX(%s:%p) or VarY(%s:%p) - WILL NOT USE",
167  GetName(), gr->GetName(), gr->GetVarX(), fGraphCoords[gr].fVarX, gr->GetVarY(), fGraphCoords[gr].fVarY);
168  }
169  }
170  // set to true if at least one grid can provide mass identification
171  SetHasMassID(mass_id);
172  // if more than one grid, need to re-order them according to [Dataset].[telescope label].GridOrder
173  if (GetListOfIDGrids()->GetEntries() > 1 && gDataSet) {
174  KVString grid_list = gDataSet->GetDataSetEnv(Form("%s.GridOrder", GetLabel()));
175  ok = kFALSE;
176  if (grid_list == "")
177  Warning("Initialize", "ID telescope %s has %d grids but no %s variable defined",
178  GetName(), GetListOfIDGrids()->GetEntries(), Form("%s.GridOrder", GetLabel()));
179  else if (grid_list.GetNValues(",") != GetListOfIDGrids()->GetEntries())
180  Warning("Initialize", "ID telescope %s has %d grids but %d grids appear in variable %s",
181  GetName(), GetListOfIDGrids()->GetEntries(), grid_list.GetNValues(","), Form("%s.GridOrder", GetLabel()));
182  else {
183  fIDGrids.Clear();
184  grid_list.Begin(",");
185  while (!grid_list.End()) {
186  auto gr_name = grid_list.Next();
187  auto gr_ob = tmp_list.FindObject(gr_name);
188  if (!gr_ob) {
189  Info("Initialize", "IDtel=%s grid %s missing", GetName(), gr_name.Data());
190  }
191  fIDGrids.Add(gr_ob);
192  }
193  ok = kTRUE;
194  }
195  }
196  if (ok) SetBit(kReadyForID);
197  }
198  }
199 
200  if (gDataSet) {
202  KVString valid;
203  if ((valid = gDataSet->GetDataSetEnv(Form("%s.MassID.Validity", GetLabel()), "")) != "") {
204  valid.ReplaceAll("Z", "_NUC_->GetZ()");
205  valid.ReplaceAll("A", "_NUC_->GetA()");
206  fMassIDValidity.reset(new KVParticleCondition(valid));
207  }
208  }
209 }
210 
211 
212 
221 
223 {
224  // Add a detector to the telescope.
225  //
226  // Detectors must be added in the order they will be hit by impinging particles,
227  // with the last detector being the one particles stopped in the telescope will stop in.
228  // i.e. dE1, dE2, ..., Eres
229  //
230  // Update name of telescope to "ID_[name of 1st detector]_[name of 2nd detector]_ ... _[name of last detector]"
231 
232  if (d) {
233  fDetectors.Add(d);
234  d->AddIDTelescope(this);
235  if (GetSize() > 1) {
236  TString old_name = GetName();
237  old_name += Form("_%s", GetDetectors()->Last()->GetName());
238  SetName(old_name);
239  }
240  else SetName(Form("ID_%s", GetDetector(1)->GetName()));
241  }
242  else {
243  Warning("AddDetector", "Called with null pointer");
244  }
245 }
246 
247 
248 
252 
254 {
255  // print out telescope structure
256  //if opt="fired" only fired detectors are printed
257 
258  TIter next(GetDetectors());
259  KVDetector* obj;
260 
261  if (!strcmp(opt, "fired")) {
262  while ((obj = (KVDetector*) next())) {
263 
264  if (obj->Fired() || obj->GetEnergy())
265  obj->Print("data");
266  }
267  }
268  else {
269  cout << "\n" << opt << "Structure of KVIDTelescope object: " <<
270  GetName() << " " << GetType() << endl;
271  cout << opt <<
272  "--------------------------------------------------------" <<
273  endl;
274  while ((obj = (KVDetector*) next())) {
275  cout << opt << "Detector: " << obj->GetName() << endl;
276  }
277  }
278 }
279 
280 
281 
282 
285 
287 {
288  // Return a pointer to the detector in the telescope with the name "name".
289 
290  KVDetector* tmp = (KVDetector*) GetDetectors()->FindObject(name);
291  if (!tmp)
292  Warning("GetDetector(const Char_t *name)",
293  "Detector %s not found in telescope %s", name, GetName());
294  return tmp;
295 }
296 
297 
298 
299 
301 
303 {
304  return fGroup;
305 }
306 
307 
308 
309 
311 
313 {
314  fGroup = kvg;
315 }
316 
317 
318 
320 
322 {
323  return (GetGroup() ? GetGroup()->GetNumber() : 0);
324 }
325 
326 
327 
328 
337 
339  Double_t Emax, Double_t Estep)
340 {
341  //For a given nucleus, generate a TGraph representing the line in the deltaE-E
342  //plane of the telescope which can be associated with nuclei of this (A,Z) with total
343  //incident energies between the two limits.
344  //NOTE: if there are other absorbers/detectors placed before this telescope,
345  //the energy losses of the particle in these will be taken into account.
346  //If the step in energy is not given, it is taken equal to 100 equal steps between min and max.
347  //The TGraph should be deleted by the user after use.
348 
349 
350  if (!Estep)
351  Estep = (Emax - Emin) / 100.;
352 
353  Int_t nsteps = 1 + (Int_t)((Emax - Emin) / Estep);
354 
355  if (nsteps < 1)
356  return 0;
357 
358  Double_t* y = new Double_t[nsteps];
359  Double_t* x = new Double_t[nsteps];
360  Int_t step = 0;
361 
362  //get list of all detectors through which particle must pass in order to reach
363  //2nd member of ID Telescope
364  TList* detectors =
366  //detectors->ls();
367  TIter next_det(detectors);
368  //cout << "nsteps =" << nsteps << endl;
369 
370  for (Double_t E = Emin; E <= Emax; E += Estep) {
371  //Set energy of nucleus
372  nuc->SetEnergy(E);
373  //cout << "Einc=" << E << endl;
374 
375  //Calculate energy loss in each member and stock in arrays x & y
376  //first member
377  KVDetector* det = 0;
378  x[step] = y[step] = -1;
379  while ((det = (KVDetector*) next_det())) {
380  //det->Print();
381  Double_t eloss = det->GetELostByParticle(nuc);
382  if (det == GetDetector(1))
383  y[step] = eloss;
384  else if (det == GetDetector(2))
385  x[step] = eloss;
386  Double_t E1 = nuc->GetEnergy() - eloss;
387  nuc->SetEnergy(E1);
388  //cout << "Eloss=" << eloss << endl;
389  //cout << "Enuc=" << nuc->GetEnergy() << endl;
390  if (E1 < 1.e-3) break;
391  }
392 
393  //cout << "step = " << step << " x = " << x[step] << " y = " << y[step] << endl;
394 
395  //make sure that some energy is lost in each member
396  //otherwise miss a step and reduce number of points in graph
397  if (x[step] > 0 && y[step] > 0) {
398  step++;
399  }
400  else {
401  nsteps--;
402  }
403 
404  //cout << "nsteps =" << nsteps << endl;
405  //reset iterator ready for next loop on detectors
406  next_det.Reset();
407  }
408  TGraph* tmp = 0;
409  if (nsteps > 1)
410  tmp = new TGraph(nsteps, x, y);
411  delete[]x;
412  delete[]y;
413  return tmp;
414 }
415 
416 
417 
418 
441 
443 {
444  //Default identification method.
445  //
446  //Works for ID telescopes for which one or more identification grids are defined, each
447  //with VARX/VARY parameters corresponding to a KVDetectorSignal or KVDetectorSignalExpression
448  //associated with one or other of the detectors constituting the telescope.
449  //
450  //For identifications using more than one grid, the default behaviour is to try identification
451  //with each grid in turn until a successful identification is obtained. The order in which
452  //the grids should be tried should be specified by a variable with the following format:
453  //
454  //~~~~~~~~~~~~~~~~
455  //[Dataset].[Array Name].[ID type].GridOrder: [Grid1],[Grid2],...
456  //~~~~~~~~~~~~~~~~
457  //
458  //where the name of each grid is given as "VARY_VARX". If no variable defining the order is found,
459  //the grids will be tried in the order they were found in the file containing the grids for this
460  //telescope.
461  //
462  // The KVIdentificationResult is first Clear()ed; then it is filled with IDtype = GetType()
463  // of this identification telescope, IDattempted = true, and the results of the identification
464  // procedure.
465 
466  idr->Clear();
467  idr->IDattempted = true;
468  idr->SetIDType(GetType());
469 
470  KVIDGraph* grid;
471  TIter it(GetListOfIDGrids());
472  while ((grid = (KVIDGraph*)it())) { //loop over grids in order given by [Dataset].[Array Name].[ID type].GridOrder:
473  Double_t de, e;
474  GetIDGridCoords(e, de, grid, x, y);
475  idr->SetGridName(grid->GetName());
476  if (grid->IsIdentifiable(e, de, &idr->Rejecting_Cut)) {
477  grid->Identify(e, de, idr);
478  if (idr->IDOK) break; // stop on first successful identification
479  }
480  else {
481  // particle rejected by cut in grid. idr->Rejecting_Cut contains its name.
482  idr->IDOK = kFALSE;
484  }
485  }
486  idr->IDcode = GetIDCode();
487 
488  return kTRUE;
489 }
490 
491 
492 
493 
529 
531 {
532  // Add an identification grid to the list of grids used by this telescope.
533  //
534  // If the grid's VARX and VARY parameters are set and contain the names of valid
535  // detector signals (see formatting rules below) they will be used by
536  // GetIDGridXCoord() and GetIDGridYCoord() to return the coordinates
537  // needed to perform particle identification using the grid.
538  //
539  // The name of the grid is set to "VARY_VARX" (just the signal names, not the detector
540  // label part - see below). This value will be stored in the
541  // KVIdentificationResult corresponding to an attempted identification of a
542  // KVReconstructedNucleus by this grid.
543  //
544  // VARX/VARY Formatting
545  //
546  // To be valid, grid VARX/Y parameters should be set as follows:
547  //
548  //~~~~~~~~~~~~~~~~~~
549  // [signal name]
550  // or [det_label]::[signal name]
551  //~~~~~~~~~~~~~~~~~~
552  //
553  // where
554  //
555  //~~~~~~~~~~~~~~~~~~
556  // [det_label] (optional) = detector label i.e. string returned by KVDetector::GetLabel()
557  // method for detector. By default, VARX is assumed to be the Eres detector
558  // or last detector and VARY the DE detector or first detector
559  // [signal_name] = name of a signal defined for the detector, possibly depending
560  // on availability of calibration
561  //
562  // To see all available signals for a detector, use
563  //
564  // KVDetector::GetListOfDetectorSignals()
565  //~~~~~~~~~~~~~~~~~~
566 
567  if (grid) {
568  fIDGrids.Add(grid);
569  KVString det_labels_x, det_labels_y;
570  KVDetectorSignal* xx = GetSignalFromGridVar(grid->GetVarX(), "X", det_labels_x);
571  KVDetectorSignal* yy = GetSignalFromGridVar(grid->GetVarY(), "Y", det_labels_y);
572  GraphCoords gc;
573  gc.fVarX = xx;
574  gc.fDetLabelsX = det_labels_x;
575  gc.fVarY = yy;
576  gc.fDetLabelsY = det_labels_y;
577  fGraphCoords[grid] = gc;
578  TString grid_name;
579  if (xx && yy) {
580  grid_name.Form("%s_%s", yy->GetName(), xx->GetName());
581  grid->SetName(grid_name);
582  }
583  }
584 }
585 
586 
587 
633 
635 {
636  // Deduce & return pointer to detector signal from grid VARX/VARY parameter
637  //
638  // To be valid, grid VARX/Y parameters should be set using mathematical expressions
639  // which use the following references to detector signals for the telescope:
640  //
641  //~~~~~~~~~~~~~~~~~~
642  // [signal name]
643  // OR [det_label]::[signal name]
644  //~~~~~~~~~~~~~~~~~~
645  //
646  // where
647  //
648  //~~~~~~~~~~~~~~~~~~
649  // [signal_name] = name of a signal defined for 1 of the detectors of the telescope
650  // [det_label] = optional detector label i.e. string returned by
651  // KVDetector::GetLabel() method for detector
652  //~~~~~~~~~~~~~~~~~~
653  //
654  // If `[det_label]` is not given, we assume for `VARX` the last (E) detector,
655  // while for `VARY` we assume the first (dE) detector. If this telescope has only
656  // one detector, we use it for both variables.
657  //
658  // To see all available signals for a detector, use
659  //
660  //~~~~~~~~~~~~~~~~~~
661  // KVDetector::GetListOfDetectorSignals()
662  //~~~~~~~~~~~~~~~~~~
663  //
664  // #### Example ####
665  // Imagine a telescope which combines 2 detectors, with labels SI and CSI (in that order, i.e. SI is the dE detector,
666  // CSI is the residual energy detector). The following cases are valid:
667  //
668  //~~~~~
669  // VARX = Energy : use 'Energy' signal of CSI detector (default for VARX)
670  // VARY = ADC : use 'ADC' signal of SI detector (default for VARY)
671  // VARY = CSI::Energy : use 'Energy' signal of CSI detector (overrides default for VARY, SI)
672  // VARX = (Q3-Q3Fast)/(0.8*Q3) : uses 'Q3' and 'Q3Fast' signals of CSI detector (default for VARX)
673  // VARY = 1.5*SI::ADC - CSI::Q3Fast/CSI::Q3 : combination of signals from both detectors
674  //~~~~~
675  //
676  // The 'det_labels' string will be filled with a comma-separated list of the labels of each detector used
677  // in the expressions.
678  //
679  // \note if any signals are not defined, they will be evaluated as zero
680 
681  if (var == "") {
682  Warning("GetSignalFromGridVar",
683  "No VAR%s defined for grid for telescope %s. KVIDTelescope-derived class handling identification must override GetIDMapX/GetIDMapY",
684  axe.Data(), GetName());
685  return nullptr;
686  }
687  KVDetector* det = nullptr;
688  KVDetectorSignal* ds(nullptr);
689  KVString sig_type;
690 
691  // check if VARX/Y is a mathematical expression
692  Bool_t is_expression = var.GetNValues("+-*/()") > 1;
693  // examine each term in expression (this will split the elements in a mathematical expression,
694  // or just take the whole expression if no math operators are present)
695  var.Begin("+-*/()");
696  bool explicit_det_ref = false;
697  bool multidetexpr = false;
698  KVString det_label;
699  while (!var.End()) {
700  KVString t = var.Next();
701  // check if we have an explicit reference to a detector
702  if (t.Contains("::")) {
703  t.Begin("::");
704  det_label = t.Next();
705  auto _det = (KVDetector*)GetDetectors()->FindObjectByLabel(det_label);
706  if (_det) {
707  explicit_det_ref = true;
708  sig_type = t.Next();
709  // check signal is defined for detector
710  if (_det->GetDetectorSignal(sig_type)) {
711  if (det_labels.Length()) {
712  if (!det_labels.Contains(det_label)) det_labels += Form(",%s", det_label.Data());
713  }
714  else
715  det_labels = det_label;
716  }
717  }
718  if (det && (_det != det)) multidetexpr = true; // more than 1 detector is referenced
719  det = _det;
720  }
721  }
722  // treat all cases with explicit detector references
723  if (explicit_det_ref) {
724  if (!multidetexpr) {
725  // only 1 detector is directly referenced
726  if (!is_expression) {
727  // it is not a mathematical expression: this is just the name of a detector signal
728  return det->GetDetectorSignal(sig_type);
729  }
730  else {
731  // math expression with explicit reference to signal(s) of 1 detector
732  sig_type = var;
733  // remove all explicit references to detector: 'det' pointer has been set
734  sig_type.ReplaceAll(Form("%s::", det_label.Data()), "");
735  // expression will be handled below
736  }
737  }
738  else {
739  // use specific constructor for explicit detector references
740  ds = new KVDetectorSignalExpression(var, var, GetDetectors());
741  if (!ds->IsValid()) {
742  delete ds;
743  ds = nullptr;
744  Error("GetSignalFromGridVar",
745  "Problem initialising ID-grid %s coordinate for telescope %s."
746  " Check definition of VAR%s for grid (=%s)",
747  axe.Data(), GetName(), axe.Data(), var.Data());
748  return ds;
749  }
751  return ds;
752  }
753  }
754  else {
755  // no explicit reference to detectors - by default, use detector (1) for Y axis,
756  // and the last detector for X axis
757  if (axe == "Y" || GetSize() == 1) det = GetDetector(1);
758  else det = GetDetector(GetSize());
759  sig_type = var;
760  det_labels = det->GetLabel();
761  }
762  ds = det->GetDetectorSignal(sig_type);
763  if (!ds) {
764  // sig_type is not the name of a known signal: assume it is an expression using known signal names
765  if (!det->AddDetectorSignalExpression(sig_type, sig_type)) {
766  Error("GetSignalFromGridVar",
767  "Problem initialising ID-grid %s coordinate for telescope %s. Request for unknown signal %s for detector %s. Check definition of VAR%s for grid (=%s)",
768  axe.Data(), GetName(), sig_type.Data(), det->GetName(), axe.Data(), var.Data());
769  ds = nullptr;
770  }
771  else
772  ds = det->GetDetectorSignal(sig_type);
773  }
774  return ds;
775 }
776 
777 
778 
780 
782 {
784  if (!cl || !cl->InheritsFrom("KVIDZAGrid")) cl = TClass::GetClass("KVIDZAGrid");
785  KVIDGrid* idgrid = (KVIDGrid*)cl->New();
786 
787  idgrid->AddIDTelescope(this);
788  idgrid->SetOnlyZId(onlyZ);
789  idgrid->SetRunList("1-10000");
790 
791  KVIDCutLine* B_line = (KVIDCutLine*)idgrid->Add("OK", "KVIDCutLine");
792  Int_t npoi_bragg = 0;
793  B_line->SetName("Bragg_line");
794  B_line->SetAcceptedDirection("right");
795 
796  return idgrid;
797 }
798 
799 
800 
808 
809 void KVIDTelescope::addLineToGrid(KVIDGrid* idgrid, int zz, int aa, int npoints)
810 {
811 
812  //loop over energy
813  //first find :
814  // ****E1 = energy at which particle passes 1st detector and starts to enter in the 2nd one****
815  // E2 = energy at which particle passes the 2nd detector
816  //then perform npoints calculations between these two energies and use these
817  //to construct a KVIDZALine
818 
819  double xfactor = 1.;
820 
821  KVNucleus part;
822 
823  KVDetector* det_de = GetDetector(1);
824  KVDetector* det_eres = GetDetector(2);
825 
826  Double_t SeuilE = 0.1;
827 
828  part.SetZ(zz);
829  part.SetA(aa);
830 
831  Double_t E1, E2;
832  //find E1
833  //go from SeuilE MeV to det_de->GetEIncOfMaxDeltaE(part.GetZ(),part.GetA()))
834  Double_t E1min = SeuilE, E1max = det_de->GetEIncOfMaxDeltaE(zz, aa);
835  E1 = (E1min + E1max) / 2.;
836 
837  while ((E1max - E1min) > SeuilE) {
838 
839  part.SetEnergy(E1);
840  det_de->Clear();
841  det_eres->Clear();
842 
843  det_de->DetectParticle(&part);
844  det_eres->DetectParticle(&part);
845  if (det_eres->GetEnergy() > SeuilE) {
846  //particle got through - decrease energy
847  E1max = E1;
848  E1 = (E1max + E1min) / 2.;
849  }
850  else {
851  //particle stopped - increase energy
852  E1min = E1;
853  E1 = (E1max + E1min) / 2.;
854  }
855  }
856 
857  //add point to Bragg line
858  Double_t dE_B = det_de->GetMaxDeltaE(zz, aa);
859  Double_t E_B = det_de->GetEIncOfMaxDeltaE(zz, aa);
860  Double_t Eres_B = det_de->GetERes(zz, aa, E_B);
861 
862  KVIDCutLine* B_line = (KVIDCutLine*)idgrid->GetCut("Bragg_line");
863  if (B_line) B_line->SetPoint(B_line->GetN(), Eres_B, dE_B);
864 
865  //find E2
866  //go from E1 MeV to maximum value where the energy loss formula is valid
867  Double_t E2min = E1, E2max = det_eres->GetEmaxValid(part.GetZ(), part.GetA());
868  E2 = (E2min + E2max) / 2.;
869 
870  while ((E2max - E2min > SeuilE)) {
871 
872  part.SetEnergy(E2);
873  det_de->Clear();
874  det_eres->Clear();
875 
876  det_de->DetectParticle(&part);
877  det_eres->DetectParticle(&part);
878  if (part.GetEnergy() > SeuilE) {
879  //particle got through - decrease energy
880  E2max = E2;
881  E2 = (E2max + E2min) / 2.;
882  }
883  else {
884  //particle stopped - increase energy
885  E2min = E2;
886  E2 = (E2max + E2min) / 2.;
887  }
888  }
889  E2 *= xfactor;
890  if ((!strcmp(det_eres->GetType(), "CSI")) && (E2 > 5000)) E2 = 5000;
891  // printf("z=%d a=%d E1=%lf E2=%lf\n",zz,aa,E1,E2);
892  KVIDZALine* line = (KVIDZALine*)idgrid->Add("ID", "KVIDZALine");
893  if (TMath::Even(zz)) line->SetLineColor(4);
894  line->SetZ(zz);
895  line->SetA(aa);
896 
897  Double_t logE1 = TMath::Log(E1);
898  Double_t logE2 = TMath::Log(E2);
899  Double_t dLog = (logE2 - logE1) / (npoints - 1.);
900 
901  for (Int_t i = 0; i < npoints; i++) {
902  // Double_t E = E1 + i*(E2-E1)/(npoints-1.);
903  Double_t E = TMath::Exp(logE1 + i * dLog);
904 
905  Double_t Eres = 0.;
906  Int_t niter = 0;
907  while (Eres < SeuilE && niter <= 20) {
908  det_de->Clear();
909  det_eres->Clear();
910 
911  part.SetEnergy(E);
912 
913  det_de->DetectParticle(&part);
914  det_eres->DetectParticle(&part);
915 
916  Eres = det_eres->GetEnergy();
917  E += SeuilE;
918  niter += 1;
919  }
920  if (!(niter > 20)) {
921  Double_t dE = det_de->GetEnergy();
922  Double_t gEres, gdE;
923  line->GetPoint(i - 1, gEres, gdE);
924  line->SetPoint(i, Eres, dE);
925 
926  }
927  }
928  //printf("sort de boucle points");
929 
930 
931 }
932 
933 
934 
943 
945 {
946  // Returns a comma-separated list of the labels of the detectors used to determine
947  // the "x" or "y" coordinates of the identification grid(s)
948  //
949  // If there is more than 1 grid and the list is not the same for all grids,
950  // prints a warning message.
951  //
952  // \param[in] axis name of grid axis i.e. "x", "X", "y" or "Y" (case insensitive)
953 
954  KVString _axis = axis;
955  _axis.ToUpper();
956 
957  if (_axis != "X" && _axis != "Y") {
958  Error("GetDetectorLabelsForGridCoord", "Called with illegal axis name '%s'", axis.Data());
959  return "";
960  }
961  KVString lab_list;
962 
963  for (auto& __gc : fGraphCoords) {
964  KVString labs;
965  if (_axis == "X") labs = __gc.second.fDetLabelsX;
966  else labs = __gc.second.fDetLabelsY;
967  if (lab_list.Length() && lab_list != labs) {
968  Error("GetDetectorLabelsForGridCoord", "Grids for telescope %s use different detector types for %s-coordinate: "
969  "%s and %s", GetName(), _axis.Data(), lab_list.Data(), labs.Data());
970  return "";
971  }
972  lab_list = labs;
973  }
974  return lab_list;
975 }
976 
977 
978 
979 
983 
985 {
986  //Return the first in the list of identification grids used by this telescope
987  //(this is for backwards compatibility with ID telescopes which had only one grid).
988  return (KVIDGraph*)GetListOfIDGrids()->First();
989 }
990 
991 
992 
993 
996 
998 {
999  //Return pointer to grid using position in list. First grid has index = 1.
1000  if (index < 1) {
1001  Error("GetIDGrid(int)", "Index must be >=1!");
1002  return nullptr;
1003  }
1004  return (KVIDGraph*)GetListOfIDGrids()->At(index - 1);
1005 }
1006 
1007 
1008 
1009 
1013 
1015 {
1016  //Return pointer to grid using "label" to search in list of grids associated
1017  //to this telescope.
1018  return (KVIDGraph*)GetListOfIDGrids()->FindObjectByLabel(label);
1019 }
1020 
1021 
1022 
1023 
1025 
1027 {
1028  AbstractMethod("GetIDMapX");
1029  return -1.;
1030 }
1031 
1032 
1033 
1038 
1040 {
1041  // Returns the pedestal associated with the 2nd detector of the telescope,
1042  // optionally depending on the given option string.
1043  // By default this returns 0, and should be overridden in specific implementations.
1044 
1045  return 0.;
1046 }
1047 
1048 
1049 
1054 
1056 {
1057  // Returns the pedestal associated with the 1st detector of the telescope,
1058  // optionally depending on the given option string.
1059  // By default this returns 0, and should be overridden in specific implementations.
1060 
1061  return 0.;
1062 }
1063 
1064 
1065 
1069 
1071 {
1072  // Return value of X coordinate to be used with the given ID grid
1073  // This corresponds to whatever was given as parameter "VARX" for the grid
1074 
1075  KVDetectorSignal* ds = fGraphCoords[g].fVarX;
1076  if (ds) return ds->GetValue();
1077  return -1;
1078 }
1079 
1080 
1081 
1085 
1087 {
1088  // Return value of Y coordinate to be used with the given ID grid
1089  // This corresponds to whatever was given as parameter "VARY" for the grid
1090 
1091  KVDetectorSignal* ds = fGraphCoords[g].fVarY;
1092  if (ds) return ds->GetValue();
1093  return -1;
1094 }
1095 
1096 
1097 
1098 
1100 
1102 {
1103  AbstractMethod("GetIDMapY");
1104  return -1.;
1105 }
1106 
1107 
1108 
1109 
1113 
1115 {
1116  //Remove all identification grids for this ID telescope
1117  //Grids are not deleted as this is handled by gIDGridManager
1118  fIDGrids.Clear();
1119  fGraphCoords.clear();
1120 }
1121 
1122 
1123 
1124 
1143 
1145 {
1146  //Static function which will create an instance of the KVIDTelescope-derived class
1147  //corresponding to 'name'
1148  //These are defined as 'Plugin' objects in the file $KVROOT/KVFiles/.kvrootrc :
1149  //~~~~~~~
1150  // # The KVMultiDetArray::GetIDTelescopes(KVDetector*de, KVDetector*e) method uses these plugins to
1151  // # create KVIDTelescope instances adapted to the specific array geometry and detector types.
1152  // # For each pair of detectors we look for a plugin with one of the following names:
1153  // # [name_of_dataset].de_detector_type[de detector thickness]-e_detector_type[de detector thickness]
1154  // # Each characteristic in [] brackets may or may not be present in the name; first we test for names
1155  // # with these characteristics, then all combinations where one or other of the characteristics is not present.
1156  // # In addition, we first test all combinations which begin with [name_of_dataset].
1157  // # The first plugin found in this order will be used.
1158  // # In addition, if for one of the two detectors there is a plugin called
1159  // # [name_of_dataset].de_detector_type[de detector thickness]
1160  // # [name_of_dataset].e_detector_type[e detector thickness]
1161  // # then we add also an instance of this 1-detector identification telescope.
1162  //~~~~~~~
1163 
1164  TPluginHandler* ph;
1165  //check and load plugin library
1166  if (!(ph = LoadPlugin("KVIDTelescope", uri)))
1167  return 0;
1168 
1169  //execute constructor/macro for identification telescope
1170  KVIDTelescope* mda = (KVIDTelescope*) ph->ExecPlugin(0);
1171  if (mda) {
1172  //set label of telescope with URI used to find plugin (minus dataset name)
1173  mda->SetLabelFromURI(uri);
1174  }
1175 
1176  return mda;
1177 }
1178 
1179 
1180 
1181 
1185 
1187 {
1188  //PRIVATE METHOD
1189  //Sets label of telescope based on URI of plugin describing child class for this telescope
1190 
1191  TString _uri(uri);
1192  if (gDataSet && _uri.BeginsWith(gDataSet->GetName())) _uri.Remove(0, strlen(gDataSet->GetName()) + 1);
1193  SetLabel(_uri.Data());
1194 }
1195 
1196 
1197 
1198 
1215 
1217 {
1218  // Initialise the identification parameters (grids, etc.) of ALL identification telescopes of this
1219  // kind (label) in the multidetector array. Therefore this method need only be called once, and not
1220  // called for each telescope. The kind/label (returned by GetLabel) of the telescope corresponds
1221  // to the URI used to find the plugin class in $KVROOT/KVFiles/.kvrootrc.
1222  // By default, this method looks for the file with name given by the environment variable
1223  //
1224  // [dataset name].IdentificationParameterList.[telescope label]: [filename]
1225  //
1226  // which is assumed to contain the list of files containing the identification grids.
1227  //
1228  // If not such envionment variable is found, the method looks for another one:
1229  //
1230  // [dataset name].IdentificationParameterFile.[telescope label]: [filename]
1231  //
1232  // which is assumed to contain identification grids.
1233 
1234  TString filename = gDataSet->GetDataSetEnv(Form("IdentificationParameterList.%s", GetLabel()));
1235 
1236  if (filename != "") {
1237  ReadIdentificationParameterFiles(filename.Data(), multidet);
1238  }
1239  else {
1240  filename = gDataSet->GetDataSetEnv(Form("IdentificationParameterFile.%s", GetLabel()));
1241 
1242  if (filename == "") {
1243  Warning("SetIdentificationParameters",
1244  "No filename defined. Should be given by %s.IdentificationParameterFile.%s or %s.IdentificationParameterFile.%s",
1246  return kFALSE;
1247  }
1248 
1249  LoadIdentificationParameters(filename, multidet);
1250  }
1251  return kTRUE;
1252 }
1253 
1254 
1255 
1256 
1264 
1266 {
1267  // In the case where the identification grids are stored in several files, this method parse
1268  // the file found with the following environment variable:
1269  //
1270  // [dataset name].IdentificationParameterList.[telescope label]: [filename]
1271  //
1272  // which contains the list of files containing the identification grids.
1273 
1274  KVFileReader fr;
1276 
1277  while (fr.IsOK()) {
1278  fr.ReadLine(0);
1279 
1280  if (fr.GetCurrentLine() != "") LoadIdentificationParameters(fr.GetCurrentLine().Data(), multidet);
1281  }
1282 
1283  fr.CloseFile();
1284 }
1285 
1286 
1287 
1288 
1291 
1293 {
1294  // This method add to the gIDGridManager list the identification grids.
1295 
1296  TString path;
1297 
1298  if ((path = gDataSet->GetFullPathToDataSetFile(filename)) == "") {
1299  Error("LoadIdentificationParameters",
1300  "File %s not found. Should be in %s",
1301  filename, gDataSet->GetDataSetDir());
1302  return;
1303  }
1304  //
1305  //Read grids from file
1306  Info("LoadIdentificationParameters", "Using file %s", path.Data());
1307  multidet->ReadGridsFromAsciiFile(path);
1308 }
1309 
1310 
1311 
1312 
1324 
1326 {
1327  //Remove identification parameters from telescope in such a way that they
1328  //can subsequently be reset e.g. with a new version.
1329  //This is used by KVMultiDetArray::UpdateIdentifications.
1330  //Child classes with specific SetIdentificationParameters methods should
1331  //also redefine this method in order to remove (destroy) cleanly the objects
1332  //created in SetIdentificationParameters.
1333  //
1334  //This default method takes the list of grids associated to the telescope,
1335  //and for each one: 1) checks if it is still in the gIDGridManager's list
1336  //2) if yes, delete the grid and remove it from gIDGridManager
1337 
1338  TIter next_grid(GetListOfIDGrids());
1339  KVIDGrid* grid;
1340  while ((grid = (KVIDGrid*)next_grid())) {
1341 
1342  if (gIDGridManager->GetGrids()->FindObject(grid)) { //this only works if KVIDTelescope uses TObject:IsEqual method (i.e. compares pointers)
1343 
1344  gIDGridManager->DeleteGrid(grid);
1345 
1346  }
1347  }
1348  //clear list of grids
1349  fIDGrids.Clear();
1351 }
1352 
1353 
1354 
1355 
1374 
1376 {
1377  // The energy of each particle is calculated as follows:
1378  //
1379  // E = dE_1 + dE_2 + ... + dE_N
1380  //
1381  // dE_1, dE_2, ... = energy losses measured in each detector through which
1382  // the particle has passed (or stopped, in the case of dE_N).
1383  // These energy losses are corrected for (Z,A)-dependent effects
1384  // such as pulse-heigth defect in silicon detectors, losses in
1385  // windows of gas detectors, etc.
1386  //
1387  // Whenever possible, the energy loss for fired detectors which are uncalibrated
1388  // or not functioning is calculated. In this case the status returned by GetCalibStatus()
1389  // will be KVIDTelescope::kCalibStatus_Calculated.
1390  // If none of the detectors is calibrated, the particle's energy cannot be calculated &
1391  // the status will be KVIDTelescope::kCalibStatus_NoCalibrations.
1392  // Otherwise, the status code will be KVIDTelescope::kCalibStatus_OK.
1393 
1394  //status code
1396 
1397  UInt_t z = nuc->GetZ();
1398  //uncharged particles
1399  if (z == 0) return;
1400 
1401  KVDetector* d1 = GetDetector(1);
1402  KVDetector* d2 = (GetSize() > 1 ? GetDetector(2) : 0);
1403  Bool_t d1_cal = d1->IsCalibrated();
1404  Bool_t d2_cal = (d2 ? d2->IsCalibrated() : kFALSE);
1405 
1406  //no calibrations
1407  if (!d1_cal && !d2)
1408  return;
1409  if ((d1 && d2) && !d1_cal && !d2_cal)
1410  return;
1411 
1412  //status code
1414 
1415  UInt_t a = nuc->GetA();
1416 
1417  // particles stopped in first member of telescope
1418  if (nuc->GetStatus() == 3) {
1419  if (d1_cal) {
1420  nuc->SetEnergy(d1->GetCorrectedEnergy(nuc, -1, kFALSE)); //N.B.: transmission=kFALSE because particle stop in d1
1421  }
1422  return;
1423  }
1424 
1425  Double_t e1, e2, einc;
1426  e1 = e2 = einc = 0.0;
1427 
1428  if (!d1_cal) {//1st detector not calibrated - calculate from residual energy in 2nd detector
1429 
1430  //second detector must exist and have all acquisition parameters fired with above-pedestal value
1431  if (d2 && d2->Fired("Pall")) e2 = d2->GetCorrectedEnergy(nuc, -1, kFALSE); //N.B.: transmission=kFALSE because particle stop in d2
1432  if (e2 <= 0.0) {
1433  // zero energy loss in 2nd detector ? can't do anything...
1435  return;
1436  }
1437  //calculate & set energy loss in dE detector
1438  //N.B. using e2 for the residual energy after detector 1 means
1439  //that we are assuming the particle stops in detector 2.
1440  //if this is not true, we will underestimate the energy of the particle.
1441  e1 = d1->GetDeltaEFromERes(z, a, e2);
1442  if (e1 < 0.0) e1 = 0.0;
1443  else {
1444  d1->SetEnergyLoss(e1);
1445  d1->SetEResAfterDetector(e2);
1446  e1 = d1->GetCorrectedEnergy(nuc);
1447  //status code
1449  }
1450  }
1451  else { //1st detector is calibrated too: get corrected energy loss
1452 
1453  e1 = d1->GetCorrectedEnergy(nuc);
1454 
1455  }
1456 
1457  if (d2 && !d2_cal) {//2nd detector not calibrated - calculate from energy loss in 1st detector
1458 
1459  e1 = d1->GetCorrectedEnergy(nuc);
1460  if (e1 <= 0.0) {
1461  // zero energy loss in 1st detector ? can't do anything...
1463  return;
1464  }
1465  //calculate & set energy loss in 2nd detector
1466  e2 = d1->GetEResFromDeltaE(z, a);
1467  if (e2 < 0.0) e2 = 0.0;
1468  else {
1469  e2 = d2->GetDeltaE(z, a, e2);
1470  d2->SetEnergyLoss(e2);
1471  e2 = d2->GetCorrectedEnergy(nuc);
1472  //status code
1474  }
1475  }
1476  else if (d2) { //2nd detector is calibrated too: get corrected energy loss
1477 
1478  e2 = d2->GetCorrectedEnergy(nuc, -1, kFALSE);//N.B.: transmission=kFALSE because particle assumed to stop in d2
1479  // recalculate corrected energy in first stage using info on Eres
1480  d1->SetEResAfterDetector(e2);
1481  e1 = d1->GetCorrectedEnergy(nuc);
1482  }
1483 
1484  //incident energy of particle (before 1st member of telescope)
1485  einc = e1 + e2;
1486 
1487  Double_t coherence_tolerance = gEnv->GetValue("KVIDTelescope.CoherencyTolerance", 1.05);
1488  if (coherence_tolerance < 1) coherence_tolerance += 1.00;
1489 
1490  //Now we have to work our way up the list of detectors from which the particle was
1491  //reconstructed. For each fired & calibrated detector which is only associated with
1492  //one particle in the events, we add the corrected measured energy loss
1493  //to the particle. For uncalibrated, unfired detectors and detectors through which
1494  //more than one particle has passed, we calculate the corrected energy loss and add it
1495  //to the particle.
1496  int ndets = nuc->GetNumDet();
1497  if (ndets > (int)GetSize()) { //particle passed through other detectors before this idtelesocpe
1498  //look at detectors not in this id telescope
1499  int idet = GetSize();//next detector after delta-e member of IDTelescope (stopping detector = 0)
1500  while (idet < ndets) {
1501 
1502  KVDetector* det = nuc->GetDetector(idet);
1503  if (det->Fired() && det->IsCalibrated() && det->GetNHits() == 1) {
1504  Double_t dE = det->GetEnergy();
1505  //in order to check if particle was really the only one to
1506  //hit each detector, we calculate the particle's energy loss
1507  //from its residual energy. if the measured energy loss is
1508  //significantly larger, there may be a second particle.
1509  e1 = det->GetDeltaEFromERes(z, a, einc);
1510  if (e1 < 0.0) e1 = 0.0;
1511  det->SetEResAfterDetector(einc);
1512  dE = det->GetCorrectedEnergy(nuc);
1513  einc += dE;
1514  }
1515  else {
1516  // Uncalibrated/unfired/multihit detector. Calculate energy loss.
1517  //calculate energy of particle before detector from energy after detector
1518  e1 = det->GetDeltaEFromERes(z, a, einc);
1519  if (e1 < 0.0) e1 = 0.0;
1520  if (det->GetNHits() > 1) {
1521  //Info("CalculateParticleEnergy",
1522  // "Detector %s was hit by %d particles. Calculated energy loss for particle %f MeV",
1523  // det->GetName(), det->GetNHits(), e1);
1524  if (!(det->Fired() && det->IsCalibrated())) {
1525  det->SetEnergyLoss(e1 + det->GetEnergy());// sum up calculated energy losses in uncalibrated detector
1526  }
1527  //status code
1529  }
1530  else if (!det->Fired() || !det->IsCalibrated()) {
1531  //Info("CalculateParticleEnergy",
1532  // "Detector %s uncalibrated/not fired. Calculated energy loss for particle %f MeV",
1533  // det->GetName(), e1);
1534  det->SetEnergyLoss(e1);
1535  //status code
1537  }
1538  det->SetEResAfterDetector(einc);
1539  e1 = det->GetCorrectedEnergy(nuc, e1);
1540  einc += e1;
1541  }
1542  idet++;
1543  }
1544  }
1545  //einc is now the energy of the particle before crossing the first detector
1546  nuc->SetEnergy(einc);
1547 }
1548 
1549 
1550 
1551 
1561 
1563 {
1564  // Returns name of default ID grid class for this ID telescope.
1565  // This is defined in a .kvrootrc configuration file by one of the following:
1566  // KVIDTelescope.DefaultGrid:
1567  // KVIDTelescope.DefaultGrid.[type]:
1568  // where [type] is the type of this identification telescope (which is given
1569  // by the character string returned by method GetLabel()... sorry :( )
1570  // If no default grid is defined for the specific type of this telescope,
1571  // the default defined by KVIDTelescope.DefaultGrid is used.
1572 
1573  TString specific;
1574  specific.Form("KVIDTelescope.DefaultGrid.%s", GetLabel());
1575  return gEnv->GetValue(specific.Data(), gEnv->GetValue("KVIDTelescope.DefaultGrid", "KVIDGraph"));
1576 }
1577 
1578 
1579 
1585 
1587 {
1588  //Create a dE-E grid (energy loss in detector 1 versus residual energy in detector 2) for a given list of isotopes
1589  // - AperZ : list of A for each Z. (ex: "1=1-3,2=3-4 5,3=6-8,4=7 9-12"...)
1590  // - npoints : number of points in each generated line
1591  // - xfactor : scales the detector 2 thickness to prolongate the lines
1592 
1593  if (GetSize() <= 1) return 0;
1594  if (!GetDetector(1) || !GetDetector(2)) return 0;
1595  double thickness = GetDetector(2)->GetThickness();
1596  GetDetector(2)->SetThickness(thickness * xfactor);
1597 
1598  Info("CalculateDeltaE_EGrid", "called with KVNameValueList");
1599 
1600  KVIDGrid* idgrid = newGrid(0);
1601 
1602  for (auto par : AperZ) {
1603  KVString tmp = par.GetName();
1604  int zz = tmp.Atoi();
1605  KVNumberList alist = par.GetString();
1606  for (auto aa : alist) {
1607  addLineToGrid(idgrid, zz, aa, npoints);
1608  }
1609  }
1610 
1611  GetDetector(2)->SetThickness(thickness);
1612 
1613  return idgrid;
1614 
1615 }
1616 
1617 
1618 
1626 
1627 KVIDGrid* KVIDTelescope::CalculateDeltaE_EGrid(const KVNumberList& Zrange, Int_t deltaA, Int_t npoints, Double_t lifetime, UChar_t massformula, Double_t xfactor)
1628 {
1629  //Create a dE-E grid (energy loss in detector 1 versus residual energy in detector 2) for a given list of isotopes
1630  // - Zrange : list of element for which we create lines
1631  // - deltaA : number of isotopes generated for each Z around massformula (ex: deltaA=1, Aref-1 Aref Aref+1)
1632  // - npoints : number of points in each generated line
1633  // - lifetime: remove isotopes with lifetime lower than this value
1634  // - xfactor : scales the detector 2 thickness to prolongate the lines
1635 
1636  if (GetSize() <= 1) return 0;
1637  if (!GetDetector(1) || !GetDetector(2)) return 0;
1638  double thickness = GetDetector(2)->GetThickness();
1639  GetDetector(2)->SetThickness(thickness * xfactor);
1640 
1641  Info("CalculateDeltaE_EGrid", "called with KVNumberList");
1642 
1643  KVIDGrid* idgrid = newGrid(0);
1644  KVNucleus part;
1645 
1646  Zrange.Begin();
1647  while (!Zrange.End()) {
1648  Int_t zz = Zrange.Next();
1649  part.SetZ(zz, massformula);
1650  Int_t aref = part.GetA();
1651  for (Int_t aa = aref - deltaA; aa <= aref + deltaA; aa += 1) {
1652  part.SetA(aa);
1653  if (part.IsKnown() && (part.GetLifeTime() > lifetime)) {
1654  addLineToGrid(idgrid, zz, aa, npoints);
1655 
1656  }
1657  }
1658  }
1659  GetDetector(2)->SetThickness(thickness);
1660  return idgrid;
1661 }
1662 
1663 
1664 
1670 
1672 {
1673  //Create a dE-E grid (energy loss in detector 1 versus residual energy in detector 2) for a given list of isotopes
1674  // - haa_zz : lines will be generated for A,Z filled with 1 in this histogram
1675  // - Zonly : if true, generate only one line per Z with the <A>(Z) of the histogram
1676  // - npoints : number of points in each generated line
1677 
1678  if (GetSize() <= 1) return 0;
1679  if (!GetDetector(1) || !GetDetector(2)) return 0;
1680  double thickness = GetDetector(2)->GetThickness();
1681 
1682  Info("CalculateDeltaE_EGrid", "called with TH2");
1683 
1684  KVIDGrid* idgrid = newGrid(0);
1685  KVNucleus part;
1686 
1687  for (Int_t nx = 1; nx <= haa_zz->GetNbinsX(); nx += 1) {
1688 
1689  Int_t zz = TMath::Nint(haa_zz->GetXaxis()->GetBinCenter(nx));
1690  KVNumberList nlA;
1691  Double_t sumA = 0, sum = 0;
1692  for (Int_t ny = 1; ny <= haa_zz->GetNbinsY(); ny += 1) {
1693  Double_t stat = haa_zz->GetBinContent(nx, ny);
1694  if (stat > 0) {
1695  Double_t val = haa_zz->GetYaxis()->GetBinCenter(ny);
1696  nlA.Add(TMath::Nint(val));
1697  sumA += val * stat;
1698  sum += stat;
1699  }
1700  }
1701  sumA /= sum;
1702  Int_t nA = nlA.GetNValues();
1703  if (nA == 0) {
1704  Warning("CalculateDeltaE_EGrid", "no count for Z=%d", zz);
1705  }
1706  else {
1707  if (Zonly) {
1708  nlA.Clear();
1709  nlA.Add(TMath::Nint(sumA));
1710  }
1711  else {
1712  if (nA == 1) {
1713  Int_t aref = nlA.Last();
1714  nlA.Add(aref - 1);
1715  nlA.Add(aref + 1);
1716  }
1717  }
1718  part.SetZ(zz);
1719  nlA.Begin();
1720  while (!nlA.End()) {
1721  Int_t aa = nlA.Next();
1722  part.SetA(aa);
1723  if (part.IsKnown()) {
1724  addLineToGrid(idgrid, zz, aa, npoints);
1725  }
1726  }
1727  }
1728  }
1729  return idgrid;
1730 }
1731 
1732 
1733 
1749 
1751 {
1752  // Returns the Y-axis value in the 2D identification map containing isotope (Z,A)
1753  // corresponding to either the given X-axis/Eres value or the current X-axis value given by GetIDGridXCoord()
1754  // If no mass information is available, just give Z.
1755  //
1756  // In this (default) implementation this means scanning the ID grids associated with
1757  // this telescope until we find an identification line Z or (Z,A), and then interpolating
1758  // the Y-coordinate for the current X-coordinate value.
1759  //
1760  // Status variable can take one of following values:
1761  //
1762  // KVIDTelescope::kMeanDE_OK all OK
1763  // KVIDTelescope::kMeanDE_XtooSmall X-coordinate is smaller than smallest X-coordinate of ID line
1764  // KVIDTelescope::kMeanDE_XtooLarge X-coordinate is larger than largest X-coordinate of ID line
1765  // KVIDTelescope::kMeanDE_NoIdentifie No identifier found for Z or (Z,A)
1766 
1767  status = kMeanDE_OK;
1768  // loop over grids
1769  TIter next(GetListOfIDGrids());
1770  KVIDGrid* grid;
1771  KVIDLine* idline = 0;
1772  while ((grid = (KVIDGrid*)next())) {
1773  idline = (KVIDLine*)grid->GetIdentifier(Z, A);
1774  if (idline) break;
1775  }
1776  if (!idline) {
1777  status = kMeanDE_NoIdentifier;
1778  return -1.;
1779  }
1780  Double_t x, x1, y1, x2, y2;
1781  x = (Eres < 0 ? GetIDGridXCoord(grid) : Eres);
1782  idline->GetEndPoint(x2, y2);
1783  if (x > x2) {
1784  status = kMeanDE_XtooLarge;
1785  return -1;
1786  }
1787  idline->GetStartPoint(x1, y1);
1788  if (x < x1) {
1789  status = kMeanDE_XtooSmall;
1790  return -1.;
1791  }
1792  return idline->Eval(x);
1793 }
1794 
1795 
1796 
1797 
1806 
1808 {
1809  // Return kTRUE if energy of ION is > minimum incident energy required for identification
1810  // This theoretical limit is defined here to be the incident energy for which the
1811  // dE in the first detector of a dE-E telescope is maximum.
1812  // If EINC>0 it is assumed to be the energy of the ion just before the first detector
1813  // (case where ion would have to pass other detectors before reaching this telescope).
1814  //
1815  // If this is not a dE-E telescope, we return kTRUE by default.
1816 
1817  if (GetSize() < 2) return kTRUE;
1818 
1819  KVDetector* dEdet = GetDetector(1);
1820  Double_t emin = dEdet->GetEIncOfMaxDeltaE(ION->GetZ(), ION->GetA());
1821  if (EINC > 0.0) return (EINC > emin);
1822  return (ION->GetEnergy() > emin);
1823 }
1824 
1825 
1826 
1853 
1855 {
1856  // For filtering simulations
1857  // Set the n->IsZMeasured() and n->IsAMeasured() status of the particle
1858  // In principle this depends on whether this telescope provides mass
1859  // identification or not, but this may depend on the particle's energy.
1860  // If A was not measured, it will be replaced with a value calculated
1861  // from whatever mass formula is used for the particle.
1862  //
1863  // In order to enable mass identification for certain telescopes without a dedicated
1864  // implementation (e.g. for simulating array response), put the following lines
1865  // in your .kvrootrc:
1866  //
1867  // [dataset].[telescope label].MassID: yes
1868  //
1869  // If you want to limit mass identification to certain values of Z and/or A,
1870  // add the following line:
1871  //
1872  // [dataset].[telescope label].MassID.Validity: [expression]
1873  //
1874  // where [expression] is some valid C++ boolean expression involving Z and/or A,
1875  // for example
1876  //
1877  // [dataset].[telescope label].MassID.Validity: (Z>3)&&(A<20)
1878  //
1879  // Then this expression will be tested here in order to determine particle
1880  // identification status
1881 
1882  n->SetZMeasured();
1883  if (!HasMassID()) {
1884  n->SetAMeasured(kFALSE);
1885  // beware - changing particle's mass changes its KE (momentum is conserved)
1886  double e = n->GetE();
1887  n->SetZ(n->GetZ());// use mass formula for A
1888  n->SetE(e);
1889  }
1890  else {
1891  if (fMassIDValidity) n->SetAMeasured(fMassIDValidity->Test(n)); // test expression for mass ID validity
1892  else n->SetAMeasured(); // no expression set; all nuclei are identified in mass
1893  if (!n->IsAMeasured()) {
1894  // beware - changing particle's mass changes its KE (momentum is conserved)
1895  double e = n->GetE();
1896  n->SetZ(n->GetZ()); // use mass formula for A
1897  n->SetE(e);
1898  }
1899  }
1900 }
1901 
1902 
1903 
1906 
1908 {
1909  // Open IdentificationBilan.dat file with given path
1910 
1912  fgIdentificationBilan = new TEnv(path);
1913 }
1914 
1915 
1916 
1919 
1921 {
1922  // Set status of ID Telescope for given system
1923  if (!(fgIdentificationBilan->GetValue(Form("%s.%s", system.Data(), GetName()), kTRUE))) ResetBit(kReadyForID);
1924 }
1925 
1926 
int Int_t
unsigned int UInt_t
KVDataSet * gDataSet
Definition: KVDataSet.cpp:29
KVIDGridManager * gIDGridManager
KVTemplateParticleCondition< KVNucleus > KVParticleCondition
ClassImp(KVPartitionList) void KVPartitionList
Initialisation.
#define d(i)
#define e(i)
unsigned char UChar_t
char Char_t
const Bool_t kFALSE
bool Bool_t
double Double_t
const Bool_t kTRUE
const char Option_t
R__EXTERN TEnv * gEnv
char * Form(const char *fmt,...)
void SetLabel(const Char_t *lab)
Definition: KVBase.h:194
virtual const Char_t * GetType() const
Definition: KVBase.h:176
const Char_t * GetLabel() const
Definition: KVBase.h:198
static TPluginHandler * LoadPlugin(const Char_t *base, const Char_t *uri="0")
Definition: KVBase.cpp:793
UInt_t GetNumber() const
Definition: KVBase.h:219
const Char_t * GetDataSetDir() const
Definition: KVDataSet.cpp:720
Bool_t HasCalibIdentInfos() const
Definition: KVDataSet.h:403
const Char_t * GetDataSetEnv(const Char_t *type, const Char_t *defval="") const
Definition: KVDataSet.cpp:758
TString GetFullPathToDataSetFile(const Char_t *filename)
Definition: KVDataSet.cpp:1887
Signal output from a mathematical combination of other signals.
Base class for output signal data produced by a detector.
virtual Bool_t IsValid() const
virtual Double_t GetValue(const KVNameValueList &params="") const
Base class for detector geometry description.
Definition: KVDetector.h:159
virtual Bool_t IsOK() const
Definition: KVDetector.h:681
virtual Double_t GetMaxDeltaE(Int_t Z, Int_t A)
virtual Double_t GetERes(Int_t Z, Int_t A, Double_t Einc)
virtual Double_t GetELostByParticle(KVNucleus *, TVector3 *norm=0)
Definition: KVDetector.cpp:276
virtual void SetEnergyLoss(Double_t e) const
Definition: KVDetector.h:372
virtual Double_t GetEnergy() const
Definition: KVDetector.h:348
Int_t GetNHits() const
Return the number of particles hitting this detector in an event.
Definition: KVDetector.h:435
virtual Double_t GetEIncOfMaxDeltaE(Int_t Z, Int_t A)
virtual void Clear(Option_t *opt="")
Definition: KVDetector.cpp:596
void SetThickness(Double_t thick)
virtual Bool_t Fired(Option_t *opt="any") const
Definition: KVDetector.h:450
virtual KVDetectorSignal * GetDetectorSignal(const KVString &type) const
Definition: KVDetector.h:532
Bool_t IsCalibrated() const
Definition: KVDetector.h:389
Bool_t AddDetectorSignalExpression(const KVString &type, const KVString &_expr)
virtual void SetEResAfterDetector(Double_t e)
Definition: KVDetector.h:629
virtual TList * GetAlignedDetectors(UInt_t direction=1)
virtual Double_t GetDeltaE(Int_t Z, Int_t A, Double_t Einc)
virtual Double_t GetDeltaEFromERes(Int_t Z, Int_t A, Double_t Eres)
virtual void DetectParticle(KVNucleus *, TVector3 *norm=0)
Definition: KVDetector.cpp:203
virtual void Print(Option_t *option="") const
Definition: KVDetector.cpp:364
virtual Double_t GetCorrectedEnergy(KVNucleus *, Double_t e=-1., Bool_t transmission=kTRUE)
Definition: KVDetector.cpp:810
Handle reading columns of numeric data in text files.
Definition: KVFileReader.h:119
KVString GetCurrentLine()
Definition: KVFileReader.h:319
void CloseFile()
Definition: KVFileReader.h:236
ReadStatus ReadLine(const KVString &pattern="")
Definition: KVFileReader.h:242
Bool_t IsOK()
Definition: KVFileReader.h:230
Bool_t OpenFileToRead(const KVString &filename)
Definition: KVFileReader.h:209
Group of detectors which can be treated independently of all others in array.
Definition: KVGroup.h:19
@ kForwards
Definition: KVGroup.h:32
Line in ID grid used to delimit regions where no identification is possible.
Definition: KVIDCutLine.h:22
virtual void SetName(const char *name)
This is redeclared to make it appear in context menus for KVIDCutLines.
Definition: KVIDCutLine.h:83
virtual void SetAcceptedDirection(const Char_t *dir)
Definition: KVIDCutLine.cpp:74
Base class for particle identification in a 2D map.
Definition: KVIDGraph.h:31
void Add(TString, KVIDentifier *)
Definition: KVIDGraph.cpp:838
void SetRunList(const char *runlist)
Definition: KVIDGraph.h:155
virtual void SetName(const char *name)
Definition: KVIDGraph.h:139
void AddIDTelescope(KVBase *t)
Definition: KVIDGraph.h:406
KVIDentifier * GetCut(const Char_t *name) const
Definition: KVIDGraph.h:279
virtual void Identify(Double_t, Double_t, KVIdentificationResult *) const =0
KVIDentifier * GetIdentifier(Int_t Z, Int_t A) const
Definition: KVIDGraph.cpp:310
virtual Bool_t IsIdentifiable(Double_t, Double_t, TString *rejected_by=nullptr) const
Definition: KVIDGraph.cpp:1269
const Char_t * GetName() const
Definition: KVIDGraph.cpp:1332
virtual void SetOnlyZId(Bool_t yes=kTRUE)
Definition: KVIDGraph.cpp:1496
KVList * GetGrids()
void DeleteGrid(KVIDGraph *, Bool_t update=kTRUE)
Abstract base class for 2D identification grids in e.g. (dE,E) maps.
Definition: KVIDGrid.h:73
Base class for lines/cuts used for particle identification in 2D data maps.
Definition: KVIDLine.h:142
void GetStartPoint(Double_t &x, Double_t &y) const
void GetEndPoint(Double_t &x, Double_t &y) const
Base class for all detectors or associations of detectors in array which can identify charged particl...
Definition: KVIDTelescope.h:83
KVIDGrid * newGrid(bool onlyZ)
void LoadIdentificationParameters(const Char_t *filename, const KVMultiDetArray *multidet)
This method add to the gIDGridManager list the identification grids.
virtual Double_t GetIDMapY(Option_t *opt="")
void init()
default init
KVIDGrid * CalculateDeltaE_EGrid(const KVNameValueList &AperZ, Int_t npoints=30, Double_t xfactor=1.)
void SetGroup(KVGroup *kvg)
KVList fMultiDetExpressions
used to clean up any multi-detector signal expressions generated to calculate X/Y coordinates
void SetLabelFromURI(const Char_t *uri)
@ kCalibStatus_NoCalibrations
Double_t GetIDGridYCoord(KVIDGraph *) const
virtual Bool_t Identify(KVIdentificationResult *, Double_t x=-1., Double_t y=-1.)
virtual Double_t GetIDMapX(Option_t *opt="")
static KVIDTelescope * MakeIDTelescope(const Char_t *name)
virtual Double_t GetPedestalY(Option_t *opt="")
virtual Bool_t SetIdentificationParameters(const KVMultiDetArray *)
void SetIDGrid(KVIDGraph *)
Bool_t HasMassID() const
KVDetector * GetDetector(UInt_t n) const
void ReadIdentificationParameterFiles(const Char_t *filename, const KVMultiDetArray *multidet)
virtual Bool_t CheckTheoreticalIdentificationThreshold(KVNucleus *, Double_t=0.0)
UInt_t GetGroupNumber()
void addLineToGrid(KVIDGrid *gg, int zz, int aa, int npoints)
KVGroup * GetGroup() const
Int_t fCalibStatus
temporary variable holding status code for last call to Calibrate(KVReconstructedNucleus*)
Definition: KVIDTelescope.h:97
virtual Double_t GetPedestalX(Option_t *opt="")
virtual void CalculateParticleEnergy(KVReconstructedNucleus *nuc)
virtual void AddDetector(KVDetector *d)
virtual Double_t GetMeanDEFromID(Int_t &status, Int_t Z, Int_t A=-1, Double_t Eres=-1.0)
const KVList * GetDetectors() const
KVDetectorSignal * GetSignalFromGridVar(const KVString &var, const KVString &axe, KVString &det_labels)
KVIDGraph * GetIDGrid()
virtual void Print(Option_t *opt="") const
std::unordered_map< KVIDGraph *, GraphCoords > fGraphCoords
X/Y coordinates from detector signals for ID maps.
virtual UShort_t GetIDCode()
void CheckIdentificationBilan(const TString &system)
Set status of ID Telescope for given system.
virtual void RemoveIdentificationParameters()
static void OpenIdentificationBilan(const TString &path)
Open IdentificationBilan.dat file with given path.
void SetHasMassID(Bool_t yes=kTRUE)
virtual void Initialize(void)
Double_t GetIDGridXCoord(KVIDGraph *) const
const Char_t * GetDefaultIDGridClass()
virtual void SetIdentificationStatus(KVReconstructedNucleus *)
void GetIDGridCoords(Double_t &X, Double_t &Y, KVIDGraph *grid, Double_t x=-1, Double_t y=-1)
UInt_t GetSize() const
KVUnownedList fIDGrids
identification grid(s)
Definition: KVIDTelescope.h:92
KVUnownedList fDetectors
list of detectors in telescope
Definition: KVIDTelescope.h:90
KVString GetDetectorLabelsForGridCoord(const KVString &axis) const
const KVList * GetListOfIDGrids() const
std::unique_ptr< KVParticleCondition > fMassIDValidity
may be used to limit mass identification to certain Z and/or A range
KVGroup * fGroup
group to which telescope belongs
Definition: KVIDTelescope.h:91
static TEnv * fgIdentificationBilan
Definition: KVIDTelescope.h:86
virtual void RemoveGrids()
virtual TGraph * MakeIDLine(KVNucleus *nuc, Double_t Emin, Double_t Emax, Double_t Estep=0.0)
Base class for identification ridge lines corresponding to different nuclear species.
Definition: KVIDZALine.h:32
Full result of one attempted particle identification.
Bool_t IDattempted
=kTRUE if identification was attempted
Bool_t IDOK
general quality of identification, =kTRUE if acceptable identification made
void SetGridName(const Char_t *n)
void Clear(Option_t *opt="")
Reset to initial values.
TString Rejecting_Cut
name of cut in grid which rejected particle for identification
Int_t IDquality
specific quality code returned by identification procedure
Int_t IDcode
a general identification code for this type of identification
void SetIDType(const Char_t *t)
virtual Double_t GetThickness() const
Definition: KVMaterial.cpp:487
Double_t GetEmaxValid(Int_t Z, Int_t A)
virtual Double_t GetEResFromDeltaE(Int_t Z, Int_t A, Double_t dE=-1.0, enum SolType type=kEmax)
Base class for describing the geometry of a detector array.
Bool_t ReadGridsFromAsciiFile(const Char_t *) const
Handles lists of named parameters with different types, a list of KVNamedParameter objects.
Description of properties and kinematics of atomic nuclei.
Definition: KVNucleus.h:125
Bool_t IsKnown(int z=-1, int a=-1) const
Definition: KVNucleus.cpp:1279
Int_t GetA() const
Definition: KVNucleus.cpp:799
void SetA(Int_t a)
Definition: KVNucleus.cpp:655
void SetZ(Int_t z, Char_t mt=-1)
Definition: KVNucleus.cpp:704
Int_t GetZ() const
Return the number of proton / atomic number.
Definition: KVNucleus.cpp:770
Double_t GetLifeTime(Int_t z=-1, Int_t a=-1) const
Definition: KVNucleus.cpp:1040
Strings used to represent a set of ranges of values.
Definition: KVNumberList.h:83
Bool_t End(void) const
Definition: KVNumberList.h:197
Int_t GetNValues() const
void Begin(void) const
void Add(Int_t)
Add value 'n' to the list.
void Clear(Option_t *="")
Empty number list, reset it to initial state.
Int_t Last() const
Returns largest number included in list.
Int_t Next(void) const
Double_t GetEnergy() const
Definition: KVParticle.h:623
void SetEnergy(Double_t e)
Definition: KVParticle.h:601
Nuclei reconstructed from data measured by a detector array .
KVDetector * GetDetector(const TString &label) const
virtual TObject * FindObjectByLabel(const Char_t *) const
virtual void Clear(Option_t *option="")
virtual TObject * At(Int_t idx) const
virtual TObject * First() const
virtual void SetCleanup(Bool_t enable=kTRUE)
virtual void Add(TObject *obj)
virtual TObject * FindObject(const char *name) const
Extension of ROOT TString class which allows backwards compatibility with ROOT v3....
Definition: KVString.h:72
void Begin(TString delim) const
Definition: KVString.cpp:565
Bool_t End() const
Definition: KVString.cpp:634
KVString Next(Bool_t strip_whitespace=kFALSE) const
Definition: KVString.cpp:695
Int_t GetNValues(TString delim) const
Definition: KVString.cpp:886
Optimised list in which named objects can only be placed once.
virtual void Add(TObject *obj)
virtual void SetLineColor(Color_t lcolor)
virtual Double_t GetBinCenter(Int_t bin) const
static TClass * GetClass(Bool_t load=kTRUE, Bool_t silent=kFALSE)
const char * GetVarY() const
const char * GetVarX() const
virtual const char * GetValue(const char *name, const char *dflt) const
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Int_t GetN() const
virtual Double_t Eval(Double_t x, TSpline *spline=nullptr, Option_t *option="") const
virtual Int_t GetNbinsY() const
TAxis * GetXaxis()
TAxis * GetYaxis()
virtual Int_t GetNbinsX() const
virtual Double_t GetBinContent(Int_t bin) const
void Reset()
virtual const char * GetName() const
virtual void SetName(const char *name)
void AbstractMethod(const char *method) const
void SetBit(UInt_t f)
virtual void Warning(const char *method, const char *msgfmt,...) const
virtual void Error(const char *method, const char *msgfmt,...) const
void ResetBit(UInt_t f)
virtual void Info(const char *method, const char *msgfmt,...) const
Longptr_t ExecPlugin(int nargs, const T &... params)
Ssiz_t Length() const
Int_t Atoi() const
void ToUpper()
Bool_t BeginsWith(const char *s, ECaseCompare cmp=kExact) const
const char * Data() const
void Form(const char *fmt,...)
TString & Remove(EStripType s, char c)
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
TString & ReplaceAll(const char *s1, const char *s2)
TLine * line
Double_t y[n]
Double_t x[n]
const Int_t n
TGraphErrors * gr
const long double cl
Definition: KVUnits.h:85
const long double g
masses
Definition: KVUnits.h:72
Int_t Nint(T x)
Double_t Exp(Double_t x)
constexpr Double_t E()
Double_t Log(Double_t x)
Bool_t Even(Long_t a)
KVDetectorSignal * fVarY
KVDetectorSignal * fVarX
auto * a