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