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