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