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