In AliMUONAlignment
[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     if (fNskip > 0) {
130       // Skip the first fNskip events
131       SkipEvents();
132     }
133 //
134 //
135 //
136 //    Read file with gas pressure values
137     char *name = 0;
138     if (fRunPeriod < 5) {
139         name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
140         fGPASize = 21;
141         fG1 = new Float_t[fGPASize];
142         fG2 = new Float_t[fGPASize];
143         fZ1 = new Float_t[fGPASize];
144         fZ2 = new Float_t[fGPASize];
145     } else if (fRunPeriod == 5) {
146         name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_startup.dat");
147         fGPASize = 18853;
148         fG1 = new Float_t[fGPASize];
149         fZ1 = new Float_t[fGPASize];
150     } else if (fRunPeriod ==6) {
151         name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_conditioned.dat");
152         fGPASize = 12719;
153         fG1 = new Float_t[fGPASize];
154         fZ1 = new Float_t[fGPASize];
155     } else {
156         Fatal("Init()", "No gas pressure file for given run period !");
157     }
158
159
160     FILE* file = fopen(name, "r");
161     Float_t z;
162     Int_t i;
163     Float_t p[5];    
164
165     const Float_t kCrossSection = 0.094e-28;      // m^2
166     const Float_t kFlux         = 1.e11 / 25.e-9; // protons/s
167     Float_t pFlux[5] = {0.2, 0.2, 0.3, 0.3, 1.0};
168
169     if (fRunPeriod < 5) {
170 //
171 //  Ring 1   
172 // 
173
174         for (i = 0; i < fGPASize; i++)
175         {
176             fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
177             fG1[i] = p[fRunPeriod];
178             
179             if (i > 0) {
180                 fZ1[i] = fZ1[i-1] + z;
181             } else {
182                 fZ1[i] = 20.;
183             }
184         }
185 //
186 // Ring 2
187 //
188         for (i = 0; i < fGPASize; i++)
189         {
190             fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
191             fG2[i] = p[fRunPeriod];
192             if (i > 0) {
193                 fZ2[i] = fZ2[i-1] + z;
194             } else {
195                 fZ2[i] = 20.;
196             }
197         }
198 //
199 // Interaction rates
200 //
201         for (i = 0; i < fGPASize; i++)  
202         {
203             fG1[i] = fG1[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s 
204             fG2[i] = fG2[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
205         }
206
207     } else {
208         for (i = 0; i < fGPASize; i++) 
209         {
210             fscanf(file, "%f %e %e %e %e %e", &z, &p[0], &p[1], &p[2], &p[3], &p[4]);
211             z /= 1000.;
212             fG1[i] = p[4] * kCrossSection * kFlux;             // 1/m/s
213             // 1/3 of nominal intensity at startup
214             if (fRunPeriod ==  kLHCPR674Startup) fG1[i] /= 3.;
215             fZ1[i] = z;
216         }
217     }
218     
219
220
221     
222 //
223 //  Transform into interaction rates
224 //
225
226
227     
228
229     Float_t sum1 = 0.;
230     Float_t sum2 = 0.;
231     
232     for (Int_t iz = 0; iz < 300; iz++) {
233         Float_t zpos = 20. + iz * 1.;
234         zpos *= 100;
235         Float_t wgt1 = GasPressureWeight( zpos);
236         Float_t wgt2 = GasPressureWeight(-zpos);
237         sum1 += wgt1;
238         sum2 += wgt2;
239     }
240     sum1/=250.;
241     sum2/=250.;
242     printf("\n %f %f \n \n", sum1, sum2);
243 }
244
245 //____________________________________________________________
246 void AliGenHalo::Generate()
247 {
248 // Generate by reading particles from input file
249  
250   Float_t polar[3]= {0., 0., 0.};
251   Float_t origin[3];
252   Float_t p[3], p0;
253   Float_t tz, txy;
254   Float_t mass;
255   //
256   Int_t nt;
257   static Bool_t first = kTRUE;
258   static Int_t  oldID = -1;
259 //
260
261   if (first && (fNskip == 0)) ReadNextParticle();
262   first = kFALSE;
263   oldID = fLossID;
264   
265   while(1) {
266       // Push particle to stack
267       mass = TDatabasePDG::Instance()->GetParticle(fPdg)->Mass();
268       p0  = TMath::Sqrt(fEkin * fEkin + 2. * mass * fEkin);
269       txy = TMath::Sqrt(fDX * fDX + fDY * fDY);
270       if (txy > 1.) {
271           tz = 0.;
272       } else {
273           tz = - TMath::Sqrt(1. - txy);
274       }
275  
276       p[0] =  p0 * fDX;
277       p[1] =  p0 * fDY;
278       p[2] =  p0 * tz;
279
280       origin[0] = fXS;
281       origin[1] = fYS;
282       origin[2] = 1950.;
283
284       PushTrack(fTrackIt , -1, fPdg , p, origin, polar, fTS - 1950./2.9979e10, kPNoProcess, nt, fWS);
285       
286       Int_t nc = ReadNextParticle();
287       
288       if (fLossID != oldID || nc == 0) {
289           oldID = fLossID;
290           break;
291       }
292   }
293   SetHighWaterMark(nt);
294 }
295  
296
297 Float_t AliGenHalo::GasPressureWeight(Float_t zPrimary)
298 {
299 //
300 // Return z-dependent gasspressure weight = interaction rate [1/m/s].
301 //
302     Float_t weight = 0.;
303     zPrimary /= 100.;        // m
304     if (fRunPeriod < 5) {
305         Float_t zAbs = TMath::Abs(zPrimary);
306         if (zPrimary > 0.) 
307         {
308             if (zAbs > fZ1[20]) {
309                 weight = 2.e4;
310             } else {
311                 for (Int_t i = 1; i < 21; i++) {
312                     if (zAbs < fZ1[i]) {
313                         weight = fG1[i];
314                         break;
315                     }
316                 }
317             }
318         } else {
319             if (zAbs > fZ2[20]) {
320                 weight = 2.e4;
321             } else {
322                 for (Int_t i = 1; i < 21; i++) {
323                     if (zAbs < fZ2[i]) {
324                     weight = fG2[i];
325                     break;
326                     }
327                 }
328             }
329         }
330     } else {
331         Int_t index = TMath::BinarySearch(fGPASize, fZ1, zPrimary);
332         weight = fG1[index];
333     }
334     return weight;
335 }
336
337 void AliGenHalo::Draw(Option_t *)
338 {
339 // Draws the gas pressure distribution
340     Float_t z[400];
341     Float_t p[400];
342     
343     for (Int_t i = 0; i < 400; i++)
344     {
345         z[i] = -20000. + Float_t(i) * 100;
346         p[i] = GasPressureWeight(z[i]);
347     }
348     
349     TGraph*  gr = new TGraph(400, z, p);   
350     TCanvas* c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
351     c1->cd();
352     gr->Draw("AL");
353 }
354
355 Int_t AliGenHalo::ReadNextParticle()
356 {
357     // Read the next particle from the file
358     Int_t ncols = fscanf(fFile,"%d %f %f %d %f %d %f %f %f %f %f %f %f %f",
359                    &fLossID, &fLossT0, &fLossZ, &fLossA, &fLossW, &fPdg, &fXS, &fYS, &fZS, &fDX, &fDY, &fEkin, &fTS, &fWS);
360     fLossZ /= 10.;
361     fXS    /= 10.;
362     fYS    /= 10.; 
363     fZS    /= 10.;   
364     fTS    *= 1.e-9;
365     return (ncols);
366 }
367
368 void AliGenHalo::SkipEvents()
369 {
370   //
371   // Skip the first fNskip events
372   Int_t skip = fNskip;
373   Int_t oldID = -1;
374
375   while (skip >= 0)
376     {
377       ReadNextParticle();
378       if (oldID != fLossID) {
379         oldID = fLossID;
380         skip--;
381       }
382     } 
383 }
384 void AliGenHalo::CountEvents()
385 {
386     Int_t nev = 0;
387     Int_t oldID = -1;
388     Int_t nc = 1;
389     while (nc != -1)
390     {
391         nc = ReadNextParticle();
392         if (oldID != fLossID) {
393             oldID = fLossID;
394             nev++;
395             printf("Number of events %10d %10d \n", nev, oldID);
396         }
397     }
398
399
400     rewind(fFile);
401 }
402
403