* provided "as is" without express or implied warranty. *
**************************************************************************/
-/*
-$Log$
-Revision 1.8 2002/03/22 13:00:25 morsch
-Initialize sum to 0. (Jiri Chudoba).
-
-Revision 1.7 2002/02/21 16:09:45 morsch
-Move SetHighwaterMark() after last possible SetTrack()
-
-Revision 1.6 2002/02/01 15:28:22 morsch
-Fastidious print statements removed.
-
-Revision 1.5 2002/02/01 14:12:28 morsch
-Include new gas pressure estimates (LHC Pproject Note 273)
-
-Revision 1.4 2001/12/19 16:30:24 morsch
-Some bugs corrected and skip method added. (Rachid Guernane)
-
-Revision 1.3 2001/07/27 17:09:36 morsch
-Use local SetTrack, KeepTrack and SetHighWaterMark methods
-to delegate either to local stack or to stack owned by AliRun.
-(Piotr Skowronski, A.M.)
-
-Revision 1.2 2001/06/14 12:15:27 morsch
-Bugs corrected. SetSide() method added.
-
-Revision 1.1 2001/01/23 15:04:33 morsch
-Generator to read beam halo file from Protvino group.
-
-*/
+/* $Id$ */
// Read background particles from a boundary source
// Very specialized generator to simulate background from beam halo.
// for this purpose.
// Author: andreas.morsch@cern.ch
-#include "AliGenHaloProtvino.h"
-#include "AliRun.h"
-#include "AliMC.h"
-#include "AliPDG.h"
+#include <stdlib.h>
#include <TDatabasePDG.h>
+#include <TCanvas.h>
+#include <TGraph.h>
+#include <TPDGCode.h>
#include <TSystem.h>
-#include <stdlib.h>
- ClassImp(AliGenHaloProtvino)
- AliGenHaloProtvino::AliGenHaloProtvino()
- :AliGenerator(-1)
+#include "AliLog.h"
+#include "AliGenHaloProtvino.h"
+#include "AliRun.h"
+
+ClassImp(AliGenHaloProtvino)
+
+AliGenHaloProtvino::AliGenHaloProtvino()
+ :AliGenerator(-1),
+ fFile(0),
+ fFileName(0),
+ fSide(1),
+ fRunPeriod(kY3D90),
+ fTimePerEvent(1.e-4),
+ fNskip(0),
+ fZ1(0),
+ fZ2(0),
+ fG1(0),
+ fG2(0),
+ fGPASize(0)
{
// Constructor
//
// Read all particles
fNpart = -1;
- fFile = 0;
- fSide = 1;
-//
- SetRunPeriod();
- SetTimePerEvent();
SetAnalog(0);
}
AliGenHaloProtvino::AliGenHaloProtvino(Int_t npart)
- :AliGenerator(npart)
+ :AliGenerator(npart),
+ fFile(0),
+ fFileName(0),
+ fSide(1),
+ fRunPeriod(kY3D90),
+ fTimePerEvent(1.e-4),
+ fNskip(0),
+ fZ1(0),
+ fZ2(0),
+ fG1(0),
+ fG2(0),
+ fGPASize(0)
{
// Constructor
fName = "Halo";
fTitle= "Halo from LHC Tunnel";
//
fNpart = npart;
- fFile = 0;
- fSide = 1;
//
- SetRunPeriod();
- SetTimePerEvent();
SetAnalog(0);
}
-AliGenHaloProtvino::AliGenHaloProtvino(const AliGenHaloProtvino & HaloProtvino)
-{
-// copy constructor
-}
-
-
//____________________________________________________________
AliGenHaloProtvino::~AliGenHaloProtvino()
{
// Initialisation
fFile = fopen(fFileName,"r");
if (fFile) {
- printf("\n File %s opened for reading, %p ! \n ", fFileName.Data(), fFile);
+ printf("\n File %s opened for reading, %p ! \n ", fFileName.Data(), (void*)fFile);
} else {
- printf("\n Opening of file %s failed, %p ! \n ", fFileName.Data(), fFile);
+ printf("\n Opening of file %s failed, %p ! \n ", fFileName.Data(), (void*)fFile);
}
//
//
//
// Read file with gas pressure values
- char *name;
-
- name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
- FILE* file = fopen(name, "r");
+ char *name = 0;
+ if (fRunPeriod < 5) {
+ name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
+ fGPASize = 21;
+ fG1 = new Float_t[fGPASize];
+ fG2 = new Float_t[fGPASize];
+ fZ1 = new Float_t[fGPASize];
+ fZ2 = new Float_t[fGPASize];
+ } else if (fRunPeriod == 5) {
+ name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_startup.dat");
+ fGPASize = 18853;
+ fG1 = new Float_t[fGPASize];
+ fZ1 = new Float_t[fGPASize];
+ } else if (fRunPeriod ==6) {
+ name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_conditioned.dat");
+ fGPASize = 12719;
+ fG1 = new Float_t[fGPASize];
+ fZ1 = new Float_t[fGPASize];
+ } else {
+ Fatal("Init()", "No gas pressure file for given run period !");
+ }
+
+ FILE* file = 0;
+ if (name) file = fopen(name, "r");
+ if (!file) {
+ AliError("No gas pressure file");
+ return;
+ }
+
Float_t z;
- Int_t i, j;
-
+ Int_t i;
+ Float_t p[5];
+
+ const Float_t kCrossSection = 0.094e-28; // m^2
+ const Float_t kFlux = 1.e11 / 25.e-9; // protons/s
+ Float_t pFlux[5] = {0.2, 0.2, 0.3, 0.3, 1.0};
+
+ Int_t ncols = 0;
+ if (fRunPeriod < 5) {
//
// Ring 1
//
- for (i = 0; i < 21; i++)
- {
- fscanf(file, "%f %f %f %f %f %f", &z,
- &fG1[i][0], &fG1[i][1], &fG1[i][2], &fG1[i][3], &fG1[i][4]);
- if (i > 0) {
- fZ1[i] = fZ1[i-1] + z;
- } else {
- fZ1[i] = 20.;
+
+ for (i = 0; i < fGPASize; i++)
+ {
+ ncols = fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
+ if (ncols<0) break;
+
+ fG1[i] = p[fRunPeriod];
+
+ if (i > 0) {
+ fZ1[i] = fZ1[i-1] + z;
+ } else {
+ fZ1[i] = 20.;
+ }
}
- }
//
// Ring 2
//
- for (i = 0; i < 21; i++)
- {
- fscanf(file, "%f %f %f %f %f %f", &z,
- &fG2[i][0], &fG2[i][1], &fG2[i][2], &fG2[i][3], &fG2[i][4]);
- if (i > 0) {
- fZ2[i] = fZ2[i-1] + z;
- } else {
- fZ2[i] = 20.;
+ for (i = 0; i < fGPASize; i++)
+ {
+ ncols = fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
+ if (ncols<0) break;
+
+ fG2[i] = p[fRunPeriod];
+ if (i > 0) {
+ fZ2[i] = fZ2[i-1] + z;
+ } else {
+ fZ2[i] = 20.;
+ }
}
- }
//
-// Transform into interaction rates
+// Interaction rates
//
- const Float_t crossSection = 0.094e-28; // m^2
- const Float_t pFlux = 1.e11/25.e-9; // 1/s
+ for (i = 0; i < fGPASize; i++)
+ {
+ fG1[i] = fG1[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
+ fG2[i] = fG2[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
+ }
- for (j = 0; j < 5; j++) {
- for (i = 0; i < 21; i++)
+ } else {
+ for (i = 0; i < fGPASize; i++)
{
- fG1[i][j] = fG1[i][j] * crossSection * pFlux; // 1/m/s
- fG2[i][j] = fG2[i][j] * crossSection * pFlux; // 1/m/s
+ ncols = fscanf(file, "%f %e %e %e %e %e", &z, &p[0], &p[1], &p[2], &p[3], &p[4]);
+ if (ncols<0) break;
+
+ z /= 1000.;
+ fG1[i] = p[4] * kCrossSection * kFlux; // 1/m/s
+ // 1/3 of nominal intensity at startup
+ if (fRunPeriod == kLHCPR674Startup) fG1[i] /= 3.;
+ fZ1[i] = z;
}
}
+
+
+//
+// Transform into interaction rates
+//
+
+
+
+
Float_t sum1 = 0.;
Float_t sum2 = 0.;
- for (Int_t i = 0; i < 300; i++) {
- Float_t z = 20.+i*1.;
- z*=100;
- Float_t wgt1 = GassPressureWeight(z);
- Float_t wgt2 = GassPressureWeight(-z);
-// printf("weight: %f %f %f %f %f \n", z, wgt1, wgt2, fZ1[20], fZ2[20]);
+ for (Int_t iz = 0; iz < 300; iz++) {
+ Float_t zpos = 20. + iz * 1.;
+ zpos *= 100;
+ Float_t wgt1 = GasPressureWeight( zpos);
+ Float_t wgt2 = GasPressureWeight(-zpos);
sum1 += wgt1;
sum2 += wgt2;
}
sum1/=250.;
sum2/=250.;
printf("\n %f %f \n \n", sum1, sum2);
+ delete file;
}
//____________________________________________________________
Float_t amass;
//
Int_t ncols, nt;
- Int_t nskip = 0;
+ static Int_t nskip = 0;
Int_t nread = 0;
Float_t* zPrimary = new Float_t [fNpart];
Float_t zVertexOld = -1.e10;
Int_t nInt = 0; // Counts number of interactions
+ Float_t wwgt = 0.;
while(1) {
//
}
// Count tracks
nread++;
- if (fNpart !=-1 && nread > fNpart) break;
+ if (fNpart !=-1 && nread >= fNpart) break;
}
//
// Mean time between interactions
//
- Float_t dT = fTimePerEvent/nInt; // sec
+
+ Float_t dT = 0.; // sec
+ if (nInt > 0)
+ dT = fTimePerEvent/nInt;
Float_t t = 0; // sec
//
// Get statistical weight according to local gas-pressure
//
- fParentWeight=wgt[nprim]*GassPressureWeight(zPrimary[nprim]);
+ fParentWeight=wgt[nprim]*GasPressureWeight(zPrimary[nprim]);
if (!fAnalog || gRandom->Rndm() < fParentWeight) {
// Pass parent particle
//
- SetTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
+ PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
KeepTrack(ntP);
- SetTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
+ PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
}
-
//
// Both sides are considered
//
if (fSide > 1) {
- fParentWeight=wgt[nprim]*GassPressureWeight(-zPrimary[nprim]);
+ fParentWeight=wgt[nprim]*GasPressureWeight(-zPrimary[nprim]);
if (!fAnalog || gRandom->Rndm() < fParentWeight) {
origin[2] = -origin[2];
originP[2] = -originP[2];
p[2]=-p[2];
- SetTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
+ PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
KeepTrack(ntP);
- SetTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
+ PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
}
}
+ wwgt += fParentWeight;
+
SetHighWaterMark(nt);
}
delete [] zPrimary;
delete [] vx;
delete [] vy;
delete [] tx;
- delete [] ty;
+ delete [] ty;
+ printf("Total weight %f\n\n", wwgt);
+
}
-AliGenHaloProtvino& AliGenHaloProtvino::operator=(const AliGenHaloProtvino& rhs)
-{
-// Assignment operator
- return *this;
-}
-
-
-
-Float_t AliGenHaloProtvino::GassPressureWeight(Float_t zPrimary)
+Float_t AliGenHaloProtvino::GasPressureWeight(Float_t zPrimary)
{
//
// Return z-dependent gasspressure weight = interaction rate [1/m/s].
//
Float_t weight = 0.;
- zPrimary/=100.; // m
- Float_t zAbs = TMath::Abs(zPrimary);
- if (zPrimary > 0.)
- {
- if (zAbs > fZ1[20]) {
- weight = 2.e4;
- } else {
- for (Int_t i = 1; i < 21; i++) {
- if (zAbs < fZ1[i]) {
- weight = fG1[i][fRunPeriod];
- break;
+ zPrimary /= 100.; // m
+ if (fRunPeriod < 5) {
+ Float_t zAbs = TMath::Abs(zPrimary);
+ if (zPrimary > 0.)
+ {
+ if (zAbs > fZ1[20]) {
+ weight = 2.e4;
+ } else {
+ for (Int_t i = 1; i < 21; i++) {
+ if (zAbs < fZ1[i]) {
+ weight = fG1[i];
+ break;
+ }
}
}
- }
- } else {
- if (zAbs > fZ2[20]) {
- weight = 2.e4;
} else {
- for (Int_t i = 1; i < 21; i++) {
- if (zAbs < fZ2[i]) {
- weight = fG2[i][fRunPeriod];
+ if (zAbs > fZ2[20]) {
+ weight = 2.e4;
+ } else {
+ for (Int_t i = 1; i < 21; i++) {
+ if (zAbs < fZ2[i]) {
+ weight = fG2[i];
break;
+ }
}
}
}
+ } else {
+ Int_t index = TMath::BinarySearch(fGPASize, fZ1, zPrimary);
+ weight = fG1[index];
}
return weight;
}
+void AliGenHaloProtvino::Draw(Option_t *)
+{
+// Draws the gas pressure distribution
+ Float_t z[400];
+ Float_t p[400];
+
+ for (Int_t i = 0; i < 400; i++)
+ {
+ z[i] = -20000. + Float_t(i) * 100;
+ p[i] = GasPressureWeight(z[i]);
+ }
+
+ TGraph* gr = new TGraph(400, z, p);
+ TCanvas* c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
+ c1->cd();
+ gr->Draw("AL");
+}
+
+
/*
# Title: README file for the sources of IR8 machine induced background
# Author: Vadim Talanov <Vadim.Talanov@cern.ch>
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff