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     Int_t ncols = 0;
 143     if (fRunPeriod < 5) {
 144 //
 145 //  Ring 1
 146 //
 147
 148         for (i = 0; i < fGPASize; i++)
 149         {
 150             ncols = fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
 151             if (ncols<0) break;
 152
 153             fG1[i] = p[fRunPeriod];
 154
 155             if (i > 0) {
 156                 fZ1[i] = fZ1[i-1] + z;
 157             } else {
 158                 fZ1[i] = 20.;
 159             }
 160         }
 161 //
 162 // Ring 2
 163 //
 164         for (i = 0; i < fGPASize; i++)
 165         {
 166             ncols = fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
 167             if (ncols<0) break;
 168
 169             fG2[i] = p[fRunPeriod];
 170             if (i > 0) {
 171                 fZ2[i] = fZ2[i-1] + z;
 172             } else {
 173                 fZ2[i] = 20.;
 174             }
 175         }
 176 //
 177 // Interaction rates
 178 //
 179         for (i = 0; i < fGPASize; i++)
 180         {
 181             fG1[i] = fG1[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
 182             fG2[i] = fG2[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
 183         }
 184
 185     } else {
 186         for (i = 0; i < fGPASize; i++)
 187         {
 188             ncols = fscanf(file, "%f %e %e %e %e %e", &z, &p[0], &p[1], &p[2], &p[3], &p[4]);
 189             if (ncols<0) break;
 190
 191             z /= 1000.;
 192             fG1[i] = p[4] * kCrossSection * kFlux;             // 1/m/s
 193             // 1/3 of nominal intensity at startup
 194             if (fRunPeriod ==  kLHCPR674Startup) fG1[i] /= 3.;
 195             fZ1[i] = z;
 196         }
 197     }
 198
 199
 200
 201
 202 //
 203 //  Transform into interaction rates
 204 //
 205
 206
 207
 208
 209     Float_t sum1 = 0.;
 210     Float_t sum2 = 0.;
 211
 212     for (Int_t iz = 0; iz < 300; iz++) {
 213         Float_t zpos = 20. + iz * 1.;
 214         zpos *= 100;
 215         Float_t wgt1 = GasPressureWeight( zpos);
 216         Float_t wgt2 = GasPressureWeight(-zpos);
 217         sum1 += wgt1;
 218         sum2 += wgt2;
 219     }
 220     sum1/=250.;
 221     sum2/=250.;
 222     printf("\n %f %f \n \n", sum1, sum2);
 223     delete file;
 224 }
 225
 226 //____________________________________________________________
 227 void AliGenHaloProtvino::Generate()
 228 {
 229 // Generate from input file
 230
 231   Float_t polar[3]= {0,0,0};
 232   Float_t origin[3];
 233   Float_t p[3], p0;
 234   Float_t tz, txy;
 235   Float_t amass;
 236   //
 237   Int_t ncols, nt;
 238   static Int_t nskip = 0;
 239   Int_t nread = 0;
 240
 241   Float_t* zPrimary = new Float_t [fNpart];
 242   Int_t  * inuc     = new Int_t   [fNpart];
 243   Int_t  * ipart    = new Int_t   [fNpart];
 244   Float_t* wgt      = new Float_t [fNpart];
 245   Float_t* ekin     = new Float_t [fNpart];
 246   Float_t* vx       = new Float_t [fNpart];
 247   Float_t* vy       = new Float_t [fNpart];
 248   Float_t* tx       = new Float_t [fNpart];
 249   Float_t* ty       = new Float_t [fNpart];
 250
 251   Float_t zVertexOld = -1.e10;
 252   Int_t   nInt       = 0;        // Counts number of interactions
 253   Float_t wwgt = 0.;
 254
 255   while(1) {
 256 //
 257 // Load event into array
 258 //
 259       ncols = fscanf(fFile,"%f %d %d %f %f %f %f %f %f",
 260                      &zPrimary[nread], &inuc[nread], &ipart[nread], &wgt[nread],
 261                      &ekin[nread], &vx[nread], &vy[nread],
 262                      &tx[nread], &ty[nread]);
 263
 264       if (ncols < 0) break;
 265 // Skip fNskip events
 266       nskip++;
 267       if (fNpart !=-1 && nskip <= fNskip) continue;
 268 // Count interactions
 269       if (zPrimary[nread] != zVertexOld) {
 270           nInt++;
 271           zVertexOld = zPrimary[nread];
 272       }
 273 // Count tracks
 274       nread++;
 275       if (fNpart !=-1 && nread >= fNpart) break;
 276   }
 277 //
 278 // Mean time between interactions
 279 //
 280
 281   Float_t dT = 0.;   // sec
 282   if (nInt > 0)
 283       dT = fTimePerEvent/nInt;
 284   Float_t t  = 0;                    // sec
 285
 286 //
 287 // Loop over primaries
 288 //
 289   zVertexOld   = -1.e10;
 290   Double_t arg = 0.;
 291
 292   for (Int_t nprim = 0; nprim < fNpart; nprim++)
 293   {
 294       amass = TDatabasePDG::Instance()->GetParticle(ipart[nprim])->Mass();
 295
 296       //
 297       // Momentum vector
 298       //
 299       p0=sqrt(ekin[nprim]*ekin[nprim] + 2.*amass*ekin[nprim]);
 300
 301       txy=TMath::Sqrt(tx[nprim]*tx[nprim]+ty[nprim]*ty[nprim]);
 302       if (txy == 1.) {
 303           tz=0;
 304       } else {
 305           tz=-TMath::Sqrt(1.-txy);
 306       }
 307
 308       p[0] = p0*tx[nprim];
 309       p[1] = p0*ty[nprim];
 310       p[2] =-p0*tz;
 311
 312       origin[0] = vx[nprim];
 313       origin[1] = vy[nprim];
 314       origin[2] = -2196.5;
 315
 316       //
 317       //
 318       // Particle weight
 319
 320       Float_t originP[3] = {0., 0., 0.};
 321       originP[2] = zPrimary[nprim];
 322
 323       Float_t pP[3] = {0., 0., 0.};
 324       Int_t ntP;
 325
 326       if (fSide == -1) {
 327           originP[2] = -zPrimary[nprim];
 328           origin[2]  = -origin[2];
 329           p[2]       = -p[2];
 330       }
 331
 332       //
 333       // Time
 334       //
 335       if (zPrimary[nprim] != zVertexOld) {
 336           while(arg==0.) arg = gRandom->Rndm();
 337           t -= dT*TMath::Log(arg);              // (sec)
 338           zVertexOld = zPrimary[nprim];
 339       }
 340
 341 //    Get statistical weight according to local gas-pressure
 342 //
 343       fParentWeight=wgt[nprim]*GasPressureWeight(zPrimary[nprim]);
 344
 345       if (!fAnalog || gRandom->Rndm() < fParentWeight) {
 346 //    Pass parent particle
 347 //
 348           PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
 349           KeepTrack(ntP);
 350           PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
 351       }
 352       //
 353       // Both sides are considered
 354       //
 355
 356       if (fSide > 1) {
 357           fParentWeight=wgt[nprim]*GasPressureWeight(-zPrimary[nprim]);
 358           if (!fAnalog || gRandom->Rndm() < fParentWeight) {
 359               origin[2]  = -origin[2];
 360               originP[2] = -originP[2];
 361               p[2]=-p[2];
 362               PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
 363               KeepTrack(ntP);
 364               PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
 365           }
 366       }
 367       wwgt += fParentWeight;
 368
 369       SetHighWaterMark(nt);
 370   }
 371   delete [] zPrimary;
 372   delete [] inuc;
 373   delete [] ipart;
 374   delete [] wgt;
 375   delete [] ekin;
 376   delete [] vx;
 377   delete [] vy;
 378   delete [] tx;
 379   delete [] ty;
 380   printf("Total weight %f\n\n", wwgt);
 381
 382 }
 383
 384
 385 Float_t AliGenHaloProtvino::GasPressureWeight(Float_t zPrimary)
 386 {
 387 //
 388 // Return z-dependent gasspressure weight = interaction rate [1/m/s].
 389 //
 390     Float_t weight = 0.;
 391     zPrimary /= 100.;        // m
 392     if (fRunPeriod < 5) {
 393         Float_t zAbs = TMath::Abs(zPrimary);
 394         if (zPrimary > 0.)
 395         {
 396             if (zAbs > fZ1[20]) {
 397                 weight = 2.e4;
 398             } else {
 399                 for (Int_t i = 1; i < 21; i++) {
 400                     if (zAbs < fZ1[i]) {
 401                         weight = fG1[i];
 402                         break;
 403                     }
 404                 }
 405             }
 406         } else {
 407             if (zAbs > fZ2[20]) {
 408                 weight = 2.e4;
 409             } else {
 410                 for (Int_t i = 1; i < 21; i++) {
 411                     if (zAbs < fZ2[i]) {
 412                     weight = fG2[i];
 413                     break;
 414                     }
 415                 }
 416             }
 417         }
 418     } else {
 419         Int_t index = TMath::BinarySearch(fGPASize, fZ1, zPrimary);
 420         weight = fG1[index];
 421     }
 422     return weight;
 423 }
 424
 425 void AliGenHaloProtvino::Draw(Option_t *)
 426 {
 427 // Draws the gas pressure distribution
 428     Float_t z[400];
 429     Float_t p[400];
 430
 431     for (Int_t i = 0; i < 400; i++)
 432     {
 433         z[i] = -20000. + Float_t(i) * 100;
 434         p[i] = GasPressureWeight(z[i]);
 435     }
 436
 437     TGraph*  gr = new TGraph(400, z, p);
 438     TCanvas* c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
 439     c1->cd();
 440     gr->Draw("AL");
 441 }
 442
 443
 444 /*
 445 # Title:    README file for the sources of IR8 machine induced background
 446 # Author:   Vadim Talanov <Vadim.Talanov@cern.ch>
 447 # Modified: 12-12-2000
 448
 449 0. Overview
 450
 451         There are three files, named ring.one.beta.[01,10,50].m, which
 452         contain the lists of background particles, induced by proton losses
 453         upstream of IP8 in the LHC ring one, for the beta* values of 1, 10
 454         and 50 m, respectively.
 455
 456 1. File contents
 457
 458         Each line in the files contains the coordinates of particle track
 459         crossing with the infinite plane, positioned at z=-1m, together with
 460         the physical properties of corresponding particle, namely:
 461
 462         S  - S coordinate of the primary interaction vertex, cm;
 463         N  - type of the gas nuclei at interaction, 1 is H, 2 - C and 3 - O;
 464         I  - particle ID in PDG particle numbering scheme;
 465         W  - particle weight;
 466         E  - particle kinetic energy, GeV;
 467         X  - x coordinate of the crossing point, cm;
 468         Y  - y coordinate of the crossing point, cm;
 469         Dx - x direction cosine;
 470         Dy - y direction cosine.
 471
 472 2. Normalisation
 473
 474         Each file is given per unity of linear density of proton inelastic
 475         interactions with the gas nuclei, [1 inelastic interaction/m].
 476
 477 # ~/vtalanov/public/README.mib: the end.
 478
 479 */
 480
 481
 482
 483