]>
Commit | Line | Data |
---|---|---|
d10d0e42 | 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 "Fpart.h" | |
16 | //#include "Fcaslim.h" //(CASLIM) fluka common | |
17 | ||
18 | //Virutal MC | |
19 | #include "AliGenerator.h" | |
20 | #include "AliStack.h" | |
21 | #include "../THijing/AliGenHijing.h" | |
22 | ||
23 | #include "TVirtualMCStack.h" | |
24 | #include "TParticle.h" | |
25 | #include "TVector3.h" | |
26 | #include "TRandom.h" | |
27 | ||
28 | //Other | |
29 | #include <Riostream.h> | |
30 | ||
31 | #ifndef WIN32 | |
32 | # define source source_ | |
33 | # define geocrs geocrs_ | |
34 | # define georeg georeg_ | |
35 | # define geohsm geohsm_ | |
36 | # define soevsv soevsv_ | |
37 | # define mcihad mcihad_ | |
38 | # define source_pbpb source_pbpb__ | |
39 | #else | |
40 | # define source SOURCE | |
41 | # define geocrs GEOCRS | |
42 | # define georeg GEOREG | |
43 | # define geohsm GEOHSM | |
44 | # define soevsv SOEVSV | |
45 | # define mcihad MCIHAD | |
46 | # define source_pbpb SOURCE_PBPB | |
47 | #endif | |
48 | ||
49 | extern "C" { | |
50 | // | |
51 | // Prototypes for FLUKA functions | |
52 | // | |
53 | void type_of_call geocrs(Double_t &, Double_t &, Double_t &); | |
54 | void type_of_call georeg(Double_t &, Double_t &, Double_t &, | |
55 | Int_t &, Int_t &); | |
56 | void type_of_call geohsm(Int_t &, Int_t &, Int_t &, Int_t &); | |
57 | void type_of_call soevsv(); | |
58 | int type_of_call mcihad(const int&); | |
59 | /* | |
60 | *----------------------------------------------------------------------* | |
61 | * * | |
62 | * Created on 07 january 1990 by Alfredo Ferrari & Paola Sala * | |
63 | * Infn - Milan * | |
64 | * * | |
65 | * Last change on 21-jun-98 by Alfredo Ferrari * | |
66 | * * | |
67 | * C++ version on 27-sep-02 by Isidro Gonzalez * | |
68 | * * | |
69 | * This is just an example of a possible user written source routine. * | |
70 | * note that the beam card still has some meaning - in the scoring the * | |
71 | * maximum momentum used in deciding the binning is taken from the * | |
72 | * beam momentum. Other beam card parameters are obsolete. * | |
73 | * * | |
74 | *----------------------------------------------------------------------*/ | |
75 | ||
76 | void source_pbpb(Int_t& nomore) { | |
77 | ||
78 | static Bool_t lfirst = true; | |
79 | static Int_t part = 0; | |
80 | static AliGenHijing* gener = 0; | |
81 | static AliStack* stack = 0; | |
82 | static Int_t npart; | |
83 | ||
84 | nomore = 0; | |
85 | TParticle* particle; | |
86 | Int_t itrack = -1; | |
87 | if (lfirst) { | |
88 | EPISOR.tkesum = zerzer; | |
89 | lfirst = false; | |
90 | EPISOR.lussrc = true; | |
91 | // | |
92 | // The generator | |
93 | // | |
94 | gener = new AliGenHijing(-1); | |
95 | // beam energy | |
96 | gener->SetEnergyCMS(5500.); | |
97 | // reference frame | |
98 | gener->SetReferenceFrame("CMS"); | |
99 | // projectile | |
100 | gener->SetProjectile("A", 208, 82); | |
101 | gener->SetTarget ("A", 208, 82); | |
102 | gener->SetBoostLHC(0); | |
103 | gener->SetImpactParameterRange(0., 0.1); | |
104 | ||
105 | // tell hijing to keep the full parent child chain | |
106 | gener->KeepFullEvent(); | |
107 | // enable jet quenching | |
108 | gener->SetJetQuenching(1); | |
109 | // enable shadowing | |
110 | gener->SetShadowing(1); | |
111 | // neutral pion and heavy particle decays switched off | |
112 | gener->SetDecaysOff(1); | |
113 | // Don't track spectators | |
114 | gener->SetSpectators(0); | |
115 | // The particle stack | |
116 | stack = new AliStack(1000); | |
117 | gener->SetStack(stack); | |
118 | gener->Init(); | |
119 | // | |
120 | // Generate event | |
121 | stack->Reset(); | |
122 | gener->Generate(); | |
123 | npart = stack->GetNprimary(); | |
124 | printf("source_pbpb generated %d particles \n", npart); | |
125 | part = 0; | |
126 | } | |
127 | if (part >= npart) { | |
128 | nomore = 1; | |
129 | return; | |
130 | } | |
131 | ||
132 | // Vertex | |
133 | Float_t za = 0.; | |
134 | // Direction | |
135 | Float_t dir = 1.; | |
136 | Int_t ic = 1; | |
137 | ||
138 | while (ic != -1) { | |
139 | particle = stack->Particle(part++); | |
140 | ic = particle->GetFirstDaughter(); | |
141 | } | |
142 | ||
143 | ||
144 | Int_t pdg = particle->GetPdgCode(); | |
145 | Int_t intfluka = mcihad(pdg); | |
146 | Int_t ifl = GetFlukaKPTOIP(intfluka); | |
147 | TVector3 polarisation; | |
148 | particle->GetPolarisation(polarisation); | |
149 | ||
150 | STACK.lstack++; | |
151 | ||
152 | // printf("Particle %5d %5d %5d %10s %10.3f %10.3f %10.3f \n", part, pdg, ifl, | |
153 | // particle->GetName(), particle->Px(), particle->Py(), particle->Pz()); | |
154 | ||
155 | ||
156 | ||
157 | /* Wt is the weight of the particle*/ | |
158 | STACK.wt[STACK.lstack] = oneone; | |
159 | STARS.weipri += STACK.wt[STACK.lstack]; | |
160 | ||
161 | STACK.ilo[STACK.lstack] = ifl; | |
162 | /* From this point ..... | |
163 | * Particle generation (1 for primaries) | |
164 | */ | |
165 | STACK.lo[STACK.lstack] = 1; | |
166 | ||
167 | /* User dependent flag:*/ | |
168 | STACK.louse[STACK.lstack] = 0; | |
169 | ||
170 | /* User dependent spare variables:*/ | |
171 | Int_t ispr = 0; | |
172 | for (ispr = 0; ispr < mkbmx1; ispr++) | |
173 | STACK.sparek[STACK.lstack][ispr] = zerzer; | |
174 | ||
175 | /* User dependent spare flags:*/ | |
176 | for (ispr = 0; ispr < mkbmx2; ispr++) | |
177 | STACK.ispark[STACK.lstack][ispr] = 0; | |
178 | ||
179 | /* Save the track number of the stack particle:*/ | |
180 | STACK.ispark[STACK.lstack][mkbmx2-1] = itrack; | |
181 | STACK.nparma++; | |
182 | STACK.numpar[STACK.lstack] = STACK.nparma; | |
183 | STACK.nevent[STACK.lstack] = 0; | |
184 | STACK.dfnear[STACK.lstack] = +zerzer; | |
185 | ||
186 | /* Particle age (s)*/ | |
187 | STACK.agestk[STACK.lstack] = +zerzer; | |
188 | STACK.aknshr[STACK.lstack] = -twotwo; | |
189 | ||
190 | /* Group number for "low" energy neutrons, set to 0 anyway*/ | |
191 | STACK.igroup[STACK.lstack] = 0; | |
192 | ||
193 | /* Kinetic energy */ | |
194 | STACK.tke[STACK.lstack] = particle->Energy() - particle->GetMass(); | |
195 | ||
196 | ||
197 | /* Particle momentum*/ | |
198 | STACK.pmom [STACK.lstack] = particle->P(); | |
199 | ||
200 | /* Cosines (tx,ty,tz)*/ | |
201 | Double_t cosx = particle->Px()/particle->P(); | |
202 | Double_t cosy = particle->Py()/particle->P(); | |
203 | Double_t cosz = TMath::Sqrt(oneone - cosx*cosx - cosy*cosy); | |
204 | if (particle->Pz() < 0.) cosz = -cosz; | |
205 | cosz *= dir; | |
206 | ||
207 | STACK.tx [STACK.lstack] = cosx; | |
208 | STACK.ty [STACK.lstack] = cosy; | |
209 | STACK.tz [STACK.lstack] = cosz; | |
210 | ||
211 | /* Polarization cosines:*/ | |
212 | if (polarisation.Mag()) { | |
213 | Double_t cospolx = polarisation.Px()/polarisation.Mag(); | |
214 | Double_t cospoly = polarisation.Py()/polarisation.Mag(); | |
215 | Double_t cospolz = sqrt(oneone - cospolx*cospolx - cospoly*cospoly); | |
216 | STACK.tx [STACK.lstack] = cospolx; | |
217 | STACK.ty [STACK.lstack] = cospoly; | |
218 | STACK.tz [STACK.lstack] = cospolz; | |
219 | } | |
220 | else { | |
221 | STACK.txpol [STACK.lstack] = -twotwo; | |
222 | STACK.typol [STACK.lstack] = +zerzer; | |
223 | STACK.tzpol [STACK.lstack] = +zerzer; | |
224 | } | |
225 | ||
226 | /* Particle coordinates*/ | |
227 | STACK.xa [STACK.lstack] = particle->Vx(); | |
228 | STACK.ya [STACK.lstack] = particle->Vy(); | |
229 | STACK.za [STACK.lstack] = za; | |
230 | ||
231 | // printf("Particle Vertex %10.3f %10.3f %10.3f %10.3f \n", | |
232 | // STACK.xa [STACK.lstack], STACK.ya [STACK.lstack], STACK.za [STACK.lstack], dir); | |
233 | ||
234 | ||
235 | ||
236 | /* Calculate the total kinetic energy of the primaries: don't change*/ | |
237 | Int_t st_ilo = STACK.ilo[STACK.lstack]; | |
238 | if ( st_ilo != 0 ) | |
239 | EPISOR.tkesum += | |
240 | ((STACK.tke[STACK.lstack] + PAPROP.amdisc[st_ilo+6]) | |
241 | * STACK.wt[STACK.lstack]); | |
242 | else | |
243 | EPISOR.tkesum += (STACK.tke[STACK.lstack] * STACK.wt[STACK.lstack]); | |
244 | ||
245 | /* Here we ask for the region number of the hitting point. | |
246 | * NREG (LSTACK) = ... | |
247 | * The following line makes the starting region search much more | |
248 | * robust if particles are starting very close to a boundary: | |
249 | */ | |
250 | geocrs( STACK.tx[STACK.lstack], | |
251 | STACK.ty[STACK.lstack], | |
252 | STACK.tz[STACK.lstack] ); | |
253 | ||
254 | Int_t idisc; | |
255 | ||
256 | georeg ( STACK.xa[STACK.lstack], | |
257 | STACK.ya[STACK.lstack], | |
258 | STACK.za[STACK.lstack], | |
259 | STACK.nreg[STACK.lstack], | |
260 | idisc);//<-- dummy return variable not used | |
261 | /* Do not change these cards:*/ | |
262 | Int_t igeohsm1 = 1; | |
263 | Int_t igeohsm2 = -11; | |
264 | geohsm ( STACK.nhspnt[STACK.lstack], igeohsm1, igeohsm2, LTCLCM.mlattc ); | |
265 | STACK.nlattc[STACK.lstack] = LTCLCM.mlattc; | |
266 | soevsv(); | |
267 | } | |
268 | } | |
269 | ||
270 |