Introducing Header instead of Log
[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$
acc86a24 18Revision 1.11 2002/10/14 14:55:35 hristov
19Merging the VirtualMC branch to the main development branch (HEAD)
20
b9d0a01d 21Revision 1.8.4.1 2002/07/24 08:56:28 alibrary
22Updating EVGEN on TVirtulaMC
23
24Revision 1.10 2002/05/28 14:28:12 morsch
25Correct LHC nominal current profile (1st year 20%, 2nd year 30%, 3rd year 100 %)
26
718e5512 27Revision 1.9 2002/05/28 13:49:17 morsch
28- Udates for new pressure table
29- calculate time
30- first provisions for real events.
31
7f98ef21 32Revision 1.8 2002/03/22 13:00:25 morsch
33Initialize sum to 0. (Jiri Chudoba).
34
c416706e 35Revision 1.7 2002/02/21 16:09:45 morsch
36Move SetHighwaterMark() after last possible SetTrack()
37
5bfdba97 38Revision 1.6 2002/02/01 15:28:22 morsch
39Fastidious print statements removed.
40
69818e67 41Revision 1.5 2002/02/01 14:12:28 morsch
42Include new gas pressure estimates (LHC Pproject Note 273)
43
eac2c8a3 44Revision 1.4 2001/12/19 16:30:24 morsch
45Some bugs corrected and skip method added. (Rachid Guernane)
46
05df024b 47Revision 1.3 2001/07/27 17:09:36 morsch
48Use local SetTrack, KeepTrack and SetHighWaterMark methods
49to delegate either to local stack or to stack owned by AliRun.
50(Piotr Skowronski, A.M.)
51
a99cf51f 52Revision 1.2 2001/06/14 12:15:27 morsch
53Bugs corrected. SetSide() method added.
54
154d524b 55Revision 1.1 2001/01/23 15:04:33 morsch
56Generator to read beam halo file from Protvino group.
57
14ab5373 58*/
59
60// Read background particles from a boundary source
61// Very specialized generator to simulate background from beam halo.
62// The input file is a text file specially prepared
63// for this purpose.
64// Author: andreas.morsch@cern.ch
65
66#include "AliGenHaloProtvino.h"
67#include "AliRun.h"
14ab5373 68#include "AliPDG.h"
69
70#include <TDatabasePDG.h>
eac2c8a3 71#include <TSystem.h>
14ab5373 72#include <stdlib.h>
73
74 ClassImp(AliGenHaloProtvino)
75 AliGenHaloProtvino::AliGenHaloProtvino()
76 :AliGenerator(-1)
77{
78// Constructor
154d524b 79
80 fName = "HaloProtvino";
81 fTitle = "Halo from LHC Tunnel";
14ab5373 82//
83// Read all particles
154d524b 84 fNpart = -1;
85 fFile = 0;
86 fSide = 1;
eac2c8a3 87//
88 SetRunPeriod();
7f98ef21 89 SetTimePerEvent();
90 SetAnalog(0);
14ab5373 91}
92
93AliGenHaloProtvino::AliGenHaloProtvino(Int_t npart)
94 :AliGenerator(npart)
95{
96// Constructor
154d524b 97 fName = "Halo";
98 fTitle= "Halo from LHC Tunnel";
14ab5373 99//
154d524b 100 fNpart = npart;
101 fFile = 0;
102 fSide = 1;
eac2c8a3 103//
104 SetRunPeriod();
7f98ef21 105 SetTimePerEvent();
106 SetAnalog(0);
14ab5373 107}
108
109AliGenHaloProtvino::AliGenHaloProtvino(const AliGenHaloProtvino & HaloProtvino)
110{
111// copy constructor
112}
113
114
115//____________________________________________________________
116AliGenHaloProtvino::~AliGenHaloProtvino()
117{
118// Destructor
119}
120
121//____________________________________________________________
122void AliGenHaloProtvino::Init()
123{
124// Initialisation
154d524b 125 fFile = fopen(fFileName,"r");
126 if (fFile) {
127 printf("\n File %s opened for reading, %p ! \n ", fFileName.Data(), fFile);
128 } else {
129 printf("\n Opening of file %s failed, %p ! \n ", fFileName.Data(), fFile);
130 }
eac2c8a3 131//
132//
133//
134// Read file with gas pressure values
135 char *name;
136
137 name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
138 FILE* file = fopen(name, "r");
139 Float_t z;
140 Int_t i, j;
141
142//
143// Ring 1
144//
7f98ef21 145 for (i = 0; i < 21; i++)
eac2c8a3 146 {
147 fscanf(file, "%f %f %f %f %f %f", &z,
148 &fG1[i][0], &fG1[i][1], &fG1[i][2], &fG1[i][3], &fG1[i][4]);
149 if (i > 0) {
150 fZ1[i] = fZ1[i-1] + z;
151 } else {
152 fZ1[i] = 20.;
153 }
154 }
155//
156// Ring 2
157//
158 for (i = 0; i < 21; i++)
159 {
160 fscanf(file, "%f %f %f %f %f %f", &z,
161 &fG2[i][0], &fG2[i][1], &fG2[i][2], &fG2[i][3], &fG2[i][4]);
162 if (i > 0) {
163 fZ2[i] = fZ2[i-1] + z;
164 } else {
165 fZ2[i] = 20.;
166 }
167 }
168//
169// Transform into interaction rates
170//
171 const Float_t crossSection = 0.094e-28; // m^2
718e5512 172 Float_t pFlux[5] = {0.2, 0.2, 0.3, 0.3, 1.0};
eac2c8a3 173
174 for (j = 0; j < 5; j++) {
718e5512 175 pFlux[j] *= 1.e11/25.e-9;
eac2c8a3 176 for (i = 0; i < 21; i++)
177 {
718e5512 178 fG1[i][j] = fG1[i][j] * crossSection * pFlux[j]; // 1/m/s
179 fG2[i][j] = fG2[i][j] * crossSection * pFlux[j]; // 1/m/s
eac2c8a3 180 }
181 }
7f98ef21 182
eac2c8a3 183
7f98ef21 184 Float_t sum1 = 0.;
185 Float_t sum2 = 0.;
eac2c8a3 186
7f98ef21 187 for (Int_t i = 0; i < 300; i++) {
eac2c8a3 188 Float_t z = 20.+i*1.;
189 z*=100;
7f98ef21 190 Float_t wgt1 = GassPressureWeight(z);
191 Float_t wgt2 = GassPressureWeight(-z);
192// printf("weight: %f %f %f %f %f \n", z, wgt1, wgt2, fZ1[20], fZ2[20]);
193 sum1 += wgt1;
194 sum2 += wgt2;
eac2c8a3 195 }
7f98ef21 196 sum1/=250.;
197 sum2/=250.;
198 printf("\n %f %f \n \n", sum1, sum2);
14ab5373 199}
200
201//____________________________________________________________
202void AliGenHaloProtvino::Generate()
203{
204// Generate from input file
14ab5373 205
206 Float_t polar[3]= {0,0,0};
207 Float_t origin[3];
208 Float_t p[3], p0;
7f98ef21 209 Float_t tz, txy;
14ab5373 210 Float_t amass;
14ab5373 211 //
7f98ef21 212 Int_t ncols, nt;
05df024b 213 Int_t nskip = 0;
154d524b 214 Int_t nread = 0;
7f98ef21 215
216 Float_t* zPrimary = new Float_t [fNpart];
217 Int_t * inuc = new Int_t [fNpart];
218 Int_t * ipart = new Int_t [fNpart];
219 Float_t* wgt = new Float_t [fNpart];
220 Float_t* ekin = new Float_t [fNpart];
221 Float_t* vx = new Float_t [fNpart];
222 Float_t* vy = new Float_t [fNpart];
223 Float_t* tx = new Float_t [fNpart];
224 Float_t* ty = new Float_t [fNpart];
225
226 Float_t zVertexOld = -1.e10;
227 Int_t nInt = 0; // Counts number of interactions
228
14ab5373 229 while(1) {
7f98ef21 230//
231// Load event into array
232//
154d524b 233 ncols = fscanf(fFile,"%f %d %d %f %f %f %f %f %f",
7f98ef21 234 &zPrimary[nread], &inuc[nread], &ipart[nread], &wgt[nread],
235 &ekin[nread], &vx[nread], &vy[nread],
236 &tx[nread], &ty[nread]);
05df024b 237
14ab5373 238 if (ncols < 0) break;
7f98ef21 239// Skip fNskip events
05df024b 240 nskip++;
241 if (fNpart !=-1 && nskip <= fNskip) continue;
7f98ef21 242// Count interactions
243 if (zPrimary[nread] != zVertexOld) {
244 nInt++;
245 zVertexOld = zPrimary[nread];
246 }
247// Count tracks
05df024b 248 nread++;
14ab5373 249 if (fNpart !=-1 && nread > fNpart) break;
7f98ef21 250 }
251//
252// Mean time between interactions
253//
254 Float_t dT = fTimePerEvent/nInt; // sec
255 Float_t t = 0; // sec
256
257//
258// Loop over primaries
259//
260 zVertexOld = -1.e10;
261 Double_t arg = 0.;
262
263 for (Int_t nprim = 0; nprim < fNpart; nprim++)
264 {
265 amass = TDatabasePDG::Instance()->GetParticle(ipart[nprim])->Mass();
14ab5373 266
267 //
268 // Momentum vector
269 //
7f98ef21 270 p0=sqrt(ekin[nprim]*ekin[nprim] + 2.*amass*ekin[nprim]);
14ab5373 271
7f98ef21 272 txy=TMath::Sqrt(tx[nprim]*tx[nprim]+ty[nprim]*ty[nprim]);
14ab5373 273 if (txy == 1.) {
274 tz=0;
275 } else {
276 tz=-TMath::Sqrt(1.-txy);
277 }
278
7f98ef21 279 p[0] = p0*tx[nprim];
280 p[1] = p0*ty[nprim];
281 p[2] =-p0*tz;
05df024b 282
7f98ef21 283 origin[0] = vx[nprim];
284 origin[1] = vy[nprim];
05df024b 285 origin[2] = -2196.5;
14ab5373 286
287 //
288 //
289 // Particle weight
290
154d524b 291 Float_t originP[3] = {0., 0., 0.};
7f98ef21 292 originP[2] = zPrimary[nprim];
154d524b 293
294 Float_t pP[3] = {0., 0., 0.};
295 Int_t ntP;
296
297 if (fSide == -1) {
7f98ef21 298 originP[2] = -zPrimary[nprim];
154d524b 299 origin[2] = -origin[2];
300 p[2] = -p[2];
301 }
05df024b 302
14ab5373 303 //
7f98ef21 304 // Time
305 //
306 if (zPrimary[nprim] != zVertexOld) {
307 while(arg==0.) arg = gRandom->Rndm();
308 t -= dT*TMath::Log(arg); // (sec)
309 zVertexOld = zPrimary[nprim];
310 }
311
312// Get statistical weight according to local gas-pressure
313//
314 fParentWeight=wgt[nprim]*GassPressureWeight(zPrimary[nprim]);
315
316 if (!fAnalog || gRandom->Rndm() < fParentWeight) {
317// Pass parent particle
318//
319 SetTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
320 KeepTrack(ntP);
321 SetTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
322 }
323
324 //
325 // Both sides are considered
326 //
154d524b 327
05df024b 328 if (fSide > 1) {
7f98ef21 329 fParentWeight=wgt[nprim]*GassPressureWeight(-zPrimary[nprim]);
330 if (!fAnalog || gRandom->Rndm() < fParentWeight) {
331 origin[2] = -origin[2];
332 originP[2] = -originP[2];
333 p[2]=-p[2];
334 SetTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
335 KeepTrack(ntP);
336 SetTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
337 }
05df024b 338 }
5bfdba97 339 SetHighWaterMark(nt);
14ab5373 340 }
7f98ef21 341 delete [] zPrimary;
342 delete [] inuc;
343 delete [] ipart;
344 delete [] wgt;
345 delete [] ekin;
346 delete [] vx;
347 delete [] vy;
348 delete [] tx;
349 delete [] ty;
14ab5373 350}
351
352
353AliGenHaloProtvino& AliGenHaloProtvino::operator=(const AliGenHaloProtvino& rhs)
354{
355// Assignment operator
356 return *this;
357}
358
359
eac2c8a3 360
14ab5373 361Float_t AliGenHaloProtvino::GassPressureWeight(Float_t zPrimary)
362{
eac2c8a3 363//
364// Return z-dependent gasspressure weight = interaction rate [1/m/s].
365//
366 Float_t weight = 0.;
367 zPrimary/=100.; // m
368 Float_t zAbs = TMath::Abs(zPrimary);
7f98ef21 369 if (zPrimary > 0.)
eac2c8a3 370 {
7f98ef21 371 if (zAbs > fZ1[20]) {
eac2c8a3 372 weight = 2.e4;
373 } else {
7f98ef21 374 for (Int_t i = 1; i < 21; i++) {
eac2c8a3 375 if (zAbs < fZ1[i]) {
376 weight = fG1[i][fRunPeriod];
377 break;
378 }
379 }
380 }
381 } else {
382 if (zAbs > fZ2[20]) {
383 weight = 2.e4;
384 } else {
385 for (Int_t i = 1; i < 21; i++) {
386 if (zAbs < fZ2[i]) {
387 weight = fG2[i][fRunPeriod];
388 break;
389 }
390 }
391 }
392 }
eac2c8a3 393 return weight;
14ab5373 394}
395
396/*
397# Title: README file for the sources of IR8 machine induced background
398# Author: Vadim Talanov <Vadim.Talanov@cern.ch>
399# Modified: 12-12-2000
400
4010. Overview
402
403 There are three files, named ring.one.beta.[01,10,50].m, which
404 contain the lists of background particles, induced by proton losses
405 upstream of IP8 in the LHC ring one, for the beta* values of 1, 10
406 and 50 m, respectively.
407
4081. File contents
409
410 Each line in the files contains the coordinates of particle track
411 crossing with the infinite plane, positioned at z=-1m, together with
412 the physical properties of corresponding particle, namely:
413
414 S - S coordinate of the primary interaction vertex, cm;
415 N - type of the gas nuclei at interaction, 1 is H, 2 - C and 3 - O;
416 I - particle ID in PDG particle numbering scheme;
417 W - particle weight;
418 E - particle kinetic energy, GeV;
419 X - x coordinate of the crossing point, cm;
420 Y - y coordinate of the crossing point, cm;
421 Dx - x direction cosine;
422 Dy - y direction cosine.
423
4242. Normalisation
425
426 Each file is given per unity of linear density of proton inelastic
427 interactions with the gas nuclei, [1 inelastic interaction/m].
428
429# ~/vtalanov/public/README.mib: the end.
430
431*/
432
433
434
435