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