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