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