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