11 #include "FairLogger.h"
12 #include "FairMCEventHeader.h"
13 #include "FairPrimaryGenerator.h"
14 #include "FairRunSim.h"
28 : FairGenerator(
"UnigenGenerator",
"CBM generator")
42 LOG(debug) << GetName() <<
": Constructor";
44 LOG(fatal) << GetName() <<
": When choosing "
45 <<
"mode kRotateFixed, a rotation angle must be specified!";
46 if (fileName[0] ==
'\0')
return;
48 LOG(fatal) << GetName() <<
": Error in intialisation! Aborting...";
57 : FairGenerator(
"UnigenGenerator",
"CBM generator")
71 LOG(debug) << GetName() <<
": Constructor";
72 if (fileName[0] ==
'\0')
return;
74 LOG(fatal) << GetName() <<
": Error in intialisation! Aborting...";
81 LOG(debug) << GetName() <<
": Destructor";
91 const TVector3& momentum) {
93 pdgCode, momentum.Px(), momentum.Py(), momentum.Pz(), 0., 0., 0.);
102 LOG(INFO) << GetName() <<
": Closing input file " <<
fFileName;
114 LOG(info) << GetName()
115 <<
": End of input tree - restarting with first entry";
119 LOG(info) << GetName() <<
": End of input tree reached";
127 LOG(ERROR) << GetName() <<
": Error reading entry " <<
fCurrentEntry
128 <<
" (returns " << result <<
")!";
142 LOG(warning) << GetName() <<
": Already initialised!";
146 LOG(info) << GetName() <<
": Initialising...";
147 std::stringstream ss;
149 case kStandard: ss <<
" standard; rotate event plane to zero";
break;
151 ss <<
" rotate event plane by fixed angle " <<
fPhi;
153 case kNoRotation: ss <<
" no event plane rotation";
break;
155 ss <<
" re-use events if necessary; random event plane angle";
157 default: ss <<
"unkown";
break;
159 LOG(INFO) << GetName() <<
": Mode " << ss.str();
164 if (!
fFile->IsOpen()) {
165 LOG(error) << GetName() <<
": Could not open input file " <<
fFileName;
168 LOG(info) << GetName() <<
": Open input file " <<
fFileName;
172 if (run ==
nullptr) {
173 LOG(error) << GetName() <<
": No run description in input file!";
177 LOG(info) << GetName() <<
": Run description";
181 Double_t mProt = 0.938272;
184 Double_t eTarg = TMath::Sqrt(pProj * pProj + mProt * mProt);
185 Double_t eProj = TMath::Sqrt(pTarg * pTarg + mProt * mProt);
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: "
195 fTree =
dynamic_cast<TTree*
>(
fFile->Get(
"events"));
196 if (
fTree ==
nullptr) {
197 LOG(error) << GetName() <<
": No event tree in input file!";
202 LOG(info) << GetName() <<
": " <<
fTree->GetEntries()
203 <<
" events in input tree";
207 LOG(info) << GetName() <<
": " << nIons <<
" ions registered.";
218 const TVector3& momentum) {
220 assert(pdgCode >= 1e9);
224 Int_t ionPdg = pdgCode;
227 LOG(warn) << GetName() <<
": Replacing hypernucleus (PDG " << pdgCode
228 <<
") by PDG " << ionPdg;
238 TVector3 neutronMom = momentum * (1. / Double_t(mass));
239 for (Int_t iNeutron = 0; iNeutron < mass; iNeutron++)
241 LOG(warn) << GetName() <<
": Neutral fragment with PDG " << ionPdg
242 <<
" is split into " << mass <<
" neutrons";
268 case kStandard: phi2 = -1. * phi1;
break;
272 default: phi2 = 0.;
break;
276 for (Int_t iPart = 0; iPart <
fEvent->
GetNpa(); iPart++) {
281 TVector3 momentum(particle->
Px(), particle->
Py(), pz);
284 if (phi2 != 0.) momentum.RotateZ(phi2);
287 if (particle->
GetPdg() < 1e9)
295 FairMCEventHeader* eventHeader = primGen->GetEvent();
297 if (eventHeader->IsSet()) {
298 LOG(warning) << GetName() <<
": Event header is already set; "
299 <<
" event info will not be stored";
301 eventHeader->SetEventID(eventId);
303 eventHeader->SetRotZ(phi1 + phi2);
304 eventHeader->MarkSet(kTRUE);
308 LOG(info) << GetName() <<
": Event ID " << eventId <<
", particles " << nPart
310 <<
" fm, phi (source) = " << phi1
311 <<
" rad , phi (generated) = " << phi2 <<
" rad";
321 LOG(debug) << GetName() <<
": Registering ions";
327 for (Int_t iEvent = 0; iEvent <
fTree->GetEntries(); ++iEvent) {
328 fTree->GetEntry(iEvent);
331 for (Int_t iParticle = 0; iParticle <
fEvent->
GetNpa(); ++iParticle) {
334 Int_t pdgCode = particle->
GetPdg();
345 if (charge == 0)
continue;
348 TString ionName = Form(
"Ion_%d_%d_%d", mass, charge, nLambda);
350 fIonMap[ionName] =
new FairIon(ionName, charge, mass, charge);
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;
362 FairRunSim* run = FairRunSim::Instance();
363 for (
const auto& entry :
fIonMap)
364 run->AddNewIon(entry.second);