In gudcay(): call ForceDecay() before Decay()
[u/mrichter/AliRoot.git] / TGeant3 / AliGeant3.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  2000/09/06 14:56:34  morsch
19 gudcay() method implemented.
20 Decays are performed using  the AliDecayer interface. The pointer to the instance of AliDecayer
21 is obtained from geant3 (will be gMC once it has been implemented for Geant4).
22
23 Revision 1.7  2000/07/13 16:19:10  fca
24 Mainly coding conventions + some small bug fixes
25
26 Revision 1.6  2000/07/11 18:24:59  fca
27 Coding convention corrections + few minor bug fixes
28
29 Revision 1.5  2000/05/20 14:49:48  fca
30 Call gdebug at the end of gustep
31
32 Revision 1.4  2000/04/26 10:17:32  fca
33 Changes in Lego for G4 compatibility
34
35 Revision 1.3  2000/04/07 11:12:35  fca
36 G4 compatibility changes
37
38 Revision 1.2  2000/02/29 19:11:17  fca
39 Move gucode into AliGeant3.cxx
40
41 Revision 1.1  2000/02/23 16:25:25  fca
42 AliVMC and AliGeant3 classes introduced
43 ReadEuclid moved from AliRun to AliModule
44
45 */
46
47 #include <TParticle.h>
48
49 #include "AliDecayer.h"
50 #include "AliGeant3.h"
51 #include "AliRun.h"
52 #include "TGeant3.h"
53 #include "AliCallf77.h" 
54
55 #ifndef WIN32
56
57 # define rxgtrak rxgtrak_
58 # define rxstrak rxstrak_
59 # define rxkeep  rxkeep_ 
60 # define rxouth  rxouth_
61 #else
62
63 # define rxgtrak RXGTRAK 
64 # define rxstrak RXSTRAK 
65 # define rxkeep  RXKEEP  
66 # define rxouth  RXOUTH
67 #endif
68
69 ClassImp(AliGeant3)
70
71 AliGeant3::AliGeant3(const char *title) : 
72   TGeant3(title) {}
73
74 //____________________________________________________________________________
75 void AliGeant3::FinishGeometry()
76 {
77   //
78   // Finalise geometry construction
79   //
80   TGeant3::FinishGeometry();
81   //Create the color table
82   SetColors();
83 }
84
85 //____________________________________________________________________________
86 void AliGeant3::Init()
87 {
88   //
89   //=================Create Materials and geometry
90   //
91   TObjArray *modules = gAlice->Modules();
92   TIter next(modules);
93   AliModule *detector;
94   while((detector = (AliModule*)next())) {
95     // Initialise detector materials and geometry
96     detector->CreateMaterials();
97     detector->CreateGeometry();
98     detector->BuildGeometry();
99     detector->Init();
100   }
101
102   //Terminate building of geometry
103   FinishGeometry();
104 }
105
106 //____________________________________________________________________________
107 void AliGeant3::ProcessRun(Int_t nevent)
108 {
109   //
110   // Process the run
111   //
112   Int_t todo = TMath::Abs(nevent);
113   for (Int_t i=0; i<todo; i++) {
114   // Process one run (one run = one event)
115      gAlice->BeginEvent();
116      ProcessEvent();
117      gAlice->FinishEvent();
118   }
119 }
120
121 //_____________________________________________________________________________
122 void AliGeant3::ProcessEvent()
123 {
124   //
125   // Process one event
126   //
127   Gtrigi();
128   Gtrigc();
129   Gtrig();
130 }
131
132 //_____________________________________________________________________________
133 void AliGeant3::SetColors()
134 {
135   //
136   // Set the colors for all the volumes
137   // this is done sequentially for all volumes
138   // based on the number of their medium
139   //
140   Int_t kv, icol;
141   Int_t jvolum=fGclink->jvolum;
142   //Int_t jtmed=fGclink->jtmed;
143   //Int_t jmate=fGclink->jmate;
144   Int_t nvolum=fGcnum->nvolum;
145   char name[5];
146   //
147   //    Now for all the volumes
148   for(kv=1;kv<=nvolum;kv++) {
149     //     Get the tracking medium
150     Int_t itm=Int_t (fZq[fZlq[jvolum-kv]+4]);
151     //     Get the material
152     //Int_t ima=Int_t (fZq[fZlq[jtmed-itm]+6]);
153     //     Get z
154     //Float_t z=fZq[fZlq[jmate-ima]+7];
155     //     Find color number
156     //icol = Int_t(z)%6+2;
157     //icol = 17+Int_t(z*150./92.);
158     //icol = kv%6+2;
159     icol = itm%6+2;
160     strncpy(name,(char*)&fZiq[jvolum+kv],4);
161     name[4]='\0';
162     Gsatt(name,"COLO",icol);
163   }
164 }
165
166 //_____________________________________________________________________________
167 //
168 //                 Interfaces to Fortran
169 //
170 //_____________________________________________________________________________
171
172 extern "C" void type_of_call  rxgtrak (Int_t &mtrack, Int_t &ipart, Float_t *pmom, 
173                                        Float_t &e, Float_t *vpos, Float_t *polar,
174                                        Float_t &tof)
175 {
176   //
177   //     Fetches next track from the ROOT stack for transport. Called by the
178   //     modified version of GTREVE.
179   //
180   //              Track number in the ROOT stack. If MTRACK=0 no
181   //      mtrack  more tracks are left in the stack to be
182   //              transported.
183   //      ipart   Particle code in the GEANT conventions.
184   //      pmom[3] Particle momentum in GeV/c
185   //      e       Particle energy in GeV
186   //      vpos[3] Particle position
187   //      tof     Particle time of flight in seconds
188   //
189   Int_t pdg;
190   gAlice->GetNextTrack(mtrack, pdg, pmom, e, vpos, polar, tof);
191   ipart = gMC->IdFromPDG(pdg);
192   mtrack++;
193 }
194
195 //_____________________________________________________________________________
196 extern "C" void type_of_call 
197 #ifndef WIN32
198 rxstrak (Int_t &keep, Int_t &parent, Int_t &ipart, Float_t *pmom, 
199                Float_t *vpos, Float_t &tof, const char* cmech, Int_t &ntr, const int cmlen)
200 #else
201 rxstrak (Int_t &keep, Int_t &parent, Int_t &ipart, Float_t *pmom,
202          Float_t *vpos, Float_t &tof, const char* cmech, const int cmlen,
203          Int_t &ntr)
204 #endif
205 {
206   //
207   //     Fetches next track from the ROOT stack for transport. Called by GUKINE
208   //     and GUSTEP.
209   //
210   //              Status of the track. If keep=0 the track is put
211   //      keep    on the ROOT stack but it is not fetched for
212   //              transport.
213   //      parent  Parent track. If parent=0 the track is a primary.
214   //              In GUSTEP the routine is normally called to store
215   //              secondaries generated by the current track whose
216   //              ROOT stack number is MTRACK (common SCKINE.
217   //      ipart   Particle code in the GEANT conventions.
218   //      pmom[3] Particle momentum in GeV/c
219   //      vpos[3] Particle position
220   //      tof     Particle time of flight in seconds
221   //
222   //      cmech   (CHARACTER*10) Particle origin. This field is user
223   //              defined and it is not used inside the GALICE code.
224   //      ntr     Number assigned to the particle in the ROOT stack.
225   //
226   char mecha[11];
227   Float_t polar[3]={0.,0.,0.};
228   for(int i=0; i<10 && i<cmlen; i++) mecha[i]=cmech[i];
229   mecha[10]=0;
230   Int_t pdg=gMC->PDGFromId(ipart);
231   gAlice->SetTrack(keep, parent-1, pdg, pmom, vpos, polar, tof, mecha, ntr);
232   ntr++;
233 }
234
235 //_____________________________________________________________________________
236 extern "C" void type_of_call  rxkeep(const Int_t &n)
237 {
238   if( NULL==gAlice ) exit(1);
239   
240   if( n<=0 || n>gAlice->Particles()->GetEntries() )
241     {
242       printf("  Bad index n=%d must be 0<n<=%d\n",
243              n,gAlice->Particles()->GetEntries());
244       exit(1);
245     }
246   
247   ((TParticle*)(gAlice->Particles()->UncheckedAt(n-1)))->SetBit(kKeepBit);
248 }
249
250 //_____________________________________________________________________________
251 extern "C" void type_of_call  rxouth ()
252 {
253   //
254   // Called by Gtreve at the end of each primary track
255   //
256   gAlice->FinishPrimary();
257 }
258
259
260 #ifndef WIN32
261 #  define gudigi gudigi_
262 #  define guhadr guhadr_
263 #  define guout  guout_
264 #  define guphad guphad_
265 #  define gudcay gudcay_
266 #  define guiget guiget_
267 #  define guinme guinme_
268 #  define guinti guinti_
269 #  define gunear gunear_
270 #  define guskip guskip_
271 #  define guview guview_
272 #  define gupara gupara_
273 #  define gudtim gudtim_
274 #  define guplsh guplsh_
275 #  define gutrev gutrev_
276 #  define gutrak gutrak_
277 #  define guswim guswim_
278 #  define gufld  gufld_
279 #  define gustep gustep_
280 #  define gukine gukine_
281 #  define uglast uglast_
282
283 #  define gheish gheish_
284 #  define flufin flufin_
285 #  define gfmfin gfmfin_
286 #  define gpghei gpghei_
287 #  define fldist fldist_
288 #  define gfmdis gfmdis_
289 #  define ghelx3 ghelx3_
290 #  define ghelix ghelix_
291 #  define grkuta grkuta_
292 #  define gtrack gtrack_
293 #  define gtreveroot gtreveroot_
294 #  define glast  glast_
295
296 #else
297 #  define gudigi GUDIGI
298 #  define guhadr GUHADR
299 #  define guout  GUOUT
300 #  define guphad GUPHAD
301 #  define gudcay GUDCAY
302 #  define guiget GUIGET
303 #  define guinme GUINME
304 #  define guinti GUINTI
305 #  define gunear GUNEAR
306 #  define guskip GUSKIP
307 #  define guview GUVIEW
308 #  define gupara GUPARA
309 #  define gudtim GUDTIM
310 #  define guplsh GUPLSH
311 #  define gutrev GUTREV
312 #  define gutrak GUTRAK
313 #  define guswim GUSWIM
314 #  define gufld  GUFLD
315 #  define gustep GUSTEP
316 #  define gukine GUKINE
317 #  define uglast UGLAST
318
319 #  define gheish GHEISH
320 #  define flufin FLUFIN
321 #  define gfmfin GFMFIN
322 #  define gpghei GPGHEI
323 #  define fldist FLDIST
324 #  define gfmdis GFMDIS
325 #  define ghelx3 GHELX3
326 #  define ghelix GHELIX
327 #  define grkuta GRKUTA
328 #  define gtrack GTRACK
329 #  define gtreveroot GTREVEROOT
330 #  define glast  GLAST
331
332 #endif
333
334 extern "C" type_of_call void gheish();
335 extern "C" type_of_call void flufin();
336 extern "C" type_of_call void gfmfin();
337 extern "C" type_of_call void gpghei();
338 extern "C" type_of_call void fldist();
339 extern "C" type_of_call void gfmdis();
340 extern "C" type_of_call void ghelx3(Float_t&, Float_t&, Float_t*, Float_t*);
341 extern "C" type_of_call void ghelix(Float_t&, Float_t&, Float_t*, Float_t*);
342 extern "C" type_of_call void grkuta(Float_t&, Float_t&, Float_t*, Float_t*);
343 extern "C" type_of_call void gtrack();
344 extern "C" type_of_call void gtreveroot();
345 extern "C" type_of_call void glast();
346
347 extern "C" type_of_call {
348
349 //______________________________________________________________________
350 void gudigi() 
351 {
352 //
353 //    ******************************************************************
354 //    *                                                                *
355 //    *       User routine to digitize one event                       *
356 //    *                                                                *
357 //    *    ==>Called by : GTRIG                                        *
358 //    *                                                                *
359 //    ******************************************************************
360
361 }
362
363
364 //______________________________________________________________________
365 void guhadr()
366 {
367 //
368 //    ******************************************************************
369 //    *                                                                *
370 //    *       User routine to generate one hadronic interaction        *
371 //    *                                                                *
372 //    *    ==>Called by : GTHADR,GTNEUT                                *
373 //    *                                                                *
374 //    ******************************************************************
375 //
376 //
377 //    ------------------------------------------------------------------
378 //
379       TGeant3* geant3 = (TGeant3*) gMC;
380       Int_t ihadr=geant3->Gcphys()->ihadr;
381       if (ihadr<4)       gheish();
382       else if (ihadr==4) flufin();
383       else               gfmfin();
384 }
385
386 //______________________________________________________________________
387 void guout()
388 {
389 //
390 //    ******************************************************************
391 //    *                                                                *
392 //    *       User routine called at the end of each event             *
393 //    *                                                                *
394 //    *    ==>Called by : GTRIG                                        *
395 //    *                                                                *
396 //    ******************************************************************
397 //
398 //
399 //    ------------------------------------------------------------------
400 //
401 }
402
403 //______________________________________________________________________
404 void guphad()
405 {
406 //
407 //    ******************************************************************
408 //    *                                                                *
409 //    *       User routine to compute Hadron. inter. probabilities     *
410 //    *                                                                *
411 //    *    ==>Called by : GTHADR,GTNEUT                                *
412 //    *                                                                *
413 //    ******************************************************************
414 //
415 //
416 //    ------------------------------------------------------------------
417 //
418       TGeant3* geant3 = (TGeant3*) gMC;
419       Int_t ihadr=geant3->Gcphys()->ihadr;
420       if (ihadr<4)       gpghei();
421       else if (ihadr==4) fldist();
422       else               gfmdis();
423 }
424
425 //______________________________________________________________________
426 void gudcay()
427 {
428 //
429 //    ******************************************************************
430 //    *                                                                *
431 //    *       User routine to decay particles                          *
432 //    *                                                                *
433 //    *    ==>Called by : GDECAY                                       *
434 //    *                                                                *
435 //    ******************************************************************
436 //
437 //
438 //    ------------------------------------------------------------------
439 //
440     
441     TGeant3* geant3=(TGeant3*) gMC;
442   // set decay table
443   gMC->Decayer()->ForceDecay();
444
445 // Initialize 4-momentum vector    
446     Int_t ipart = geant3->Gckine()->ipart;
447     TLorentzVector p;
448     
449     p[0]=geant3->Gctrak()->vect[3];
450     p[1]=geant3->Gctrak()->vect[4];
451     p[2]=geant3->Gctrak()->vect[5];
452     p[3]=geant3->Gctrak()->vect[6];    
453     
454 // Convert from geant to lund particle code
455     Int_t iplund=gMC->PDGFromId(ipart);
456 // Particle list
457     static TClonesArray *particles;
458     if(!particles) particles=new TClonesArray("TParticle",1000);
459 // Decay
460     gMC->Decayer()->Decay(iplund, &p);
461     
462 // Fetch Particles
463     Int_t np = geant3->Decayer()->ImportParticles(particles);
464     if (np <=1) return;
465     
466     for (Int_t i=0; i< np; i++) 
467     {
468         TParticle *  iparticle = (TParticle *) particles->At(i);
469         Int_t ks = iparticle->GetStatusCode();
470 // Final state particles only
471         if (ks < 1 || ks >  10) continue;
472 // Skip neutrinos
473         Int_t  kf=iparticle->GetPdgCode();
474         if (kf==12 || kf ==-12) continue;
475         if (kf==14 || kf ==-14) continue;
476         if (kf==16 || kf ==-16) continue;
477         
478         Int_t index=geant3->Gcking()->ngkine;
479 // Put particle on geant stack
480 // momentum vector
481
482         (geant3->Gcking()->gkin[index][0]) = iparticle->Px();
483         (geant3->Gcking()->gkin[index][1]) = iparticle->Py();
484         (geant3->Gcking()->gkin[index][2]) = iparticle->Pz();
485         (geant3->Gcking()->gkin[index][3]) = iparticle->Energy();
486         Int_t ilu = gMC->IdFromPDG(kf);
487
488 // particle type        
489         (geant3->Gcking()->gkin[index][4]) = Float_t(ilu);
490 // position
491         (geant3->Gckin3()->gpos[index][0]) = geant3->Gctrak()->vect[0];
492         (geant3->Gckin3()->gpos[index][1]) = geant3->Gctrak()->vect[1];
493         (geant3->Gckin3()->gpos[index][2]) = geant3->Gctrak()->vect[2];
494 // time of flight offset (mm)
495         (geant3->Gcking()->tofd[index])    = 0.;
496 //      (geant3->Gcking()->tofd[index])    = iparticle->T()/(10*kSpeedOfLight);
497 // increase stack counter
498         (geant3->Gcking()->ngkine)=index+1;
499     }
500 }
501
502 //______________________________________________________________________
503 void guiget(Int_t&, Int_t&, Int_t&)
504 {
505 //
506 //    ******************************************************************
507 //    *                                                                *
508 //    *       User routine for interactive control of GEANT            *
509 //    *                                                                *
510 //    *    ==>Called by : <GXINT>, GINCOM                              *
511 //    *                                                                *
512 //    ******************************************************************
513 //
514 //
515 //    ------------------------------------------------------------------
516 //
517 }
518
519 //______________________________________________________________________
520 void guinme(Float_t*, Int_t&, Float_t*, Int_t& IYES)
521 {
522 //
523 //    **********************************************
524 //    *                                            *
525 //    *    USER ROUTINE TO PROVIDE GINME FUNCTION  *
526 //    *    FOR ALL USER SHAPES IDENTIFIED BY THE   *
527 //    *    SHAPE NUMBER SH. POINT IS GIVEN IN X    *
528 //    *    THE PARAMETERS ARE GIVEN IN P. IYES IS  *
529 //    *    RETURNED 1 IF POINT IS IN, 0 IF POINT   *
530 //    *    IS OUT AND LESS THAN ZERO IF SHAPE      *
531 //    *    NUMBER IS NOT SUPPORTED.                *
532 //    *                                            *
533 //    *    ==>Called by : GINME                    *
534 //    *                                            *
535 //    **********************************************
536 //
537       IYES=-1;
538 }
539
540 //______________________________________________________________________
541 void guinti()
542 {
543 //
544 //    ******************************************************************
545 //    *                                                                *
546 //    *       User routine for interactive version                     *
547 //    *                                                                *
548 //    *    ==>Called by : <GXINT>,  GINTRI                             *
549 //    *                                                                *
550 //    ******************************************************************
551 //
552 //
553 //    ------------------------------------------------------------------
554 //
555 }
556
557 //______________________________________________________________________
558 void gunear(Int_t&, Int_t&, Float_t*, Int_t&)
559 {
560 //
561 //    ******************************************************************
562 //    *                                                                *
563 //    *    User search                                                 *
564 //    *       ISEARC to identify the given volume                      *
565 //    *       ICALL  to identify the calling routine                   *
566 //    *              1 GMEDIA like                                     *
567 //    *              2 GNEXT like                                      *
568 //    *       X      coordinates (+direction for ICALL=2)              *
569 //    *       JNEAR  address of default list of neighbours             *
570 //    *              (list to be overwriten by user)                   *
571 //    *                                                                *
572 //    *    Called by : GFTRAC, GINVOL, GTMEDI, GTNEXT, GNEXT, GMEDIA   *
573 //    *                                                                *
574 //    ******************************************************************
575 //
576 //
577 //    ------------------------------------------------------------------
578 //
579 }
580
581 //______________________________________________________________________
582 void guskip(Int_t& ISKIP)
583 {
584 //
585 //    ******************************************************************
586 //    *                                                                *
587 //    *   User routine to skip unwanted tracks                         *
588 //    *                                                                *
589 //    *   Called by : GSSTAK                                           *
590 //    *   Author    : F.Bruyant                                        *
591 //    *                                                                *
592 //    ******************************************************************
593 //
594 //
595 //    ------------------------------------------------------------------
596 //
597       ISKIP = 0;
598 }
599
600 //______________________________________________________________________
601 void guswim(Float_t& CHARGE, Float_t& STEP, Float_t* VECT, Float_t* VOUT)
602 {
603 //
604 //    ******************************************************************
605 //    *                                                                *
606 //    *       User routine to control tracking of one track            *
607 //    *       in a magnetic field                                      *
608 //    *                                                                *
609 //    *    ==>Called by : GTELEC,GTHADR,GTMUON                         *
610 //    *                                                                *
611 //    ******************************************************************
612 //
613 //
614 //    ------------------------------------------------------------------
615 //
616   TGeant3* geant3 = (TGeant3*) gMC;
617   Int_t ifield=geant3->Gctmed()->ifield;
618   Float_t fieldm=geant3->Gctmed()->fieldm;
619
620   if (ifield==3) {
621     Float_t fldcharge = fieldm*CHARGE;
622     ghelx3(fldcharge,STEP,VECT,VOUT);
623   }
624   else if (ifield==2) ghelix(CHARGE,STEP,VECT,VOUT);
625   else                grkuta(CHARGE,STEP,VECT,VOUT);
626 }
627
628 //______________________________________________________________________
629 void guview(Int_t&, Int_t&, DEFCHARD, Int_t& DEFCHARL)
630 {
631 //
632 //    ******************************************************************
633 //    *                                                                *
634 //    *       User routine for interactive version                     *
635 //    *                                                                *
636 //    *    ==>Called by : <GXINT>, GINC1                               *
637 //    *                                                                *
638 //    ******************************************************************
639 //
640 //
641 //    ------------------------------------------------------------------
642 //
643 }
644
645 //______________________________________________________________________
646 void gupara()
647 {
648 //
649 //    ******************************************************************
650 //    *                                                                *
651 //    *       User routine called every time a particle falls below    *
652 //    *       parametrization threshold. This routine should create    *
653 //    *       the parametrization stack, and, when this is full,       *
654 //    *       parametrize the shower and track the geantinos.          *
655 //    *                                                                *
656 //    *    ==>Called by : GTRACK                                       *
657 //    *                                                                *
658 //    ******************************************************************
659 //
660 //
661 //    ------------------------------------------------------------------
662 //
663 }
664
665 //______________________________________________________________________
666 Float_t gudtim(Float_t&, Float_t&, Int_t&, Int_t&)
667 {
668 //
669 //    ******************************************************************
670 //    *                                                                *
671 //    *       User function called by GCDRIF to return drift time      *
672 //    *                                                                *
673 //    *    ==>Called by : GCDRIF                                       *
674 //    *                                                                *
675 //    ******************************************************************
676 //
677 //
678 //    ------------------------------------------------------------------
679 //
680       return 0;
681 }
682
683
684 //______________________________________________________________________
685 Float_t guplsh(Int_t&, Int_t&)
686 {
687 //
688 //    ******************************************************************
689 //    *                                                                *
690 //    *                                                                *
691 //    *    ==>Called by : GLISUR                                       *
692 //    *                                                                *
693 //    ******************************************************************
694 //
695 //
696 //    ------------------------------------------------------------------
697 //
698 //
699 //*** By default this defines perfect smoothness
700       return 1;
701 }
702
703 //______________________________________________________________________
704 void gutrak()
705 {
706 //
707 //    ******************************************************************
708 //    *                                                                *
709 //    *       User routine to control tracking of one track            *
710 //    *                                                                *
711 //    *    ==>Called by : GTREVE                                       *
712 //    *                                                                *
713 //    ******************************************************************
714 //
715 //
716 //    ------------------------------------------------------------------
717 //
718      Int_t ndet = gAlice->Modules()->GetLast();
719      TObjArray &dets = *gAlice->Modules();
720      AliModule *module;
721      Int_t i;
722
723      for(i=0; i<=ndet; i++)
724        if((module = (AliModule*)dets[i]))
725          module->PreTrack();
726
727      gtrack();
728
729      for(i=0; i<=ndet; i++)
730        if((module = (AliModule*)dets[i]))
731          module->PostTrack();
732 }
733
734 //______________________________________________________________________
735 void gutrev()
736 {
737 //
738 //    ******************************************************************
739 //    *                                                                *
740 //    *       User routine to control tracking of one event            *
741 //    *                                                                *
742 //    *    ==>Called by : GTRIG                                        *
743 //    *                                                                *
744 //    ******************************************************************
745 //
746 //
747 //    ------------------------------------------------------------------
748 //
749   gtreveroot();
750 }
751
752
753 //______________________________________________________________________
754 void gufld(Float_t *x, Float_t *b)
755 {
756       if(gAlice->Field()) {
757          gAlice->Field()->Field(x,b);
758       } else {
759          printf("No mag field defined!\n");
760          b[0]=b[1]=b[2]=0;
761       }
762 }
763
764 //______________________________________________________________________
765 void gustep()
766 {
767 //
768 //    ******************************************************************
769 //    *                                                                *
770 //    *       User routine called at the end of each tracking step     *
771 //    *       INWVOL is different from 0 when the track has reached    *
772 //    *              a volume boundary                                 *
773 //    *       ISTOP is different from 0 if the track has stopped       *
774 //    *                                                                *
775 //    *    ==>Called by : GTRACK                                       *
776 //    *                                                                *
777 //    ******************************************************************
778 //
779
780
781   TLorentzVector x;
782   Float_t r;
783   Int_t ipp, jk, id, nt;
784   Float_t polar[3]={0,0,0};
785   Float_t mom[3];
786   const char *kChproc;
787
788   
789   TGeant3* geant3 = (TGeant3*) gMC;
790   //     Stop particle if outside user defined tracking region 
791   gMC->TrackPosition(x);
792   r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
793   if (r > gAlice->TrackingRmax() || TMath::Abs(x[2]) > gAlice->TrackingZmax()) {
794         gMC->StopTrack();
795   }
796
797   // --- Add new created particles 
798   if (gMC->NSecondaries() > 0) {
799     kChproc=gMC->ProdProcess();
800     for (jk = 0; jk < geant3->Gcking()->ngkine; ++jk) {
801       ipp = Int_t (geant3->Gcking()->gkin[jk][4]+0.5);
802       // --- Skip neutrinos! 
803       if (ipp != 4) {
804         gAlice->SetTrack(1,gAlice->CurrentTrack(),gMC->PDGFromId(ipp), geant3->Gcking()->gkin[jk], 
805                          geant3->Gckin3()->gpos[jk], polar,geant3->Gctrak()->tofg, kChproc, nt);
806       }
807     }
808   }
809   // Cherenkov photons here
810   if ( geant3->Gckin2()->ngphot ) {
811     for (jk = 0; jk < geant3->Gckin2()->ngphot; ++jk) {
812       mom[0]=geant3->Gckin2()->xphot[jk][3]*geant3->Gckin2()->xphot[jk][6];
813       mom[1]=geant3->Gckin2()->xphot[jk][4]*geant3->Gckin2()->xphot[jk][6];
814       mom[2]=geant3->Gckin2()->xphot[jk][5]*geant3->Gckin2()->xphot[jk][6];
815       gAlice->SetTrack(1, gAlice->CurrentTrack(), gMC->PDGFromId(50),
816                        mom,                             //momentum
817                        geant3->Gckin2()->xphot[jk],     //position
818                        &geant3->Gckin2()->xphot[jk][7], //polarisation
819                        geant3->Gckin2()->xphot[jk][10], //time of flight
820                        "Cherenkov", nt);
821       }
822   }
823   // --- Particle leaving the setup ?
824   if (!gMC->IsTrackOut()) 
825     if ((id=gAlice->DetFromMate(geant3->Gctmed()->numed)) >= 0) gAlice->StepManager(id);
826
827   // --- Standard GEANT debug routine 
828   if(geant3->Gcflag()->idebug) geant3->Gdebug();
829 }
830
831 //______________________________________________________________________
832 void gukine ()
833 {
834 //
835 //    ******************************************************************
836 //    *                                                                *
837 //    *       Read or Generates Kinematics for primary tracks          *
838 //    *                                                                *
839 //    *    ==>Called by : GTRIG                                        *
840 //    *                                                                *
841 //    ******************************************************************
842 //
843 //
844 //    ------------------------------------------------------------------
845 //
846   gAlice->Generator()->Generate();
847 }
848
849
850 //______________________________________________________________________
851 void uglast()
852 {
853 //
854 //    ******************************************************************
855 //    *                                                                *
856 //    *       User routine called at the end of the run                *
857 //    *                                                                *
858 //    *    ==>Called by : GRUN                                         *
859 //    *                                                                *
860 //    ******************************************************************
861 //
862 //
863 }
864 }
865