1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
26 #include <TDatabasePDG.h>
32 #include "AliGenHalo.h"
33 #include "AliGenEventHeader.h"
39 AliGenHalo::AliGenHalo()
70 fTitle = "Halo from LHC Beam";
77 AliGenHalo::AliGenHalo(Int_t npart)
107 fTitle= "Halo from LHC Beam";
114 //____________________________________________________________
115 AliGenHalo::~AliGenHalo()
120 //____________________________________________________________
121 void AliGenHalo::Init()
124 fFile = fopen(fFileName,"r");
129 printf("\n File %s opened for reading, %p ! \n ", fFileName.Data(), (void*)fFile);
131 printf("\n Opening of file %s failed, %p ! \n ", fFileName.Data(), (void*)fFile);
136 // Skip the first fNskip events
142 // Read file with gas pressure values
144 if (fRunPeriod < 5) {
145 name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
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");
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");
159 fG1 = new Float_t[fGPASize];
160 fZ1 = new Float_t[fGPASize];
162 Fatal("Init()", "No gas pressure file for given run period !");
167 if (name) file = fopen(name, "r");
169 AliError("No gas pressure file");
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};
181 if (fRunPeriod < 5) {
186 for (i = 0; i < fGPASize; i++)
188 ir = fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
191 fG1[i] = p[fRunPeriod];
194 fZ1[i] = fZ1[i-1] + z;
202 for (i = 0; i < fGPASize; i++)
204 ir = fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
207 fG2[i] = p[fRunPeriod];
209 fZ2[i] = fZ2[i-1] + z;
217 for (i = 0; i < fGPASize; i++)
219 fG1[i] = fG1[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
220 fG2[i] = fG2[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
224 for (i = 0; i < fGPASize; i++)
226 ir = fscanf(file, "%f %e %e %e %e %e", &z, &p[0], &p[1], &p[2], &p[3], &p[4]);
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.;
240 // Transform into interaction rates
249 for (Int_t iz = 0; iz < 300; iz++) {
250 Float_t zpos = 20. + iz * 1.;
252 Float_t wgt1 = GasPressureWeight( zpos);
253 Float_t wgt2 = GasPressureWeight(-zpos);
259 printf("\n %f %f \n \n", sum1, sum2);
263 //____________________________________________________________
264 void AliGenHalo::Generate()
266 // Generate by reading particles from input file
268 Float_t polar[3]= {0., 0., 0.};
275 static Bool_t first = kTRUE;
276 static Int_t oldID = -1;
279 if (first && (fNskip == 0)) ReadNextParticle();
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);
292 tz = - TMath::Sqrt(1. - txy);
303 PushTrack(fTrackIt , -1, fPdg , p, origin, polar, fTS - 1950./2.9979e10, kPNoProcess, nt, fWS);
305 Int_t nc = ReadNextParticle();
307 if (fLossID != oldID || nc == 0) {
313 SetHighWaterMark(nt);
314 AliGenEventHeader* header = new AliGenEventHeader("HALO");
315 header->SetNProduced(np);
316 // Passes header either to the container or to gAlice
318 header->SetName(fName);
319 fContainer->AddHeader(header);
321 gAlice->SetGenEventHeader(header);
326 Float_t AliGenHalo::GasPressureWeight(Float_t zPrimary)
329 // Return z-dependent gasspressure weight = interaction rate [1/m/s].
332 zPrimary /= 100.; // m
333 if (fRunPeriod < 5) {
334 Float_t zAbs = TMath::Abs(zPrimary);
337 if (zAbs > fZ1[20]) {
340 for (Int_t i = 1; i < 21; i++) {
348 if (zAbs > fZ2[20]) {
351 for (Int_t i = 1; i < 21; i++) {
360 Int_t index = TMath::BinarySearch(fGPASize, fZ1, zPrimary);
366 void AliGenHalo::Draw(Option_t *)
368 // Draws the gas pressure distribution
372 for (Int_t i = 0; i < 400; i++)
374 z[i] = -20000. + Float_t(i) * 100;
375 p[i] = GasPressureWeight(z[i]);
378 TGraph* gr = new TGraph(400, z, p);
379 TCanvas* c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
384 Int_t AliGenHalo::ReadNextParticle()
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);
397 void AliGenHalo::SkipEvents()
400 // Skip the first fNskip events
407 if (oldID != fLossID) {
413 void AliGenHalo::CountEvents()
415 // Count total number of events
421 nc = ReadNextParticle();
422 if (oldID != fLossID) {
425 printf("Number of events %10d %10d \n", nev, oldID);