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