Assymmetry due TDI taken into account.
[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
16/*
17$Log$
c416706e 18Revision 1.7 2002/02/21 16:09:45 morsch
19Move SetHighwaterMark() after last possible SetTrack()
20
5bfdba97 21Revision 1.6 2002/02/01 15:28:22 morsch
22Fastidious print statements removed.
23
69818e67 24Revision 1.5 2002/02/01 14:12:28 morsch
25Include new gas pressure estimates (LHC Pproject Note 273)
26
eac2c8a3 27Revision 1.4 2001/12/19 16:30:24 morsch
28Some bugs corrected and skip method added. (Rachid Guernane)
29
05df024b 30Revision 1.3 2001/07/27 17:09:36 morsch
31Use local SetTrack, KeepTrack and SetHighWaterMark methods
32to delegate either to local stack or to stack owned by AliRun.
33(Piotr Skowronski, A.M.)
34
a99cf51f 35Revision 1.2 2001/06/14 12:15:27 morsch
36Bugs corrected. SetSide() method added.
37
154d524b 38Revision 1.1 2001/01/23 15:04:33 morsch
39Generator to read beam halo file from Protvino group.
40
14ab5373 41*/
42
43// Read background particles from a boundary source
44// Very specialized generator to simulate background from beam halo.
45// The input file is a text file specially prepared
46// for this purpose.
47// Author: andreas.morsch@cern.ch
48
49#include "AliGenHaloProtvino.h"
50#include "AliRun.h"
51#include "AliMC.h"
52#include "AliPDG.h"
53
54#include <TDatabasePDG.h>
eac2c8a3 55#include <TSystem.h>
14ab5373 56#include <stdlib.h>
57
58 ClassImp(AliGenHaloProtvino)
59 AliGenHaloProtvino::AliGenHaloProtvino()
60 :AliGenerator(-1)
61{
62// Constructor
154d524b 63
64 fName = "HaloProtvino";
65 fTitle = "Halo from LHC Tunnel";
14ab5373 66//
67// Read all particles
154d524b 68 fNpart = -1;
69 fFile = 0;
70 fSide = 1;
eac2c8a3 71//
72 SetRunPeriod();
14ab5373 73}
74
75AliGenHaloProtvino::AliGenHaloProtvino(Int_t npart)
76 :AliGenerator(npart)
77{
78// Constructor
154d524b 79 fName = "Halo";
80 fTitle= "Halo from LHC Tunnel";
14ab5373 81//
154d524b 82 fNpart = npart;
83 fFile = 0;
84 fSide = 1;
eac2c8a3 85//
86 SetRunPeriod();
14ab5373 87}
88
89AliGenHaloProtvino::AliGenHaloProtvino(const AliGenHaloProtvino & HaloProtvino)
90{
91// copy constructor
92}
93
94
95//____________________________________________________________
96AliGenHaloProtvino::~AliGenHaloProtvino()
97{
98// Destructor
99}
100
101//____________________________________________________________
102void AliGenHaloProtvino::Init()
103{
104// Initialisation
154d524b 105 fFile = fopen(fFileName,"r");
106 if (fFile) {
107 printf("\n File %s opened for reading, %p ! \n ", fFileName.Data(), fFile);
108 } else {
109 printf("\n Opening of file %s failed, %p ! \n ", fFileName.Data(), fFile);
110 }
eac2c8a3 111//
112//
113//
114// Read file with gas pressure values
115 char *name;
116
117 name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
118 FILE* file = fopen(name, "r");
119 Float_t z;
120 Int_t i, j;
121
122//
123// Ring 1
124//
125 for (i = 0; i < 20; i++)
126 {
127 fscanf(file, "%f %f %f %f %f %f", &z,
128 &fG1[i][0], &fG1[i][1], &fG1[i][2], &fG1[i][3], &fG1[i][4]);
129 if (i > 0) {
130 fZ1[i] = fZ1[i-1] + z;
131 } else {
132 fZ1[i] = 20.;
133 }
134 }
135//
136// Ring 2
137//
138 for (i = 0; i < 21; i++)
139 {
140 fscanf(file, "%f %f %f %f %f %f", &z,
141 &fG2[i][0], &fG2[i][1], &fG2[i][2], &fG2[i][3], &fG2[i][4]);
142 if (i > 0) {
143 fZ2[i] = fZ2[i-1] + z;
144 } else {
145 fZ2[i] = 20.;
146 }
147 }
148//
149// Transform into interaction rates
150//
151 const Float_t crossSection = 0.094e-28; // m^2
152 const Float_t pFlux = 1.e11/25.e-9; // 1/s
153
154 for (j = 0; j < 5; j++) {
155 for (i = 0; i < 21; i++)
156 {
157 if (i <20)
158 fG1[i][j] = fG1[i][j] * crossSection * pFlux; // 1/m/s
159 fG2[i][j] = fG2[i][j] * crossSection * pFlux; // 1/m/s
160 }
161 }
162
163
c416706e 164 Float_t sum = 0.;
eac2c8a3 165
166 for (Int_t i = 0; i < 250; i++) {
167 Float_t z = 20.+i*1.;
168 z*=100;
169 Float_t wgt = GassPressureWeight(z);
eac2c8a3 170 sum+=wgt;
171 }
172 sum/=250.;
173 printf("\n %f \n \n", sum);
14ab5373 174}
175
176//____________________________________________________________
177void AliGenHaloProtvino::Generate()
178{
179// Generate from input file
14ab5373 180
181 Float_t polar[3]= {0,0,0};
182 Float_t origin[3];
183 Float_t p[3], p0;
184 Float_t ekin, wgt, tx, ty, tz, txy;
185 Float_t zPrimary;
186 Float_t amass;
187 Int_t inuc;
188 //
189 Int_t ipart, ncols, nt;
05df024b 190 Int_t nskip = 0;
154d524b 191 Int_t nread = 0;
14ab5373 192 while(1) {
154d524b 193 ncols = fscanf(fFile,"%f %d %d %f %f %f %f %f %f",
14ab5373 194 &zPrimary, &inuc, &ipart, &wgt,
195 &ekin, &origin[0], &origin[1],
196 &tx, &ty);
154d524b 197/*
14ab5373 198 printf(" \n %f %d %d %f %f %f %f %f %f",
05df024b 199 zPrimary, inuc, ipart, wgt,
200 ekin, origin[0], origin[1],
201 tx, ty);
154d524b 202*/
05df024b 203
14ab5373 204 if (ncols < 0) break;
154d524b 205
05df024b 206 nskip++;
207 if (fNpart !=-1 && nskip <= fNskip) continue;
208
209 nread++;
14ab5373 210 if (fNpart !=-1 && nread > fNpart) break;
211
14ab5373 212 amass = TDatabasePDG::Instance()->GetParticle(ipart)->Mass();
213
214 //
215 // Momentum vector
216 //
05df024b 217 p0=sqrt(ekin*ekin + 2.*amass*ekin);
14ab5373 218
219 txy=TMath::Sqrt(tx*tx+ty*ty);
220 if (txy == 1.) {
221 tz=0;
222 } else {
223 tz=-TMath::Sqrt(1.-txy);
224 }
225
226 p[0]=p0*tx;
227 p[1]=p0*ty;
154d524b 228 p[2]=-p0*tz;
05df024b 229
230 origin[2] = -2196.5;
14ab5373 231
232 //
233 //
234 // Particle weight
235
154d524b 236 Float_t originP[3] = {0., 0., 0.};
237 originP[2] = zPrimary;
238
239 Float_t pP[3] = {0., 0., 0.};
240 Int_t ntP;
241
242 if (fSide == -1) {
243 originP[2] = -zPrimary;
244 origin[2] = -origin[2];
245 p[2] = -p[2];
246 }
05df024b 247
a99cf51f 248 SetTrack(0,-1,kProton,pP,originP,polar,0,kPNoProcess,ntP);
249 KeepTrack(ntP);
14ab5373 250 fParentWeight=wgt*GassPressureWeight(zPrimary);
a99cf51f 251 SetTrack(fTrackIt,ntP,ipart,p,origin,polar,0,kPNoProcess,nt,fParentWeight);
14ab5373 252 //
253 // Assume particles come from two directions with same probability
154d524b 254
05df024b 255 // Both Side are considered
256 if (fSide > 1) {
257 origin[2]=-origin[2];
258 p[2]=-p[2];
259 fParentWeight=wgt*GassPressureWeight(-zPrimary);
260 SetTrack(fTrackIt,ntP,ipart,p,origin,polar,0,kPNoProcess,nt,fParentWeight);
261 }
5bfdba97 262 SetHighWaterMark(nt);
14ab5373 263 }
264}
265
266
267AliGenHaloProtvino& AliGenHaloProtvino::operator=(const AliGenHaloProtvino& rhs)
268{
269// Assignment operator
270 return *this;
271}
272
273
eac2c8a3 274
14ab5373 275Float_t AliGenHaloProtvino::GassPressureWeight(Float_t zPrimary)
276{
eac2c8a3 277//
278// Return z-dependent gasspressure weight = interaction rate [1/m/s].
279//
280 Float_t weight = 0.;
281 zPrimary/=100.; // m
282 Float_t zAbs = TMath::Abs(zPrimary);
283 zPrimary = TMath::Abs(zPrimary);
154d524b 284
eac2c8a3 285 if (zPrimary < 0.)
286 {
287 if (zAbs > fZ1[19]) {
288 weight = 2.e4;
289 } else {
290 for (Int_t i = 1; i < 20; i++) {
291 if (zAbs < fZ1[i]) {
292 weight = fG1[i][fRunPeriod];
293 break;
294 }
295 }
296 }
297 } else {
298 if (zAbs > fZ2[20]) {
299 weight = 2.e4;
300 } else {
301 for (Int_t i = 1; i < 21; i++) {
302 if (zAbs < fZ2[i]) {
303 weight = fG2[i][fRunPeriod];
304 break;
305 }
306 }
307 }
308 }
154d524b 309
eac2c8a3 310 return weight;
14ab5373 311}
312
313/*
314# Title: README file for the sources of IR8 machine induced background
315# Author: Vadim Talanov <Vadim.Talanov@cern.ch>
316# Modified: 12-12-2000
317
3180. Overview
319
320 There are three files, named ring.one.beta.[01,10,50].m, which
321 contain the lists of background particles, induced by proton losses
322 upstream of IP8 in the LHC ring one, for the beta* values of 1, 10
323 and 50 m, respectively.
324
3251. File contents
326
327 Each line in the files contains the coordinates of particle track
328 crossing with the infinite plane, positioned at z=-1m, together with
329 the physical properties of corresponding particle, namely:
330
331 S - S coordinate of the primary interaction vertex, cm;
332 N - type of the gas nuclei at interaction, 1 is H, 2 - C and 3 - O;
333 I - particle ID in PDG particle numbering scheme;
334 W - particle weight;
335 E - particle kinetic energy, GeV;
336 X - x coordinate of the crossing point, cm;
337 Y - y coordinate of the crossing point, cm;
338 Dx - x direction cosine;
339 Dy - y direction cosine.
340
3412. Normalisation
342
343 Each file is given per unity of linear density of proton inelastic
344 interactions with the gas nuclei, [1 inelastic interaction/m].
345
346# ~/vtalanov/public/README.mib: the end.
347
348*/
349
350
351
352