]>
Commit | Line | Data |
---|---|---|
14ab5373 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
88cb7938 | 16 | /* $Id$ */ |
14ab5373 | 17 | |
18 | // Read background particles from a boundary source | |
19 | // Very specialized generator to simulate background from beam halo. | |
20 | // The input file is a text file specially prepared | |
21 | // for this purpose. | |
22 | // Author: andreas.morsch@cern.ch | |
23 | ||
116cbefd | 24 | #include <stdlib.h> |
14ab5373 | 25 | |
26 | #include <TDatabasePDG.h> | |
d41a82ac | 27 | #include <TCanvas.h> |
28 | #include <TGraph.h> | |
116cbefd | 29 | #include <TPDGCode.h> |
eac2c8a3 | 30 | #include <TSystem.h> |
116cbefd | 31 | |
32 | #include "AliGenHaloProtvino.h" | |
33 | #include "AliRun.h" | |
14ab5373 | 34 | |
198bb1c7 | 35 | ClassImp(AliGenHaloProtvino) |
36 | ||
37 | AliGenHaloProtvino::AliGenHaloProtvino() | |
1c56e311 | 38 | :AliGenerator(-1), |
39 | fFile(0), | |
40 | fFileName(0), | |
41 | fSide(1), | |
42 | fRunPeriod(kY3D90), | |
43 | fTimePerEvent(1.e-4), | |
44 | fNskip(0), | |
45 | fZ1(0), | |
46 | fZ2(0), | |
47 | fG1(0), | |
48 | fG2(0), | |
49 | fGPASize(0) | |
14ab5373 | 50 | { |
51 | // Constructor | |
154d524b | 52 | |
53 | fName = "HaloProtvino"; | |
54 | fTitle = "Halo from LHC Tunnel"; | |
14ab5373 | 55 | // |
56 | // Read all particles | |
154d524b | 57 | fNpart = -1; |
7f98ef21 | 58 | SetAnalog(0); |
14ab5373 | 59 | } |
60 | ||
61 | AliGenHaloProtvino::AliGenHaloProtvino(Int_t npart) | |
1c56e311 | 62 | :AliGenerator(npart), |
63 | fFile(0), | |
64 | fFileName(0), | |
65 | fSide(1), | |
66 | fRunPeriod(kY3D90), | |
67 | fTimePerEvent(1.e-4), | |
68 | fNskip(0), | |
69 | fZ1(0), | |
70 | fZ2(0), | |
71 | fG1(0), | |
72 | fG2(0), | |
73 | fGPASize(0) | |
14ab5373 | 74 | { |
75 | // Constructor | |
154d524b | 76 | fName = "Halo"; |
77 | fTitle= "Halo from LHC Tunnel"; | |
14ab5373 | 78 | // |
154d524b | 79 | fNpart = npart; |
eac2c8a3 | 80 | // |
7f98ef21 | 81 | SetAnalog(0); |
14ab5373 | 82 | } |
83 | ||
198bb1c7 | 84 | AliGenHaloProtvino::AliGenHaloProtvino(const AliGenHaloProtvino & HaloProtvino): |
1c56e311 | 85 | AliGenerator(HaloProtvino), |
86 | fFile(0), | |
87 | fFileName(0), | |
88 | fSide(1), | |
89 | fRunPeriod(kY3D90), | |
90 | fTimePerEvent(1.e-4), | |
91 | fNskip(0), | |
92 | fZ1(0), | |
93 | fZ2(0), | |
94 | fG1(0), | |
95 | fG2(0), | |
96 | fGPASize(0) | |
14ab5373 | 97 | { |
198bb1c7 | 98 | // Copy constructor |
99 | HaloProtvino.Copy(*this); | |
14ab5373 | 100 | } |
101 | ||
102 | ||
103 | //____________________________________________________________ | |
104 | AliGenHaloProtvino::~AliGenHaloProtvino() | |
105 | { | |
106 | // Destructor | |
107 | } | |
108 | ||
109 | //____________________________________________________________ | |
110 | void AliGenHaloProtvino::Init() | |
111 | { | |
112 | // Initialisation | |
154d524b | 113 | fFile = fopen(fFileName,"r"); |
114 | if (fFile) { | |
e81f4d4f | 115 | printf("\n File %s opened for reading, %p ! \n ", fFileName.Data(), (void*)fFile); |
154d524b | 116 | } else { |
e81f4d4f | 117 | printf("\n Opening of file %s failed, %p ! \n ", fFileName.Data(), (void*)fFile); |
154d524b | 118 | } |
eac2c8a3 | 119 | // |
120 | // | |
121 | // | |
122 | // Read file with gas pressure values | |
d41a82ac | 123 | char *name = 0; |
124 | if (fRunPeriod < 5) { | |
125 | name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" ); | |
126 | fGPASize = 21; | |
127 | fG1 = new Float_t[fGPASize]; | |
128 | fG2 = new Float_t[fGPASize]; | |
129 | fZ1 = new Float_t[fGPASize]; | |
130 | fZ2 = new Float_t[fGPASize]; | |
131 | } else if (fRunPeriod == 5) { | |
132 | name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_startup.dat"); | |
133 | fGPASize = 18853; | |
134 | fG1 = new Float_t[fGPASize]; | |
135 | fZ1 = new Float_t[fGPASize]; | |
136 | } else if (fRunPeriod ==6) { | |
137 | name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_conditioned.dat"); | |
138 | fGPASize = 12719; | |
139 | fG1 = new Float_t[fGPASize]; | |
140 | fZ1 = new Float_t[fGPASize]; | |
141 | } else { | |
142 | Fatal("Init()", "No gas pressure file for given run period !"); | |
143 | } | |
144 | ||
145 | ||
eac2c8a3 | 146 | FILE* file = fopen(name, "r"); |
147 | Float_t z; | |
d41a82ac | 148 | Int_t i; |
149 | Float_t p[5]; | |
150 | ||
151 | const Float_t kCrossSection = 0.094e-28; // m^2 | |
152 | const Float_t kFlux = 1.e11 / 25.e-9; // protons/s | |
153 | Float_t pFlux[5] = {0.2, 0.2, 0.3, 0.3, 1.0}; | |
154 | ||
155 | if (fRunPeriod < 5) { | |
eac2c8a3 | 156 | // |
157 | // Ring 1 | |
158 | // | |
d41a82ac | 159 | |
160 | for (i = 0; i < fGPASize; i++) | |
161 | { | |
162 | fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]); | |
163 | fG1[i] = p[fRunPeriod]; | |
164 | ||
165 | if (i > 0) { | |
166 | fZ1[i] = fZ1[i-1] + z; | |
167 | } else { | |
168 | fZ1[i] = 20.; | |
169 | } | |
eac2c8a3 | 170 | } |
eac2c8a3 | 171 | // |
172 | // Ring 2 | |
173 | // | |
d41a82ac | 174 | for (i = 0; i < fGPASize; i++) |
175 | { | |
176 | fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]); | |
177 | fG2[i] = p[fRunPeriod]; | |
178 | if (i > 0) { | |
179 | fZ2[i] = fZ2[i-1] + z; | |
180 | } else { | |
181 | fZ2[i] = 20.; | |
182 | } | |
eac2c8a3 | 183 | } |
eac2c8a3 | 184 | // |
d41a82ac | 185 | // Interaction rates |
eac2c8a3 | 186 | // |
d41a82ac | 187 | for (i = 0; i < fGPASize; i++) |
188 | { | |
189 | fG1[i] = fG1[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s | |
190 | fG2[i] = fG2[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s | |
191 | } | |
eac2c8a3 | 192 | |
d41a82ac | 193 | } else { |
194 | for (i = 0; i < fGPASize; i++) | |
eac2c8a3 | 195 | { |
d41a82ac | 196 | fscanf(file, "%f %e %e %e %e %e", &z, &p[0], &p[1], &p[2], &p[3], &p[4]); |
197 | z /= 1000.; | |
198 | fG1[i] = p[4] * kCrossSection * kFlux; // 1/m/s | |
199 | // 1/3 of nominal intensity at startup | |
200 | if (fRunPeriod == kLHCPR674Startup) fG1[i] /= 3.; | |
201 | fZ1[i] = z; | |
eac2c8a3 | 202 | } |
203 | } | |
7f98ef21 | 204 | |
eac2c8a3 | 205 | |
d41a82ac | 206 | |
207 | ||
208 | // | |
209 | // Transform into interaction rates | |
210 | // | |
211 | ||
212 | ||
213 | ||
214 | ||
7f98ef21 | 215 | Float_t sum1 = 0.; |
216 | Float_t sum2 = 0.; | |
eac2c8a3 | 217 | |
7f98ef21 | 218 | for (Int_t i = 0; i < 300; i++) { |
eac2c8a3 | 219 | Float_t z = 20.+i*1.; |
220 | z*=100; | |
d41a82ac | 221 | Float_t wgt1 = GasPressureWeight(z); |
222 | Float_t wgt2 = GasPressureWeight(-z); | |
7f98ef21 | 223 | // printf("weight: %f %f %f %f %f \n", z, wgt1, wgt2, fZ1[20], fZ2[20]); |
224 | sum1 += wgt1; | |
225 | sum2 += wgt2; | |
eac2c8a3 | 226 | } |
7f98ef21 | 227 | sum1/=250.; |
228 | sum2/=250.; | |
229 | printf("\n %f %f \n \n", sum1, sum2); | |
14ab5373 | 230 | } |
231 | ||
232 | //____________________________________________________________ | |
233 | void AliGenHaloProtvino::Generate() | |
234 | { | |
235 | // Generate from input file | |
14ab5373 | 236 | |
237 | Float_t polar[3]= {0,0,0}; | |
238 | Float_t origin[3]; | |
239 | Float_t p[3], p0; | |
7f98ef21 | 240 | Float_t tz, txy; |
14ab5373 | 241 | Float_t amass; |
14ab5373 | 242 | // |
7f98ef21 | 243 | Int_t ncols, nt; |
3969be0c | 244 | static Int_t nskip = 0; |
154d524b | 245 | Int_t nread = 0; |
7f98ef21 | 246 | |
247 | Float_t* zPrimary = new Float_t [fNpart]; | |
248 | Int_t * inuc = new Int_t [fNpart]; | |
249 | Int_t * ipart = new Int_t [fNpart]; | |
250 | Float_t* wgt = new Float_t [fNpart]; | |
251 | Float_t* ekin = new Float_t [fNpart]; | |
252 | Float_t* vx = new Float_t [fNpart]; | |
253 | Float_t* vy = new Float_t [fNpart]; | |
254 | Float_t* tx = new Float_t [fNpart]; | |
255 | Float_t* ty = new Float_t [fNpart]; | |
256 | ||
257 | Float_t zVertexOld = -1.e10; | |
258 | Int_t nInt = 0; // Counts number of interactions | |
ac3faee4 | 259 | Float_t wwgt = 0.; |
7f98ef21 | 260 | |
14ab5373 | 261 | while(1) { |
7f98ef21 | 262 | // |
263 | // Load event into array | |
264 | // | |
154d524b | 265 | ncols = fscanf(fFile,"%f %d %d %f %f %f %f %f %f", |
7f98ef21 | 266 | &zPrimary[nread], &inuc[nread], &ipart[nread], &wgt[nread], |
267 | &ekin[nread], &vx[nread], &vy[nread], | |
268 | &tx[nread], &ty[nread]); | |
05df024b | 269 | |
14ab5373 | 270 | if (ncols < 0) break; |
7f98ef21 | 271 | // Skip fNskip events |
05df024b | 272 | nskip++; |
273 | if (fNpart !=-1 && nskip <= fNskip) continue; | |
7f98ef21 | 274 | // Count interactions |
275 | if (zPrimary[nread] != zVertexOld) { | |
276 | nInt++; | |
277 | zVertexOld = zPrimary[nread]; | |
278 | } | |
279 | // Count tracks | |
05df024b | 280 | nread++; |
14ab5373 | 281 | if (fNpart !=-1 && nread > fNpart) break; |
7f98ef21 | 282 | } |
283 | // | |
284 | // Mean time between interactions | |
285 | // | |
286 | Float_t dT = fTimePerEvent/nInt; // sec | |
287 | Float_t t = 0; // sec | |
288 | ||
289 | // | |
290 | // Loop over primaries | |
291 | // | |
292 | zVertexOld = -1.e10; | |
293 | Double_t arg = 0.; | |
294 | ||
295 | for (Int_t nprim = 0; nprim < fNpart; nprim++) | |
296 | { | |
297 | amass = TDatabasePDG::Instance()->GetParticle(ipart[nprim])->Mass(); | |
14ab5373 | 298 | |
299 | // | |
300 | // Momentum vector | |
301 | // | |
7f98ef21 | 302 | p0=sqrt(ekin[nprim]*ekin[nprim] + 2.*amass*ekin[nprim]); |
14ab5373 | 303 | |
7f98ef21 | 304 | txy=TMath::Sqrt(tx[nprim]*tx[nprim]+ty[nprim]*ty[nprim]); |
14ab5373 | 305 | if (txy == 1.) { |
306 | tz=0; | |
307 | } else { | |
308 | tz=-TMath::Sqrt(1.-txy); | |
309 | } | |
310 | ||
7f98ef21 | 311 | p[0] = p0*tx[nprim]; |
312 | p[1] = p0*ty[nprim]; | |
313 | p[2] =-p0*tz; | |
05df024b | 314 | |
7f98ef21 | 315 | origin[0] = vx[nprim]; |
316 | origin[1] = vy[nprim]; | |
05df024b | 317 | origin[2] = -2196.5; |
14ab5373 | 318 | |
319 | // | |
320 | // | |
321 | // Particle weight | |
322 | ||
154d524b | 323 | Float_t originP[3] = {0., 0., 0.}; |
7f98ef21 | 324 | originP[2] = zPrimary[nprim]; |
154d524b | 325 | |
326 | Float_t pP[3] = {0., 0., 0.}; | |
327 | Int_t ntP; | |
328 | ||
329 | if (fSide == -1) { | |
7f98ef21 | 330 | originP[2] = -zPrimary[nprim]; |
154d524b | 331 | origin[2] = -origin[2]; |
332 | p[2] = -p[2]; | |
333 | } | |
05df024b | 334 | |
14ab5373 | 335 | // |
7f98ef21 | 336 | // Time |
337 | // | |
338 | if (zPrimary[nprim] != zVertexOld) { | |
339 | while(arg==0.) arg = gRandom->Rndm(); | |
340 | t -= dT*TMath::Log(arg); // (sec) | |
341 | zVertexOld = zPrimary[nprim]; | |
342 | } | |
343 | ||
344 | // Get statistical weight according to local gas-pressure | |
345 | // | |
d41a82ac | 346 | fParentWeight=wgt[nprim]*GasPressureWeight(zPrimary[nprim]); |
7f98ef21 | 347 | |
348 | if (!fAnalog || gRandom->Rndm() < fParentWeight) { | |
349 | // Pass parent particle | |
350 | // | |
642f15cf | 351 | PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight); |
7f98ef21 | 352 | KeepTrack(ntP); |
642f15cf | 353 | PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight); |
7f98ef21 | 354 | } |
7f98ef21 | 355 | // |
356 | // Both sides are considered | |
357 | // | |
154d524b | 358 | |
05df024b | 359 | if (fSide > 1) { |
d41a82ac | 360 | fParentWeight=wgt[nprim]*GasPressureWeight(-zPrimary[nprim]); |
7f98ef21 | 361 | if (!fAnalog || gRandom->Rndm() < fParentWeight) { |
362 | origin[2] = -origin[2]; | |
363 | originP[2] = -originP[2]; | |
364 | p[2]=-p[2]; | |
642f15cf | 365 | PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight); |
7f98ef21 | 366 | KeepTrack(ntP); |
642f15cf | 367 | PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight); |
7f98ef21 | 368 | } |
05df024b | 369 | } |
ac3faee4 | 370 | wwgt += fParentWeight; |
d41a82ac | 371 | |
5bfdba97 | 372 | SetHighWaterMark(nt); |
14ab5373 | 373 | } |
7f98ef21 | 374 | delete [] zPrimary; |
375 | delete [] inuc; | |
376 | delete [] ipart; | |
377 | delete [] wgt; | |
378 | delete [] ekin; | |
379 | delete [] vx; | |
380 | delete [] vy; | |
381 | delete [] tx; | |
d41a82ac | 382 | delete [] ty; |
ac3faee4 | 383 | printf("Total weight %f\n\n", wwgt); |
d41a82ac | 384 | |
14ab5373 | 385 | } |
386 | ||
387 | ||
388 | AliGenHaloProtvino& AliGenHaloProtvino::operator=(const AliGenHaloProtvino& rhs) | |
389 | { | |
390 | // Assignment operator | |
198bb1c7 | 391 | rhs.Copy(*this); |
14ab5373 | 392 | return *this; |
393 | } | |
394 | ||
395 | ||
eac2c8a3 | 396 | |
d41a82ac | 397 | Float_t AliGenHaloProtvino::GasPressureWeight(Float_t zPrimary) |
14ab5373 | 398 | { |
eac2c8a3 | 399 | // |
400 | // Return z-dependent gasspressure weight = interaction rate [1/m/s]. | |
401 | // | |
402 | Float_t weight = 0.; | |
d41a82ac | 403 | zPrimary /= 100.; // m |
404 | if (fRunPeriod < 5) { | |
405 | Float_t zAbs = TMath::Abs(zPrimary); | |
406 | if (zPrimary > 0.) | |
407 | { | |
408 | if (zAbs > fZ1[20]) { | |
409 | weight = 2.e4; | |
410 | } else { | |
411 | for (Int_t i = 1; i < 21; i++) { | |
412 | if (zAbs < fZ1[i]) { | |
413 | weight = fG1[i]; | |
414 | break; | |
415 | } | |
eac2c8a3 | 416 | } |
417 | } | |
eac2c8a3 | 418 | } else { |
d41a82ac | 419 | if (zAbs > fZ2[20]) { |
420 | weight = 2.e4; | |
421 | } else { | |
422 | for (Int_t i = 1; i < 21; i++) { | |
423 | if (zAbs < fZ2[i]) { | |
424 | weight = fG2[i]; | |
eac2c8a3 | 425 | break; |
d41a82ac | 426 | } |
eac2c8a3 | 427 | } |
428 | } | |
429 | } | |
d41a82ac | 430 | } else { |
431 | Int_t index = TMath::BinarySearch(fGPASize, fZ1, zPrimary); | |
432 | weight = fG1[index]; | |
eac2c8a3 | 433 | } |
eac2c8a3 | 434 | return weight; |
14ab5373 | 435 | } |
436 | ||
d64dff25 | 437 | void AliGenHaloProtvino::Draw(Option_t *) |
d41a82ac | 438 | { |
439 | // Draws the gas pressure distribution | |
440 | Float_t z[400]; | |
441 | Float_t p[400]; | |
442 | ||
443 | for (Int_t i = 0; i < 400; i++) | |
444 | { | |
445 | z[i] = -20000. + Float_t(i) * 100; | |
446 | p[i] = GasPressureWeight(z[i]); | |
447 | } | |
448 | ||
449 | TGraph* gr = new TGraph(400, z, p); | |
450 | TCanvas* c1 = new TCanvas("c1","Canvas 1",400,10,600,700); | |
451 | c1->cd(); | |
452 | gr->Draw("AL"); | |
453 | } | |
454 | ||
455 | ||
dc1d768c | 456 | void AliGenHaloProtvino::Copy(TObject&) const |
198bb1c7 | 457 | { |
458 | // | |
459 | // Copy | |
460 | // | |
461 | Fatal("Copy","Not implemented!\n"); | |
462 | } | |
463 | ||
464 | ||
14ab5373 | 465 | /* |
466 | # Title: README file for the sources of IR8 machine induced background | |
467 | # Author: Vadim Talanov <Vadim.Talanov@cern.ch> | |
468 | # Modified: 12-12-2000 | |
469 | ||
470 | 0. Overview | |
471 | ||
472 | There are three files, named ring.one.beta.[01,10,50].m, which | |
473 | contain the lists of background particles, induced by proton losses | |
474 | upstream of IP8 in the LHC ring one, for the beta* values of 1, 10 | |
475 | and 50 m, respectively. | |
476 | ||
477 | 1. File contents | |
478 | ||
479 | Each line in the files contains the coordinates of particle track | |
480 | crossing with the infinite plane, positioned at z=-1m, together with | |
481 | the physical properties of corresponding particle, namely: | |
482 | ||
483 | S - S coordinate of the primary interaction vertex, cm; | |
484 | N - type of the gas nuclei at interaction, 1 is H, 2 - C and 3 - O; | |
485 | I - particle ID in PDG particle numbering scheme; | |
486 | W - particle weight; | |
487 | E - particle kinetic energy, GeV; | |
488 | X - x coordinate of the crossing point, cm; | |
489 | Y - y coordinate of the crossing point, cm; | |
490 | Dx - x direction cosine; | |
491 | Dy - y direction cosine. | |
492 | ||
493 | 2. Normalisation | |
494 | ||
495 | Each file is given per unity of linear density of proton inelastic | |
496 | interactions with the gas nuclei, [1 inelastic interaction/m]. | |
497 | ||
498 | # ~/vtalanov/public/README.mib: the end. | |
499 | ||
500 | */ | |
501 | ||
502 | ||
503 | ||
504 |