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   Float_t dT = fTimePerEvent/nInt;   // sec
 274   Float_t t  = 0;                    // sec
 275
 276 //
 277 // Loop over primaries
 278 //
 279   zVertexOld   = -1.e10;
 280   Double_t arg = 0.;
 281
 282   for (Int_t nprim = 0; nprim < fNpart; nprim++)
 283   {
 284       amass = TDatabasePDG::Instance()->GetParticle(ipart[nprim])->Mass();
 285
 286       //
 287       // Momentum vector
 288       //
 289       p0=sqrt(ekin[nprim]*ekin[nprim] + 2.*amass*ekin[nprim]);
 290
 291       txy=TMath::Sqrt(tx[nprim]*tx[nprim]+ty[nprim]*ty[nprim]);
 292       if (txy == 1.) {
 293           tz=0;
 294       } else {
 295           tz=-TMath::Sqrt(1.-txy);
 296       }
 297
 298       p[0] = p0*tx[nprim];
 299       p[1] = p0*ty[nprim];
 300       p[2] =-p0*tz;
 301
 302       origin[0] = vx[nprim];
 303       origin[1] = vy[nprim];
 304       origin[2] = -2196.5;
 305
 306       //
 307       //
 308       // Particle weight
 309
 310       Float_t originP[3] = {0., 0., 0.};
 311       originP[2] = zPrimary[nprim];
 312
 313       Float_t pP[3] = {0., 0., 0.};
 314       Int_t ntP;
 315
 316       if (fSide == -1) {
 317           originP[2] = -zPrimary[nprim];
 318           origin[2]  = -origin[2];
 319           p[2]       = -p[2];
 320       }
 321
 322       //
 323       // Time
 324       //
 325       if (zPrimary[nprim] != zVertexOld) {
 326           while(arg==0.) arg = gRandom->Rndm();
 327           t -= dT*TMath::Log(arg);              // (sec)
 328           zVertexOld = zPrimary[nprim];
 329       }
 330
 331 //    Get statistical weight according to local gas-pressure
 332 //
 333       fParentWeight=wgt[nprim]*GasPressureWeight(zPrimary[nprim]);
 334
 335       if (!fAnalog || gRandom->Rndm() < fParentWeight) {
 336 //    Pass parent particle
 337 //
 338           PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
 339           KeepTrack(ntP);
 340           PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
 341       }
 342       //
 343       // Both sides are considered
 344       //
 345
 346       if (fSide > 1) {
 347           fParentWeight=wgt[nprim]*GasPressureWeight(-zPrimary[nprim]);
 348           if (!fAnalog || gRandom->Rndm() < fParentWeight) {
 349               origin[2]  = -origin[2];
 350               originP[2] = -originP[2];
 351               p[2]=-p[2];
 352               PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
 353               KeepTrack(ntP);
 354               PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
 355           }
 356       }
 357       wwgt += fParentWeight;
 358
 359       SetHighWaterMark(nt);
 360   }
 361   delete [] zPrimary;
 362   delete [] inuc;
 363   delete [] ipart;
 364   delete [] wgt;
 365   delete [] ekin;
 366   delete [] vx;
 367   delete [] vy;
 368   delete [] tx;
 369   delete [] ty;
 370   printf("Total weight %f\n\n", wwgt);
 371
 372 }
 373
 374
 375 Float_t AliGenHaloProtvino::GasPressureWeight(Float_t zPrimary)
 376 {
 377 //
 378 // Return z-dependent gasspressure weight = interaction rate [1/m/s].
 379 //
 380     Float_t weight = 0.;
 381     zPrimary /= 100.;        // m
 382     if (fRunPeriod < 5) {
 383         Float_t zAbs = TMath::Abs(zPrimary);
 384         if (zPrimary > 0.)
 385         {
 386             if (zAbs > fZ1[20]) {
 387                 weight = 2.e4;
 388             } else {
 389                 for (Int_t i = 1; i < 21; i++) {
 390                     if (zAbs < fZ1[i]) {
 391                         weight = fG1[i];
 392                         break;
 393                     }
 394                 }
 395             }
 396         } else {
 397             if (zAbs > fZ2[20]) {
 398                 weight = 2.e4;
 399             } else {
 400                 for (Int_t i = 1; i < 21; i++) {
 401                     if (zAbs < fZ2[i]) {
 402                     weight = fG2[i];
 403                     break;
 404                     }
 405                 }
 406             }
 407         }
 408     } else {
 409         Int_t index = TMath::BinarySearch(fGPASize, fZ1, zPrimary);
 410         weight = fG1[index];
 411     }
 412     return weight;
 413 }
 414
 415 void AliGenHaloProtvino::Draw(Option_t *)
 416 {
 417 // Draws the gas pressure distribution
 418     Float_t z[400];
 419     Float_t p[400];
 420
 421     for (Int_t i = 0; i < 400; i++)
 422     {
 423         z[i] = -20000. + Float_t(i) * 100;
 424         p[i] = GasPressureWeight(z[i]);
 425     }
 426
 427     TGraph*  gr = new TGraph(400, z, p);
 428     TCanvas* c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
 429     c1->cd();
 430     gr->Draw("AL");
 431 }
 432
 433
 434 /*
 435 # Title:    README file for the sources of IR8 machine induced background
 436 # Author:   Vadim Talanov <Vadim.Talanov@cern.ch>
 437 # Modified: 12-12-2000
 438
 439 0. Overview
 440
 441         There are three files, named ring.one.beta.[01,10,50].m, which
 442         contain the lists of background particles, induced by proton losses
 443         upstream of IP8 in the LHC ring one, for the beta* values of 1, 10
 444         and 50 m, respectively.
 445
 446 1. File contents
 447
 448         Each line in the files contains the coordinates of particle track
 449         crossing with the infinite plane, positioned at z=-1m, together with
 450         the physical properties of corresponding particle, namely:
 451
 452         S  - S coordinate of the primary interaction vertex, cm;
 453         N  - type of the gas nuclei at interaction, 1 is H, 2 - C and 3 - O;
 454         I  - particle ID in PDG particle numbering scheme;
 455         W  - particle weight;
 456         E  - particle kinetic energy, GeV;
 457         X  - x coordinate of the crossing point, cm;
 458         Y  - y coordinate of the crossing point, cm;
 459         Dx - x direction cosine;
 460         Dy - y direction cosine.
 461
 462 2. Normalisation
 463
 464         Each file is given per unity of linear density of proton inelastic
 465         interactions with the gas nuclei, [1 inelastic interaction/m].
 466
 467 # ~/vtalanov/public/README.mib: the end.
 468
 469 */
 470
 471
 472
 473