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