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