[u/mrichter/AliRoot.git] / EVGEN / AliGenHaloProtvino.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/*
$Log$
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.

*/

// Read background particles from a boundary source
// Very specialized generator to simulate background from beam halo.
// The input file is a text file specially prepared 
// for this purpose.
// Author: andreas.morsch@cern.ch

#include "AliGenHaloProtvino.h"
#include "AliRun.h"
#include "AliMC.h"
#include "AliPDG.h"

#include <TDatabasePDG.h>
#include <TSystem.h>
#include <stdlib.h>

 ClassImp(AliGenHaloProtvino)
     AliGenHaloProtvino::AliGenHaloProtvino()
	 :AliGenerator(-1)
{
// Constructor
    
    fName  = "HaloProtvino";
    fTitle = "Halo from LHC Tunnel";
//
//  Read all particles
    fNpart = -1;
    fFile  =  0;
    fSide  =  1;
//
    SetRunPeriod();
    SetTimePerEvent();
    SetAnalog(0);
}

AliGenHaloProtvino::AliGenHaloProtvino(Int_t npart)
    :AliGenerator(npart)
{
// 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()
{
// Destructor
}

//____________________________________________________________
void AliGenHaloProtvino::Init() 
{
// Initialisation
    fFile = fopen(fFileName,"r");
    if (fFile) {
	printf("\n File %s opened for reading, %p ! \n ",  fFileName.Data(), fFile);
    } else {
	printf("\n Opening of file %s failed,  %p ! \n ",  fFileName.Data(), fFile);
    }
//
//
//
//    Read file with gas pressure values
    char *name;
    
    name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
    FILE* file = fopen(name, "r");
    Float_t z;
    Int_t i, j;
    
//
//  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.;
	}
  }
//
// 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.;
	}
    }
//
//  Transform into 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 (j = 0; j <  5; j++) {
	pFlux[j] *= 1.e11/25.e-9;
	for (i = 0; i < 21; 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
	}
    }
    

    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]);
	sum1 += wgt1;
	sum2 += wgt2;
    }
    sum1/=250.;
    sum2/=250.;
    printf("\n %f %f \n \n", sum1, sum2);
}

//____________________________________________________________
void AliGenHaloProtvino::Generate()
{
// Generate from input file
 
  Float_t polar[3]= {0,0,0};
  Float_t origin[3];
  Float_t p[3], p0;
  Float_t tz, txy;
  Float_t amass;
  //
  Int_t ncols, nt;
  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
  
  while(1) {
//
// Load event into array
//
      ncols = fscanf(fFile,"%f %d %d %f %f %f %f %f %f",
		     &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;
  }
//
// Mean time between interactions
//
  Float_t dT = fTimePerEvent/nInt;   // sec 
  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[nprim]*ekin[nprim] + 2.*amass*ekin[nprim]);
      
      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[nprim];
      p[1] = p0*ty[nprim];
      p[2] =-p0*tz;
      
      origin[0] = vx[nprim];
      origin[1] = vy[nprim];
      origin[2] = -2196.5;

      //
      //
      // Particle weight

      Float_t originP[3] = {0., 0., 0.};
      originP[2] = zPrimary[nprim];
      
      Float_t pP[3] = {0., 0., 0.};
      Int_t ntP;
      
      if (fSide == -1) {
	  originP[2] = -zPrimary[nprim];
	  origin[2]  = -origin[2];
	  p[2]       = -p[2];
      }

      //
      // 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]*GassPressureWeight(zPrimary[nprim]);

      if (!fAnalog || gRandom->Rndm() < fParentWeight) {
//    Pass parent particle
//
	  SetTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
	  KeepTrack(ntP);
	  SetTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
      }

      //
      // Both sides are considered
      //

      if (fSide > 1) {
	  fParentWeight=wgt[nprim]*GassPressureWeight(-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);
	      KeepTrack(ntP);
	      SetTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
	  }
      }
      SetHighWaterMark(nt);
  }
  delete [] zPrimary;
  delete [] inuc;    
  delete [] ipart;   
  delete [] wgt;     
  delete [] ekin;    
  delete [] vx;      
  delete [] vy;      
  delete [] tx;      
  delete [] ty;      
}
 

AliGenHaloProtvino& AliGenHaloProtvino::operator=(const  AliGenHaloProtvino& rhs)
{
// Assignment operator
    return *this;
}


Float_t AliGenHaloProtvino::GassPressureWeight(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;
		}
	    }
	}
    } 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];
		    break;
		}
	    }
	}
    }
    return weight;
}

/*
# Title:    README file for the sources of IR8 machine induced background
# Author:   Vadim Talanov <Vadim.Talanov@cern.ch>
# Modified: 12-12-2000 

0. Overview

	There are three files, named ring.one.beta.[01,10,50].m, which
	contain the lists of background particles, induced by proton losses
	upstream of IP8 in the LHC ring one, for the beta* values of 1, 10
	and 50 m, respectively.

1. File contents

	Each line in the files contains the coordinates of particle track
	crossing with the infinite plane, positioned at z=-1m, together with
	the physical properties of corresponding particle, namely:

	S  - S coordinate of the primary interaction vertex, cm;
	N  - type of the gas nuclei at interaction, 1 is H, 2 - C and 3 - O;
	I  - particle ID in PDG particle numbering scheme;
	W  - particle weight;
	E  - particle kinetic energy, GeV;
	X  - x coordinate of the crossing point, cm;
	Y  - y coordinate of the crossing point, cm;
	Dx - x direction cosine;
	Dy - y direction cosine.

2. Normalisation

	Each file is given per unity of linear density of proton inelastic
	interactions with the gas nuclei, [1 inelastic interaction/m].

# ~/vtalanov/public/README.mib: the end.

*/
Commit	Line	Data
14ab5373	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/*
	17	$Log$
718e5512	18	Revision 1.9 2002/05/28 13:49:17 morsch
	19	- Udates for new pressure table
	20	- calculate time
	21	- first provisions for real events.
	22
7f98ef21	23	Revision 1.8 2002/03/22 13:00:25 morsch
	24	Initialize sum to 0. (Jiri Chudoba).
	25
c416706e	26	Revision 1.7 2002/02/21 16:09:45 morsch
	27	Move SetHighwaterMark() after last possible SetTrack()
	28
5bfdba97	29	Revision 1.6 2002/02/01 15:28:22 morsch
	30	Fastidious print statements removed.
	31
69818e67	32	Revision 1.5 2002/02/01 14:12:28 morsch
	33	Include new gas pressure estimates (LHC Pproject Note 273)
	34
eac2c8a3	35	Revision 1.4 2001/12/19 16:30:24 morsch
	36	Some bugs corrected and skip method added. (Rachid Guernane)
	37
05df024b	38	Revision 1.3 2001/07/27 17:09:36 morsch
	39	Use local SetTrack, KeepTrack and SetHighWaterMark methods
	40	to delegate either to local stack or to stack owned by AliRun.
	41	(Piotr Skowronski, A.M.)
	42
a99cf51f	43	Revision 1.2 2001/06/14 12:15:27 morsch
	44	Bugs corrected. SetSide() method added.
	45
154d524b	46	Revision 1.1 2001/01/23 15:04:33 morsch
	47	Generator to read beam halo file from Protvino group.
	48
14ab5373	49	*/
	50
	51	// Read background particles from a boundary source
	52	// Very specialized generator to simulate background from beam halo.
	53	// The input file is a text file specially prepared
	54	// for this purpose.
	55	// Author: andreas.morsch@cern.ch
	56
	57	#include "AliGenHaloProtvino.h"
	58	#include "AliRun.h"
	59	#include "AliMC.h"
	60	#include "AliPDG.h"
	61
	62	#include <TDatabasePDG.h>
eac2c8a3	63	#include <TSystem.h>
14ab5373	64	#include <stdlib.h>
	65
	66	ClassImp(AliGenHaloProtvino)
	67	AliGenHaloProtvino::AliGenHaloProtvino()
	68	:AliGenerator(-1)
	69	{
	70	// Constructor
154d524b	71
	72	fName = "HaloProtvino";
	73	fTitle = "Halo from LHC Tunnel";
14ab5373	74	//
14ab5373	75	// Read all particles
154d524b	76	fNpart = -1;
	77	fFile = 0;
	78	fSide = 1;
eac2c8a3	79	//
eac2c8a3	80	SetRunPeriod();
7f98ef21	81	SetTimePerEvent();
7f98ef21	82	SetAnalog(0);
14ab5373	83	}
	84
	85	AliGenHaloProtvino::AliGenHaloProtvino(Int_t npart)
	86	:AliGenerator(npart)
	87	{
	88	// Constructor
154d524b	89	fName = "Halo";
154d524b	90	fTitle= "Halo from LHC Tunnel";
14ab5373	91	//
154d524b	92	fNpart = npart;
	93	fFile = 0;
	94	fSide = 1;
eac2c8a3	95	//
eac2c8a3	96	SetRunPeriod();
7f98ef21	97	SetTimePerEvent();
7f98ef21	98	SetAnalog(0);
14ab5373	99	}
	100
	101	AliGenHaloProtvino::AliGenHaloProtvino(const AliGenHaloProtvino & HaloProtvino)
	102	{
	103	// copy constructor
	104	}
	105
	106
	107	//____________________________________________________________
	108	AliGenHaloProtvino::~AliGenHaloProtvino()
	109	{
	110	// Destructor
	111	}
	112
	113	//____________________________________________________________
	114	void AliGenHaloProtvino::Init()
	115	{
	116	// Initialisation
154d524b	117	fFile = fopen(fFileName,"r");
	118	if (fFile) {
	119	printf("\n File %s opened for reading, %p ! \n ", fFileName.Data(), fFile);
	120	} else {
	121	printf("\n Opening of file %s failed, %p ! \n ", fFileName.Data(), fFile);
	122	}
eac2c8a3	123	//
	124	//
	125	//
	126	// Read file with gas pressure values
	127	char *name;
	128
	129	name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
	130	FILE* file = fopen(name, "r");
	131	Float_t z;
	132	Int_t i, j;
	133
	134	//
	135	// Ring 1
	136	//
7f98ef21	137	for (i = 0; i < 21; i++)
eac2c8a3	138	{
	139	fscanf(file, "%f %f %f %f %f %f", &z,
	140	&fG1[i][0], &fG1[i][1], &fG1[i][2], &fG1[i][3], &fG1[i][4]);
	141	if (i > 0) {
	142	fZ1[i] = fZ1[i-1] + z;
	143	} else {
	144	fZ1[i] = 20.;
	145	}
	146	}
	147	//
	148	// Ring 2
	149	//
	150	for (i = 0; i < 21; i++)
	151	{
	152	fscanf(file, "%f %f %f %f %f %f", &z,
	153	&fG2[i][0], &fG2[i][1], &fG2[i][2], &fG2[i][3], &fG2[i][4]);
	154	if (i > 0) {
	155	fZ2[i] = fZ2[i-1] + z;
	156	} else {
	157	fZ2[i] = 20.;
	158	}
	159	}
	160	//
	161	// Transform into interaction rates
	162	//
	163	const Float_t crossSection = 0.094e-28; // m^2
718e5512	164	Float_t pFlux[5] = {0.2, 0.2, 0.3, 0.3, 1.0};
eac2c8a3	165
eac2c8a3	166	for (j = 0; j < 5; j++) {
718e5512	167	pFlux[j] *= 1.e11/25.e-9;
eac2c8a3	168	for (i = 0; i < 21; i++)
eac2c8a3	169	{
718e5512	170	fG1[i][j] = fG1[i][j] * crossSection * pFlux[j]; // 1/m/s
718e5512	171	fG2[i][j] = fG2[i][j] * crossSection * pFlux[j]; // 1/m/s
eac2c8a3	172	}
eac2c8a3	173	}
7f98ef21	174
eac2c8a3	175
7f98ef21	176	Float_t sum1 = 0.;
7f98ef21	177	Float_t sum2 = 0.;
eac2c8a3	178
7f98ef21	179	for (Int_t i = 0; i < 300; i++) {
eac2c8a3	180	Float_t z = 20.+i*1.;
eac2c8a3	181	z*=100;
7f98ef21	182	Float_t wgt1 = GassPressureWeight(z);
	183	Float_t wgt2 = GassPressureWeight(-z);
	184	// printf("weight: %f %f %f %f %f \n", z, wgt1, wgt2, fZ1[20], fZ2[20]);
	185	sum1 += wgt1;
	186	sum2 += wgt2;
eac2c8a3	187	}
7f98ef21	188	sum1/=250.;
	189	sum2/=250.;
	190	printf("\n %f %f \n \n", sum1, sum2);
14ab5373	191	}
	192
	193	//____________________________________________________________
	194	void AliGenHaloProtvino::Generate()
	195	{
	196	// Generate from input file
14ab5373	197
	198	Float_t polar[3]= {0,0,0};
	199	Float_t origin[3];
	200	Float_t p[3], p0;
7f98ef21	201	Float_t tz, txy;
14ab5373	202	Float_t amass;
14ab5373	203	//
7f98ef21	204	Int_t ncols, nt;
05df024b	205	Int_t nskip = 0;
154d524b	206	Int_t nread = 0;
7f98ef21	207
	208	Float_t* zPrimary = new Float_t [fNpart];
	209	Int_t * inuc = new Int_t [fNpart];
	210	Int_t * ipart = new Int_t [fNpart];
	211	Float_t* wgt = new Float_t [fNpart];
	212	Float_t* ekin = new Float_t [fNpart];
	213	Float_t* vx = new Float_t [fNpart];
	214	Float_t* vy = new Float_t [fNpart];
	215	Float_t* tx = new Float_t [fNpart];
	216	Float_t* ty = new Float_t [fNpart];
	217
	218	Float_t zVertexOld = -1.e10;
	219	Int_t nInt = 0; // Counts number of interactions
	220
14ab5373	221	while(1) {
7f98ef21	222	//
	223	// Load event into array
	224	//
154d524b	225	ncols = fscanf(fFile,"%f %d %d %f %f %f %f %f %f",
7f98ef21	226	&zPrimary[nread], &inuc[nread], &ipart[nread], &wgt[nread],
	227	&ekin[nread], &vx[nread], &vy[nread],
	228	&tx[nread], &ty[nread]);
05df024b	229
14ab5373	230	if (ncols < 0) break;
7f98ef21	231	// Skip fNskip events
05df024b	232	nskip++;
05df024b	233	if (fNpart !=-1 && nskip <= fNskip) continue;
7f98ef21	234	// Count interactions
	235	if (zPrimary[nread] != zVertexOld) {
	236	nInt++;
	237	zVertexOld = zPrimary[nread];
	238	}
	239	// Count tracks
05df024b	240	nread++;
14ab5373	241	if (fNpart !=-1 && nread > fNpart) break;
7f98ef21	242	}
	243	//
	244	// Mean time between interactions
	245	//
	246	Float_t dT = fTimePerEvent/nInt; // sec
	247	Float_t t = 0; // sec
	248
	249	//
	250	// Loop over primaries
	251	//
	252	zVertexOld = -1.e10;
	253	Double_t arg = 0.;
	254
	255	for (Int_t nprim = 0; nprim < fNpart; nprim++)
	256	{
	257	amass = TDatabasePDG::Instance()->GetParticle(ipart[nprim])->Mass();
14ab5373	258
	259	//
	260	// Momentum vector
	261	//
7f98ef21	262	p0=sqrt(ekin[nprim]ekin[nprim] + 2.amass*ekin[nprim]);
14ab5373	263
7f98ef21	264	txy=TMath::Sqrt(tx[nprim]tx[nprim]+ty[nprim]ty[nprim]);
14ab5373	265	if (txy == 1.) {
	266	tz=0;
	267	} else {
	268	tz=-TMath::Sqrt(1.-txy);
	269	}
	270
7f98ef21	271	p[0] = p0*tx[nprim];
	272	p[1] = p0*ty[nprim];
	273	p[2] =-p0*tz;
05df024b	274
7f98ef21	275	origin[0] = vx[nprim];
7f98ef21	276	origin[1] = vy[nprim];
05df024b	277	origin[2] = -2196.5;
14ab5373	278
	279	//
	280	//
	281	// Particle weight
	282
154d524b	283	Float_t originP[3] = {0., 0., 0.};
7f98ef21	284	originP[2] = zPrimary[nprim];
154d524b	285
	286	Float_t pP[3] = {0., 0., 0.};
	287	Int_t ntP;
	288
	289	if (fSide == -1) {
7f98ef21	290	originP[2] = -zPrimary[nprim];
154d524b	291	origin[2] = -origin[2];
	292	p[2] = -p[2];
	293	}
05df024b	294
14ab5373	295	//
7f98ef21	296	// Time
	297	//
	298	if (zPrimary[nprim] != zVertexOld) {
	299	while(arg==0.) arg = gRandom->Rndm();
	300	t -= dT*TMath::Log(arg); // (sec)
	301	zVertexOld = zPrimary[nprim];
	302	}
	303
	304	// Get statistical weight according to local gas-pressure
	305	//
	306	fParentWeight=wgt[nprim]*GassPressureWeight(zPrimary[nprim]);
	307
	308	if (!fAnalog \|\| gRandom->Rndm() < fParentWeight) {
	309	// Pass parent particle
	310	//
	311	SetTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
	312	KeepTrack(ntP);
	313	SetTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
	314	}
	315
	316	//
	317	// Both sides are considered
	318	//
154d524b	319
05df024b	320	if (fSide > 1) {
7f98ef21	321	fParentWeight=wgt[nprim]*GassPressureWeight(-zPrimary[nprim]);
	322	if (!fAnalog \|\| gRandom->Rndm() < fParentWeight) {
	323	origin[2] = -origin[2];
	324	originP[2] = -originP[2];
	325	p[2]=-p[2];
	326	SetTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
	327	KeepTrack(ntP);
	328	SetTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
	329	}
05df024b	330	}
5bfdba97	331	SetHighWaterMark(nt);
14ab5373	332	}
7f98ef21	333	delete [] zPrimary;
	334	delete [] inuc;
	335	delete [] ipart;
	336	delete [] wgt;
	337	delete [] ekin;
	338	delete [] vx;
	339	delete [] vy;
	340	delete [] tx;
	341	delete [] ty;
14ab5373	342	}
	343
	344
	345	AliGenHaloProtvino& AliGenHaloProtvino::operator=(const AliGenHaloProtvino& rhs)
	346	{
	347	// Assignment operator
	348	return *this;
	349	}
	350
	351
eac2c8a3	352
14ab5373	353	Float_t AliGenHaloProtvino::GassPressureWeight(Float_t zPrimary)
14ab5373	354	{
eac2c8a3	355	//
	356	// Return z-dependent gasspressure weight = interaction rate [1/m/s].
	357	//
	358	Float_t weight = 0.;
	359	zPrimary/=100.; // m
	360	Float_t zAbs = TMath::Abs(zPrimary);
7f98ef21	361	if (zPrimary > 0.)
eac2c8a3	362	{
7f98ef21	363	if (zAbs > fZ1[20]) {
eac2c8a3	364	weight = 2.e4;
eac2c8a3	365	} else {
7f98ef21	366	for (Int_t i = 1; i < 21; i++) {
eac2c8a3	367	if (zAbs < fZ1[i]) {
	368	weight = fG1[i][fRunPeriod];
	369	break;
	370	}
	371	}
	372	}
	373	} else {
	374	if (zAbs > fZ2[20]) {
	375	weight = 2.e4;
	376	} else {
	377	for (Int_t i = 1; i < 21; i++) {
	378	if (zAbs < fZ2[i]) {
	379	weight = fG2[i][fRunPeriod];
	380	break;
	381	}
	382	}
	383	}
	384	}
eac2c8a3	385	return weight;
14ab5373	386	}
	387
	388	/*
	389	# Title: README file for the sources of IR8 machine induced background
	390	# Author: Vadim Talanov <Vadim.Talanov@cern.ch>
	391	# Modified: 12-12-2000
	392
	393	0. Overview
	394
	395	There are three files, named ring.one.beta.[01,10,50].m, which
	396	contain the lists of background particles, induced by proton losses
	397	upstream of IP8 in the LHC ring one, for the beta* values of 1, 10
	398	and 50 m, respectively.
	399
	400	1. File contents
	401
	402	Each line in the files contains the coordinates of particle track
	403	crossing with the infinite plane, positioned at z=-1m, together with
	404	the physical properties of corresponding particle, namely:
	405
	406	S - S coordinate of the primary interaction vertex, cm;
	407	N - type of the gas nuclei at interaction, 1 is H, 2 - C and 3 - O;
	408	I - particle ID in PDG particle numbering scheme;
	409	W - particle weight;
	410	E - particle kinetic energy, GeV;
	411	X - x coordinate of the crossing point, cm;
	412	Y - y coordinate of the crossing point, cm;
	413	Dx - x direction cosine;
	414	Dy - y direction cosine.
	415
	416	2. Normalisation
	417
	418	Each file is given per unity of linear density of proton inelastic
	419	interactions with the gas nuclei, [1 inelastic interaction/m].
	420
	421	# ~/vtalanov/public/README.mib: the end.
	422
	423	*/
	424
	425
	426
	427