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.5  2002/02/01 14:12:28  morsch
  19 Include new gas pressure estimates (LHC Pproject Note 273)
  20
  21 Revision 1.4  2001/12/19 16:30:24  morsch
  22 Some bugs corrected and skip method added. (Rachid Guernane)
  23
  24 Revision 1.3  2001/07/27 17:09:36  morsch
  25 Use local SetTrack, KeepTrack and SetHighWaterMark methods
  26 to delegate either to local stack or to stack owned by AliRun.
  27 (Piotr Skowronski, A.M.)
  28
  29 Revision 1.2  2001/06/14 12:15:27  morsch
  30 Bugs corrected. SetSide() method added.
  31
  32 Revision 1.1  2001/01/23 15:04:33  morsch
  33 Generator to read beam halo file from Protvino group.
  34
  35 */
  36
  37 // Read background particles from a boundary source
  38 // Very specialized generator to simulate background from beam halo.
  39 // The input file is a text file specially prepared
  40 // for this purpose.
  41 // Author: andreas.morsch@cern.ch
  42
  43 #include "AliGenHaloProtvino.h"
  44 #include "AliRun.h"
  45 #include "AliMC.h"
  46 #include "AliPDG.h"
  47
  48 #include <TDatabasePDG.h>
  49 #include <TSystem.h>
  50 #include <stdlib.h>
  51
  52  ClassImp(AliGenHaloProtvino)
  53      AliGenHaloProtvino::AliGenHaloProtvino()
  54          :AliGenerator(-1)
  55 {
  56 // Constructor
  57
  58     fName  = "HaloProtvino";
  59     fTitle = "Halo from LHC Tunnel";
  60 //
  61 //  Read all particles
  62     fNpart = -1;
  63     fFile  =  0;
  64     fSide  =  1;
  65 //
  66     SetRunPeriod();
  67 }
  68
  69 AliGenHaloProtvino::AliGenHaloProtvino(Int_t npart)
  70     :AliGenerator(npart)
  71 {
  72 // Constructor
  73     fName = "Halo";
  74     fTitle= "Halo from LHC Tunnel";
  75 //
  76     fNpart   = npart;
  77     fFile    = 0;
  78     fSide    = 1;
  79 //
  80     SetRunPeriod();
  81 }
  82
  83 AliGenHaloProtvino::AliGenHaloProtvino(const AliGenHaloProtvino & HaloProtvino)
  84 {
  85 // copy constructor
  86 }
  87
  88
  89 //____________________________________________________________
  90 AliGenHaloProtvino::~AliGenHaloProtvino()
  91 {
  92 // Destructor
  93 }
  94
  95 //____________________________________________________________
  96 void AliGenHaloProtvino::Init()
  97 {
  98 // Initialisation
  99     fFile = fopen(fFileName,"r");
 100     if (fFile) {
 101         printf("\n File %s opened for reading, %p ! \n ",  fFileName.Data(), fFile);
 102     } else {
 103         printf("\n Opening of file %s failed,  %p ! \n ",  fFileName.Data(), fFile);
 104     }
 105 //
 106 //
 107 //
 108 //    Read file with gas pressure values
 109     char *name;
 110
 111     name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
 112     FILE* file = fopen(name, "r");
 113     Float_t z;
 114     Int_t i, j;
 115
 116 //
 117 //  Ring 1
 118 //
 119     for (i = 0; i < 20; i++)
 120     {
 121         fscanf(file, "%f %f %f %f %f %f", &z,
 122                &fG1[i][0], &fG1[i][1], &fG1[i][2], &fG1[i][3], &fG1[i][4]);
 123         if (i > 0) {
 124             fZ1[i] = fZ1[i-1] + z;
 125         } else {
 126             fZ1[i] = 20.;
 127         }
 128   }
 129 //
 130 // Ring 2
 131 //
 132     for (i = 0; i < 21; i++)
 133     {
 134         fscanf(file, "%f %f %f %f %f %f", &z,
 135                &fG2[i][0], &fG2[i][1], &fG2[i][2], &fG2[i][3], &fG2[i][4]);
 136         if (i > 0) {
 137             fZ2[i] = fZ2[i-1] + z;
 138         } else {
 139             fZ2[i] = 20.;
 140         }
 141     }
 142 //
 143 //  Transform into interaction rates
 144 //
 145     const Float_t crossSection = 0.094e-28;     // m^2
 146     const Float_t pFlux        = 1.e11/25.e-9;  // 1/s
 147
 148     for (j = 0; j <  5; j++) {
 149         for (i = 0; i < 21; i++)
 150         {
 151             if (i <20)
 152                 fG1[i][j] = fG1[i][j] * crossSection * pFlux; // 1/m/s
 153                 fG2[i][j] = fG2[i][j] * crossSection * pFlux; // 1/m/s
 154         }
 155     }
 156
 157
 158     Float_t sum;
 159
 160     for (Int_t i = 0; i < 250; i++) {
 161         Float_t z = 20.+i*1.;
 162         z*=100;
 163         Float_t wgt = GassPressureWeight(z);
 164 //      printf("\n %f %f", z, wgt);
 165         sum+=wgt;
 166     }
 167     sum/=250.;
 168     printf("\n %f \n \n", sum);
 169 }
 170
 171 //____________________________________________________________
 172 void AliGenHaloProtvino::Generate()
 173 {
 174 // Generate from input file
 175
 176   Float_t polar[3]= {0,0,0};
 177   Float_t origin[3];
 178   Float_t p[3], p0;
 179   Float_t ekin, wgt, tx, ty, tz, txy;
 180   Float_t zPrimary;
 181   Float_t amass;
 182   Int_t   inuc;
 183   //
 184   Int_t ipart, ncols, nt;
 185   Int_t nskip = 0;
 186   Int_t nread = 0;
 187   while(1) {
 188       ncols = fscanf(fFile,"%f %d %d %f %f %f %f %f %f",
 189                      &zPrimary, &inuc, &ipart, &wgt,
 190                      &ekin, &origin[0], &origin[1],
 191                      &tx, &ty);
 192 /*
 193       printf(" \n %f %d %d %f %f %f %f %f %f",
 194              zPrimary, inuc, ipart, wgt,
 195              ekin, origin[0], origin[1],
 196              tx, ty);
 197 */
 198
 199       if (ncols < 0) break;
 200
 201       nskip++;
 202       if (fNpart !=-1 && nskip <= fNskip) continue;
 203
 204       nread++;
 205       if (fNpart !=-1 && nread > fNpart) break;
 206
 207       amass = TDatabasePDG::Instance()->GetParticle(ipart)->Mass();
 208
 209       //
 210       // Momentum vector
 211       //
 212       p0=sqrt(ekin*ekin + 2.*amass*ekin);
 213
 214       txy=TMath::Sqrt(tx*tx+ty*ty);
 215       if (txy == 1.) {
 216           tz=0;
 217       } else {
 218           tz=-TMath::Sqrt(1.-txy);
 219       }
 220
 221       p[0]=p0*tx;
 222       p[1]=p0*ty;
 223       p[2]=-p0*tz;
 224
 225       origin[2] = -2196.5;
 226
 227       //
 228       //
 229       // Particle weight
 230
 231       Float_t originP[3] = {0., 0., 0.};
 232       originP[2] = zPrimary;
 233
 234       Float_t pP[3] = {0., 0., 0.};
 235       Int_t ntP;
 236
 237       if (fSide == -1) {
 238           originP[2] = -zPrimary;
 239           origin[2]  = -origin[2];
 240           p[2]       = -p[2];
 241       }
 242
 243       SetTrack(0,-1,kProton,pP,originP,polar,0,kPNoProcess,ntP);
 244       KeepTrack(ntP);
 245       fParentWeight=wgt*GassPressureWeight(zPrimary);
 246       SetTrack(fTrackIt,ntP,ipart,p,origin,polar,0,kPNoProcess,nt,fParentWeight);
 247       SetHighWaterMark(nt);
 248
 249       //
 250       // Assume particles come from two directions with same probability
 251
 252       // Both Side are considered
 253       if (fSide > 1) {
 254           origin[2]=-origin[2];
 255           p[2]=-p[2];
 256           fParentWeight=wgt*GassPressureWeight(-zPrimary);
 257           SetTrack(fTrackIt,ntP,ipart,p,origin,polar,0,kPNoProcess,nt,fParentWeight);
 258       }
 259
 260   }
 261 }
 262
 263
 264 AliGenHaloProtvino& AliGenHaloProtvino::operator=(const  AliGenHaloProtvino& rhs)
 265 {
 266 // Assignment operator
 267     return *this;
 268 }
 269
 270
 271
 272 Float_t AliGenHaloProtvino::GassPressureWeight(Float_t zPrimary)
 273 {
 274 //
 275 // Return z-dependent gasspressure weight = interaction rate [1/m/s].
 276 //
 277     Float_t weight = 0.;
 278     zPrimary/=100.;        // m
 279     Float_t zAbs = TMath::Abs(zPrimary);
 280     zPrimary = TMath::Abs(zPrimary);
 281
 282     if (zPrimary < 0.)
 283     {
 284         if (zAbs > fZ1[19]) {
 285             weight = 2.e4;
 286         } else {
 287             for (Int_t i = 1; i < 20; i++) {
 288                 if (zAbs < fZ1[i]) {
 289                     weight = fG1[i][fRunPeriod];
 290                     break;
 291                 }
 292             }
 293         }
 294     } else {
 295         if (zAbs > fZ2[20]) {
 296             weight = 2.e4;
 297         } else {
 298             for (Int_t i = 1; i < 21; i++) {
 299                 if (zAbs < fZ2[i]) {
 300                     weight = fG2[i][fRunPeriod];
 301                     break;
 302                 }
 303             }
 304         }
 305     }
 306
 307     return weight;
 308 }
 309
 310 /*
 311 # Title:    README file for the sources of IR8 machine induced background
 312 # Author:   Vadim Talanov <Vadim.Talanov@cern.ch>
 313 # Modified: 12-12-2000
 314
 315 0. Overview
 316
 317         There are three files, named ring.one.beta.[01,10,50].m, which
 318         contain the lists of background particles, induced by proton losses
 319         upstream of IP8 in the LHC ring one, for the beta* values of 1, 10
 320         and 50 m, respectively.
 321
 322 1. File contents
 323
 324         Each line in the files contains the coordinates of particle track
 325         crossing with the infinite plane, positioned at z=-1m, together with
 326         the physical properties of corresponding particle, namely:
 327
 328         S  - S coordinate of the primary interaction vertex, cm;
 329         N  - type of the gas nuclei at interaction, 1 is H, 2 - C and 3 - O;
 330         I  - particle ID in PDG particle numbering scheme;
 331         W  - particle weight;
 332         E  - particle kinetic energy, GeV;
 333         X  - x coordinate of the crossing point, cm;
 334         Y  - y coordinate of the crossing point, cm;
 335         Dx - x direction cosine;
 336         Dy - y direction cosine.
 337
 338 2. Normalisation
 339
 340         Each file is given per unity of linear density of proton inelastic
 341         interactions with the gas nuclei, [1 inelastic interaction/m].
 342
 343 # ~/vtalanov/public/README.mib: the end.
 344
 345 */
 346
 347
 348
 349