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