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

    Int_t ncols = 0;
    if (fRunPeriod < 5) {
//
//  Ring 1   
// 

	for (i = 0; i < fGPASize; i++)
	{
	    ncols = fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
	    if (ncols<0) break;

	    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++)
	{
	    ncols = fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
	    if (ncols<0) break;

	    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++) 
	{
	    ncols = fscanf(file, "%f %e %e %e %e %e", &z, &p[0], &p[1], &p[2], &p[3], &p[4]);
	    if (ncols<0) break;

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