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