099c8abee4406f8ebc93ae40e23bb6ae2c2d5fe6
[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   Float_t dT = fTimePerEvent/nInt;   // sec 
274   Float_t t  = 0;                    // sec
275   
276 //
277 // Loop over primaries
278 //
279   zVertexOld   = -1.e10;
280   Double_t arg = 0.;
281   
282   for (Int_t nprim = 0; nprim < fNpart; nprim++) 
283   {
284       amass = TDatabasePDG::Instance()->GetParticle(ipart[nprim])->Mass();
285
286       //
287       // Momentum vector
288       //
289       p0=sqrt(ekin[nprim]*ekin[nprim] + 2.*amass*ekin[nprim]);
290       
291       txy=TMath::Sqrt(tx[nprim]*tx[nprim]+ty[nprim]*ty[nprim]);
292       if (txy == 1.) {
293           tz=0;
294       } else {
295           tz=-TMath::Sqrt(1.-txy);
296       }
297     
298       p[0] = p0*tx[nprim];
299       p[1] = p0*ty[nprim];
300       p[2] =-p0*tz;
301       
302       origin[0] = vx[nprim];
303       origin[1] = vy[nprim];
304       origin[2] = -2196.5;
305
306       //
307       //
308       // Particle weight
309
310       Float_t originP[3] = {0., 0., 0.};
311       originP[2] = zPrimary[nprim];
312       
313       Float_t pP[3] = {0., 0., 0.};
314       Int_t ntP;
315       
316       if (fSide == -1) {
317           originP[2] = -zPrimary[nprim];
318           origin[2]  = -origin[2];
319           p[2]       = -p[2];
320       }
321
322       //
323       // Time
324       //
325       if (zPrimary[nprim] != zVertexOld) {
326           while(arg==0.) arg = gRandom->Rndm();
327           t -= dT*TMath::Log(arg);              // (sec)   
328           zVertexOld = zPrimary[nprim];
329       }
330
331 //    Get statistical weight according to local gas-pressure
332 //
333       fParentWeight=wgt[nprim]*GasPressureWeight(zPrimary[nprim]);
334
335       if (!fAnalog || gRandom->Rndm() < fParentWeight) {
336 //    Pass parent particle
337 //
338           PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
339           KeepTrack(ntP);
340           PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
341       }
342       //
343       // Both sides are considered
344       //
345
346       if (fSide > 1) {
347           fParentWeight=wgt[nprim]*GasPressureWeight(-zPrimary[nprim]);
348           if (!fAnalog || gRandom->Rndm() < fParentWeight) {
349               origin[2]  = -origin[2];
350               originP[2] = -originP[2];
351               p[2]=-p[2];
352               PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
353               KeepTrack(ntP);
354               PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
355           }
356       }
357       wwgt += fParentWeight;
358       
359       SetHighWaterMark(nt);
360   }
361   delete [] zPrimary;
362   delete [] inuc;    
363   delete [] ipart;   
364   delete [] wgt;     
365   delete [] ekin;    
366   delete [] vx;      
367   delete [] vy;      
368   delete [] tx;      
369   delete [] ty;     
370   printf("Total weight %f\n\n", wwgt);
371   
372 }
373  
374
375 Float_t AliGenHaloProtvino::GasPressureWeight(Float_t zPrimary)
376 {
377 //
378 // Return z-dependent gasspressure weight = interaction rate [1/m/s].
379 //
380     Float_t weight = 0.;
381     zPrimary /= 100.;        // m
382     if (fRunPeriod < 5) {
383         Float_t zAbs = TMath::Abs(zPrimary);
384         if (zPrimary > 0.) 
385         {
386             if (zAbs > fZ1[20]) {
387                 weight = 2.e4;
388             } else {
389                 for (Int_t i = 1; i < 21; i++) {
390                     if (zAbs < fZ1[i]) {
391                         weight = fG1[i];
392                         break;
393                     }
394                 }
395             }
396         } else {
397             if (zAbs > fZ2[20]) {
398                 weight = 2.e4;
399             } else {
400                 for (Int_t i = 1; i < 21; i++) {
401                     if (zAbs < fZ2[i]) {
402                     weight = fG2[i];
403                     break;
404                     }
405                 }
406             }
407         }
408     } else {
409         Int_t index = TMath::BinarySearch(fGPASize, fZ1, zPrimary);
410         weight = fG1[index];
411     }
412     return weight;
413 }
414
415 void AliGenHaloProtvino::Draw(Option_t *)
416 {
417 // Draws the gas pressure distribution
418     Float_t z[400];
419     Float_t p[400];
420     
421     for (Int_t i = 0; i < 400; i++)
422     {
423         z[i] = -20000. + Float_t(i) * 100;
424         p[i] = GasPressureWeight(z[i]);
425     }
426     
427     TGraph*  gr = new TGraph(400, z, p);   
428     TCanvas* c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
429     c1->cd();
430     gr->Draw("AL");
431 }
432
433
434 /*
435 # Title:    README file for the sources of IR8 machine induced background
436 # Author:   Vadim Talanov <Vadim.Talanov@cern.ch>
437 # Modified: 12-12-2000 
438
439 0. Overview
440
441         There are three files, named ring.one.beta.[01,10,50].m, which
442         contain the lists of background particles, induced by proton losses
443         upstream of IP8 in the LHC ring one, for the beta* values of 1, 10
444         and 50 m, respectively.
445
446 1. File contents
447
448         Each line in the files contains the coordinates of particle track
449         crossing with the infinite plane, positioned at z=-1m, together with
450         the physical properties of corresponding particle, namely:
451
452         S  - S coordinate of the primary interaction vertex, cm;
453         N  - type of the gas nuclei at interaction, 1 is H, 2 - C and 3 - O;
454         I  - particle ID in PDG particle numbering scheme;
455         W  - particle weight;
456         E  - particle kinetic energy, GeV;
457         X  - x coordinate of the crossing point, cm;
458         Y  - y coordinate of the crossing point, cm;
459         Dx - x direction cosine;
460         Dy - y direction cosine.
461
462 2. Normalisation
463
464         Each file is given per unity of linear density of proton inelastic
465         interactions with the gas nuclei, [1 inelastic interaction/m].
466
467 # ~/vtalanov/public/README.mib: the end.
468
469 */
470
471
472
473