Class to read new beam halo boundary source.
[u/mrichter/AliRoot.git] / EVGEN / AliGenHalo.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 beam halo background particles from a boundary source
19 // Boundary source is in the LHCb format http://www.hep.manchester.ac.uk/u/robert/LHC_backgrounds/Note-MIBStudies.pdf
20 // and has been provided by Robert Appleby
21 // Author: andreas.morsch@cern.ch
22
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 "AliGenHalo.h"
33 #include "AliRun.h"
34
35 ClassImp(AliGenHalo)
36
37 AliGenHalo::AliGenHalo()
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      fLossID(0),   
51      fLossA(0),   
52      fPdg(0),     
53      fLossT0(0),
54      fLossZ(0), 
55      fLossW(0), 
56      fXS(0),    
57      fYS(0),    
58      fZS(0),    
59      fDX(0),    
60      fDY(0),    
61      fEkin(0),  
62      fTS(0),    
63      fWS(0)     
64 {
65 // Constructor
66     
67     fName  = "Halo";
68     fTitle = "Halo from LHC Beam";
69 //
70 //  Read all particles
71     fNpart = -1;
72     SetAnalog(0);
73 }
74
75 AliGenHalo::AliGenHalo(Int_t npart)
76     :AliGenerator(npart),
77      fFile(0),
78      fFileName(0),
79      fSide(1),
80      fRunPeriod(kY3D90),
81      fTimePerEvent(1.e-4),
82      fNskip(0),
83      fZ1(0),
84      fZ2(0),
85      fG1(0),
86      fG2(0),
87      fGPASize(0),
88      fLossID(0),   
89      fLossA(0),   
90      fPdg(0),     
91      fLossT0(0),
92      fLossZ(0), 
93      fLossW(0), 
94      fXS(0),    
95      fYS(0),    
96      fZS(0),    
97      fDX(0),    
98      fDY(0),    
99      fEkin(0),  
100      fTS(0),    
101      fWS(0)     
102 {
103 // Constructor
104     fName = "Halo";
105     fTitle= "Halo from LHC Beam";
106 //
107     fNpart   = npart;
108 //
109     SetAnalog(0);
110 }
111
112 //____________________________________________________________
113 AliGenHalo::~AliGenHalo()
114 {
115 // Destructor
116 }
117
118 //____________________________________________________________
119 void AliGenHalo::Init() 
120 {
121 // Initialisation
122     fFile = fopen(fFileName,"r");
123     if (fFile) {
124         printf("\n File %s opened for reading, %p ! \n ",  fFileName.Data(), (void*)fFile);
125     } else {
126         printf("\n Opening of file %s failed,  %p ! \n ",  fFileName.Data(), (void*)fFile);
127     }
128 //
129 //
130 //
131 //    Read file with gas pressure values
132     char *name = 0;
133     if (fRunPeriod < 5) {
134         name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
135         fGPASize = 21;
136         fG1 = new Float_t[fGPASize];
137         fG2 = new Float_t[fGPASize];
138         fZ1 = new Float_t[fGPASize];
139         fZ2 = new Float_t[fGPASize];
140     } else if (fRunPeriod == 5) {
141         name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_startup.dat");
142         fGPASize = 18853;
143         fG1 = new Float_t[fGPASize];
144         fZ1 = new Float_t[fGPASize];
145     } else if (fRunPeriod ==6) {
146         name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_conditioned.dat");
147         fGPASize = 12719;
148         fG1 = new Float_t[fGPASize];
149         fZ1 = new Float_t[fGPASize];
150     } else {
151         Fatal("Init()", "No gas pressure file for given run period !");
152     }
153
154
155     FILE* file = fopen(name, "r");
156     Float_t z;
157     Int_t i;
158     Float_t p[5];    
159
160     const Float_t kCrossSection = 0.094e-28;      // m^2
161     const Float_t kFlux         = 1.e11 / 25.e-9; // protons/s
162     Float_t pFlux[5] = {0.2, 0.2, 0.3, 0.3, 1.0};
163
164     if (fRunPeriod < 5) {
165 //
166 //  Ring 1   
167 // 
168
169         for (i = 0; i < fGPASize; i++)
170         {
171             fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
172             fG1[i] = p[fRunPeriod];
173             
174             if (i > 0) {
175                 fZ1[i] = fZ1[i-1] + z;
176             } else {
177                 fZ1[i] = 20.;
178             }
179         }
180 //
181 // Ring 2
182 //
183         for (i = 0; i < fGPASize; i++)
184         {
185             fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
186             fG2[i] = p[fRunPeriod];
187             if (i > 0) {
188                 fZ2[i] = fZ2[i-1] + z;
189             } else {
190                 fZ2[i] = 20.;
191             }
192         }
193 //
194 // Interaction rates
195 //
196         for (i = 0; i < fGPASize; i++)  
197         {
198             fG1[i] = fG1[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s 
199             fG2[i] = fG2[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
200         }
201
202     } else {
203         for (i = 0; i < fGPASize; i++) 
204         {
205             fscanf(file, "%f %e %e %e %e %e", &z, &p[0], &p[1], &p[2], &p[3], &p[4]);
206             z /= 1000.;
207             fG1[i] = p[4] * kCrossSection * kFlux;             // 1/m/s
208             // 1/3 of nominal intensity at startup
209             if (fRunPeriod ==  kLHCPR674Startup) fG1[i] /= 3.;
210             fZ1[i] = z;
211         }
212     }
213     
214
215
216     
217 //
218 //  Transform into interaction rates
219 //
220
221
222     
223
224     Float_t sum1 = 0.;
225     Float_t sum2 = 0.;
226     
227     for (Int_t iz = 0; iz < 300; iz++) {
228         Float_t zpos = 20. + iz * 1.;
229         zpos *= 100;
230         Float_t wgt1 = GasPressureWeight( zpos);
231         Float_t wgt2 = GasPressureWeight(-zpos);
232         sum1 += wgt1;
233         sum2 += wgt2;
234     }
235     sum1/=250.;
236     sum2/=250.;
237     printf("\n %f %f \n \n", sum1, sum2);
238 }
239
240 //____________________________________________________________
241 void AliGenHalo::Generate()
242 {
243 // Generate by reading particles from input file
244  
245   Float_t polar[3]= {0., 0., 0.};
246   Float_t origin[3];
247   Float_t p[3], p0;
248   Float_t tz, txy;
249   Float_t mass;
250   //
251   Int_t nt;
252   static Bool_t first = kTRUE;
253   static Int_t  oldID = -1;
254 //
255   if (first) {
256       ReadNextParticle();
257       first = kFALSE;
258       oldID = fLossID;
259   }
260   
261   
262   
263   while(1) {
264       // Push particle to stack
265       mass = TDatabasePDG::Instance()->GetParticle(fPdg)->Mass();
266       p0  = TMath::Sqrt(fEkin * fEkin + 2.* mass * fEkin);
267       txy = TMath::Sqrt(fDX * fDX + fDY * fDY);
268       if (txy == 1.) {
269           tz = 0;
270       } else {
271           tz = - TMath::Sqrt(1. - txy);
272       }
273       p[0] =  p0 * fDX;
274       p[1] =  p0 * fDY;
275       p[2] =  p0 * tz;
276
277       origin[0] = fXS;
278       origin[1] = fYS;
279       origin[2] = 1950.;
280
281       PushTrack(fTrackIt , -1, fPdg , p, origin, polar, fTS, kPNoProcess, nt, fWS);
282
283       Int_t nc = ReadNextParticle();
284       
285       if (fLossID != oldID || nc == 0) {
286           oldID = fLossID;
287           break;
288       }
289   }
290   SetHighWaterMark(nt);
291 }
292  
293
294 Float_t AliGenHalo::GasPressureWeight(Float_t zPrimary)
295 {
296 //
297 // Return z-dependent gasspressure weight = interaction rate [1/m/s].
298 //
299     Float_t weight = 0.;
300     zPrimary /= 100.;        // m
301     if (fRunPeriod < 5) {
302         Float_t zAbs = TMath::Abs(zPrimary);
303         if (zPrimary > 0.) 
304         {
305             if (zAbs > fZ1[20]) {
306                 weight = 2.e4;
307             } else {
308                 for (Int_t i = 1; i < 21; i++) {
309                     if (zAbs < fZ1[i]) {
310                         weight = fG1[i];
311                         break;
312                     }
313                 }
314             }
315         } else {
316             if (zAbs > fZ2[20]) {
317                 weight = 2.e4;
318             } else {
319                 for (Int_t i = 1; i < 21; i++) {
320                     if (zAbs < fZ2[i]) {
321                     weight = fG2[i];
322                     break;
323                     }
324                 }
325             }
326         }
327     } else {
328         Int_t index = TMath::BinarySearch(fGPASize, fZ1, zPrimary);
329         weight = fG1[index];
330     }
331     return weight;
332 }
333
334 void AliGenHalo::Draw(Option_t *)
335 {
336 // Draws the gas pressure distribution
337     Float_t z[400];
338     Float_t p[400];
339     
340     for (Int_t i = 0; i < 400; i++)
341     {
342         z[i] = -20000. + Float_t(i) * 100;
343         p[i] = GasPressureWeight(z[i]);
344     }
345     
346     TGraph*  gr = new TGraph(400, z, p);   
347     TCanvas* c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
348     c1->cd();
349     gr->Draw("AL");
350 }
351
352 Int_t AliGenHalo::ReadNextParticle()
353 {
354     // Read the next particle from the file
355     Int_t ncols = fscanf(fFile,"%d %f %f %d %f %d %f %f %f %f %f %f %f %f",
356                    &fLossID, &fLossT0, &fLossZ, &fLossA, &fLossW, &fPdg, &fXS, &fYS, &fZS, &fDX, &fDY, &fEkin, &fTS, &fWS);
357     fLossZ /= 10.;
358     fXS    /= 10.;
359     fYS    /= 10.; 
360     fZS    /= 10.;   
361     fTS    *= 1.e-9;
362     return (ncols);
363 }