/* $Id$ */
-// Read background particles from a boundary source
-// Very specialized generator to simulate background from beam halo.
-// The input file is a text file specially prepared
-// for this purpose.
+// Read beam halo background particles from a boundary source
+// Boundary source is in the LHCb format http://www.hep.manchester.ac.uk/u/robert/LHC_backgrounds/Note-MIBStudies.pdf
+// and has been provided by Robert Appleby
// Author: andreas.morsch@cern.ch
+
#include <stdlib.h>
#include <TDatabasePDG.h>
+#include <TCanvas.h>
+#include <TGraph.h>
#include <TPDGCode.h>
+#include <TSystem.h>
#include "AliGenHalo.h"
#include "AliRun.h"
ClassImp(AliGenHalo)
AliGenHalo::AliGenHalo()
- :AliGenerator(-1),
- fp(0),
- fFileName(0)
+ :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),
+ fLossID(0),
+ fLossA(0),
+ fPdg(0),
+ fLossT0(0),
+ fLossZ(0),
+ fLossW(0),
+ fXS(0),
+ fYS(0),
+ fZS(0),
+ fDX(0),
+ fDY(0),
+ fEkin(0),
+ fTS(0),
+ fWS(0)
{
// Constructor
- fName="Halo";
- fTitle="Halo from LHC Tunnel";
+
+ fName = "Halo";
+ fTitle = "Halo from LHC Beam";
//
// Read all particles
- fNpart=-1;
+ fNpart = -1;
+ SetAnalog(0);
}
AliGenHalo::AliGenHalo(Int_t npart)
:AliGenerator(npart),
- fp(0),
- fFileName(0)
+ fFile(0),
+ fFileName(0),
+ fSide(1),
+ fRunPeriod(kY3D90),
+ fTimePerEvent(1.e-4),
+ fNskip(0),
+ fZ1(0),
+ fZ2(0),
+ fG1(0),
+ fG2(0),
+ fGPASize(0),
+ fLossID(0),
+ fLossA(0),
+ fPdg(0),
+ fLossT0(0),
+ fLossZ(0),
+ fLossW(0),
+ fXS(0),
+ fYS(0),
+ fZS(0),
+ fDX(0),
+ fDY(0),
+ fEkin(0),
+ fTS(0),
+ fWS(0)
{
// Constructor
- fName="Halo";
- fTitle="Halo from LHC Tunnel";
+ fName = "Halo";
+ fTitle= "Halo from LHC Beam";
//
-// Read all particles
- fNpart=-1;
+ fNpart = npart;
+//
+ SetAnalog(0);
}
//____________________________________________________________
void AliGenHalo::Init()
{
// Initialisation
+ fFile = fopen(fFileName,"r");
+ if (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(), (void*)fFile);
+ }
+//
+//
+//
+// Read file with gas pressure values
+ 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 = fopen(name, "r");
+ Float_t z;
+ 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};
+
+ if (fRunPeriod < 5) {
+//
+// Ring 1
+//
+
+ for (i = 0; i < fGPASize; i++)
+ {
+ fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
+ fG1[i] = p[fRunPeriod];
+
+ if (i > 0) {
+ fZ1[i] = fZ1[i-1] + z;
+ } else {
+ fZ1[i] = 20.;
+ }
+ }
+//
+// Ring 2
+//
+ for (i = 0; i < fGPASize; i++)
+ {
+ fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
+ fG2[i] = p[fRunPeriod];
+ if (i > 0) {
+ fZ2[i] = fZ2[i-1] + z;
+ } else {
+ fZ2[i] = 20.;
+ }
+ }
+//
+// Interaction rates
+//
+ 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
+ }
+
+ } else {
+ for (i = 0; i < fGPASize; i++)
+ {
+ fscanf(file, "%f %e %e %e %e %e", &z, &p[0], &p[1], &p[2], &p[3], &p[4]);
+ 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 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);
}
//____________________________________________________________
void AliGenHalo::Generate()
{
-// Generate from input file
- fp = fopen(fFileName,"r");
- if (fp) {
- printf("\n File %s opened for reading ! \n ", (char*) &fFileName);
- } else {
- printf("\n Opening of file %s failed ! \n ", (char*) &fFileName);
- }
-//
-// MARS particle codes
- const Int_t kmars[12]={0,kProton,kNeutron,kPiPlus,kPiMinus,kKPlus,kKMinus,
- kMuonPlus,kMuonMinus,kGamma,kElectron,kPositron};
+// Generate by reading particles from input file
- Float_t polar[3]= {0,0,0};
+ Float_t polar[3]= {0., 0., 0.};
Float_t origin[3];
Float_t p[3], p0;
- Float_t ekin, wgt, tx, ty, tz, txy;
- Float_t amass;
+ Float_t tz, txy;
+ Float_t mass;
//
- Int_t ipart, ncols, nt;
+ Int_t nt;
+ static Bool_t first = kTRUE;
+ static Int_t oldID = -1;
+//
+ if (first) {
+ ReadNextParticle();
+ first = kFALSE;
+ oldID = fLossID;
+ }
+
- Int_t nread=0;
- origin[2]=2650;
while(1) {
- ncols = fscanf(fp,"%i %f %f %f %f %f %f",
- &ipart, &ekin, &wgt,
- &origin[0], &origin[1],
- &tx, &ty);
- if (ncols < 0) break;
- nread++;
- if (fNpart !=-1 && nread > fNpart) break;
- ipart = kmars[ipart];
- amass = TDatabasePDG::Instance()->GetParticle(ipart)->Mass();
- p0=sqrt(ekin*ekin + 2.*amass);
-
- txy=TMath::Sqrt(tx*tx+ty*ty);
+ // Push particle to stack
+ mass = TDatabasePDG::Instance()->GetParticle(fPdg)->Mass();
+ p0 = TMath::Sqrt(fEkin * fEkin + 2.* mass * fEkin);
+ txy = TMath::Sqrt(fDX * fDX + fDY * fDY);
if (txy == 1.) {
- tz=0;
+ tz = 0;
} else {
- tz=-TMath::Sqrt(1.-txy);
+ tz = - TMath::Sqrt(1. - txy);
}
- p[0]=p0*tx;
- p[1]=p0*ty;
- p[2]=p0*tz;
- fParentWeight=wgt;
- PushTrack(fTrackIt,-1,ipart,p,origin,polar,0,kPNoProcess,nt,fParentWeight);
-// PushTrack(fTrackIt,-1,ipart,p,origin,polar,0,"Halo+",nt,fParentWeight);
- origin[2]=-origin[2];
- p[2]=-p[2];
- PushTrack(fTrackIt,-1,ipart,p,origin,polar,0,kPNoProcess,nt,fParentWeight);
-// PushTrack(fTrackIt,-1,ipart,p,origin,polar,0,"Halo-",nt,fParentWeight);
- origin[2]=-origin[2];
- p[2]=-p[2];
- }
-}
-
-
+ p[0] = p0 * fDX;
+ p[1] = p0 * fDY;
+ p[2] = p0 * tz;
+ origin[0] = fXS;
+ origin[1] = fYS;
+ origin[2] = 1950.;
+ PushTrack(fTrackIt , -1, fPdg , p, origin, polar, fTS, kPNoProcess, nt, fWS);
+ Int_t nc = ReadNextParticle();
+
+ if (fLossID != oldID || nc == 0) {
+ oldID = fLossID;
+ break;
+ }
+ }
+ SetHighWaterMark(nt);
+}
+
+Float_t AliGenHalo::GasPressureWeight(Float_t zPrimary)
+{
+//
+// Return z-dependent gasspressure weight = interaction rate [1/m/s].
+//
+ Float_t weight = 0.;
+ 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];
+ break;
+ }
+ }
+ }
+ }
+ } else {
+ Int_t index = TMath::BinarySearch(fGPASize, fZ1, zPrimary);
+ weight = fG1[index];
+ }
+ return weight;
+}
+void AliGenHalo::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");
+}
+Int_t AliGenHalo::ReadNextParticle()
+{
+ // Read the next particle from the file
+ Int_t ncols = fscanf(fFile,"%d %f %f %d %f %d %f %f %f %f %f %f %f %f",
+ &fLossID, &fLossT0, &fLossZ, &fLossA, &fLossW, &fPdg, &fXS, &fYS, &fZS, &fDX, &fDY, &fEkin, &fTS, &fWS);
+ fLossZ /= 10.;
+ fXS /= 10.;
+ fYS /= 10.;
+ fZS /= 10.;
+ fTS *= 1.e-9;
+ return (ncols);
+}
#include "AliGenerator.h"
#include <TString.h>
-// Read background particles from a boundary source
-// Very specialized generator to simulate background from beam halo.
+// Read beam halo background particles from a boundary source
+// Boundary source is in the LHCb format
+// http://www.hep.manchester.ac.uk/u/robert/LHC_backgrounds/Note-MIBStudies.pdf
+// and has been provided by Robert Appleby
// Author: andreas.morsch@cern.ch
class AliGenHalo : public AliGenerator
{
- public:
+public:
+ enum constants{kY1Day0, kY1Day70, kY2D0, kY2D10, kY3D90, kLHCPR674Startup, kLHCPR674Conditioned};
AliGenHalo();
AliGenHalo(Int_t npart);
virtual ~AliGenHalo();
virtual void Init();
virtual void SetFileName(TString filename) {fFileName=TString(filename);}
virtual void Generate();
+ virtual Float_t GasPressureWeight(Float_t zPrimary);
+ virtual void SetSide(Int_t flag = 1) {fSide = flag;}
+ virtual void SetNskip(Int_t nskip) {fNskip = nskip;}
+ virtual void SetRunPeriod(Int_t t = kY3D90) {fRunPeriod = t;}
+ virtual void SetTimePerEvent(Float_t t = 1.e-4) {fTimePerEvent = t;}
+ virtual void Draw(Option_t * opt="");
+ private:
+ virtual Int_t ReadNextParticle();
protected:
- FILE *fp; // ! Pointer to file
+ FILE* fFile; // ! Pointer to file
TString fFileName; // Choose the file
+ Int_t fSide; // Muon arm side (1) / Castor side (-1)
+ Int_t fRunPeriod; // LHC Running Period
+ Float_t fTimePerEvent; // Time corresponding to one event [s]
+ Int_t fNskip; // Number of entries to skip
+ Float_t* fZ1; // ! z-positions for gas pressure tables
+ Float_t* fZ2; // ! z-positions for gas pressure tables
+ Float_t* fG1; // ! gas pressures
+ Float_t* fG2; // ! gas pressures
+ Int_t fGPASize; // ! Size of arrays
+ Int_t fLossID; // ! unique loss ID
+ Int_t fLossA; // ! atomic number of scatterer
+ Int_t fPdg; // ! pdg code
+ Float_t fLossT0; // ! relative time
+ Float_t fLossZ; // ! scaterring position
+ Float_t fLossW; // ! weight of proton loss
+ Float_t fXS; // ! x-position on scoring plane
+ Float_t fYS; // ! y-position on scoring plane
+ Float_t fZS; // ! z-position on scoring plane
+ Float_t fDX; // ! direction cosine x
+ Float_t fDY; // ! direction cosine y
+ Float_t fEkin; // ! kinetic energy
+ Float_t fTS; // ! relative arrival time
+ Float_t fWS; // ! weight
private:
AliGenHalo(const AliGenHalo &Halo);
AliGenHalo & operator=(const AliGenHalo & rhs);
-
- ClassDef(AliGenHalo,1) // LHC background boundary source (MARS input)
+
+ ClassDef(AliGenHalo,1) // LHC beam halo boundary source
};
#endif