newly implemented model for slow nucleon production
[u/mrichter/AliRoot.git] / EVGEN / AliGenHaloProtvino.cxx
index 2eb40d2..4466a19 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/*
-$Log$
-Revision 1.12  2002/11/21 16:22:35  alibrary
-Removing AliMCProcess and AliMC
-
-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.
@@ -69,15 +24,30 @@ Generator to read beam halo file from Protvino group.
 #include <stdlib.h>
 
 #include <TDatabasePDG.h>
+#include <TCanvas.h>
+#include <TGraph.h>
 #include <TPDGCode.h>
 #include <TSystem.h>
 
+#include "AliLog.h"
 #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
     
@@ -86,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()
 {
@@ -128,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;
 }
 
 //____________________________________________________________
@@ -214,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];
@@ -229,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) {
 //
@@ -250,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
   
 //
@@ -315,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;
@@ -350,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 <Vadim.Talanov@cern.ch>