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