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