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