From d41a82acbec38799a7e16bf413a64baef52d2695 Mon Sep 17 00:00:00 2001 From: morsch Date: Mon, 1 Mar 2004 11:09:15 +0000 Subject: [PATCH] Option to read gas pressure tables from LHC PR 674. --- EVGEN/AliGenHaloProtvino.cxx | 185 ++++++++++++++++++++++++----------- EVGEN/AliGenHaloProtvino.h | 14 +-- 2 files changed, 137 insertions(+), 62 deletions(-) diff --git a/EVGEN/AliGenHaloProtvino.cxx b/EVGEN/AliGenHaloProtvino.cxx index 8ebb365cac4..f7aec29e99c 100644 --- a/EVGEN/AliGenHaloProtvino.cxx +++ b/EVGEN/AliGenHaloProtvino.cxx @@ -24,6 +24,8 @@ #include #include +#include +#include #include #include @@ -94,63 +96,106 @@ void AliGenHaloProtvino::Init() // // // Read file with gas pressure values - char *name; - - name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" ); + 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, 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}; + + 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++) + { + 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 < 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++) + { + 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.; + } } - } // -// Transform into interaction rates +// Interaction rates // - const Float_t kCrossSection = 0.094e-28; // m^2 - Float_t pFlux[5] = {0.2, 0.2, 0.3, 0.3, 1.0}; + 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++) { - pFlux[j] *= 1.e11/25.e-9; - for (i = 0; i < 21; i++) + } else { + for (i = 0; i < fGPASize; i++) { - fG1[i][j] = fG1[i][j] * kCrossSection * pFlux[j]; // 1/m/s - fG2[i][j] = fG2[i][j] * kCrossSection * pFlux[j]; // 1/m/s + 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 i = 0; i < 300; i++) { Float_t z = 20.+i*1.; z*=100; - Float_t wgt1 = GassPressureWeight(z); - Float_t wgt2 = GassPressureWeight(-z); + Float_t wgt1 = GasPressureWeight(z); + Float_t wgt2 = GasPressureWeight(-z); // printf("weight: %f %f %f %f %f \n", z, wgt1, wgt2, fZ1[20], fZ2[20]); sum1 += wgt1; sum2 += wgt2; @@ -187,6 +232,7 @@ void AliGenHaloProtvino::Generate() Float_t zVertexOld = -1.e10; Int_t nInt = 0; // Counts number of interactions + Float_t Wgt = 0.; while(1) { // @@ -273,7 +319,7 @@ void AliGenHaloProtvino::Generate() // 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 @@ -282,13 +328,12 @@ void AliGenHaloProtvino::Generate() KeepTrack(ntP); 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]; @@ -298,6 +343,8 @@ void AliGenHaloProtvino::Generate() PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight); } } + Wgt += fParentWeight; + SetHighWaterMark(nt); } delete [] zPrimary; @@ -308,7 +355,9 @@ void AliGenHaloProtvino::Generate() delete [] vx; delete [] vy; delete [] tx; - delete [] ty; + delete [] ty; + printf("Total weight %f\n\n", Wgt); + } @@ -321,41 +370,65 @@ AliGenHaloProtvino& AliGenHaloProtvino::operator=(const AliGenHaloProtvino& rhs -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() +{ +// 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"); +} + + void AliGenHaloProtvino::Copy(TObject&) const { // diff --git a/EVGEN/AliGenHaloProtvino.h b/EVGEN/AliGenHaloProtvino.h index c1b9804f326..4961d728da0 100644 --- a/EVGEN/AliGenHaloProtvino.h +++ b/EVGEN/AliGenHaloProtvino.h @@ -16,7 +16,7 @@ class AliGenHaloProtvino : public AliGenerator { public: - enum constants{kY1Day0, kY1Day70, kY2D0, kY2D10, kY3D90}; + enum constants{kY1Day0, kY1Day70, kY2D0, kY2D10, kY3D90, kLHCPR674Startup, kLHCPR674Conditioned}; AliGenHaloProtvino(); AliGenHaloProtvino(Int_t npart); AliGenHaloProtvino(const AliGenHaloProtvino &HaloProtvino); @@ -24,13 +24,12 @@ public: virtual void Init(); virtual void SetFileName(TString filename) {fFileName=TString(filename);} virtual void Generate(); - virtual Float_t GassPressureWeight(Float_t zPrimary); + 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(); AliGenHaloProtvino & operator=(const AliGenHaloProtvino & rhs); protected: @@ -40,8 +39,11 @@ protected: 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[21], fZ2[21]; // ! z-positions for gas pressure tables - Float_t fG1[21][5], fG2[21][5]; // ! gas pressures + 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 private: void Copy(TObject&) const; ClassDef(AliGenHaloProtvino,1) // LHC background boundary source (Protvino Group results) -- 2.43.0