[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)
{
// Constructor
    
    fName  = "HaloProtvino";
    fTitle = "Halo from LHC Tunnel";
//
//  Read all particles
    fNpart = -1;
    fFile  =  0;
    fSide  =  1;
//
    SetRunPeriod();
    SetTimePerEvent();
    SetAnalog(0);
}

AliGenHaloProtvino::AliGenHaloProtvino(Int_t npart)
    :AliGenerator(npart)
{
// Constructor
    fName = "Halo";
    fTitle= "Halo from LHC Tunnel";
//
    fNpart   = npart;
    fFile    = 0;
    fSide    = 1;
//
    SetRunPeriod();
    SetTimePerEvent();
    SetAnalog(0);
}

AliGenHaloProtvino::AliGenHaloProtvino(const AliGenHaloProtvino & HaloProtvino):
    AliGenerator(HaloProtvino)
{
// Copy constructor
    HaloProtvino.Copy(*this);
}


//____________________________________________________________
AliGenHaloProtvino::~AliGenHaloProtvino()
{
// Destructor
}

//____________________________________________________________
void AliGenHaloProtvino::Init() 
{
// Initialisation
    fFile = fopen(fFileName,"r");
    if (fFile) {
	printf("\n File %s opened for reading, %p ! \n ",  fFileName.Data(), (void*)fFile);
    } else {
	printf("\n Opening of file %s failed,  %p ! \n ",  fFileName.Data(), (void*)fFile);
    }
//
//
//
//    Read file with gas pressure values
    char *name = 0;
    if (fRunPeriod < 5) {
	name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
	fGPASize = 21;
	fG1 = new Float_t[fGPASize];
	fG2 = new Float_t[fGPASize];
	fZ1 = new Float_t[fGPASize];
	fZ2 = new Float_t[fGPASize];
    } else if (fRunPeriod == 5) {
	name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_startup.dat");
	fGPASize = 18853;
	fG1 = new Float_t[fGPASize];
	fZ1 = new Float_t[fGPASize];
    } else if (fRunPeriod ==6) {
	name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_conditioned.dat");
	fGPASize = 12719;
	fG1 = new Float_t[fGPASize];
	fZ1 = new Float_t[fGPASize];
    } else {
	Fatal("Init()", "No gas pressure file for given run period !");
    }


    FILE* file = fopen(name, "r");
    Float_t z;
    Int_t i;
    Float_t p[5];    

    const Float_t kCrossSection = 0.094e-28;      // m^2
    const Float_t kFlux         = 1.e11 / 25.e-9; // protons/s
    Float_t pFlux[5] = {0.2, 0.2, 0.3, 0.3, 1.0};

    if (fRunPeriod < 5) {
//
//  Ring 1   
// 

	for (i = 0; i < fGPASize; i++)
	{
	    fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
	    fG1[i] = p[fRunPeriod];
	    
	    if (i > 0) {
		fZ1[i] = fZ1[i-1] + z;
	    } else {
		fZ1[i] = 20.;
	    }
	}
//
// Ring 2
//
	for (i = 0; i < fGPASize; i++)
	{
	    fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
	    fG2[i] = p[fRunPeriod];
	    if (i > 0) {
		fZ2[i] = fZ2[i-1] + z;
	    } else {
		fZ2[i] = 20.;
	    }
	}
//
// Interaction rates
//
	for (i = 0; i < fGPASize; i++)  
	{
	    fG1[i] = fG1[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s 
	    fG2[i] = fG2[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
	}

    } else {
	for (i = 0; i < fGPASize; i++) 
	{
	    fscanf(file, "%f %e %e %e %e %e", &z, &p[0], &p[1], &p[2], &p[3], &p[4]);
	    z /= 1000.;
	    fG1[i] = p[4] * kCrossSection * kFlux;             // 1/m/s
	    // 1/3 of nominal intensity at startup
	    if (fRunPeriod ==  kLHCPR674Startup) fG1[i] /= 3.;
	    fZ1[i] = z;
	}
    }
    

//
//  Transform into interaction rates
//


    Float_t sum1 = 0.;
    Float_t sum2 = 0.;
    
    for (Int_t i = 0; i < 300; i++) {
	Float_t z = 20.+i*1.;
	z*=100;
	Float_t wgt1 = GasPressureWeight(z);
	Float_t wgt2 = GasPressureWeight(-z);
//	printf("weight: %f %f %f %f %f \n", z, wgt1, wgt2, fZ1[20], fZ2[20]);
	sum1 += wgt1;
	sum2 += wgt2;
    }
    sum1/=250.;
    sum2/=250.;
    printf("\n %f %f \n \n", sum1, sum2);
}

//____________________________________________________________
void AliGenHaloProtvino::Generate()
{
// Generate from input file
 
  Float_t polar[3]= {0,0,0};
  Float_t origin[3];
  Float_t p[3], p0;
  Float_t tz, txy;
  Float_t amass;
  //
  Int_t ncols, nt;
  static Int_t nskip = 0;
  Int_t nread = 0;

  Float_t* zPrimary = new Float_t [fNpart];
  Int_t  * inuc     = new Int_t   [fNpart];
  Int_t  * ipart    = new Int_t   [fNpart];
  Float_t* wgt      = new Float_t [fNpart]; 
  Float_t* ekin     = new Float_t [fNpart];
  Float_t* vx       = new Float_t [fNpart];
  Float_t* vy       = new Float_t [fNpart];
  Float_t* tx       = new Float_t [fNpart];
  Float_t* ty       = new Float_t [fNpart];
  
  Float_t zVertexOld = -1.e10;
  Int_t   nInt       = 0;        // Counts number of interactions
  Float_t wwgt = 0.;
  
  while(1) {
//
// Load event into array
//
      ncols = fscanf(fFile,"%f %d %d %f %f %f %f %f %f",
		     &zPrimary[nread], &inuc[nread], &ipart[nread], &wgt[nread], 
		     &ekin[nread], &vx[nread], &vy[nread],
		     &tx[nread], &ty[nread]);
      
      if (ncols < 0) break;
// Skip fNskip events
      nskip++;
      if (fNpart !=-1 && nskip <= fNskip) continue;
// Count interactions
      if (zPrimary[nread] != zVertexOld) {
	  nInt++;
	  zVertexOld = zPrimary[nread];
      }
// Count tracks      
      nread++;
      if (fNpart !=-1 && nread > fNpart) break;
  }
//
// Mean time between interactions
//
  Float_t dT = fTimePerEvent/nInt;   // sec 
  Float_t t  = 0;                    // sec
  
//
// Loop over primaries
//
  zVertexOld   = -1.e10;
  Double_t arg = 0.;
  
  for (Int_t nprim = 0; nprim < fNpart; nprim++) 
  {
      amass = TDatabasePDG::Instance()->GetParticle(ipart[nprim])->Mass();

      //
      // Momentum vector
      //
      p0=sqrt(ekin[nprim]*ekin[nprim] + 2.*amass*ekin[nprim]);
      
      txy=TMath::Sqrt(tx[nprim]*tx[nprim]+ty[nprim]*ty[nprim]);
      if (txy == 1.) {
	  tz=0;
      } else {
	  tz=-TMath::Sqrt(1.-txy);
      }
    
      p[0] = p0*tx[nprim];
      p[1] = p0*ty[nprim];
      p[2] =-p0*tz;
      
      origin[0] = vx[nprim];
      origin[1] = vy[nprim];
      origin[2] = -2196.5;

      //
      //
      // Particle weight

      Float_t originP[3] = {0., 0., 0.};
      originP[2] = zPrimary[nprim];
      
      Float_t pP[3] = {0., 0., 0.};
      Int_t ntP;
      
      if (fSide == -1) {
	  originP[2] = -zPrimary[nprim];
	  origin[2]  = -origin[2];
	  p[2]       = -p[2];
      }

      //
      // Time
      //
      if (zPrimary[nprim] != zVertexOld) {
	  while(arg==0.) arg = gRandom->Rndm();
	  t -= dT*TMath::Log(arg);              // (sec)   
	  zVertexOld = zPrimary[nprim];
      }

//    Get statistical weight according to local gas-pressure
//
      fParentWeight=wgt[nprim]*GasPressureWeight(zPrimary[nprim]);

      if (!fAnalog || gRandom->Rndm() < fParentWeight) {
//    Pass parent particle
//
	  PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
	  KeepTrack(ntP);
	  PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
      }
      //
      // Both sides are considered
      //

      if (fSide > 1) {
	  fParentWeight=wgt[nprim]*GasPressureWeight(-zPrimary[nprim]);
	  if (!fAnalog || gRandom->Rndm() < fParentWeight) {
	      origin[2]  = -origin[2];
	      originP[2] = -originP[2];
	      p[2]=-p[2];
	      PushTrack(0,-1,kProton,pP,originP,polar,t,kPNoProcess,ntP, fParentWeight);
	      KeepTrack(ntP);
	      PushTrack(fTrackIt,ntP,ipart[nprim],p,origin,polar,t,kPNoProcess,nt,fParentWeight);
	  }
      }
      wwgt += fParentWeight;
      
      SetHighWaterMark(nt);
  }
  delete [] zPrimary;
  delete [] inuc;    
  delete [] ipart;   
  delete [] wgt;     
  delete [] ekin;    
  delete [] vx;      
  delete [] vy;      
  delete [] tx;      
  delete [] ty;     
  printf("Total weight %f\n\n", wwgt);
  
}
 

AliGenHaloProtvino& AliGenHaloProtvino::operator=(const  AliGenHaloProtvino& rhs)
{
// Assignment operator
    rhs.Copy(*this);
    return *this;
}


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


void AliGenHaloProtvino::Copy(TObject&) const
{
    //
    // Copy 
    //
    Fatal("Copy","Not implemented!\n");
}


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