FairShip
Loading...
Searching...
No Matches
DPPythia8Generator.cxx
Go to the documentation of this file.
1// SPDX-License-Identifier: LGPL-3.0-or-later
2// SPDX-FileCopyrightText: Copyright CERN for the benefit of the SHiP
3// Collaboration
4
5#include "DPPythia8Generator.h"
6
7#include <cmath>
8
9#include "FairPrimaryGenerator.h"
10#include "Pythia8/Pythia.h"
11#include "TDatabasePDG.h" // for TDatabasePDG
12#include "TMath.h"
13#include "TROOT.h"
14#include "TSystem.h"
15const Double_t cm = 10.; // pythia units are mm
16const Double_t c_light = 2.99792458e+10; // speed of light in cm/sec (c_light
17 // = 2.99792458e+8 * m/s)
18const Int_t debug = 1;
19
20// ----- Default constructor -------------------------------------------
21DPPythia8Generator::DPPythia8Generator() {
22 // fHadDecay = false;
23 fId = 2212; // proton
24 fMom = 400; // proton
25 fDP = 9900015; // DP pdg code
26 fLmin = 5000. * cm; // mm minimum decay position z ROOT units !
27 fLmax = 12000. * cm; // mm maximum decay position z
28 fFDs = 7.7 / 10.4; // correction for Pythia6 to match measured Ds production
29 fpbrem = kFALSE;
30 fpbremPDF = 0;
31 fdy = kFALSE;
32 fDPminM = 0.5;
33 fInputFile = nullptr;
34 fnRetries = 0;
35 fnDPtot = 0;
36 fShipEventNr = 0;
37 fPythia = new Pythia8::Pythia();
38 // fPythiaHadDecay = new Pythia8::Pythia();
39}
40// -------------------------------------------------------------------------
41
42// ----- Default constructor -------------------------------------------
43Bool_t DPPythia8Generator::Init() {
44 if (fUseRandom1) fRandomEngine = std::make_shared<PyTr1Rng>();
45 if (fUseRandom3) fRandomEngine = std::make_shared<PyTr3Rng>();
46 fPythia->setRndmEnginePtr(fRandomEngine);
47 // fPythiaHadDecay->setRndmEnginePtr(fRandomEngine);
48 fn = 0;
49 if (fextFile) {
50 /*if (0 == strncmp("/eos",fextFile,4) ) {
51 if (0 == strncmp("/eos",fextFile,4) ) {
52 TString tmp = gSystem->Getenv("EOSSHIP");
53 tmp+=fextFile;
54 fInputFile = TFile::Open(tmp);
55 fLogger->Info(MESSAGE_ORIGIN,""Open external file with charm or beauty
56 hadrons on eos: %s",tmp->Data()); if (!fInputFile) {
57 fLogger->Fatal(MESSAGE_ORIGIN, "Error opening input file. You may have
58 forgotten to provide a krb5 token. Try kinit username@lxplus.cern.ch");
59 return kFALSE; }
60 }else{
61 fLogger->Info(MESSAGE_ORIGIN,"Open external file with charm or beauty
62 hadrons: %s",fextFile); fInputFile = new TFile(fextFile->c_str()); if
63 (!fInputFile) { fLogger->Fatal(MESSAGE_ORIGIN, "Error opening input file");
64 return kFALSE;
65 }
66 }
67 if (fInputFile->IsZombie()) {
68 fLogger->Fatal(MESSAGE_ORIGIN, "File is corrupted");
69 return kFALSE; }
70 fTree = (TTree *)fInputFile->Get("pythia6");
71 fNevents = fTree->GetEntries();
72 fn = firstEvent;
73 fTree->SetBranchAddress("id",&hid); // particle id
74 fTree->SetBranchAddress("px",&hpx); // momentum
75 fTree->SetBranchAddress("py",&hpy);
76 fTree->SetBranchAddress("pz",&hpz);
77 fTree->SetBranchAddress("E",&hE);
78 fTree->SetBranchAddress("M",&hM);
79 fTree->SetBranchAddress("mid",&mid); // mother
80 fTree->SetBranchAddress("mpx",&mpx); // momentum
81 fTree->SetBranchAddress("mpy",&mpy);
82 fTree->SetBranchAddress("mpz",&mpz);
83 fTree->SetBranchAddress("mE",&mE);
84 */
85 } else if (!fpbrem) {
86 if (debug) {
87 std::cout << "Beam Momentum " << fMom << std::endl;
88 }
89 fPythia->settings.mode("Beams:idA", fId);
90 fPythia->settings.mode("Beams:idB", 2212);
91 fPythia->settings.mode("Beams:frameType", 2);
92 fPythia->settings.parm("Beams:eA", fMom); // codespell:ignore parm
93 fPythia->settings.parm("Beams:eB", 0.); // codespell:ignore parm
94
95 if (fdy)
96 fPythia->settings.parm("PhaseSpace:mHatMin", // codespell:ignore parm
97 fDPminM);
98
99 } else {
100 if (!fpbremPDF) {
101 // std::cout << " Failed in retrieving dark photon PDF for production by
102 // proton bremstrahlung! Exiting..." << std::endl;
103 LOG(fatal) << "Failed in retrieving dark photon PDF for production by "
104 "proton bremstrahlung!";
105 return kFALSE;
106 }
107 }
108 /*if (fHadDecay) {
109 std::cout << " ******************************** " << std::endl
110 << " ** Initialise Pythia for e+e-->hadrons " << std::endl
111 << " ******************************** " << std::endl
112 << " Mass of the A: " << fPythia->particleData.m0(fDP) << " GeV"
113 << std::endl; fPythiaHadDecay->settings.mode("Beams:idA", 11);
114 fPythiaHadDecay->settings.mode("Beams:idB", -11);
115 fPythiaHadDecay->settings.mode("Beams:frameType", 1);
116 fPythiaHadDecay->settings.parm("Beams:eCM",10); // codespell:ignore parm
117 fPythiaHadDecay->readString("WeakSingleBoson:ffbar2ffbar(s:gm) = on");
118 //fPythiaHadDecay->readString("23:isResonance = true")
119 //force to hadrons
120 fPythiaHadDecay->readString("23:onMode = off");
121 fPythiaHadDecay->readString("23:onIfAny = 1 2 3 4 5");
122 }*/
123 TDatabasePDG* pdgBase = TDatabasePDG::Instance();
124 Double_t root_ctau = pdgBase->GetParticle(fDP)->Lifetime();
125 // fPythia->particleData.readString("4900023:useBreitWigner = false");
126 if (debug) {
127 std::cout << "Final particle parameters for PDGID " << fDP << ":"
128 << std::endl;
129 List(fDP);
130 }
131 if (debug) {
132 std::cout << "tau root PDG database " << root_ctau
133 << "[s] ctau root = " << root_ctau * 3e10 << "[cm]" << std::endl;
134 }
135 fctau = fPythia->particleData.tau0(fDP); //* 3.3333e-12
136 if (debug) {
137 std::cout << "ctau pythia " << fctau << "[mm]" << std::endl;
138 }
139 int initPass = fPythia->init();
140 if (debug) {
141 std::cout << "Pythia initialisation bool: " << initPass << std::endl;
142 }
143
144 if (!initPass) {
145 LOG(fatal) << "Pythia initialisation failed";
146 return kFALSE;
147 }
148
149 return kTRUE;
150 // if (fHadDecay) fPythiaHadDecay->init();
151 // return kTRUE;
152}
153// -------------------------------------------------------------------------
154
155// ----- Destructor ----------------------------------------------------
156DPPythia8Generator::~DPPythia8Generator() {}
157// -------------------------------------------------------------------------
158
159// ----- Passing the event ---------------------------------------------
160Bool_t DPPythia8Generator::ReadEvent(FairPrimaryGenerator* cpg) {
161 Double_t tp, tS, zp, xp, yp, zS, xS, yS, pz, px, py, e, w;
162 Double_t tm, zm, xm, ym, pmz, pmx, pmy, em;
163 // Double_t td,zd,xd,yd;
164 Int_t im;
165 // take DP decay of Pythia, move it to the SHiP decay region
166 // recalculate decay time
167 // weight event with exp(-t_ship/tau_DP) / exp(-t_pythia/tau_DP)
168
169 int iDP =
170 0; // index of the chosen DP, also ensures that at least 1 DP is produced
171 std::vector<int>
172 dec_chain; // pythia indices of the particles to be stored on the stack
173 std::vector<int> dpvec; // pythia indices of DP particles
174 do {
175 if (fextFile) {
176 // take charm or beauty hadron from external file
177 // correct for too much Ds produced by pythia6
178 /*bool x = true;
179 while(x){
180 if (fn==fNevents) {fLogger->Warning(MESSAGE_ORIGIN, "End of input file.
181 Rewind.");} fTree->GetEntry(fn%fNevents); fn++; if ( int(fabs(hid[0]) )
182 != 431){ x = false; } else { Double_t rnr = gRandom->Uniform(0,1); if(
183 rnr<fFDs ) { x = false; };
184 //std::cout<<"what is x "<<x<<" id "<<int(fabs(hid[0]))<<" rnr " << rnr
185 <<" "<< fFDs <<std::endl ;
186 }
187 }
188 fPythia->event.reset();
189 fPythia->event.append( (Int_t)hid[0], 1, 0, 0, hpx[0], hpy[0], hpz[0],
190 hE[0], hM[0], 0., 9. );
191 */
192 }
193 // bit for proton brem
194 if (fpbrem) {
195 fPythia->event.reset();
196 double dpmom = 0;
197 double thetain = 0;
198 fpbremPDF->GetRandom2(dpmom, thetain);
199 double dpm = fPythia->particleData.m0(fDP);
200 double dpe = sqrt(dpmom * dpmom + dpm * dpm);
201 double phiin = 2. * M_PI * gRandom->Rndm();
202
203 if (debug > 1) {
204 std::cout << " Adding DP gun with p " << dpmom << " m " << dpm << " e "
205 << dpe << " theta,phi " << thetain << "," << phiin
206 << std::endl
207 << std::flush;
208 }
209 fPythia->event.append(fDP, 1, 0, 0, dpmom * sin(thetain) * cos(phiin),
210 dpmom * sin(thetain) * sin(phiin),
211 dpmom * cos(thetain), dpe, dpm);
212 }
213
214 if (!fPythia->next()) LOG(fatal) << "fPythia->next() failed";
215
216 // fPythia->event.list();
217 for (int i = 0; i < fPythia->event.size(); i++) {
218 // find all DP
219 if (abs(fPythia->event[i].id()) == fDP) {
220 dpvec.push_back(i);
221 }
222 }
223 iDP = dpvec.size();
224 fnDPtot += iDP;
225 if (iDP == 0) {
226 // fLogger->Info(MESSAGE_ORIGIN,"Event without DP. Retry.");
227 // fPythia->event.list();
228 fnRetries +=
229 1; // can happen if phasespace does not allow meson to decay to DP
230 } else {
231 // for mesons, could have more than one ... but for DY prod, need to take
232 // the last one... int r = int( gRandom->Uniform(0,iDP) );
233 int r = iDP - 1;
234 // std::cout << " ----> debug 2 " << r << std::endl;
235 int i = dpvec[r];
236 // production vertex
237 zp = fPythia->event[i].zProd();
238 xp = fPythia->event[i].xProd();
239 yp = fPythia->event[i].yProd();
240 tp = fPythia->event[i].tProd();
241 // momentum
242 pz = fPythia->event[i].pz();
243 px = fPythia->event[i].px();
244 py = fPythia->event[i].py();
245 e = fPythia->event[i].e();
246 // old decay vertex
247 // Int_t ida =fPythia->event[i].daughter1();
248 std::cout << " Debug: decay product of A: "
249 << fPythia->event[fPythia->event[i].daughter1()].id() << " "
250 << fPythia->event[fPythia->event[i].daughter2()].id()
251 << std::endl;
252 // hadDecay = abs(fPythia->event[ida].id())!=11 &&
253 // abs(fPythia->event[ida].id())!=13 && abs(fPythia->event[ida].id())!=15;
254
255 // zd =fPythia->event[ida].zProd();
256 // xd =fPythia->event[ida].xProd();
257 // yd =fPythia->event[ida].yProd();
258 // td =fPythia->event[ida].tProd();
259 // new decay vertex
260 Double_t LS = gRandom->Uniform(fLmin, fLmax); // mm, G4 and Pythia8 units
261 Double_t p = TMath::Sqrt(px * px + py * py + pz * pz);
262 Double_t lam = LS / p;
263 xS = xp + lam * px;
264 yS = yp + lam * py;
265 zS = zp + lam * pz;
266 Double_t gam = e / TMath::Sqrt(e * e - p * p);
267 Double_t beta = p / e;
268 tS = tp + LS / beta; // units ? [mm/c] + [mm/beta] (beta is dimensionless
269 // speed, and c=1 here)
270 // if one would use [s], then tS = tp/(cm*c_light) +
271 // (LS/cm)/(beta*c_light) = tS/(cm*c_light) i.e. units look consistent
272 w = TMath::Exp(-LS / (beta * gam * fctau)) *
273 ((fLmax - fLmin) / (beta * gam * fctau));
274 im = (Int_t)fPythia->event[i].mother1();
275 zm = fPythia->event[im].zProd();
276 xm = fPythia->event[im].xProd();
277 ym = fPythia->event[im].yProd();
278 pmz = fPythia->event[im].pz();
279 pmx = fPythia->event[im].px();
280 pmy = fPythia->event[im].py();
281 em = fPythia->event[im].e();
282 tm = fPythia->event[im].tProd();
283 // foresee finite beam size
284 Double_t dx = 0;
285 Double_t dy = 0;
286 if (fsmearBeam > 0) {
287 Double_t test = fsmearBeam * fsmearBeam;
288 Double_t Rsq = test + 1.;
289 while (Rsq > test) {
290 dx = gRandom->Uniform(-1., 1.) * fsmearBeam;
291 dy = gRandom->Uniform(-1., 1.) * fsmearBeam;
292 Rsq = dx * dx + dy * dy;
293 }
294 }
295 if (fextFile) {
296 // take grand mother particle from input file, to know if primary or
297 // secondary production
298 // cpg->AddTrack((Int_t)mid[0],mpx[0],mpy[0],mpz[0],xm/cm+dx,ym/cm+dy,zm/cm,-1,false,mE[0],0.,1.);
299 // cpg->AddTrack((Int_t)fPythia->event[im].id(),pmx,pmy,pmz,xm/cm+dx,ym/cm+dy,zm/cm,0,false,em,tm/cm/c_light,w);
300 // // convert pythia's (x,y,z[mm], t[mm/c]) to ([cm], [s])
301 // cpg->AddTrack(fDP, px, py, pz, xp/cm+dx,yp/cm+dy,zp/cm,
302 // 1,false,e,tp/cm/c_light,w);
303 } else {
304 cpg->AddTrack(
305 (Int_t)fPythia->event[im].id(), pmx, pmy, pmz, xm / cm + dx,
306 ym / cm + dy, zm / cm, -1, false, em, tm / cm / c_light,
307 w); // convert pythia's (x,y,z[mm], t[mm/c]) to ([cm], [s])
308 cpg->AddTrack(fDP, px, py, pz, xp / cm + dx, yp / cm + dy, zp / cm, 0,
309 false, e, tp / cm / c_light, w);
310 }
311 // bookkeep the indices of stored particles
312 dec_chain.push_back(im);
313 dec_chain.push_back(i);
314 if (debug > 1)
315 std::cout << std::endl
316 << " insert mother id " << im
317 << " pdg=" << fPythia->event[im].id() << " pmz = " << pmz
318 << " [GeV], zm = " << zm << " [mm] tm = " << tm << " [mm/c]"
319 << std::endl;
320 if (debug > 1)
321 std::cout << " ----> insert DP id " << i << " pdg=" << fDP
322 << " pz = " << pz << " [GeV] zp = " << zp
323 << " [mm] tp = " << tp << " [mm/c]" << std::endl;
324 iDP = i;
325 }
326 } while (iDP ==
327 0); // ----------- avoid rare empty events w/o any DP's produced
328
329 if (fShipEventNr % 100 == 0) {
330 LOGF(info, "ship event %i / pythia event-nr %i", fShipEventNr, fn);
331 }
332 fShipEventNr += 1;
333 // fill a container with pythia indices of the DP decay chain
334 if (debug > 1) std::cout << "Filling daughter particles" << std::endl;
335 // if (!hadDecay){
336 for (int k = 0; k < fPythia->event.size(); k++) {
337 // if daughter of DP, copy
338 if (debug > 1)
339 std::cout << k << " pdg =" << fPythia->event[k].id() << " mum "
340 << fPythia->event[k].mother1() << std::endl;
341 im = fPythia->event[k].mother1();
342 while (im > 0) {
343 if (im == iDP) {
344 break;
345 } // pick the decay products of only 1 chosen DP
346 // if ( abs(fPythia->event[im].id())==fDP && im == iDP ){break;}
347 else {
348 im = fPythia->event[im].mother1();
349 }
350 }
351 if (im < 1) {
352 if (debug > 1) std::cout << "reject" << std::endl;
353 continue;
354 }
355 if (debug > 1) std::cout << "accept" << std::endl;
356 dec_chain.push_back(k);
357 }
358 //}
359 /*else {
360 //get decay from e+e- collision.....
361 fPythiaHadDecay->settings.parm("Beams:eCM",20); // codespell:ignore parm
362 fPythiaHadDecay->next();
363 for (int k=0; k<fPythiaHadDecay->event.size(); k++){
364 fPythia->event.append(
365 fPythiaHadDecay->event[k].id(),fPythiaHadDecay->event[k].status()
366 ,fPythiaHadDecay->event[k].mother1() , fPythiaHadDecay->event[k].mother2(),
367 fPythiaHadDecay->event[k].daughter1(), fPythiaHadDecay->event[k].daughter2(),
368 fPythiaHadDecay->event[k].col(), fPythiaHadDecay->event[k].acol(),
369 fPythiaHadDecay->event[k].px(), fPythiaHadDecay->event[k].py(),
370 fPythiaHadDecay->event[k].pz(), fPythiaHadDecay->event[k].e(),
371 fPythiaHadDecay->event[k].m(), 0., 9. ); dec_chain.push_back(
372 fPythia->event.size()-1 ); std::cout << " Adding decay product: " << k << " "
373 << fPythiaHadDecay->event[k].id() << " "
374 << fPythiaHadDecay->event[k].status() << " "
375 << fPythiaHadDecay->event[k].mother1() << " "
376 << fPythiaHadDecay->event[k].mother2() << " "
377 << std::endl;
378 }
379 }*/
380
381 // go over daughters and store them on the stack, starting from 2 to account
382 // for DP and its mother
383 for (std::vector<int>::iterator it = dec_chain.begin() + 2;
384 it != dec_chain.end(); ++it) {
385 // pythia index of the particle to store
386 int k = *it;
387 // find mother position on the output stack: impy -> im
388 int impy = fPythia->event[k].mother1();
389 std::vector<int>::iterator itm =
390 std::find(dec_chain.begin(), dec_chain.end(), impy);
391 im = -1; // safety
392 if (itm != dec_chain.end())
393 im = itm - dec_chain.begin(); // convert iterator into sequence number
394
395 Bool_t wanttracking = false;
396 if (fPythia->event[k].isFinal()) {
397 wanttracking = true;
398 }
399 pz = fPythia->event[k].pz();
400 px = fPythia->event[k].px();
401 py = fPythia->event[k].py();
402 e = fPythia->event[k].e();
403 if (fextFile) {
404 im += 1;
405 };
406 cpg->AddTrack((Int_t)fPythia->event[k].id(), px, py, pz, xS / cm, yS / cm,
407 zS / cm, im, wanttracking, e, tS / cm / c_light, w);
408 // std::cout <<k<< " insert pdg =" <<fPythia->event[k].id() << " pz = " <<
409 // pz << " [GeV] zS = " << zS << " [mm] tS = " << tS << "[mm/c]" <<
410 // std::endl;
411 }
412 return kTRUE;
413}
414// -------------------------------------------------------------------------
415void DPPythia8Generator::SetParameters(char* par) {
416 // Set Parameters
417 fPythia->readString(par);
418 if (debug) {
419 std::cout << "fPythia->readString(\"" << par << "\")" << std::endl;
420 }
421}
422
423// -------------------------------------------------------------------------
debug
const Int_t debug
Definition: DPPythia8Generator.cxx:18
cm
const Double_t cm
Definition: DPPythia8Generator.cxx:15
c_light
const Double_t c_light
Definition: DPPythia8Generator.cxx:16
cm
Double_t cm
Definition: ShipBellField.cxx:24
DPPythia8Generator::fPythia
Pythia8::Pythia * fPythia
Definition: DPPythia8Generator.h:89
DPPythia8Generator::ReadEvent
Bool_t ReadEvent(FairPrimaryGenerator *) override
Definition: DPPythia8Generator.cxx:160
DPPythia8Generator::fRandomEngine
std::shared_ptr< Pythia8::RndmEngine > fRandomEngine
Definition: DPPythia8Generator.h:92
DPPythia8Generator::Init
Bool_t Init() override
Definition: DPPythia8Generator.cxx:43
SHiP::Generator::fextFile
std::optional< std::string > fextFile
Definition: Generator.h:47