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

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
199	for (Int_t i = 0; i < 300; i++) {
200	Float_t z = 20.+i*1.;
201	z*=100;
202	Float_t wgt1 = GasPressureWeight(z);
203	Float_t wgt2 = GasPressureWeight(-z);
204	// printf("weight: %f %f %f %f %f \n", z, wgt1, wgt2, fZ1[20], fZ2[20]);
205	sum1 += wgt1;
206	sum2 += wgt2;
207	}
208	sum1/=250.;
209	sum2/=250.;
210	printf("\n %f %f \n \n", sum1, sum2);
211	}
212
213	//____________________________________________________________
214	void AliGenHaloProtvino::Generate()
215	{
216	// Generate from input file
217
218	Float_t polar[3]= {0,0,0};
219	Float_t origin[3];
220	Float_t p[3], p0;
221	Float_t tz, txy;
222	Float_t amass;
223	//
224	Int_t ncols, nt;
225	static Int_t nskip = 0;
226	Int_t nread = 0;
227
228	Float_t* zPrimary = new Float_t [fNpart];
229	Int_t * inuc = new Int_t [fNpart];
230	Int_t * ipart = new Int_t [fNpart];
231	Float_t* wgt = new Float_t [fNpart];
232	Float_t* ekin = new Float_t [fNpart];
233	Float_t* vx = new Float_t [fNpart];
234	Float_t* vy = new Float_t [fNpart];
235	Float_t* tx = new Float_t [fNpart];
236	Float_t* ty = new Float_t [fNpart];
237
238	Float_t zVertexOld = -1.e10;
239	Int_t nInt = 0; // Counts number of interactions
240	Float_t wwgt = 0.;
241
242	while(1) {
243	//
244	// Load event into array
245	//
246	ncols = fscanf(fFile,"%f %d %d %f %f %f %f %f %f",
247	&zPrimary[nread], &inuc[nread], &ipart[nread], &wgt[nread],
248	&ekin[nread], &vx[nread], &vy[nread],
249	&tx[nread], &ty[nread]);
250
251	if (ncols < 0) break;
252	// Skip fNskip events
253	nskip++;
254	if (fNpart !=-1 && nskip <= fNskip) continue;
255	// Count interactions
256	if (zPrimary[nread] != zVertexOld) {
257	nInt++;
258	zVertexOld = zPrimary[nread];
259	}
260	// Count tracks
261	nread++;
262	if (fNpart !=-1 && nread > fNpart) break;
263	}
264	//
265	// Mean time between interactions
266	//
267	Float_t dT = fTimePerEvent/nInt; // sec
268	Float_t t = 0; // sec
269
270	//
271	// Loop over primaries
272	//
273	zVertexOld = -1.e10;
274	Double_t arg = 0.;
275
276	for (Int_t nprim = 0; nprim < fNpart; nprim++)
277	{
278	amass = TDatabasePDG::Instance()->GetParticle(ipart[nprim])->Mass();
279
280	//
281	// Momentum vector
282	//
283	p0=sqrt(ekin[nprim]ekin[nprim] + 2.amass*ekin[nprim]);
284
285	txy=TMath::Sqrt(tx[nprim]tx[nprim]+ty[nprim]ty[nprim]);
286	if (txy == 1.) {
287	tz=0;
288	} else {
289	tz=-TMath::Sqrt(1.-txy);
290	}
291
292	p[0] = p0*tx[nprim];
293	p[1] = p0*ty[nprim];
294	p[2] =-p0*tz;
295
296	origin[0] = vx[nprim];
297	origin[1] = vy[nprim];
298	origin[2] = -2196.5;
299
300	//
301	//
302	// Particle weight
303
304	Float_t originP[3] = {0., 0., 0.};
305	originP[2] = zPrimary[nprim];
306
307	Float_t pP[3] = {0., 0., 0.};
308	Int_t ntP;
309
310	if (fSide == -1) {
311	originP[2] = -zPrimary[nprim];
312	origin[2] = -origin[2];
313	p[2] = -p[2];
314	}
315
316	//
317	// Time
318	//
319	if (zPrimary[nprim] != zVertexOld) {
320	while(arg==0.) arg = gRandom->Rndm();
321	t -= dT*TMath::Log(arg); // (sec)
322	zVertexOld = zPrimary[nprim];
323	}
324
325	// Get statistical weight according to local gas-pressure
326	//
327	fParentWeight=wgt[nprim]*GasPressureWeight(zPrimary[nprim]);
328
329	if (!fAnalog \|\| gRandom->Rndm() < fParentWeight) {
330	// Pass parent particle
331	//
332	PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
333	KeepTrack(ntP);
334	PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
335	}
336	//
337	// Both sides are considered
338	//
339
340	if (fSide > 1) {
341	fParentWeight=wgt[nprim]*GasPressureWeight(-zPrimary[nprim]);
342	if (!fAnalog \|\| gRandom->Rndm() < fParentWeight) {
343	origin[2] = -origin[2];
344	originP[2] = -originP[2];
345	p[2]=-p[2];
346	PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
347	KeepTrack(ntP);
348	PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
349	}
350	}
351	wwgt += fParentWeight;
352
353	SetHighWaterMark(nt);
354	}
355	delete [] zPrimary;
356	delete [] inuc;
357	delete [] ipart;
358	delete [] wgt;
359	delete [] ekin;
360	delete [] vx;
361	delete [] vy;
362	delete [] tx;
363	delete [] ty;
364	printf("Total weight %f\n\n", wwgt);
365
366	}
367
368
369	Float_t AliGenHaloProtvino::GasPressureWeight(Float_t zPrimary)
370	{
371	//
372	// Return z-dependent gasspressure weight = interaction rate [1/m/s].
373	//
374	Float_t weight = 0.;
375	zPrimary /= 100.; // m
376	if (fRunPeriod < 5) {
377	Float_t zAbs = TMath::Abs(zPrimary);
378	if (zPrimary > 0.)
379	{
380	if (zAbs > fZ1[20]) {
381	weight = 2.e4;
382	} else {
383	for (Int_t i = 1; i < 21; i++) {
384	if (zAbs < fZ1[i]) {
385	weight = fG1[i];
386	break;
387	}
388	}
389	}
390	} else {
391	if (zAbs > fZ2[20]) {
392	weight = 2.e4;
393	} else {
394	for (Int_t i = 1; i < 21; i++) {
395	if (zAbs < fZ2[i]) {
396	weight = fG2[i];
397	break;
398	}
399	}
400	}
401	}
402	} else {
403	Int_t index = TMath::BinarySearch(fGPASize, fZ1, zPrimary);
404	weight = fG1[index];
405	}
406	return weight;
407	}
408
409	void AliGenHaloProtvino::Draw(Option_t *)
410	{
411	// Draws the gas pressure distribution
412	Float_t z[400];
413	Float_t p[400];
414
415	for (Int_t i = 0; i < 400; i++)
416	{
417	z[i] = -20000. + Float_t(i) * 100;
418	p[i] = GasPressureWeight(z[i]);
419	}
420
421	TGraph* gr = new TGraph(400, z, p);
422	TCanvas* c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
423	c1->cd();
424	gr->Draw("AL");
425	}
426
198bb1c7	427
14ab5373	428	/*
	429	# Title: README file for the sources of IR8 machine induced background
	430	# Author: Vadim Talanov <Vadim.Talanov@cern.ch>
	431	# Modified: 12-12-2000
	432
	433	0. Overview
	434
	435	There are three files, named ring.one.beta.[01,10,50].m, which
	436	contain the lists of background particles, induced by proton losses
	437	upstream of IP8 in the LHC ring one, for the beta* values of 1, 10
	438	and 50 m, respectively.
	439
	440	1. File contents
	441
	442	Each line in the files contains the coordinates of particle track
	443	crossing with the infinite plane, positioned at z=-1m, together with
	444	the physical properties of corresponding particle, namely:
	445
	446	S - S coordinate of the primary interaction vertex, cm;
	447	N - type of the gas nuclei at interaction, 1 is H, 2 - C and 3 - O;
	448	I - particle ID in PDG particle numbering scheme;
	449	W - particle weight;
	450	E - particle kinetic energy, GeV;
	451	X - x coordinate of the crossing point, cm;
	452	Y - y coordinate of the crossing point, cm;
	453	Dx - x direction cosine;
	454	Dy - y direction cosine.
	455
	456	2. Normalisation
	457
	458	Each file is given per unity of linear density of proton inelastic
	459	interactions with the gas nuclei, [1 inelastic interaction/m].
	460
	461	# ~/vtalanov/public/README.mib: the end.
	462
	463	*/
	464
	465
	466
	467