X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EVGEN%2FAliGenHaloProtvino.cxx;h=4466a194e37e2ee9f556687316d333efa80215ca;hb=0300efeba6711342232ca4eed23219dab37af3f8;hp=c287244be0f2f249f7b22182c376533215386cba;hpb=69818e67d6bfaf109e02c9d9df084eb3010e31ed;p=u%2Fmrichter%2FAliRoot.git diff --git a/EVGEN/AliGenHaloProtvino.cxx b/EVGEN/AliGenHaloProtvino.cxx index c287244be0f..4466a194e37 100644 --- a/EVGEN/AliGenHaloProtvino.cxx +++ b/EVGEN/AliGenHaloProtvino.cxx @@ -13,26 +13,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -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. @@ -40,18 +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 "AliMC.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 @@ -60,32 +56,32 @@ Generator to read beam halo file from Protvino group. // // Read all particles fNpart = -1; - fFile = 0; - fSide = 1; -// - SetRunPeriod(); + 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(); -} - -AliGenHaloProtvino::AliGenHaloProtvino(const AliGenHaloProtvino & HaloProtvino) -{ -// copy constructor + SetAnalog(0); } - //____________________________________________________________ AliGenHaloProtvino::~AliGenHaloProtvino() { @@ -98,74 +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 < 20; 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++) { - if (i <20) - 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 sum; - for (Int_t i = 0; i < 250; i++) { - Float_t z = 20.+i*1.; - z*=100; - Float_t wgt = GassPressureWeight(z); -// printf("\n %f %f", z, wgt); - sum+=wgt; + + 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; } - sum/=250.; - printf("\n %f \n \n", sum); + sum1/=250.; + sum2/=250.; + printf("\n %f %f \n \n", sum1, sum2); + delete file; } //____________________________________________________________ @@ -176,52 +231,86 @@ void AliGenHaloProtvino::Generate() 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 zPrimary; + Float_t tz, txy; Float_t amass; - Int_t inuc; // - Int_t ipart, ncols, nt; - Int_t nskip = 0; + Int_t ncols, nt; + static Int_t nskip = 0; Int_t nread = 0; + + Float_t* zPrimary = new Float_t [fNpart]; + Int_t * inuc = new Int_t [fNpart]; + Int_t * ipart = new Int_t [fNpart]; + Float_t* wgt = new Float_t [fNpart]; + Float_t* ekin = new Float_t [fNpart]; + Float_t* vx = new Float_t [fNpart]; + Float_t* vy = new Float_t [fNpart]; + Float_t* tx = new Float_t [fNpart]; + Float_t* ty = new Float_t [fNpart]; + + Float_t zVertexOld = -1.e10; + Int_t nInt = 0; // Counts number of interactions + Float_t wwgt = 0.; + while(1) { +// +// Load event into array +// ncols = fscanf(fFile,"%f %d %d %f %f %f %f %f %f", - &zPrimary, &inuc, &ipart, &wgt, - &ekin, &origin[0], &origin[1], - &tx, &ty); -/* - printf(" \n %f %d %d %f %f %f %f %f %f", - zPrimary, inuc, ipart, wgt, - ekin, origin[0], origin[1], - tx, ty); -*/ + &zPrimary[nread], &inuc[nread], &ipart[nread], &wgt[nread], + &ekin[nread], &vx[nread], &vy[nread], + &tx[nread], &ty[nread]); if (ncols < 0) break; - +// Skip fNskip events nskip++; if (fNpart !=-1 && nskip <= fNskip) continue; - +// Count interactions + if (zPrimary[nread] != zVertexOld) { + nInt++; + zVertexOld = zPrimary[nread]; + } +// Count tracks nread++; - if (fNpart !=-1 && nread > fNpart) break; + if (fNpart !=-1 && nread >= fNpart) break; + } +// +// Mean time between interactions +// - amass = TDatabasePDG::Instance()->GetParticle(ipart)->Mass(); + Float_t dT = 0.; // sec + if (nInt > 0) + dT = fTimePerEvent/nInt; + Float_t t = 0; // sec + +// +// Loop over primaries +// + zVertexOld = -1.e10; + Double_t arg = 0.; + + for (Int_t nprim = 0; nprim < fNpart; nprim++) + { + amass = TDatabasePDG::Instance()->GetParticle(ipart[nprim])->Mass(); // // Momentum vector // - p0=sqrt(ekin*ekin + 2.*amass*ekin); + p0=sqrt(ekin[nprim]*ekin[nprim] + 2.*amass*ekin[nprim]); - txy=TMath::Sqrt(tx*tx+ty*ty); + txy=TMath::Sqrt(tx[nprim]*tx[nprim]+ty[nprim]*ty[nprim]); if (txy == 1.) { tz=0; } else { tz=-TMath::Sqrt(1.-txy); } - p[0]=p0*tx; - p[1]=p0*ty; - p[2]=-p0*tz; + p[0] = p0*tx[nprim]; + p[1] = p0*ty[nprim]; + p[2] =-p0*tz; + origin[0] = vx[nprim]; + origin[1] = vy[nprim]; origin[2] = -2196.5; // @@ -229,84 +318,129 @@ void AliGenHaloProtvino::Generate() // Particle weight Float_t originP[3] = {0., 0., 0.}; - originP[2] = zPrimary; + originP[2] = zPrimary[nprim]; Float_t pP[3] = {0., 0., 0.}; Int_t ntP; if (fSide == -1) { - originP[2] = -zPrimary; + originP[2] = -zPrimary[nprim]; origin[2] = -origin[2]; p[2] = -p[2]; } - SetTrack(0,-1,kProton,pP,originP,polar,0,kPNoProcess,ntP); - KeepTrack(ntP); - fParentWeight=wgt*GassPressureWeight(zPrimary); - SetTrack(fTrackIt,ntP,ipart,p,origin,polar,0,kPNoProcess,nt,fParentWeight); - SetHighWaterMark(nt); - // - // Assume particles come from two directions with same probability + // Time + // + if (zPrimary[nprim] != zVertexOld) { + while(arg==0.) arg = gRandom->Rndm(); + t -= dT*TMath::Log(arg); // (sec) + zVertexOld = zPrimary[nprim]; + } + +// Get statistical weight according to local gas-pressure +// + fParentWeight=wgt[nprim]*GasPressureWeight(zPrimary[nprim]); + + if (!fAnalog || gRandom->Rndm() < fParentWeight) { +// Pass parent particle +// + PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight); + KeepTrack(ntP); + PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight); + } + // + // Both sides are considered + // - // Both Side are considered if (fSide > 1) { - origin[2]=-origin[2]; - p[2]=-p[2]; - fParentWeight=wgt*GassPressureWeight(-zPrimary); - SetTrack(fTrackIt,ntP,ipart,p,origin,polar,0,kPNoProcess,nt,fParentWeight); + fParentWeight=wgt[nprim]*GasPressureWeight(-zPrimary[nprim]); + if (!fAnalog || gRandom->Rndm() < fParentWeight) { + origin[2] = -origin[2]; + originP[2] = -originP[2]; + p[2]=-p[2]; + PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight); + KeepTrack(ntP); + PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight); + } } + wwgt += fParentWeight; + SetHighWaterMark(nt); } + delete [] zPrimary; + delete [] inuc; + delete [] ipart; + delete [] wgt; + delete [] ekin; + delete [] vx; + delete [] vy; + delete [] tx; + 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); - zPrimary = TMath::Abs(zPrimary); - - if (zPrimary < 0.) - { - if (zAbs > fZ1[19]) { - weight = 2.e4; - } else { - for (Int_t i = 1; i < 20; 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