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