Option to read gas pressure tables from LHC PR 674.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 1 Mar 2004 11:09:15 +0000 (11:09 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 1 Mar 2004 11:09:15 +0000 (11:09 +0000)
EVGEN/AliGenHaloProtvino.cxx
EVGEN/AliGenHaloProtvino.h

index 8ebb365..f7aec29 100644 (file)
@@ -24,6 +24,8 @@
 #include <stdlib.h>
 
 #include <TDatabasePDG.h>
+#include <TCanvas.h>
+#include <TGraph.h>
 #include <TPDGCode.h>
 #include <TSystem.h>
 
@@ -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
 {
     //
index c1b9804..4961d72 100644 (file)
@@ -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)