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