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