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