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