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