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