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