Shadowing of variabled corrected.
[u/mrichter/AliRoot.git] / EVGEN / AliGenHaloProtvino.cxx
CommitLineData
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 35ClassImp(AliGenHaloProtvino)
36
37AliGenHaloProtvino::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
61AliGenHaloProtvino::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
2ddeccf2 84//____________________________________________________________
85AliGenHaloProtvino::~AliGenHaloProtvino()
86{
87// Destructor
88}
89
90//____________________________________________________________
91void AliGenHaloProtvino::Init()
92{
93// Initialisation
94 fFile = fopen(fFileName,"r");
95 if (fFile) {
96 printf("\n File %s opened for reading, %p ! \n ", fFileName.Data(), (void*)fFile);
97 } else {
98 printf("\n Opening of file %s failed, %p ! \n ", fFileName.Data(), (void*)fFile);
99 }
100//
101//
102//
103// Read file with gas pressure values
104 char *name = 0;
105 if (fRunPeriod < 5) {
106 name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
107 fGPASize = 21;
108 fG1 = new Float_t[fGPASize];
109 fG2 = new Float_t[fGPASize];
110 fZ1 = new Float_t[fGPASize];
111 fZ2 = new Float_t[fGPASize];
112 } else if (fRunPeriod == 5) {
113 name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_startup.dat");
114 fGPASize = 18853;
115 fG1 = new Float_t[fGPASize];
116 fZ1 = new Float_t[fGPASize];
117 } else if (fRunPeriod ==6) {
118 name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_conditioned.dat");
119 fGPASize = 12719;
120 fG1 = new Float_t[fGPASize];
121 fZ1 = new Float_t[fGPASize];
122 } else {
123 Fatal("Init()", "No gas pressure file for given run period !");
124 }
125
126
127 FILE* file = fopen(name, "r");
128 Float_t z;
129 Int_t i;
130 Float_t p[5];
131
132 const Float_t kCrossSection = 0.094e-28; // m^2
133 const Float_t kFlux = 1.e11 / 25.e-9; // protons/s
134 Float_t pFlux[5] = {0.2, 0.2, 0.3, 0.3, 1.0};
135
136 if (fRunPeriod < 5) {
137//
138// Ring 1
139//
140
141 for (i = 0; i < fGPASize; i++)
142 {
143 fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
144 fG1[i] = p[fRunPeriod];
145
146 if (i > 0) {
147 fZ1[i] = fZ1[i-1] + z;
148 } else {
149 fZ1[i] = 20.;
150 }
151 }
152//
153// Ring 2
154//
155 for (i = 0; i < fGPASize; i++)
156 {
157 fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
158 fG2[i] = p[fRunPeriod];
159 if (i > 0) {
160 fZ2[i] = fZ2[i-1] + z;
161 } else {
162 fZ2[i] = 20.;
163 }
164 }
165//
166// Interaction rates
167//
168 for (i = 0; i < fGPASize; i++)
169 {
170 fG1[i] = fG1[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
171 fG2[i] = fG2[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
172 }
173
174 } else {
175 for (i = 0; i < fGPASize; i++)
176 {
177 fscanf(file, "%f %e %e %e %e %e", &z, &p[0], &p[1], &p[2], &p[3], &p[4]);
178 z /= 1000.;
179 fG1[i] = p[4] * kCrossSection * kFlux; // 1/m/s
180 // 1/3 of nominal intensity at startup
181 if (fRunPeriod == kLHCPR674Startup) fG1[i] /= 3.;
182 fZ1[i] = z;
183 }
184 }
185
186
187
188
189//
190// Transform into interaction rates
191//
192
193
194
195
196 Float_t sum1 = 0.;
197 Float_t sum2 = 0.;
198
ea79897e 199 for (Int_t iz = 0; iz < 300; iz++) {
200 Float_t zpos = 20. + iz * 1.;
201 zpos *= 100;
202 Float_t wgt1 = GasPressureWeight( zpos);
203 Float_t wgt2 = GasPressureWeight(-zpos);
2ddeccf2 204 sum1 += wgt1;
205 sum2 += wgt2;
206 }
207 sum1/=250.;
208 sum2/=250.;
209 printf("\n %f %f \n \n", sum1, sum2);
210}
211
212//____________________________________________________________
213void AliGenHaloProtvino::Generate()
214{
215// Generate from input file
216
217 Float_t polar[3]= {0,0,0};
218 Float_t origin[3];
219 Float_t p[3], p0;
220 Float_t tz, txy;
221 Float_t amass;
222 //
223 Int_t ncols, nt;
224 static Int_t nskip = 0;
225 Int_t nread = 0;
226
227 Float_t* zPrimary = new Float_t [fNpart];
228 Int_t * inuc = new Int_t [fNpart];
229 Int_t * ipart = new Int_t [fNpart];
230 Float_t* wgt = new Float_t [fNpart];
231 Float_t* ekin = new Float_t [fNpart];
232 Float_t* vx = new Float_t [fNpart];
233 Float_t* vy = new Float_t [fNpart];
234 Float_t* tx = new Float_t [fNpart];
235 Float_t* ty = new Float_t [fNpart];
236
237 Float_t zVertexOld = -1.e10;
238 Int_t nInt = 0; // Counts number of interactions
239 Float_t wwgt = 0.;
240
241 while(1) {
242//
243// Load event into array
244//
245 ncols = fscanf(fFile,"%f %d %d %f %f %f %f %f %f",
246 &zPrimary[nread], &inuc[nread], &ipart[nread], &wgt[nread],
247 &ekin[nread], &vx[nread], &vy[nread],
248 &tx[nread], &ty[nread]);
249
250 if (ncols < 0) break;
251// Skip fNskip events
252 nskip++;
253 if (fNpart !=-1 && nskip <= fNskip) continue;
254// Count interactions
255 if (zPrimary[nread] != zVertexOld) {
256 nInt++;
257 zVertexOld = zPrimary[nread];
258 }
259// Count tracks
260 nread++;
261 if (fNpart !=-1 && nread > fNpart) break;
262 }
263//
264// Mean time between interactions
265//
266 Float_t dT = fTimePerEvent/nInt; // sec
267 Float_t t = 0; // sec
268
269//
270// Loop over primaries
271//
272 zVertexOld = -1.e10;
273 Double_t arg = 0.;
274
275 for (Int_t nprim = 0; nprim < fNpart; nprim++)
276 {
277 amass = TDatabasePDG::Instance()->GetParticle(ipart[nprim])->Mass();
278
279 //
280 // Momentum vector
281 //
282 p0=sqrt(ekin[nprim]*ekin[nprim] + 2.*amass*ekin[nprim]);
283
284 txy=TMath::Sqrt(tx[nprim]*tx[nprim]+ty[nprim]*ty[nprim]);
285 if (txy == 1.) {
286 tz=0;
287 } else {
288 tz=-TMath::Sqrt(1.-txy);
289 }
290
291 p[0] = p0*tx[nprim];
292 p[1] = p0*ty[nprim];
293 p[2] =-p0*tz;
294
295 origin[0] = vx[nprim];
296 origin[1] = vy[nprim];
297 origin[2] = -2196.5;
298
299 //
300 //
301 // Particle weight
302
303 Float_t originP[3] = {0., 0., 0.};
304 originP[2] = zPrimary[nprim];
305
306 Float_t pP[3] = {0., 0., 0.};
307 Int_t ntP;
308
309 if (fSide == -1) {
310 originP[2] = -zPrimary[nprim];
311 origin[2] = -origin[2];
312 p[2] = -p[2];
313 }
314
315 //
316 // Time
317 //
318 if (zPrimary[nprim] != zVertexOld) {
319 while(arg==0.) arg = gRandom->Rndm();
320 t -= dT*TMath::Log(arg); // (sec)
321 zVertexOld = zPrimary[nprim];
322 }
323
324// Get statistical weight according to local gas-pressure
325//
326 fParentWeight=wgt[nprim]*GasPressureWeight(zPrimary[nprim]);
327
328 if (!fAnalog || gRandom->Rndm() < fParentWeight) {
329// Pass parent particle
330//
331 PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
332 KeepTrack(ntP);
333 PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
334 }
335 //
336 // Both sides are considered
337 //
338
339 if (fSide > 1) {
340 fParentWeight=wgt[nprim]*GasPressureWeight(-zPrimary[nprim]);
341 if (!fAnalog || gRandom->Rndm() < fParentWeight) {
342 origin[2] = -origin[2];
343 originP[2] = -originP[2];
344 p[2]=-p[2];
345 PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
346 KeepTrack(ntP);
347 PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
348 }
349 }
350 wwgt += fParentWeight;
351
352 SetHighWaterMark(nt);
353 }
354 delete [] zPrimary;
355 delete [] inuc;
356 delete [] ipart;
357 delete [] wgt;
358 delete [] ekin;
359 delete [] vx;
360 delete [] vy;
361 delete [] tx;
362 delete [] ty;
363 printf("Total weight %f\n\n", wwgt);
364
365}
366
367
368Float_t AliGenHaloProtvino::GasPressureWeight(Float_t zPrimary)
369{
370//
371// Return z-dependent gasspressure weight = interaction rate [1/m/s].
372//
373 Float_t weight = 0.;
374 zPrimary /= 100.; // m
375 if (fRunPeriod < 5) {
376 Float_t zAbs = TMath::Abs(zPrimary);
377 if (zPrimary > 0.)
378 {
379 if (zAbs > fZ1[20]) {
380 weight = 2.e4;
381 } else {
382 for (Int_t i = 1; i < 21; i++) {
383 if (zAbs < fZ1[i]) {
384 weight = fG1[i];
385 break;
386 }
387 }
388 }
389 } else {
390 if (zAbs > fZ2[20]) {
391 weight = 2.e4;
392 } else {
393 for (Int_t i = 1; i < 21; i++) {
394 if (zAbs < fZ2[i]) {
395 weight = fG2[i];
396 break;
397 }
398 }
399 }
400 }
401 } else {
402 Int_t index = TMath::BinarySearch(fGPASize, fZ1, zPrimary);
403 weight = fG1[index];
404 }
405 return weight;
406}
407
408void AliGenHaloProtvino::Draw(Option_t *)
409{
410// Draws the gas pressure distribution
411 Float_t z[400];
412 Float_t p[400];
413
414 for (Int_t i = 0; i < 400; i++)
415 {
416 z[i] = -20000. + Float_t(i) * 100;
417 p[i] = GasPressureWeight(z[i]);
418 }
419
420 TGraph* gr = new TGraph(400, z, p);
421 TCanvas* c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
422 c1->cd();
423 gr->Draw("AL");
424}
425
198bb1c7 426
14ab5373 427/*
428# Title: README file for the sources of IR8 machine induced background
429# Author: Vadim Talanov <Vadim.Talanov@cern.ch>
430# Modified: 12-12-2000
431
4320. Overview
433
434 There are three files, named ring.one.beta.[01,10,50].m, which
435 contain the lists of background particles, induced by proton losses
436 upstream of IP8 in the LHC ring one, for the beta* values of 1, 10
437 and 50 m, respectively.
438
4391. File contents
440
441 Each line in the files contains the coordinates of particle track
442 crossing with the infinite plane, positioned at z=-1m, together with
443 the physical properties of corresponding particle, namely:
444
445 S - S coordinate of the primary interaction vertex, cm;
446 N - type of the gas nuclei at interaction, 1 is H, 2 - C and 3 - O;
447 I - particle ID in PDG particle numbering scheme;
448 W - particle weight;
449 E - particle kinetic energy, GeV;
450 X - x coordinate of the crossing point, cm;
451 Y - y coordinate of the crossing point, cm;
452 Dx - x direction cosine;
453 Dy - y direction cosine.
454
4552. Normalisation
456
457 Each file is given per unity of linear density of proton inelastic
458 interactions with the gas nuclei, [1 inelastic interaction/m].
459
460# ~/vtalanov/public/README.mib: the end.
461
462*/
463
464
465
466