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