[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), 
     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  = "HaloProtvino";
    fTitle = "Halo from LHC Tunnel";
//
//  Read all particles
    fNpart = -1;
    SetAnalog(0);
}

AliGenHaloProtvino::AliGenHaloProtvino(Int_t 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;
//
    SetAnalog(0);
}

//____________________________________________________________
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 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);
}

//____________________________________________________________
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);
  
}
 

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");
}


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