]>
Commit | Line | Data |
---|---|---|
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" | |
20 | #include "TParticle.h" | |
21 | #include "TVector3.h" | |
22 | ||
23 | //Other | |
eae0fe66 | 24 | #include <Riostream.h> |
b9d0a01d | 25 | |
26 | #ifndef WIN32 | |
27 | # define source source_ | |
28 | # define geocrs geocrs_ | |
29 | # define georeg georeg_ | |
30 | # define geohsm geohsm_ | |
31 | # define soevsv soevsv_ | |
32 | #else | |
33 | # define source SOURCE | |
34 | # define geocrs GEOCRS | |
35 | # define georeg GEOREG | |
36 | # define geohsm GEOHSM | |
37 | # define soevsv SOEVSV | |
38 | #endif | |
39 | ||
40 | ||
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 | |
104 | Int_t itrack = -1; | |
105 | TParticle* particle = cppstack->GetNextTrack(itrack); | |
106 | ||
107 | //Exit if itrack is negative (-1). Set lsouit to false to mark last track for | |
108 | //this event | |
109 | if (itrack<0) { | |
110 | nomore = 1; | |
111 | EPISOR.lsouit = false; | |
112 | cout << "\t* EPISOR.lsouit = " << (EPISOR.lsouit?'T':'F') << endl; | |
113 | cout << "\t* No more particles. Exiting..." << endl; | |
114 | #ifdef METHODDEBUG | |
115 | cout << "<== source(" << nomore << ")" << endl; | |
116 | #endif | |
117 | return; | |
118 | } | |
119 | ||
120 | //Get some info about the particle and print it | |
121 | TVector3 polarisation; | |
122 | particle->GetPolarisation(polarisation); | |
123 | cout << "\t* Particle " << itrack << " retrieved..." << endl; | |
124 | cout << "\t\t+ Name = " << particle->GetName() << endl; | |
125 | cout << "\t\t+ PDG/Fluka code = " << particle->GetPdgCode() | |
126 | << " / " << fluka->IdFromPDG(particle->GetPdgCode()) << endl; | |
127 | cout << "\t\t+ E = " << particle->Energy() << " GeV" << endl; | |
128 | cout << "\t\t+ P = (" | |
129 | << particle->Px() << " , " | |
130 | << particle->Py() << " , " | |
131 | << particle->Pz() << " ) --> " | |
132 | << particle->P() << " GeV" << endl; | |
133 | cout << "\t\t+ M = " << particle->GetMass() << " GeV" << endl; | |
134 | cout << "\t\t+ Initial point = ( " | |
135 | << particle->Vx() << " , " | |
136 | << particle->Vy() << " , " | |
137 | << particle->Vz() << " )" | |
138 | << endl; | |
139 | cout << "\t\t+ Polarisation = ( " | |
140 | << polarisation.Px() << " , " | |
141 | << polarisation.Py() << " , " | |
142 | << polarisation.Pz() << " )" | |
143 | << endl; | |
144 | /* Lstack is the stack counter: of course any time source is called it | |
145 | * must be =0 | |
146 | */ | |
147 | STACK.lstack++; | |
148 | cout << "\t* Storing particle parameters in the stack, lstack = " | |
149 | << STACK.lstack << endl; | |
150 | /* Wt is the weight of the particle*/ | |
151 | STACK.wt[STACK.lstack] = oneone; | |
152 | STARS.weipri += STACK.wt[STACK.lstack]; | |
153 | /* Particle type (1=proton.....). Ijbeam is the type set by the BEAM | |
154 | * card | |
155 | */ | |
156 | //STACK.ilo[STACK.lstack] = BEAM.ijbeam; | |
157 | STACK.ilo[STACK.lstack] = fluka-> IdFromPDG(particle->GetPdgCode()); | |
158 | /* From this point ..... | |
159 | * Particle generation (1 for primaries) | |
160 | */ | |
161 | STACK.lo[STACK.lstack] = 1; | |
162 | /* User dependent flag:*/ | |
163 | STACK.louse[STACK.lstack] = 0; | |
164 | /* User dependent spare variables:*/ | |
165 | for (Int_t ispr = 0; ispr < mkbmx1; ispr++) | |
166 | STACK.sparek[STACK.lstack][ispr] = zerzer; | |
167 | /* User dependent spare flags:*/ | |
168 | for (Int_t ispr = 0; ispr < mkbmx2; ispr++) | |
169 | STACK.ispark[STACK.lstack][ispr] = 0; | |
170 | /* Save the track number of the stack particle:*/ | |
171 | STACK.ispark[STACK.lstack][mkbmx2-1] = STACK.lstack; | |
172 | STACK.nparma++; | |
173 | STACK.numpar[STACK.lstack] = STACK.nparma; | |
174 | STACK.nevent[STACK.lstack] = 0; | |
175 | STACK.dfnear[STACK.lstack] = +zerzer; | |
176 | /* ... to this point: don't change anything | |
177 | * Particle age (s) | |
178 | */ | |
179 | STACK.agestk[STACK.lstack] = +zerzer; | |
180 | STACK.aknshr[STACK.lstack] = -twotwo; | |
181 | /* Group number for "low" energy neutrons, set to 0 anyway*/ | |
182 | STACK.igroup[STACK.lstack] = 0; | |
183 | /* Kinetic energy of the particle (GeV)*/ | |
184 | //STACK.tke[STACK.lstack] = | |
185 | //sqrt( BEAM.pbeam*BEAM.pbeam + | |
186 | // PAPROP.am[BEAM.ijbeam+6]*PAPROP.am[BEAM.ijbeam+6] ) | |
187 | //- PAPROP.am[BEAM.ijbeam+6]; | |
188 | STACK.tke[STACK.lstack] = particle->Energy() - particle->GetMass(); | |
189 | ||
190 | /* Particle momentum*/ | |
191 | //STACK.pmom [STACK.lstack] = BEAM.pbeam; | |
192 | STACK.pmom [STACK.lstack] = particle->P(); | |
193 | ||
194 | /* PMOM (lstack) = SQRT ( TKE (stack) * ( TKE (lstack) + TWOTWO | |
195 | * & * AM (ILO(lstack)) ) ) | |
196 | * Cosines (tx,ty,tz) | |
197 | */ | |
198 | //STACK.tx [STACK.lstack] = BEAM.tinx; | |
199 | //STACK.ty [STACK.lstack] = BEAM.tiny; | |
200 | //STACK.tz [STACK.lstack] = BEAM.tinz; | |
201 | Double_t cosx = particle->Px()/particle->P(); | |
202 | Double_t cosy = particle->Py()/particle->P(); | |
203 | Double_t cosz = sqrt(oneone - cosx*cosx - cosy*cosy); | |
204 | STACK.tx [STACK.lstack] = cosx; | |
205 | STACK.ty [STACK.lstack] = cosy; | |
206 | STACK.tz [STACK.lstack] = cosz; | |
207 | ||
208 | /* Polarization cosines: | |
209 | */ | |
210 | //STACK.txpol [STACK.lstack] = -twotwo; | |
211 | //STACK.typol [STACK.lstack] = +zerzer; | |
212 | //STACK.tzpol [STACK.lstack] = +zerzer; | |
213 | if (polarisation.Mag()) { | |
214 | Double_t cospolx = polarisation.Px()/polarisation.Mag(); | |
215 | Double_t cospoly = polarisation.Py()/polarisation.Mag(); | |
216 | Double_t cospolz = sqrt(oneone - cospolx*cospolx - cospoly*cospoly); | |
217 | STACK.tx [STACK.lstack] = cospolx; | |
218 | STACK.ty [STACK.lstack] = cospoly; | |
219 | STACK.tz [STACK.lstack] = cospolz; | |
220 | } | |
221 | else { | |
222 | STACK.txpol [STACK.lstack] = -twotwo; | |
223 | STACK.typol [STACK.lstack] = +zerzer; | |
224 | STACK.tzpol [STACK.lstack] = +zerzer; | |
225 | } | |
226 | ||
227 | /* Particle coordinates*/ | |
228 | //STACK.xa [STACK.lstack] = BEAM.xina; | |
229 | //STACK.ya [STACK.lstack] = BEAM.yina; | |
230 | //STACK.za [STACK.lstack] = BEAM.zina | |
231 | //Vertext coordinates; | |
232 | STACK.xa [STACK.lstack] = particle->Vx(); | |
233 | STACK.ya [STACK.lstack] = particle->Vy(); | |
234 | STACK.za [STACK.lstack] = particle->Vz(); | |
235 | ||
236 | // Some printout | |
237 | cout << "\t* Particle information transfered to stack..." << endl; | |
238 | ||
239 | /* Calculate the total kinetic energy of the primaries: don't change*/ | |
240 | Int_t st_ilo = STACK.ilo[STACK.lstack]; | |
241 | if ( st_ilo != 0 ) | |
242 | EPISOR.tkesum += | |
243 | ((STACK.tke[STACK.lstack] + PAPROP.amdisc[st_ilo+6]) | |
244 | * STACK.wt[STACK.lstack]); | |
245 | else | |
246 | EPISOR.tkesum += (STACK.tke[STACK.lstack] * STACK.wt[STACK.lstack]); | |
247 | ||
248 | /* Here we ask for the region number of the hitting point. | |
249 | * NREG (LSTACK) = ... | |
250 | * The following line makes the starting region search much more | |
251 | * robust if particles are starting very close to a boundary: | |
252 | */ | |
253 | geocrs( STACK.tx[STACK.lstack], | |
254 | STACK.ty[STACK.lstack], | |
255 | STACK.tz[STACK.lstack] ); | |
256 | Int_t idisc; | |
257 | georeg ( STACK.xa[STACK.lstack], | |
258 | STACK.ya[STACK.lstack], | |
259 | STACK.za[STACK.lstack], | |
260 | STACK.nreg[STACK.lstack], | |
261 | idisc);//<-- dummy return variable not used | |
262 | ||
263 | /* Do not change these cards:*/ | |
264 | Int_t igeohsm1 = 1; | |
265 | Int_t igeohsm2 = -11; | |
266 | geohsm ( STACK.nhspnt[STACK.lstack], igeohsm1, igeohsm2, LTCLCM.mlattc ); | |
267 | STACK.nlattc[STACK.lstack] = LTCLCM.mlattc; | |
268 | soevsv(); | |
269 | ||
270 | cout << "\t* EPISOR.lsouit = " << (EPISOR.lsouit?'T':'F') << endl; | |
271 | cout << "\t* " << STACK.lstack << " particles in the event" << endl; | |
272 | ||
273 | #ifdef METHODDEBUG | |
274 | cout << "<== source(" << nomore << ")" << endl; | |
275 | #endif | |
276 | } | |
277 | } |