Resolving all symbols in the library
[u/mrichter/AliRoot.git] / EVGEN / AliGenHalo.cxx
CommitLineData
4c039060 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
803d1ab0 16/* $Id$ */
4c039060 17
d12ebe74 18// Read beam halo background particles from a boundary source
19// Boundary source is in the LHCb format http://www.hep.manchester.ac.uk/u/robert/LHC_backgrounds/Note-MIBStudies.pdf
20// and has been provided by Robert Appleby
675e9664 21// Author: andreas.morsch@cern.ch
22
d12ebe74 23
116cbefd 24#include <stdlib.h>
c5970876 25
1578254f 26#include <TDatabasePDG.h>
d12ebe74 27#include <TCanvas.h>
28#include <TGraph.h>
116cbefd 29#include <TPDGCode.h>
d12ebe74 30#include <TSystem.h>
116cbefd 31
32#include "AliGenHalo.h"
2a0b4912 33#include "AliGenEventHeader.h"
116cbefd 34#include "AliRun.h"
418e091f 35#include "AliLog.h"
f87cfe57 36
198bb1c7 37ClassImp(AliGenHalo)
38
39AliGenHalo::AliGenHalo()
d12ebe74 40 :AliGenerator(-1),
41 fFile(0),
42 fFileName(0),
43 fSide(1),
44 fRunPeriod(kY3D90),
45 fTimePerEvent(1.e-4),
46 fNskip(0),
47 fZ1(0),
48 fZ2(0),
49 fG1(0),
50 fG2(0),
51 fGPASize(0),
52 fLossID(0),
53 fLossA(0),
54 fPdg(0),
55 fLossT0(0),
56 fLossZ(0),
57 fLossW(0),
58 fXS(0),
59 fYS(0),
60 fZS(0),
61 fDX(0),
62 fDY(0),
63 fEkin(0),
64 fTS(0),
65 fWS(0)
fe4da5cc 66{
f87cfe57 67// Constructor
d12ebe74 68
69 fName = "Halo";
70 fTitle = "Halo from LHC Beam";
fe4da5cc 71//
72// Read all particles
d12ebe74 73 fNpart = -1;
74 SetAnalog(0);
fe4da5cc 75}
76
77AliGenHalo::AliGenHalo(Int_t npart)
1c56e311 78 :AliGenerator(npart),
d12ebe74 79 fFile(0),
80 fFileName(0),
81 fSide(1),
82 fRunPeriod(kY3D90),
83 fTimePerEvent(1.e-4),
84 fNskip(0),
85 fZ1(0),
86 fZ2(0),
87 fG1(0),
88 fG2(0),
89 fGPASize(0),
90 fLossID(0),
91 fLossA(0),
92 fPdg(0),
93 fLossT0(0),
94 fLossZ(0),
95 fLossW(0),
96 fXS(0),
97 fYS(0),
98 fZS(0),
99 fDX(0),
100 fDY(0),
101 fEkin(0),
102 fTS(0),
103 fWS(0)
fe4da5cc 104{
f87cfe57 105// Constructor
d12ebe74 106 fName = "Halo";
107 fTitle= "Halo from LHC Beam";
fe4da5cc 108//
d12ebe74 109 fNpart = npart;
110//
111 SetAnalog(0);
fe4da5cc 112}
113
114//____________________________________________________________
115AliGenHalo::~AliGenHalo()
116{
f87cfe57 117// Destructor
fe4da5cc 118}
119
120//____________________________________________________________
121void AliGenHalo::Init()
f87cfe57 122{
123// Initialisation
d12ebe74 124 fFile = fopen(fFileName,"r");
1701cff6 125 Int_t ir = 0;
126
127
d12ebe74 128 if (fFile) {
129 printf("\n File %s opened for reading, %p ! \n ", fFileName.Data(), (void*)fFile);
130 } else {
131 printf("\n Opening of file %s failed, %p ! \n ", fFileName.Data(), (void*)fFile);
418e091f 132 return;
d12ebe74 133 }
c321951c 134
135 if (fNskip > 0) {
136 // Skip the first fNskip events
137 SkipEvents();
138 }
d12ebe74 139//
140//
141//
142// Read file with gas pressure values
143 char *name = 0;
144 if (fRunPeriod < 5) {
145 name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
146 fGPASize = 21;
147 fG1 = new Float_t[fGPASize];
148 fG2 = new Float_t[fGPASize];
149 fZ1 = new Float_t[fGPASize];
150 fZ2 = new Float_t[fGPASize];
151 } else if (fRunPeriod == 5) {
152 name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_startup.dat");
153 fGPASize = 18853;
154 fG1 = new Float_t[fGPASize];
155 fZ1 = new Float_t[fGPASize];
156 } else if (fRunPeriod ==6) {
157 name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_conditioned.dat");
158 fGPASize = 12719;
159 fG1 = new Float_t[fGPASize];
160 fZ1 = new Float_t[fGPASize];
161 } else {
162 Fatal("Init()", "No gas pressure file for given run period !");
163 }
164
418e091f 165
166 FILE* file = 0;
167 if (name) file = fopen(name, "r");
168 if (!file) {
169 AliError("No gas pressure file");
170 return;
171 }
172
d12ebe74 173 Float_t z;
174 Int_t i;
175 Float_t p[5];
176
177 const Float_t kCrossSection = 0.094e-28; // m^2
178 const Float_t kFlux = 1.e11 / 25.e-9; // protons/s
179 Float_t pFlux[5] = {0.2, 0.2, 0.3, 0.3, 1.0};
180
181 if (fRunPeriod < 5) {
182//
183// Ring 1
184//
185
186 for (i = 0; i < fGPASize; i++)
187 {
1701cff6 188 ir = fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
189 if (ir == 0) break;
190
d12ebe74 191 fG1[i] = p[fRunPeriod];
192
193 if (i > 0) {
194 fZ1[i] = fZ1[i-1] + z;
195 } else {
196 fZ1[i] = 20.;
197 }
198 }
199//
200// Ring 2
201//
202 for (i = 0; i < fGPASize; i++)
203 {
1701cff6 204 ir = fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
205 if (ir == 0) break;
206
d12ebe74 207 fG2[i] = p[fRunPeriod];
208 if (i > 0) {
209 fZ2[i] = fZ2[i-1] + z;
210 } else {
211 fZ2[i] = 20.;
212 }
213 }
214//
215// Interaction rates
216//
217 for (i = 0; i < fGPASize; i++)
218 {
219 fG1[i] = fG1[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
220 fG2[i] = fG2[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
221 }
222
223 } else {
224 for (i = 0; i < fGPASize; i++)
225 {
1701cff6 226 ir = fscanf(file, "%f %e %e %e %e %e", &z, &p[0], &p[1], &p[2], &p[3], &p[4]);
227 if (ir == 0) break;
d12ebe74 228 z /= 1000.;
229 fG1[i] = p[4] * kCrossSection * kFlux; // 1/m/s
230 // 1/3 of nominal intensity at startup
231 if (fRunPeriod == kLHCPR674Startup) fG1[i] /= 3.;
232 fZ1[i] = z;
233 }
234 }
235
236
237
238
239//
240// Transform into interaction rates
241//
242
243
244
245
246 Float_t sum1 = 0.;
247 Float_t sum2 = 0.;
248
249 for (Int_t iz = 0; iz < 300; iz++) {
250 Float_t zpos = 20. + iz * 1.;
251 zpos *= 100;
252 Float_t wgt1 = GasPressureWeight( zpos);
253 Float_t wgt2 = GasPressureWeight(-zpos);
254 sum1 += wgt1;
255 sum2 += wgt2;
256 }
257 sum1/=250.;
258 sum2/=250.;
259 printf("\n %f %f \n \n", sum1, sum2);
de4ab12f 260 delete file;
f87cfe57 261}
fe4da5cc 262
263//____________________________________________________________
264void AliGenHalo::Generate()
265{
d12ebe74 266// Generate by reading particles from input file
fe4da5cc 267
d12ebe74 268 Float_t polar[3]= {0., 0., 0.};
fe4da5cc 269 Float_t origin[3];
270 Float_t p[3], p0;
d12ebe74 271 Float_t tz, txy;
272 Float_t mass;
fe4da5cc 273 //
d12ebe74 274 Int_t nt;
275 static Bool_t first = kTRUE;
276 static Int_t oldID = -1;
277//
c321951c 278
279 if (first && (fNskip == 0)) ReadNextParticle();
280 first = kFALSE;
281 oldID = fLossID;
2a0b4912 282 Int_t np = 0;
fe4da5cc 283
284 while(1) {
d12ebe74 285 // Push particle to stack
286 mass = TDatabasePDG::Instance()->GetParticle(fPdg)->Mass();
c321951c 287 p0 = TMath::Sqrt(fEkin * fEkin + 2. * mass * fEkin);
d12ebe74 288 txy = TMath::Sqrt(fDX * fDX + fDY * fDY);
73d273a4 289 if (txy > 1.) {
290 tz = 0.;
fe4da5cc 291 } else {
d12ebe74 292 tz = - TMath::Sqrt(1. - txy);
fe4da5cc 293 }
c321951c 294
d12ebe74 295 p[0] = p0 * fDX;
296 p[1] = p0 * fDY;
297 p[2] = p0 * tz;
fe4da5cc 298
d12ebe74 299 origin[0] = fXS;
300 origin[1] = fYS;
301 origin[2] = 1950.;
fe4da5cc 302
dc3747df 303 PushTrack(fTrackIt , -1, fPdg , p, origin, polar, fTS - 1950./2.9979e10, kPNoProcess, nt, fWS);
2a0b4912 304 np++;
d12ebe74 305 Int_t nc = ReadNextParticle();
306
307 if (fLossID != oldID || nc == 0) {
308 oldID = fLossID;
309 break;
310 }
311 }
2a0b4912 312
d12ebe74 313 SetHighWaterMark(nt);
2a0b4912 314 AliGenEventHeader* header = new AliGenEventHeader("HALO");
315 header->SetNProduced(np);
316 // Passes header either to the container or to gAlice
317 if (fContainer) {
e052713d 318 header->SetName(fName);
2a0b4912 319 fContainer->AddHeader(header);
320 } else {
321 gAlice->SetGenEventHeader(header);
322 }
d12ebe74 323}
2a0b4912 324
fe4da5cc 325
d12ebe74 326Float_t AliGenHalo::GasPressureWeight(Float_t zPrimary)
327{
328//
329// Return z-dependent gasspressure weight = interaction rate [1/m/s].
330//
331 Float_t weight = 0.;
332 zPrimary /= 100.; // m
333 if (fRunPeriod < 5) {
334 Float_t zAbs = TMath::Abs(zPrimary);
335 if (zPrimary > 0.)
336 {
337 if (zAbs > fZ1[20]) {
338 weight = 2.e4;
339 } else {
340 for (Int_t i = 1; i < 21; i++) {
341 if (zAbs < fZ1[i]) {
342 weight = fG1[i];
343 break;
344 }
345 }
346 }
347 } else {
348 if (zAbs > fZ2[20]) {
349 weight = 2.e4;
350 } else {
351 for (Int_t i = 1; i < 21; i++) {
352 if (zAbs < fZ2[i]) {
353 weight = fG2[i];
354 break;
355 }
356 }
357 }
358 }
359 } else {
360 Int_t index = TMath::BinarySearch(fGPASize, fZ1, zPrimary);
361 weight = fG1[index];
362 }
363 return weight;
364}
fe4da5cc 365
d12ebe74 366void AliGenHalo::Draw(Option_t *)
367{
368// Draws the gas pressure distribution
369 Float_t z[400];
370 Float_t p[400];
371
372 for (Int_t i = 0; i < 400; i++)
373 {
374 z[i] = -20000. + Float_t(i) * 100;
375 p[i] = GasPressureWeight(z[i]);
376 }
377
378 TGraph* gr = new TGraph(400, z, p);
379 TCanvas* c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
380 c1->cd();
381 gr->Draw("AL");
382}
fe4da5cc 383
d12ebe74 384Int_t AliGenHalo::ReadNextParticle()
385{
386 // Read the next particle from the file
387 Int_t ncols = fscanf(fFile,"%d %f %f %d %f %d %f %f %f %f %f %f %f %f",
388 &fLossID, &fLossT0, &fLossZ, &fLossA, &fLossW, &fPdg, &fXS, &fYS, &fZS, &fDX, &fDY, &fEkin, &fTS, &fWS);
389 fLossZ /= 10.;
390 fXS /= 10.;
391 fYS /= 10.;
392 fZS /= 10.;
393 fTS *= 1.e-9;
394 return (ncols);
395}
c321951c 396
397void AliGenHalo::SkipEvents()
398{
399 //
400 // Skip the first fNskip events
401 Int_t skip = fNskip;
402 Int_t oldID = -1;
403
404 while (skip >= 0)
405 {
406 ReadNextParticle();
407 if (oldID != fLossID) {
408 oldID = fLossID;
409 skip--;
410 }
411 }
412}
73d273a4 413void AliGenHalo::CountEvents()
414{
4a33c50d 415 // Count total number of events
73d273a4 416 Int_t nev = 0;
417 Int_t oldID = -1;
418 Int_t nc = 1;
419 while (nc != -1)
420 {
421 nc = ReadNextParticle();
422 if (oldID != fLossID) {
423 oldID = fLossID;
424 nev++;
425 printf("Number of events %10d %10d \n", nev, oldID);
426 }
427 }
428
429
430 rewind(fFile);
431}
432
433