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