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