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