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