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