X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EVGEN%2FAliGenHaloProtvino.cxx;h=4466a194e37e2ee9f556687316d333efa80215ca;hb=5b64c50285c4ceb169ad1d72481f03e471f52714;hp=e7fdcbb17bc96f7fab3b658b01744f20bd93a4ca;hpb=acc86a241efbb2413429532a3006c66e29f2ab8c;p=u%2Fmrichter%2FAliRoot.git diff --git a/EVGEN/AliGenHaloProtvino.cxx b/EVGEN/AliGenHaloProtvino.cxx index e7fdcbb17bc..4466a194e37 100644 --- a/EVGEN/AliGenHaloProtvino.cxx +++ b/EVGEN/AliGenHaloProtvino.cxx @@ -13,49 +13,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.11 2002/10/14 14:55:35 hristov -Merging the VirtualMC branch to the main development branch (HEAD) - -Revision 1.8.4.1 2002/07/24 08:56:28 alibrary -Updating EVGEN on TVirtulaMC - -Revision 1.10 2002/05/28 14:28:12 morsch -Correct LHC nominal current profile (1st year 20%, 2nd year 30%, 3rd year 100 %) - -Revision 1.9 2002/05/28 13:49:17 morsch -- Udates for new pressure table -- calculate time -- first provisions for real events. - -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. @@ -63,17 +21,33 @@ Generator to read beam halo file from Protvino group. // for this purpose. // Author: andreas.morsch@cern.ch -#include "AliGenHaloProtvino.h" -#include "AliRun.h" -#include "AliPDG.h" +#include #include +#include +#include +#include #include -#include - 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 @@ -82,36 +56,32 @@ Generator to read beam halo file from Protvino group. // // 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() { @@ -124,78 +94,133 @@ void AliGenHaloProtvino::Init() // 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 - 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] * crossSection * pFlux[j]; // 1/m/s - fG2[i][j] = fG2[i][j] * crossSection * pFlux[j]; // 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; } //____________________________________________________________ @@ -210,7 +235,7 @@ void AliGenHaloProtvino::Generate() 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]; @@ -225,6 +250,7 @@ void AliGenHaloProtvino::Generate() Float_t zVertexOld = -1.e10; Int_t nInt = 0; // Counts number of interactions + Float_t wwgt = 0.; while(1) { // @@ -246,12 +272,15 @@ void AliGenHaloProtvino::Generate() } // 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 // @@ -311,31 +340,32 @@ 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 // - 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; @@ -346,53 +376,71 @@ void AliGenHaloProtvino::Generate() 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