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