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