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