-
Notifications
You must be signed in to change notification settings - Fork 496
Expand file tree
/
Copy pathPythia8Generator.cxx
More file actions
196 lines (176 loc) · 6.78 KB
/
Pythia8Generator.cxx
File metadata and controls
196 lines (176 loc) · 6.78 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
// Copyright CERN and copyright holders of ALICE O2. This software is
// distributed under the terms of the GNU General Public License v3 (GPL
// Version 3), copied verbatim in the file "COPYING".
//
// See http://alice-o2.web.cern.ch/license for full licensing information.
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
/********************************************************************************
* Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH *
* *
* This software is distributed under the terms of the *
* GNU Lesser General Public Licence version 3 (LGPL) version 3, *
* copied verbatim in the file "LICENSE" *
********************************************************************************/
// -------------------------------------------------------------------------
// ----- M. Al-Turany June 2014 -----
// -------------------------------------------------------------------------
#include <cmath>
#include "TROOT.h"
#include "Pythia8/Basics.h" // for RndmEngine
#include "Pythia8/Pythia.h"
#include "FairPrimaryGenerator.h"
#include "FairGenerator.h"
#include "TRandom.h" // for TRandom
#include "TRandom1.h" // for TRandom1
#include "TRandom3.h" // for TRandom3, gRandom
#include "Generators/Pythia8Generator.h"
using namespace Pythia8;
namespace o2
{
namespace eventgen
{
class PyTr1Rng : public RndmEngine
{
public:
PyTr1Rng() { mRng = std::make_unique<TRandom1>(gRandom->GetSeed()); };
~PyTr1Rng() override = default;
Double_t flat() override { return mRng->TRandom1::Rndm(); };
private:
std::unique_ptr<TRandom1> mRng; //!
};
class PyTr3Rng : public RndmEngine
{
public:
PyTr3Rng() { mRng = std::make_unique<TRandom3>(gRandom->GetSeed()); };
~PyTr3Rng() override = default;
Double_t flat() override { return mRng->TRandom3::Rndm(); };
private:
std::unique_ptr<TRandom3> mRng; //!
};
// ----- Default constructor -------------------------------------------
Pythia8Generator::Pythia8Generator()
{
mUseRandom1 = kFALSE;
mUseRandom3 = kTRUE;
mId = 2212; // proton
mMom = 400; // proton
mHNL = 0; // HNL if set to !=0, for example 9900014, only track
mPythia = std::make_unique<Pythia>();
}
// -------------------------------------------------------------------------
// ----- Default constructor -------------------------------------------
Bool_t Pythia8Generator::Init()
{
if (mUseRandom1) {
mRandomEngine = std::make_unique<PyTr1Rng>();
}
if (mUseRandom3) {
mRandomEngine = std::make_unique<PyTr3Rng>();
}
mPythia->setRndmEnginePtr(mRandomEngine.get());
/** commenting these lines
as they would override external settings **/
/**
cout<<"Beam Momentum "<<fMom<<endl;
// Set arguments in Settings database.
mPythia.settings.mode("Beams:idA", fId);
mPythia.settings.mode("Beams:idB", 2212);
mPythia.settings.mode("Beams:frameType", 3);
mPythia.settings.parm("Beams:pxA", 0.);
mPythia.settings.parm("Beams:pyA", 0.);
mPythia.settings.parm("Beams:pzA", fMom);
mPythia.settings.parm("Beams:pxB", 0.);
mPythia.settings.parm("Beams:pyB", 0.);
mPythia.settings.parm("Beams:pzB", 0.);
**/
mPythia->init();
return kTRUE;
}
// -------------------------------------------------------------------------
// ----- Destructor ----------------------------------------------------
Pythia8Generator::~Pythia8Generator() = default;
// -------------------------------------------------------------------------
// ----- Passing the event ---------------------------------------------
Bool_t Pythia8Generator::ReadEvent(FairPrimaryGenerator* cpg)
{
const double mm2cm = 0.1;
const double clight = 2.997924580e10; //cm/c
Int_t npart = 0;
while (npart == 0) {
mPythia->next();
for (int i = 0; i < mPythia->event.size(); i++) {
if (mPythia->event[i].isFinal()) {
// only send HNL decay products to G4
if (mHNL != 0) {
Int_t im = mPythia->event[i].mother1();
if (mPythia->event[im].id() == mHNL) {
// for the moment, hardcode 110m is maximum decay length
Double_t z = mPythia->event[i].zProd();
Double_t x = abs(mPythia->event[i].xProd());
Double_t y = abs(mPythia->event[i].yProd());
// cout<<"debug HNL decay pos "<<x<<" "<< y<<" "<< z <<endl;
if (z < 11000. && z > 7000. && x < 250. && y < 250.) {
npart++;
}
}
} else {
npart++;
}
};
};
// happens if a charm particle being produced which does decay without producing a HNL. Try another event.
// if (npart == 0){ mPythia->event.list();}
};
// cout<<"debug p8 event 0 " << mPythia->event[0].id()<< " "<< mPythia->event[1].id()<< " "
// << mPythia->event[2].id()<< " "<< npart <<endl;
for (Int_t ii = 0; ii < mPythia->event.size(); ii++) {
Bool_t wanttracking = true;
if (!mPythia->event[ii].isFinal())
wanttracking = false;
Double_t z = mPythia->event[ii].zProd();
Double_t x = mPythia->event[ii].xProd();
Double_t y = mPythia->event[ii].yProd();
Double_t pz = mPythia->event[ii].pz();
Double_t px = mPythia->event[ii].px();
Double_t py = mPythia->event[ii].py();
Double_t t = mPythia->event[ii].tProd();
x *= mm2cm;
y *= mm2cm;
z *= mm2cm;
t *= mm2cm / clight;
Double_t e = -9e9;
// handle system entry separately
if (!wanttracking)
e = mPythia->event[ii].e();
//
printf("AddTrack from Pythia8 %5d \n", (Int_t)mPythia->event[ii].id());
cpg->AddTrack((Int_t)mPythia->event[ii].id(), px, py, pz, x, y, z,
(Int_t)mPythia->event[ii].mother1(), wanttracking, e, t);
} // particle lloop
return kTRUE;
}
// -------------------------------------------------------------------------
void Pythia8Generator::SetParameters(const char* par)
{
// Set Parameters
mPythia->readString(par);
cout << R"(mPythia->readString(")" << par << R"("))" << endl;
}
// -------------------------------------------------------------------------
void Pythia8Generator::Print()
{
mPythia->settings.listAll();
}
// -------------------------------------------------------------------------
void Pythia8Generator::GetPythiaInstance(int arg)
{
mPythia->particleData.list(arg);
cout << "canDecay " << mPythia->particleData.canDecay(arg) << " " << mPythia->particleData.mayDecay(arg) << endl;
}
// -------------------------------------------------------------------------
} // namespace eventgen
} // namespace o2
ClassImp(o2::eventgen::Pythia8Generator);