Mess from last merge cleaned up.
[u/mrichter/AliRoot.git] / TFluka / TFluka.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 /*
17 $Log$
18 Revision 1.8  2002/12/06 12:28:44  morsch
19 Region to media mapping corrected and improved.
20
21 Revision 1.7  2002/12/06 12:21:32  morsch
22 User stepping methods added (E. Futo)
23
24 Revision 1.6  2002/11/21 18:40:06  iglez2
25 Media handling added
26
27 Revision 1.5  2002/11/07 17:59:10  iglez2
28 Included the geometry through geant4_vmc/FLUGG
29
30 Revision 1.4  2002/11/04 16:00:46  iglez2
31 The conversion between ID and PDG now uses Fluka routines and arrays which is more consistent.
32
33 Revision 1.3  2002/10/22 15:12:14  alibrary
34 Introducing Riostream.h
35
36 Revision 1.2  2002/10/14 14:57:40  hristov
37 Merging the VirtualMC branch to the main development branch (HEAD)
38
39 Revision 1.1.2.8  2002/10/08 16:33:17  iglez2
40 LSOUIT is set to true before the second call to flukam.
41
42 Revision 1.1.2.7  2002/10/08 09:30:37  iglez2
43 Solved stupid missing ;
44
45 Revision 1.1.2.6  2002/10/07 13:40:22  iglez2
46 First implementations of the PDG <--> Fluka Id conversion routines
47
48 Revision 1.1.2.5  2002/09/26 16:26:03  iglez2
49 Added verbosity
50 Call to gAlice->Generator()->Generate()
51
52 Revision 1.1.2.4  2002/09/26 13:22:23  iglez2
53 Naive implementation of ProcessRun and ProcessEvent
54 Opening/Closing of input file (fInputFileName) with FORTRAN unit 5 before/after the first call to flukam inside Init()
55
56 Revision 1.1.2.3  2002/09/20 15:35:51  iglez2
57 Modification of LFDRTR. Value is passed to FLUKA !!!
58
59 Revision 1.1.2.2  2002/09/18 14:34:44  iglez2
60 Revised version with all pure virtual methods implemented
61
62 Revision 1.1.2.1  2002/07/24 08:49:41  alibrary
63 Adding TFluka to VirtualMC
64
65 Revision 1.1  2002/07/05 13:10:07  morsch
66 First commit of Fluka interface.
67
68 */
69
70 #include <Riostream.h>
71
72 #include "TClonesArray.h"
73 #include "TFluka.h"
74 #include "TCallf77.h"      //For the fortran calls
75 #include "Fdblprc.h"       //(DBLPRC) fluka common
76 #include "Fepisor.h"       //(EPISOR) fluka common
77 #include "Ffinuc.h"        //(FINUC) fluka common
78 #include "Fiounit.h"       //(IOUNIT) fluka common
79 #include "Fpaprop.h"       //(PAPROP) fluka common
80 #include "Fpart.h"         //(PART)   fluka common
81 #include "Ftrackr.h"       //(TRACKR) fluka common
82 #include "Fpaprop.h"       //(PAPROP) fluka common
83 #include "Ffheavy.h"       //(FHEAVY) fluka common
84
85 #include "TVirtualMC.h"
86 #include "TG4GeometryManager.h" //For the geometry management
87 #include "TG4DetConstruction.h" //For the detector construction
88
89 #include "FGeometryInit.hh"
90 #include "TLorentzVector.h"
91 #include "FlukaVolume.h"
92
93 // Fluka methods that may be needed.
94 #ifndef WIN32 
95 # define flukam  flukam_
96 # define fluka_openinp fluka_openinp_
97 # define fluka_closeinp fluka_closeinp_
98 # define mcihad mcihad_
99 # define mpdgha mpdgha_
100 #else 
101 # define flukam  FLUKAM
102 # define fluka_openinp FLUKA_OPENINP
103 # define fluka_closeinp FLUKA_CLOSEINP
104 # define mcihad MCIHAD
105 # define mpdgha MPDGHA
106 #endif
107
108 extern "C" 
109 {
110   //
111   // Prototypes for FLUKA functions
112   //
113   void type_of_call flukam(const int&);
114   void type_of_call fluka_openinp(const int&, DEFCHARA);
115   void type_of_call fluka_closeinp(const int&);
116   int  type_of_call mcihad(const int&);
117   int  type_of_call mpdgha(const int&);
118 }
119
120 //
121 // Class implementation for ROOT
122 //
123 ClassImp(TFluka)
124
125 //
126 //----------------------------------------------------------------------------
127 // TFluka constructors and destructors.
128 //____________________________________________________________________________ 
129 TFluka::TFluka()
130   :TVirtualMC(),
131    fVerbosityLevel(0),
132    fInputFileName(""),
133    fDetector(0),
134    fCurrentFlukaRegion(-1)
135
136   //
137   // Default constructor
138   //
139
140  
141 TFluka::TFluka(const char *title, Int_t verbosity)
142   :TVirtualMC("TFluka",title),
143    fVerbosityLevel(verbosity),
144    fInputFileName(""),
145    fDetector(0),
146    fCurrentFlukaRegion(-1)
147 {
148   if (fVerbosityLevel >=3)
149     cout << "==> TFluka::TFluka(" << title << ") constructor called." << endl;
150
151   
152   // create geometry manager
153   if (fVerbosityLevel >=2)
154     cout << "\t* Creating G4 Geometry manager..." << endl;
155   fGeometryManager = new TG4GeometryManager();
156   if (fVerbosityLevel >=2)
157     cout << "\t* Creating G4 Detector..." << endl;
158   fDetector = new TG4DetConstruction();
159   FGeometryInit* geominit = FGeometryInit::GetInstance();
160   if (geominit)
161     geominit->setDetConstruction(fDetector);
162   else {
163     cerr << "ERROR: Could not create FGeometryInit!" << endl;
164     cerr << "       Exiting!!!" << endl;
165     abort();
166   }
167
168   if (fVerbosityLevel >=3)
169     cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
170
171   fVolumeMediaMap = new TClonesArray("FlukaVolume",1000);
172   fNVolumes      = 0;
173   fMediaByRegion = 0;
174 }
175
176 TFluka::~TFluka() {
177   if (fVerbosityLevel >=3)
178     cout << "==> TFluka::~TFluka() destructor called." << endl;
179
180   delete fGeometryManager;
181   fVolumeMediaMap->Delete();
182   delete  fVolumeMediaMap;
183   
184
185   if (fVerbosityLevel >=3)
186     cout << "<== TFluka::~TFluka() destructor called." << endl;
187 }
188
189 //
190 //_____________________________________________________________________________
191 // TFluka control methods
192 //____________________________________________________________________________ 
193 void TFluka::Init() {
194   if (fVerbosityLevel >=3)
195     cout << "==> TFluka::Init() called." << endl;
196
197   if (fVerbosityLevel >=2)
198     cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F')
199          << ") in fluka..." << endl;
200   GLOBAL.lfdrtr = true;
201
202   if (fVerbosityLevel >=2)
203     cout << "\t* Opening file " << fInputFileName << endl;
204   const char* fname = fInputFileName;
205   fluka_openinp(lunin, PASSCHARA(fname));
206
207   if (fVerbosityLevel >=2)
208     cout << "\t* Calling flukam..." << endl;
209   flukam(1);
210
211   if (fVerbosityLevel >=2)
212     cout << "\t* Closing file " << fInputFileName << endl;
213   fluka_closeinp(lunin);
214
215   if (fVerbosityLevel >=3)
216     cout << "<== TFluka::Init() called." << endl;
217
218   FinishGeometry();
219
220 }
221
222 void TFluka::FinishGeometry() {
223 //
224 // Build-up table with region to medium correspondance
225 //
226     char tmp[5];
227     
228   if (fVerbosityLevel >=3)
229     cout << "==> TFluka::FinishGeometry() called." << endl;
230
231 //  fGeometryManager->Ggclos();
232
233   FGeometryInit* flugg = FGeometryInit::GetInstance();  
234   
235   fMediaByRegion = new Int_t[fNVolumes+2];
236   for (Int_t i = 0; i < fNVolumes; i++)
237   {
238       FlukaVolume* vol = dynamic_cast<FlukaVolume*>((*fVolumeMediaMap)[i]);
239       TString volName = vol->GetName();
240       Int_t   media   = vol->GetMedium();
241       printf("Finish Geometry: volName, media %d %s %d \n", i, volName.Data(), media);
242       strcpy(tmp, volName.Data());
243       tmp[4] = '\0';
244       flugg->SetMediumFromName(tmp, media);
245   }
246
247   flugg->BuildMediaMap();
248   
249   if (fVerbosityLevel >=3)
250     cout << "<== TFluka::FinishGeometry() called." << endl;
251
252
253 void TFluka::BuildPhysics() {
254   if (fVerbosityLevel >=3)
255     cout << "==> TFluka::BuildPhysics() called." << endl;
256
257
258   if (fVerbosityLevel >=3)
259     cout << "<== TFluka::BuildPhysics() called." << endl;
260 }  
261
262 void TFluka::ProcessEvent() {
263   if (fVerbosityLevel >=3)
264     cout << "==> TFluka::ProcessEvent() called." << endl;
265
266   if (fVerbosityLevel >=3)
267     cout << "<== TFluka::ProcessEvent() called." << endl;
268 }
269
270
271 void TFluka::ProcessRun(Int_t nevent) {
272   if (fVerbosityLevel >=3)
273     cout << "==> TFluka::ProcessRun(" << nevent << ") called." 
274          << endl;
275
276   if (fVerbosityLevel >=2) {
277     cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
278     cout << "\t* Calling flukam again..." << endl;
279   }
280   fApplication->GeneratePrimaries();
281   EPISOR.lsouit = true;
282   flukam(1);
283
284   if (fVerbosityLevel >=3)
285     cout << "<== TFluka::ProcessRun(" << nevent << ") called." 
286          << endl;
287 }
288
289 //_____________________________________________________________________________
290 // methods for building/management of geometry
291 //____________________________________________________________________________ 
292 // functions from GCONS 
293 void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,  
294                     Float_t &dens, Float_t &radl, Float_t &absl,
295                     Float_t* ubuf, Int_t& nbuf) {
296 //
297   fGeometryManager->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf);
298
299
300 void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,  
301                     Double_t &dens, Double_t &radl, Double_t &absl,
302                     Double_t* ubuf, Int_t& nbuf) {
303 //
304   fGeometryManager->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf);
305
306
307 // detector composition
308 void TFluka::Material(Int_t& kmat, const char* name, Double_t a, 
309                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
310                       Float_t* buf, Int_t nwbuf) {
311 //
312   fGeometryManager
313     ->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf); 
314
315 void TFluka::Material(Int_t& kmat, const char* name, Double_t a, 
316                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
317                       Double_t* buf, Int_t nwbuf) {
318 //
319   fGeometryManager
320     ->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf); 
321
322
323 void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a, 
324                      Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
325 //
326    fGeometryManager
327      ->Mixture(kmat, name, a, z, dens, nlmat, wmat); 
328
329 void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a, 
330                      Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
331 //
332    fGeometryManager
333      ->Mixture(kmat, name, a, z, dens, nlmat, wmat); 
334
335
336 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat, 
337                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd, 
338                     Double_t stemax, Double_t deemax, Double_t epsil, 
339                     Double_t stmin, Float_t* ubuf, Int_t nbuf) { 
340   //
341   fGeometryManager
342     ->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, 
343              epsil, stmin, ubuf, nbuf);
344
345 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat, 
346                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd, 
347                     Double_t stemax, Double_t deemax, Double_t epsil, 
348                     Double_t stmin, Double_t* ubuf, Int_t nbuf) { 
349   //
350   fGeometryManager
351     ->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, 
352              epsil, stmin, ubuf, nbuf);
353
354
355 void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX, 
356                     Double_t thetaY, Double_t phiY, Double_t thetaZ, 
357                     Double_t phiZ) {
358 //                   
359   fGeometryManager
360     ->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ); 
361
362
363 void TFluka::Gstpar(Int_t itmed, const char *param, Double_t parval) {
364 //
365   fGeometryManager->Gstpar(itmed, param, parval); 
366 }    
367
368 // functions from GGEOM 
369 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,  
370                      Float_t *upar, Int_t np)  {
371 //
372 //  fVolumeMediaMap[TString(name)] = nmed;
373     TClonesArray &lvols = *fVolumeMediaMap;
374     new(lvols[fNVolumes++]) 
375         FlukaVolume(name, nmed);
376     return fGeometryManager->Gsvolu(name, shape, nmed, upar, np); 
377 }
378 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,  
379                      Double_t *upar, Int_t np)  {
380 //
381     TClonesArray &lvols = *fVolumeMediaMap;
382     new(lvols[fNVolumes++]) 
383         FlukaVolume(name, nmed);
384
385     return fGeometryManager->Gsvolu(name, shape, nmed, upar, np); 
386 }
387  
388 void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv, 
389                    Int_t iaxis) {
390 //
391     fGeometryManager->Gsdvn(name, mother, ndiv, iaxis); 
392
393
394 void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv, 
395                     Int_t iaxis, Double_t c0i, Int_t numed) {
396 //
397     TClonesArray &lvols = *fVolumeMediaMap;
398     new(lvols[fNVolumes++]) 
399         FlukaVolume(name, numed);
400     fGeometryManager->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed); 
401
402
403 void TFluka::Gsdvt(const char *name, const char *mother, Double_t step, 
404                    Int_t iaxis, Int_t numed, Int_t ndvmx) {
405 //      
406     TClonesArray &lvols = *fVolumeMediaMap;
407     new(lvols[fNVolumes++]) 
408         FlukaVolume(name, numed);               
409     fGeometryManager->Gsdvt(name, mother, step, iaxis, numed, ndvmx); 
410
411
412 void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step, 
413                     Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) { 
414 //
415     TClonesArray &lvols = *fVolumeMediaMap;
416     new(lvols[fNVolumes++]) 
417         FlukaVolume(name, numed);
418     fGeometryManager->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx); 
419
420
421 void TFluka::Gsord(const char *name, Int_t iax) {
422 //
423   fGeometryManager->Gsord(name, iax); 
424
425
426 void TFluka::Gspos(const char *name, Int_t nr, const char *mother,  
427                    Double_t x, Double_t y, Double_t z, Int_t irot, 
428                    const char *konly) {
429 //
430   fGeometryManager->Gspos(name, nr, mother, x, y, z, irot, konly); 
431
432
433 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,  
434                     Double_t x, Double_t y, Double_t z, Int_t irot,
435                     const char *konly, Float_t *upar, Int_t np)  {
436   //
437   fGeometryManager->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np); 
438
439 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,  
440                     Double_t x, Double_t y, Double_t z, Int_t irot,
441                     const char *konly, Double_t *upar, Int_t np)  {
442   //
443   fGeometryManager->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np); 
444
445
446 void TFluka::Gsbool(const char* onlyVolName, const char* manyVolName) {
447 //
448   fGeometryManager->Gsbool(onlyVolName, manyVolName);
449 }
450
451 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t *ppckov,
452                          Float_t *absco, Float_t *effic, Float_t *rindex) {
453 //
454   fGeometryManager->SetCerenkov(itmed, npckov, ppckov, absco, effic, rindex);
455 }  
456 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t *ppckov,
457                          Double_t *absco, Double_t *effic, Double_t *rindex) {
458 //
459   fGeometryManager->SetCerenkov(itmed, npckov, ppckov, absco, effic, rindex);
460 }  
461     
462 // Euclid
463 void TFluka::WriteEuclid(const char* fileName, const char* topVol, 
464                           Int_t number, Int_t nlevel) {
465 //
466   fGeometryManager->WriteEuclid(fileName, topVol, number, nlevel); 
467
468
469
470
471 //_____________________________________________________________________________
472 // methods needed by the stepping
473 //____________________________________________________________________________ 
474
475 Int_t TFluka::GetMedium() const {
476     FGeometryInit* flugg = FGeometryInit::GetInstance();  
477     return flugg->GetMedium(fCurrentFlukaRegion);
478 }
479
480
481
482 //____________________________________________________________________________ 
483 // ID <--> PDG transformations
484 //_____________________________________________________________________________
485 Int_t TFluka::IdFromPDG(Int_t pdg) const 
486 {
487   //
488   // Return Fluka code from PDG and pseudo ENDF code
489
490   // MCIHAD() goes from pdg to fluka internal.
491   Int_t intfluka = mcihad(pdg);
492   // KPTOIP array goes from internal to official
493   return GetFlukaKPTOIP(intfluka);
494 }
495
496 Int_t TFluka::PDGFromId(Int_t id) const 
497 {
498   //
499   // Return PDG code and pseudo ENDF code from Fluka code
500
501   //IPTOKP array goes from official to internal
502   Int_t intfluka = GetFlukaIPTOKP(id);
503   //MPKDHA() goes from internal to PDG
504   return mpdgha(intfluka);
505 }
506
507 //_____________________________________________________________________________
508 // methods for step management
509 //____________________________________________________________________________ 
510 //
511 // dynamic properties
512 //
513 void TFluka::TrackPosition(TLorentzVector& position) const
514 {
515 // Return the current position in the master reference frame of the
516 // track being transported
517 // TRACKR.atrack = age of the particle
518 // TRACKR.xtrack = x-position of the last point
519 // TRACKR.ytrack = y-position of the last point
520 // TRACKR.ztrack = z-position of the last point
521   position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
522   position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
523   position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
524   position.SetT(TRACKR.atrack);
525 }
526
527 void TFluka::TrackMomentum(TLorentzVector& momentum) const
528 {
529 // Return the direction and the momentum (GeV/c) of the track
530 // currently being transported
531 // TRACKR.ptrack = momentum of the particle (not always defined, if
532 //               < 0 must be obtained from etrack) 
533 // TRACKR.cx,y,ztrck = direction cosines of the current particle
534 // TRACKR.etrack = total energy of the particle
535 // TRACKR.jtrack = identity number of the particle
536 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
537   if (TRACKR.ptrack >= 0) {
538     momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
539     momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
540     momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
541     momentum.SetE(TRACKR.etrack);
542     return;
543   }
544   else {
545     Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack]*PAPROP.am[TRACKR.jtrack]);
546     momentum.SetPx(p*TRACKR.cxtrck);
547     momentum.SetPy(p*TRACKR.cytrck);
548     momentum.SetPz(p*TRACKR.cztrck);
549     momentum.SetE(TRACKR.etrack);
550     return;
551   }
552 }
553
554 Double_t TFluka::TrackStep() const
555 {
556 // Return the length in centimeters of the current step
557 // TRACKR.ctrack = total curved path
558     return TRACKR.ctrack;
559 }
560
561 Double_t TFluka::TrackLength() const
562 {
563 // It is wrong
564 // should be the sum of all steps starting from the beginning of the track
565 // for the time being returns only the length in centimeters of the current step
566     return TRACKR.ctrack;
567 }
568
569 Double_t TFluka::TrackTime() const
570 {
571 // Return the current time of flight of the track being transported
572 // TRACKR.atrack = age of the particle
573   return TRACKR.atrack;
574 }
575
576 Double_t TFluka::Edep() const
577 {
578 // Energy deposition
579 // if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
580 // -->local energy deposition (the value and the point are not recorded in TRACKR)
581 //    but in the variable "rull" of the procedure "endraw.cxx"
582 // if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
583 // -->no energy loss along the track
584 // if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
585 // -->energy loss distributed along the track
586 // TRACKR.dtrack = energy deposition of the jth deposition even
587   if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
588     return fRull;
589   else {
590     Double_t sum = 0;
591     for ( Int_t j=0;j<TRACKR.mtrack;j++) {
592       sum +=TRACKR.dtrack[j];  
593     }
594     return sum;
595   }
596 }
597
598 Int_t TFluka::TrackPid() const
599 {
600 // Return the id of the particle transported
601 // TRACKR.jtrack = identity number of the particle
602   return PDGFromId(TRACKR.jtrack);
603 }
604
605 Double_t TFluka::TrackCharge() const
606 {
607 // Return charge of the track currently transported
608 // PAPROP.ichrge = electric charge of the particle
609   return PAPROP.ichrge[TRACKR.jtrack+6];
610 }
611
612 Double_t TFluka::TrackMass() const
613 {
614 // PAPROP.am = particle mass in GeV
615   return PAPROP.am[TRACKR.jtrack+6];
616 }
617
618 Double_t TFluka::Etot() const
619 {
620 // TRACKR.etrack = total energy of the particle
621   return TRACKR.etrack;
622 }
623
624 //
625 // track status
626 //
627 Bool_t   TFluka::IsNewTrack() const
628 {
629 // ???????????????,
630 // True if the track is not at the boundary of the current volume
631 // Not true in some cases in bxdraw - to be solved
632   return 1;
633 }
634
635 Bool_t   TFluka::IsTrackInside() const
636 {
637 // True if the track is not at the boundary of the current volume
638 // In Fluka a step is always inside one kind of material
639 // If the step would go behind the region of one material,
640 // it will be shortened to reach only the boundary.
641 // Therefore IsTrackInside() is always true.
642 // Not true in some cases in bxdraw - to be solved
643   return 1;
644 }
645
646 Bool_t   TFluka::IsTrackEntering() const
647 {
648 // True if this is the first step of the track in the current volume
649 // Boundary- (X) crossing
650 // Icode = 19: boundary crossing - call from Kaskad
651 // Icode = 29: boundary crossing - call from Emfsco
652 // Icode = 39: boundary crossing - call from Kasneu
653 // Icode = 49: boundary crossing - call from Kashea
654 // Icode = 59: boundary crossing - call from Kasoph
655   if (fIcode == 19 ||
656       fIcode == 29 ||
657       fIcode == 39 ||
658       fIcode == 49 ||
659       fIcode == 59) return 1;
660   else return 0;
661 }
662
663 Bool_t   TFluka::IsTrackExiting() const
664 {
665 // True if this is the last step of the track in the current volume
666 // Boundary- (X) crossing
667 // Icode = 19: boundary crossing - call from Kaskad
668 // Icode = 29: boundary crossing - call from Emfsco
669 // Icode = 39: boundary crossing - call from Kasneu
670 // Icode = 49: boundary crossing - call from Kashea
671 // Icode = 59: boundary crossing - call from Kasoph
672   if (fIcode == 19 ||
673       fIcode == 29 ||
674       fIcode == 39 ||
675       fIcode == 49 ||
676       fIcode == 59) return 1;
677   else return 0;
678 }
679
680 Bool_t   TFluka::IsTrackOut() const
681 {
682 // True if the track is out of the setup
683 // means escape
684 // Icode = 14: escape - call from Kaskad
685 // Icode = 23: escape - call from Emfsco
686 // Icode = 32: escape - call from Kasneu
687 // Icode = 40: escape - call from Kashea
688 // Icode = 51: escape - call from Kasoph
689   if (fIcode == 14 ||
690       fIcode == 23 ||
691       fIcode == 32 ||
692       fIcode == 40 ||
693       fIcode == 51) return 1;
694   else return 0;
695 }
696
697 Bool_t   TFluka::IsTrackDisappeared() const
698 {
699 // means all inelastic interactions and decays
700 // fIcode from usdraw
701   if (fIcode == 101 || // inelastic interaction
702       fIcode == 102 || // particle decay
703       fIcode == 214 || // in-flight annihilation
704       fIcode == 215 || // annihilation at rest
705       fIcode == 217 || // pair production
706       fIcode == 221) return 1;
707   else return 0;
708 }
709
710 Bool_t   TFluka::IsTrackStop() const
711 {
712 // True if the track energy has fallen below the threshold
713 // means stopped by signal or below energy threshold
714 // Icode = 12: stopping particle       - call from Kaskad
715 // Icode = 15: time kill               - call from Kaskad
716 // Icode = 21: below threshold, iarg=1 - call from Emfsco
717 // Icode = 22: below threshold, iarg=2 - call from Emfsco
718 // Icode = 24: time kill               - call from Emfsco
719 // Icode = 31: below threshold         - call from Kasneu
720 // Icode = 33: time kill               - call from Kasneu
721 // Icode = 41: time kill               - call from Kashea
722 // Icode = 52: time kill               - call from Kasoph
723   if (fIcode == 12 ||
724       fIcode == 15 ||
725       fIcode == 21 ||
726       fIcode == 22 ||
727       fIcode == 24 ||
728       fIcode == 31 ||
729       fIcode == 33 ||
730       fIcode == 41 ||
731       fIcode == 52) return 1;
732   else return 0;
733 }
734
735 Bool_t   TFluka::IsTrackAlive() const
736 {
737 // means not disappeared or not out
738   if (IsTrackDisappeared() || IsTrackOut() ) return 0;
739   else return 1;
740 }
741
742 //
743 // secondaries
744 //
745
746 Int_t TFluka::NSecondaries() const
747 // Number of secondary particles generated in the current step
748 // fIcode = 100 = elastic interaction
749 // fIcode = 101 = inelastic interaction
750 // fIcode = 102 = particle decay
751 // fIcode = 103 = delta ray
752 // fIcode = 104 = pair production
753 // fIcode = 105 = bremsstrahlung
754 {
755   if (fIcode >= 100 && fIcode <= 105)
756     return FINUC.np + FHEAVY.npheav;
757   else 
758     return -1;
759 }
760
761 void     TFluka::GetSecondary(Int_t isec, Int_t& particleId,
762                 TLorentzVector& position, TLorentzVector& momentum)
763 // 
764 // fIcode = 100 = elastic interaction
765 // fIcode = 101 = inelastic interaction
766 // fIcode = 102 = particle decay
767 // fIcode = 103 = delta ray
768 // fIcode = 104 = pair production
769 // fIcode = 105 = bremsstrahlung
770 {
771   if (fIcode >= 100 && fIcode <= 105) {
772     if (isec >= 0 && isec < FINUC.np) {
773       particleId = PDGFromId(FINUC.kpart[isec]);
774       position.SetX(fXsco);
775       position.SetY(fYsco);
776       position.SetZ(fZsco);
777       position.SetT(FINUC.agesec[isec]);
778       momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
779       momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
780       momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
781       momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
782     }
783     if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
784       Int_t jsec = isec - FINUC.np;
785       particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
786       position.SetX(fXsco);
787       position.SetY(fYsco);
788       position.SetZ(fZsco);
789       position.SetT(FHEAVY.agheav[jsec]);
790       momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
791       momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
792       momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
793       if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6) 
794         momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
795       else if (FHEAVY.tkheav[jsec] > 6)
796         momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
797     }
798   }
799 }
800
801 //TMCProcess ProdProcess(Int_t isec) const
802 // Name of the process that has produced the secondary particles
803 // in the current step
804 //{
805 // will come from FINUC when called from USDRAW
806 //}
807
808 //Int_t StepProcesses(TArrayI &proc) const
809 // Return processes active in the current step
810 //{
811 //ck = total energy of the particl ???????????????? 
812 //}
813
814
815 // ===============================================================
816 void TFluka::FutoTest() 
817 {
818   Int_t icode, mreg, newreg, particleId;
819 //  Int_t medium;
820   Double_t rull, xsco, ysco, zsco;
821   TLorentzVector position, momentum;
822   icode = GetIcode();
823   if (icode == 0) {
824     cout << " icode=" << icode << endl;
825     /*
826     cout << "TLorentzVector positionX=" << position.X()
827        << "positionY=" << position.Y()
828        << "positionZ=" << position.Z()
829        << "timeT=" << position.T() << endl;
830     cout << "TLorentzVector momentumX=" << momentum.X()
831        << "momentumY=" << momentum.Y()
832        << "momentumZ=" << momentum.Z()
833        << "energyE=" << momentum.E() << endl;
834     cout << "TrackPid=" << TrackPid() << endl;
835     */
836   }
837
838   else if (icode > 0 && icode <= 5) {
839 // mgdraw
840     mreg = GetMreg();
841 //    medium = GetMedium();
842     cout << " icode=" << icode
843          << " mreg=" << mreg
844 //       << " medium=" << medium
845          << endl;
846   TrackPosition(position);
847   TrackMomentum(momentum);
848   cout << "TLorentzVector positionX=" << position.X()
849        << "positionY=" << position.Y()
850        << "positionZ=" << position.Z()
851        << "timeT=" << position.T() << endl;
852   cout << "TLorentzVector momentumX=" << momentum.X()
853        << "momentumY=" << momentum.Y()
854        << "momentumZ=" << momentum.Z()
855        << "energyE=" << momentum.E() << endl;
856   cout << "TrackStep=" << TrackStep() << endl;
857   cout << "TrackLength=" << TrackLength() << endl;
858   cout << "TrackTime=" << TrackTime() << endl;
859   cout << "Edep=" << Edep() << endl;
860   cout << "TrackPid=" << TrackPid() << endl;
861   cout << "TrackCharge=" << TrackCharge() << endl;
862   cout << "TrackMass=" << TrackMass() << endl;
863   cout << "Etot=" << Etot() << endl;
864   cout << "IsNewTrack=" << IsNewTrack() << endl;
865   cout << "IsTrackInside=" << IsTrackInside() << endl;
866   cout << "IsTrackEntering=" << IsTrackEntering() << endl;
867   cout << "IsTrackExiting=" << IsTrackExiting() << endl;
868   cout << "IsTrackOut=" << IsTrackOut() << endl;
869   cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl;
870   cout << "IsTrackAlive=" << IsTrackAlive() << endl;
871   }
872
873   else if((icode >= 10 && icode <= 15) ||
874           (icode >= 20 && icode <= 24) ||
875           (icode >= 30 && icode <= 33) ||
876           (icode >= 40 && icode <= 41) ||
877           (icode >= 50 && icode <= 52)) {
878 // endraw
879     mreg = GetMreg();
880 //    medium = GetMedium();
881     rull = GetRull();
882     xsco = GetXsco();
883     ysco = GetYsco();
884     zsco = GetZsco();
885     cout << " icode=" << icode
886          << " mreg=" << mreg
887 //       << " medium=" << medium
888          << " rull=" << rull
889          << " xsco=" << xsco
890          << " ysco=" << ysco
891          << " zsco=" << zsco << endl;
892   TrackPosition(position);
893   TrackMomentum(momentum);
894   cout << "Edep=" << Edep() << endl;
895   cout << "Etot=" << Etot() << endl;
896   cout << "TrackPid=" << TrackPid() << endl;
897   cout << "TrackCharge=" << TrackCharge() << endl;
898   cout << "TrackMass=" << TrackMass() << endl;
899   cout << "IsTrackOut=" << IsTrackOut() << endl;
900   cout << "IsTrackDisappeared=" << IsTrackDisappeared() << endl;
901   cout << "IsTrackStop=" << IsTrackStop() << endl;
902   cout << "IsTrackAlive=" << IsTrackAlive() << endl;
903   }
904
905   else if((icode >= 100 && icode <= 105) ||
906            (icode == 208) ||
907            (icode == 210) ||
908            (icode == 212) ||
909            (icode >= 214 && icode <= 215) ||
910            (icode == 217) ||
911            (icode == 219) ||
912            (icode == 221) ||
913            (icode == 225) ||
914            (icode == 300) ||
915            (icode == 400)) {
916 // usdraw
917     mreg = GetMreg();
918 //    medium = GetMedium();
919     xsco = GetXsco();
920     ysco = GetYsco();
921     zsco = GetZsco();
922     cout << " icode=" << icode
923          << " mreg=" << mreg
924 //       << " medium=" << medium
925          << " xsco=" << xsco
926          << " ysco=" << ysco
927          << " zsco=" << zsco << endl;
928     cout << "TrackPid=" << TrackPid() << endl;
929     cout << "NSecondaries=" << NSecondaries() << endl;
930     for (Int_t isec=0; isec< NSecondaries(); isec++) {
931 //void     TFluka::GetSecondary(Int_t isec, Int_t& particleId,
932 //                 TLorentzVector& position, TLorentzVector& momentum)
933       TFluka::GetSecondary(isec, particleId, position, momentum);
934       cout << "TLorentzVector positionX=" << position.X()
935            << "positionY=" << position.Y()
936            << "positionZ=" << position.Z()
937            << "timeT=" << position.T() << endl;
938       cout << "TLorentzVector momentumX=" << momentum.X()
939            << "momentumY=" << momentum.Y()
940            << "momentumZ=" << momentum.Z()
941            << "energyE=" << momentum.E() << endl;
942       cout << "TrackPid=" << particleId << endl;
943
944     }
945   }
946
947   else if((icode == 19) ||
948           (icode == 29) ||
949           (icode == 39) ||
950           (icode == 49) ||
951           (icode == 59)) {
952     mreg = GetMreg();
953 //    medium = GetMedium();
954     newreg = GetNewreg();
955     xsco = GetXsco();
956     ysco = GetYsco();
957     zsco = GetZsco();
958     cout << " icode=" << icode
959          << " mreg=" << mreg
960 //       << " medium=" << medium
961          << " newreg=" << newreg
962          << " xsco=" << xsco
963          << " ysco=" << ysco
964          << " zsco=" << zsco << endl;
965   }
966 //
967 // ====================================================================
968 //
969
970   
971
972 } // end of FutoTest