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