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