[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.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.

*/

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