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 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 :
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 "AliLog.h"
33 : #include "AliGenHaloProtvino.h"
34 : #include "AliRun.h"
35 :
36 6 : ClassImp(AliGenHaloProtvino)
37 :
38 : AliGenHaloProtvino::AliGenHaloProtvino()
39 0 : :AliGenerator(-1),
40 0 : fFile(0),
41 0 : fFileName(0),
42 0 : fSide(1),
43 0 : fRunPeriod(kY3D90),
44 0 : fTimePerEvent(1.e-4),
45 0 : fNskip(0),
46 0 : fZ1(0),
47 0 : fZ2(0),
48 0 : fG1(0),
49 0 : fG2(0),
50 0 : fGPASize(0)
51 0 : {
52 : // Constructor
53 :
54 0 : fName = "HaloProtvino";
55 0 : fTitle = "Halo from LHC Tunnel";
56 : //
57 : // Read all particles
58 0 : fNpart = -1;
59 0 : SetAnalog(0);
60 0 : }
61 :
62 : AliGenHaloProtvino::AliGenHaloProtvino(Int_t npart)
63 0 : :AliGenerator(npart),
64 0 : fFile(0),
65 0 : fFileName(0),
66 0 : fSide(1),
67 0 : fRunPeriod(kY3D90),
68 0 : fTimePerEvent(1.e-4),
69 0 : fNskip(0),
70 0 : fZ1(0),
71 0 : fZ2(0),
72 0 : fG1(0),
73 0 : fG2(0),
74 0 : fGPASize(0)
75 0 : {
76 : // Constructor
77 0 : fName = "Halo";
78 0 : fTitle= "Halo from LHC Tunnel";
79 : //
80 0 : fNpart = npart;
81 : //
82 0 : SetAnalog(0);
83 0 : }
84 :
85 : //____________________________________________________________
86 0 : AliGenHaloProtvino::~AliGenHaloProtvino()
87 0 : {
88 : // Destructor
89 0 : }
90 :
91 : //____________________________________________________________
92 : void AliGenHaloProtvino::Init()
93 : {
94 : // Initialisation
95 0 : fFile = fopen(fFileName,"r");
96 0 : if (fFile) {
97 0 : printf("\n File %s opened for reading, %p ! \n ", fFileName.Data(), (void*)fFile);
98 0 : } else {
99 0 : 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 0 : if (fRunPeriod < 5) {
107 0 : name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
108 0 : fGPASize = 21;
109 0 : fG1 = new Float_t[fGPASize];
110 0 : fG2 = new Float_t[fGPASize];
111 0 : fZ1 = new Float_t[fGPASize];
112 0 : fZ2 = new Float_t[fGPASize];
113 0 : } else if (fRunPeriod == 5) {
114 0 : name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_startup.dat");
115 0 : fGPASize = 18853;
116 0 : fG1 = new Float_t[fGPASize];
117 0 : fZ1 = new Float_t[fGPASize];
118 0 : } else if (fRunPeriod ==6) {
119 0 : name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_conditioned.dat");
120 0 : fGPASize = 12719;
121 0 : fG1 = new Float_t[fGPASize];
122 0 : fZ1 = new Float_t[fGPASize];
123 0 : } else {
124 0 : Fatal("Init()", "No gas pressure file for given run period !");
125 : }
126 :
127 : FILE* file = 0;
128 0 : if (name) file = fopen(name, "r");
129 0 : if (!file) {
130 0 : AliError("No gas pressure file");
131 0 : return;
132 : }
133 :
134 0 : Float_t z;
135 : Int_t i;
136 0 : 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 0 : Float_t pFlux[5] = {0.2, 0.2, 0.3, 0.3, 1.0};
141 :
142 : Int_t ncols = 0;
143 0 : if (fRunPeriod < 5) {
144 : //
145 : // Ring 1
146 : //
147 :
148 0 : for (i = 0; i < fGPASize; i++)
149 : {
150 0 : ncols = fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
151 0 : if (ncols<0) break;
152 :
153 0 : fG1[i] = p[fRunPeriod];
154 :
155 0 : if (i > 0) {
156 0 : fZ1[i] = fZ1[i-1] + z;
157 0 : } else {
158 0 : fZ1[i] = 20.;
159 : }
160 : }
161 : //
162 : // Ring 2
163 : //
164 0 : for (i = 0; i < fGPASize; i++)
165 : {
166 0 : ncols = fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
167 0 : if (ncols<0) break;
168 :
169 0 : fG2[i] = p[fRunPeriod];
170 0 : if (i > 0) {
171 0 : fZ2[i] = fZ2[i-1] + z;
172 0 : } else {
173 0 : fZ2[i] = 20.;
174 : }
175 : }
176 : //
177 : // Interaction rates
178 : //
179 0 : for (i = 0; i < fGPASize; i++)
180 : {
181 0 : fG1[i] = fG1[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
182 0 : fG2[i] = fG2[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
183 : }
184 :
185 : } else {
186 0 : for (i = 0; i < fGPASize; i++)
187 : {
188 0 : ncols = fscanf(file, "%f %e %e %e %e %e", &z, &p[0], &p[1], &p[2], &p[3], &p[4]);
189 0 : if (ncols<0) break;
190 :
191 0 : z /= 1000.;
192 0 : fG1[i] = p[4] * kCrossSection * kFlux; // 1/m/s
193 : // 1/3 of nominal intensity at startup
194 0 : if (fRunPeriod == kLHCPR674Startup) fG1[i] /= 3.;
195 0 : 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 :
212 0 : for (Int_t iz = 0; iz < 300; iz++) {
213 0 : Float_t zpos = 20. + iz * 1.;
214 0 : zpos *= 100;
215 0 : Float_t wgt1 = GasPressureWeight( zpos);
216 0 : Float_t wgt2 = GasPressureWeight(-zpos);
217 0 : sum1 += wgt1;
218 0 : sum2 += wgt2;
219 : }
220 0 : sum1/=250.;
221 0 : sum2/=250.;
222 0 : printf("\n %f %f \n \n", sum1, sum2);
223 0 : delete file;
224 0 : }
225 :
226 : //____________________________________________________________
227 : void AliGenHaloProtvino::Generate()
228 : {
229 : // Generate from input file
230 :
231 0 : Float_t polar[3]= {0,0,0};
232 0 : Float_t origin[3];
233 0 : Float_t p[3], p0;
234 : Float_t tz, txy;
235 : Float_t amass;
236 : //
237 0 : Int_t ncols, nt;
238 : static Int_t nskip = 0;
239 : Int_t nread = 0;
240 :
241 0 : Float_t* zPrimary = new Float_t [fNpart];
242 0 : Int_t * inuc = new Int_t [fNpart];
243 0 : Int_t * ipart = new Int_t [fNpart];
244 0 : Float_t* wgt = new Float_t [fNpart];
245 0 : Float_t* ekin = new Float_t [fNpart];
246 0 : Float_t* vx = new Float_t [fNpart];
247 0 : Float_t* vy = new Float_t [fNpart];
248 0 : Float_t* tx = new Float_t [fNpart];
249 0 : 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 0 : while(1) {
256 : //
257 : // Load event into array
258 : //
259 0 : ncols = fscanf(fFile,"%f %d %d %f %f %f %f %f %f",
260 0 : &zPrimary[nread], &inuc[nread], &ipart[nread], &wgt[nread],
261 0 : &ekin[nread], &vx[nread], &vy[nread],
262 0 : &tx[nread], &ty[nread]);
263 :
264 0 : if (ncols < 0) break;
265 : // Skip fNskip events
266 0 : nskip++;
267 0 : if (fNpart !=-1 && nskip <= fNskip) continue;
268 : // Count interactions
269 0 : if (zPrimary[nread] != zVertexOld) {
270 0 : nInt++;
271 : zVertexOld = zPrimary[nread];
272 0 : }
273 : // Count tracks
274 0 : nread++;
275 0 : if (fNpart !=-1 && nread >= fNpart) break;
276 : }
277 : //
278 : // Mean time between interactions
279 : //
280 :
281 : Float_t dT = 0.; // sec
282 0 : if (nInt > 0)
283 0 : dT = fTimePerEvent/nInt;
284 : Float_t t = 0; // sec
285 :
286 : //
287 : // Loop over primaries
288 : //
289 : zVertexOld = -1.e10;
290 : Double_t arg = 0.;
291 :
292 0 : for (Int_t nprim = 0; nprim < fNpart; nprim++)
293 : {
294 0 : amass = TDatabasePDG::Instance()->GetParticle(ipart[nprim])->Mass();
295 :
296 : //
297 : // Momentum vector
298 : //
299 0 : p0=sqrt(ekin[nprim]*ekin[nprim] + 2.*amass*ekin[nprim]);
300 :
301 0 : txy=TMath::Sqrt(tx[nprim]*tx[nprim]+ty[nprim]*ty[nprim]);
302 0 : if (txy == 1.) {
303 : tz=0;
304 0 : } else {
305 0 : tz=-TMath::Sqrt(1.-txy);
306 : }
307 :
308 0 : p[0] = p0*tx[nprim];
309 0 : p[1] = p0*ty[nprim];
310 0 : p[2] =-p0*tz;
311 :
312 0 : origin[0] = vx[nprim];
313 0 : origin[1] = vy[nprim];
314 0 : origin[2] = -2196.5;
315 :
316 : //
317 : //
318 : // Particle weight
319 :
320 0 : Float_t originP[3] = {0., 0., 0.};
321 0 : originP[2] = zPrimary[nprim];
322 :
323 0 : Float_t pP[3] = {0., 0., 0.};
324 0 : Int_t ntP;
325 :
326 0 : if (fSide == -1) {
327 0 : originP[2] = -zPrimary[nprim];
328 0 : origin[2] = -origin[2];
329 0 : p[2] = -p[2];
330 0 : }
331 :
332 : //
333 : // Time
334 : //
335 0 : if (zPrimary[nprim] != zVertexOld) {
336 0 : while(arg==0.) arg = gRandom->Rndm();
337 0 : t -= dT*TMath::Log(arg); // (sec)
338 0 : zVertexOld = zPrimary[nprim];
339 0 : }
340 :
341 : // Get statistical weight according to local gas-pressure
342 : //
343 0 : fParentWeight=wgt[nprim]*GasPressureWeight(zPrimary[nprim]);
344 :
345 0 : if (!fAnalog || gRandom->Rndm() < fParentWeight) {
346 : // Pass parent particle
347 : //
348 0 : PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
349 0 : KeepTrack(ntP);
350 0 : PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
351 0 : }
352 : //
353 : // Both sides are considered
354 : //
355 :
356 0 : if (fSide > 1) {
357 0 : fParentWeight=wgt[nprim]*GasPressureWeight(-zPrimary[nprim]);
358 0 : if (!fAnalog || gRandom->Rndm() < fParentWeight) {
359 0 : origin[2] = -origin[2];
360 0 : originP[2] = -originP[2];
361 0 : p[2]=-p[2];
362 0 : PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
363 0 : KeepTrack(ntP);
364 0 : PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
365 0 : }
366 : }
367 0 : wwgt += fParentWeight;
368 :
369 0 : SetHighWaterMark(nt);
370 0 : }
371 0 : delete [] zPrimary;
372 0 : delete [] inuc;
373 0 : delete [] ipart;
374 0 : delete [] wgt;
375 0 : delete [] ekin;
376 0 : delete [] vx;
377 0 : delete [] vy;
378 0 : delete [] tx;
379 0 : delete [] ty;
380 0 : printf("Total weight %f\n\n", wwgt);
381 :
382 0 : }
383 :
384 :
385 : Float_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 0 : zPrimary /= 100.; // m
392 0 : if (fRunPeriod < 5) {
393 0 : Float_t zAbs = TMath::Abs(zPrimary);
394 0 : if (zPrimary > 0.)
395 : {
396 0 : if (zAbs > fZ1[20]) {
397 : weight = 2.e4;
398 0 : } else {
399 0 : for (Int_t i = 1; i < 21; i++) {
400 0 : if (zAbs < fZ1[i]) {
401 0 : weight = fG1[i];
402 0 : break;
403 : }
404 : }
405 : }
406 : } else {
407 0 : if (zAbs > fZ2[20]) {
408 : weight = 2.e4;
409 0 : } else {
410 0 : for (Int_t i = 1; i < 21; i++) {
411 0 : if (zAbs < fZ2[i]) {
412 0 : weight = fG2[i];
413 0 : break;
414 : }
415 : }
416 : }
417 : }
418 0 : } else {
419 0 : Int_t index = TMath::BinarySearch(fGPASize, fZ1, zPrimary);
420 0 : weight = fG1[index];
421 : }
422 0 : return weight;
423 : }
424 :
425 : void AliGenHaloProtvino::Draw(Option_t *)
426 : {
427 : // Draws the gas pressure distribution
428 0 : Float_t z[400];
429 0 : Float_t p[400];
430 :
431 0 : for (Int_t i = 0; i < 400; i++)
432 : {
433 0 : z[i] = -20000. + Float_t(i) * 100;
434 0 : p[i] = GasPressureWeight(z[i]);
435 : }
436 :
437 0 : TGraph* gr = new TGraph(400, z, p);
438 0 : TCanvas* c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
439 0 : c1->cd();
440 0 : gr->Draw("AL");
441 0 : }
442 :
443 :
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 :
449 : 0. 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 :
456 : 1. 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 :
472 : 2. 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 :
|
