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