X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EVGEN%2FAliGenHaloProtvino.cxx;h=5234b8d15e90ffe36726b38b2b0d57356bb23ef0;hb=de4ab12f58c7ef589671d8aded3b26afa972501a;hp=5b2a60d8536bbbe8803816bf2c29b9a65cad64f9;hpb=88cb7938ca21d4a80991d4e7aa564008c29340f7;p=u%2Fmrichter%2FAliRoot.git diff --git a/EVGEN/AliGenHaloProtvino.cxx b/EVGEN/AliGenHaloProtvino.cxx index 5b2a60d8536..5234b8d15e9 100644 --- a/EVGEN/AliGenHaloProtvino.cxx +++ b/EVGEN/AliGenHaloProtvino.cxx @@ -24,15 +24,29 @@ #include #include +#include +#include #include #include #include "AliGenHaloProtvino.h" #include "AliRun.h" - ClassImp(AliGenHaloProtvino) - AliGenHaloProtvino::AliGenHaloProtvino() - :AliGenerator(-1) +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 @@ -41,36 +55,32 @@ // // 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() { @@ -83,78 +93,121 @@ 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" ); + 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 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 + 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); -// 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; } //____________________________________________________________ @@ -169,7 +222,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]; @@ -184,6 +237,7 @@ void AliGenHaloProtvino::Generate() Float_t zVertexOld = -1.e10; Int_t nInt = 0; // Counts number of interactions + Float_t wwgt = 0.; while(1) { // @@ -270,31 +324,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; @@ -305,53 +360,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