EVGEN/AliGenHalo.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 beam halo background particles from a boundary source
  19 // Boundary source is in the LHCb format http://www.hep.manchester.ac.uk/u/robert/LHC_backgrounds/Note-MIBStudies.pdf
  20 // and has been provided by Robert Appleby
  21 // Author: andreas.morsch@cern.ch
  22
  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 "AliGenHalo.h"
  33 #include "AliRun.h"
  34
  35 ClassImp(AliGenHalo)
  36
  37 AliGenHalo::AliGenHalo()
  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      fLossID(0),
  51      fLossA(0),
  52      fPdg(0),
  53      fLossT0(0),
  54      fLossZ(0),
  55      fLossW(0),
  56      fXS(0),
  57      fYS(0),
  58      fZS(0),
  59      fDX(0),
  60      fDY(0),
  61      fEkin(0),
  62      fTS(0),
  63      fWS(0)
  64 {
  65 // Constructor
  66
  67     fName  = "Halo";
  68     fTitle = "Halo from LHC Beam";
  69 //
  70 //  Read all particles
  71     fNpart = -1;
  72     SetAnalog(0);
  73 }
  74
  75 AliGenHalo::AliGenHalo(Int_t npart)
  76     :AliGenerator(npart),
  77      fFile(0),
  78      fFileName(0),
  79      fSide(1),
  80      fRunPeriod(kY3D90),
  81      fTimePerEvent(1.e-4),
  82      fNskip(0),
  83      fZ1(0),
  84      fZ2(0),
  85      fG1(0),
  86      fG2(0),
  87      fGPASize(0),
  88      fLossID(0),
  89      fLossA(0),
  90      fPdg(0),
  91      fLossT0(0),
  92      fLossZ(0),
  93      fLossW(0),
  94      fXS(0),
  95      fYS(0),
  96      fZS(0),
  97      fDX(0),
  98      fDY(0),
  99      fEkin(0),
 100      fTS(0),
 101      fWS(0)
 102 {
 103 // Constructor
 104     fName = "Halo";
 105     fTitle= "Halo from LHC Beam";
 106 //
 107     fNpart   = npart;
 108 //
 109     SetAnalog(0);
 110 }
 111
 112 //____________________________________________________________
 113 AliGenHalo::~AliGenHalo()
 114 {
 115 // Destructor
 116 }
 117
 118 //____________________________________________________________
 119 void AliGenHalo::Init()
 120 {
 121 // Initialisation
 122     fFile = fopen(fFileName,"r");
 123     if (fFile) {
 124         printf("\n File %s opened for reading, %p ! \n ",  fFileName.Data(), (void*)fFile);
 125     } else {
 126         printf("\n Opening of file %s failed,  %p ! \n ",  fFileName.Data(), (void*)fFile);
 127     }
 128
 129     if (fNskip > 0) {
 130       // Skip the first fNskip events
 131       SkipEvents();
 132     }
 133 //
 134 //
 135 //
 136 //    Read file with gas pressure values
 137     char *name = 0;
 138     if (fRunPeriod < 5) {
 139         name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
 140         fGPASize = 21;
 141         fG1 = new Float_t[fGPASize];
 142         fG2 = new Float_t[fGPASize];
 143         fZ1 = new Float_t[fGPASize];
 144         fZ2 = new Float_t[fGPASize];
 145     } else if (fRunPeriod == 5) {
 146         name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_startup.dat");
 147         fGPASize = 18853;
 148         fG1 = new Float_t[fGPASize];
 149         fZ1 = new Float_t[fGPASize];
 150     } else if (fRunPeriod ==6) {
 151         name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_conditioned.dat");
 152         fGPASize = 12719;
 153         fG1 = new Float_t[fGPASize];
 154         fZ1 = new Float_t[fGPASize];
 155     } else {
 156         Fatal("Init()", "No gas pressure file for given run period !");
 157     }
 158
 159
 160     FILE* file = fopen(name, "r");
 161     Float_t z;
 162     Int_t i;
 163     Float_t p[5];
 164
 165     const Float_t kCrossSection = 0.094e-28;      // m^2
 166     const Float_t kFlux         = 1.e11 / 25.e-9; // protons/s
 167     Float_t pFlux[5] = {0.2, 0.2, 0.3, 0.3, 1.0};
 168
 169     if (fRunPeriod < 5) {
 170 //
 171 //  Ring 1
 172 //
 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             fG1[i] = p[fRunPeriod];
 178
 179             if (i > 0) {
 180                 fZ1[i] = fZ1[i-1] + z;
 181             } else {
 182                 fZ1[i] = 20.;
 183             }
 184         }
 185 //
 186 // Ring 2
 187 //
 188         for (i = 0; i < fGPASize; i++)
 189         {
 190             fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
 191             fG2[i] = p[fRunPeriod];
 192             if (i > 0) {
 193                 fZ2[i] = fZ2[i-1] + z;
 194             } else {
 195                 fZ2[i] = 20.;
 196             }
 197         }
 198 //
 199 // Interaction rates
 200 //
 201         for (i = 0; i < fGPASize; i++)
 202         {
 203             fG1[i] = fG1[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
 204             fG2[i] = fG2[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
 205         }
 206
 207     } else {
 208         for (i = 0; i < fGPASize; i++)
 209         {
 210             fscanf(file, "%f %e %e %e %e %e", &z, &p[0], &p[1], &p[2], &p[3], &p[4]);
 211             z /= 1000.;
 212             fG1[i] = p[4] * kCrossSection * kFlux;             // 1/m/s
 213             // 1/3 of nominal intensity at startup
 214             if (fRunPeriod ==  kLHCPR674Startup) fG1[i] /= 3.;
 215             fZ1[i] = z;
 216         }
 217     }
 218
 219
 220
 221
 222 //
 223 //  Transform into interaction rates
 224 //
 225
 226
 227
 228
 229     Float_t sum1 = 0.;
 230     Float_t sum2 = 0.;
 231
 232     for (Int_t iz = 0; iz < 300; iz++) {
 233         Float_t zpos = 20. + iz * 1.;
 234         zpos *= 100;
 235         Float_t wgt1 = GasPressureWeight( zpos);
 236         Float_t wgt2 = GasPressureWeight(-zpos);
 237         sum1 += wgt1;
 238         sum2 += wgt2;
 239     }
 240     sum1/=250.;
 241     sum2/=250.;
 242     printf("\n %f %f \n \n", sum1, sum2);
 243 }
 244
 245 //____________________________________________________________
 246 void AliGenHalo::Generate()
 247 {
 248 // Generate by reading particles from input file
 249
 250   Float_t polar[3]= {0., 0., 0.};
 251   Float_t origin[3];
 252   Float_t p[3], p0;
 253   Float_t tz, txy;
 254   Float_t mass;
 255   //
 256   Int_t nt;
 257   static Bool_t first = kTRUE;
 258   static Int_t  oldID = -1;
 259 //
 260
 261   if (first && (fNskip == 0)) ReadNextParticle();
 262   first = kFALSE;
 263   oldID = fLossID;
 264
 265   while(1) {
 266       // Push particle to stack
 267       mass = TDatabasePDG::Instance()->GetParticle(fPdg)->Mass();
 268       p0  = TMath::Sqrt(fEkin * fEkin + 2. * mass * fEkin);
 269       txy = TMath::Sqrt(fDX * fDX + fDY * fDY);
 270       if (txy > 1.) {
 271           tz = 0.;
 272       } else {
 273           tz = - TMath::Sqrt(1. - txy);
 274       }
 275
 276       p[0] =  p0 * fDX;
 277       p[1] =  p0 * fDY;
 278       p[2] =  p0 * tz;
 279
 280       origin[0] = fXS;
 281       origin[1] = fYS;
 282       origin[2] = 1950.;
 283
 284       PushTrack(fTrackIt , -1, fPdg , p, origin, polar, fTS - 1950./2.9979e10, kPNoProcess, nt, fWS);
 285
 286       Int_t nc = ReadNextParticle();
 287
 288       if (fLossID != oldID || nc == 0) {
 289           oldID = fLossID;
 290           break;
 291       }
 292   }
 293   SetHighWaterMark(nt);
 294 }
 295
 296
 297 Float_t AliGenHalo::GasPressureWeight(Float_t zPrimary)
 298 {
 299 //
 300 // Return z-dependent gasspressure weight = interaction rate [1/m/s].
 301 //
 302     Float_t weight = 0.;
 303     zPrimary /= 100.;        // m
 304     if (fRunPeriod < 5) {
 305         Float_t zAbs = TMath::Abs(zPrimary);
 306         if (zPrimary > 0.)
 307         {
 308             if (zAbs > fZ1[20]) {
 309                 weight = 2.e4;
 310             } else {
 311                 for (Int_t i = 1; i < 21; i++) {
 312                     if (zAbs < fZ1[i]) {
 313                         weight = fG1[i];
 314                         break;
 315                     }
 316                 }
 317             }
 318         } else {
 319             if (zAbs > fZ2[20]) {
 320                 weight = 2.e4;
 321             } else {
 322                 for (Int_t i = 1; i < 21; i++) {
 323                     if (zAbs < fZ2[i]) {
 324                     weight = fG2[i];
 325                     break;
 326                     }
 327                 }
 328             }
 329         }
 330     } else {
 331         Int_t index = TMath::BinarySearch(fGPASize, fZ1, zPrimary);
 332         weight = fG1[index];
 333     }
 334     return weight;
 335 }
 336
 337 void AliGenHalo::Draw(Option_t *)
 338 {
 339 // Draws the gas pressure distribution
 340     Float_t z[400];
 341     Float_t p[400];
 342
 343     for (Int_t i = 0; i < 400; i++)
 344     {
 345         z[i] = -20000. + Float_t(i) * 100;
 346         p[i] = GasPressureWeight(z[i]);
 347     }
 348
 349     TGraph*  gr = new TGraph(400, z, p);
 350     TCanvas* c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
 351     c1->cd();
 352     gr->Draw("AL");
 353 }
 354
 355 Int_t AliGenHalo::ReadNextParticle()
 356 {
 357     // Read the next particle from the file
 358     Int_t ncols = fscanf(fFile,"%d %f %f %d %f %d %f %f %f %f %f %f %f %f",
 359                    &fLossID, &fLossT0, &fLossZ, &fLossA, &fLossW, &fPdg, &fXS, &fYS, &fZS, &fDX, &fDY, &fEkin, &fTS, &fWS);
 360     fLossZ /= 10.;
 361     fXS    /= 10.;
 362     fYS    /= 10.;
 363     fZS    /= 10.;
 364     fTS    *= 1.e-9;
 365     return (ncols);
 366 }
 367
 368 void AliGenHalo::SkipEvents()
 369 {
 370   //
 371   // Skip the first fNskip events
 372   Int_t skip = fNskip;
 373   Int_t oldID = -1;
 374
 375   while (skip >= 0)
 376     {
 377       ReadNextParticle();
 378       if (oldID != fLossID) {
 379         oldID = fLossID;
 380         skip--;
 381       }
 382     }
 383 }
 384 void AliGenHalo::CountEvents()
 385 {
 386     Int_t nev = 0;
 387     Int_t oldID = -1;
 388     Int_t nc = 1;
 389     while (nc != -1)
 390     {
 391         nc = ReadNextParticle();
 392         if (oldID != fLossID) {
 393             oldID = fLossID;
 394             nev++;
 395             printf("Number of events %10d %10d \n", nev, oldID);
 396         }
 397     }
 398
 399
 400     rewind(fFile);
 401 }
 402
 403