Transition to NewIO
[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>
116cbefd 27#include <TPDGCode.h>
eac2c8a3 28#include <TSystem.h>
116cbefd 29
30#include "AliGenHaloProtvino.h"
31#include "AliRun.h"
14ab5373 32
33 ClassImp(AliGenHaloProtvino)
34 AliGenHaloProtvino::AliGenHaloProtvino()
35 :AliGenerator(-1)
36{
37// Constructor
154d524b 38
39 fName = "HaloProtvino";
40 fTitle = "Halo from LHC Tunnel";
14ab5373 41//
42// Read all particles
154d524b 43 fNpart = -1;
44 fFile = 0;
45 fSide = 1;
eac2c8a3 46//
47 SetRunPeriod();
7f98ef21 48 SetTimePerEvent();
49 SetAnalog(0);
14ab5373 50}
51
52AliGenHaloProtvino::AliGenHaloProtvino(Int_t npart)
53 :AliGenerator(npart)
54{
55// Constructor
154d524b 56 fName = "Halo";
57 fTitle= "Halo from LHC Tunnel";
14ab5373 58//
154d524b 59 fNpart = npart;
60 fFile = 0;
61 fSide = 1;
eac2c8a3 62//
63 SetRunPeriod();
7f98ef21 64 SetTimePerEvent();
65 SetAnalog(0);
14ab5373 66}
67
68AliGenHaloProtvino::AliGenHaloProtvino(const AliGenHaloProtvino & HaloProtvino)
69{
70// copy constructor
71}
72
73
74//____________________________________________________________
75AliGenHaloProtvino::~AliGenHaloProtvino()
76{
77// Destructor
78}
79
80//____________________________________________________________
81void AliGenHaloProtvino::Init()
82{
83// Initialisation
154d524b 84 fFile = fopen(fFileName,"r");
85 if (fFile) {
86 printf("\n File %s opened for reading, %p ! \n ", fFileName.Data(), fFile);
87 } else {
88 printf("\n Opening of file %s failed, %p ! \n ", fFileName.Data(), fFile);
89 }
eac2c8a3 90//
91//
92//
93// Read file with gas pressure values
94 char *name;
95
96 name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
97 FILE* file = fopen(name, "r");
98 Float_t z;
99 Int_t i, j;
100
101//
102// Ring 1
103//
7f98ef21 104 for (i = 0; i < 21; i++)
eac2c8a3 105 {
106 fscanf(file, "%f %f %f %f %f %f", &z,
107 &fG1[i][0], &fG1[i][1], &fG1[i][2], &fG1[i][3], &fG1[i][4]);
108 if (i > 0) {
109 fZ1[i] = fZ1[i-1] + z;
110 } else {
111 fZ1[i] = 20.;
112 }
113 }
114//
115// Ring 2
116//
117 for (i = 0; i < 21; i++)
118 {
119 fscanf(file, "%f %f %f %f %f %f", &z,
120 &fG2[i][0], &fG2[i][1], &fG2[i][2], &fG2[i][3], &fG2[i][4]);
121 if (i > 0) {
122 fZ2[i] = fZ2[i-1] + z;
123 } else {
124 fZ2[i] = 20.;
125 }
126 }
127//
128// Transform into interaction rates
129//
130 const Float_t crossSection = 0.094e-28; // m^2
718e5512 131 Float_t pFlux[5] = {0.2, 0.2, 0.3, 0.3, 1.0};
eac2c8a3 132
133 for (j = 0; j < 5; j++) {
718e5512 134 pFlux[j] *= 1.e11/25.e-9;
eac2c8a3 135 for (i = 0; i < 21; i++)
136 {
718e5512 137 fG1[i][j] = fG1[i][j] * crossSection * pFlux[j]; // 1/m/s
138 fG2[i][j] = fG2[i][j] * crossSection * pFlux[j]; // 1/m/s
eac2c8a3 139 }
140 }
7f98ef21 141
eac2c8a3 142
7f98ef21 143 Float_t sum1 = 0.;
144 Float_t sum2 = 0.;
eac2c8a3 145
7f98ef21 146 for (Int_t i = 0; i < 300; i++) {
eac2c8a3 147 Float_t z = 20.+i*1.;
148 z*=100;
7f98ef21 149 Float_t wgt1 = GassPressureWeight(z);
150 Float_t wgt2 = GassPressureWeight(-z);
151// printf("weight: %f %f %f %f %f \n", z, wgt1, wgt2, fZ1[20], fZ2[20]);
152 sum1 += wgt1;
153 sum2 += wgt2;
eac2c8a3 154 }
7f98ef21 155 sum1/=250.;
156 sum2/=250.;
157 printf("\n %f %f \n \n", sum1, sum2);
14ab5373 158}
159
160//____________________________________________________________
161void AliGenHaloProtvino::Generate()
162{
163// Generate from input file
14ab5373 164
165 Float_t polar[3]= {0,0,0};
166 Float_t origin[3];
167 Float_t p[3], p0;
7f98ef21 168 Float_t tz, txy;
14ab5373 169 Float_t amass;
14ab5373 170 //
7f98ef21 171 Int_t ncols, nt;
05df024b 172 Int_t nskip = 0;
154d524b 173 Int_t nread = 0;
7f98ef21 174
175 Float_t* zPrimary = new Float_t [fNpart];
176 Int_t * inuc = new Int_t [fNpart];
177 Int_t * ipart = new Int_t [fNpart];
178 Float_t* wgt = new Float_t [fNpart];
179 Float_t* ekin = new Float_t [fNpart];
180 Float_t* vx = new Float_t [fNpart];
181 Float_t* vy = new Float_t [fNpart];
182 Float_t* tx = new Float_t [fNpart];
183 Float_t* ty = new Float_t [fNpart];
184
185 Float_t zVertexOld = -1.e10;
186 Int_t nInt = 0; // Counts number of interactions
187
14ab5373 188 while(1) {
7f98ef21 189//
190// Load event into array
191//
154d524b 192 ncols = fscanf(fFile,"%f %d %d %f %f %f %f %f %f",
7f98ef21 193 &zPrimary[nread], &inuc[nread], &ipart[nread], &wgt[nread],
194 &ekin[nread], &vx[nread], &vy[nread],
195 &tx[nread], &ty[nread]);
05df024b 196
14ab5373 197 if (ncols < 0) break;
7f98ef21 198// Skip fNskip events
05df024b 199 nskip++;
200 if (fNpart !=-1 && nskip <= fNskip) continue;
7f98ef21 201// Count interactions
202 if (zPrimary[nread] != zVertexOld) {
203 nInt++;
204 zVertexOld = zPrimary[nread];
205 }
206// Count tracks
05df024b 207 nread++;
14ab5373 208 if (fNpart !=-1 && nread > fNpart) break;
7f98ef21 209 }
210//
211// Mean time between interactions
212//
213 Float_t dT = fTimePerEvent/nInt; // sec
214 Float_t t = 0; // sec
215
216//
217// Loop over primaries
218//
219 zVertexOld = -1.e10;
220 Double_t arg = 0.;
221
222 for (Int_t nprim = 0; nprim < fNpart; nprim++)
223 {
224 amass = TDatabasePDG::Instance()->GetParticle(ipart[nprim])->Mass();
14ab5373 225
226 //
227 // Momentum vector
228 //
7f98ef21 229 p0=sqrt(ekin[nprim]*ekin[nprim] + 2.*amass*ekin[nprim]);
14ab5373 230
7f98ef21 231 txy=TMath::Sqrt(tx[nprim]*tx[nprim]+ty[nprim]*ty[nprim]);
14ab5373 232 if (txy == 1.) {
233 tz=0;
234 } else {
235 tz=-TMath::Sqrt(1.-txy);
236 }
237
7f98ef21 238 p[0] = p0*tx[nprim];
239 p[1] = p0*ty[nprim];
240 p[2] =-p0*tz;
05df024b 241
7f98ef21 242 origin[0] = vx[nprim];
243 origin[1] = vy[nprim];
05df024b 244 origin[2] = -2196.5;
14ab5373 245
246 //
247 //
248 // Particle weight
249
154d524b 250 Float_t originP[3] = {0., 0., 0.};
7f98ef21 251 originP[2] = zPrimary[nprim];
154d524b 252
253 Float_t pP[3] = {0., 0., 0.};
254 Int_t ntP;
255
256 if (fSide == -1) {
7f98ef21 257 originP[2] = -zPrimary[nprim];
154d524b 258 origin[2] = -origin[2];
259 p[2] = -p[2];
260 }
05df024b 261
14ab5373 262 //
7f98ef21 263 // Time
264 //
265 if (zPrimary[nprim] != zVertexOld) {
266 while(arg==0.) arg = gRandom->Rndm();
267 t -= dT*TMath::Log(arg); // (sec)
268 zVertexOld = zPrimary[nprim];
269 }
270
271// Get statistical weight according to local gas-pressure
272//
273 fParentWeight=wgt[nprim]*GassPressureWeight(zPrimary[nprim]);
274
275 if (!fAnalog || gRandom->Rndm() < fParentWeight) {
276// Pass parent particle
277//
278 SetTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
279 KeepTrack(ntP);
280 SetTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
281 }
282
283 //
284 // Both sides are considered
285 //
154d524b 286
05df024b 287 if (fSide > 1) {
7f98ef21 288 fParentWeight=wgt[nprim]*GassPressureWeight(-zPrimary[nprim]);
289 if (!fAnalog || gRandom->Rndm() < fParentWeight) {
290 origin[2] = -origin[2];
291 originP[2] = -originP[2];
292 p[2]=-p[2];
293 SetTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
294 KeepTrack(ntP);
295 SetTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
296 }
05df024b 297 }
5bfdba97 298 SetHighWaterMark(nt);
14ab5373 299 }
7f98ef21 300 delete [] zPrimary;
301 delete [] inuc;
302 delete [] ipart;
303 delete [] wgt;
304 delete [] ekin;
305 delete [] vx;
306 delete [] vy;
307 delete [] tx;
308 delete [] ty;
14ab5373 309}
310
311
312AliGenHaloProtvino& AliGenHaloProtvino::operator=(const AliGenHaloProtvino& rhs)
313{
314// Assignment operator
315 return *this;
316}
317
318
eac2c8a3 319
14ab5373 320Float_t AliGenHaloProtvino::GassPressureWeight(Float_t zPrimary)
321{
eac2c8a3 322//
323// Return z-dependent gasspressure weight = interaction rate [1/m/s].
324//
325 Float_t weight = 0.;
326 zPrimary/=100.; // m
327 Float_t zAbs = TMath::Abs(zPrimary);
7f98ef21 328 if (zPrimary > 0.)
eac2c8a3 329 {
7f98ef21 330 if (zAbs > fZ1[20]) {
eac2c8a3 331 weight = 2.e4;
332 } else {
7f98ef21 333 for (Int_t i = 1; i < 21; i++) {
eac2c8a3 334 if (zAbs < fZ1[i]) {
335 weight = fG1[i][fRunPeriod];
336 break;
337 }
338 }
339 }
340 } else {
341 if (zAbs > fZ2[20]) {
342 weight = 2.e4;
343 } else {
344 for (Int_t i = 1; i < 21; i++) {
345 if (zAbs < fZ2[i]) {
346 weight = fG2[i][fRunPeriod];
347 break;
348 }
349 }
350 }
351 }
eac2c8a3 352 return weight;
14ab5373 353}
354
355/*
356# Title: README file for the sources of IR8 machine induced background
357# Author: Vadim Talanov <Vadim.Talanov@cern.ch>
358# Modified: 12-12-2000
359
3600. Overview
361
362 There are three files, named ring.one.beta.[01,10,50].m, which
363 contain the lists of background particles, induced by proton losses
364 upstream of IP8 in the LHC ring one, for the beta* values of 1, 10
365 and 50 m, respectively.
366
3671. File contents
368
369 Each line in the files contains the coordinates of particle track
370 crossing with the infinite plane, positioned at z=-1m, together with
371 the physical properties of corresponding particle, namely:
372
373 S - S coordinate of the primary interaction vertex, cm;
374 N - type of the gas nuclei at interaction, 1 is H, 2 - C and 3 - O;
375 I - particle ID in PDG particle numbering scheme;
376 W - particle weight;
377 E - particle kinetic energy, GeV;
378 X - x coordinate of the crossing point, cm;
379 Y - y coordinate of the crossing point, cm;
380 Dx - x direction cosine;
381 Dy - y direction cosine.
382
3832. Normalisation
384
385 Each file is given per unity of linear density of proton inelastic
386 interactions with the gas nuclei, [1 inelastic interaction/m].
387
388# ~/vtalanov/public/README.mib: the end.
389
390*/
391
392
393
394