0a6041e434412212f402629c6a41eb351c4c861d
[u/mrichter/AliRoot.git] / EVGEN / AliGenHaloProtvino.cxx
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 ClassImp(AliGenHaloProtvino)
37
38 AliGenHaloProtvino::AliGenHaloProtvino()
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 {
52 // Constructor
53     
54     fName  = "HaloProtvino";
55     fTitle = "Halo from LHC Tunnel";
56 //
57 //  Read all particles
58     fNpart = -1;
59     SetAnalog(0);
60 }
61
62 AliGenHaloProtvino::AliGenHaloProtvino(Int_t npart)
63     :AliGenerator(npart),
64      fFile(0),
65      fFileName(0),
66      fSide(1),
67      fRunPeriod(kY3D90),
68      fTimePerEvent(1.e-4),
69      fNskip(0),
70      fZ1(0),
71      fZ2(0),
72      fG1(0),
73      fG2(0),
74      fGPASize(0)
75 {
76 // Constructor
77     fName = "Halo";
78     fTitle= "Halo from LHC Tunnel";
79 //
80     fNpart   = npart;
81 //
82     SetAnalog(0);
83 }
84
85 //____________________________________________________________
86 AliGenHaloProtvino::~AliGenHaloProtvino()
87 {
88 // Destructor
89 }
90
91 //____________________________________________________________
92 void AliGenHaloProtvino::Init() 
93 {
94 // Initialisation
95     fFile = fopen(fFileName,"r");
96     if (fFile) {
97         printf("\n File %s opened for reading, %p ! \n ",  fFileName.Data(), (void*)fFile);
98     } else {
99         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     if (fRunPeriod < 5) {
107         name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
108         fGPASize = 21;
109         fG1 = new Float_t[fGPASize];
110         fG2 = new Float_t[fGPASize];
111         fZ1 = new Float_t[fGPASize];
112         fZ2 = new Float_t[fGPASize];
113     } else if (fRunPeriod == 5) {
114         name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_startup.dat");
115         fGPASize = 18853;
116         fG1 = new Float_t[fGPASize];
117         fZ1 = new Float_t[fGPASize];
118     } else if (fRunPeriod ==6) {
119         name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_conditioned.dat");
120         fGPASize = 12719;
121         fG1 = new Float_t[fGPASize];
122         fZ1 = new Float_t[fGPASize];
123     } else {
124         Fatal("Init()", "No gas pressure file for given run period !");
125     }
126
127     FILE* file = 0;
128     if (name) file = fopen(name, "r");
129     if (!file) {
130         AliError("No gas pressure file");
131         return;
132     }
133
134     Float_t z;
135     Int_t i;
136     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     Float_t pFlux[5] = {0.2, 0.2, 0.3, 0.3, 1.0};
141
142     if (fRunPeriod < 5) {
143 //
144 //  Ring 1   
145 // 
146
147         for (i = 0; i < fGPASize; i++)
148         {
149             fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
150             fG1[i] = p[fRunPeriod];
151             
152             if (i > 0) {
153                 fZ1[i] = fZ1[i-1] + z;
154             } else {
155                 fZ1[i] = 20.;
156             }
157         }
158 //
159 // Ring 2
160 //
161         for (i = 0; i < fGPASize; i++)
162         {
163             fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
164             fG2[i] = p[fRunPeriod];
165             if (i > 0) {
166                 fZ2[i] = fZ2[i-1] + z;
167             } else {
168                 fZ2[i] = 20.;
169             }
170         }
171 //
172 // Interaction rates
173 //
174         for (i = 0; i < fGPASize; i++)  
175         {
176             fG1[i] = fG1[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s 
177             fG2[i] = fG2[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
178         }
179
180     } else {
181         for (i = 0; i < fGPASize; i++) 
182         {
183             fscanf(file, "%f %e %e %e %e %e", &z, &p[0], &p[1], &p[2], &p[3], &p[4]);
184             z /= 1000.;
185             fG1[i] = p[4] * kCrossSection * kFlux;             // 1/m/s
186             // 1/3 of nominal intensity at startup
187             if (fRunPeriod ==  kLHCPR674Startup) fG1[i] /= 3.;
188             fZ1[i] = z;
189         }
190     }
191     
192
193
194     
195 //
196 //  Transform into interaction rates
197 //
198
199
200     
201
202     Float_t sum1 = 0.;
203     Float_t sum2 = 0.;
204     
205     for (Int_t iz = 0; iz < 300; iz++) {
206         Float_t zpos = 20. + iz * 1.;
207         zpos *= 100;
208         Float_t wgt1 = GasPressureWeight( zpos);
209         Float_t wgt2 = GasPressureWeight(-zpos);
210         sum1 += wgt1;
211         sum2 += wgt2;
212     }
213     sum1/=250.;
214     sum2/=250.;
215     printf("\n %f %f \n \n", sum1, sum2);
216     delete file;
217 }
218
219 //____________________________________________________________
220 void AliGenHaloProtvino::Generate()
221 {
222 // Generate from input file
223  
224   Float_t polar[3]= {0,0,0};
225   Float_t origin[3];
226   Float_t p[3], p0;
227   Float_t tz, txy;
228   Float_t amass;
229   //
230   Int_t ncols, nt;
231   static Int_t nskip = 0;
232   Int_t nread = 0;
233
234   Float_t* zPrimary = new Float_t [fNpart];
235   Int_t  * inuc     = new Int_t   [fNpart];
236   Int_t  * ipart    = new Int_t   [fNpart];
237   Float_t* wgt      = new Float_t [fNpart]; 
238   Float_t* ekin     = new Float_t [fNpart];
239   Float_t* vx       = new Float_t [fNpart];
240   Float_t* vy       = new Float_t [fNpart];
241   Float_t* tx       = new Float_t [fNpart];
242   Float_t* ty       = new Float_t [fNpart];
243   
244   Float_t zVertexOld = -1.e10;
245   Int_t   nInt       = 0;        // Counts number of interactions
246   Float_t wwgt = 0.;
247   
248   while(1) {
249 //
250 // Load event into array
251 //
252       ncols = fscanf(fFile,"%f %d %d %f %f %f %f %f %f",
253                      &zPrimary[nread], &inuc[nread], &ipart[nread], &wgt[nread], 
254                      &ekin[nread], &vx[nread], &vy[nread],
255                      &tx[nread], &ty[nread]);
256       
257       if (ncols < 0) break;
258 // Skip fNskip events
259       nskip++;
260       if (fNpart !=-1 && nskip <= fNskip) continue;
261 // Count interactions
262       if (zPrimary[nread] != zVertexOld) {
263           nInt++;
264           zVertexOld = zPrimary[nread];
265       }
266 // Count tracks      
267       nread++;
268       if (fNpart !=-1 && nread >= fNpart) break;
269   }
270 //
271 // Mean time between interactions
272 //
273
274   Float_t dT = 0.;   // sec 
275   if (nInt > 0) 
276       dT = fTimePerEvent/nInt;   
277   Float_t t  = 0;                    // sec
278   
279 //
280 // Loop over primaries
281 //
282   zVertexOld   = -1.e10;
283   Double_t arg = 0.;
284   
285   for (Int_t nprim = 0; nprim < fNpart; nprim++) 
286   {
287       amass = TDatabasePDG::Instance()->GetParticle(ipart[nprim])->Mass();
288
289       //
290       // Momentum vector
291       //
292       p0=sqrt(ekin[nprim]*ekin[nprim] + 2.*amass*ekin[nprim]);
293       
294       txy=TMath::Sqrt(tx[nprim]*tx[nprim]+ty[nprim]*ty[nprim]);
295       if (txy == 1.) {
296           tz=0;
297       } else {
298           tz=-TMath::Sqrt(1.-txy);
299       }
300     
301       p[0] = p0*tx[nprim];
302       p[1] = p0*ty[nprim];
303       p[2] =-p0*tz;
304       
305       origin[0] = vx[nprim];
306       origin[1] = vy[nprim];
307       origin[2] = -2196.5;
308
309       //
310       //
311       // Particle weight
312
313       Float_t originP[3] = {0., 0., 0.};
314       originP[2] = zPrimary[nprim];
315       
316       Float_t pP[3] = {0., 0., 0.};
317       Int_t ntP;
318       
319       if (fSide == -1) {
320           originP[2] = -zPrimary[nprim];
321           origin[2]  = -origin[2];
322           p[2]       = -p[2];
323       }
324
325       //
326       // Time
327       //
328       if (zPrimary[nprim] != zVertexOld) {
329           while(arg==0.) arg = gRandom->Rndm();
330           t -= dT*TMath::Log(arg);              // (sec)   
331           zVertexOld = zPrimary[nprim];
332       }
333
334 //    Get statistical weight according to local gas-pressure
335 //
336       fParentWeight=wgt[nprim]*GasPressureWeight(zPrimary[nprim]);
337
338       if (!fAnalog || gRandom->Rndm() < fParentWeight) {
339 //    Pass parent particle
340 //
341           PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
342           KeepTrack(ntP);
343           PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
344       }
345       //
346       // Both sides are considered
347       //
348
349       if (fSide > 1) {
350           fParentWeight=wgt[nprim]*GasPressureWeight(-zPrimary[nprim]);
351           if (!fAnalog || gRandom->Rndm() < fParentWeight) {
352               origin[2]  = -origin[2];
353               originP[2] = -originP[2];
354               p[2]=-p[2];
355               PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
356               KeepTrack(ntP);
357               PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
358           }
359       }
360       wwgt += fParentWeight;
361       
362       SetHighWaterMark(nt);
363   }
364   delete [] zPrimary;
365   delete [] inuc;    
366   delete [] ipart;   
367   delete [] wgt;     
368   delete [] ekin;    
369   delete [] vx;      
370   delete [] vy;      
371   delete [] tx;      
372   delete [] ty;     
373   printf("Total weight %f\n\n", wwgt);
374   
375 }
376  
377
378 Float_t AliGenHaloProtvino::GasPressureWeight(Float_t zPrimary)
379 {
380 //
381 // Return z-dependent gasspressure weight = interaction rate [1/m/s].
382 //
383     Float_t weight = 0.;
384     zPrimary /= 100.;        // m
385     if (fRunPeriod < 5) {
386         Float_t zAbs = TMath::Abs(zPrimary);
387         if (zPrimary > 0.) 
388         {
389             if (zAbs > fZ1[20]) {
390                 weight = 2.e4;
391             } else {
392                 for (Int_t i = 1; i < 21; i++) {
393                     if (zAbs < fZ1[i]) {
394                         weight = fG1[i];
395                         break;
396                     }
397                 }
398             }
399         } else {
400             if (zAbs > fZ2[20]) {
401                 weight = 2.e4;
402             } else {
403                 for (Int_t i = 1; i < 21; i++) {
404                     if (zAbs < fZ2[i]) {
405                     weight = fG2[i];
406                     break;
407                     }
408                 }
409             }
410         }
411     } else {
412         Int_t index = TMath::BinarySearch(fGPASize, fZ1, zPrimary);
413         weight = fG1[index];
414     }
415     return weight;
416 }
417
418 void AliGenHaloProtvino::Draw(Option_t *)
419 {
420 // Draws the gas pressure distribution
421     Float_t z[400];
422     Float_t p[400];
423     
424     for (Int_t i = 0; i < 400; i++)
425     {
426         z[i] = -20000. + Float_t(i) * 100;
427         p[i] = GasPressureWeight(z[i]);
428     }
429     
430     TGraph*  gr = new TGraph(400, z, p);   
431     TCanvas* c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
432     c1->cd();
433     gr->Draw("AL");
434 }
435
436
437 /*
438 # Title:    README file for the sources of IR8 machine induced background
439 # Author:   Vadim Talanov <Vadim.Talanov@cern.ch>
440 # Modified: 12-12-2000 
441
442 0. Overview
443
444         There are three files, named ring.one.beta.[01,10,50].m, which
445         contain the lists of background particles, induced by proton losses
446         upstream of IP8 in the LHC ring one, for the beta* values of 1, 10
447         and 50 m, respectively.
448
449 1. File contents
450
451         Each line in the files contains the coordinates of particle track
452         crossing with the infinite plane, positioned at z=-1m, together with
453         the physical properties of corresponding particle, namely:
454
455         S  - S coordinate of the primary interaction vertex, cm;
456         N  - type of the gas nuclei at interaction, 1 is H, 2 - C and 3 - O;
457         I  - particle ID in PDG particle numbering scheme;
458         W  - particle weight;
459         E  - particle kinetic energy, GeV;
460         X  - x coordinate of the crossing point, cm;
461         Y  - y coordinate of the crossing point, cm;
462         Dx - x direction cosine;
463         Dy - y direction cosine.
464
465 2. Normalisation
466
467         Each file is given per unity of linear density of proton inelastic
468         interactions with the gas nuclei, [1 inelastic interaction/m].
469
470 # ~/vtalanov/public/README.mib: the end.
471
472 */
473
474
475
476