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