er  dev
ERTelescopePID.cxx
1 /********************************************************************************
2  * Copyright (C) Joint Institute for Nuclear Research *
3  * *
4  * This software is distributed under the terms of the *
5  * GNU Lesser General Public Licence version 3 (LGPL) version 3, *
6  * copied verbatim in the file "LICENSE" *
7  ********************************************************************************/
8 
9 #include "ERTelescopePID.h"
10 
11 #include <numeric>
12 
13 #include "G4IonTable.hh"
14 #include "G4ParticleDefinition.hh"
15 #include "G4EmCalculator.hh"
16 #include "G4NistManager.hh"
17 
18 #include "TVector3.h"
19 #include "TMath.h"
20 #include "TGeoNode.h"
21 #include "TGeoManager.h"
22 
23 #include "FairRootManager.h"
24 #include "FairRuntimeDb.h"
25 #include "FairLogger.h"
26 #include "FairEventHeader.h"
27 
28 #include "ERBeamDetTrack.h"
29 #include "ERRunAna.h"
30 #include "ERSupport.h"
31 
32 //--------------------------------------------------------------------------------------------------
34  : ERTask("ER qtelescope particle identification scheme") {
35  fAvailibleRunManagers.push_back("ERRunAna");
36 }
37 //--------------------------------------------------------------------------------------------------
39  : ERTask("ER qtelescope particle identification scheme", verbose) {
40  fAvailibleRunManagers.push_back("ERRunAna");
41 }
42 //--------------------------------------------------------------------------------------------------
43 InitStatus ERTelescopePID::Init() {
44  if (ERTask::Init() != kSUCCESS)
45  return kFATAL;
46  FairRootManager* ioman = FairRootManager::Instance();
47  if ( ! ioman )
48  Fatal("Init", "No FairRootManager");
49  TList* allbrNames = ioman->GetBranchNameList();
50  TIter nextBranch(allbrNames);
51  TObjString* bName;
52  while (bName = (TObjString*)nextBranch()) {
53  TString bFullName = bName->GetString();
54  if (bFullName.Contains("Digi") && bFullName.Contains("Telescope")) {
55  fQTelescopeDigi[bFullName] = (TClonesArray*) ioman->GetObject(bFullName);
56  LOG(DEBUG) << "Digi branch " << bFullName << FairLogger::endl;
57  }
58  if (bFullName.Contains("Track") && bFullName.Contains("Telescope")) {
59  Int_t bPrefixNameLength = bFullName.First('_');
60  TString brName(bFullName(bPrefixNameLength + 1, bFullName.Length()));
61  if (fParticleDescriptions.find(brName) == fParticleDescriptions.end())
62  continue;
63  fQTelescopeTrack[brName] = (TClonesArray*) ioman->GetObject(bFullName);
64  //Creating particles collections
65  for (auto particleDescription : fParticleDescriptions[brName]){
66  TString brParticleName;
67  const auto pdg = particleDescription.fPDG;
68  brParticleName.Form("%s_%d",brName.Data(), pdg);
69  fQTelescopeParticle[brName][pdg] = new TClonesArray("ERTelescopeParticle",
70  consts::approx_telescope_particle_number);
71  ioman->Register("TelescopeParticle_" + brParticleName, "Telescope",
72  fQTelescopeParticle[brName][pdg], kTRUE);
73  }
74  }
75  }
76  fQTelescopeSetup = ERTelescopeSetup::Instance();
77  fQTelescopeSetup->ReadGeoParamsFromParContainer();
78  return kSUCCESS;
79 }
80 //--------------------------------------------------------------------------------------------------
81 void ERTelescopePID::SetParticle(
82  const TString& trackBranchName, const PDG pdg,
83  const TString& deStation /*= ""*/, const TString& eStation /*= ""*/,
84  const Double_t deNormalizedThickness /*= 0.002*/,
85  const std::vector<TString>& stations_to_use_em_calculator_for_de_e /* = {}*/,
86  const std::vector<TString>& stations_to_use_em_calculator_for_kinetic_energy/*= {}*/) {
87  if ((deStation == "" && eStation != "") || (deStation != "" && eStation == ""))
88  LOG(FATAL) << "If one of dE or E station is defined, another should be defined."
89  << FairLogger::endl;
90  std::list<TString> eStations;
91  if (eStation != "" )
92  eStations.push_back(eStation);
93  fParticleDescriptions[trackBranchName].emplace_back(pdg, deStation, eStations,
94  deNormalizedThickness,
95  stations_to_use_em_calculator_for_de_e,
96  stations_to_use_em_calculator_for_kinetic_energy);
97 }
98 //--------------------------------------------------------------------------------------------------
99 void ERTelescopePID::SetParticle(
100  const TString& trackBranchName, const PDG pdg,
101  const TString& deStation /*= ""*/, const std::list<TString>& eStations /*= {}*/,
102  const Double_t deNormalizedThickness /*= 0.002*/,
103  const std::vector<TString>& stations_to_use_em_calculator_for_de_e /* = {}*/,
104  const std::vector<TString>& stations_to_use_em_calculator_for_kinetic_energy/*= {}*/) {
105  fParticleDescriptions[trackBranchName].emplace_back(pdg, deStation, eStations,
106  deNormalizedThickness,
107  stations_to_use_em_calculator_for_de_e,
108  stations_to_use_em_calculator_for_kinetic_energy);
109 }
110 //--------------------------------------------------------------------------------------------------
111 void ERTelescopePID::Exec(Option_t* opt) {
112  LOG(DEBUG) << "[ERTelescopePID]------------Started-----------------------------------------------"
113  << FairLogger::endl;
114  Reset();
115  for (const auto& itTrackBranches : fQTelescopeTrack) {
116  const auto& trackBranch = itTrackBranches.first;
117  const auto& tracks = itTrackBranches.second;
118  if (tracks->GetEntriesFast() == 0)
119  continue;
120  LOG(DEBUG) << "[ERTelescopePID] Work with traсks in " << trackBranch << FairLogger::endl;
121  for (Int_t iTrack(0); iTrack < tracks->GetEntriesFast(); iTrack++) {
122  const auto* track = static_cast<ERTelescopeTrack*>(tracks->At(iTrack));
123  for (const auto particleDescription : fParticleDescriptions[trackBranch]) {
124  // Get particle mass by PDG from G4IonTable
125  const auto pdg = particleDescription.fPDG;
126  const auto* particle = G4ParticleTable::GetParticleTable()->FindParticle(pdg);
127  if (!particle)
128  LOG(FATAL) << "Particle with code " << pdg << " not found in Geant database "<< FairLogger::endl;
129  const auto mass = particle->GetPDGMass(); // MeV
130  // Find geometry point to start track back propagation.
131  const auto back_propagation_start_point = FindBackPropagationStartPoint(*track);
132  // Calc particle energy deposites in setup using back propagation from start point
133  // to middle of target volume. We assume that kinetic energy
134  // is fully detected by the setup, so calculate it as sum of
135  // energy deposites in digi and energy deposites in passive volumes (using G4EmCalculator).
136  std::list<DigiOnTrack> digisOnTrack;
137  const auto energyDeposites = CalcEnergyDeposites(
138  *track, back_propagation_start_point, *particle, digisOnTrack,
139  particleDescription.stations_to_use_em_calculator_for_de_e_,
140  particleDescription.stations_to_use_em_calculator_for_kinetic_energy_);
141  Double_t edepInThickStation = -1., edepInThinStation = -1.,
142  edepInThickStationCorrected = -1., edepInThinStationCorrected = -1.;
143  ERChannel channelOfThinStation = consts::undefined_channel;
144  ERChannel channelOfThickStation = consts::undefined_channel;
145  FindEnergiesForDeEAnalysis(trackBranch, digisOnTrack,
146  particleDescription.fEStations, particleDescription.fDeStation,
147  particleDescription.fNormalizedThickness, edepInThickStation, edepInThinStation,
148  edepInThickStationCorrected, edepInThinStationCorrected,
149  channelOfThinStation, channelOfThickStation);
150  const auto digisDeposite = energyDeposites.first;
151  const auto deadDeposite = energyDeposites.second;
152  const auto kinetic_energy = digisDeposite + deadDeposite;
153  const Double_t fullEnergy = mass + kinetic_energy;
154  // Calculate Lorentz vector of particle at the exit of the reaction in lab system.
155  const Double_t momentumMag = sqrt(pow(fullEnergy, 2) - pow(mass, 2) );
156  TVector3 direction = track->GetDirection();
157  const auto momentum = momentumMag * direction;
158  const TLorentzVector lvTarget (momentum, fullEnergy);
159  // Add particle to collection.
160  auto& particleCollection = *fQTelescopeParticle[trackBranch][pdg];
161  AddParticle(lvTarget, kinetic_energy, deadDeposite, edepInThickStation, edepInThinStation,
162  edepInThickStationCorrected, edepInThinStationCorrected,
163  channelOfThinStation, channelOfThickStation, particleCollection);
164  }
165  }
166  }
167  LOG(DEBUG) << "[ERTelescopePID]------------Finished----------------------------------------------"
168  << FairLogger::endl;
169 }
170 //--------------------------------------------------------------------------------------------------
172  for (const auto itTrackBranches : fQTelescopeParticle) {
173  for (const auto itParticleBranches : itTrackBranches.second)
174  if (itParticleBranches.second) {
175  itParticleBranches.second->Delete();
176  }
177  }
178 }
179 //--------------------------------------------------------------------------------------------------
181 AddParticle(const TLorentzVector& lvInteraction, const Double_t kinetic_energy, const Double_t deadEloss,
182  Double_t edepInThickStation, Double_t edepInThinStation,
183  Double_t edepInThickStationCorrected, Double_t edepInThinStationCorrected,
184  const ERChannel channelOfThinStation, const ERChannel channelOfThickStation,
185  TClonesArray& col) {
186  return new(col[col.GetEntriesFast()]) ERTelescopeParticle(lvInteraction, kinetic_energy, deadEloss,
187  edepInThickStation, edepInThinStation,
188  edepInThickStationCorrected, edepInThinStationCorrected,
189  channelOfThinStation, channelOfThickStation);
190 }
191 //--------------------------------------------------------------------------------------------------
192 TVector3 ERTelescopePID::FindBackPropagationStartPoint(const ERTelescopeTrack& track) {
193  // Return geometry point at which the track exit last sensetive volume
194  // on its direct propagation. If sensetive volume was not found, return
195  // track telescope vertex.
196  const TVector3 telescopeVertex = track.GetXStationVertex();
197  TVector3 back_propagation_start_point = telescopeVertex;
198  const auto direction = track.GetDirection();
199  TGeoNode* currentNode = gGeoManager->InitTrack(
200  telescopeVertex.X(), telescopeVertex.Y(), telescopeVertex.Z(),
201  direction.X(), direction.Y(), direction.Z());
202  TGeoNode* lastSensetiveNode = nullptr;
203  TString lastSensetivePath;
204  TVector3 lastSensetivePosition;
205  while(!gGeoManager->IsOutside()) {
206  bool inSensetiveVolume = false;
207  const TString path = gGeoManager->GetPath();
208  if (path.Contains("Sensitive")) {
209  // Enter sensetive volume. Next step will be in sensetive volume.
210  inSensetiveVolume = true;
211  lastSensetiveNode = currentNode;
212  lastSensetivePath = path;
213  }
214  currentNode = gGeoManager->FindNextBoundary();
215  currentNode = gGeoManager->Step();
216  if (inSensetiveVolume) {
217  // position, when track exit sensetve volume
218  lastSensetivePosition = TVector3(gGeoManager->GetCurrentPoint());
219  }
220  }
221  if (lastSensetiveNode) {
222  back_propagation_start_point = lastSensetivePosition;
223  LOG(DEBUG) << "[FindBackPropagationStartPoint] Last sensetive volume for track "
224  << lastSensetivePath << FairLogger::endl;
225  } else {
226  LOG(DEBUG) << "[FindBackPropagationStartPoint] Last sensetive volume for track not found. "
227  << "Track telescope vertex will be used as start point for back propagation\n";
228  }
229  LOG(DEBUG) << "[FindBackPropagationStartPoint] Back propagation start point "
230  << "(" << back_propagation_start_point.X() << "," << back_propagation_start_point.Y()
231  << "," << back_propagation_start_point.Z() << ")" << FairLogger::endl;
232  return back_propagation_start_point;
233 }
234 //--------------------------------------------------------------------------------------------------
235 std::pair<Double_t, Double_t> ERTelescopePID::
236 CalcEnergyDeposites(const ERTelescopeTrack& track, const TVector3& start_point,
237  const G4ParticleDefinition& particle,
238  std::list<DigiOnTrack>& digisOnTrack,
239  const std::vector<TString>& stations_to_use_em_calculator_for_de_e,
240  const std::vector<TString>& stations_to_use_em_calculator_for_kinetic_energy) {
241  // Calc paritcle energy deposites in setup using back track propagation from start point.
242  // Return pair: first - sum of energy deposites in digi(sensetive volumes);
243  // second - sum of energy deposites in passive volumes (or dead energy deposite).
244  Double_t digiDepositesSum = 0.;
245  Double_t deadDepositesSum = 0.;
246  // We start with kinetic energy equal zero, because we assume
247  // that start point is an exit point of the last sensetive volume,
248  // which track has passed in setup. If the track stopped earlier,
249  // it will be taken into account automatically, because particle
250  // with zero kinetic energy can not loss energy ;)
251  Double_t kinetic_energy = 0.;
252  // Init track in back direction.
253  auto back_direction = track.GetBackDirection();
254  const TString log_prefix = "[ERTelescopePID] [CalcEnergyDeposites] ";
255  LOG(DEBUG) << log_prefix << "Energy deposites calculation for particle "
256  << particle.GetParticleName() << "; start point = (" << start_point.X() << ","
257  << start_point.Y() << "," << start_point.Z() << "; and direction = "
258  << back_direction.X() << "," << back_direction.Y() << "," << back_direction.Z()
259  << FairLogger::endl;
260  const auto calc_edep_using_em_calulator =
261  [&log_prefix, &particle](TGeoNode* node, float step, float T) {
262  const TString materialName = node->GetMedium()->GetMaterial()->GetName();
263  const auto* material = G4NistManager::Instance()->FindOrBuildMaterial(materialName.Data());
264  LOG(DEBUG) << log_prefix <<" Calc energy deposite for range " << step << " in material "
265  << materialName << " with kinetic energy " << T << FairLogger::endl;
266  return CalcElossIntegralVolStep(T, particle, *material, step);
267  };
268  TGeoNode* node = gGeoManager->InitTrack(start_point.X(), start_point.Y(), start_point.Z(),
269  back_direction.X(), back_direction.Y(), back_direction.Z());
270  // While track not in interaction volume or outside the setup,
271  // accumulate energy deposites and kinetic energy.
272  while(!gGeoManager->IsOutside()) {
273  gGeoManager->FindNextBoundary();
274  const TVector3 current_position(gGeoManager->GetCurrentPoint());
275  LOG(DEBUG) << log_prefix << "track position (" << current_position.X() << ", "
276  << current_position.Y() << ", " << current_position.Z() << ")" << FairLogger::endl;
277  LOG(DEBUG) << log_prefix << "path = "
278  << gGeoManager->GetPath() << FairLogger::endl;
279  // If track in sensetive volume, try to find digi
280  const TString nodePath = TString(gGeoManager->GetPath());
281  bool is_last_step = false;
282  if (nodePath.Contains("Sensitive")){
283  const Bool_t use_em_calcaulator_for_de_e =
284  std::find_if(stations_to_use_em_calculator_for_de_e.begin(),
285  stations_to_use_em_calculator_for_de_e.end(),
286  [&nodePath](const TString& station) {return nodePath.Contains(station);})
287  != stations_to_use_em_calculator_for_de_e.end();
288  const Bool_t use_em_calcaulator_for_kinetic_energy =
289  std::find_if(stations_to_use_em_calculator_for_kinetic_energy.begin(),
290  stations_to_use_em_calculator_for_kinetic_energy.end(),
291  [&nodePath](const TString& station) {return nodePath.Contains(station);})
292  != stations_to_use_em_calculator_for_kinetic_energy.end();
293  LOG(DEBUG) << log_prefix <<" Finding digis..." << FairLogger::endl;
294  auto branchAndDigis = FindDigisByNode(*node, gGeoManager->GetPath());
295  Bool_t digis_for_node_are_found = !branchAndDigis.empty();
296  float edep_using_em_calculator = -1.;
297  if (!digis_for_node_are_found
298  || use_em_calcaulator_for_de_e
299  || use_em_calcaulator_for_kinetic_energy) {
300  edep_using_em_calculator = calc_edep_using_em_calulator(node, gGeoManager->GetStep(),
301  kinetic_energy);
302  }
303  for (auto& branchAndDigi : branchAndDigis) {
304  const auto digiBranch = branchAndDigi.first;
305  auto* digi = branchAndDigi.second;
306  if (use_em_calcaulator_for_de_e)
307  digi->SetEdep(edep_using_em_calculator);
308  digisOnTrack.emplace_back(digiBranch, digi, gGeoManager->GetStep());
309  }
310  const Double_t edepInStation = digis_for_node_are_found && !use_em_calcaulator_for_kinetic_energy
311  ? ApplyEdepAccountingStrategy(branchAndDigis)
312  : edep_using_em_calculator;
313  kinetic_energy += edepInStation;
314  digiDepositesSum += edepInStation;
315  } else {
316  // track in passive volume
317  auto step = gGeoManager->GetStep();
318  const auto target_vertex = track.GetTargetVertex();
319  const double step_to_target_vertex = (target_vertex - current_position).Mag();
320  is_last_step = step_to_target_vertex <= step;
321  step = is_last_step ? step_to_target_vertex : step;
322  const auto deadDeposite = calc_edep_using_em_calulator(node, step, kinetic_energy);
323  kinetic_energy += deadDeposite;
324  deadDepositesSum += deadDeposite;
325  }
326  LOG(DEBUG) << log_prefix << "Current kinetic Energy = " << kinetic_energy << FairLogger::endl;
327  if (is_last_step)
328  break;
329  node = gGeoManager->Step();
330  }
331  LOG(DEBUG) << log_prefix << "Finish deposite calculation with Kinetic energy = "
332  << kinetic_energy << " = digis sum deposite : " << digiDepositesSum
333  << " + dead deposites sum : " << deadDepositesSum << FairLogger::endl;
334  return {digiDepositesSum, deadDepositesSum};
335 }
336 //--------------------------------------------------------------------------------------------------
337 std::map<TString, ERDigi*> ERTelescopePID::FindDigisByNode(const TGeoNode& node, const TString& nodePath) {
338  // example SSD: /cave_1/QTelescopeTmp_0/Telescope_4_5/Telescope_4_SingleSi_SSD_V_4_X_11/SensitiveSingleSiStrip_X_2
339  // example DSD: /cave_1/QTelescopeTmp_0/CT_3/CT_DoubleSi_DSD_XY_0/doubleSiStrip_XY_0/SensitiveDoubleSiBox_XY_16
340  // example CsI: /cave_1/QTelescopeTmp_0/Central_telescope_3/Central_telescope_CsI_0/CsIBoxShell_8/SensitiveCsIBox_1
341  // example NonUniform SSD: /cave_1/QTelescopeTmp_0/Telescope_1_1/Telescope_1_SingleSi_SSD20_1_X_0/pseudoSiStrip_4_4/SensitivePixelSiBox_X4_Y1_0
342  // example R SSD: /cave_1/QTelescopeTmp_0/Telescope_1_1/Telescope_1_SingleSi_Phi1_X_0/r_station_0/Sensitivestrip_12
343  // example R DSD: /cave_1/QTelescopeTmp_0/Telescope_1_1/Telescope_1_DoubleSi_R_XY_0/r_station_0/Segment_3/SensitiveSegment_1
344  // example R CsI: /cave_1/QTelescopeTmp_0/Telescope_1_1/Telescope_1_CsI_Radial_0/CrystalShell_7/SensitiveCrystal_0
345  std::map<TString, ERDigi*> resultDigis;
346  const bool nodeOfDoubleSiStation = nodePath.Contains("DoubleSi");
347  //@TODO here Setup->GetComponent(path) interface should be used in future.
348  const auto getDigiBranchSubString = [&node, &nodePath, nodeOfDoubleSiStation]() -> TString {
349  TString digiBranchSubString = nodePath;
350  digiBranchSubString.Remove(digiBranchSubString.Last('/'), digiBranchSubString.Length());
351  if (nodeOfDoubleSiStation || nodePath.Contains("CsI") || nodePath.Contains("pseudo"))
352  digiBranchSubString.Remove(digiBranchSubString.Last('/'), digiBranchSubString.Length());
353  if (digiBranchSubString.EndsWith("r_station_0"))
354  digiBranchSubString.Remove(digiBranchSubString.Last('/'), digiBranchSubString.Length());
355  digiBranchSubString.Remove(0, digiBranchSubString.Last('/') + 1);
356  // remove node id
357  digiBranchSubString.Remove(digiBranchSubString.Last('_'), digiBranchSubString.Length());
358  return digiBranchSubString; // for ex: Telescope_4_SingleSi_SSD_V_4_X
359  };
360  const auto digiBranchSubstring = getDigiBranchSubString();
361  LOG(DEBUG) <<" [ERTelescopePID][CalcEnergyDeposites] Digi branch substring "
362  << digiBranchSubstring << FairLogger::endl;
363  // Multiple digi collections can correspond to one geometry node.
364  // DoubleSi station is example.
365  //@TODO here component->GetBranchNames() interface should be used in future.
366  std::list<TString> digiBranchNames;
367  for (const auto& digiCollectionPair : fQTelescopeDigi) {
368  if (digiCollectionPair.first.Contains(digiBranchSubstring))
369  digiBranchNames.push_back(digiCollectionPair.first);
370  }
371  if (digiBranchNames.empty()) {
372  LOG(ERROR) << " [ERTelescopePID][CalcEnergyDeposites] Branch with substring " << digiBranchSubstring
373  << " not found." << FairLogger::endl;
374  return resultDigis;
375  }
376  //@TODO here component->GetChannelFromSensetiveNode() interface should be used in future.
377  const auto getChannels = [&node, &nodePath, nodeOfDoubleSiStation]() -> std::vector<Int_t> {
378  TString pathWithChannelPostfix = nodePath;
379  if (nodePath.Contains("CsI") || nodePath.Contains("pseudo"))
380  pathWithChannelPostfix.Remove(pathWithChannelPostfix.Last('/'), pathWithChannelPostfix.Length());
381  const TString channelStr(pathWithChannelPostfix(pathWithChannelPostfix.Last('_') + 1,
382  pathWithChannelPostfix.Length()));
383  if (!nodeOfDoubleSiStation)
384  return {channelStr.Atoi()};
385  pathWithChannelPostfix.Remove(pathWithChannelPostfix.Last('/'), pathWithChannelPostfix.Length());
386  const TString secondChannelStr(pathWithChannelPostfix(pathWithChannelPostfix.Last('_') + 1,
387  pathWithChannelPostfix.Length()));
388  return {channelStr.Atoi(), secondChannelStr.Atoi()};
389  };
390  const auto channels = getChannels();
391  for (const auto& digiBranchName : digiBranchNames) {
392  auto* digis = fQTelescopeDigi[digiBranchName];
393  // @BUG Code only for XY doubleSi stations.
394  Int_t channel = (nodeOfDoubleSiStation && digiBranchName.EndsWith("_X")) ? channels[1] : channels[0];
395  LOG(DEBUG) << " [ERTelescopePID][CalcEnergyDeposites] Finding digi with channel number "
396  << channel << " in branch " << digiBranchName << FairLogger::endl;
397  Bool_t digiFound = false;
398  for (Int_t iDigi = 0; iDigi < digis->GetEntriesFast(); iDigi++){
399  auto* digi = static_cast<ERDigi*>(digis->At(iDigi));
400  LOG(DEBUG) << " [ERTelescopePID][CalcEnergyDeposites] Collection has digi with channel "
401  << digi->Channel() << " and edep " << digi->Edep() << FairLogger::endl;
402  if (digi->Channel() != channel)
403  continue;
404  digiFound = true;
405  LOG(DEBUG) << " [ERTelescopePID][CalcEnergyDeposites] Found digi with edep "
406  << digi->Edep() << FairLogger::endl;
407  resultDigis[digiBranchName] = digi;
408  break;
409  }
410  if (!digiFound)
411  LOG(ERROR) << " [ERTelescopePID][CalcEnergyDeposites] Digi with channel number "
412  << channel << " not found in collection " << digiBranchName << FairLogger::endl;
413  }
414  return resultDigis;
415 }
416 //--------------------------------------------------------------------------------------------------
417 void ERTelescopePID::FindEnergiesForDeEAnalysis(const TString& trackBranch,
418  const std::list<DigiOnTrack>& digisOnTrack, const std::list<TString>& eStations,
419  const TString& deStation, const Double_t normalizedThickness, Double_t& edepInThickStation,
420  Double_t& edepInThinStation, Double_t& edepInThickStationCorrected,
421  Double_t& edepInThinStationCorrected, ERChannel& channelOfThinStation,
422  ERChannel& channelOfThickStation) {
423  const TString log_prefix = "[ERTelescopePID][FindEnergiesForDeEAnalysis] ";
424  if (eStations.empty() || deStation == "")
425  return;
426  const auto getDigisOnTrack = [&digisOnTrack](const TString& stationName) -> std::list<DigiOnTrack> {
427  std::list<DigiOnTrack> digisOnTrackInStation;
428  for (const auto digiOnTrack : digisOnTrack) {
429  if (digiOnTrack.fBranch.Contains(stationName)) {
430  digisOnTrackInStation.push_back(digiOnTrack);
431  }
432  }
433  return digisOnTrackInStation;
434  };
435  const auto deDigisOnTrack = getDigisOnTrack(deStation);
436  if (deDigisOnTrack.empty()) {
437  LOG(WARNING) << log_prefix << "Digi for station " << deStation
438  << " not found on track from " << trackBranch << " path." << FairLogger::endl;
439  return;
440  }
441  channelOfThinStation = deDigisOnTrack.front().fDigi->Channel();
442  TString eStations_string;
443  for (const auto& eStation : eStations) {
444  const auto eDigisOnTrack = getDigisOnTrack(eStation);
445  if (eDigisOnTrack.empty()) {
446  LOG(WARNING) << log_prefix << "Digi for station " << eStation
447  << " not found on track from " << trackBranch << " path" << FairLogger::endl;
448  continue;
449  }
450  if (edepInThickStation == -1.) {
451  edepInThickStation = ApplyEdepAccountingStrategy(eDigisOnTrack);
452  eStations_string += eStation;
453  } else {
454  edepInThickStation += ApplyEdepAccountingStrategy(eDigisOnTrack);
455  eStations_string += TString(", ") + eStation;
456  }
457  //TODO: Works only for single e station
458  channelOfThickStation = eDigisOnTrack.front().fDigi->Channel();
459  }
460  if (edepInThickStation == -1.) {
461  LOG(WARNING) << log_prefix << "No digis found for E stations in de-E analysis." << FairLogger::endl;
462  return;
463  }
464  edepInThinStation = ApplyEdepAccountingStrategy(deDigisOnTrack);
465  const auto deSensetiveThickness = deDigisOnTrack.front().fSensetiveThickness;
466  // de and E correction
467  edepInThinStationCorrected = edepInThinStation * normalizedThickness / deSensetiveThickness;
468  LOG(DEBUG) << log_prefix << "Digi from station " << deStation
469  << " with edep = " << edepInThinStation << " [MeV] corrected to edep = "
470  << edepInThinStationCorrected << ". Normalized thickness = "
471  << normalizedThickness << ", step in sensetive volume = " << deSensetiveThickness
472  << FairLogger::endl;
473  edepInThickStationCorrected = edepInThickStation + edepInThinStation - edepInThinStationCorrected;
474  LOG(DEBUG) << log_prefix << "Digis from stations " << eStations_string
475  << " with edep = " << edepInThickStation
476  << " corrected to " << edepInThickStationCorrected << FairLogger::endl;
477 }
478 //--------------------------------------------------------------------------------------------------
479 Double_t ERTelescopePID::ApplyEdepAccountingStrategy(const std::map<TString, ERDigi*>& digisByBranchName) {
480  if (digisByBranchName.empty())
481  return 0;
482  // @TODO setup->GetComponent(branchName) and component->ComponentName() should be used here
483  const auto getStrategyForStation = [this](const TString& branchName) -> EdepAccountingStrategy {
484  const auto itStationStrategy = std::find_if(
485  fEdepAccountingStrategies.begin(), fEdepAccountingStrategies.end(),
486  [&branchName](const std::pair<TString, EdepAccountingStrategy> stationAndStrategy) {
487  return branchName.Contains(stationAndStrategy.first);
488  });
489  if (itStationStrategy != fEdepAccountingStrategies.end())
490  return itStationStrategy->second;
491  if (branchName.Contains("_X"))
492  return EdepFromXChannel;
493  else if (branchName.Contains("_Y"))
494  return EdepFromYChannel;
495  if (branchName.Contains("CsI"))
496  return SummarizedEdep;
497  LOG(FATAL) << "Unable to select edep accaunting strategy by branch " << branchName << FairLogger::endl;
498  };
499  const auto strategy = getStrategyForStation(digisByBranchName.begin()->first);
500  Double_t summarizedEdep = 0.;
501  for (const auto& branchAndDigi : digisByBranchName) {
502  const auto branchName = branchAndDigi.first;
503  const auto* digi = branchAndDigi.second;
504  // @TODO component->GetBranchName(...) should be used here
505  if ((strategy == EdepFromXChannel && branchName.Contains("_X"))
506  || (strategy == EdepFromYChannel && branchName.Contains("_Y"))) {
507  return digi->Edep();
508  }
509  summarizedEdep += digi->Edep();
510  }
511  if (strategy == AverageEdep)
512  return summarizedEdep / digisByBranchName.size();
513  else if (strategy == SummarizedEdep)
514  return summarizedEdep;
515  LOG(FATAL) << "[ERTelescopePID][ApplyEdepAccountingStrategy] Unable to apply edep accounting strategy." << FairLogger::endl;
516  return -1;
517 }
518 //--------------------------------------------------------------------------------------------------
519 Double_t ERTelescopePID::ApplyEdepAccountingStrategy(const std::list<DigiOnTrack>& digisOnTrack) {
520  std::map<TString, ERDigi*> digisByBranchName;
521  for (auto& digiOnTrack : digisOnTrack) {
522  digisByBranchName[digiOnTrack.fBranch] = digiOnTrack.fDigi;
523  }
524  return ApplyEdepAccountingStrategy(digisByBranchName);
525 }
526 //--------------------------------------------------------------------------------------------------
527 ClassImp(ERTelescopePID)
Definition: ERDigi.h:15
virtual void Exec(Option_t *opt)
Defines the transformation actions for all input data (Digi) to output data (Track) for each event...
virtual void Reset()
Resets all output data.
ERTelescopeSetup * fQTelescopeSetup
access to ERTelescopeSetup class instance
virtual InitStatus Init()
Defines all input and output object colletions participates in track finding.
virtual InitStatus Init()
Definition: ERTask.cxx:31
ERTelescopePID()
Default constructor.
ERTelescopeParticle * AddParticle(const TLorentzVector &lvInteraction, Double_t kineticEnergy, Double_t deadEloss, Double_t edepInThickStation, Double_t edepInThinStation, Double_t edepInThickStationCorrected, Double_t edepInThinStationCorrected, ERChannel channelOfThinStaion, ERChannel channelOfThickStation, TClonesArray &col)
Adds a ERTelescopeParticles to the output Collection.
std::vector< TString > fAvailibleRunManagers
Run managers that are availible for this task.
Definition: ERTask.h:64
Base abstract class for all tasks in er.
Definition: ERTask.h:27