[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 = 0;
    if (name) file = fopen(name, "r");
    if (!file) {
	AliError("No gas pressure file");
	return;
    }

    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);
    delete file;
}

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