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