Coding conventions.
[u/mrichter/AliRoot.git] / TFluka / source_ray.cxx
1 // Fortran 
2 #include "TCallf77.h"
3
4 // Fluka commons
5 #include "Fdblprc.h"  //(DBLPRC) fluka common
6 #include "Fdimpar.h"  //(DIMPAR) fluka parameters
7 #include "Fepisor.h"  //(EPISOR) fluka common
8 #include "Fstack.h"   //(STACK)  fluka common
9 #include "Fstars.h"   //(STARS)  fluka common
10 #include "Fbeam.h"    //(BEAM)   fluka common
11 #include "Fpaprop.h"  //(PAPROP) fluka common
12 #include "Fltclcm.h"  //(LTCLCM) fluka common
13 #include "Fpart.h"  
14 //#include "Fcaslim.h"  //(CASLIM) fluka common
15
16 #include "TVector3.h"
17 #include "TRandom.h"
18 #include "TParticle.h"
19 #include "../EVGEN/AliGenScan.h"
20 #include "AliStack.h"
21
22 //Other
23 #include <Riostream.h>
24
25 #ifndef WIN32
26 # define source source_
27 # define geocrs geocrs_
28 # define georeg georeg_
29 # define geohsm geohsm_
30 # define soevsv soevsv_
31 # define mcihad mcihad_
32 # define source_ray source_ray__
33 #else
34 # define source SOURCE
35 # define geocrs GEOCRS
36 # define georeg GEOREG
37 # define geohsm GEOHSM
38 # define soevsv SOEVSV
39 # define mcihad MCIHAD
40 # define source_ray SOURCE_RAY
41 #endif
42
43 extern "C" {
44   //
45   // Prototypes for FLUKA functions
46   //
47   void type_of_call geocrs(Double_t &, Double_t &, Double_t &);
48   void type_of_call georeg(Double_t &, Double_t &, Double_t &, 
49                            Int_t &, Int_t &);
50   void type_of_call geohsm(Int_t &, Int_t &, Int_t &, Int_t &);
51   void type_of_call soevsv();
52   int  type_of_call mcihad(const int&);
53  /*
54    *----------------------------------------------------------------------*
55    *                                                                      *
56    *     Created on 07 january 1990   by    Alfredo Ferrari & Paola Sala  *
57    *                                                   Infn - Milan       *
58    *                                                                      *
59    *     Last change on 21-jun-98     by    Alfredo Ferrari               *
60    *                                                                      *
61    *     C++ version on 27-sep-02     by    Isidro Gonzalez               *
62    *                                                                      *
63    *  This is just an example of a possible user written source routine.  *
64    *  note that the beam card still has some meaning - in the scoring the *
65    *  maximum momentum used in deciding the binning is taken from the     *
66    *  beam momentum.  Other beam card parameters are obsolete.            *
67    *                                                                      *
68    *----------------------------------------------------------------------*/
69
70   void source_ray(Int_t& nomore) {
71
72       static Bool_t lfirst       = true;
73       static AliGenScan* gener   = 0;
74       static AliStack* stack     = 0;   
75
76       nomore = 0;
77       TParticle* particle = 0;
78       Int_t itrack        = -1;
79       if (lfirst) {
80           EPISOR.tkesum = zerzer;
81           lfirst = false;
82           EPISOR.lussrc = true;
83 //
84 // The generator
85 //
86           gener  = new AliGenScan(0);
87           gener->SetRange(40, -200., 200., 40, -200., 200., 1, -700., -699.);
88           gener->SetPart(0);
89           gener->SetThetaRange(180., 180.);
90           gener->SetPhiRange(0.,0.);
91           
92           
93 // Stack 
94           stack = new AliStack(10000);
95           gener->SetStack(stack);
96           gener->Init();
97
98       } else {
99           //
100           // Generate event
101           stack->Reset();
102           gener->Generate();
103           Int_t npart = stack->GetNprimary();
104           for (Int_t part = 0; part < npart; part++) {
105               STACK.lstack++;
106               particle = stack->Particle(part);
107               
108               printf("Particle %5d %10s %10.3f %10.3f %10.3f \n", 
109                      STACK.lstack, particle->GetName(), 
110                      particle->Px(), particle->Py(), particle->Pz());
111               
112
113               
114               /* Wt is the weight of the particle*/
115               STACK.wt[STACK.lstack] = oneone;
116               STARS.weipri += STACK.wt[STACK.lstack];
117               
118               // ray
119               STACK.ilo[STACK.lstack] = 0;
120               /* From this point .....
121                * Particle generation (1 for primaries)
122                */
123               STACK.lo[STACK.lstack] = 1;
124               
125               /* User dependent flag:*/
126               STACK.louse[STACK.lstack] = 0;
127               
128               /* User dependent spare variables:*/
129               Int_t ispr = 0;
130               for (ispr = 0; ispr < mkbmx1; ispr++)
131                   STACK.sparek[STACK.lstack][ispr] = zerzer;
132               
133               /* User dependent spare flags:*/
134               for (ispr = 0; ispr < mkbmx2; ispr++)
135                   STACK.ispark[STACK.lstack][ispr] = 0;
136               
137               /* Save the track number of the stack particle:*/
138               STACK.ispark[STACK.lstack][mkbmx2-1] = itrack;
139               STACK.nparma++;
140               STACK.numpar[STACK.lstack] = STACK.nparma;
141               STACK.nevent[STACK.lstack] = 0;
142               STACK.dfnear[STACK.lstack] = +zerzer;
143               
144               /* Particle age (s)*/
145               STACK.agestk[STACK.lstack] = +zerzer;
146               STACK.aknshr[STACK.lstack] = -twotwo;
147               
148               /* Group number for "low" energy neutrons, set to 0 anyway*/
149               STACK.igroup[STACK.lstack] = 0;
150               
151               /* Kinetic energy */
152               STACK.tke[STACK.lstack] = particle->Energy() - particle->GetMass();
153               
154               
155               /* Particle momentum*/
156               STACK.pmom [STACK.lstack] = particle->P();
157               
158               /* Cosines (tx,ty,tz)*/
159               Double_t cosx = particle->Px()/particle->P();
160               Double_t cosy = particle->Py()/particle->P();
161               Double_t cosz = TMath::Sqrt(oneone - cosx*cosx - cosy*cosy);
162               if (particle->Pz() < 0.) cosz = -cosz;
163               
164               STACK.tx [STACK.lstack] = cosx;
165               STACK.ty [STACK.lstack] = cosy;
166               STACK.tz [STACK.lstack] = cosz;
167               
168               /* Polarization cosines:*/
169               STACK.txpol [STACK.lstack] = -twotwo;
170               STACK.typol [STACK.lstack] = +zerzer;
171               STACK.tzpol [STACK.lstack] = +zerzer;
172               
173               /* Particle coordinates*/
174               STACK.xa [STACK.lstack] =   particle->Vx();
175               STACK.ya [STACK.lstack] =   particle->Vy();
176               STACK.za [STACK.lstack] =   particle->Vz();
177               
178               printf("Particle Vertex %10.3f %10.3f %10.3f  \n",  
179                      STACK.xa [STACK.lstack],  STACK.ya [STACK.lstack],  STACK.za [STACK.lstack]);
180               
181               
182               
183               /*  Calculate the total kinetic energy of the primaries: don't change*/
184               Int_t st_ilo =  STACK.ilo[STACK.lstack];
185               if ( st_ilo != 0 )
186                   EPISOR.tkesum += 
187                       ((STACK.tke[STACK.lstack] + PAPROP.amdisc[st_ilo+6])
188                        * STACK.wt[STACK.lstack]);
189               else
190                   EPISOR.tkesum += (STACK.tke[STACK.lstack] * STACK.wt[STACK.lstack]);
191               
192               /*  Here we ask for the region number of the hitting point.
193                *     NREG (LSTACK) = ...
194                *  The following line makes the starting region search much more
195                *  robust if particles are starting very close to a boundary:
196                */
197               geocrs( STACK.tx[STACK.lstack], 
198                       STACK.ty[STACK.lstack], 
199                       STACK.tz[STACK.lstack] );
200               
201               Int_t idisc;
202               
203               georeg ( STACK.xa[STACK.lstack], 
204                        STACK.ya[STACK.lstack], 
205                        STACK.za[STACK.lstack],
206                        STACK.nreg[STACK.lstack], 
207                        idisc);//<-- dummy return variable not used
208               /*  Do not change these cards:*/
209               Int_t igeohsm1 = 1;
210               Int_t igeohsm2 = -11;
211               geohsm ( STACK.nhspnt[STACK.lstack], igeohsm1, igeohsm2, LTCLCM.mlattc );
212               STACK.nlattc[STACK.lstack] = LTCLCM.mlattc;
213               soevsv();
214           }
215       }
216   }
217 }
218
219