ERDecayEXP1803.cxx

Ivan Muzalevsky, 03/09/2018 05:28 PM

Download (9.5 KB)

 
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 "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
      // std::cout << "Start reation in target. Defined pos: " << fTargetReactZ << ", current pos: " << curPos.Z() << endl;
113
      // 6He + 2H → 3He + 5H
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
        // cout  << "exc 123 " << exc << endl;
143
        // LOG(DEBUG) << "Ion H5 mass in reaction " << mass5H << FairLogger::endl;
144

    
145
        reactMasses[0] = f3He->Mass();
146
        // cout << "H3 mass " << reactMasses[0] << endl;
147
        reactMasses[1] = mass5H;
148
        // cout << "mass5H mass " << reactMasses[1] << endl;
149
       
150
        // cout << "reaction decay momentum" << lvReaction.Px() << " " << lvReaction.Py() << " " << lvReaction.Pz() << endl; 
151
        if (fReactionPhaseSpace->SetDecay(lvReaction, 2, reactMasses) == kFALSE){
152
          // cout << " not work " << endl;
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
      //5H → f3H + n +n.
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)