CbmRoot
CbmUnigenGenerator.cxx
Go to the documentation of this file.
1 
9 #include "CbmUnigenGenerator.h"
10 
11 #include "FairLogger.h"
12 #include "FairMCEventHeader.h"
13 #include "FairPrimaryGenerator.h"
14 #include "FairRunSim.h"
15 #include "TFile.h"
16 #include "TRandom.h"
17 #include "TTree.h"
18 #include "TVector3.h"
19 #include "UEvent.h"
20 #include "UParticle.h"
21 #include "URun.h"
22 #include <cassert>
23 #include <sstream>
24 
25 
26 // ----- Constructor ----------------------------------------------------
27 CbmUnigenGenerator::CbmUnigenGenerator(const char* fileName, EMode mode)
28  : FairGenerator("UnigenGenerator", "CBM generator")
29  , fFileName(fileName)
30  , fMode(mode)
31  , fPhi(0.)
32  , fIsInit(kFALSE)
33  , fFile(nullptr)
34  , fTree(nullptr)
35  , fCurrentEntry(-1)
36  , fEvent(new UEvent())
37  , fNofPrimaries(0)
38  , fNofEvents(0)
39  , fBetaCM(0.)
40  , fGammaCM(1.)
41  , fIonMap() {
42  LOG(debug) << GetName() << ": Constructor";
43  if (mode == kRotateFixed)
44  LOG(fatal) << GetName() << ": When choosing "
45  << "mode kRotateFixed, a rotation angle must be specified!";
46  if (fileName[0] == '\0') return;
47  if (!Init())
48  LOG(fatal) << GetName() << ": Error in intialisation! Aborting...";
49 }
50 // --------------------------------------------------------------------------
51 
52 
53 // ----- Constructor ----------------------------------------------------
55  EMode mode,
56  Double_t phi)
57  : FairGenerator("UnigenGenerator", "CBM generator")
58  , fFileName(fileName)
59  , fMode(mode)
60  , fPhi(phi)
61  , fIsInit(kFALSE)
62  , fFile(nullptr)
63  , fTree(nullptr)
64  , fCurrentEntry(-1)
65  , fEvent(new UEvent())
66  , fNofPrimaries(0)
67  , fNofEvents(0)
68  , fBetaCM(0.)
69  , fGammaCM(1.)
70  , fIonMap() {
71  LOG(debug) << GetName() << ": Constructor";
72  if (fileName[0] == '\0') return;
73  if (!Init())
74  LOG(fatal) << GetName() << ": Error in intialisation! Aborting...";
75 }
76 // --------------------------------------------------------------------------
77 
78 
79 // ----- Destructor -----------------------------------------------------
81  LOG(debug) << GetName() << ": Destructor";
82  CloseInput();
83  if (fEvent) delete fEvent;
84 }
85 // --------------------------------------------------------------------------
86 
87 
88 // ----- Add a primary to the event generator ---------------------------
89 void CbmUnigenGenerator::AddPrimary(FairPrimaryGenerator* primGen,
90  Int_t pdgCode,
91  const TVector3& momentum) {
92  primGen->AddTrack(
93  pdgCode, momentum.Px(), momentum.Py(), momentum.Pz(), 0., 0., 0.);
94  fNofPrimaries++;
95 }
96 // --------------------------------------------------------------------------
97 
98 
99 // ----- Close the input file -------------------------------------------
101  if (!fFile) return;
102  LOG(INFO) << GetName() << ": Closing input file " << fFileName;
103  fFile->Close();
104  fFile = nullptr;
105 }
106 // --------------------------------------------------------------------------
107 
108 
109 // ----- Read next entry from tree --------------------------------------
111 
112  if (++fCurrentEntry >= fTree->GetEntries()) {
113  if (fMode == kReuseEvents) {
114  LOG(info) << GetName()
115  << ": End of input tree - restarting with first entry";
116  fCurrentEntry = 0;
117  } //? Re-use entries
118  else {
119  LOG(info) << GetName() << ": End of input tree reached";
120  return kFALSE;
121  } //? Do not re-use entries
122  } //? End of input tree
123 
124  // Read entry
125  Int_t result = fTree->GetEntry(fCurrentEntry);
126  if (result <= 0) {
127  LOG(ERROR) << GetName() << ": Error reading entry " << fCurrentEntry
128  << " (returns " << result << ")!";
129  return kFALSE;
130  }
131 
132  return kTRUE;
133 }
134 // --------------------------------------------------------------------------
135 
136 
137 // ----- Open the input file --------------------------------------------
139 
140  // --- Check for file being already initialised
141  if (fIsInit) {
142  LOG(warning) << GetName() << ": Already initialised!";
143  return kTRUE;
144  }
145 
146  LOG(info) << GetName() << ": Initialising...";
147  std::stringstream ss;
148  switch (fMode) {
149  case kStandard: ss << " standard; rotate event plane to zero"; break;
150  case kRotateFixed:
151  ss << " rotate event plane by fixed angle " << fPhi;
152  break;
153  case kNoRotation: ss << " no event plane rotation"; break;
154  case kReuseEvents:
155  ss << " re-use events if necessary; random event plane angle";
156  break;
157  default: ss << "unkown"; break;
158  }
159  LOG(INFO) << GetName() << ": Mode " << ss.str();
160 
161 
162  // --- Try to open file
163  fFile = new TFile(fFileName, "READ");
164  if (!fFile->IsOpen()) {
165  LOG(error) << GetName() << ": Could not open input file " << fFileName;
166  return kFALSE;
167  }
168  LOG(info) << GetName() << ": Open input file " << fFileName;
169 
170  // --- Get and print run description
171  URun* run = dynamic_cast<URun*>(fFile->Get("run"));
172  if (run == nullptr) {
173  LOG(error) << GetName() << ": No run description in input file!";
174  fFile->Close();
175  return kFALSE;
176  }
177  LOG(info) << GetName() << ": Run description";
178  run->Print();
179 
180  // --- Lorentz transformation to the lab (target) system
181  Double_t mProt = 0.938272;
182  Double_t pTarg = run->GetPTarg(); // target momentum per nucleon
183  Double_t pProj = run->GetPProj(); // projectile momentum per nucleon
184  Double_t eTarg = TMath::Sqrt(pProj * pProj + mProt * mProt);
185  Double_t eProj = TMath::Sqrt(pTarg * pTarg + mProt * mProt);
186  fBetaCM = pTarg / eTarg;
187  fGammaCM = 1. / TMath::Sqrt(1. - fBetaCM * fBetaCM);
188  Double_t pBeam = fGammaCM * (pProj - fBetaCM * eProj);
189  LOG(info) << GetName() << ": sqrt(s_NN) = " << run->GetNNSqrtS()
190  << " GeV, p_beam = " << pBeam << " GeV/u";
191  LOG(info) << GetName() << ": Lorentz transformation to lab system: "
192  << " beta " << fBetaCM << ", gamma " << fGammaCM;
193 
194  // --- Get input tree and connect event object to its branch
195  fTree = dynamic_cast<TTree*>(fFile->Get("events"));
196  if (fTree == nullptr) {
197  LOG(error) << GetName() << ": No event tree in input file!";
198  fFile->Close();
199  return kFALSE;
200  }
201  fTree->SetBranchAddress("event", &fEvent);
202  LOG(info) << GetName() << ": " << fTree->GetEntries()
203  << " events in input tree";
204 
205  // --- Register ions found in the input file
206  Int_t nIons = RegisterIons();
207  LOG(info) << GetName() << ": " << nIons << " ions registered.";
208 
209  fIsInit = kTRUE;
210  return kTRUE;
211 }
212 // --------------------------------------------------------------------------
213 
214 
215 // ----- Process a composite particle -----------------------------------
216 void CbmUnigenGenerator::ProcessIon(FairPrimaryGenerator* primGen,
217  Int_t pdgCode,
218  const TVector3& momentum) {
219 
220  assert(pdgCode >= 1e9); // Should be an ion PDG
221 
222  // Since hyper-nuclei are not (yet) supported by FairRoot, their PDG
223  // is replaced by that of the non-strange analogue.
224  Int_t ionPdg = pdgCode;
225  if (GetIonLambdas(pdgCode)) {
226  ionPdg = pdgCode - GetIonLambdas(pdgCode) * kPdgLambda;
227  LOG(warn) << GetName() << ": Replacing hypernucleus (PDG " << pdgCode
228  << ") by PDG " << ionPdg;
229  } //? Lambdas in ion
230 
231  // Charged ions can be registered
232  if (GetIonCharge(ionPdg)) AddPrimary(primGen, ionPdg, momentum);
233 
234  // Neutral ions are not supported by GEANT4.
235  // They are thus decomposed into neutrons (as an approximation)
236  else {
237  Int_t mass = GetIonMass(ionPdg);
238  TVector3 neutronMom = momentum * (1. / Double_t(mass));
239  for (Int_t iNeutron = 0; iNeutron < mass; iNeutron++)
240  AddPrimary(primGen, 2112, neutronMom);
241  LOG(warn) << GetName() << ": Neutral fragment with PDG " << ionPdg
242  << " is split into " << mass << " neutrons";
243  } //? Neutral ion
244 }
245 // --------------------------------------------------------------------------
246 
247 
248 // ----- Read input event -----------------------------------------------
249 Bool_t CbmUnigenGenerator::ReadEvent(FairPrimaryGenerator* primGen) {
250 
251  // Init must have been called before
252  assert(fIsInit);
253 
254  // Get next entry from tree
255  fNofEvents++;
256  if (!GetNextEntry()) return kFALSE;
257 
258  // Event properties
259  Int_t eventId = fEvent->GetEventNr(); // Generator event ID
260  Int_t nPart = fEvent->GetNpa(); // Number of particles in event
261  Double_t impact = fEvent->GetB(); // Impact parameter
262  Double_t phi1 = fEvent->GetPhi(); // Event plane angle (generator) [rad]
263  fNofPrimaries = 0; // Number of registered primaries
264 
265  // Event rotation angle to be applied on the input momentum vectors
266  Double_t phi2 = 0.;
267  switch (fMode) {
268  case kStandard: phi2 = -1. * phi1; break;
269  case kRotateFixed: phi2 = fPhi; break;
270  case kNoRotation: phi2 = 0.; break;
271  case kReuseEvents: phi2 = gRandom->Uniform(0., 2 * TMath::Pi()); break;
272  default: phi2 = 0.; break;
273  } //? mode
274 
275  // Loop over particles in current event
276  for (Int_t iPart = 0; iPart < fEvent->GetNpa(); iPart++) {
277  UParticle* particle = fEvent->GetParticle(iPart);
278 
279  // Lorentz transformation into lab (target) frame
280  Double_t pz = fGammaCM * (particle->Pz() - fBetaCM * particle->E());
281  TVector3 momentum(particle->Px(), particle->Py(), pz);
282 
283  // Apply event rotation if needed
284  if (phi2 != 0.) momentum.RotateZ(phi2);
285 
286  // Normal particle
287  if (particle->GetPdg() < 1e9) // elementary particle
288  AddPrimary(primGen, particle->GetPdg(), momentum);
289  else
290  ProcessIon(primGen, particle->GetPdg(), momentum);
291 
292  } //# Particles
293 
294  // Set event information in the MC event header
295  FairMCEventHeader* eventHeader = primGen->GetEvent();
296  assert(eventHeader);
297  if (eventHeader->IsSet()) {
298  LOG(warning) << GetName() << ": Event header is already set; "
299  << " event info will not be stored";
300  } else {
301  eventHeader->SetEventID(eventId); // Generator event ID
302  eventHeader->SetB(fEvent->GetB()); // Impact parameter
303  eventHeader->SetRotZ(phi1 + phi2); // Event plane angle
304  eventHeader->MarkSet(kTRUE);
305  }
306 
307  // Info to screen
308  LOG(info) << GetName() << ": Event ID " << eventId << ", particles " << nPart
309  << ", primaries " << fNofPrimaries << ", b = " << impact
310  << " fm, phi (source) = " << phi1
311  << " rad , phi (generated) = " << phi2 << " rad";
312 
313  return kTRUE;
314 }
315 // --------------------------------------------------------------------------
316 
317 
318 // ----- Register ions to the run ---------------------------------------
320 
321  LOG(debug) << GetName() << ": Registering ions";
322  UParticle* particle {nullptr};
323  Int_t nIons {0};
324  fIonMap.clear();
325 
326  // --- Event loop
327  for (Int_t iEvent = 0; iEvent < fTree->GetEntries(); ++iEvent) {
328  fTree->GetEntry(iEvent);
329 
330  // --- Particle loop
331  for (Int_t iParticle = 0; iParticle < fEvent->GetNpa(); ++iParticle) {
332  particle = fEvent->GetParticle(iParticle);
333 
334  Int_t pdgCode = particle->GetPdg();
335  if (pdgCode > 1e9) { // ion
336 
337  // For ions the PDG code is +-10LZZZAAAI,
338  // where A = n_Lambda + n_proton + n_neutron, Z = n_proton, L = n_Lambda
339  // and I = 0 - ground state, I>0 - excitations
340  const Int_t nLambda = GetIonLambdas(pdgCode);
341  const Int_t charge = GetIonCharge(pdgCode);
342  const Int_t mass = GetIonMass(pdgCode);
343 
344  // Neutral ions are skipped (not present in G4IonTable)
345  if (charge == 0) continue;
346 
347  // Add ion to ion map
348  TString ionName = Form("Ion_%d_%d_%d", mass, charge, nLambda);
349  if (fIonMap.find(ionName) == fIonMap.end()) { // new ion
350  fIonMap[ionName] = new FairIon(ionName, charge, mass, charge);
351  nIons++;
352  LOG(debug) << GetName() << ": Adding new ion " << ionName << ", PDG "
353  << pdgCode << " (mass " << mass << ", charge " << charge
354  << ", nLambda " << nLambda << ") from event " << iEvent
355  << " at index " << iParticle;
356  } //? new ion
357 
358  } //? ion
359  } //# particles
360  } //# events
361 
362  FairRunSim* run = FairRunSim::Instance();
363  for (const auto& entry : fIonMap)
364  run->AddNewIon(entry.second);
365 
366  return fIonMap.size();
367 }
368 // ------------------------------------------------------------------------
369 
CbmUnigenGenerator::EMode
EMode
Mode enumerator.
Definition: CbmUnigenGenerator.h:48
CbmUnigenGenerator::ProcessIon
void ProcessIon(FairPrimaryGenerator *primGen, Int_t pdgCode, const TVector3 &momentum)
Treat a composite particle (ion)
Definition: CbmUnigenGenerator.cxx:216
UParticle::Pz
Double_t Pz() const
Definition: UParticle.h:79
UEvent::GetPhi
Double_t GetPhi() const
Definition: UEvent.h:33
CbmUnigenGenerator::GetIonLambdas
Int_t GetIonLambdas(Int_t pdgCode) const
Number of Lambdas in an ion.
Definition: CbmUnigenGenerator.h:145
CbmUnigenGenerator::fBetaCM
Double_t fBetaCM
CM velocity in the lab frame.
Definition: CbmUnigenGenerator.h:101
CbmUnigenGenerator::AddPrimary
void AddPrimary(FairPrimaryGenerator *primGen, Int_t pdgCode, const TVector3 &momentum)
Add a primary particle to the event generator.
Definition: CbmUnigenGenerator.cxx:89
UEvent::GetB
Double_t GetB() const
Definition: UEvent.h:32
UParticle
Definition: UParticle.h:10
UParticle::GetPdg
Int_t GetPdg() const
Definition: UParticle.h:69
CbmUnigenGenerator::fEvent
UEvent * fEvent
Current input event.
Definition: CbmUnigenGenerator.h:98
UEvent::GetNpa
Int_t GetNpa() const
Definition: UEvent.h:37
CbmUnigenGenerator::Init
Bool_t Init()
Initialisation.
Definition: CbmUnigenGenerator.cxx:138
CbmUnigenGenerator::fIonMap
std::map< TString, FairIon * > fIonMap
Map from ion name to FairIon.
Definition: CbmUnigenGenerator.h:103
CbmUnigenGenerator::fFile
TFile * fFile
Input ROOT file.
Definition: CbmUnigenGenerator.h:95
CbmUnigenGenerator
Generates input to transport simulation from files in Unigen format.
Definition: CbmUnigenGenerator.h:44
ClassImp
ClassImp(CbmUnigenGenerator)
CbmUnigenGenerator::fNofEvents
Int_t fNofEvents
Number of processed events.
Definition: CbmUnigenGenerator.h:100
CbmUnigenGenerator::kStandard
@ kStandard
Rotate events to zero event plane (default)
Definition: CbmUnigenGenerator.h:49
URun::GetNNSqrtS
Double_t GetNNSqrtS()
Definition: URun.cxx:153
UParticle.h
CbmUnigenGenerator::kRotateFixed
@ kRotateFixed
Rotate events around z by a fixed angle.
Definition: CbmUnigenGenerator.h:51
UEvent
Definition: UEvent.h:12
UEvent.h
URun::GetPProj
Double_t GetPProj() const
Definition: URun.h:54
CbmUnigenGenerator::ReadEvent
virtual Bool_t ReadEvent(FairPrimaryGenerator *primGen)
Read one event from the input file.
Definition: CbmUnigenGenerator.cxx:249
CbmUnigenGenerator::fNofPrimaries
Int_t fNofPrimaries
Number of primaries registered in current event.
Definition: CbmUnigenGenerator.h:99
CbmUnigenGenerator::fFileName
TString fFileName
Input file name.
Definition: CbmUnigenGenerator.h:91
CbmUnigenGenerator::kReuseEvents
@ kReuseEvents
Reuse events if more are requested than present.
Definition: CbmUnigenGenerator.h:52
UEvent::GetParticle
UParticle * GetParticle(Int_t index) const
Definition: UEvent.cxx:92
CbmUnigenGenerator::RegisterIons
Int_t RegisterIons()
Register ions to the simulation.
Definition: CbmUnigenGenerator.cxx:319
CbmUnigenGenerator::kPdgLambda
static const Int_t kPdgLambda
Decomposition of ion PDG code.
Definition: CbmUnigenGenerator.h:109
CbmUnigenGenerator::fMode
EMode fMode
Rotation mode.
Definition: CbmUnigenGenerator.h:92
xMath::Pi
double Pi()
Definition: xMath.h:5
CbmUnigenGenerator::fPhi
Double_t fPhi
Event plane rotation angle.
Definition: CbmUnigenGenerator.h:93
CbmUnigenGenerator.h
CbmUnigenGenerator::kNoRotation
@ kNoRotation
No event rotation.
Definition: CbmUnigenGenerator.h:50
UParticle::E
Double_t E() const
Definition: UParticle.h:80
CbmUnigenGenerator::fGammaCM
Double_t fGammaCM
Gamma factor of CM in lab frame.
Definition: CbmUnigenGenerator.h:102
URun
Definition: URun.h:8
CbmUnigenGenerator::GetNextEntry
Bool_t GetNextEntry()
Get next entry from input tree.
Definition: CbmUnigenGenerator.cxx:110
CbmUnigenGenerator::fTree
TTree * fTree
Input ROOT tree.
Definition: CbmUnigenGenerator.h:96
CbmUnigenGenerator::CbmUnigenGenerator
CbmUnigenGenerator(const char *fileName="", EMode mode=kStandard)
Default constructor.
Definition: CbmUnigenGenerator.cxx:27
UParticle::Px
Double_t Px() const
Definition: UParticle.h:77
URun::GetPTarg
Double_t GetPTarg() const
Definition: URun.h:57
CbmUnigenGenerator::fCurrentEntry
Int_t fCurrentEntry
Current entry number.
Definition: CbmUnigenGenerator.h:97
URun.h
CbmUnigenGenerator::CloseInput
void CloseInput()
Close the input file.
Definition: CbmUnigenGenerator.cxx:100
CbmUnigenGenerator::fIsInit
Bool_t fIsInit
Flag whether generator is initialised.
Definition: CbmUnigenGenerator.h:94
CbmUnigenGenerator::~CbmUnigenGenerator
virtual ~CbmUnigenGenerator()
Destructor.
Definition: CbmUnigenGenerator.cxx:80
URun::Print
void Print(Option_t *="") const
Definition: URun.cxx:84
CbmUnigenGenerator::GetIonCharge
Int_t GetIonCharge(Int_t pdgCode) const
Charge number of an ion.
Definition: CbmUnigenGenerator.h:134
CbmUnigenGenerator::GetIonMass
Int_t GetIonMass(Int_t pdgCode) const
Mass number of an ion.
Definition: CbmUnigenGenerator.h:156
UParticle::Py
Double_t Py() const
Definition: UParticle.h:78
UEvent::GetEventNr
Int_t GetEventNr() const
Definition: UEvent.h:31