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