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

/* $Id$ */

// 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 <stdlib.h>

#include <TDatabasePDG.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TPDGCode.h>
#include <TSystem.h>

#include "AliGenHaloProtvino.h"
#include "AliRun.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):
    AliGenerator(HaloProtvino)
{
// Copy constructor
    HaloProtvino.Copy(*this);
}


//____________________________________________________________
AliGenHaloProtvino::~AliGenHaloProtvino()
{
// Destructor
}

//____________________________________________________________
void AliGenHaloProtvino::Init() 
{
// Initialisation
    fFile = fopen(fFileName,"r");
    if (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(), (void*)fFile);
    }
//
//
//
//    Read file with gas pressure values
    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;
    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 < 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 < 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.;
	    }
	}
//
// Interaction rates
//
	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
	}

    } else {
	for (i = 0; i < fGPASize; i++) 
	{
	    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 = 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;
    }
    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;
  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[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]*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
      //

      if (fSide > 1) {
	  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
    rhs.Copy(*this);
    return *this;
}


Float_t AliGenHaloProtvino::GasPressureWeight(Float_t zPrimary)
{
//
// Return z-dependent gasspressure weight = interaction rate [1/m/s].
//
    Float_t weight = 0.;
    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];
		    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");
}


void AliGenHaloProtvino::Copy(TObject&) const
{
    //
    // Copy 
    //
    Fatal("Copy","Not implemented!\n");
}


/*
# 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
	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	/* $Id$ */
	17
	18	// Read background particles from a boundary source
	19	// Very specialized generator to simulate background from beam halo.
	20	// The input file is a text file specially prepared
	21	// for this purpose.
	22	// Author: andreas.morsch@cern.ch
	23
	24	#include <stdlib.h>
	25
	26	#include <TDatabasePDG.h>
	27	#include <TCanvas.h>
	28	#include <TGraph.h>
	29	#include <TPDGCode.h>
	30	#include <TSystem.h>
	31
	32	#include "AliGenHaloProtvino.h"
	33	#include "AliRun.h"
	34
	35	ClassImp(AliGenHaloProtvino)
	36
	37	AliGenHaloProtvino::AliGenHaloProtvino()
	38	:AliGenerator(-1)
	39	{
	40	// Constructor
	41
	42	fName = "HaloProtvino";
	43	fTitle = "Halo from LHC Tunnel";
	44	//
	45	// Read all particles
	46	fNpart = -1;
	47	fFile = 0;
	48	fSide = 1;
	49	//
	50	SetRunPeriod();
	51	SetTimePerEvent();
	52	SetAnalog(0);
	53	}
	54
	55	AliGenHaloProtvino::AliGenHaloProtvino(Int_t npart)
	56	:AliGenerator(npart)
	57	{
	58	// Constructor
	59	fName = "Halo";
	60	fTitle= "Halo from LHC Tunnel";
	61	//
	62	fNpart = npart;
	63	fFile = 0;
	64	fSide = 1;
	65	//
	66	SetRunPeriod();
	67	SetTimePerEvent();
	68	SetAnalog(0);
	69	}
	70
	71	AliGenHaloProtvino::AliGenHaloProtvino(const AliGenHaloProtvino & HaloProtvino):
	72	AliGenerator(HaloProtvino)
	73	{
	74	// Copy constructor
	75	HaloProtvino.Copy(*this);
	76	}
	77
	78
	79	//____________________________________________________________
	80	AliGenHaloProtvino::~AliGenHaloProtvino()
	81	{
	82	// Destructor
	83	}
	84
	85	//____________________________________________________________
	86	void AliGenHaloProtvino::Init()
	87	{
	88	// Initialisation
	89	fFile = fopen(fFileName,"r");
	90	if (fFile) {
	91	printf("\n File %s opened for reading, %p ! \n ", fFileName.Data(), (void*)fFile);
	92	} else {
	93	printf("\n Opening of file %s failed, %p ! \n ", fFileName.Data(), (void*)fFile);
	94	}
	95	//
	96	//
	97	//
	98	// Read file with gas pressure values
	99	char *name = 0;
	100	if (fRunPeriod < 5) {
	101	name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
	102	fGPASize = 21;
	103	fG1 = new Float_t[fGPASize];
	104	fG2 = new Float_t[fGPASize];
	105	fZ1 = new Float_t[fGPASize];
	106	fZ2 = new Float_t[fGPASize];
	107	} else if (fRunPeriod == 5) {
	108	name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_startup.dat");
	109	fGPASize = 18853;
	110	fG1 = new Float_t[fGPASize];
	111	fZ1 = new Float_t[fGPASize];
	112	} else if (fRunPeriod ==6) {
	113	name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_conditioned.dat");
	114	fGPASize = 12719;
	115	fG1 = new Float_t[fGPASize];
	116	fZ1 = new Float_t[fGPASize];
	117	} else {
	118	Fatal("Init()", "No gas pressure file for given run period !");
	119	}
	120
	121
	122	FILE* file = fopen(name, "r");
	123	Float_t z;
	124	Int_t i;
	125	Float_t p[5];
	126
	127	const Float_t kCrossSection = 0.094e-28; // m^2
	128	const Float_t kFlux = 1.e11 / 25.e-9; // protons/s
	129	Float_t pFlux[5] = {0.2, 0.2, 0.3, 0.3, 1.0};
	130
	131	if (fRunPeriod < 5) {
	132	//
	133	// Ring 1
	134	//
	135
	136	for (i = 0; i < fGPASize; i++)
	137	{
	138	fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
	139	fG1[i] = p[fRunPeriod];
	140
	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 < fGPASize; i++)
	151	{
	152	fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
	153	fG2[i] = p[fRunPeriod];
	154	if (i > 0) {
	155	fZ2[i] = fZ2[i-1] + z;
	156	} else {
	157	fZ2[i] = 20.;
	158	}
	159	}
	160	//
	161	// Interaction rates
	162	//
	163	for (i = 0; i < fGPASize; i++)
	164	{
	165	fG1[i] = fG1[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
	166	fG2[i] = fG2[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
	167	}
	168
	169	} else {
	170	for (i = 0; i < fGPASize; i++)
	171	{
	172	fscanf(file, "%f %e %e %e %e %e", &z, &p[0], &p[1], &p[2], &p[3], &p[4]);
	173	z /= 1000.;
	174	fG1[i] = p[4] * kCrossSection * kFlux; // 1/m/s
	175	// 1/3 of nominal intensity at startup
	176	if (fRunPeriod == kLHCPR674Startup) fG1[i] /= 3.;
	177	fZ1[i] = z;
	178	}
	179	}
	180
	181
	182
	183
	184	//
	185	// Transform into interaction rates
	186	//
	187
	188
	189
	190
	191	Float_t sum1 = 0.;
	192	Float_t sum2 = 0.;
	193
	194	for (Int_t i = 0; i < 300; i++) {
	195	Float_t z = 20.+i*1.;
	196	z*=100;
	197	Float_t wgt1 = GasPressureWeight(z);
	198	Float_t wgt2 = GasPressureWeight(-z);
	199	// printf("weight: %f %f %f %f %f \n", z, wgt1, wgt2, fZ1[20], fZ2[20]);
	200	sum1 += wgt1;
	201	sum2 += wgt2;
	202	}
	203	sum1/=250.;
	204	sum2/=250.;
	205	printf("\n %f %f \n \n", sum1, sum2);
	206	}
	207
	208	//____________________________________________________________
	209	void AliGenHaloProtvino::Generate()
	210	{
	211	// Generate from input file
	212
	213	Float_t polar[3]= {0,0,0};
	214	Float_t origin[3];
	215	Float_t p[3], p0;
	216	Float_t tz, txy;
	217	Float_t amass;
	218	//
	219	Int_t ncols, nt;
	220	static Int_t nskip = 0;
	221	Int_t nread = 0;
	222
	223	Float_t* zPrimary = new Float_t [fNpart];
	224	Int_t * inuc = new Int_t [fNpart];
	225	Int_t * ipart = new Int_t [fNpart];
	226	Float_t* wgt = new Float_t [fNpart];
	227	Float_t* ekin = new Float_t [fNpart];
	228	Float_t* vx = new Float_t [fNpart];
	229	Float_t* vy = new Float_t [fNpart];
	230	Float_t* tx = new Float_t [fNpart];
	231	Float_t* ty = new Float_t [fNpart];
	232
	233	Float_t zVertexOld = -1.e10;
	234	Int_t nInt = 0; // Counts number of interactions
	235	Float_t wwgt = 0.;
	236
	237	while(1) {
	238	//
	239	// Load event into array
	240	//
	241	ncols = fscanf(fFile,"%f %d %d %f %f %f %f %f %f",
	242	&zPrimary[nread], &inuc[nread], &ipart[nread], &wgt[nread],
	243	&ekin[nread], &vx[nread], &vy[nread],
	244	&tx[nread], &ty[nread]);
	245
	246	if (ncols < 0) break;
	247	// Skip fNskip events
	248	nskip++;
	249	if (fNpart !=-1 && nskip <= fNskip) continue;
	250	// Count interactions
	251	if (zPrimary[nread] != zVertexOld) {
	252	nInt++;
	253	zVertexOld = zPrimary[nread];
	254	}
	255	// Count tracks
	256	nread++;
	257	if (fNpart !=-1 && nread > fNpart) break;
	258	}
	259	//
	260	// Mean time between interactions
	261	//
	262	Float_t dT = fTimePerEvent/nInt; // sec
	263	Float_t t = 0; // sec
	264
	265	//
	266	// Loop over primaries
	267	//
	268	zVertexOld = -1.e10;
	269	Double_t arg = 0.;
	270
	271	for (Int_t nprim = 0; nprim < fNpart; nprim++)
	272	{
	273	amass = TDatabasePDG::Instance()->GetParticle(ipart[nprim])->Mass();
	274
	275	//
	276	// Momentum vector
	277	//
	278	p0=sqrt(ekin[nprim]ekin[nprim] + 2.amass*ekin[nprim]);
	279
	280	txy=TMath::Sqrt(tx[nprim]tx[nprim]+ty[nprim]ty[nprim]);
	281	if (txy == 1.) {
	282	tz=0;
	283	} else {
	284	tz=-TMath::Sqrt(1.-txy);
	285	}
	286
	287	p[0] = p0*tx[nprim];
	288	p[1] = p0*ty[nprim];
	289	p[2] =-p0*tz;
	290
	291	origin[0] = vx[nprim];
	292	origin[1] = vy[nprim];
	293	origin[2] = -2196.5;
	294
	295	//
	296	//
	297	// Particle weight
	298
	299	Float_t originP[3] = {0., 0., 0.};
	300	originP[2] = zPrimary[nprim];
	301
	302	Float_t pP[3] = {0., 0., 0.};
	303	Int_t ntP;
	304
	305	if (fSide == -1) {
	306	originP[2] = -zPrimary[nprim];
	307	origin[2] = -origin[2];
	308	p[2] = -p[2];
	309	}
	310
	311	//
	312	// Time
	313	//
	314	if (zPrimary[nprim] != zVertexOld) {
	315	while(arg==0.) arg = gRandom->Rndm();
	316	t -= dT*TMath::Log(arg); // (sec)
	317	zVertexOld = zPrimary[nprim];
	318	}
	319
	320	// Get statistical weight according to local gas-pressure
	321	//
	322	fParentWeight=wgt[nprim]*GasPressureWeight(zPrimary[nprim]);
	323
	324	if (!fAnalog \|\| gRandom->Rndm() < fParentWeight) {
	325	// Pass parent particle
	326	//
	327	PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
	328	KeepTrack(ntP);
	329	PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
	330	}
	331	//
	332	// Both sides are considered
	333	//
	334
	335	if (fSide > 1) {
	336	fParentWeight=wgt[nprim]*GasPressureWeight(-zPrimary[nprim]);
	337	if (!fAnalog \|\| gRandom->Rndm() < fParentWeight) {
	338	origin[2] = -origin[2];
	339	originP[2] = -originP[2];
	340	p[2]=-p[2];
	341	PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
	342	KeepTrack(ntP);
	343	PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
	344	}
	345	}
	346	wwgt += fParentWeight;
	347
	348	SetHighWaterMark(nt);
	349	}
	350	delete [] zPrimary;
	351	delete [] inuc;
	352	delete [] ipart;
	353	delete [] wgt;
	354	delete [] ekin;
	355	delete [] vx;
	356	delete [] vy;
	357	delete [] tx;
	358	delete [] ty;
	359	printf("Total weight %f\n\n", wwgt);
	360
	361	}
	362
	363
	364	AliGenHaloProtvino& AliGenHaloProtvino::operator=(const AliGenHaloProtvino& rhs)
	365	{
	366	// Assignment operator
	367	rhs.Copy(*this);
	368	return *this;
	369	}
	370
	371
	372
	373	Float_t AliGenHaloProtvino::GasPressureWeight(Float_t zPrimary)
	374	{
	375	//
	376	// Return z-dependent gasspressure weight = interaction rate [1/m/s].
	377	//
	378	Float_t weight = 0.;
	379	zPrimary /= 100.; // m
	380	if (fRunPeriod < 5) {
	381	Float_t zAbs = TMath::Abs(zPrimary);
	382	if (zPrimary > 0.)
	383	{
	384	if (zAbs > fZ1[20]) {
	385	weight = 2.e4;
	386	} else {
	387	for (Int_t i = 1; i < 21; i++) {
	388	if (zAbs < fZ1[i]) {
	389	weight = fG1[i];
	390	break;
	391	}
	392	}
	393	}
	394	} else {
	395	if (zAbs > fZ2[20]) {
	396	weight = 2.e4;
	397	} else {
	398	for (Int_t i = 1; i < 21; i++) {
	399	if (zAbs < fZ2[i]) {
	400	weight = fG2[i];
	401	break;
	402	}
	403	}
	404	}
	405	}
	406	} else {
	407	Int_t index = TMath::BinarySearch(fGPASize, fZ1, zPrimary);
	408	weight = fG1[index];
	409	}
	410	return weight;
	411	}
	412
	413	void AliGenHaloProtvino::Draw(Option_t *)
	414	{
	415	// Draws the gas pressure distribution
	416	Float_t z[400];
	417	Float_t p[400];
	418
	419	for (Int_t i = 0; i < 400; i++)
	420	{
	421	z[i] = -20000. + Float_t(i) * 100;
	422	p[i] = GasPressureWeight(z[i]);
	423	}
	424
	425	TGraph* gr = new TGraph(400, z, p);
	426	TCanvas* c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
	427	c1->cd();
	428	gr->Draw("AL");
	429	}
	430
	431
	432	void AliGenHaloProtvino::Copy(TObject&) const
	433	{
	434	//
	435	// Copy
	436	//
	437	Fatal("Copy","Not implemented!\n");
	438	}
	439
	440
	441	/*
	442	# Title: README file for the sources of IR8 machine induced background
	443	# Author: Vadim Talanov <Vadim.Talanov@cern.ch>
	444	# Modified: 12-12-2000
	445
	446	0. Overview
	447
	448	There are three files, named ring.one.beta.[01,10,50].m, which
	449	contain the lists of background particles, induced by proton losses
	450	upstream of IP8 in the LHC ring one, for the beta* values of 1, 10
	451	and 50 m, respectively.
	452
	453	1. File contents
	454
	455	Each line in the files contains the coordinates of particle track
	456	crossing with the infinite plane, positioned at z=-1m, together with
	457	the physical properties of corresponding particle, namely:
	458
	459	S - S coordinate of the primary interaction vertex, cm;
	460	N - type of the gas nuclei at interaction, 1 is H, 2 - C and 3 - O;
	461	I - particle ID in PDG particle numbering scheme;
	462	W - particle weight;
	463	E - particle kinetic energy, GeV;
	464	X - x coordinate of the crossing point, cm;
	465	Y - y coordinate of the crossing point, cm;
	466	Dx - x direction cosine;
	467	Dy - y direction cosine.
	468
	469	2. Normalisation
	470
	471	Each file is given per unity of linear density of proton inelastic
	472	interactions with the gas nuclei, [1 inelastic interaction/m].
	473
	474	# ~/vtalanov/public/README.mib: the end.
	475
	476	*/
	477
	478
	479
	480