]>
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:*/ | |
148ba0b4 | 165 | Int_t ispr = 0; |
166 | for (ispr = 0; ispr < mkbmx1; ispr++) | |
b9d0a01d | 167 | STACK.sparek[STACK.lstack][ispr] = zerzer; |
168 | /* User dependent spare flags:*/ | |
148ba0b4 | 169 | for (ispr = 0; ispr < mkbmx2; ispr++) |
b9d0a01d | 170 | STACK.ispark[STACK.lstack][ispr] = 0; |
171 | /* Save the track number of the stack particle:*/ | |
172 | STACK.ispark[STACK.lstack][mkbmx2-1] = STACK.lstack; | |
173 | STACK.nparma++; | |
174 | STACK.numpar[STACK.lstack] = STACK.nparma; | |
175 | STACK.nevent[STACK.lstack] = 0; | |
176 | STACK.dfnear[STACK.lstack] = +zerzer; | |
177 | /* ... to this point: don't change anything | |
178 | * Particle age (s) | |
179 | */ | |
180 | STACK.agestk[STACK.lstack] = +zerzer; | |
181 | STACK.aknshr[STACK.lstack] = -twotwo; | |
182 | /* Group number for "low" energy neutrons, set to 0 anyway*/ | |
183 | STACK.igroup[STACK.lstack] = 0; | |
184 | /* Kinetic energy of the particle (GeV)*/ | |
185 | //STACK.tke[STACK.lstack] = | |
186 | //sqrt( BEAM.pbeam*BEAM.pbeam + | |
187 | // PAPROP.am[BEAM.ijbeam+6]*PAPROP.am[BEAM.ijbeam+6] ) | |
188 | //- PAPROP.am[BEAM.ijbeam+6]; | |
189 | STACK.tke[STACK.lstack] = particle->Energy() - particle->GetMass(); | |
190 | ||
191 | /* Particle momentum*/ | |
192 | //STACK.pmom [STACK.lstack] = BEAM.pbeam; | |
193 | STACK.pmom [STACK.lstack] = particle->P(); | |
194 | ||
195 | /* PMOM (lstack) = SQRT ( TKE (stack) * ( TKE (lstack) + TWOTWO | |
196 | * & * AM (ILO(lstack)) ) ) | |
197 | * Cosines (tx,ty,tz) | |
198 | */ | |
199 | //STACK.tx [STACK.lstack] = BEAM.tinx; | |
200 | //STACK.ty [STACK.lstack] = BEAM.tiny; | |
201 | //STACK.tz [STACK.lstack] = BEAM.tinz; | |
202 | Double_t cosx = particle->Px()/particle->P(); | |
203 | Double_t cosy = particle->Py()/particle->P(); | |
204 | Double_t cosz = sqrt(oneone - cosx*cosx - cosy*cosy); | |
205 | STACK.tx [STACK.lstack] = cosx; | |
206 | STACK.ty [STACK.lstack] = cosy; | |
207 | STACK.tz [STACK.lstack] = cosz; | |
208 | ||
209 | /* Polarization cosines: | |
210 | */ | |
211 | //STACK.txpol [STACK.lstack] = -twotwo; | |
212 | //STACK.typol [STACK.lstack] = +zerzer; | |
213 | //STACK.tzpol [STACK.lstack] = +zerzer; | |
214 | if (polarisation.Mag()) { | |
215 | Double_t cospolx = polarisation.Px()/polarisation.Mag(); | |
216 | Double_t cospoly = polarisation.Py()/polarisation.Mag(); | |
217 | Double_t cospolz = sqrt(oneone - cospolx*cospolx - cospoly*cospoly); | |
218 | STACK.tx [STACK.lstack] = cospolx; | |
219 | STACK.ty [STACK.lstack] = cospoly; | |
220 | STACK.tz [STACK.lstack] = cospolz; | |
221 | } | |
222 | else { | |
223 | STACK.txpol [STACK.lstack] = -twotwo; | |
224 | STACK.typol [STACK.lstack] = +zerzer; | |
225 | STACK.tzpol [STACK.lstack] = +zerzer; | |
226 | } | |
227 | ||
228 | /* Particle coordinates*/ | |
229 | //STACK.xa [STACK.lstack] = BEAM.xina; | |
230 | //STACK.ya [STACK.lstack] = BEAM.yina; | |
231 | //STACK.za [STACK.lstack] = BEAM.zina | |
232 | //Vertext coordinates; | |
233 | STACK.xa [STACK.lstack] = particle->Vx(); | |
234 | STACK.ya [STACK.lstack] = particle->Vy(); | |
235 | STACK.za [STACK.lstack] = particle->Vz(); | |
236 | ||
237 | // Some printout | |
238 | cout << "\t* Particle information transfered to stack..." << endl; | |
239 | ||
240 | /* Calculate the total kinetic energy of the primaries: don't change*/ | |
241 | Int_t st_ilo = STACK.ilo[STACK.lstack]; | |
242 | if ( st_ilo != 0 ) | |
243 | EPISOR.tkesum += | |
244 | ((STACK.tke[STACK.lstack] + PAPROP.amdisc[st_ilo+6]) | |
245 | * STACK.wt[STACK.lstack]); | |
246 | else | |
247 | EPISOR.tkesum += (STACK.tke[STACK.lstack] * STACK.wt[STACK.lstack]); | |
248 | ||
249 | /* Here we ask for the region number of the hitting point. | |
250 | * NREG (LSTACK) = ... | |
251 | * The following line makes the starting region search much more | |
252 | * robust if particles are starting very close to a boundary: | |
253 | */ | |
254 | geocrs( STACK.tx[STACK.lstack], | |
255 | STACK.ty[STACK.lstack], | |
256 | STACK.tz[STACK.lstack] ); | |
257 | Int_t idisc; | |
258 | georeg ( STACK.xa[STACK.lstack], | |
259 | STACK.ya[STACK.lstack], | |
260 | STACK.za[STACK.lstack], | |
261 | STACK.nreg[STACK.lstack], | |
262 | idisc);//<-- dummy return variable not used | |
263 | ||
264 | /* Do not change these cards:*/ | |
265 | Int_t igeohsm1 = 1; | |
266 | Int_t igeohsm2 = -11; | |
267 | geohsm ( STACK.nhspnt[STACK.lstack], igeohsm1, igeohsm2, LTCLCM.mlattc ); | |
268 | STACK.nlattc[STACK.lstack] = LTCLCM.mlattc; | |
269 | soevsv(); | |
270 | ||
271 | cout << "\t* EPISOR.lsouit = " << (EPISOR.lsouit?'T':'F') << endl; | |
272 | cout << "\t* " << STACK.lstack << " particles in the event" << endl; | |
273 | ||
274 | #ifdef METHODDEBUG | |
275 | cout << "<== source(" << nomore << ")" << endl; | |
276 | #endif | |
277 | } | |
278 | } |