1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.15 2001/03/20 06:28:49 alibrary
19 New detector loop split in 2
21 Revision 1.14 2000/12/20 08:39:39 fca
22 Support for Cerenkov and process list in Virtual MC
24 Revision 1.13 2000/11/30 07:12:54 alibrary
25 Introducing new Rndm and QA classes
27 Revision 1.12 2000/11/06 11:35:46 morsch
28 Call BuildGeometry() after Init() to be able to share common detector parameters.
30 Revision 1.11 2000/10/04 16:30:22 fca
31 Add include for exit()
33 Revision 1.10 2000/10/02 21:28:16 fca
34 Removal of useless dependecies via forward declarations
36 Revision 1.9 2000/09/12 14:36:17 morsch
37 In gudcay(): call ForceDecay() before Decay()
39 Revision 1.8 2000/09/06 14:56:34 morsch
40 gudcay() method implemented.
41 Decays are performed using the AliDecayer interface. The pointer to the instance of AliDecayer
42 is obtained from geant3 (will be gMC once it has been implemented for Geant4).
44 Revision 1.7 2000/07/13 16:19:10 fca
45 Mainly coding conventions + some small bug fixes
47 Revision 1.6 2000/07/11 18:24:59 fca
48 Coding convention corrections + few minor bug fixes
50 Revision 1.5 2000/05/20 14:49:48 fca
51 Call gdebug at the end of gustep
53 Revision 1.4 2000/04/26 10:17:32 fca
54 Changes in Lego for G4 compatibility
56 Revision 1.3 2000/04/07 11:12:35 fca
57 G4 compatibility changes
59 Revision 1.2 2000/02/29 19:11:17 fca
60 Move gucode into AliGeant3.cxx
62 Revision 1.1 2000/02/23 16:25:25 fca
63 AliVMC and AliGeant3 classes introduced
64 ReadEuclid moved from AliRun to AliModule
70 #include <TParticle.h>
71 #include <TStopwatch.h>
73 #include "AliDecayer.h"
74 #include "AliGeant3.h"
77 #include "AliCallf77.h"
78 #include "AliModule.h"
80 #include "AliGenerator.h"
84 # define rxgtrak rxgtrak_
85 # define rxouth rxouth_
89 # define rxgtrak RXGTRAK
90 # define rxouth RXOUTH
96 AliGeant3::AliGeant3(const char *title) :
99 //____________________________________________________________________________
100 void AliGeant3::FinishGeometry()
103 // Finalise geometry construction
105 TGeant3::FinishGeometry();
106 //Create the color table
110 //____________________________________________________________________________
111 void AliGeant3::Init()
114 //=================Create Materials and geometry
117 TObjArray *modules = gAlice->Modules();
120 printf("Geometry creation:\n");
121 while((detector = (AliModule*)next())) {
123 // Initialise detector materials and geometry
124 detector->CreateMaterials();
125 detector->CreateGeometry();
126 printf("%10s R:%.2fs C:%.2fs\n",
127 detector->GetName(),stw.RealTime(),stw.CpuTime());
129 //Terminate building of geometry
132 printf("Initialisation:\n");
134 while((detector = (AliModule*)next())) {
136 // Initialise detector and display geometry
138 detector->BuildGeometry();
139 printf("%10s R:%.2fs C:%.2fs\n",
140 detector->GetName(),stw.RealTime(),stw.CpuTime());
145 //____________________________________________________________________________
146 void AliGeant3::ProcessRun(Int_t nevent)
151 Int_t todo = TMath::Abs(nevent);
152 for (Int_t i=0; i<todo; i++) {
153 // Process one run (one run = one event)
154 gAlice->BeginEvent();
156 gAlice->FinishEvent();
160 //_____________________________________________________________________________
161 void AliGeant3::ProcessEvent()
171 //_____________________________________________________________________________
172 void AliGeant3::SetColors()
175 // Set the colors for all the volumes
176 // this is done sequentially for all volumes
177 // based on the number of their medium
180 Int_t jvolum=fGclink->jvolum;
181 //Int_t jtmed=fGclink->jtmed;
182 //Int_t jmate=fGclink->jmate;
183 Int_t nvolum=fGcnum->nvolum;
186 // Now for all the volumes
187 for(kv=1;kv<=nvolum;kv++) {
188 // Get the tracking medium
189 Int_t itm=Int_t (fZq[fZlq[jvolum-kv]+4]);
191 //Int_t ima=Int_t (fZq[fZlq[jtmed-itm]+6]);
193 //Float_t z=fZq[fZlq[jmate-ima]+7];
195 //icol = Int_t(z)%6+2;
196 //icol = 17+Int_t(z*150./92.);
199 strncpy(name,(char*)&fZiq[jvolum+kv],4);
201 Gsatt(name,"COLO",icol);
205 //_____________________________________________________________________________
207 // Interfaces to Fortran
209 //_____________________________________________________________________________
211 extern "C" void type_of_call rxgtrak (Int_t &mtrack, Int_t &ipart, Float_t *pmom,
212 Float_t &e, Float_t *vpos, Float_t *polar,
216 // Fetches next track from the ROOT stack for transport. Called by the
217 // modified version of GTREVE.
219 // Track number in the ROOT stack. If MTRACK=0 no
220 // mtrack more tracks are left in the stack to be
222 // ipart Particle code in the GEANT conventions.
223 // pmom[3] Particle momentum in GeV/c
224 // e Particle energy in GeV
225 // vpos[3] Particle position
226 // tof Particle time of flight in seconds
229 gAlice->GetNextTrack(mtrack, pdg, pmom, e, vpos, polar, tof);
230 ipart = gMC->IdFromPDG(pdg);
234 //_____________________________________________________________________________
235 extern "C" void type_of_call rxouth ()
238 // Called by Gtreve at the end of each primary track
240 gAlice->FinishPrimary();
243 //_____________________________________________________________________________
244 extern "C" void type_of_call rxinh ()
247 // Called by Gtreve at the beginning of each primary track
249 gAlice->BeginPrimary();
253 # define gudigi gudigi_
254 # define guhadr guhadr_
255 # define guout guout_
256 # define guphad guphad_
257 # define gudcay gudcay_
258 # define guiget guiget_
259 # define guinme guinme_
260 # define guinti guinti_
261 # define gunear gunear_
262 # define guskip guskip_
263 # define guview guview_
264 # define gupara gupara_
265 # define gudtim gudtim_
266 # define guplsh guplsh_
267 # define gutrev gutrev_
268 # define gutrak gutrak_
269 # define guswim guswim_
270 # define gufld gufld_
271 # define gustep gustep_
272 # define gukine gukine_
273 # define uglast uglast_
275 # define gheish gheish_
276 # define flufin flufin_
277 # define gfmfin gfmfin_
278 # define gpghei gpghei_
279 # define fldist fldist_
280 # define gfmdis gfmdis_
281 # define ghelx3 ghelx3_
282 # define ghelix ghelix_
283 # define grkuta grkuta_
284 # define gtrack gtrack_
285 # define gtreveroot gtreveroot_
286 # define glast glast_
289 # define gudigi GUDIGI
290 # define guhadr GUHADR
292 # define guphad GUPHAD
293 # define gudcay GUDCAY
294 # define guiget GUIGET
295 # define guinme GUINME
296 # define guinti GUINTI
297 # define gunear GUNEAR
298 # define guskip GUSKIP
299 # define guview GUVIEW
300 # define gupara GUPARA
301 # define gudtim GUDTIM
302 # define guplsh GUPLSH
303 # define gutrev GUTREV
304 # define gutrak GUTRAK
305 # define guswim GUSWIM
307 # define gustep GUSTEP
308 # define gukine GUKINE
309 # define uglast UGLAST
311 # define gheish GHEISH
312 # define flufin FLUFIN
313 # define gfmfin GFMFIN
314 # define gpghei GPGHEI
315 # define fldist FLDIST
316 # define gfmdis GFMDIS
317 # define ghelx3 GHELX3
318 # define ghelix GHELIX
319 # define grkuta GRKUTA
320 # define gtrack GTRACK
321 # define gtreveroot GTREVEROOT
326 extern "C" type_of_call void gheish();
327 extern "C" type_of_call void flufin();
328 extern "C" type_of_call void gfmfin();
329 extern "C" type_of_call void gpghei();
330 extern "C" type_of_call void fldist();
331 extern "C" type_of_call void gfmdis();
332 extern "C" type_of_call void ghelx3(Float_t&, Float_t&, Float_t*, Float_t*);
333 extern "C" type_of_call void ghelix(Float_t&, Float_t&, Float_t*, Float_t*);
334 extern "C" type_of_call void grkuta(Float_t&, Float_t&, Float_t*, Float_t*);
335 extern "C" type_of_call void gtrack();
336 extern "C" type_of_call void gtreveroot();
337 extern "C" type_of_call void glast();
339 extern "C" type_of_call {
341 //______________________________________________________________________
345 // ******************************************************************
347 // * User routine to digitize one event *
349 // * ==>Called by : GTRIG *
351 // ******************************************************************
356 //______________________________________________________________________
360 // ******************************************************************
362 // * User routine to generate one hadronic interaction *
364 // * ==>Called by : GTHADR,GTNEUT *
366 // ******************************************************************
369 // ------------------------------------------------------------------
371 TGeant3* geant3 = (TGeant3*) gMC;
372 Int_t ihadr=geant3->Gcphys()->ihadr;
373 if (ihadr<4) gheish();
374 else if (ihadr==4) flufin();
378 //______________________________________________________________________
382 // ******************************************************************
384 // * User routine called at the end of each event *
386 // * ==>Called by : GTRIG *
388 // ******************************************************************
391 // ------------------------------------------------------------------
395 //______________________________________________________________________
399 // ******************************************************************
401 // * User routine to compute Hadron. inter. probabilities *
403 // * ==>Called by : GTHADR,GTNEUT *
405 // ******************************************************************
408 // ------------------------------------------------------------------
410 TGeant3* geant3 = (TGeant3*) gMC;
411 Int_t ihadr=geant3->Gcphys()->ihadr;
412 if (ihadr<4) gpghei();
413 else if (ihadr==4) fldist();
417 //______________________________________________________________________
421 // ******************************************************************
423 // * User routine to decay particles *
425 // * ==>Called by : GDECAY *
427 // ******************************************************************
430 // ------------------------------------------------------------------
433 TGeant3* geant3=(TGeant3*) gMC;
435 gMC->Decayer()->ForceDecay();
437 // Initialize 4-momentum vector
438 Int_t ipart = geant3->Gckine()->ipart;
441 p[0]=geant3->Gctrak()->vect[3];
442 p[1]=geant3->Gctrak()->vect[4];
443 p[2]=geant3->Gctrak()->vect[5];
444 p[3]=geant3->Gctrak()->vect[6];
446 // Convert from geant to lund particle code
447 Int_t iplund=gMC->PDGFromId(ipart);
449 static TClonesArray *particles;
450 if(!particles) particles=new TClonesArray("TParticle",1000);
452 gMC->Decayer()->Decay(iplund, &p);
455 Int_t np = geant3->Decayer()->ImportParticles(particles);
458 for (Int_t i=0; i< np; i++)
460 TParticle * iparticle = (TParticle *) particles->At(i);
461 Int_t ks = iparticle->GetStatusCode();
462 // Final state particles only
463 if (ks < 1 || ks > 10) continue;
465 Int_t kf=iparticle->GetPdgCode();
466 if (kf==12 || kf ==-12) continue;
467 if (kf==14 || kf ==-14) continue;
468 if (kf==16 || kf ==-16) continue;
470 Int_t index=geant3->Gcking()->ngkine;
471 // Put particle on geant stack
474 (geant3->Gcking()->gkin[index][0]) = iparticle->Px();
475 (geant3->Gcking()->gkin[index][1]) = iparticle->Py();
476 (geant3->Gcking()->gkin[index][2]) = iparticle->Pz();
477 (geant3->Gcking()->gkin[index][3]) = iparticle->Energy();
478 Int_t ilu = gMC->IdFromPDG(kf);
481 (geant3->Gcking()->gkin[index][4]) = Float_t(ilu);
483 (geant3->Gckin3()->gpos[index][0]) = geant3->Gctrak()->vect[0];
484 (geant3->Gckin3()->gpos[index][1]) = geant3->Gctrak()->vect[1];
485 (geant3->Gckin3()->gpos[index][2]) = geant3->Gctrak()->vect[2];
486 // time of flight offset (mm)
487 (geant3->Gcking()->tofd[index]) = 0.;
488 // (geant3->Gcking()->tofd[index]) = iparticle->T()/(10*kSpeedOfLight);
489 // increase stack counter
490 (geant3->Gcking()->ngkine)=index+1;
494 //______________________________________________________________________
495 void guiget(Int_t&, Int_t&, Int_t&)
498 // ******************************************************************
500 // * User routine for interactive control of GEANT *
502 // * ==>Called by : <GXINT>, GINCOM *
504 // ******************************************************************
507 // ------------------------------------------------------------------
511 //______________________________________________________________________
512 void guinme(Float_t*, Int_t&, Float_t*, Int_t& IYES)
515 // **********************************************
517 // * USER ROUTINE TO PROVIDE GINME FUNCTION *
518 // * FOR ALL USER SHAPES IDENTIFIED BY THE *
519 // * SHAPE NUMBER SH. POINT IS GIVEN IN X *
520 // * THE PARAMETERS ARE GIVEN IN P. IYES IS *
521 // * RETURNED 1 IF POINT IS IN, 0 IF POINT *
522 // * IS OUT AND LESS THAN ZERO IF SHAPE *
523 // * NUMBER IS NOT SUPPORTED. *
525 // * ==>Called by : GINME *
527 // **********************************************
532 //______________________________________________________________________
536 // ******************************************************************
538 // * User routine for interactive version *
540 // * ==>Called by : <GXINT>, GINTRI *
542 // ******************************************************************
545 // ------------------------------------------------------------------
549 //______________________________________________________________________
550 void gunear(Int_t&, Int_t&, Float_t*, Int_t&)
553 // ******************************************************************
556 // * ISEARC to identify the given volume *
557 // * ICALL to identify the calling routine *
560 // * X coordinates (+direction for ICALL=2) *
561 // * JNEAR address of default list of neighbours *
562 // * (list to be overwriten by user) *
564 // * Called by : GFTRAC, GINVOL, GTMEDI, GTNEXT, GNEXT, GMEDIA *
566 // ******************************************************************
569 // ------------------------------------------------------------------
573 //______________________________________________________________________
574 void guskip(Int_t& ISKIP)
577 // ******************************************************************
579 // * User routine to skip unwanted tracks *
581 // * Called by : GSSTAK *
582 // * Author : F.Bruyant *
584 // ******************************************************************
587 // ------------------------------------------------------------------
592 //______________________________________________________________________
593 void guswim(Float_t& CHARGE, Float_t& STEP, Float_t* VECT, Float_t* VOUT)
596 // ******************************************************************
598 // * User routine to control tracking of one track *
599 // * in a magnetic field *
601 // * ==>Called by : GTELEC,GTHADR,GTMUON *
603 // ******************************************************************
606 // ------------------------------------------------------------------
608 TGeant3* geant3 = (TGeant3*) gMC;
609 Int_t ifield=geant3->Gctmed()->ifield;
610 Float_t fieldm=geant3->Gctmed()->fieldm;
613 Float_t fldcharge = fieldm*CHARGE;
614 ghelx3(fldcharge,STEP,VECT,VOUT);
616 else if (ifield==2) ghelix(CHARGE,STEP,VECT,VOUT);
617 else grkuta(CHARGE,STEP,VECT,VOUT);
620 //______________________________________________________________________
621 void guview(Int_t&, Int_t&, DEFCHARD, Int_t& DEFCHARL)
624 // ******************************************************************
626 // * User routine for interactive version *
628 // * ==>Called by : <GXINT>, GINC1 *
630 // ******************************************************************
633 // ------------------------------------------------------------------
637 //______________________________________________________________________
641 // ******************************************************************
643 // * User routine called every time a particle falls below *
644 // * parametrization threshold. This routine should create *
645 // * the parametrization stack, and, when this is full, *
646 // * parametrize the shower and track the geantinos. *
648 // * ==>Called by : GTRACK *
650 // ******************************************************************
653 // ------------------------------------------------------------------
657 //______________________________________________________________________
658 Float_t gudtim(Float_t&, Float_t&, Int_t&, Int_t&)
661 // ******************************************************************
663 // * User function called by GCDRIF to return drift time *
665 // * ==>Called by : GCDRIF *
667 // ******************************************************************
670 // ------------------------------------------------------------------
676 //______________________________________________________________________
677 Float_t guplsh(Int_t&, Int_t&)
680 // ******************************************************************
683 // * ==>Called by : GLISUR *
685 // ******************************************************************
688 // ------------------------------------------------------------------
691 //*** By default this defines perfect smoothness
695 //______________________________________________________________________
699 // ******************************************************************
701 // * User routine to control tracking of one track *
703 // * ==>Called by : GTREVE *
705 // ******************************************************************
708 // ------------------------------------------------------------------
717 //______________________________________________________________________
721 // ******************************************************************
723 // * User routine to control tracking of one event *
725 // * ==>Called by : GTRIG *
727 // ******************************************************************
730 // ------------------------------------------------------------------
736 //______________________________________________________________________
737 void gufld(Float_t *x, Float_t *b)
739 if(gAlice->Field()) {
740 gAlice->Field()->Field(x,b);
742 printf("No mag field defined!\n");
747 //______________________________________________________________________
751 // ******************************************************************
753 // * User routine called at the end of each tracking step *
754 // * INWVOL is different from 0 when the track has reached *
755 // * a volume boundary *
756 // * ISTOP is different from 0 if the track has stopped *
758 // * ==>Called by : GTRACK *
760 // ******************************************************************
766 Int_t ipp, jk, id, nt;
767 Float_t polar[3]={0,0,0};
772 TGeant3* geant3 = (TGeant3*) gMC;
773 // Stop particle if outside user defined tracking region
774 gMC->TrackPosition(x);
775 r=TMath::Sqrt(x[0]*x[0]+x[1]*x[1]);
776 if (r > gAlice->TrackingRmax() || TMath::Abs(x[2]) > gAlice->TrackingZmax()) {
780 // --- Add new created particles
781 if (gMC->NSecondaries() > 0) {
782 pProc=gMC->ProdProcess(0);
783 for (jk = 0; jk < geant3->Gcking()->ngkine; ++jk) {
784 ipp = Int_t (geant3->Gcking()->gkin[jk][4]+0.5);
785 // --- Skip neutrinos!
787 gAlice->SetTrack(1,gAlice->CurrentTrack(),gMC->PDGFromId(ipp), geant3->Gcking()->gkin[jk],
788 geant3->Gckin3()->gpos[jk], polar,geant3->Gctrak()->tofg, pProc, nt);
792 // Cherenkov photons here
793 if ( geant3->Gckin2()->ngphot ) {
794 for (jk = 0; jk < geant3->Gckin2()->ngphot; ++jk) {
795 mom[0]=geant3->Gckin2()->xphot[jk][3]*geant3->Gckin2()->xphot[jk][6];
796 mom[1]=geant3->Gckin2()->xphot[jk][4]*geant3->Gckin2()->xphot[jk][6];
797 mom[2]=geant3->Gckin2()->xphot[jk][5]*geant3->Gckin2()->xphot[jk][6];
798 gAlice->SetTrack(1, gAlice->CurrentTrack(), gMC->PDGFromId(50),
800 geant3->Gckin2()->xphot[jk], //position
801 &geant3->Gckin2()->xphot[jk][7], //polarisation
802 geant3->Gckin2()->xphot[jk][10], //time of flight
806 // --- Particle leaving the setup ?
807 if (!gMC->IsTrackOut())
808 if ((id=gAlice->DetFromMate(geant3->Gctmed()->numed)) >= 0) gAlice->StepManager(id);
810 // --- Standard GEANT debug routine
811 if(geant3->Gcflag()->idebug) geant3->Gdebug();
814 //______________________________________________________________________
818 // ******************************************************************
820 // * Read or Generates Kinematics for primary tracks *
822 // * ==>Called by : GTRIG *
824 // ******************************************************************
827 // ------------------------------------------------------------------
829 gAlice->Generator()->Generate();
833 //______________________________________________________________________
837 // ******************************************************************
839 // * User routine called at the end of the run *
841 // * ==>Called by : GRUN *
843 // ******************************************************************