b9d0a01d |
1 | #define METHODDEBUG |
2 | |
3 | // Fortran |
4 | #include "TCallf77.h" |
5 | |
6 | // Fluka commons |
7 | #include "Fdblprc.h" //(DBLPRC) fluka common |
8 | #include "Fdimpar.h" //(DIMPAR) fluka parameters |
9 | #include "Fepisor.h" //(EPISOR) fluka common |
10 | #include "Fstack.h" //(STACK) fluka common |
11 | #include "Fstars.h" //(STARS) fluka common |
12 | #include "Fbeam.h" //(BEAM) fluka common |
13 | #include "Fpaprop.h" //(PAPROP) fluka common |
14 | #include "Fltclcm.h" //(LTCLCM) fluka common |
15 | //#include "Fcaslim.h" //(CASLIM) fluka common |
16 | |
17 | //Virutal MC |
18 | #include "TFluka.h" |
19 | #include "TVirtualMCStack.h" |
45dc600a |
20 | #include "TVirtualMCApplication.h" |
b9d0a01d |
21 | #include "TParticle.h" |
22 | #include "TVector3.h" |
23 | |
24 | //Other |
eae0fe66 |
25 | #include <Riostream.h> |
b9d0a01d |
26 | |
27 | #ifndef WIN32 |
28 | # define source source_ |
29 | # define geocrs geocrs_ |
30 | # define georeg georeg_ |
31 | # define geohsm geohsm_ |
32 | # define soevsv soevsv_ |
33 | #else |
34 | # define source SOURCE |
35 | # define geocrs GEOCRS |
36 | # define georeg GEOREG |
37 | # define geohsm GEOHSM |
38 | # define soevsv SOEVSV |
39 | #endif |
40 | |
b9d0a01d |
41 | extern "C" { |
42 | // |
43 | // Prototypes for FLUKA functions |
44 | // |
45 | void type_of_call geocrs(Double_t &, Double_t &, Double_t &); |
46 | void type_of_call georeg(Double_t &, Double_t &, Double_t &, |
47 | Int_t &, Int_t &); |
48 | void type_of_call geohsm(Int_t &, Int_t &, Int_t &, Int_t &); |
49 | void type_of_call soevsv(); |
50 | /* |
51 | *----------------------------------------------------------------------* |
52 | * * |
53 | * Created on 07 january 1990 by Alfredo Ferrari & Paola Sala * |
54 | * Infn - Milan * |
55 | * * |
56 | * Last change on 21-jun-98 by Alfredo Ferrari * |
57 | * * |
58 | * C++ version on 27-sep-02 by Isidro Gonzalez * |
59 | * * |
60 | * This is just an example of a possible user written source routine. * |
61 | * note that the beam card still has some meaning - in the scoring the * |
62 | * maximum momentum used in deciding the binning is taken from the * |
63 | * beam momentum. Other beam card parameters are obsolete. * |
64 | * * |
65 | *----------------------------------------------------------------------*/ |
66 | |
67 | void source(Int_t& nomore) { |
68 | #ifdef METHODDEBUG |
69 | cout << "==> source(" << nomore << ")" << endl; |
70 | #endif |
71 | |
72 | cout << "\t* EPISOR.lsouit = " << (EPISOR.lsouit?'T':'F') << endl; |
73 | |
74 | static Bool_t lfirst = true; |
75 | /*======================================================================* |
76 | * * |
77 | * BASIC VERSION * |
78 | * * |
79 | *======================================================================*/ |
80 | nomore = 0; |
81 | /* +-------------------------------------------------------------------* |
82 | * | First call initializations:*/ |
83 | if (lfirst) { |
84 | |
85 | /*| *** The following 3 cards are mandatory ***/ |
86 | |
87 | EPISOR.tkesum = zerzer; |
88 | lfirst = false; |
89 | EPISOR.lussrc = true; |
90 | /*| *** User initialization ***/ |
91 | } |
92 | /* | |
93 | * +-------------------------------------------------------------------* |
94 | * Push one source particle to the stack. Note that you could as well |
95 | * push many but this way we reserve a maximum amount of space in the |
96 | * stack for the secondaries to be generated |
97 | */ |
98 | |
99 | // Get the pointer to the VMC |
100 | TVirtualMC* fluka = TFluka::GetMC(); |
101 | // Get the stack produced from the generator |
102 | TVirtualMCStack* cppstack = fluka->GetStack(); |
103 | //Get next particle |
45dc600a |
104 | if (STACK.lstack != 1) { |
105 | TVirtualMCApplication::Instance()->PostTrack(); |
106 | TVirtualMCApplication::Instance()->FinishPrimary(); |
107 | } |
b9d0a01d |
108 | Int_t itrack = -1; |
109 | TParticle* particle = cppstack->GetNextTrack(itrack); |
110 | |
111 | //Exit if itrack is negative (-1). Set lsouit to false to mark last track for |
112 | //this event |
113 | if (itrack<0) { |
114 | nomore = 1; |
115 | EPISOR.lsouit = false; |
116 | cout << "\t* EPISOR.lsouit = " << (EPISOR.lsouit?'T':'F') << endl; |
117 | cout << "\t* No more particles. Exiting..." << endl; |
118 | #ifdef METHODDEBUG |
119 | cout << "<== source(" << nomore << ")" << endl; |
120 | #endif |
121 | return; |
122 | } |
123 | |
124 | //Get some info about the particle and print it |
125 | TVector3 polarisation; |
126 | particle->GetPolarisation(polarisation); |
127 | cout << "\t* Particle " << itrack << " retrieved..." << endl; |
128 | cout << "\t\t+ Name = " << particle->GetName() << endl; |
129 | cout << "\t\t+ PDG/Fluka code = " << particle->GetPdgCode() |
130 | << " / " << fluka->IdFromPDG(particle->GetPdgCode()) << endl; |
131 | cout << "\t\t+ E = " << particle->Energy() << " GeV" << endl; |
132 | cout << "\t\t+ P = (" |
133 | << particle->Px() << " , " |
134 | << particle->Py() << " , " |
135 | << particle->Pz() << " ) --> " |
136 | << particle->P() << " GeV" << endl; |
137 | cout << "\t\t+ M = " << particle->GetMass() << " GeV" << endl; |
138 | cout << "\t\t+ Initial point = ( " |
139 | << particle->Vx() << " , " |
140 | << particle->Vy() << " , " |
141 | << particle->Vz() << " )" |
142 | << endl; |
143 | cout << "\t\t+ Polarisation = ( " |
144 | << polarisation.Px() << " , " |
145 | << polarisation.Py() << " , " |
146 | << polarisation.Pz() << " )" |
147 | << endl; |
148 | /* Lstack is the stack counter: of course any time source is called it |
149 | * must be =0 |
150 | */ |
45dc600a |
151 | |
b9d0a01d |
152 | STACK.lstack++; |
153 | cout << "\t* Storing particle parameters in the stack, lstack = " |
154 | << STACK.lstack << endl; |
155 | /* Wt is the weight of the particle*/ |
156 | STACK.wt[STACK.lstack] = oneone; |
157 | STARS.weipri += STACK.wt[STACK.lstack]; |
158 | /* Particle type (1=proton.....). Ijbeam is the type set by the BEAM |
159 | * card |
160 | */ |
161 | //STACK.ilo[STACK.lstack] = BEAM.ijbeam; |
162 | STACK.ilo[STACK.lstack] = fluka-> IdFromPDG(particle->GetPdgCode()); |
163 | /* From this point ..... |
164 | * Particle generation (1 for primaries) |
165 | */ |
166 | STACK.lo[STACK.lstack] = 1; |
167 | /* User dependent flag:*/ |
168 | STACK.louse[STACK.lstack] = 0; |
169 | /* User dependent spare variables:*/ |
148ba0b4 |
170 | Int_t ispr = 0; |
171 | for (ispr = 0; ispr < mkbmx1; ispr++) |
b9d0a01d |
172 | STACK.sparek[STACK.lstack][ispr] = zerzer; |
173 | /* User dependent spare flags:*/ |
148ba0b4 |
174 | for (ispr = 0; ispr < mkbmx2; ispr++) |
b9d0a01d |
175 | STACK.ispark[STACK.lstack][ispr] = 0; |
176 | /* Save the track number of the stack particle:*/ |
177 | STACK.ispark[STACK.lstack][mkbmx2-1] = STACK.lstack; |
178 | STACK.nparma++; |
179 | STACK.numpar[STACK.lstack] = STACK.nparma; |
180 | STACK.nevent[STACK.lstack] = 0; |
181 | STACK.dfnear[STACK.lstack] = +zerzer; |
182 | /* ... to this point: don't change anything |
183 | * Particle age (s) |
184 | */ |
185 | STACK.agestk[STACK.lstack] = +zerzer; |
186 | STACK.aknshr[STACK.lstack] = -twotwo; |
187 | /* Group number for "low" energy neutrons, set to 0 anyway*/ |
188 | STACK.igroup[STACK.lstack] = 0; |
189 | /* Kinetic energy of the particle (GeV)*/ |
190 | //STACK.tke[STACK.lstack] = |
191 | //sqrt( BEAM.pbeam*BEAM.pbeam + |
192 | // PAPROP.am[BEAM.ijbeam+6]*PAPROP.am[BEAM.ijbeam+6] ) |
193 | //- PAPROP.am[BEAM.ijbeam+6]; |
194 | STACK.tke[STACK.lstack] = particle->Energy() - particle->GetMass(); |
195 | |
196 | /* Particle momentum*/ |
197 | //STACK.pmom [STACK.lstack] = BEAM.pbeam; |
198 | STACK.pmom [STACK.lstack] = particle->P(); |
199 | |
200 | /* PMOM (lstack) = SQRT ( TKE (stack) * ( TKE (lstack) + TWOTWO |
201 | * & * AM (ILO(lstack)) ) ) |
202 | * Cosines (tx,ty,tz) |
203 | */ |
204 | //STACK.tx [STACK.lstack] = BEAM.tinx; |
205 | //STACK.ty [STACK.lstack] = BEAM.tiny; |
206 | //STACK.tz [STACK.lstack] = BEAM.tinz; |
207 | Double_t cosx = particle->Px()/particle->P(); |
208 | Double_t cosy = particle->Py()/particle->P(); |
209 | Double_t cosz = sqrt(oneone - cosx*cosx - cosy*cosy); |
210 | STACK.tx [STACK.lstack] = cosx; |
211 | STACK.ty [STACK.lstack] = cosy; |
212 | STACK.tz [STACK.lstack] = cosz; |
213 | |
214 | /* Polarization cosines: |
215 | */ |
216 | //STACK.txpol [STACK.lstack] = -twotwo; |
217 | //STACK.typol [STACK.lstack] = +zerzer; |
218 | //STACK.tzpol [STACK.lstack] = +zerzer; |
219 | if (polarisation.Mag()) { |
220 | Double_t cospolx = polarisation.Px()/polarisation.Mag(); |
221 | Double_t cospoly = polarisation.Py()/polarisation.Mag(); |
222 | Double_t cospolz = sqrt(oneone - cospolx*cospolx - cospoly*cospoly); |
223 | STACK.tx [STACK.lstack] = cospolx; |
224 | STACK.ty [STACK.lstack] = cospoly; |
225 | STACK.tz [STACK.lstack] = cospolz; |
226 | } |
227 | else { |
228 | STACK.txpol [STACK.lstack] = -twotwo; |
229 | STACK.typol [STACK.lstack] = +zerzer; |
230 | STACK.tzpol [STACK.lstack] = +zerzer; |
231 | } |
232 | |
233 | /* Particle coordinates*/ |
234 | //STACK.xa [STACK.lstack] = BEAM.xina; |
235 | //STACK.ya [STACK.lstack] = BEAM.yina; |
236 | //STACK.za [STACK.lstack] = BEAM.zina |
237 | //Vertext coordinates; |
238 | STACK.xa [STACK.lstack] = particle->Vx(); |
239 | STACK.ya [STACK.lstack] = particle->Vy(); |
240 | STACK.za [STACK.lstack] = particle->Vz(); |
241 | |
242 | // Some printout |
243 | cout << "\t* Particle information transfered to stack..." << endl; |
244 | |
245 | /* Calculate the total kinetic energy of the primaries: don't change*/ |
246 | Int_t st_ilo = STACK.ilo[STACK.lstack]; |
247 | if ( st_ilo != 0 ) |
248 | EPISOR.tkesum += |
249 | ((STACK.tke[STACK.lstack] + PAPROP.amdisc[st_ilo+6]) |
250 | * STACK.wt[STACK.lstack]); |
251 | else |
252 | EPISOR.tkesum += (STACK.tke[STACK.lstack] * STACK.wt[STACK.lstack]); |
253 | |
254 | /* Here we ask for the region number of the hitting point. |
255 | * NREG (LSTACK) = ... |
256 | * The following line makes the starting region search much more |
257 | * robust if particles are starting very close to a boundary: |
258 | */ |
259 | geocrs( STACK.tx[STACK.lstack], |
260 | STACK.ty[STACK.lstack], |
261 | STACK.tz[STACK.lstack] ); |
262 | Int_t idisc; |
263 | georeg ( STACK.xa[STACK.lstack], |
264 | STACK.ya[STACK.lstack], |
265 | STACK.za[STACK.lstack], |
266 | STACK.nreg[STACK.lstack], |
267 | idisc);//<-- dummy return variable not used |
268 | |
269 | /* Do not change these cards:*/ |
270 | Int_t igeohsm1 = 1; |
271 | Int_t igeohsm2 = -11; |
272 | geohsm ( STACK.nhspnt[STACK.lstack], igeohsm1, igeohsm2, LTCLCM.mlattc ); |
273 | STACK.nlattc[STACK.lstack] = LTCLCM.mlattc; |
274 | soevsv(); |
275 | |
276 | cout << "\t* EPISOR.lsouit = " << (EPISOR.lsouit?'T':'F') << endl; |
277 | cout << "\t* " << STACK.lstack << " particles in the event" << endl; |
45dc600a |
278 | TVirtualMCApplication::Instance()->PreTrack(); |
b9d0a01d |
279 | #ifdef METHODDEBUG |
280 | cout << "<== source(" << nomore << ")" << endl; |
281 | #endif |
282 | } |
283 | } |