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