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