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