]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TFluka/TFlukaGeo.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Modifications for Mac
[u/mrichter/AliRoot.git] / TFluka / TFlukaGeo.cxx
CommitLineData
8495a208 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
18#include <Riostream.h>
19
e4b9dd05 20#include "AliModule.h"
21#include "AliRun.h"
8495a208 22#include "TClonesArray.h"
23#include "TFlukaGeo.h"
24#include "TCallf77.h" //For the fortran calls
25#include "Fdblprc.h" //(DBLPRC) fluka common
26#include "Fepisor.h" //(EPISOR) fluka common
27#include "Ffinuc.h" //(FINUC) fluka common
28#include "Fiounit.h" //(IOUNIT) fluka common
29#include "Fpaprop.h" //(PAPROP) fluka common
30#include "Fpart.h" //(PART) fluka common
31#include "Ftrackr.h" //(TRACKR) fluka common
32#include "Fpaprop.h" //(PAPROP) fluka common
33#include "Ffheavy.h" //(FHEAVY) fluka common
34
35#include "TVirtualMC.h"
36#include "TGeoManager.h"
c7a186a7 37#include "TGeoMaterial.h"
38#include "TGeoMedium.h"
8495a208 39#include "TFlukaMCGeometry.h"
c7a186a7 40#include "TFlukaCerenkov.h"
8495a208 41#include "TLorentzVector.h"
42
43// Fluka methods that may be needed.
44#ifndef WIN32
45# define flukam flukam_
46# define fluka_openinp fluka_openinp_
47# define fluka_closeinp fluka_closeinp_
48# define mcihad mcihad_
49# define mpdgha mpdgha_
50#else
51# define flukam FLUKAM
52# define fluka_openinp FLUKA_OPENINP
53# define fluka_closeinp FLUKA_CLOSEINP
54# define mcihad MCIHAD
55# define mpdgha MPDGHA
56#endif
57
58extern "C"
59{
60 //
61 // Prototypes for FLUKA functions
62 //
63 void type_of_call flukam(const int&);
64 void type_of_call fluka_openinp(const int&, DEFCHARA);
65 void type_of_call fluka_closeinp(const int&);
66 int type_of_call mcihad(const int&);
67 int type_of_call mpdgha(const int&);
68}
69
70//
71// Class implementation for ROOT
72//
73ClassImp(TFluka)
74
75//
76//----------------------------------------------------------------------------
77// TFluka constructors and destructors.
78//______________________________________________________________________________
79TFluka::TFluka()
80 :TVirtualMC(),
81 fVerbosityLevel(0),
82 sInputFileName("")
83{
84 //
85 // Default constructor
86 //
87 fNVolumes = 0;
88 fCurrentFlukaRegion = -1;
89 fGeom = 0;
c7a186a7 90 fMaterials = 0;
8495a208 91}
92
93//______________________________________________________________________________
94TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupported)
95 :TVirtualMC("TFluka",title, isRootGeometrySupported),
96 fVerbosityLevel(verbosity),
97 sInputFileName(""),
98 fTrackIsEntering(0),
fbf08100 99 fTrackIsExiting(0),
100 fTrackIsNew(0)
8495a208 101{
102 // create geometry interface
103 if (fVerbosityLevel >=3)
104 cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
105
106 fNVolumes = 0;
107 fCurrentFlukaRegion = -1;
108 fGeom = new TFlukaMCGeometry("geom", "ALICE geometry");
c7a186a7 109 fMaterials = 0;
8495a208 110}
111
112//______________________________________________________________________________
113TFluka::~TFluka() {
114 if (fVerbosityLevel >=3)
115 cout << "==> TFluka::~TFluka() destructor called." << endl;
116
117 delete fGeom;
118
119 if (fVerbosityLevel >=3)
120 cout << "<== TFluka::~TFluka() destructor called." << endl;
121}
122
123//
124//______________________________________________________________________________
125// TFluka control methods
126//______________________________________________________________________________
127void TFluka::Init() {
128
c7a186a7 129 if (fVerbosityLevel >=3)
130 cout << "==> TFluka::Init() called." << endl;
131
132 if (!gGeoManager) new TGeoManager("geom", "FLUKA geometry");
133 fApplication->ConstructGeometry();
134 TGeoVolume *top = (TGeoVolume*)gGeoManager->GetListOfVolumes()->First();
135 gGeoManager->SetTopVolume(top);
136 gGeoManager->CloseGeometry("di");
137 gGeoManager->DefaultColors(); // to be removed
138 fNVolumes = fGeom->NofVolumes();
139 printf("== Number of volumes: %i\n ==", fNVolumes);
140 fGeom->CreateFlukaMatFile("flukaMat.inp");
141 cout << "\t* InitPhysics() - Prepare input file to be called" << endl;
142 // now we have TGeo geometry created and we have to patch alice.inp
143 // with the material mapping file FlukaMat.inp
e4b9dd05 144}
145
146
147//______________________________________________________________________________
148void TFluka::FinishGeometry() {
149//
150// Build-up table with region to medium correspondance
151//
152 if (fVerbosityLevel >=3)
153 cout << "==> TFluka::FinishGeometry() called." << endl;
154
155 printf("----FinishGeometry - nothing to do with TGeo\n");
156
157 if (fVerbosityLevel >=3)
158 cout << "<== TFluka::FinishGeometry() called." << endl;
159}
160
161//______________________________________________________________________________
162void TFluka::BuildPhysics() {
163 if (fVerbosityLevel >=3)
164 cout << "==> TFluka::BuildPhysics() called." << endl;
c7a186a7 165 InitPhysics(); // prepare input file with the current physics settings
8495a208 166 cout << "\t* InitPhysics() - Prepare input file was called" << endl;
c7a186a7 167
168 if (fVerbosityLevel >=2)
169 cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F')
170 << ") in fluka..." << endl;
171 GLOBAL.lfdrtr = true;
172
173 if (fVerbosityLevel >=2)
174 cout << "\t* Opening file " << sInputFileName << endl;
175 const char* fname = sInputFileName;
176 fluka_openinp(lunin, PASSCHARA(fname));
177
178 if (fVerbosityLevel >=2)
179 cout << "\t* Calling flukam..." << endl;
180 flukam(1);
181
182 if (fVerbosityLevel >=2)
183 cout << "\t* Closing file " << sInputFileName << endl;
184 fluka_closeinp(lunin);
185
186 FinishGeometry();
187
188 if (fVerbosityLevel >=3)
189 cout << "<== TFluka::Init() called." << endl;
e4b9dd05 190
191
192 if (fVerbosityLevel >=3)
193 cout << "<== TFluka::BuildPhysics() called." << endl;
8495a208 194}
195
196//______________________________________________________________________________
197void TFluka::ProcessEvent() {
198 if (fVerbosityLevel >=3)
199 cout << "==> TFluka::ProcessEvent() called." << endl;
200 fApplication->GeneratePrimaries();
201 EPISOR.lsouit = true;
202 flukam(1);
203 if (fVerbosityLevel >=3)
204 cout << "<== TFluka::ProcessEvent() called." << endl;
205}
206
207//______________________________________________________________________________
208void TFluka::ProcessRun(Int_t nevent) {
209 if (fVerbosityLevel >=3)
210 cout << "==> TFluka::ProcessRun(" << nevent << ") called."
211 << endl;
212
213 if (fVerbosityLevel >=2) {
214 cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
215 cout << "\t* Calling flukam again..." << endl;
216 }
217 fApplication->InitGeometry();
218 fApplication->BeginEvent();
219 ProcessEvent();
220 fApplication->FinishEvent();
221 if (fVerbosityLevel >=3)
222 cout << "<== TFluka::ProcessRun(" << nevent << ") called."
223 << endl;
224
225}
226
227//_____________________________________________________________________________
228// methods for building/management of geometry
229
230// functions from GCONS
231//____________________________________________________________________________
232void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,
233 Float_t &dens, Float_t &radl, Float_t &absl,
234 Float_t* ubuf, Int_t& nbuf) {
235//
236 fGeom->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf);
237}
238
239//______________________________________________________________________________
240void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,
241 Double_t &dens, Double_t &radl, Double_t &absl,
242 Double_t* ubuf, Int_t& nbuf) {
243//
244 fGeom->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf);
245}
246
247// detector composition
248//______________________________________________________________________________
249void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
250 Double_t z, Double_t dens, Double_t radl, Double_t absl,
251 Float_t* buf, Int_t nwbuf) {
252//
253 fGeom->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf);
254}
255
256//______________________________________________________________________________
257void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
258 Double_t z, Double_t dens, Double_t radl, Double_t absl,
259 Double_t* buf, Int_t nwbuf) {
260//
261 fGeom->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf);
262}
263
264//______________________________________________________________________________
265void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a,
266 Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
267//
268 fGeom->Mixture(kmat, name, a, z, dens, nlmat, wmat);
269}
270
271//______________________________________________________________________________
272void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a,
273 Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
274//
275 fGeom->Mixture(kmat, name, a, z, dens, nlmat, wmat);
276}
277
278//______________________________________________________________________________
279void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
280 Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
281 Double_t stemax, Double_t deemax, Double_t epsil,
282 Double_t stmin, Float_t* ubuf, Int_t nbuf) {
283 //
284 fGeom->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
285 epsil, stmin, ubuf, nbuf);
286}
287
288//______________________________________________________________________________
289void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
290 Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
291 Double_t stemax, Double_t deemax, Double_t epsil,
292 Double_t stmin, Double_t* ubuf, Int_t nbuf) {
293 //
294 fGeom->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
295 epsil, stmin, ubuf, nbuf);
296}
297
298//______________________________________________________________________________
299void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX,
300 Double_t thetaY, Double_t phiY, Double_t thetaZ,
301 Double_t phiZ) {
302//
303 fGeom->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
304}
305
306//______________________________________________________________________________
e4b9dd05 307void TFluka::Gstpar(Int_t itmed, const char* param, Double_t parval) {
8495a208 308//
309// Is it needed with TGeo ??? - to clear-up
e4b9dd05 310//
311
312 printf("Gstpar called with %6d %5s %12.4e %6d\n", itmed, param, parval, fGeom->GetFlukaMaterial(itmed));
313
314 Bool_t process = kFALSE;
315 if (strncmp(param, "DCAY", 4) == 0 ||
316 strncmp(param, "PAIR", 4) == 0 ||
317 strncmp(param, "COMP", 4) == 0 ||
318 strncmp(param, "PHOT", 4) == 0 ||
319 strncmp(param, "PFIS", 4) == 0 ||
320 strncmp(param, "DRAY", 4) == 0 ||
321 strncmp(param, "ANNI", 4) == 0 ||
322 strncmp(param, "BREM", 4) == 0 ||
323 strncmp(param, "MUNU", 4) == 0 ||
324 strncmp(param, "CKOV", 4) == 0 ||
325 strncmp(param, "HADR", 4) == 0 ||
326 strncmp(param, "LOSS", 4) == 0 ||
327 strncmp(param, "MULS", 4) == 0 ||
328 strncmp(param, "RAYL", 4) == 0)
329 {
330 process = kTRUE;
331 }
332 if (process) {
333 SetProcess(param, Int_t (parval), fGeom->GetFlukaMaterial(itmed));
334 } else {
335 SetCut(param, parval, fGeom->GetFlukaMaterial(itmed));
336 }
337
338
339
8495a208 340}
341
342// functions from GGEOM
343//_____________________________________________________________________________
344void TFluka::Gsatt(const char *name, const char *att, Int_t val)
345{
346 fGeom->Gsatt(name,att, val);
347}
348
349//______________________________________________________________________________
350Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
351 Float_t *upar, Int_t np) {
352//
353 return fGeom->Gsvolu(name, shape, nmed, upar, np);
354}
355
356//______________________________________________________________________________
357Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
358 Double_t *upar, Int_t np) {
359//
360 return fGeom->Gsvolu(name, shape, nmed, upar, np);
361}
362
363//______________________________________________________________________________
364void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv,
365 Int_t iaxis) {
366//
367 fGeom->Gsdvn(name, mother, ndiv, iaxis);
368}
369
370//______________________________________________________________________________
371void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv,
372 Int_t iaxis, Double_t c0i, Int_t numed) {
373//
374 fGeom->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed);
375}
376
377//______________________________________________________________________________
378void TFluka::Gsdvt(const char *name, const char *mother, Double_t step,
379 Int_t iaxis, Int_t numed, Int_t ndvmx) {
380//
381 fGeom->Gsdvt(name, mother, step, iaxis, numed, ndvmx);
382}
383
384//______________________________________________________________________________
385void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step,
386 Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) {
387//
388 fGeom->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx);
389}
390
391//______________________________________________________________________________
392void TFluka::Gsord(const char * /*name*/, Int_t /*iax*/) {
393//
394// Nothing to do with TGeo
395}
396
397//______________________________________________________________________________
398void TFluka::Gspos(const char *name, Int_t nr, const char *mother,
399 Double_t x, Double_t y, Double_t z, Int_t irot,
400 const char *konly) {
401//
402 fGeom->Gspos(name, nr, mother, x, y, z, irot, konly);
403}
404
405//______________________________________________________________________________
406void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
407 Double_t x, Double_t y, Double_t z, Int_t irot,
408 const char *konly, Float_t *upar, Int_t np) {
409 //
410 fGeom->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
411}
412
413//______________________________________________________________________________
414void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
415 Double_t x, Double_t y, Double_t z, Int_t irot,
416 const char *konly, Double_t *upar, Int_t np) {
417 //
418 fGeom->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
419}
420
421//______________________________________________________________________________
422void TFluka::Gsbool(const char* /*onlyVolName*/, const char* /*manyVolName*/) {
423//
424// Nothing to do with TGeo
425 Warning("Gsbool", "Not implemented with TGeo");
426}
427
428//______________________________________________________________________________
c7a186a7 429void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
430 Float_t* absco, Float_t* effic, Float_t* rindex) {
8495a208 431//
c7a186a7 432// Set Cerenkov properties for medium itmed
433//
434// npckov: number of sampling points
435// ppckov: energy values
436// absco: absorption length
437// effic: quantum efficiency
438// rindex: refraction index
439//
440//
441//
442// Create object holding Cerenkov properties
443//
444 TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex);
445//
446// Pass object to medium
447 TGeoMedium* medium = gGeoManager->GetMedium(itmed);
448 medium->SetCerenkovProperties(cerenkovProperties);
8495a208 449}
450
451//______________________________________________________________________________
452void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t * /*ppckov*/,
453 Double_t * /*absco*/, Double_t * /*effic*/, Double_t * /*rindex*/) {
454//
455// Not implemented with TGeo - what G4 did ? Any FLUKA card generated?
456 Warning("SetCerenkov", "Not implemented with TGeo");
457}
458
459// Euclid
460//______________________________________________________________________________
461void TFluka::WriteEuclid(const char* /*fileName*/, const char* /*topVol*/,
462 Int_t /*number*/, Int_t /*nlevel*/) {
463//
464// Not with TGeo
465 Warning("WriteEuclid", "Not implemented with TGeo");
466}
467
468
469
470//_____________________________________________________________________________
471// methods needed by the stepping
472//____________________________________________________________________________
473
474Int_t TFluka::GetMedium() const {
475//
476// Get the medium number for the current fluka region
477//
478 return fGeom->GetMedium(); // this I need to check due to remapping !!!
479}
480
481
482
483//____________________________________________________________________________
484// particle table usage
485// ID <--> PDG transformations
486//_____________________________________________________________________________
487Int_t TFluka::IdFromPDG(Int_t pdg) const
488{
489 //
490 // Return Fluka code from PDG and pseudo ENDF code
491
492 // Catch the feedback photons
493 if (pdg == 50000051) return (-1);
494 // MCIHAD() goes from pdg to fluka internal.
495 Int_t intfluka = mcihad(pdg);
496 // KPTOIP array goes from internal to official
497 return GetFlukaKPTOIP(intfluka);
498}
499
500//______________________________________________________________________________
501Int_t TFluka::PDGFromId(Int_t id) const
502{
503 //
504 // Return PDG code and pseudo ENDF code from Fluka code
505
506 // IPTOKP array goes from official to internal
507
508 if (id == -1) {
509// Cerenkov photon
510 if (fVerbosityLevel >= 1)
511 printf("\n PDGFromId: Cerenkov Photon \n");
512 return 50000050;
513 }
514// Error id
515 if (id == 0) {
516 if (fVerbosityLevel >= 1)
517 printf("PDGFromId: Error id = 0\n");
518 return -1;
519 }
520// Good id
521 Int_t intfluka = GetFlukaIPTOKP(id);
522 if (intfluka == 0) {
523 if (fVerbosityLevel >= 1)
524 printf("PDGFromId: Error intfluka = 0: %d\n", id);
525 return -1;
526 } else if (intfluka < 0) {
527 if (fVerbosityLevel >= 1)
528 printf("PDGFromId: Error intfluka < 0: %d\n", id);
529 return -1;
530 }
531 if (fVerbosityLevel >= 3)
532 printf("mpdgha called with %d %d \n", id, intfluka);
533 // MPDGHA() goes from fluka internal to pdg.
534 return mpdgha(intfluka);
535}
536
537//_____________________________________________________________________________
538// methods for physics management
539//____________________________________________________________________________
540//
541// set methods
542//
543
aa067d5f 544void TFluka::SetProcess(const char* flagName, Int_t flagValue, Int_t imat)
e4b9dd05 545{
546 strcpy(&fProcessFlag[fNbOfProc][0],flagName);
547 fProcessValue[fNbOfProc] = flagValue;
aa067d5f 548 fProcessMaterial[fNbOfProc] = imat;
e4b9dd05 549 fNbOfProc++;
550}
551
8495a208 552//______________________________________________________________________________
553void TFluka::SetProcess(const char* flagName, Int_t flagValue)
554{
555 Int_t i;
e4b9dd05 556 if (fNbOfProc < 100) {
557 for (i=0; i<fNbOfProc; i++) {
558 if (strcmp(&fProcessFlag[i][0],flagName) == 0) {
559 fProcessValue[fNbOfProc] = flagValue;
f37b3910 560 fProcessMaterial[fNbOfProc] = -1;
561 return;
8495a208 562 }
563 }
e4b9dd05 564 strcpy(&fProcessFlag[fNbOfProc][0],flagName);
f37b3910 565 fProcessMaterial[fNbOfProc] = -1;
566 fProcessValue[fNbOfProc++] = flagValue;
567
8495a208 568 }
569 else
570 cout << "Nb of SetProcess calls exceeds 100 - ignored" << endl;
e4b9dd05 571}
572
573//______________________________________________________________________________
574void TFluka::SetCut(const char* cutName, Double_t cutValue, Int_t imed)
575{
576 strcpy(&fCutFlag[fNbOfCut][0],cutName);
577 fCutValue[fNbOfCut] = cutValue;
aa067d5f 578 fCutMaterial[fNbOfCut] = imed;
e4b9dd05 579 fNbOfCut++;
8495a208 580}
581
582//______________________________________________________________________________
583void TFluka::SetCut(const char* cutName, Double_t cutValue)
584{
585 Int_t i;
e4b9dd05 586 if (fNbOfCut < 100) {
587 for (i=0; i<fNbOfCut; i++) {
588 if (strcmp(&fCutFlag[i][0],cutName) == 0) {
589 fCutValue[fNbOfCut] = cutValue;
590 return;
8495a208 591 }
592 }
e4b9dd05 593 strcpy(&fCutFlag[fNbOfCut][0],cutName);
aa067d5f 594 fCutMaterial[fNbOfCut] = -1;
e4b9dd05 595 fCutValue[fNbOfCut++] = cutValue;
8495a208 596 }
597 else
598 cout << "Nb of SetCut calls exceeds 100 - ignored" << endl;
8495a208 599}
600
601//______________________________________________________________________________
602Double_t TFluka::Xsec(char*, Double_t, Int_t, Int_t)
603{
604 printf("WARNING: Xsec not yet implemented !\n"); return -1.;
605}
606
607
608//______________________________________________________________________________
609void TFluka::InitPhysics()
610{
8495a208 611 Int_t i, j, k;
612 Double_t fCut;
c7a186a7 613
aa067d5f 614 FILE *pAliceCoreInp, *pAliceFlukaMat, *pAliceInp;
e749c124 615
c7a186a7 616 Double_t zero = 0.0;
617 Double_t one = 1.0;
618 Double_t two = 2.0;
619 Double_t three = 3.0;
e749c124 620
65398acc 621 Float_t fLastMaterial = fGeom->GetLastMaterialIndex();
622 printf(" last FLUKA material is %g\n", fLastMaterial);
c7a186a7 623
e4b9dd05 624 // Prepare Cerenkov
625 TList *matList = gGeoManager->GetListOfMaterials();
626 Int_t nmaterial = matList->GetSize();
627 fMaterials = new Int_t[nmaterial];
628
8495a208 629// construct file names
c7a186a7 630
0c160c74 631 TString sAliceCoreInp = getenv("ALICE_ROOT");
632 sAliceCoreInp +="/TFluka/input/";
633 TString sAliceTmp = "flukaMat.inp";
634 TString sAliceInp = GetInputFileName();
8495a208 635 sAliceCoreInp += GetCoreInputFileName();
c7a186a7 636
637// open files
638
e749c124 639 if ((pAliceCoreInp = fopen(sAliceCoreInp.Data(),"r")) == NULL) {
c7a186a7 640 printf("\nCannot open file %s\n",sAliceCoreInp.Data());
641 exit(1);
e749c124 642 }
643 if ((pAliceFlukaMat = fopen(sAliceTmp.Data(),"r")) == NULL) {
c7a186a7 644 printf("\nCannot open file %s\n",sAliceTmp.Data());
645 exit(1);
e749c124 646 }
647 if ((pAliceInp = fopen(sAliceInp.Data(),"w")) == NULL) {
c7a186a7 648 printf("\nCannot open file %s\n",sAliceInp.Data());
649 exit(1);
e749c124 650 }
8495a208 651
652// copy core input file
653 Char_t sLine[255];
654 Float_t fEventsPerRun;
c7a186a7 655
e749c124 656 while ((fgets(sLine,255,pAliceCoreInp)) != NULL) {
c7a186a7 657 if (strncmp(sLine,"GEOEND",6) != 0)
658 fprintf(pAliceInp,"%s",sLine); // copy until GEOEND card
659 else {
660 fprintf(pAliceInp,"GEOEND\n"); // add GEOEND card
661 goto flukamat;
662 }
8495a208 663 } // end of while until GEOEND card
c7a186a7 664
8495a208 665
c7a186a7 666 flukamat:
e749c124 667 while ((fgets(sLine,255,pAliceFlukaMat)) != NULL) { // copy flukaMat.inp file
c7a186a7 668 fprintf(pAliceInp,"%s\n",sLine);
8495a208 669 }
c7a186a7 670
e749c124 671 while ((fgets(sLine,255,pAliceCoreInp)) != NULL) {
c7a186a7 672 if (strncmp(sLine,"START",5) != 0)
673 fprintf(pAliceInp,"%s\n",sLine);
674 else {
675 sscanf(sLine+10,"%10f",&fEventsPerRun);
8495a208 676 goto fin;
c7a186a7 677 }
8495a208 678 } //end of while until START card
c7a186a7 679
8495a208 680fin:
681// in G3 the process control values meaning can be different for
682// different processes, but for most of them is:
683// 0 process is not activated
684// 1 process is activated WITH generation of secondaries
685// 2 process is activated WITHOUT generation of secondaries
686// if process does not generate secondaries => 1 same as 2
687//
688// Exceptions:
689// MULS: also 3
690// LOSS: also 3, 4
691// RAYL: only 0,1
692// HADR: may be > 2
693//
694
c7a186a7 695// Loop over number of SetProcess calls
e749c124 696 fprintf(pAliceInp,"*----------------------------------------------------------------------------- \n");
697 fprintf(pAliceInp,"*----- The following data are generated from SetProcess and SetCut calls ----- \n");
698 fprintf(pAliceInp,"*----------------------------------------------------------------------------- \n");
8495a208 699
e4b9dd05 700 for (i = 0; i < fNbOfProc; i++) {
701 Float_t matMin = three;
702 Float_t matMax = fLastMaterial;
703 Bool_t global = kTRUE;
aa067d5f 704 if (fProcessMaterial[i] != -1) {
705 matMin = Float_t(fProcessMaterial[i]);
e4b9dd05 706 matMax = matMin;
707 global = kFALSE;
708 }
709
8495a208 710 // annihilation
711 // G3 default value: 1
712 // G4 processes: G4eplusAnnihilation/G4IeplusAnnihilation
713 // Particles: e+
714 // Physics: EM
715 // flag = 0 no annihilation
716 // flag = 1 annihilation, decays processed
717 // flag = 2 annihilation, no decay product stored
718 // gMC ->SetProcess("ANNI",1); // EMFCUT -1. 0. 0. 3. lastmat 0. ANNH-THR
e4b9dd05 719 if (strncmp(&fProcessFlag[i][0],"ANNI",4) == 0) {
720 if (fProcessValue[i] == 1 || fProcessValue[i] == 2) {
c7a186a7 721 fprintf(pAliceInp,"*\n*Kinetic energy threshold (GeV) for e+ annihilation - resets to default=0.\n");
722 fprintf(pAliceInp,"*Generated from call: SetProcess('ANNI',1) or SetProcess('ANNI',2)\n");
723 // -one = kinetic energy threshold (GeV) for e+ annihilation (resets to default=0)
724 // zero = not used
725 // zero = not used
e4b9dd05 726 // matMin = lower bound of the material indices in which the respective thresholds apply
727 // matMax = upper bound of the material indices in which the respective thresholds apply
c7a186a7 728 // one = step length in assigning indices
729 // "ANNH-THR";
e4b9dd05 730 fprintf(pAliceInp,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",-one,zero,zero,matMin,matMax,one);
c7a186a7 731 }
e4b9dd05 732 else if (fProcessValue[i] == 0) {
c7a186a7 733 fprintf(pAliceInp,"*\n*No annihilation - no FLUKA card generated\n");
e4b9dd05 734 fprintf(pAliceInp,"*Generated from call: SetProcess('ANNI',0)\n");
c7a186a7 735 }
736 else {
737 fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('ANNI',?) call.\n");
738 fprintf(pAliceInp,"*No FLUKA card generated\n");
739 }
8495a208 740 }
8495a208 741
742 // bremsstrahlung and pair production are both activated
743 // G3 default value: 1
744 // G4 processes: G4eBremsstrahlung/G4IeBremsstrahlung,
745 // G4MuBremsstrahlung/G4IMuBremsstrahlung,
746 // G4LowEnergyBremstrahlung
747 // Particles: e-/e+; mu+/mu-
748 // Physics: EM
749 // flag = 0 no bremsstrahlung
750 // flag = 1 bremsstrahlung, photon processed
751 // flag = 2 bremsstrahlung, no photon stored
752 // gMC ->SetProcess("BREM",1); // PAIRBREM 2. 0. 0. 3. lastmat
c7a186a7 753 // EMFCUT -1. 0. 0. 3. lastmat 0. ELPO-THR
8495a208 754 // G3 default value: 1
755 // G4 processes: G4GammaConversion,
756 // G4MuPairProduction/G4IMuPairProduction
757 // G4LowEnergyGammaConversion
758 // Particles: gamma, mu
759 // Physics: EM
760 // flag = 0 no delta rays
761 // flag = 1 delta rays, secondaries processed
762 // flag = 2 delta rays, no secondaries stored
763 // gMC ->SetProcess("PAIR",1); // PAIRBREM 1. 0. 0. 3. lastmat
c7a186a7 764 // EMFCUT 0. 0. -1. 3. lastmat 0. PHOT-THR
e4b9dd05 765 else if ((strncmp(&fProcessFlag[i][0],"PAIR",4) == 0) && (fProcessValue[i] == 1 || fProcessValue[i] == 2)) {
766
767 for (j=0; j<fNbOfProc; j++) {
768 if ((strncmp(&fProcessFlag[j][0],"BREM",4) == 0) &&
769 (fProcessValue[j] == 1 || fProcessValue[j] == 2) &&
aa067d5f 770 (fProcessMaterial[j] == fProcessMaterial[i])) {
e4b9dd05 771 fprintf(pAliceInp,"*\n*Bremsstrahlung and pair production by muons and charged hadrons both activated\n");
772 fprintf(pAliceInp,"*Generated from call: SetProcess('BREM',1) and SetProcess('PAIR',1)\n");
773 fprintf(pAliceInp,"*Energy threshold set by call SetCut('BCUTM',cut) or set to 0.\n");
774 fprintf(pAliceInp,"*Energy threshold set by call SetCut('PPCUTM',cut) or set to 0.\n");
775 // three = bremsstrahlung and pair production by muons and charged hadrons both are activated
776 fprintf(pAliceInp,"PAIRBREM %10.1f",three);
777 // direct pair production by muons
778 // G4 particles: "e-", "e+"
779 // G3 default value: 0.01 GeV
780 //gMC ->SetCut("PPCUTM",cut); // total energy cut for direct pair prod. by muons
781 fCut = 0.0;
782 for (k=0; k<fNbOfCut; k++) {
783 if (strncmp(&fCutFlag[k][0],"PPCUTM",6) == 0 &&
aa067d5f 784 (fCutMaterial[k] == fProcessMaterial[i])) fCut = fCutValue[k];
e4b9dd05 785 }
786 fprintf(pAliceInp,"%10.4g",fCut);
787 // fCut; = e+, e- kinetic energy threshold (in GeV) for explicit pair production.
788 // muon and hadron bremsstrahlung
789 // G4 particles: "gamma"
790 // G3 default value: CUTGAM=0.001 GeV
791 //gMC ->SetCut("BCUTM",cut); // cut for muon and hadron bremsstrahlung
792 fCut = 0.0;
793 for (k=0; k<fNbOfCut; k++) {
794 if (strncmp(&fCutFlag[k][0],"BCUTM",5) == 0 &&
aa067d5f 795 (fCutMaterial[k] == fProcessMaterial[i])) fCut = fCutValue[k];
e4b9dd05 796 }
797 fprintf(pAliceInp,"%10.4g%10.1f%10.1f\n",fCut,matMin,matMax);
798 // fCut = photon energy threshold (GeV) for explicit bremsstrahlung production
799 // matMin = lower bound of the material indices in which the respective thresholds apply
800 // matMax = upper bound of the material indices in which the respective thresholds apply
801
802 // for e+ and e-
803 fprintf(pAliceInp,"*\n*Kinetic energy threshold (GeV) for e+/e- bremsstrahlung - resets to default=0.\n");
804 fprintf(pAliceInp,"*Generated from call: SetProcess('BREM',1);\n");
805 fCut = -1.0;
806 for (k=0; k<fNbOfCut; k++) {
807 if (strncmp(&fCutFlag[k][0],"BCUTE",5) == 0 &&
aa067d5f 808 (fCutMaterial[k] == fProcessMaterial[i])) fCut = fCutValue[k];
e4b9dd05 809 }
810 //fCut = kinetic energy threshold (GeV) for e+/e- bremsstrahlung (resets to default=0)
811 // zero = not used
812 // zero = not used
813 // matMin = lower bound of the material indices in which the respective thresholds apply
814 // matMax = upper bound of the material indices in which the respective thresholds apply
815 // one = step length in assigning indices
816 // "ELPO-THR";
817 fprintf(pAliceInp,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",fCut,zero,zero,matMin,matMax,one);
818
8495a208 819 // for e+ and e-
e4b9dd05 820 fprintf(pAliceInp,"*\n*Pair production by electrons is activated\n");
821 fprintf(pAliceInp,"*Generated from call: SetProcess('PAIR',1);\n");
822 fCut = -1.0;
823 for (k=0; k<fNbOfCut; k++) {
824 if (strncmp(&fCutFlag[k][0],"CUTGAM",6) == 0 &&
aa067d5f 825 (fCutMaterial[k] == fProcessMaterial[i])) fCut = fCutValue[k];
e4b9dd05 826 }
827 // fCut = energy threshold (GeV) for gamma pair production (< 0.0 : resets to default, = 0.0 : ignored)
828 // matMin = lower bound of the material indices in which the respective thresholds apply
829 // matMax = upper bound of the material indices in which the respective thresholds apply
830 // one = step length in assigning indices
831 fprintf(pAliceInp,"EMFCUT %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",zero,zero,fCut,matMin,matMax,one);
832 goto BOTH;
833 } // end of if for BREM
834 } // end of loop for BREM
835
836 // only pair production by muons and charged hadrons is activated
837 fprintf(pAliceInp,"*\n*Pair production by muons and charged hadrons is activated\n");
838 fprintf(pAliceInp,"*Generated from call: SetProcess('PAIR',1) or SetProcess('PAIR',2)\n");
839 fprintf(pAliceInp,"*Energy threshold set by call SetCut('PPCUTM',cut) or set to 0.\n");
840 // direct pair production by muons
841 // G4 particles: "e-", "e+"
842 // G3 default value: 0.01 GeV
843 //gMC ->SetCut("PPCUTM",cut); // total energy cut for direct pair prod. by muons
844 // one = pair production by muons and charged hadrons is activated
845 // zero = e+, e- kinetic energy threshold (in GeV) for explicit pair production.
846 // zero = no explicit bremsstrahlung production is simulated
847 // matMin = lower bound of the material indices in which the respective thresholds apply
848 // matMax = upper bound of the material indices in which the respective thresholds apply
849 fprintf(pAliceInp,"PAIRBREM %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax);
850
851 // for e+ and e-
852 fprintf(pAliceInp,"*\n*Pair production by electrons is activated\n");
853 fprintf(pAliceInp,"*Generated from call: SetProcess('PAIR',1) or SetProcess('PAIR',2)\n");
854 fCut = -1.0;
855 for (j=0; j<fNbOfCut; j++) {
856 if (strncmp(&fCutFlag[j][0],"CUTGAM",6) == 0 &&
aa067d5f 857 (fCutMaterial[j] == fProcessMaterial[i])) fCut = fCutValue[j];
e4b9dd05 858 }
859 // zero = energy threshold (GeV) for Compton scattering (= 0.0 : ignored)
860 // zero = energy threshold (GeV) for Photoelectric (= 0.0 : ignored)
861 // fCut = energy threshold (GeV) for gamma pair production (< 0.0 : resets to default, = 0.0 : ignored)
862 // matMin = lower bound of the material indices in which the respective thresholds apply
863 // matMax = upper bound of the material indices in which the respective thresholds apply
864 // one = step length in assigning indices
865 fprintf(pAliceInp,"EMFCUT %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",zero,zero,fCut,matMin,matMax,one);
8495a208 866
e4b9dd05 867 BOTH:
868 k = 0;
8495a208 869 } // end of if for PAIR
c7a186a7 870
e4b9dd05 871
872
873 // bremsstrahlung
874 // G3 default value: 1
875 // G4 processes: G4eBremsstrahlung/G4IeBremsstrahlung,
876 // G4MuBremsstrahlung/G4IMuBremsstrahlung,
877 // G4LowEnergyBremstrahlung
878 // Particles: e-/e+; mu+/mu-
879 // Physics: EM
880 // flag = 0 no bremsstrahlung
881 // flag = 1 bremsstrahlung, photon processed
882 // flag = 2 bremsstrahlung, no photon stored
883 // gMC ->SetProcess("BREM",1); // PAIRBREM 2. 0. 0. 3. lastmat
884 // EMFCUT -1. 0. 0. 3. lastmat 0. ELPO-THR
885 else if (strncmp(&fProcessFlag[i][0],"BREM",4) == 0) {
886 for (j = 0; j < fNbOfProc; j++) {
887 if ((strncmp(&fProcessFlag[j][0],"PAIR",4) == 0) &&
888 fProcessValue[j] == 1 &&
aa067d5f 889 (fProcessMaterial[j] == fProcessMaterial[i])) goto NOBREM;
e4b9dd05 890 }
891 if (fProcessValue[i] == 1 || fProcessValue[i] == 2) {
892 fprintf(pAliceInp,"*\n*Bremsstrahlung by muons and charged hadrons is activated\n");
893 fprintf(pAliceInp,"*Generated from call: SetProcess('BREM',1) or SetProcess('BREM',2)\n");
894 fprintf(pAliceInp,"*Energy threshold set by call SetCut('BCUTM',cut) or set to 0.\n");
895 // two = bremsstrahlung by muons and charged hadrons is activated
896 // zero = no meaning
897 // muon and hadron bremsstrahlung
898 // G4 particles: "gamma"
899 // G3 default value: CUTGAM=0.001 GeV
900 //gMC ->SetCut("BCUTM",cut); // cut for muon and hadron bremsstrahlung
901 fCut = 0.0;
902 for (j=0; j<fNbOfCut; j++) {
903 if (strncmp(&fCutFlag[j][0],"BCUTM",5) == 0 &&
aa067d5f 904 (fCutMaterial[j] == fProcessMaterial[i])) fCut = fCutValue[j];
c7a186a7 905 }
e4b9dd05 906 // fCut = photon energy threshold (GeV) for explicit bremsstrahlung production
907 // matMin = lower bound of the material indices in which the respective thresholds apply
908 // matMax = upper bound of the material indices in which the respective thresholds apply
909 fprintf(pAliceInp,"PAIRBREM %10.1f%10.1f%10.4g%10.1f%10.1f\n",two,zero,fCut,matMin,matMax);
c7a186a7 910
e4b9dd05 911 // for e+ and e-
912 fprintf(pAliceInp,"*\n*Kinetic energy threshold (GeV) for e+/e- bremsstrahlung - resets to default=0.\n");
913 fprintf(pAliceInp,"*Generated from call: SetProcess('BREM',1);");
914 // - one = kinetic energy threshold (GeV) for e+/e- bremsstrahlung (resets to default=0)
915 // zero = not used
916 // zero = not used
917 // matMin = lower bound of the material indices in which the respective thresholds apply
918 // matMax = upper bound of the material indices in which the respective thresholds apply
919 // one = step length in assigning indices
920 //"ELPO-THR";
921 fprintf(pAliceInp,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",-one,zero,zero,matMin,matMax,one);
922 }
923 else if (fProcessValue[i] == 0) {
924 fprintf(pAliceInp,"*\n*No bremsstrahlung - no FLUKA card generated\n");
925 fprintf(pAliceInp,"*Generated from call: SetProcess('BREM',0)\n");
926 }
927 else {
928 fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('BREM',?) call.\n");
929 fprintf(pAliceInp,"*No FLUKA card generated\n");
930 }
931 NOBREM:
932 j = 0;
933 } // end of else if (strncmp(&fProcessFlag[i][0],"BREM",4) == 0)
934
935 // Cerenkov photon generation
936 // G3 default value: 0
937 // G4 process: G4Cerenkov
938 //
939 // Particles: charged
940 // Physics: Optical
941 // flag = 0 no Cerenkov photon generation
942 // flag = 1 Cerenkov photon generation
943 // flag = 2 Cerenkov photon generation with primary stopped at each step
944 //xx gMC ->SetProcess("CKOV",1); // ??? Cerenkov photon generation
945
946 else if (strncmp(&fProcessFlag[i][0],"CKOV",4) == 0) {
947 if ((fProcessValue[i] == 1 || fProcessValue[i] == 2) && global) {
948 // Write comments
949 fprintf(pAliceInp, "* \n");
950 fprintf(pAliceInp, "*Cerenkov photon generation\n");
951 fprintf(pAliceInp, "*Generated from call: SetProcess('CKOV',1) or SetProcess('CKOV',2)\n");
952 // Loop over media
953 for (Int_t im = 0; im < nmaterial; im++)
954 {
955 TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
956 Int_t idmat = material->GetIndex();
957
aa067d5f 958 if (!global && idmat != fProcessMaterial[i]) continue;
e4b9dd05 959
960 fMaterials[idmat] = im;
961 // Skip media with no Cerenkov properties
962 TFlukaCerenkov* cerenkovProp;
963 if (!(cerenkovProp = dynamic_cast<TFlukaCerenkov*>(material->GetCerenkovProperties()))) continue;
964 //
965 // This medium has Cerenkov properties
966 //
967 //
968 // Write OPT-PROD card for each medium
969 Float_t emin = cerenkovProp->GetMinimumEnergy();
970 Float_t emax = cerenkovProp->GetMaximumEnergy();
971 fprintf(pAliceInp, "OPT-PROD %10.4g%10.4g%10.4g%10.4g%10.4g%10.4gCERENKOV\n", emin, emax, 0.,
972 Float_t(idmat), Float_t(idmat), 0.);
973 //
974 // Write OPT-PROP card for each medium
975 // Forcing FLUKA to call user routines (queffc.cxx, rflctv.cxx, rfrndx.cxx)
976 //
977 fprintf(pAliceInp, "OPT-PROP %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fWV-LIMIT\n",
978 cerenkovProp->GetMinimumWavelength(),
979 cerenkovProp->GetMaximumWavelength(),
980 cerenkovProp->GetMaximumWavelength(),
c7a186a7 981 Float_t(idmat), Float_t(idmat), 0.0);
e4b9dd05 982
983 if (cerenkovProp->IsMetal()) {
984 fprintf(pAliceInp, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fMETAL\n",
985 -100., -100., -100.,
986 Float_t(idmat), Float_t(idmat), 0.0);
987 } else {
988 fprintf(pAliceInp, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",
989 -100., -100., -100.,
990 Float_t(idmat), Float_t(idmat), 0.0);
991 }
992
993
994 for (Int_t j = 0; j < 3; j++) {
995 fprintf(pAliceInp, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f&\n",
996 -100., -100., -100.,
997 Float_t(idmat), Float_t(idmat), 0.0);
998 }
999 // Photon detection efficiency user defined
1000
1001 if (cerenkovProp->IsSensitive())
1002 fprintf(pAliceInp, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fSENSITIV\n",
1003 -100., -100., -100.,
1004 Float_t(idmat), Float_t(idmat), 0.0);
1005
1006 } // materials
1007 } else if (fProcessValue[i] == 0) {
1008 fprintf(pAliceInp,"*\n*No Cerenkov photon generation\n");
1009 fprintf(pAliceInp,"*Generated from call: SetProcess('CKOV',0)\n");
1010 // zero = not used
1011 // zero = not used
1012 // zero = not used
1013 // matMin = lower bound of the material indices in which the respective thresholds apply
1014 // matMax = upper bound of the material indices in which the respective thresholds apply
1015 // one = step length in assigning indices
1016 //"CERE-OFF";
1017 fprintf(pAliceInp,"OPT-PROD %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fCERE-OFF\n",zero,zero,zero,matMin,matMax,one);
1018 }
1019 else {
1020 fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('CKOV',?) call.\n");
1021 fprintf(pAliceInp,"*No FLUKA card generated\n");
1022 }
1023 } // end of else if (strncmp(&fProcessFlag[i][0],"CKOV",4) == 0)
1024
1025 // Compton scattering
1026 // G3 default value: 1
1027 // G4 processes: G4ComptonScattering,
1028 // G4LowEnergyCompton,
1029 // G4PolarizedComptonScattering
1030 // Particles: gamma
1031 // Physics: EM
1032 // flag = 0 no Compton scattering
1033 // flag = 1 Compton scattering, electron processed
1034 // flag = 2 Compton scattering, no electron stored
1035 // gMC ->SetProcess("COMP",1); // EMFCUT -1. 0. 0. 3. lastmat 0. PHOT-THR
1036 else if (strncmp(&fProcessFlag[i][0],"COMP",4) == 0) {
1037 if (fProcessValue[i] == 1 || fProcessValue[i] == 2) {
1038 fprintf(pAliceInp,"*\n*Energy threshold (GeV) for Compton scattering - resets to default=0.\n");
1039 fprintf(pAliceInp,"*Generated from call: SetProcess('COMP',1);\n");
1040 // - one = energy threshold (GeV) for Compton scattering - resets to default=0.
1041 // zero = not used
1042 // zero = not used
1043 // matMin = lower bound of the material indices in which the respective thresholds apply
1044 // matMax = upper bound of the material indices in which the respective thresholds apply
1045 // one = step length in assigning indices
1046 //"PHOT-THR";
1047 fprintf(pAliceInp,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",-one,zero,zero,matMin,matMax,one);
1048 }
1049 else if (fProcessValue[i] == 0) {
1050 fprintf(pAliceInp,"*\n*No Compton scattering - no FLUKA card generated\n");
1051 fprintf(pAliceInp,"*Generated from call: SetProcess('COMP',0)\n");
1052 }
1053 else {
1054 fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('COMP',?) call.\n");
1055 fprintf(pAliceInp,"*No FLUKA card generated\n");
1056 }
1057 } // end of else if (strncmp(&fProcessFlag[i][0],"COMP",4) == 0)
1058
1059 // decay
1060 // G3 default value: 1
1061 // G4 process: G4Decay
1062 //
1063 // Particles: all which decay is applicable for
1064 // Physics: General
1065 // flag = 0 no decays
1066 // flag = 1 decays, secondaries processed
1067 // flag = 2 decays, no secondaries stored
1068 //gMC ->SetProcess("DCAY",1); // not available
1069 else if ((strncmp(&fProcessFlag[i][0],"DCAY",4) == 0) && fProcessValue[i] == 1)
1070 cout << "SetProcess for flag=" << &fProcessFlag[i][0] << " value=" << fProcessValue[i] << " not avaliable!" << endl;
1071
1072 // delta-ray
1073 // G3 default value: 2
1074 // !! G4 treats delta rays in different way
1075 // G4 processes: G4eIonisation/G4IeIonization,
1076 // G4MuIonisation/G4IMuIonization,
1077 // G4hIonisation/G4IhIonisation
1078 // Particles: charged
1079 // Physics: EM
1080 // flag = 0 no energy loss
1081 // flag = 1 restricted energy loss fluctuations
1082 // flag = 2 complete energy loss fluctuations
1083 // flag = 3 same as 1
1084 // flag = 4 no energy loss fluctuations
1085 // gMC ->SetProcess("DRAY",0); // DELTARAY 1.E+6 0. 0. 3. lastmat 0.
1086 else if (strncmp(&fProcessFlag[i][0],"DRAY",4) == 0) {
1087 if (fProcessValue[i] == 0 || fProcessValue[i] == 4) {
1088 fprintf(pAliceInp,"*\n*Kinetic energy threshold (GeV) for delta ray production\n");
1089 fprintf(pAliceInp,"*Generated from call: SetProcess('DRAY',0) or SetProcess('DRAY',4)\n");
1090 fprintf(pAliceInp,"*No delta ray production by muons - threshold set artificially high\n");
1091 Double_t emin = 1.0e+6; // kinetic energy threshold (GeV) for delta ray production (discrete energy transfer)
1092 // zero = ignored
1093 // zero = ignored
1094 // matMin = lower bound of the material indices in which the respective thresholds apply
1095 // matMax = upper bound of the material indices in which the respective thresholds apply
1096 // one = step length in assigning indices
1097 fprintf(pAliceInp,"DELTARAY %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n",emin,zero,zero,matMin,matMax,one);
1098 }
1099 else if (fProcessValue[i] == 1 || fProcessValue[i] == 2 || fProcessValue[i] == 3) {
1100 fprintf(pAliceInp,"*\n*Kinetic energy threshold (GeV) for delta ray production\n");
1101 fprintf(pAliceInp,"*Generated from call: SetProcess('DRAY',flag), flag=1,2,3\n");
1102 fprintf(pAliceInp,"*Delta ray production by muons switched on\n");
1103 fprintf(pAliceInp,"*Energy threshold set by call SetCut('DCUTM',cut) or set to 1.0e+6.\n");
1104 fCut = 1.0e+6;
1105 for (j = 0; j < fNbOfCut; j++) {
1106 if (strncmp(&fCutFlag[j][0],"DCUTM",5) == 0 &&
aa067d5f 1107 fCutMaterial[j] == fProcessMaterial[i]) fCut = fCutValue[j];
c7a186a7 1108 }
e4b9dd05 1109 // fCut = kinetic energy threshold (GeV) for delta ray production (discrete energy transfer)
1110 // zero = ignored
1111 // zero = ignored
1112 // matMin = lower bound of the material indices in which the respective thresholds apply
1113 // matMax = upper bound of the material indices in which the respective thresholds apply
1114 // one = step length in assigning indices
1115 fprintf(pAliceInp,"DELTARAY %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n",fCut,zero,zero,matMin,matMax,one);
1116 }
1117 else {
1118 fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('DRAY',?) call.\n");
1119 fprintf(pAliceInp,"*No FLUKA card generated\n");
1120 }
1121 } // end of else if (strncmp(&fProcessFlag[i][0],"DRAY",4) == 0)
1122
1123 // hadronic process
1124 // G3 default value: 1
1125 // G4 processes: all defined by TG4PhysicsConstructorHadron
1126 //
1127 // Particles: hadrons
1128 // Physics: Hadron
1129 // flag = 0 no multiple scattering
1130 // flag = 1 hadronic interactions, secondaries processed
1131 // flag = 2 hadronic interactions, no secondaries stored
1132 // gMC ->SetProcess("HADR",1); // ??? hadronic process
1133 //Select pure GEANH (HADR 1) or GEANH/NUCRIN (HADR 3) ?????
1134 else if (strncmp(&fProcessFlag[i][0],"HADR",4) == 0) {
1135 if (fProcessValue[i] == 1 || fProcessValue[i] == 2) {
1136 fprintf(pAliceInp,"*\n*Hadronic interaction is ON by default in FLUKA\n");
1137 fprintf(pAliceInp,"*No FLUKA card generated\n");
1138 }
1139 else if (fProcessValue[i] == 0) {
1140 fprintf(pAliceInp,"*\n*Hadronic interaction is set OFF\n");
1141 fprintf(pAliceInp,"*Generated from call: SetProcess('HADR',0);\n");
1142 // zero = ignored
1143 // three = multiple scattering for hadrons and muons is completely suppressed
1144 // zero = no spin-relativistic corrections
1145 // matMin = lower bound of the material indices in which the respective thresholds apply
1146 // matMax = upper bound of the material indices in which the respective thresholds apply
1147 fprintf(pAliceInp,"MULSOPT %10.1f%10.1f%10.1f%10.1f%10.1f\n",zero,three,zero,matMin,matMax);
1148
1149 }
1150 else {
1151 fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('HADR',?) call.\n");
1152 fprintf(pAliceInp,"*No FLUKA card generated\n");
1153 }
1154 } // end of else if (strncmp(&fProcessFlag[i][0],"HADR",4) == 0)
1155
1156
1157 // energy loss
1158 // G3 default value: 2
1159 // G4 processes: G4eIonisation/G4IeIonization,
1160 // G4MuIonisation/G4IMuIonization,
1161 // G4hIonisation/G4IhIonisation
1162 //
1163 // Particles: charged
1164 // Physics: EM
1165 // flag=0 no energy loss
1166 // flag=1 restricted energy loss fluctuations
1167 // flag=2 complete energy loss fluctuations
1168 // flag=3 same as 1
1169 // flag=4 no energy loss fluctuations
1170 // If the value ILOSS is changed, then (in G3) cross-sections and energy
1171 // loss tables must be recomputed via the command 'PHYSI'
1172 // gMC ->SetProcess("LOSS",2); // ??? IONFLUCT ? energy loss
1173 else if (strncmp(&fProcessFlag[i][0],"LOSS",4) == 0) {
1174 if (fProcessValue[i] == 2) { // complete energy loss fluctuations
1175 fprintf(pAliceInp,"*\n*Complete energy loss fluctuations do not exist in FLUKA\n");
1176 fprintf(pAliceInp,"*Generated from call: SetProcess('LOSS',2);\n");
1177 fprintf(pAliceInp,"*flag=2=complete energy loss fluctuations\n");
1178 fprintf(pAliceInp,"*No FLUKA card generated\n");
1179 }
1180 else if (fProcessValue[i] == 1 || fProcessValue[i] == 3) { // restricted energy loss fluctuations
1181 fprintf(pAliceInp,"*\n*Restricted energy loss fluctuations\n");
1182 fprintf(pAliceInp,"*Generated from call: SetProcess('LOSS',1) or SetProcess('LOSS',3)\n");
1183 // one = restricted energy loss fluctuations (for hadrons and muons) switched on
1184 // one = restricted energy loss fluctuations (for e+ and e-) switched on
1185 // one = minimal accuracy
1186 // matMin = lower bound of the material indices in which the respective thresholds apply
1187 // upper bound of the material indices in which the respective thresholds apply
1188 fprintf(pAliceInp,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,one,one,matMin,matMax);
1189 }
1190 else if (fProcessValue[i] == 4) { // no energy loss fluctuations
1191 fprintf(pAliceInp,"*\n*No energy loss fluctuations\n");
1192 fprintf(pAliceInp,"*\n*Generated from call: SetProcess('LOSS',4)\n");
1193 // - one = restricted energy loss fluctuations (for hadrons and muons) switched off
1194 // - one = restricted energy loss fluctuations (for e+ and e-) switched off
1195 // one = minimal accuracy
1196 // matMin = lower bound of the material indices in which the respective thresholds apply
1197 // matMax = upper bound of the material indices in which the respective thresholds apply
1198 fprintf(pAliceInp,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n",-one,-one,one,matMin,matMax);
1199 }
1200 else {
1201 fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('LOSS',?) call.\n");
1202 fprintf(pAliceInp,"*No FLUKA card generated\n");
1203 }
1204 } // end of else if (strncmp(&fProcessFlag[i][0],"LOSS",4) == 0)
1205
1206
1207 // multiple scattering
1208 // G3 default value: 1
1209 // G4 process: G4MultipleScattering/G4IMultipleScattering
1210 //
1211 // Particles: charged
1212 // Physics: EM
1213 // flag = 0 no multiple scattering
1214 // flag = 1 Moliere or Coulomb scattering
1215 // flag = 2 Moliere or Coulomb scattering
1216 // flag = 3 Gaussian scattering
1217 // gMC ->SetProcess("MULS",1); // MULSOPT multiple scattering
1218 else if (strncmp(&fProcessFlag[i][0],"MULS",4) == 0) {
1219 if (fProcessValue[i] == 1 || fProcessValue[i] == 2 || fProcessValue[i] == 3) {
1220 fprintf(pAliceInp,"*\n*Multiple scattering is ON by default for e+e- and for hadrons/muons\n");
1221 fprintf(pAliceInp,"*No FLUKA card generated\n");
1222 }
1223 else if (fProcessValue[i] == 0) {
1224 fprintf(pAliceInp,"*\n*Multiple scattering is set OFF\n");
1225 fprintf(pAliceInp,"*Generated from call: SetProcess('MULS',0);\n");
1226 // zero = ignored
1227 // three = multiple scattering for hadrons and muons is completely suppressed
1228 // three = multiple scattering for e+ and e- is completely suppressed
1229 // matMin = lower bound of the material indices in which the respective thresholds apply
1230 // matMax = upper bound of the material indices in which the respective thresholds apply
1231 fprintf(pAliceInp,"MULSOPT %10.1f%10.1f%10.1f%10.1f%10.1f\n",zero,three,three,matMin,matMax);
1232 }
1233 else {
1234 fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('MULS',?) call.\n");
1235 fprintf(pAliceInp,"*No FLUKA card generated\n");
1236 }
1237 } // end of else if (strncmp(&fProcessFlag[i][0],"MULS",4) == 0)
1238
2b8052b0 1239
e4b9dd05 1240 // muon nuclear interaction
1241 // G3 default value: 0
1242 // G4 processes: G4MuNuclearInteraction,
1243 // G4MuonMinusCaptureAtRest
1244 //
1245 // Particles: mu
1246 // Physics: Not set
1247 // flag = 0 no muon-nuclear interaction
1248 // flag = 1 nuclear interaction, secondaries processed
1249 // flag = 2 nuclear interaction, secondaries not processed
1250 // gMC ->SetProcess("MUNU",1); // MUPHOTON 1. 0. 0. 3. lastmat
1251 else if (strncmp(&fProcessFlag[i][0],"MUNU",4) == 0) {
1252 if (fProcessValue[i] == 1) {
1253 fprintf(pAliceInp,"*\n*Muon nuclear interactions with production of secondary hadrons\n");
1254 fprintf(pAliceInp,"*\n*Generated from call: SetProcess('MUNU',1);\n");
1255 // one = full simulation of muon nuclear interactions and production of secondary hadrons
1256 // zero = ratio of longitudinal to transverse virtual photon cross-section - Default = 0.25.
1257 // zero = fraction of rho-like interactions ( must be < 1) - Default = 0.75.
1258 // matMin = lower bound of the material indices in which the respective thresholds apply
1259 // matMax = upper bound of the material indices in which the respective thresholds apply
1260 fprintf(pAliceInp,"MUPHOTON %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax);
1261 }
1262 else if (fProcessValue[i] == 2) {
1263 fprintf(pAliceInp,"*\n*Muon nuclear interactions without production of secondary hadrons\n");
1264 fprintf(pAliceInp,"*Generated from call: SetProcess('MUNU',2);\n");
1265 // two = full simulation of muon nuclear interactions and production of secondary hadrons
1266 // zero = ratio of longitudinal to transverse virtual photon cross-section - Default = 0.25.
1267 // zero = fraction of rho-like interactions ( must be < 1) - Default = 0.75.
1268 // matMin = lower bound of the material indices in which the respective thresholds apply
1269 // matMax = upper bound of the material indices in which the respective thresholds apply
1270 fprintf(pAliceInp,"MUPHOTON %10.1f%10.1f%10.1f%10.1f%10.1f\n",two,zero,zero,matMin,matMax);
1271 }
1272 else if (fProcessValue[i] == 0) {
1273 fprintf(pAliceInp,"*\n*No muon nuclear interaction - no FLUKA card generated\n");
1274 fprintf(pAliceInp,"*Generated from call: SetProcess('MUNU',0)\n");
1275 }
1276 else {
1277 fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('MUNU',?) call.\n");
1278 fprintf(pAliceInp,"*No FLUKA card generated\n");
1279 }
1280 } // end of else if (strncmp(&fProcessFlag[i][0],"MUNU",4) == 0)
c7a186a7 1281
e4b9dd05 1282
1283 // photofission
1284 // G3 default value: 0
1285 // G4 process: ??
1286 //
1287 // Particles: gamma
1288 // Physics: ??
1289 // gMC ->SetProcess("PFIS",0); // PHOTONUC -1. 0. 0. 3. lastmat 0.
1290 // flag = 0 no photon fission
1291 // flag = 1 photon fission, secondaries processed
1292 // flag = 2 photon fission, no secondaries stored
1293 else if (strncmp(&fProcessFlag[i][0],"PFIS",4) == 0) {
1294 if (fProcessValue[i] == 0) {
1295 fprintf(pAliceInp,"*\n*No photonuclear interactions\n");
1296 fprintf(pAliceInp,"*Generated from call: SetProcess('PFIS',0);\n");
1297 // - one = no photonuclear interactions
1298 // zero = not used
1299 // zero = not used
1300 // matMin = lower bound of the material indices in which the respective thresholds apply
1301 // matMax = upper bound of the material indices in which the respective thresholds apply
1302 fprintf(pAliceInp,"PHOTONUC %10.1f%10.1f%10.1f%10.1f%10.1f\n",-one,zero,zero,matMin,matMax);
1303 }
1304 else if (fProcessValue[i] == 1) {
1305 fprintf(pAliceInp,"*\n*Photon nuclear interactions are activated at all energies\n");
1306 fprintf(pAliceInp,"*Generated from call: SetProcess('PFIS',1);\n");
1307 // one = photonuclear interactions are activated at all energies
1308 // zero = not used
1309 // zero = not used
1310 // matMin = lower bound of the material indices in which the respective thresholds apply
1311 // matMax = upper bound of the material indices in which the respective thresholds apply
1312 fprintf(pAliceInp,"PHOTONUC %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax);
1313 }
1314 else if (fProcessValue[i] == 0) {
1315 fprintf(pAliceInp,"*\n*No photofission - no FLUKA card generated\n");
1316 fprintf(pAliceInp,"*Generated from call: SetProcess('PFIS',0)\n");
1317 }
1318 else {
1319 fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('PFIS',?) call.\n");
1320 fprintf(pAliceInp,"*No FLUKA card generated\n");
1321 }
8495a208 1322 }
8495a208 1323
e4b9dd05 1324
1325 // photo electric effect
1326 // G3 default value: 1
1327 // G4 processes: G4PhotoElectricEffect
1328 // G4LowEnergyPhotoElectric
1329 // Particles: gamma
1330 // Physics: EM
1331 // flag = 0 no photo electric effect
1332 // flag = 1 photo electric effect, electron processed
1333 // flag = 2 photo electric effect, no electron stored
1334 // gMC ->SetProcess("PHOT",1); // EMFCUT 0. -1. 0. 3. lastmat 0. PHOT-THR
1335 else if (strncmp(&fProcessFlag[i][0],"PHOT",4) == 0) {
1336 if (fProcessValue[i] == 1 || fProcessValue[i] == 2) {
1337 fprintf(pAliceInp,"*\n*Photo electric effect is activated\n");
1338 fprintf(pAliceInp,"*Generated from call: SetProcess('PHOT',1);\n");
1339 // zero = ignored
1340 // - one = resets to default=0.
1341 // zero = ignored
1342 // matMin = lower bound of the material indices in which the respective thresholds apply
1343 // matMax = upper bound of the material indices in which the respective thresholds apply
1344 // one = step length in assigning indices
1345 //"PHOT-THR";
1346 fprintf(pAliceInp,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",zero,-one,zero,matMin,matMax,one);
1347 }
1348 else if (fProcessValue[i] == 0) {
1349 fprintf(pAliceInp,"*\n*No photo electric effect - no FLUKA card generated\n");
1350 fprintf(pAliceInp,"*Generated from call: SetProcess('PHOT',0)\n");
1351 }
1352 else {
1353 fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('PHOT',?) call.\n");
1354 fprintf(pAliceInp,"*No FLUKA card generated\n");
1355 }
1356 } // else if (strncmp(&fProcessFlag[i][0],"PHOT",4) == 0)
8495a208 1357
e4b9dd05 1358
1359 // Rayleigh scattering
1360 // G3 default value: 0
1361 // G4 process: G4OpRayleigh
1362 //
1363 // Particles: optical photon
1364 // Physics: Optical
1365 // flag = 0 Rayleigh scattering off
1366 // flag = 1 Rayleigh scattering on
1367 //xx gMC ->SetProcess("RAYL",1);
1368 else if (strncmp(&fProcessFlag[i][0],"RAYL",4) == 0) {
1369 if (fProcessValue[i] == 1) {
1370 fprintf(pAliceInp,"*\n*Rayleigh scattering is ON by default in FLUKA\n");
1371 fprintf(pAliceInp,"*No FLUKA card generated\n");
1372 }
1373 else if (fProcessValue[i] == 0) {
1374 fprintf(pAliceInp,"*\n*Rayleigh scattering is set OFF\n");
1375 fprintf(pAliceInp,"*Generated from call: SetProcess('RAYL',0);\n");
1376 // - one = no Rayleigh scattering and no binding corrections for Compton
1377 // matMin = lower bound of the material indices in which the respective thresholds apply
1378 // matMax = upper bound of the material indices in which the respective thresholds apply
1379 fprintf(pAliceInp,"EMFRAY %10.1f%10.1f%10.1f%10.1f\n",-one,three,matMin,matMax);
1380 }
1381 else {
1382 fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('RAYL',?) call.\n");
1383 fprintf(pAliceInp,"*No FLUKA card generated\n");
1384 }
1385 } // end of else if (strncmp(&fProcessFlag[i][0],"RAYL",4) == 0)
1386
1387
1388 // synchrotron radiation in magnetic field
1389 // G3 default value: 0
1390 // G4 process: G4SynchrotronRadiation
1391 //
1392 // Particles: ??
1393 // Physics: Not set
1394 // flag = 0 no synchrotron radiation
1395 // flag = 1 synchrotron radiation
1396 //xx gMC ->SetProcess("SYNC",1); // synchrotron radiation generation
1397 else if (strncmp(&fProcessFlag[i][0],"SYNC",4) == 0) {
1398 fprintf(pAliceInp,"*\n*Synchrotron radiation generation is NOT implemented in FLUKA\n");
1399 fprintf(pAliceInp,"*No FLUKA card generated\n");
8495a208 1400 }
e4b9dd05 1401
1402
1403 // Automatic calculation of tracking medium parameters
1404 // flag = 0 no automatic calculation
1405 // flag = 1 automatic calculation
1406 //xx gMC ->SetProcess("AUTO",1); // ??? automatic computation of the tracking medium parameters
1407 else if (strncmp(&fProcessFlag[i][0],"AUTO",4) == 0) {
1408 fprintf(pAliceInp,"*\n*Automatic calculation of tracking medium parameters is always ON in FLUKA\n");
1409 fprintf(pAliceInp,"*No FLUKA card generated\n");
8495a208 1410 }
e4b9dd05 1411
1412
1413 // To control energy loss fluctuation model
1414 // flag = 0 Urban model
1415 // flag = 1 PAI model
1416 // flag = 2 PAI+ASHO model (not active at the moment)
1417 //xx gMC ->SetProcess("STRA",1); // ??? energy fluctuation model
1418 else if (strncmp(&fProcessFlag[i][0],"STRA",4) == 0) {
1419 if (fProcessValue[i] == 0 || fProcessValue[i] == 2 || fProcessValue[i] == 3) {
1420 fprintf(pAliceInp,"*\n*Ionization energy losses calculation is activated\n");
1421 fprintf(pAliceInp,"*Generated from call: SetProcess('STRA',n);, n=0,1,2\n");
1422 // one = restricted energy loss fluctuations (for hadrons and muons) switched on
1423 // one = restricted energy loss fluctuations (for e+ and e-) switched on
1424 // one = minimal accuracy
1425 // matMin = lower bound of the material indices in which the respective thresholds apply
1426 // matMax = upper bound of the material indices in which the respective thresholds apply
1427 fprintf(pAliceInp,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,one,one,matMin,matMax);
1428 }
1429 else {
1430 fprintf(pAliceInp,"*\n*Illegal flag value in SetProcess('STRA',?) call.\n");
1431 fprintf(pAliceInp,"*No FLUKA card generated\n");
1432 }
1433 } // else if (strncmp(&fProcessFlag[i][0],"STRA",4) == 0)
1434
8495a208 1435
8495a208 1436
1437
e4b9dd05 1438 else { // processes not yet treated
1439
1440 // light photon absorption (Cerenkov photons)
1441 // it is turned on when Cerenkov process is turned on
1442 // G3 default value: 0
1443 // G4 process: G4OpAbsorption, G4OpBoundaryProcess
1444 //
1445 // Particles: optical photon
1446 // Physics: Optical
1447 // flag = 0 no absorption of Cerenkov photons
1448 // flag = 1 absorption of Cerenkov photons
1449 // gMC ->SetProcess("LABS",2); // ??? Cerenkov light absorption
1450
1451
1452
1453 cout << "SetProcess for flag=" << &fProcessFlag[i][0] << " value=" << fProcessValue[i] << " not yet implemented!" << endl;
8495a208 1454 }
e4b9dd05 1455 } //end of loop number of SetProcess calls
1456
8495a208 1457
e4b9dd05 1458// Loop over number of SetCut calls
1459 for (Int_t i = 0; i < fNbOfCut; i++) {
1460 Float_t matMin = three;
1461 Float_t matMax = fLastMaterial;
1462 Bool_t global = kTRUE;
aa067d5f 1463 if (fCutMaterial[i] != -1) {
1464 matMin = Float_t(fCutMaterial[i]);
1465 matMax = matMin;
1466 global = kFALSE;
8495a208 1467 }
aa067d5f 1468
1469 // cuts handled in SetProcess calls
e4b9dd05 1470 if (strncmp(&fCutFlag[i][0],"BCUTM",5) == 0) continue;
1471 else if (strncmp(&fCutFlag[i][0],"BCUTE",5) == 0) continue;
1472 else if (strncmp(&fCutFlag[i][0],"DCUTM",5) == 0) continue;
1473 else if (strncmp(&fCutFlag[i][0],"PPCUTM",6) == 0) continue;
1474
aa067d5f 1475 // delta-rays by electrons
1476 // G4 particles: "e-"
1477 // G3 default value: 10**4 GeV
1478 // gMC ->SetCut("DCUTE",cut); // cut for deltarays by electrons
1479 else if (strncmp(&fCutFlag[i][0],"DCUTE",5) == 0) {
1480 fprintf(pAliceInp,"*\n*Cut for delta rays by electrons\n");
1481 fprintf(pAliceInp,"*Generated from call: SetCut('DCUTE',cut);\n");
1482 // -fCutValue[i];
1483 // zero = ignored
1484 // zero = ignored
1485 // matMin = lower bound of the material indices in which the respective thresholds apply
1486 // matMax = upper bound of the material indices in which the respective thresholds apply
1487 // loop over materials for EMFCUT FLUKA cards
1488 for (j=0; j < matMax-matMin+1; j++) {
1489 Int_t nreg, imat, *reglist;
1490 Float_t ireg;
1491 imat = (Int_t) matMin + j;
1492 reglist = fGeom->GetMaterialList(imat, nreg);
1493 // loop over regions of a given material
1494 for (k=0; k<nreg; k++) {
1495 ireg = reglist[k];
1496 fprintf(pAliceInp,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f\n",-fCutValue[i],zero,zero,ireg,ireg);
1497 }
1498 }
1499 fprintf(pAliceInp,"DELTARAY %10.4g%10.3f%10.3f%10.1f%10.1f%10.1f\n",fCutValue[i], 100., 1.03, matMin, matMax, 1.0);
1500 fprintf(pAliceInp,"STEPSIZE %10.4g%10.3f%10.3f%10.1f%10.1f\n", 0.1, 1.0, 1.00,
1501 Float_t(gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1), 1.0);
1502 } // end of if for delta-rays by electrons
1503
1504
e4b9dd05 1505 // gammas
1506 // G4 particles: "gamma"
1507 // G3 default value: 0.001 GeV
1508 // gMC ->SetCut("CUTGAM",cut); // cut for gammas
1509
1510 else if (strncmp(&fCutFlag[i][0],"CUTGAM",6) == 0 && global) {
aa067d5f 1511 fprintf(pAliceInp,"*\n*Cut for gamma\n");
1512 fprintf(pAliceInp,"*Generated from call: SetCut('CUTGAM',cut);\n");
1513 // -fCutValue[i];
1514 // 7.0 = lower bound of the particle id-numbers to which the cut-off
1515 fprintf(pAliceInp,"PART-THR %10.4g%10.1f\n",-fCutValue[i],7.0);
8495a208 1516 }
e4b9dd05 1517 else if (strncmp(&fCutFlag[i][0],"CUTGAM",6) == 0 && !global) {
aa067d5f 1518 fprintf(pAliceInp,"*\n*Cut specific to material for gamma\n");
1519 fprintf(pAliceInp,"*Generated from call: SetCut('CUTGAM',cut);\n");
1520 // fCutValue[i];
1521 // loop over materials for EMFCUT FLUKA cards
1522 for (j=0; j < matMax-matMin+1; j++) {
1523 Int_t nreg, imat, *reglist;
1524 Float_t ireg;
1525 imat = (Int_t) matMin + j;
1526 reglist = fGeom->GetMaterialList(imat, nreg);
1527 // loop over regions of a given material
1528 for (Int_t k=0; k<nreg; k++) {
1529 ireg = reglist[k];
1530 fprintf(pAliceInp,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", zero, fCutValue[i], zero, ireg, ireg, one);
1531 }
1532 }
1533 } // end of else if for gamma
1534
1535
e4b9dd05 1536 // electrons
1537 // G4 particles: "e-"
1538 // ?? positrons
1539 // G3 default value: 0.001 GeV
1540 //gMC ->SetCut("CUTELE",cut); // cut for e+,e-
1541 else if (strncmp(&fCutFlag[i][0],"CUTELE",6) == 0 && global) {
aa067d5f 1542 fprintf(pAliceInp,"*\n*Cut for electrons\n");
1543 fprintf(pAliceInp,"*Generated from call: SetCut('CUTELE',cut);\n");
1544 // -fCutValue[i];
1545 // three = lower bound of the particle id-numbers to which the cut-off
1546 // 4.0 = upper bound of the particle id-numbers to which the cut-off
1547 // one = step length in assigning numbers
1548 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f%10.1f\n",-fCutValue[i],three,4.0,one);
8495a208 1549 }
e4b9dd05 1550 else if (strncmp(&fCutFlag[i][0],"CUTELE",6) == 0 && !global) {
aa067d5f 1551 fprintf(pAliceInp,"*\n*Cut specific to material for electrons\n");
1552 fprintf(pAliceInp,"*Generated from call: SetCut('CUTELE',cut);\n");
1553 // -fCutValue[i];
1554 // loop over materials for EMFCUT FLUKA cards
1555 for (j=0; j < matMax-matMin+1; j++) {
1556 Int_t nreg, imat, *reglist;
1557 Float_t ireg;
1558 imat = (Int_t) matMin + j;
1559 reglist = fGeom->GetMaterialList(imat, nreg);
1560 // loop over regions of a given material
1561 for (k=0; k<nreg; k++) {
1562 ireg = reglist[k];
1563 fprintf(pAliceInp,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", -fCutValue[i], zero, zero, ireg, ireg, one);
1564 }
1565 }
1566 } // end of else if for electrons
1567
1568
e4b9dd05 1569 // neutral hadrons
1570 // G4 particles: of type "baryon", "meson", "nucleus" with zero charge
1571 // G3 default value: 0.01 GeV
1572 //gMC ->SetCut("CUTNEU",cut); // cut for neutral hadrons
1573 else if (strncmp(&fCutFlag[i][0],"CUTNEU",6) == 0 && global) {
aa067d5f 1574 fprintf(pAliceInp,"*\n*Cut for neutral hadrons\n");
1575 fprintf(pAliceInp,"*Generated from call: SetCut('CUTNEU',cut);\n");
e4b9dd05 1576
aa067d5f 1577 // 8.0 = Neutron
1578 // 9.0 = Antineutron
1579 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],8.0,9.0);
e4b9dd05 1580
aa067d5f 1581 // 12.0 = Kaon zero long
1582 // 12.0 = Kaon zero long
1583 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],12.0,12.0);
e4b9dd05 1584
aa067d5f 1585 // 17.0 = Lambda, 18.0 = Antilambda
1586 // 19.0 = Kaon zero short
1587 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],17.0,19.0);
e4b9dd05 1588
aa067d5f 1589 // 22.0 = Sigma zero, Pion zero, Kaon zero
1590 // 25.0 = Antikaon zero
1591 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],22.0,25.0);
e4b9dd05 1592
aa067d5f 1593 // 32.0 = Antisigma zero
1594 // 32.0 = Antisigma zero
1595 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],32.0,32.0);
e4b9dd05 1596
aa067d5f 1597 // 34.0 = Xi zero
1598 // 35.0 = AntiXi zero
1599 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],34.0,35.0);
e4b9dd05 1600
aa067d5f 1601 // 47.0 = D zero
1602 // 48.0 = AntiD zero
1603 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],47.0,48.0);
e4b9dd05 1604
aa067d5f 1605 // 53.0 = Xi_c zero
1606 // 53.0 = Xi_c zero
1607 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],53.0,53.0);
e4b9dd05 1608
aa067d5f 1609 // 55.0 = Xi'_c zero
1610 // 56.0 = Omega_c zero
1611 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],55.0,56.0);
e4b9dd05 1612
aa067d5f 1613 // 59.0 = AntiXi_c zero
1614 // 59.0 = AntiXi_c zero
1615 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],59.0,59.0);
e4b9dd05 1616
aa067d5f 1617 // 61.0 = AntiXi'_c zero
1618 // 62.0 = AntiOmega_c zero
1619 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],61.0,62.0);
8495a208 1620 }
e4b9dd05 1621
1622 // charged hadrons
1623 // G4 particles: of type "baryon", "meson", "nucleus" with non-zero charge
1624 // G3 default value: 0.01 GeV
1625 //gMC ->SetCut("CUTHAD",cut); // cut for charged hadrons
1626 else if (strncmp(&fCutFlag[i][0],"CUTHAD",6) == 0 && global) {
aa067d5f 1627 fprintf(pAliceInp,"*\n*Cut for charged hadrons\n");
1628 fprintf(pAliceInp,"*Generated from call: SetCut('CUTHAD',cut);\n");
e4b9dd05 1629
aa067d5f 1630 // 1.0 = Proton
1631 // 2.0 = Antiproton
1632 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],1.0,2.0);
e4b9dd05 1633
aa067d5f 1634 // 13.0 = Positive Pion, Negative Pion, Positive Kaon
1635 // 16.0 = Negative Kaon
1636 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],13.0,16.0);
e4b9dd05 1637
aa067d5f 1638 // 20.0 = Negative Sigma
1639 // 21.0 = Positive Sigma
1640 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],20.0,21.0);
e4b9dd05 1641
aa067d5f 1642 // 31.0 = Antisigma minus
1643 // 33.0 = Antisigma plus
1644 // 2.0 = step length
1645 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f%10.1f\n",-fCutValue[i],31.0,33.0,2.0);
e4b9dd05 1646
aa067d5f 1647 // 36.0 = Negative Xi, Positive Xi, Omega minus
1648 // 39.0 = Antiomega
1649 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],36.0,39.0);
e4b9dd05 1650
aa067d5f 1651 // 45.0 = D plus
1652 // 46.0 = D minus
1653 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],45.0,46.0);
e4b9dd05 1654
aa067d5f 1655 // 49.0 = D_s plus, D_s minus, Lambda_c plus
1656 // 52.0 = Xi_c plus
1657 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],49.0,52.0);
e4b9dd05 1658
aa067d5f 1659 // 54.0 = Xi'_c plus
1660 // 60.0 = AntiXi'_c minus
1661 // 6.0 = step length
1662 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f%10.1f\n",-fCutValue[i],54.0,60.0,6.0);
e4b9dd05 1663
aa067d5f 1664 // 57.0 = Antilambda_c minus
1665 // 58.0 = AntiXi_c minus
1666 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],57.0,58.0);
8495a208 1667 }
e4b9dd05 1668
1669 // muons
1670 // G4 particles: "mu+", "mu-"
1671 // G3 default value: 0.01 GeV
1672 //gMC ->SetCut("CUTMUO",cut); // cut for mu+, mu-
1673 else if (strncmp(&fCutFlag[i][0],"CUTMUO",6)== 0 && global) {
aa067d5f 1674 fprintf(pAliceInp,"*\n*Cut for muons\n");
1675 fprintf(pAliceInp,"*Generated from call: SetCut('CUTMUO',cut);\n");
1676 // 10.0 = Muon+
1677 // 11.0 = Muon-
1678 fprintf(pAliceInp,"PART-THR %10.4g%10.1f%10.1f\n",-fCutValue[i],10.0,11.0);
8495a208 1679 }
e4b9dd05 1680
e4b9dd05 1681 //
1682 // time of flight cut in seconds
1683 // G4 particles: all
1684 // G3 default value: 0.01 GeV
1685 //gMC ->SetCut("TOFMAX",tofmax); // time of flight cuts in seconds
1686 else if (strncmp(&fCutFlag[i][0],"TOFMAX",6) == 0) {
aa067d5f 1687 fprintf(pAliceInp,"*\n*Time of flight cuts in seconds\n");
1688 fprintf(pAliceInp,"*Generated from call: SetCut('TOFMAX',tofmax);\n");
1689 // zero = ignored
1690 // zero = ignored
1691 // -6.0 = lower bound of the particle numbers for which the transport time cut-off and/or the start signal is to be applied
1692 // 64.0 = upper bound of the particle numbers for which the transport time cut-off and/or the start signal is to be applied
1693 fprintf(pAliceInp,"TIME-CUT %10.4g%10.1f%10.1f%10.1f%10.1f\n",fCutValue[i]*1.e9,zero,zero,-6.0,64.0);
8495a208 1694 }
e4b9dd05 1695
1696 else if (global){
aa067d5f 1697 cout << "SetCut for flag=" << &fCutFlag[i][0] << " value=" << fCutValue[i] << " not yet implemented!" << endl;
8495a208 1698 }
1699 else {
aa067d5f 1700 cout << "SetCut for flag=" << &fCutFlag[i][0] << " value=" << fCutValue[i] << " (material specific) not yet implemented!" << endl;
8495a208 1701 }
e4b9dd05 1702
e749c124 1703 } //end of loop over SetCut calls
e4b9dd05 1704
8495a208 1705// Add START and STOP card
e4b9dd05 1706 fprintf(pAliceInp,"START %10.1f\n",fEventsPerRun);
1707 fprintf(pAliceInp,"STOP \n");
2b8052b0 1708
e4b9dd05 1709
2b8052b0 1710// Close files
e4b9dd05 1711
2b8052b0 1712 fclose(pAliceCoreInp);
1713 fclose(pAliceFlukaMat);
1714 fclose(pAliceInp);
1715
8495a208 1716} // end of InitPhysics
1717
1718
1719//______________________________________________________________________________
1720void TFluka::SetMaxStep(Double_t)
1721{
1722// SetMaxStep is dummy procedure in TFluka !
1723 if (fVerbosityLevel >=3)
1724 cout << "SetMaxStep is dummy procedure in TFluka !" << endl;
1725}
1726
1727//______________________________________________________________________________
1728void TFluka::SetMaxNStep(Int_t)
1729{
1730// SetMaxNStep is dummy procedure in TFluka !
1731 if (fVerbosityLevel >=3)
1732 cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
1733}
1734
1735//______________________________________________________________________________
1736void TFluka::SetUserDecay(Int_t)
1737{
1738// SetUserDecay is dummy procedure in TFluka !
1739 if (fVerbosityLevel >=3)
1740 cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
1741}
1742
1743//
1744// dynamic properties
1745//
1746//______________________________________________________________________________
1747void TFluka::TrackPosition(TLorentzVector& position) const
1748{
1749// Return the current position in the master reference frame of the
1750// track being transported
1751// TRACKR.atrack = age of the particle
1752// TRACKR.xtrack = x-position of the last point
1753// TRACKR.ytrack = y-position of the last point
1754// TRACKR.ztrack = z-position of the last point
1755 Int_t caller = GetCaller();
fbf08100 1756 if (caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw
8495a208 1757 position.SetX(GetXsco());
1758 position.SetY(GetYsco());
1759 position.SetZ(GetZsco());
1760 position.SetT(TRACKR.atrack);
1761 }
1762 else if (caller == 4) { // mgdraw
1763 position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1764 position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1765 position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1766 position.SetT(TRACKR.atrack);
1767 }
1768 else if (caller == 5) { // sodraw
1769 position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1770 position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1771 position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1772 position.SetT(0);
1773 }
1774 else
1775 Warning("TrackPosition","position not available");
1776}
1777
1778//______________________________________________________________________________
1779void TFluka::TrackPosition(Double_t& x, Double_t& y, Double_t& z) const
1780{
1781// Return the current position in the master reference frame of the
1782// track being transported
1783// TRACKR.atrack = age of the particle
1784// TRACKR.xtrack = x-position of the last point
1785// TRACKR.ytrack = y-position of the last point
1786// TRACKR.ztrack = z-position of the last point
1787 Int_t caller = GetCaller();
fbf08100 1788 if (caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw
8495a208 1789 x = GetXsco();
1790 y = GetYsco();
1791 z = GetZsco();
1792 }
e8f0734b 1793 else if (caller == 4 || caller == 5) { // mgdraw, sodraw
8495a208 1794 x = TRACKR.xtrack[TRACKR.ntrack];
1795 y = TRACKR.ytrack[TRACKR.ntrack];
1796 z = TRACKR.ztrack[TRACKR.ntrack];
1797 }
1798 else
1799 Warning("TrackPosition","position not available");
1800}
1801
1802//______________________________________________________________________________
1803void TFluka::TrackMomentum(TLorentzVector& momentum) const
1804{
1805// Return the direction and the momentum (GeV/c) of the track
1806// currently being transported
1807// TRACKR.ptrack = momentum of the particle (not always defined, if
1808// < 0 must be obtained from etrack)
1809// TRACKR.cx,y,ztrck = direction cosines of the current particle
1810// TRACKR.etrack = total energy of the particle
1811// TRACKR.jtrack = identity number of the particle
1812// PAPROP.am[TRACKR.jtrack] = particle mass in gev
1813 Int_t caller = GetCaller();
1814 if (caller != 2) { // not eedraw
1815 if (TRACKR.ptrack >= 0) {
1816 momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
1817 momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
1818 momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
1819 momentum.SetE(TRACKR.etrack);
1820 return;
1821 }
1822 else {
1823 Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
1824 momentum.SetPx(p*TRACKR.cxtrck);
1825 momentum.SetPy(p*TRACKR.cytrck);
1826 momentum.SetPz(p*TRACKR.cztrck);
1827 momentum.SetE(TRACKR.etrack);
1828 return;
1829 }
1830 }
1831 else
1832 Warning("TrackMomentum","momentum not available");
1833}
1834
1835//______________________________________________________________________________
1836void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e) const
1837{
1838// Return the direction and the momentum (GeV/c) of the track
1839// currently being transported
1840// TRACKR.ptrack = momentum of the particle (not always defined, if
1841// < 0 must be obtained from etrack)
1842// TRACKR.cx,y,ztrck = direction cosines of the current particle
1843// TRACKR.etrack = total energy of the particle
1844// TRACKR.jtrack = identity number of the particle
1845// PAPROP.am[TRACKR.jtrack] = particle mass in gev
1846 Int_t caller = GetCaller();
1847 if (caller != 2) { // not eedraw
1848 if (TRACKR.ptrack >= 0) {
1849 px = TRACKR.ptrack*TRACKR.cxtrck;
1850 py = TRACKR.ptrack*TRACKR.cytrck;
1851 pz = TRACKR.ptrack*TRACKR.cztrck;
1852 e = TRACKR.etrack;
1853 return;
1854 }
1855 else {
1856 Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
1857 px = p*TRACKR.cxtrck;
1858 py = p*TRACKR.cytrck;
1859 pz = p*TRACKR.cztrck;
1860 e = TRACKR.etrack;
1861 return;
1862 }
1863 }
1864 else
1865 Warning("TrackMomentum","momentum not available");
1866}
1867
1868//______________________________________________________________________________
1869Double_t TFluka::TrackStep() const
1870{
1871// Return the length in centimeters of the current step
1872// TRACKR.ctrack = total curved path
1873 Int_t caller = GetCaller();
fbf08100 1874 if (caller == 11 || caller==12 || caller == 3 || caller == 6) //bxdraw,endraw,usdraw
8495a208 1875 return 0.0;
1876 else if (caller == 4) //mgdraw
1877 return TRACKR.ctrack;
1878 else
1879 return -1.0;
1880}
1881
1882//______________________________________________________________________________
1883Double_t TFluka::TrackLength() const
1884{
1885// TRACKR.cmtrck = cumulative curved path since particle birth
1886 Int_t caller = GetCaller();
fbf08100 1887 if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
8495a208 1888 return TRACKR.cmtrck;
1889 else
1890 return -1.0;
1891}
1892
1893//______________________________________________________________________________
1894Double_t TFluka::TrackTime() const
1895{
1896// Return the current time of flight of the track being transported
1897// TRACKR.atrack = age of the particle
1898 Int_t caller = GetCaller();
fbf08100 1899 if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
8495a208 1900 return TRACKR.atrack;
1901 else
1902 return -1;
1903}
1904
1905//______________________________________________________________________________
1906Double_t TFluka::Edep() const
1907{
1908// Energy deposition
1909// if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
1910// -->local energy deposition (the value and the point are not recorded in TRACKR)
1911// but in the variable "rull" of the procedure "endraw.cxx"
1912// if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
1913// -->no energy loss along the track
1914// if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
1915// -->energy loss distributed along the track
1916// TRACKR.dtrack = energy deposition of the jth deposition even
fbf08100 1917
1918 // If coming from bxdraw we have 2 steps of 0 length and 0 edep
1919 Int_t caller = GetCaller();
1920 if (caller == 11 || caller==12) return 0.0;
8495a208 1921 Double_t sum = 0;
1922 for ( Int_t j=0;j<TRACKR.mtrack;j++) {
1923 sum +=TRACKR.dtrack[j];
1924 }
1925 if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
1926 return fRull + sum;
1927 else {
1928 return sum;
1929 }
1930}
1931
1932//______________________________________________________________________________
1933Int_t TFluka::TrackPid() const
1934{
1935// Return the id of the particle transported
1936// TRACKR.jtrack = identity number of the particle
1937 Int_t caller = GetCaller();
1938 if (caller != 2) // not eedraw
1939 return PDGFromId(TRACKR.jtrack);
1940 else
1941 return -1000;
1942}
1943
1944//______________________________________________________________________________
1945Double_t TFluka::TrackCharge() const
1946{
1947// Return charge of the track currently transported
1948// PAPROP.ichrge = electric charge of the particle
1949// TRACKR.jtrack = identity number of the particle
1950 Int_t caller = GetCaller();
1951 if (caller != 2) // not eedraw
1952 return PAPROP.ichrge[TRACKR.jtrack+6];
1953 else
1954 return -1000.0;
1955}
1956
1957//______________________________________________________________________________
1958Double_t TFluka::TrackMass() const
1959{
1960// PAPROP.am = particle mass in GeV
1961// TRACKR.jtrack = identity number of the particle
1962 Int_t caller = GetCaller();
1963 if (caller != 2) // not eedraw
1964 return PAPROP.am[TRACKR.jtrack+6];
1965 else
1966 return -1000.0;
1967}
1968
1969//______________________________________________________________________________
1970Double_t TFluka::Etot() const
1971{
1972// TRACKR.etrack = total energy of the particle
1973 Int_t caller = GetCaller();
1974 if (caller != 2) // not eedraw
1975 return TRACKR.etrack;
1976 else
1977 return -1000.0;
1978}
1979
1980//
1981// track status
1982//
1983//______________________________________________________________________________
1984Bool_t TFluka::IsNewTrack() const
1985{
fbf08100 1986// Return true for the first call of Stepping()
fbf08100 1987 return fTrackIsNew;
8495a208 1988}
1989
1990//______________________________________________________________________________
1991Bool_t TFluka::IsTrackInside() const
1992{
1993// True if the track is not at the boundary of the current volume
1994// In Fluka a step is always inside one kind of material
1995// If the step would go behind the region of one material,
1996// it will be shortened to reach only the boundary.
1997// Therefore IsTrackInside() is always true.
1998 Int_t caller = GetCaller();
fbf08100 1999 if (caller == 11 || caller==12) // bxdraw
8495a208 2000 return 0;
2001 else
2002 return 1;
2003}
2004
2005//______________________________________________________________________________
2006Bool_t TFluka::IsTrackEntering() const
2007{
2008// True if this is the first step of the track in the current volume
2009
2010 Int_t caller = GetCaller();
2011 if (caller == 11) // bxdraw entering
2012 return 1;
2013 else return 0;
2014}
2015
2016//______________________________________________________________________________
2017Bool_t TFluka::IsTrackExiting() const
2018{
2019 Int_t caller = GetCaller();
2020 if (caller == 12) // bxdraw exiting
2021 return 1;
2022 else return 0;
2023}
2024
2025//______________________________________________________________________________
2026Bool_t TFluka::IsTrackOut() const
2027{
2028// True if the track is out of the setup
2029// means escape
2030// Icode = 14: escape - call from Kaskad
2031// Icode = 23: escape - call from Emfsco
2032// Icode = 32: escape - call from Kasneu
2033// Icode = 40: escape - call from Kashea
2034// Icode = 51: escape - call from Kasoph
2035 if (fIcode == 14 ||
2036 fIcode == 23 ||
2037 fIcode == 32 ||
2038 fIcode == 40 ||
2039 fIcode == 51) return 1;
2040 else return 0;
2041}
2042
2043//______________________________________________________________________________
2044Bool_t TFluka::IsTrackDisappeared() const
2045{
2046// means all inelastic interactions and decays
2047// fIcode from usdraw
2048 if (fIcode == 101 || // inelastic interaction
2049 fIcode == 102 || // particle decay
2050 fIcode == 214 || // in-flight annihilation
2051 fIcode == 215 || // annihilation at rest
2052 fIcode == 217 || // pair production
2053 fIcode == 221) return 1;
2054 else return 0;
2055}
2056
2057//______________________________________________________________________________
2058Bool_t TFluka::IsTrackStop() const
2059{
2060// True if the track energy has fallen below the threshold
2061// means stopped by signal or below energy threshold
2062// Icode = 12: stopping particle - call from Kaskad
2063// Icode = 15: time kill - call from Kaskad
2064// Icode = 21: below threshold, iarg=1 - call from Emfsco
2065// Icode = 22: below threshold, iarg=2 - call from Emfsco
2066// Icode = 24: time kill - call from Emfsco
2067// Icode = 31: below threshold - call from Kasneu
2068// Icode = 33: time kill - call from Kasneu
2069// Icode = 41: time kill - call from Kashea
2070// Icode = 52: time kill - call from Kasoph
2071 if (fIcode == 12 ||
2072 fIcode == 15 ||
2073 fIcode == 21 ||
2074 fIcode == 22 ||
2075 fIcode == 24 ||
2076 fIcode == 31 ||
2077 fIcode == 33 ||
2078 fIcode == 41 ||
2079 fIcode == 52) return 1;
2080 else return 0;
2081}
2082
2083//______________________________________________________________________________
2084Bool_t TFluka::IsTrackAlive() const
2085{
2086// means not disappeared or not out
2087 if (IsTrackDisappeared() || IsTrackOut() ) return 0;
2088 else return 1;
2089}
2090
2091//
2092// secondaries
2093//
2094
2095//______________________________________________________________________________
2096Int_t TFluka::NSecondaries() const
2097// Number of secondary particles generated in the current step
2098// FINUC.np = number of secondaries except light and heavy ions
2099// FHEAVY.npheav = number of secondaries for light and heavy secondary ions
2100{
2101 Int_t caller = GetCaller();
2102 if (caller == 6) // valid only after usdraw
2103 return FINUC.np + FHEAVY.npheav;
2104 else
2105 return 0;
2106} // end of NSecondaries
2107
2108//______________________________________________________________________________
2109void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
2110 TLorentzVector& position, TLorentzVector& momentum)
2111{
2112 Int_t caller = GetCaller();
2113 if (caller == 6) { // valid only after usdraw
2114 if (isec >= 0 && isec < FINUC.np) {
2115 particleId = PDGFromId(FINUC.kpart[isec]);
2116 position.SetX(fXsco);
2117 position.SetY(fYsco);
2118 position.SetZ(fZsco);
2119 position.SetT(TRACKR.atrack);
8495a208 2120 momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
2121 momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
2122 momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
2123 momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
2124 }
2125 else if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
2126 Int_t jsec = isec - FINUC.np;
2127 particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
2128 position.SetX(fXsco);
2129 position.SetY(fYsco);
2130 position.SetZ(fZsco);
2131 position.SetT(TRACKR.atrack);
8495a208 2132 momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
2133 momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
2134 momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
2135 if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
2136 momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
2137 else if (FHEAVY.tkheav[jsec] > 6)
2138 momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
2139 }
2140 else
2141 Warning("GetSecondary","isec out of range");
2142 }
2143 else
2144 Warning("GetSecondary","no secondaries available");
2145} // end of GetSecondary
2146
2147//______________________________________________________________________________
2148TMCProcess TFluka::ProdProcess(Int_t) const
2149// Name of the process that has produced the secondary particles
2150// in the current step
2151{
2152 const TMCProcess kIpNoProc = kPNoProcess;
2153 const TMCProcess kIpPDecay = kPDecay;
2154 const TMCProcess kIpPPair = kPPair;
2155// const TMCProcess kIpPPairFromPhoton = kPPairFromPhoton;
2156// const TMCProcess kIpPPairFromVirtualPhoton = kPPairFromVirtualPhoton;
2157 const TMCProcess kIpPCompton = kPCompton;
2158 const TMCProcess kIpPPhotoelectric = kPPhotoelectric;
2159 const TMCProcess kIpPBrem = kPBrem;
2160// const TMCProcess kIpPBremFromHeavy = kPBremFromHeavy;
2161// const TMCProcess kIpPBremFromElectronOrPositron = kPBremFromElectronOrPositron;
2162 const TMCProcess kIpPDeltaRay = kPDeltaRay;
2163// const TMCProcess kIpPMoller = kPMoller;
2164// const TMCProcess kIpPBhabha = kPBhabha;
2165 const TMCProcess kIpPAnnihilation = kPAnnihilation;
2166// const TMCProcess kIpPAnnihilInFlight = kPAnnihilInFlight;
2167// const TMCProcess kIpPAnnihilAtRest = kPAnnihilAtRest;
2168 const TMCProcess kIpPHadronic = kPHadronic;
2169 const TMCProcess kIpPMuonNuclear = kPMuonNuclear;
2170 const TMCProcess kIpPPhotoFission = kPPhotoFission;
2171 const TMCProcess kIpPRayleigh = kPRayleigh;
2172// const TMCProcess kIpPCerenkov = kPCerenkov;
2173// const TMCProcess kIpPSynchrotron = kPSynchrotron;
2174
2175 Int_t mugamma = TRACKR.jtrack == 7 || TRACKR.jtrack == 10 || TRACKR.jtrack == 11;
2176 if (fIcode == 102) return kIpPDecay;
2177 else if (fIcode == 104 || fIcode == 217) return kIpPPair;
2178// else if (fIcode == 104) return kIpPairFromPhoton;
2179// else if (fIcode == 217) return kIpPPairFromVirtualPhoton;
2180 else if (fIcode == 219) return kIpPCompton;
2181 else if (fIcode == 221) return kIpPPhotoelectric;
2182 else if (fIcode == 105 || fIcode == 208) return kIpPBrem;
2183// else if (fIcode == 105) return kIpPBremFromHeavy;
2184// else if (fIcode == 208) return kPBremFromElectronOrPositron;
2185 else if (fIcode == 103 || fIcode == 400) return kIpPDeltaRay;
2186 else if (fIcode == 210 || fIcode == 212) return kIpPDeltaRay;
2187// else if (fIcode == 210) return kIpPMoller;
2188// else if (fIcode == 212) return kIpPBhabha;
2189 else if (fIcode == 214 || fIcode == 215) return kIpPAnnihilation;
2190// else if (fIcode == 214) return kIpPAnnihilInFlight;
2191// else if (fIcode == 215) return kIpPAnnihilAtRest;
2192 else if (fIcode == 101) return kIpPHadronic;
2193 else if (fIcode == 101) {
2194 if (!mugamma) return kIpPHadronic;
2195 else if (TRACKR.jtrack == 7) return kIpPPhotoFission;
2196 else return kIpPMuonNuclear;
2197 }
2198 else if (fIcode == 225) return kIpPRayleigh;
2199// Fluka codes 100, 300 and 400 still to be investigasted
2200 else return kIpNoProc;
2201}
2202
2203//Int_t StepProcesses(TArrayI &proc) const
2204// Return processes active in the current step
2205//{
2206//ck = total energy of the particl ????????????????
2207//}
2208
2209
2210//______________________________________________________________________________
2211Int_t TFluka::VolId2Mate(Int_t id) const
2212{
2213//
2214// Returns the material number for a given volume ID
2215//
2216 return fGeom->VolId2Mate(id);
2217}
2218
2219//______________________________________________________________________________
2220const char* TFluka::VolName(Int_t id) const
2221{
2222//
2223// Returns the volume name for a given volume ID
2224//
2225 return fGeom->VolName(id);
2226}
2227
2228//______________________________________________________________________________
2229Int_t TFluka::VolId(const Text_t* volName) const
2230{
2231//
2232// Converts from volume name to volume ID.
2233// Time consuming. (Only used during set-up)
2234// Could be replaced by hash-table
2235//
2236 return fGeom->VolId(volName);
2237}
2238
2239//______________________________________________________________________________
2240Int_t TFluka::CurrentVolID(Int_t& copyNo) const
2241{
2242//
2243// Return the logical id and copy number corresponding to the current fluka region
2244//
2245 return fGeom->CurrentVolID(copyNo);
2246}
2247
2248//______________________________________________________________________________
2249Int_t TFluka::CurrentVolOffID(Int_t off, Int_t& copyNo) const
2250{
2251//
2252// Return the logical id and copy number of off'th mother
2253// corresponding to the current fluka region
2254//
2255 return fGeom->CurrentVolOffID(off, copyNo);
2256}
2257
2258//______________________________________________________________________________
2259const char* TFluka::CurrentVolName() const
2260{
2261//
2262// Return the current volume name
2263//
2264 return fGeom->CurrentVolName();
2265}
2266
2267//______________________________________________________________________________
2268const char* TFluka::CurrentVolOffName(Int_t off) const
2269{
2270//
2271// Return the volume name of the off'th mother of the current volume
2272//
2273 return fGeom->CurrentVolOffName(off);
2274}
2275
2276//______________________________________________________________________________
2277Int_t TFluka::CurrentMaterial(Float_t & /*a*/, Float_t & /*z*/,
2278 Float_t & /*dens*/, Float_t & /*radl*/, Float_t & /*absl*/) const
2279{
2280//
2281// Return the current medium number ??? what about material properties
2282//
2283 Int_t copy;
2284 Int_t id = TFluka::CurrentVolID(copy);
2285 Int_t med = TFluka::VolId2Mate(id);
2286 return med;
2287}
2288
2289//______________________________________________________________________________
2290void TFluka::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
2291{
2292// Transforms a position from the world reference frame
2293// to the current volume reference frame.
2294//
2295// Geant3 desription:
2296// ==================
2297// Computes coordinates XD (in DRS)
2298// from known coordinates XM in MRS
2299// The local reference system can be initialized by
2300// - the tracking routines and GMTOD used in GUSTEP
2301// - a call to GMEDIA(XM,NUMED)
2302// - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
2303// (inverse routine is GDTOM)
2304//
2305// If IFLAG=1 convert coordinates
2306// IFLAG=2 convert direction cosinus
2307//
2308// ---
2309 fGeom->Gmtod(xm,xd,iflag);
2310}
2311
2312//______________________________________________________________________________
2313void TFluka::Gmtod(Double_t* xm, Double_t* xd, Int_t iflag)
2314{
2315// Transforms a position from the world reference frame
2316// to the current volume reference frame.
2317//
2318// Geant3 desription:
2319// ==================
2320// Computes coordinates XD (in DRS)
2321// from known coordinates XM in MRS
2322// The local reference system can be initialized by
2323// - the tracking routines and GMTOD used in GUSTEP
2324// - a call to GMEDIA(XM,NUMED)
2325// - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
2326// (inverse routine is GDTOM)
2327//
2328// If IFLAG=1 convert coordinates
2329// IFLAG=2 convert direction cosinus
2330//
2331// ---
2332 fGeom->Gmtod(xm,xd,iflag);
2333}
2334
2335//______________________________________________________________________________
2336void TFluka::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
2337{
2338// Transforms a position from the current volume reference frame
2339// to the world reference frame.
2340//
2341// Geant3 desription:
2342// ==================
2343// Computes coordinates XM (Master Reference System
2344// knowing the coordinates XD (Detector Ref System)
2345// The local reference system can be initialized by
2346// - the tracking routines and GDTOM used in GUSTEP
2347// - a call to GSCMED(NLEVEL,NAMES,NUMBER)
2348// (inverse routine is GMTOD)
2349//
2350// If IFLAG=1 convert coordinates
2351// IFLAG=2 convert direction cosinus
2352//
2353// ---
2354 fGeom->Gdtom(xd,xm,iflag);
2355}
2356
2357//______________________________________________________________________________
2358void TFluka::Gdtom(Double_t* xd, Double_t* xm, Int_t iflag)
2359{
2360// Transforms a position from the current volume reference frame
2361// to the world reference frame.
2362//
2363// Geant3 desription:
2364// ==================
2365// Computes coordinates XM (Master Reference System
2366// knowing the coordinates XD (Detector Ref System)
2367// The local reference system can be initialized by
2368// - the tracking routines and GDTOM used in GUSTEP
2369// - a call to GSCMED(NLEVEL,NAMES,NUMBER)
2370// (inverse routine is GMTOD)
2371//
2372// If IFLAG=1 convert coordinates
2373// IFLAG=2 convert direction cosinus
2374//
2375// ---
2376 fGeom->Gdtom(xd,xm,iflag);
2377}
fbf08100 2378//______________________________________________________________________________
2379void TFluka::SetMreg(Int_t l)
2380{
2381// Set current fluka region
2382 fCurrentFlukaRegion = l;
2383 fGeom->SetMreg(l);
2384}
8495a208 2385
c7a186a7 2386
2387
2388