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