Transition to NewIO
[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 <TPDGCode.h>
28 #include <TSystem.h>
29
30 #include "AliGenHaloProtvino.h"
31 #include "AliRun.h"
32
33  ClassImp(AliGenHaloProtvino)
34      AliGenHaloProtvino::AliGenHaloProtvino()
35          :AliGenerator(-1)
36 {
37 // Constructor
38     
39     fName  = "HaloProtvino";
40     fTitle = "Halo from LHC Tunnel";
41 //
42 //  Read all particles
43     fNpart = -1;
44     fFile  =  0;
45     fSide  =  1;
46 //
47     SetRunPeriod();
48     SetTimePerEvent();
49     SetAnalog(0);
50 }
51
52 AliGenHaloProtvino::AliGenHaloProtvino(Int_t npart)
53     :AliGenerator(npart)
54 {
55 // Constructor
56     fName = "Halo";
57     fTitle= "Halo from LHC Tunnel";
58 //
59     fNpart   = npart;
60     fFile    = 0;
61     fSide    = 1;
62 //
63     SetRunPeriod();
64     SetTimePerEvent();
65     SetAnalog(0);
66 }
67
68 AliGenHaloProtvino::AliGenHaloProtvino(const AliGenHaloProtvino & HaloProtvino)
69 {
70 // copy constructor
71 }
72
73
74 //____________________________________________________________
75 AliGenHaloProtvino::~AliGenHaloProtvino()
76 {
77 // Destructor
78 }
79
80 //____________________________________________________________
81 void AliGenHaloProtvino::Init() 
82 {
83 // Initialisation
84     fFile = fopen(fFileName,"r");
85     if (fFile) {
86         printf("\n File %s opened for reading, %p ! \n ",  fFileName.Data(), fFile);
87     } else {
88         printf("\n Opening of file %s failed,  %p ! \n ",  fFileName.Data(), fFile);
89     }
90 //
91 //
92 //
93 //    Read file with gas pressure values
94     char *name;
95     
96     name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
97     FILE* file = fopen(name, "r");
98     Float_t z;
99     Int_t i, j;
100     
101 //
102 //  Ring 1   
103 // 
104     for (i = 0; i < 21; i++)
105     {
106         fscanf(file, "%f %f %f %f %f %f", &z, 
107                &fG1[i][0], &fG1[i][1], &fG1[i][2], &fG1[i][3], &fG1[i][4]);
108         if (i > 0) {
109             fZ1[i] = fZ1[i-1] + z;
110         } else {
111             fZ1[i] = 20.;
112         }
113   }
114 //
115 // Ring 2
116 //
117     for (i = 0; i < 21; i++)
118     {
119         fscanf(file, "%f %f %f %f %f %f", &z, 
120                &fG2[i][0], &fG2[i][1], &fG2[i][2], &fG2[i][3], &fG2[i][4]);
121         if (i > 0) {
122             fZ2[i] = fZ2[i-1] + z;
123         } else {
124             fZ2[i] = 20.;
125         }
126     }
127 //
128 //  Transform into interaction rates
129 //
130     const Float_t crossSection = 0.094e-28;     // m^2
131     Float_t pFlux[5] = {0.2, 0.2, 0.3, 0.3, 1.0};
132
133     for (j = 0; j <  5; j++) {
134         pFlux[j] *= 1.e11/25.e-9;
135         for (i = 0; i < 21; i++)  
136         {
137             fG1[i][j] = fG1[i][j] * crossSection * pFlux[j]; // 1/m/s 
138             fG2[i][j] = fG2[i][j] * crossSection * pFlux[j]; // 1/m/s
139         }
140     }
141     
142
143     Float_t sum1 = 0.;
144     Float_t sum2 = 0.;
145     
146     for (Int_t i = 0; i < 300; i++) {
147         Float_t z = 20.+i*1.;
148         z*=100;
149         Float_t wgt1 = GassPressureWeight(z);
150         Float_t wgt2 = GassPressureWeight(-z);
151 //      printf("weight: %f %f %f %f %f \n", z, wgt1, wgt2, fZ1[20], fZ2[20]);
152         sum1 += wgt1;
153         sum2 += wgt2;
154     }
155     sum1/=250.;
156     sum2/=250.;
157     printf("\n %f %f \n \n", sum1, sum2);
158 }
159
160 //____________________________________________________________
161 void AliGenHaloProtvino::Generate()
162 {
163 // Generate from input file
164  
165   Float_t polar[3]= {0,0,0};
166   Float_t origin[3];
167   Float_t p[3], p0;
168   Float_t tz, txy;
169   Float_t amass;
170   //
171   Int_t ncols, nt;
172   Int_t nskip = 0;
173   Int_t nread = 0;
174
175   Float_t* zPrimary = new Float_t [fNpart];
176   Int_t  * inuc     = new Int_t   [fNpart];
177   Int_t  * ipart    = new Int_t   [fNpart];
178   Float_t* wgt      = new Float_t [fNpart]; 
179   Float_t* ekin     = new Float_t [fNpart];
180   Float_t* vx       = new Float_t [fNpart];
181   Float_t* vy       = new Float_t [fNpart];
182   Float_t* tx       = new Float_t [fNpart];
183   Float_t* ty       = new Float_t [fNpart];
184   
185   Float_t zVertexOld = -1.e10;
186   Int_t   nInt       = 0;        // Counts number of interactions
187   
188   while(1) {
189 //
190 // Load event into array
191 //
192       ncols = fscanf(fFile,"%f %d %d %f %f %f %f %f %f",
193                      &zPrimary[nread], &inuc[nread], &ipart[nread], &wgt[nread], 
194                      &ekin[nread], &vx[nread], &vy[nread],
195                      &tx[nread], &ty[nread]);
196       
197       if (ncols < 0) break;
198 // Skip fNskip events
199       nskip++;
200       if (fNpart !=-1 && nskip <= fNskip) continue;
201 // Count interactions
202       if (zPrimary[nread] != zVertexOld) {
203           nInt++;
204           zVertexOld = zPrimary[nread];
205       }
206 // Count tracks      
207       nread++;
208       if (fNpart !=-1 && nread > fNpart) break;
209   }
210 //
211 // Mean time between interactions
212 //
213   Float_t dT = fTimePerEvent/nInt;   // sec 
214   Float_t t  = 0;                    // sec
215   
216 //
217 // Loop over primaries
218 //
219   zVertexOld   = -1.e10;
220   Double_t arg = 0.;
221   
222   for (Int_t nprim = 0; nprim < fNpart; nprim++) 
223   {
224       amass = TDatabasePDG::Instance()->GetParticle(ipart[nprim])->Mass();
225
226       //
227       // Momentum vector
228       //
229       p0=sqrt(ekin[nprim]*ekin[nprim] + 2.*amass*ekin[nprim]);
230       
231       txy=TMath::Sqrt(tx[nprim]*tx[nprim]+ty[nprim]*ty[nprim]);
232       if (txy == 1.) {
233           tz=0;
234       } else {
235           tz=-TMath::Sqrt(1.-txy);
236       }
237     
238       p[0] = p0*tx[nprim];
239       p[1] = p0*ty[nprim];
240       p[2] =-p0*tz;
241       
242       origin[0] = vx[nprim];
243       origin[1] = vy[nprim];
244       origin[2] = -2196.5;
245
246       //
247       //
248       // Particle weight
249
250       Float_t originP[3] = {0., 0., 0.};
251       originP[2] = zPrimary[nprim];
252       
253       Float_t pP[3] = {0., 0., 0.};
254       Int_t ntP;
255       
256       if (fSide == -1) {
257           originP[2] = -zPrimary[nprim];
258           origin[2]  = -origin[2];
259           p[2]       = -p[2];
260       }
261
262       //
263       // Time
264       //
265       if (zPrimary[nprim] != zVertexOld) {
266           while(arg==0.) arg = gRandom->Rndm();
267           t -= dT*TMath::Log(arg);              // (sec)   
268           zVertexOld = zPrimary[nprim];
269       }
270
271 //    Get statistical weight according to local gas-pressure
272 //
273       fParentWeight=wgt[nprim]*GassPressureWeight(zPrimary[nprim]);
274
275       if (!fAnalog || gRandom->Rndm() < fParentWeight) {
276 //    Pass parent particle
277 //
278           SetTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
279           KeepTrack(ntP);
280           SetTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
281       }
282
283       //
284       // Both sides are considered
285       //
286
287       if (fSide > 1) {
288           fParentWeight=wgt[nprim]*GassPressureWeight(-zPrimary[nprim]);
289           if (!fAnalog || gRandom->Rndm() < fParentWeight) {
290               origin[2]  = -origin[2];
291               originP[2] = -originP[2];
292               p[2]=-p[2];
293               SetTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
294               KeepTrack(ntP);
295               SetTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
296           }
297       }
298       SetHighWaterMark(nt);
299   }
300   delete [] zPrimary;
301   delete [] inuc;    
302   delete [] ipart;   
303   delete [] wgt;     
304   delete [] ekin;    
305   delete [] vx;      
306   delete [] vy;      
307   delete [] tx;      
308   delete [] ty;      
309 }
310  
311
312 AliGenHaloProtvino& AliGenHaloProtvino::operator=(const  AliGenHaloProtvino& rhs)
313 {
314 // Assignment operator
315     return *this;
316 }
317
318
319
320 Float_t AliGenHaloProtvino::GassPressureWeight(Float_t zPrimary)
321 {
322 //
323 // Return z-dependent gasspressure weight = interaction rate [1/m/s].
324 //
325     Float_t weight = 0.;
326     zPrimary/=100.;        // m
327     Float_t zAbs = TMath::Abs(zPrimary);
328     if (zPrimary > 0.) 
329     {
330         if (zAbs > fZ1[20]) {
331             weight = 2.e4;
332         } else {
333             for (Int_t i = 1; i < 21; i++) {
334                 if (zAbs < fZ1[i]) {
335                     weight = fG1[i][fRunPeriod];
336                     break;
337                 }
338             }
339         }
340     } else {
341         if (zAbs > fZ2[20]) {
342             weight = 2.e4;
343         } else {
344             for (Int_t i = 1; i < 21; i++) {
345                 if (zAbs < fZ2[i]) {
346                     weight = fG2[i][fRunPeriod];
347                     break;
348                 }
349             }
350         }
351     }
352     return weight;
353 }
354
355 /*
356 # Title:    README file for the sources of IR8 machine induced background
357 # Author:   Vadim Talanov <Vadim.Talanov@cern.ch>
358 # Modified: 12-12-2000 
359
360 0. Overview
361
362         There are three files, named ring.one.beta.[01,10,50].m, which
363         contain the lists of background particles, induced by proton losses
364         upstream of IP8 in the LHC ring one, for the beta* values of 1, 10
365         and 50 m, respectively.
366
367 1. File contents
368
369         Each line in the files contains the coordinates of particle track
370         crossing with the infinite plane, positioned at z=-1m, together with
371         the physical properties of corresponding particle, namely:
372
373         S  - S coordinate of the primary interaction vertex, cm;
374         N  - type of the gas nuclei at interaction, 1 is H, 2 - C and 3 - O;
375         I  - particle ID in PDG particle numbering scheme;
376         W  - particle weight;
377         E  - particle kinetic energy, GeV;
378         X  - x coordinate of the crossing point, cm;
379         Y  - y coordinate of the crossing point, cm;
380         Dx - x direction cosine;
381         Dy - y direction cosine.
382
383 2. Normalisation
384
385         Each file is given per unity of linear density of proton inelastic
386         interactions with the gas nuclei, [1 inelastic interaction/m].
387
388 # ~/vtalanov/public/README.mib: the end.
389
390 */
391
392
393
394