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 /*
  17 $Log$
  18 Revision 1.8.4.1  2002/07/24 08:56:28  alibrary
  19 Updating EVGEN on TVirtulaMC
  20
  21 Revision 1.10  2002/05/28 14:28:12  morsch
  22 Correct LHC nominal current profile (1st year 20%, 2nd year 30%, 3rd year 100 %)
  23
  24 Revision 1.9  2002/05/28 13:49:17  morsch
  25 - Udates for new pressure table
  26 - calculate time
  27 - first provisions for real events.
  28
  29 Revision 1.8  2002/03/22 13:00:25  morsch
  30 Initialize sum to 0. (Jiri Chudoba).
  31
  32 Revision 1.7  2002/02/21 16:09:45  morsch
  33 Move SetHighwaterMark() after last possible SetTrack()
  34
  35 Revision 1.6  2002/02/01 15:28:22  morsch
  36 Fastidious print statements removed.
  37
  38 Revision 1.5  2002/02/01 14:12:28  morsch
  39 Include new gas pressure estimates (LHC Pproject Note 273)
  40
  41 Revision 1.4  2001/12/19 16:30:24  morsch
  42 Some bugs corrected and skip method added. (Rachid Guernane)
  43
  44 Revision 1.3  2001/07/27 17:09:36  morsch
  45 Use local SetTrack, KeepTrack and SetHighWaterMark methods
  46 to delegate either to local stack or to stack owned by AliRun.
  47 (Piotr Skowronski, A.M.)
  48
  49 Revision 1.2  2001/06/14 12:15:27  morsch
  50 Bugs corrected. SetSide() method added.
  51
  52 Revision 1.1  2001/01/23 15:04:33  morsch
  53 Generator to read beam halo file from Protvino group.
  54
  55 */
  56
  57 // Read background particles from a boundary source
  58 // Very specialized generator to simulate background from beam halo.
  59 // The input file is a text file specially prepared
  60 // for this purpose.
  61 // Author: andreas.morsch@cern.ch
  62
  63 #include "AliGenHaloProtvino.h"
  64 #include "AliRun.h"
  65 #include "AliMC.h"
  66 #include "AliPDG.h"
  67
  68 #include <TDatabasePDG.h>
  69 #include <TSystem.h>
  70 #include <stdlib.h>
  71
  72  ClassImp(AliGenHaloProtvino)
  73      AliGenHaloProtvino::AliGenHaloProtvino()
  74          :AliGenerator(-1)
  75 {
  76 // Constructor
  77
  78     fName  = "HaloProtvino";
  79     fTitle = "Halo from LHC Tunnel";
  80 //
  81 //  Read all particles
  82     fNpart = -1;
  83     fFile  =  0;
  84     fSide  =  1;
  85 //
  86     SetRunPeriod();
  87     SetTimePerEvent();
  88     SetAnalog(0);
  89 }
  90
  91 AliGenHaloProtvino::AliGenHaloProtvino(Int_t npart)
  92     :AliGenerator(npart)
  93 {
  94 // Constructor
  95     fName = "Halo";
  96     fTitle= "Halo from LHC Tunnel";
  97 //
  98     fNpart   = npart;
  99     fFile    = 0;
 100     fSide    = 1;
 101 //
 102     SetRunPeriod();
 103     SetTimePerEvent();
 104     SetAnalog(0);
 105 }
 106
 107 AliGenHaloProtvino::AliGenHaloProtvino(const AliGenHaloProtvino & HaloProtvino)
 108 {
 109 // copy constructor
 110 }
 111
 112
 113 //____________________________________________________________
 114 AliGenHaloProtvino::~AliGenHaloProtvino()
 115 {
 116 // Destructor
 117 }
 118
 119 //____________________________________________________________
 120 void AliGenHaloProtvino::Init()
 121 {
 122 // Initialisation
 123     fFile = fopen(fFileName,"r");
 124     if (fFile) {
 125         printf("\n File %s opened for reading, %p ! \n ",  fFileName.Data(), fFile);
 126     } else {
 127         printf("\n Opening of file %s failed,  %p ! \n ",  fFileName.Data(), fFile);
 128     }
 129 //
 130 //
 131 //
 132 //    Read file with gas pressure values
 133     char *name;
 134
 135     name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
 136     FILE* file = fopen(name, "r");
 137     Float_t z;
 138     Int_t i, j;
 139
 140 //
 141 //  Ring 1
 142 //
 143     for (i = 0; i < 21; i++)
 144     {
 145         fscanf(file, "%f %f %f %f %f %f", &z,
 146                &fG1[i][0], &fG1[i][1], &fG1[i][2], &fG1[i][3], &fG1[i][4]);
 147         if (i > 0) {
 148             fZ1[i] = fZ1[i-1] + z;
 149         } else {
 150             fZ1[i] = 20.;
 151         }
 152   }
 153 //
 154 // Ring 2
 155 //
 156     for (i = 0; i < 21; i++)
 157     {
 158         fscanf(file, "%f %f %f %f %f %f", &z,
 159                &fG2[i][0], &fG2[i][1], &fG2[i][2], &fG2[i][3], &fG2[i][4]);
 160         if (i > 0) {
 161             fZ2[i] = fZ2[i-1] + z;
 162         } else {
 163             fZ2[i] = 20.;
 164         }
 165     }
 166 //
 167 //  Transform into interaction rates
 168 //
 169     const Float_t crossSection = 0.094e-28;     // m^2
 170     Float_t pFlux[5] = {0.2, 0.2, 0.3, 0.3, 1.0};
 171
 172     for (j = 0; j <  5; j++) {
 173         pFlux[j] *= 1.e11/25.e-9;
 174         for (i = 0; i < 21; i++)
 175         {
 176             fG1[i][j] = fG1[i][j] * crossSection * pFlux[j]; // 1/m/s
 177             fG2[i][j] = fG2[i][j] * crossSection * pFlux[j]; // 1/m/s
 178         }
 179     }
 180
 181
 182     Float_t sum1 = 0.;
 183     Float_t sum2 = 0.;
 184
 185     for (Int_t i = 0; i < 300; i++) {
 186         Float_t z = 20.+i*1.;
 187         z*=100;
 188         Float_t wgt1 = GassPressureWeight(z);
 189         Float_t wgt2 = GassPressureWeight(-z);
 190 //      printf("weight: %f %f %f %f %f \n", z, wgt1, wgt2, fZ1[20], fZ2[20]);
 191         sum1 += wgt1;
 192         sum2 += wgt2;
 193     }
 194     sum1/=250.;
 195     sum2/=250.;
 196     printf("\n %f %f \n \n", sum1, sum2);
 197 }
 198
 199 //____________________________________________________________
 200 void AliGenHaloProtvino::Generate()
 201 {
 202 // Generate from input file
 203
 204   Float_t polar[3]= {0,0,0};
 205   Float_t origin[3];
 206   Float_t p[3], p0;
 207   Float_t tz, txy;
 208   Float_t amass;
 209   //
 210   Int_t ncols, nt;
 211   Int_t nskip = 0;
 212   Int_t nread = 0;
 213
 214   Float_t* zPrimary = new Float_t [fNpart];
 215   Int_t  * inuc     = new Int_t   [fNpart];
 216   Int_t  * ipart    = new Int_t   [fNpart];
 217   Float_t* wgt      = new Float_t [fNpart];
 218   Float_t* ekin     = new Float_t [fNpart];
 219   Float_t* vx       = new Float_t [fNpart];
 220   Float_t* vy       = new Float_t [fNpart];
 221   Float_t* tx       = new Float_t [fNpart];
 222   Float_t* ty       = new Float_t [fNpart];
 223
 224   Float_t zVertexOld = -1.e10;
 225   Int_t   nInt       = 0;        // Counts number of interactions
 226
 227   while(1) {
 228 //
 229 // Load event into array
 230 //
 231       ncols = fscanf(fFile,"%f %d %d %f %f %f %f %f %f",
 232                      &zPrimary[nread], &inuc[nread], &ipart[nread], &wgt[nread],
 233                      &ekin[nread], &vx[nread], &vy[nread],
 234                      &tx[nread], &ty[nread]);
 235
 236       if (ncols < 0) break;
 237 // Skip fNskip events
 238       nskip++;
 239       if (fNpart !=-1 && nskip <= fNskip) continue;
 240 // Count interactions
 241       if (zPrimary[nread] != zVertexOld) {
 242           nInt++;
 243           zVertexOld = zPrimary[nread];
 244       }
 245 // Count tracks
 246       nread++;
 247       if (fNpart !=-1 && nread > fNpart) break;
 248   }
 249 //
 250 // Mean time between interactions
 251 //
 252   Float_t dT = fTimePerEvent/nInt;   // sec
 253   Float_t t  = 0;                    // sec
 254
 255 //
 256 // Loop over primaries
 257 //
 258   zVertexOld   = -1.e10;
 259   Double_t arg = 0.;
 260
 261   for (Int_t nprim = 0; nprim < fNpart; nprim++)
 262   {
 263       amass = TDatabasePDG::Instance()->GetParticle(ipart[nprim])->Mass();
 264
 265       //
 266       // Momentum vector
 267       //
 268       p0=sqrt(ekin[nprim]*ekin[nprim] + 2.*amass*ekin[nprim]);
 269
 270       txy=TMath::Sqrt(tx[nprim]*tx[nprim]+ty[nprim]*ty[nprim]);
 271       if (txy == 1.) {
 272           tz=0;
 273       } else {
 274           tz=-TMath::Sqrt(1.-txy);
 275       }
 276
 277       p[0] = p0*tx[nprim];
 278       p[1] = p0*ty[nprim];
 279       p[2] =-p0*tz;
 280
 281       origin[0] = vx[nprim];
 282       origin[1] = vy[nprim];
 283       origin[2] = -2196.5;
 284
 285       //
 286       //
 287       // Particle weight
 288
 289       Float_t originP[3] = {0., 0., 0.};
 290       originP[2] = zPrimary[nprim];
 291
 292       Float_t pP[3] = {0., 0., 0.};
 293       Int_t ntP;
 294
 295       if (fSide == -1) {
 296           originP[2] = -zPrimary[nprim];
 297           origin[2]  = -origin[2];
 298           p[2]       = -p[2];
 299       }
 300
 301       //
 302       // Time
 303       //
 304       if (zPrimary[nprim] != zVertexOld) {
 305           while(arg==0.) arg = gRandom->Rndm();
 306           t -= dT*TMath::Log(arg);              // (sec)
 307           zVertexOld = zPrimary[nprim];
 308       }
 309
 310 //    Get statistical weight according to local gas-pressure
 311 //
 312       fParentWeight=wgt[nprim]*GassPressureWeight(zPrimary[nprim]);
 313
 314       if (!fAnalog || gRandom->Rndm() < fParentWeight) {
 315 //    Pass parent particle
 316 //
 317           SetTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
 318           KeepTrack(ntP);
 319           SetTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
 320       }
 321
 322       //
 323       // Both sides are considered
 324       //
 325
 326       if (fSide > 1) {
 327           fParentWeight=wgt[nprim]*GassPressureWeight(-zPrimary[nprim]);
 328           if (!fAnalog || gRandom->Rndm() < fParentWeight) {
 329               origin[2]  = -origin[2];
 330               originP[2] = -originP[2];
 331               p[2]=-p[2];
 332               SetTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
 333               KeepTrack(ntP);
 334               SetTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
 335           }
 336       }
 337       SetHighWaterMark(nt);
 338   }
 339   delete [] zPrimary;
 340   delete [] inuc;
 341   delete [] ipart;
 342   delete [] wgt;
 343   delete [] ekin;
 344   delete [] vx;
 345   delete [] vy;
 346   delete [] tx;
 347   delete [] ty;
 348 }
 349
 350
 351 AliGenHaloProtvino& AliGenHaloProtvino::operator=(const  AliGenHaloProtvino& rhs)
 352 {
 353 // Assignment operator
 354     return *this;
 355 }
 356
 357
 358
 359 Float_t AliGenHaloProtvino::GassPressureWeight(Float_t zPrimary)
 360 {
 361 //
 362 // Return z-dependent gasspressure weight = interaction rate [1/m/s].
 363 //
 364     Float_t weight = 0.;
 365     zPrimary/=100.;        // m
 366     Float_t zAbs = TMath::Abs(zPrimary);
 367     if (zPrimary > 0.)
 368     {
 369         if (zAbs > fZ1[20]) {
 370             weight = 2.e4;
 371         } else {
 372             for (Int_t i = 1; i < 21; i++) {
 373                 if (zAbs < fZ1[i]) {
 374                     weight = fG1[i][fRunPeriod];
 375                     break;
 376                 }
 377             }
 378         }
 379     } else {
 380         if (zAbs > fZ2[20]) {
 381             weight = 2.e4;
 382         } else {
 383             for (Int_t i = 1; i < 21; i++) {
 384                 if (zAbs < fZ2[i]) {
 385                     weight = fG2[i][fRunPeriod];
 386                     break;
 387                 }
 388             }
 389         }
 390     }
 391     return weight;
 392 }
 393
 394 /*
 395 # Title:    README file for the sources of IR8 machine induced background
 396 # Author:   Vadim Talanov <Vadim.Talanov@cern.ch>
 397 # Modified: 12-12-2000
 398
 399 0. Overview
 400
 401         There are three files, named ring.one.beta.[01,10,50].m, which
 402         contain the lists of background particles, induced by proton losses
 403         upstream of IP8 in the LHC ring one, for the beta* values of 1, 10
 404         and 50 m, respectively.
 405
 406 1. File contents
 407
 408         Each line in the files contains the coordinates of particle track
 409         crossing with the infinite plane, positioned at z=-1m, together with
 410         the physical properties of corresponding particle, namely:
 411
 412         S  - S coordinate of the primary interaction vertex, cm;
 413         N  - type of the gas nuclei at interaction, 1 is H, 2 - C and 3 - O;
 414         I  - particle ID in PDG particle numbering scheme;
 415         W  - particle weight;
 416         E  - particle kinetic energy, GeV;
 417         X  - x coordinate of the crossing point, cm;
 418         Y  - y coordinate of the crossing point, cm;
 419         Dx - x direction cosine;
 420         Dy - y direction cosine.
 421
 422 2. Normalisation
 423
 424         Each file is given per unity of linear density of proton inelastic
 425         interactions with the gas nuclei, [1 inelastic interaction/m].
 426
 427 # ~/vtalanov/public/README.mib: the end.
 428
 429 */
 430
 431
 432
 433