KaliVeda
Toolkit for HIC analysis
KVINDRADSTConverter.cpp
1 #include "KVDataAnalyser.h"
2 #define KVDSTTreeConverter_cxx
3 
4 #include "KVINDRADSTConverter.h"
5 #include "KVUserAnalysisOptionList.h"
6 #include "KVINDRADB.h"
7 #include <KVReconNucTrajectory.h>
8 
9 
14 
16 {
17  // The Begin() function is called at the start of the query.
18  // When running with PROOF Begin() is only called on the client.
19  // The tree argument is deprecated (on PROOF 0 is passed).
20 
21  Info("Begin", "Reading %lld events from file", total_events);
22  events_read = 0;
23  KVUserAnalysisOptionList option_list;
24  option_list.ParseOptions(GetOption());
25  run_number = option_list.GetOpt("run_number").Atoi();
27  if (!gIndra)
29  else
30  gIndra->SetParameters(run_number);
31 
32  with_phoswich = gIndra->GetListOfPhoswich() && !gIndra->GetListOfPhoswich()->IsEmpty();
33  camp1 = gDataSet->IsCalled("INDRA_camp1");
34  camp2 = gDataSet->IsCalled("INDRA_camp2");
35  camp4 = gDataSet->IsCalled("INDRA_camp4");
36 
38 
39  output_file = gDataSet->NewRunfile("recon", run_index);
40  output_tree = new TTree("ReconEvents", Form("%s : %s : root events converted from DST",
41  gIndraDB->GetDBRunFile(run_index).GetName(),
42  gIndraDB->GetDBRunFile(run_index).GetTitle()));
44 }
45 
46 
47 
52 
54 {
55  // The SlaveBegin() function is called after the Begin() function.
56  // When running with PROOF SlaveBegin() is called on each slave server.
57  // The tree argument is deprecated (on PROOF 0 is passed).
58 
60 }
61 
62 
63 
65 
67 {
68  fReader.SetLocalEntry(entry);
69 
70  the_event->Clear();
71 
72  for (int i = 0; i < *npart_traite; ++i) {
73  if (idcode[i] >= KVINDRA::NO_IDENTIFICATION) continue;
74 
75  if (idcode[i] == KVINDRA::ID_GAMMA) {
76  KVDetector* CSI = nullptr;
77  if (with_phoswich && icou[i] == 1)
78  CSI = get_phos(i);
79  else
80  CSI = get_csi(i);
81  if (!CSI) {
82  Error("Process", "Got no CSI for cod=%d cou=%d mod=%d", idcode[i], icou[i], imod[i]);
83  }
84  else {
85  the_event->GetParameters()->IncrementValue("INDRA_GAMMA_MULT", 1);
86  the_event->GetParameters()->IncrementValue("INDRA_GAMMA_DETS", CSI->GetName());
87  }
88  continue;
89  }
90  KVINDRADetector* det = nullptr;
91 
92  if (stopped_in_csi(i)) det = get_csi(i);
93  else if (stopped_in_si(i)) det = get_si(i);
94  else if (stopped_in_ci(i)) det = get_ci(i);
95  else if (stopped_in_phos(i)) det = get_phos(i);
96  else if (stopped_in_si75(i)) det = get_si75(i);
97  else if (stopped_in_sili(i)) det = get_sili(i);
98  if (!det) {
99  // nothing fired? try direct approach?
100  switch (idcode[i]) {
101  case KVINDRA::ID_CSI_PSA:
103  if (!with_phoswich || icou[i] > 1)
104  det = gIndra->GetDetectorByType((UInt_t)icou[i], (UInt_t)imod[i], CsI_T);
105  else
106  det = gIndra->GetDetectorByType((UInt_t)icou[i], (UInt_t)imod[i], Phos_T);
107  break;
108 
109  default:
110  det = nullptr;
111  }
112  }
113 
114  auto idtelescope_code = idcode[i];
115  switch (idcode[i]) {
117  if (!with_phoswich && icou[i] == 1)
118  idtelescope_code = KVINDRA::ID_SI_CSI;
119  else if (icou[i] < 10)
120  idtelescope_code = KVINDRA::ID_CI_SI;
121  else if (icou[i] >= 10) {
122  if (si75_fired(i)) idtelescope_code = KVINDRA::ID_CI_SI75;
123  else idtelescope_code = KVINDRA::ID_CI_CSI;
124  }
125  break;
126 
128  if (icou[i] > 1 && icou[i] < 10) {
129  idtelescope_code = KVINDRA::ID_CI_SI;
130  }
131  else if (icou[i] > 9) {
132  idtelescope_code = KVINDRA::ID_CI_CSI;
133  det = get_csi(i);
134  }
135  break;
136 
137  case KVINDRA::ID_NEUTRON:
138  if (icou[i] < 10) {
139  idtelescope_code = KVINDRA::ID_SI_CSI;
140  }
141  break;
142 
144  if (icou[i] > 1 && icou[i] < 10) {
145  idtelescope_code = KVINDRA::ID_CI_SI;
146  }
147  else if (icou[i] > 9) {
148  idtelescope_code = KVINDRA::ID_CI_CSI;
149  }
150  break;
151 
153  if (icou[i] > 1 && icou[i] < 10) {
154  idtelescope_code = KVINDRA::ID_CI_SI;
155  det = get_si(i);
156  }
157  else if (icou[i] > 9) {
158  idtelescope_code = KVINDRA::ID_CI_CSI;
159  det = get_csi(i);
160  }
161  break;
162 
164  if (with_phoswich && icou[i] == 1) {
165  idtelescope_code = KVINDRA::ID_PHOSWICH;
166  }
167  else {
168  idtelescope_code = KVINDRA::ID_CSI_PSA;
169  }
170  break;
171 
173  if (icou[i] > 1) {
174  idtelescope_code = KVINDRA::ID_CSI_PSA;
175  }
176  break;
177  }
178  if (!det) {
179  Warning("Process",
180  "Got no detector with id=%d cou=%d mod=%d z=%d a=%d ecod=%d ener=%f",
181  idcode[i], icou[i], imod[i], z[i], a[i], ecode[i], ener[i]);
182  continue;
183  }
184  auto node = det->GetNode();
185  auto traj = (KVGeoDNTrajectory*)node->GetTrajectories()->First();
186  bool etalon = sili_fired(i) || si75_fired(i);
187  auto chio = det->GetChIo();
188  if (node->GetNTrajForwards() > 1) {
189  // Stopped in CsI behind etalon telescope - we have a particle stopped in CsI which may have gone CI-CsI or CI-Si75-SiLi-CsI
190  if (etalon)
191  traj = (KVGeoDNTrajectory*)node->GetTrajectories()->FindObjectWithMethod("4", "GetN");
192  else
193  traj = (KVGeoDNTrajectory*)node->GetTrajectories()->FindObjectWithMethod("2", "GetN");
194  }
195  // correct idtelescope code in case of code=6 or code=10 stopped in SiLi or Si75 (1st campaign data)
197  if (stopped_in_sili(i))
198  idtelescope_code = KVINDRA::ID_SI75_SILI;
199  else if (stopped_in_si75(i))
200  idtelescope_code = KVINDRA::ID_CI_SI75;
201  }
202 
203  if (!node->GetTrajectories()) {
204  Warning("Process", "node %s has no trajectories", node->GetName());
205  }
206  auto rtraj = (KVReconNucTrajectory*)det->GetGroup()->GetTrajectoryForReconstruction(traj, node);
207  TIter next_id(rtraj->GetIDTelescopes());
208  KVIDTelescope* idt = nullptr;
209  while ((idt = (KVIDTelescope*)next_id())) {
210  if (idt->GetIDCode() == idtelescope_code) break;
211  }
212  if (!idt || idt->GetIDCode() != idtelescope_code) {
213  idt = (KVIDTelescope*)rtraj->GetIDTelescopes()->First();
214  }
215  // idt may still = nullptr here if e.g. particle stopped in ChIo (idcode=stopped_in_first_stage)
216  // i.e. normally only in cases where no identification is available
217 
218  auto nuc = the_event->AddParticle();
219  nuc->SetParameter("ARRAY", "INDRA");
220  nuc->SetMassFormula(KVNucleus::kVedaMass);
221  nuc->SetReconstructionTrajectory(rtraj);
222  if (phos_fired(i)) {
223  the_event->SetParameter(Form("ACQPAR.INDRA.PHOS_%02d_R", imod[i]), ci_gg[i]);
224  the_event->SetParameter(Form("ACQPAR.INDRA.PHOS_%02d_L", imod[i]), ci_pg[i]);
225  }
226  if (ci_fired(i)) {
227  if (chio) {
228  the_event->SetParameter(Form("ACQPAR.INDRA.%s_GG", chio->GetName()), ci_gg[i]);
229  the_event->SetParameter(Form("ACQPAR.INDRA.%s_PG", chio->GetName()), ci_pg[i]);
230  }
231  };
232  if (si_fired(i)) {
233  the_event->SetParameter(Form("ACQPAR.INDRA.SI_%02d%02d_GG", icou[i], imod[i]), si_gg[i]);
234  the_event->SetParameter(Form("ACQPAR.INDRA.SI_%02d%02d_PG", icou[i], imod[i]), si_pg[i]);
235  }
236  if (csi_fired(i)) {
237  the_event->SetParameter(Form("ACQPAR.INDRA.CSI_%02d%02d_R", icou[i], imod[i]), csi_r[i]);
238  the_event->SetParameter(Form("ACQPAR.INDRA.CSI_%02d%02d_L", icou[i], imod[i]), csi_l[i]);
239  }
240 
241  // set up identification result
242  auto idr = nuc->GetIdentificationResult(1);
243  idr->IDattempted = true;
244  if (idt) {
245  idr->SetIDType(idt->GetType());
246  idr->IDOK = true;
247  idr->Zident = true;
248  }
249  idr->IDcode = idcode[i];
250  // correction for 1st campaign data: ALL masses are positive!!!
251  if (camp1) {
252  // after inspection of data, only CsI R/L provides isotopic identification (for Z<5)
253  if (idcode[i] != 2) {
254  a_indra[i] = -a_indra[i];
255  a[i] = -a[i];
256  }
257  }
258  idr->Aident = (a[i] > 0);
259  if (idr->Aident) idr->PID = a_indra[i];
260  else idr->PID = z_indra[i];
261  idr->Z = z[i];
262  idr->A = TMath::Abs(a[i]);
268  nuc->SetIsIdentified();
269  nuc->SetIdentification(idr, idt);
270 
271  if (ecode[i] > 0 && ecode[i] < 4 && ener[i] > 0) {
272  nuc->SetIsCalibrated();
273  nuc->SetECode(ecode[i]);
274  nuc->SetEnergy(ener[i]);
275  nuc->GetAnglesFromReconstructionTrajectory();
276  if (phos_fired(i)) {
277  nuc->SetParameter("INDRA.EPHOS", de1[i]);
278  }
279  if (ci_fired(i)) {
280  nuc->SetParameter("INDRA.ECHIO", de1[i]);
281  nuc->SetParameter("INDRA.DE_MYLAR", de_mylar[i]);
282  };
283  if (si_fired(i)) {
284  nuc->SetParameter("INDRA.ESI", de2[i]);
285  }
286  if (csi_fired(i)) {
287  nuc->SetParameter("INDRA.ECSI", de3[i]);
288  }
289  if (si75_fired(i)) nuc->SetParameter("INDRA.ESI75", de4[i]);
290  if (sili_fired(i)) nuc->SetParameter("INDRA.ESILI", de5[i]);
291  }
292  else
293  nuc->SetIsUncalibrated();
294  }
295 
296  output_tree->Fill();
297 
298  ++events_read;
299  if (!(events_read % 1000)) std::cout << "Read " << events_read << "/" << total_events << " events...\n";
300 
301  return kTRUE;
302 }
303 
304 
305 
310 
312 {
313  // The SlaveTerminate() function is called after all entries or objects
314  // have been processed. When running with PROOF SlaveTerminate() is called
315  // on each slave server.
316 
317 }
318 
319 
320 
325 
327 {
328  // The Terminate() function is the last function to be called during
329  // a query. It always runs on the client, it can be used to present
330  // the results graphically or save the results to file.
331 
332  output_file->cd();
333  gDataAnalyser->WriteBatchInfo(output_tree);
334  output_tree->Write(); //write tree to file
335 
336  //add new file to repository
337  gDataSet->CommitRunfile("recon", run_index, output_file);
338 }
339 
340 
unsigned int UInt_t
bool Bool_t
constexpr Bool_t kTRUE
Option_t Option_t option
char * Form(const char *fmt,...)
virtual const Char_t * GetType() const
Definition: KVBase.h:176
virtual Bool_t IsCalled(const Char_t *name) const
Definition: KVBase.h:189
void WriteBatchInfo(TTree *)
void CommitRunfile(const KVString &type, const run_index_t &run, TFile *file)
Definition: KVDataSet.cpp:1372
TFile * NewRunfile(const KVString &type, const run_index_t &run)
Definition: KVDataSet.cpp:1112
KVNameValueList * GetParameters() const
Definition: KVEvent.h:179
static void MakeEventBranch(TTree *tree, const TString &branchname, T &event, Int_t bufsize=10000000)
Definition: KVEvent.h:210
void Clear(Option_t *opt="") override
Definition: KVEvent.h:238
void SetParameter(const Char_t *name, ValType value) const
Definition: KVEvent.h:197
const KVDBRunFile & GetDBRunFile(const run_index_t &r) const
Definition: KVExpDB.h:93
Path taken by particles through multidetector geometry.
Base class for all detectors or associations of detectors in array which can identify charged particl...
Definition: KVIDTelescope.h:84
virtual UShort_t GetIDCode()
TTreeReaderArray< Float_t > si_gg
TTreeReaderArray< Int_t > z
TTreeReaderArray< Float_t > de1
TTreeReaderArray< Float_t > z_indra
KVINDRADetector * get_sili(int i)
TTreeReaderValue< Int_t > npart_traite
pointer to the analyzed TTree or TChain
TTreeReaderArray< Float_t > a_indra
KVINDRADetector * get_ci(int i)
void SlaveBegin(TTree *tree) override
TTreeReaderArray< Float_t > de5
TTreeReaderArray< Float_t > de3
TTreeReaderArray< Float_t > de4
TTreeReaderArray< Float_t > ci_gg
TTreeReaderArray< Float_t > si_pg
TTreeReaderArray< Int_t > a
TTreeReaderArray< Float_t > csi_r
KVINDRADetector * get_phos(int i)
KVINDRADetector * get_si75(int i)
Bool_t Process(Long64_t entry) override
TTreeReaderArray< Float_t > csi_l
TTreeReaderArray< Float_t > ci_pg
TTreeReaderArray< Int_t > icou
KVINDRADetector * get_csi(int i)
TTreeReaderArray< Float_t > de_mylar
TTreeReaderArray< Int_t > imod
KVReconstructedEvent * the_event
TTreeReaderArray< Int_t > idcode
TTreeReaderArray< Float_t > ener
TTreeReader fReader
the tree reader
void SlaveTerminate() override
KVINDRADetector * get_si(int i)
TTreeReaderArray< Float_t > de2
TTreeReaderArray< Int_t > ecode
void Begin(TTree *tree) override
Base class for detectors of INDRA array.
KVINDRADetector * GetChIo() const
@ ID_CI_SI75
particle identified in ChIo-Si75 etalon telescope
Definition: KVINDRA.h:136
@ ID_SI_CSI
particle identified in Si-CsI telescope
Definition: KVINDRA.h:131
@ ID_NEUTRON
'neutron' discriminated by coherency between CsI and Si-CsI identifications
Definition: KVINDRA.h:128
@ ID_CI_SI
particle identified in ChIo-Si telescope
Definition: KVINDRA.h:134
@ ID_CI_SI_COHERENCY
particle identified in ChIo-Si telescope in coincidence with light particle identified in CsI
Definition: KVINDRA.h:137
@ ID_STOPPED_IN_FIRST_STAGE
particle stopped in first detector of telescope, only minimum Z can be estimated
Definition: KVINDRA.h:126
@ ID_CSI_PSA
particle identified in CsI detector by pulse shape analysis
Definition: KVINDRA.h:130
@ ID_PHOSWICH
particle identified in phoswich (campaigns 1-3)
Definition: KVINDRA.h:129
@ NO_IDENTIFICATION
no identification either attempted or available for particle
Definition: KVINDRA.h:125
@ ID_CI_COHERENCY
particle stopped in ChIo revealed by coherency tests (Zmin)
Definition: KVINDRA.h:138
@ ID_CI_MULTIHIT
particles stopped in multiple Si (ring<10) or CsI (ring>9) behind same ChIo, bad identification
Definition: KVINDRA.h:139
@ ID_CSI_FRAGMENT
particle partially identified in CsI detector, with Z greater than identifiable
Definition: KVINDRA.h:140
@ ID_GAMMA
'gamma' particle detected in CsI
Definition: KVINDRA.h:127
@ ID_CSI_MASS_OUT_OF_RANGE
particle partially identified in CsI detector, mass out of range of apparent Z (pile-up?...
Definition: KVINDRA.h:141
@ ID_SI75_SILI
particle identified in Si75-SiLi etalon telescope
Definition: KVINDRA.h:132
@ ID_CI_CSI
particle identified in ChIo-CsI telescope
Definition: KVINDRA.h:135
virtual KVINDRADetector * GetDetectorByType(UInt_t cou, UInt_t mod, UInt_t type) const
Definition: KVINDRA.cpp:240
const KVSeqCollection * GetListOfPhoswich() const
Definition: KVINDRA.h:247
virtual void SetParameters(UInt_t n, Bool_t physics_parameters_only=kFALSE)
static KVMultiDetArray * MakeMultiDetector(const Char_t *dataset_name, Int_t run=-1, TString classname="KVMultiDetArray")
void IncrementValue(const Char_t *name, value_type value)
void SetParameter(const Char_t *name, ValType value) const
Definition: KVParticle.h:822
Path through detector array used to reconstruct detected particle.
Event containing KVReconstructedNucleus nuclei reconstructed from hits in detectors.
Particle * AddParticle()
Handle list of options and input parameters for user analyis .
TString GetOpt(const Char_t *opt) const
void ParseOptions(const KVString &optlist)
virtual Bool_t IsEmpty() const
Bool_t cd() override
const char * GetName() const override
const char * GetTitle() const override
virtual void Warning(const char *method, const char *msgfmt,...) const
virtual void Error(const char *method, const char *msgfmt,...) const
virtual void Info(const char *method, const char *msgfmt,...) const
const char * GetOption() const override
Int_t Atoi() const
EEntryStatus SetLocalEntry(Long64_t entry)
virtual Int_t Fill()
Int_t Write(const char *name=nullptr, Int_t option=0, Int_t bufsize=0) const override
void set_run(int r)
Definition: run_index.h:65
long long Long64_t
Double_t Abs(Double_t d)