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