Add the number of local boards
[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 "AliGenHaloProtvino.h"
33 #include "AliRun.h"
34
35 ClassImp(AliGenHaloProtvino)
36
37 AliGenHaloProtvino::AliGenHaloProtvino()
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 {
51 // Constructor
52     
53     fName  = "HaloProtvino";
54     fTitle = "Halo from LHC Tunnel";
55 //
56 //  Read all particles
57     fNpart = -1;
58     SetAnalog(0);
59 }
60
61 AliGenHaloProtvino::AliGenHaloProtvino(Int_t npart)
62     :AliGenerator(npart),
63      fFile(0),
64      fFileName(0),
65      fSide(1),
66      fRunPeriod(kY3D90),
67      fTimePerEvent(1.e-4),
68      fNskip(0),
69      fZ1(0),
70      fZ2(0),
71      fG1(0),
72      fG2(0),
73      fGPASize(0)
74 {
75 // Constructor
76     fName = "Halo";
77     fTitle= "Halo from LHC Tunnel";
78 //
79     fNpart   = npart;
80 //
81     SetAnalog(0);
82 }
83
84 //____________________________________________________________
85 AliGenHaloProtvino::~AliGenHaloProtvino()
86 {
87 // Destructor
88 }
89
90 //____________________________________________________________
91 void AliGenHaloProtvino::Init() 
92 {
93 // Initialisation
94     fFile = fopen(fFileName,"r");
95     if (fFile) {
96         printf("\n File %s opened for reading, %p ! \n ",  fFileName.Data(), (void*)fFile);
97     } else {
98         printf("\n Opening of file %s failed,  %p ! \n ",  fFileName.Data(), (void*)fFile);
99     }
100 //
101 //
102 //
103 //    Read file with gas pressure values
104     char *name = 0;
105     if (fRunPeriod < 5) {
106         name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
107         fGPASize = 21;
108         fG1 = new Float_t[fGPASize];
109         fG2 = new Float_t[fGPASize];
110         fZ1 = new Float_t[fGPASize];
111         fZ2 = new Float_t[fGPASize];
112     } else if (fRunPeriod == 5) {
113         name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_startup.dat");
114         fGPASize = 18853;
115         fG1 = new Float_t[fGPASize];
116         fZ1 = new Float_t[fGPASize];
117     } else if (fRunPeriod ==6) {
118         name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_conditioned.dat");
119         fGPASize = 12719;
120         fG1 = new Float_t[fGPASize];
121         fZ1 = new Float_t[fGPASize];
122     } else {
123         Fatal("Init()", "No gas pressure file for given run period !");
124     }
125
126
127     FILE* file = fopen(name, "r");
128     Float_t z;
129     Int_t i;
130     Float_t p[5];    
131
132     const Float_t kCrossSection = 0.094e-28;      // m^2
133     const Float_t kFlux         = 1.e11 / 25.e-9; // protons/s
134     Float_t pFlux[5] = {0.2, 0.2, 0.3, 0.3, 1.0};
135
136     if (fRunPeriod < 5) {
137 //
138 //  Ring 1   
139 // 
140
141         for (i = 0; i < fGPASize; i++)
142         {
143             fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
144             fG1[i] = p[fRunPeriod];
145             
146             if (i > 0) {
147                 fZ1[i] = fZ1[i-1] + z;
148             } else {
149                 fZ1[i] = 20.;
150             }
151         }
152 //
153 // Ring 2
154 //
155         for (i = 0; i < fGPASize; i++)
156         {
157             fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
158             fG2[i] = p[fRunPeriod];
159             if (i > 0) {
160                 fZ2[i] = fZ2[i-1] + z;
161             } else {
162                 fZ2[i] = 20.;
163             }
164         }
165 //
166 // Interaction rates
167 //
168         for (i = 0; i < fGPASize; i++)  
169         {
170             fG1[i] = fG1[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s 
171             fG2[i] = fG2[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
172         }
173
174     } else {
175         for (i = 0; i < fGPASize; i++) 
176         {
177             fscanf(file, "%f %e %e %e %e %e", &z, &p[0], &p[1], &p[2], &p[3], &p[4]);
178             z /= 1000.;
179             fG1[i] = p[4] * kCrossSection * kFlux;             // 1/m/s
180             // 1/3 of nominal intensity at startup
181             if (fRunPeriod ==  kLHCPR674Startup) fG1[i] /= 3.;
182             fZ1[i] = z;
183         }
184     }
185     
186
187
188     
189 //
190 //  Transform into interaction rates
191 //
192
193
194     
195
196     Float_t sum1 = 0.;
197     Float_t sum2 = 0.;
198     
199     for (Int_t i = 0; i < 300; i++) {
200         Float_t z = 20.+i*1.;
201         z*=100;
202         Float_t wgt1 = GasPressureWeight(z);
203         Float_t wgt2 = GasPressureWeight(-z);
204 //      printf("weight: %f %f %f %f %f \n", z, wgt1, wgt2, fZ1[20], fZ2[20]);
205         sum1 += wgt1;
206         sum2 += wgt2;
207     }
208     sum1/=250.;
209     sum2/=250.;
210     printf("\n %f %f \n \n", sum1, sum2);
211 }
212
213 //____________________________________________________________
214 void AliGenHaloProtvino::Generate()
215 {
216 // Generate from input file
217  
218   Float_t polar[3]= {0,0,0};
219   Float_t origin[3];
220   Float_t p[3], p0;
221   Float_t tz, txy;
222   Float_t amass;
223   //
224   Int_t ncols, nt;
225   static Int_t nskip = 0;
226   Int_t nread = 0;
227
228   Float_t* zPrimary = new Float_t [fNpart];
229   Int_t  * inuc     = new Int_t   [fNpart];
230   Int_t  * ipart    = new Int_t   [fNpart];
231   Float_t* wgt      = new Float_t [fNpart]; 
232   Float_t* ekin     = new Float_t [fNpart];
233   Float_t* vx       = new Float_t [fNpart];
234   Float_t* vy       = new Float_t [fNpart];
235   Float_t* tx       = new Float_t [fNpart];
236   Float_t* ty       = new Float_t [fNpart];
237   
238   Float_t zVertexOld = -1.e10;
239   Int_t   nInt       = 0;        // Counts number of interactions
240   Float_t wwgt = 0.;
241   
242   while(1) {
243 //
244 // Load event into array
245 //
246       ncols = fscanf(fFile,"%f %d %d %f %f %f %f %f %f",
247                      &zPrimary[nread], &inuc[nread], &ipart[nread], &wgt[nread], 
248                      &ekin[nread], &vx[nread], &vy[nread],
249                      &tx[nread], &ty[nread]);
250       
251       if (ncols < 0) break;
252 // Skip fNskip events
253       nskip++;
254       if (fNpart !=-1 && nskip <= fNskip) continue;
255 // Count interactions
256       if (zPrimary[nread] != zVertexOld) {
257           nInt++;
258           zVertexOld = zPrimary[nread];
259       }
260 // Count tracks      
261       nread++;
262       if (fNpart !=-1 && nread > fNpart) break;
263   }
264 //
265 // Mean time between interactions
266 //
267   Float_t dT = fTimePerEvent/nInt;   // sec 
268   Float_t t  = 0;                    // sec
269   
270 //
271 // Loop over primaries
272 //
273   zVertexOld   = -1.e10;
274   Double_t arg = 0.;
275   
276   for (Int_t nprim = 0; nprim < fNpart; nprim++) 
277   {
278       amass = TDatabasePDG::Instance()->GetParticle(ipart[nprim])->Mass();
279
280       //
281       // Momentum vector
282       //
283       p0=sqrt(ekin[nprim]*ekin[nprim] + 2.*amass*ekin[nprim]);
284       
285       txy=TMath::Sqrt(tx[nprim]*tx[nprim]+ty[nprim]*ty[nprim]);
286       if (txy == 1.) {
287           tz=0;
288       } else {
289           tz=-TMath::Sqrt(1.-txy);
290       }
291     
292       p[0] = p0*tx[nprim];
293       p[1] = p0*ty[nprim];
294       p[2] =-p0*tz;
295       
296       origin[0] = vx[nprim];
297       origin[1] = vy[nprim];
298       origin[2] = -2196.5;
299
300       //
301       //
302       // Particle weight
303
304       Float_t originP[3] = {0., 0., 0.};
305       originP[2] = zPrimary[nprim];
306       
307       Float_t pP[3] = {0., 0., 0.};
308       Int_t ntP;
309       
310       if (fSide == -1) {
311           originP[2] = -zPrimary[nprim];
312           origin[2]  = -origin[2];
313           p[2]       = -p[2];
314       }
315
316       //
317       // Time
318       //
319       if (zPrimary[nprim] != zVertexOld) {
320           while(arg==0.) arg = gRandom->Rndm();
321           t -= dT*TMath::Log(arg);              // (sec)   
322           zVertexOld = zPrimary[nprim];
323       }
324
325 //    Get statistical weight according to local gas-pressure
326 //
327       fParentWeight=wgt[nprim]*GasPressureWeight(zPrimary[nprim]);
328
329       if (!fAnalog || gRandom->Rndm() < fParentWeight) {
330 //    Pass parent particle
331 //
332           PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
333           KeepTrack(ntP);
334           PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
335       }
336       //
337       // Both sides are considered
338       //
339
340       if (fSide > 1) {
341           fParentWeight=wgt[nprim]*GasPressureWeight(-zPrimary[nprim]);
342           if (!fAnalog || gRandom->Rndm() < fParentWeight) {
343               origin[2]  = -origin[2];
344               originP[2] = -originP[2];
345               p[2]=-p[2];
346               PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
347               KeepTrack(ntP);
348               PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
349           }
350       }
351       wwgt += fParentWeight;
352       
353       SetHighWaterMark(nt);
354   }
355   delete [] zPrimary;
356   delete [] inuc;    
357   delete [] ipart;   
358   delete [] wgt;     
359   delete [] ekin;    
360   delete [] vx;      
361   delete [] vy;      
362   delete [] tx;      
363   delete [] ty;     
364   printf("Total weight %f\n\n", wwgt);
365   
366 }
367  
368
369 Float_t AliGenHaloProtvino::GasPressureWeight(Float_t zPrimary)
370 {
371 //
372 // Return z-dependent gasspressure weight = interaction rate [1/m/s].
373 //
374     Float_t weight = 0.;
375     zPrimary /= 100.;        // m
376     if (fRunPeriod < 5) {
377         Float_t zAbs = TMath::Abs(zPrimary);
378         if (zPrimary > 0.) 
379         {
380             if (zAbs > fZ1[20]) {
381                 weight = 2.e4;
382             } else {
383                 for (Int_t i = 1; i < 21; i++) {
384                     if (zAbs < fZ1[i]) {
385                         weight = fG1[i];
386                         break;
387                     }
388                 }
389             }
390         } else {
391             if (zAbs > fZ2[20]) {
392                 weight = 2.e4;
393             } else {
394                 for (Int_t i = 1; i < 21; i++) {
395                     if (zAbs < fZ2[i]) {
396                     weight = fG2[i];
397                     break;
398                     }
399                 }
400             }
401         }
402     } else {
403         Int_t index = TMath::BinarySearch(fGPASize, fZ1, zPrimary);
404         weight = fG1[index];
405     }
406     return weight;
407 }
408
409 void AliGenHaloProtvino::Draw(Option_t *)
410 {
411 // Draws the gas pressure distribution
412     Float_t z[400];
413     Float_t p[400];
414     
415     for (Int_t i = 0; i < 400; i++)
416     {
417         z[i] = -20000. + Float_t(i) * 100;
418         p[i] = GasPressureWeight(z[i]);
419     }
420     
421     TGraph*  gr = new TGraph(400, z, p);   
422     TCanvas* c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
423     c1->cd();
424     gr->Draw("AL");
425 }
426
427
428 /*
429 # Title:    README file for the sources of IR8 machine induced background
430 # Author:   Vadim Talanov <Vadim.Talanov@cern.ch>
431 # Modified: 12-12-2000 
432
433 0. Overview
434
435         There are three files, named ring.one.beta.[01,10,50].m, which
436         contain the lists of background particles, induced by proton losses
437         upstream of IP8 in the LHC ring one, for the beta* values of 1, 10
438         and 50 m, respectively.
439
440 1. File contents
441
442         Each line in the files contains the coordinates of particle track
443         crossing with the infinite plane, positioned at z=-1m, together with
444         the physical properties of corresponding particle, namely:
445
446         S  - S coordinate of the primary interaction vertex, cm;
447         N  - type of the gas nuclei at interaction, 1 is H, 2 - C and 3 - O;
448         I  - particle ID in PDG particle numbering scheme;
449         W  - particle weight;
450         E  - particle kinetic energy, GeV;
451         X  - x coordinate of the crossing point, cm;
452         Y  - y coordinate of the crossing point, cm;
453         Dx - x direction cosine;
454         Dy - y direction cosine.
455
456 2. Normalisation
457
458         Each file is given per unity of linear density of proton inelastic
459         interactions with the gas nuclei, [1 inelastic interaction/m].
460
461 # ~/vtalanov/public/README.mib: the end.
462
463 */
464
465
466
467