1
|
|
2
|
|
3
|
|
4
|
|
5
|
|
6
|
|
7
|
|
8
|
|
9
|
#include "ERDecayEXP1803.h"
|
10
|
|
11
|
#include <iostream>
|
12
|
|
13
|
#include "TVirtualMC.h"
|
14
|
#include "TLorentzVector.h"
|
15
|
#include "TMCProcess.h"
|
16
|
#include "TRandom.h"
|
17
|
|
18
|
#include "FairRunSim.h"
|
19
|
#include "FairLogger.h"
|
20
|
|
21
|
#include "ERDecayMCEventHeader.h"
|
22
|
#include "ERMCEventHeader.h"
|
23
|
|
24
|
using namespace std;
|
25
|
|
26
|
ERDecayEXP1803::ERDecayEXP1803():
|
27
|
ERDecay("EXP1803"),
|
28
|
fDecayFinish(kFALSE),
|
29
|
fTargetReactZ(0.),
|
30
|
fMinStep(0.01),
|
31
|
f6He(NULL),
|
32
|
f2H (NULL),
|
33
|
f3He(NULL),
|
34
|
f5H (NULL),
|
35
|
f3H (NULL),
|
36
|
fn (NULL),
|
37
|
fIon3He(NULL),
|
38
|
f5HMass(0.),
|
39
|
fIs5HUserMassSet(false),
|
40
|
fIs5HExcitationSet(false)
|
41
|
{
|
42
|
fRnd = new TRandom3();
|
43
|
fReactionPhaseSpace = new TGenPhaseSpace();
|
44
|
fDecayPhaseSpace = new TGenPhaseSpace();
|
45
|
FairRunSim* run = FairRunSim::Instance();
|
46
|
fUnstableIon5H = new FairIon("5H", 1, 5, 1);
|
47
|
fIon3He = new FairIon("3He", 2, 3, 2);
|
48
|
run->AddNewIon(fUnstableIon5H);
|
49
|
run->AddNewIon(fIon3He);
|
50
|
}
|
51
|
|
52
|
ERDecayEXP1803::~ERDecayEXP1803() {
|
53
|
}
|
54
|
|
55
|
void ERDecayEXP1803::SetH5Exitation(Double_t excMean, Double_t fwhm, Double_t distibWeight) {
|
56
|
f5HExcitationMean.push_back(excMean);
|
57
|
f5HExcitationSigma.push_back(fwhm / 2.355);
|
58
|
if (!fIs5HExcitationSet) {
|
59
|
f5HExcitationWeight.push_back(distibWeight);
|
60
|
fIs5HExcitationSet = true;
|
61
|
return ;
|
62
|
}
|
63
|
f5HExcitationWeight.push_back(f5HExcitationWeight.back() + distibWeight);
|
64
|
}
|
65
|
|
66
|
Bool_t ERDecayEXP1803::Init(){
|
67
|
f6He = TDatabasePDG::Instance()->GetParticle("6He");
|
68
|
if ( ! f6He ) {
|
69
|
std::cerr << "-W- ERDecayEXP1803: Ion 6He not found in database!" << endl;
|
70
|
return kFALSE;
|
71
|
}
|
72
|
|
73
|
f2H = TDatabasePDG::Instance()->GetParticle("Deuteron");
|
74
|
if ( ! f2H ) {
|
75
|
std::cerr << "-W- ERDecayEXP1803: Ion Deuteron not found in database!" << endl;
|
76
|
return kFALSE;
|
77
|
}
|
78
|
|
79
|
f5H = TDatabasePDG::Instance()->GetParticle("5H");
|
80
|
if ( ! f5H ) {
|
81
|
std::cerr << "-W- ERDecayEXP1803: Ion 5H not found in database!" << endl;
|
82
|
return kFALSE;
|
83
|
}
|
84
|
|
85
|
f3He = TDatabasePDG::Instance()->GetParticle(fIon3He->GetName());
|
86
|
if ( ! f3He ) {
|
87
|
std::cerr << "-W- ERDecayEXP1803: Ion 3He not found in database!" << endl;
|
88
|
return kFALSE;
|
89
|
}
|
90
|
|
91
|
f3H = TDatabasePDG::Instance()->GetParticle("Triton");
|
92
|
if ( ! f3H ) {
|
93
|
std::cerr << "-W- ERDecayEXP1803: Ion Triton not found in database!" << endl;
|
94
|
return kFALSE;
|
95
|
}
|
96
|
|
97
|
fn = TDatabasePDG::Instance()->GetParticle("neutron");
|
98
|
if ( ! fn ) {
|
99
|
std::cerr << "-W- ERDecayEXP1803: Particle neutron not found in database!" << endl;
|
100
|
return kFALSE;
|
101
|
}
|
102
|
|
103
|
return kTRUE;
|
104
|
}
|
105
|
|
106
|
Bool_t ERDecayEXP1803::Stepping() {
|
107
|
if(!fDecayFinish && gMC->TrackPid() == 1000020060 && TString(gMC->CurrentVolName()).Contains(fVolumeName)){
|
108
|
gMC->SetMaxStep(fMinStep);
|
109
|
TLorentzVector curPos;
|
110
|
gMC->TrackPosition(curPos);
|
111
|
if (curPos.Z() > fTargetReactZ){
|
112
|
|
113
|
|
114
|
TLorentzVector lv6He;
|
115
|
gMC->TrackMomentum(lv6He);
|
116
|
cout << "He6 momentum " << lv6He.Px() << " " << lv6He.Py() << " " << lv6He.Pz() << endl;
|
117
|
|
118
|
TLorentzVector lv2H(0., 0., 0., f2H->Mass());
|
119
|
cout << "H2 momentum " << lv2H.Px() << " " << lv2H.Py() << " " << lv2H.Pz() << endl;
|
120
|
TLorentzVector lvReaction;
|
121
|
lvReaction = lv6He + lv2H;
|
122
|
cout << "lvReaction momentum" << lvReaction.Px() << " " << lvReaction.Py() << " " << lvReaction.Pz() << endl;
|
123
|
|
124
|
Double_t mass5H = (fIs5HUserMassSet) ? f5HMass : f5H->Mass();
|
125
|
Double_t exc = 0;
|
126
|
Double_t randWeight,reactMasses[2];
|
127
|
Int_t distribNum;
|
128
|
Bool_t threshold1 = kFALSE;
|
129
|
while (1) {
|
130
|
if (fIs5HExcitationSet) {
|
131
|
randWeight = gRandom->Uniform(0., f5HExcitationWeight.back());
|
132
|
distribNum = 0;
|
133
|
for (; distribNum < f5HExcitationWeight.size(); distribNum++) {
|
134
|
if (randWeight < f5HExcitationWeight[distribNum]) {
|
135
|
break;
|
136
|
}
|
137
|
}
|
138
|
exc = gRandom->Gaus(f5HExcitationMean[distribNum], f5HExcitationSigma[distribNum]);
|
139
|
fUnstableIon5H->SetExcEnergy(exc);
|
140
|
}
|
141
|
mass5H += exc;
|
142
|
|
143
|
|
144
|
|
145
|
reactMasses[0] = f3He->Mass();
|
146
|
|
147
|
reactMasses[1] = mass5H;
|
148
|
|
149
|
|
150
|
|
151
|
if (fReactionPhaseSpace->SetDecay(lvReaction, 2, reactMasses) == kFALSE){
|
152
|
|
153
|
mass5H -= exc;
|
154
|
}
|
155
|
else break;
|
156
|
}
|
157
|
fReactionPhaseSpace->Generate();
|
158
|
|
159
|
TLorentzVector *lv3He = fReactionPhaseSpace->GetDecay(0);
|
160
|
cout << "h3 decay momentum " << lv3He->Px() << " " << lv3He->Py() << " " << lv3He->Pz() << endl;
|
161
|
TLorentzVector *lv5H = fReactionPhaseSpace->GetDecay(1);
|
162
|
cout << "H5 decay momentum " << lv5H->Px() << " " << lv5H->Py() << " " << lv5H->Pz() << endl;
|
163
|
|
164
|
Double_t decayMasses[3];
|
165
|
decayMasses[0] = f3H->Mass();
|
166
|
decayMasses[1] = fn->Mass();
|
167
|
decayMasses[2] = fn->Mass();
|
168
|
|
169
|
fDecayPhaseSpace->SetDecay(*lv5H, 3, decayMasses);
|
170
|
fDecayPhaseSpace->Generate();
|
171
|
|
172
|
TLorentzVector *lv3H = fDecayPhaseSpace->GetDecay(0);
|
173
|
TLorentzVector *lvn1 = fDecayPhaseSpace->GetDecay(1);
|
174
|
TLorentzVector *lvn2 = fDecayPhaseSpace->GetDecay(2);
|
175
|
|
176
|
|
177
|
Int_t He6TrackNb, H5TrackNb, He3TrackNb, H3TrackNb, n1TrackNb, n2TrackNb;
|
178
|
|
179
|
He6TrackNb = gMC->GetStack()->GetCurrentTrackNumber();
|
180
|
|
181
|
cout << "h5 pdg " << f5H->PdgCode() << endl;
|
182
|
gMC->GetStack()->PushTrack(1, He6TrackNb, f5H->PdgCode(),
|
183
|
lv5H->Px(),lv5H->Py(),lv5H->Pz(),
|
184
|
lv5H->E(), curPos.X(), curPos.Y(), curPos.Z(),
|
185
|
gMC->TrackTime(), 0., 0., 0.,
|
186
|
kPDecay, H5TrackNb, mass5H, 0);
|
187
|
gMC->GetStack()->PushTrack(1, He6TrackNb, f3He->PdgCode(),
|
188
|
lv3He->Px(),lv3He->Py(),lv3He->Pz(),
|
189
|
lv3He->E(), curPos.X(), curPos.Y(), curPos.Z(),
|
190
|
gMC->TrackTime(), 0., 0., 0.,
|
191
|
kPDecay, He3TrackNb, f3He->Mass(), 0);
|
192
|
gMC->GetStack()->PushTrack(1, He6TrackNb, f3H->PdgCode(),
|
193
|
lv3H->Px(),lv3H->Py(),lv3H->Pz(),
|
194
|
lv3H->E(), curPos.X(), curPos.Y(), curPos.Z(),
|
195
|
gMC->TrackTime(), 0., 0., 0.,
|
196
|
kPDecay, H3TrackNb, f3H->Mass(), 0);
|
197
|
gMC->GetStack()->PushTrack(1, He6TrackNb, fn->PdgCode(),
|
198
|
lvn1->Px(),lvn1->Py(),lvn1->Pz(),
|
199
|
lvn1->E(), curPos.X(), curPos.Y(), curPos.Z(),
|
200
|
gMC->TrackTime(), 0., 0., 0.,
|
201
|
kPDecay, n1TrackNb, fn->Mass(), 0);
|
202
|
gMC->GetStack()->PushTrack(1, He6TrackNb, fn->PdgCode(),
|
203
|
lvn2->Px(),lvn2->Py(),lvn2->Pz(),
|
204
|
lvn2->E(), curPos.X(), curPos.Y(), curPos.Z(),
|
205
|
gMC->TrackTime(), 0., 0., 0.,
|
206
|
kPDecay, n2TrackNb, fn->Mass(), 0);
|
207
|
LOG(INFO) << He6TrackNb << " " << H5TrackNb << " " << He3TrackNb << FairLogger::endl;
|
208
|
gMC->StopTrack();
|
209
|
fDecayFinish = kTRUE;
|
210
|
gMC->SetMaxStep(100.);
|
211
|
|
212
|
FairRunSim* run = FairRunSim::Instance();
|
213
|
if (TString(run->GetMCEventHeader()->ClassName()).Contains("ERDecayMCEventHeader")){
|
214
|
ERDecayMCEventHeader* header = (ERDecayMCEventHeader*)run->GetMCEventHeader();
|
215
|
header->SetDecayPos(curPos.Vect());
|
216
|
header->SetInputIon(He6TrackNb);
|
217
|
header->AddOutputParticle(H5TrackNb);
|
218
|
header->AddOutputParticle(He3TrackNb);
|
219
|
header->AddOutputParticle(H3TrackNb);
|
220
|
header->AddOutputParticle(n1TrackNb);
|
221
|
header->AddOutputParticle(n2TrackNb);
|
222
|
}
|
223
|
}
|
224
|
}
|
225
|
return kTRUE;
|
226
|
}
|
227
|
|
228
|
void ERDecayEXP1803::BeginEvent() {
|
229
|
fDecayFinish = kFALSE;
|
230
|
fTargetReactZ = fRnd->Uniform(-fTargetThickness / 2, fTargetThickness / 2);
|
231
|
FairRunSim* run = FairRunSim::Instance();
|
232
|
}
|
233
|
|
234
|
void ERDecayEXP1803::FinishEvent() {
|
235
|
}
|
236
|
|
237
|
ClassImp(ERDecayEXP1803)
|