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