[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.                  *
 **************************************************************************/

/*
$Log$
Revision 1.6  2002/02/01 15:28:22  morsch
Fastidious print statements removed.

Revision 1.5  2002/02/01 14:12:28  morsch
Include new gas pressure estimates (LHC Pproject Note 273)

Revision 1.4  2001/12/19 16:30:24  morsch
Some bugs corrected and skip method added. (Rachid Guernane)

Revision 1.3  2001/07/27 17:09:36  morsch
Use local SetTrack, KeepTrack and SetHighWaterMark methods
to delegate either to local stack or to stack owned by AliRun.
(Piotr Skowronski, A.M.)

Revision 1.2  2001/06/14 12:15:27  morsch
Bugs corrected. SetSide() method added.

Revision 1.1  2001/01/23 15:04:33  morsch
Generator to read beam halo file from Protvino group.

*/

// 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 "AliGenHaloProtvino.h"
#include "AliRun.h"
#include "AliMC.h"
#include "AliPDG.h"

#include <TDatabasePDG.h>
#include <TSystem.h>
#include <stdlib.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();
}

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

AliGenHaloProtvino::AliGenHaloProtvino(const AliGenHaloProtvino & HaloProtvino)
{
// copy constructor
}


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

//____________________________________________________________
void AliGenHaloProtvino::Init() 
{
// Initialisation
    fFile = fopen(fFileName,"r");
    if (fFile) {
	printf("\n File %s opened for reading, %p ! \n ",  fFileName.Data(), fFile);
    } else {
	printf("\n Opening of file %s failed,  %p ! \n ",  fFileName.Data(), fFile);
    }
//
//
//
//    Read file with gas pressure values
    char *name;
    
    name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
    FILE* file = fopen(name, "r");
    Float_t z;
    Int_t i, j;
    
//
//  Ring 1   
// 
    for (i = 0; i < 20; i++)
    {
	fscanf(file, "%f %f %f %f %f %f", &z, 
	       &fG1[i][0], &fG1[i][1], &fG1[i][2], &fG1[i][3], &fG1[i][4]);
	if (i > 0) {
	    fZ1[i] = fZ1[i-1] + z;
	} else {
	    fZ1[i] = 20.;
	}
  }
//
// Ring 2
//
    for (i = 0; i < 21; i++)
    {
	fscanf(file, "%f %f %f %f %f %f", &z, 
	       &fG2[i][0], &fG2[i][1], &fG2[i][2], &fG2[i][3], &fG2[i][4]);
	if (i > 0) {
	    fZ2[i] = fZ2[i-1] + z;
	} else {
	    fZ2[i] = 20.;
	}
    }
//
//  Transform into interaction rates
//
    const Float_t crossSection = 0.094e-28;     // m^2
    const Float_t pFlux        = 1.e11/25.e-9;  // 1/s

    for (j = 0; j <  5; j++) {
	for (i = 0; i < 21; i++)  
	{
	    if (i <20)
		fG1[i][j] = fG1[i][j] * crossSection * pFlux; // 1/m/s 
	        fG2[i][j] = fG2[i][j] * crossSection * pFlux; // 1/m/s 
	}
    }


    Float_t sum;
    
    for (Int_t i = 0; i < 250; i++) {
	Float_t z = 20.+i*1.;
	z*=100;
	Float_t wgt = GassPressureWeight(z);
//	printf("\n %f %f", z, wgt);
	sum+=wgt;
    }
    sum/=250.;
    printf("\n %f \n \n", sum);
}

//____________________________________________________________
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 ekin, wgt, tx, ty, tz, txy;
  Float_t zPrimary;
  Float_t amass;
  Int_t   inuc;
  //
  Int_t ipart, ncols, nt;
  Int_t nskip = 0;
  Int_t nread = 0;
  while(1) {
      ncols = fscanf(fFile,"%f %d %d %f %f %f %f %f %f",
		     &zPrimary, &inuc, &ipart, &wgt, 
		     &ekin, &origin[0], &origin[1],
		     &tx, &ty);
/*
      printf(" \n %f %d %d %f %f %f %f %f %f",
             zPrimary, inuc, ipart, wgt, 
             ekin, origin[0], origin[1],
             tx, ty);
*/
      
      if (ncols < 0) break;

      nskip++;
      if (fNpart !=-1 && nskip <= fNskip) continue;
      
      nread++;
      if (fNpart !=-1 && nread > fNpart) break;

      amass = TDatabasePDG::Instance()->GetParticle(ipart)->Mass();

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

      //
      //
      // Particle weight

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

      SetTrack(0,-1,kProton,pP,originP,polar,0,kPNoProcess,ntP);
      KeepTrack(ntP);
      fParentWeight=wgt*GassPressureWeight(zPrimary);
      SetTrack(fTrackIt,ntP,ipart,p,origin,polar,0,kPNoProcess,nt,fParentWeight);
      //
      // Assume particles come from two directions with same probability

      // Both Side are considered
      if (fSide > 1) {
          origin[2]=-origin[2];
          p[2]=-p[2];
          fParentWeight=wgt*GassPressureWeight(-zPrimary);
          SetTrack(fTrackIt,ntP,ipart,p,origin,polar,0,kPNoProcess,nt,fParentWeight);
      }
      SetHighWaterMark(nt);
  }
}
 

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


Float_t AliGenHaloProtvino::GassPressureWeight(Float_t zPrimary)
{
//
// Return z-dependent gasspressure weight = interaction rate [1/m/s].
//
    Float_t weight = 0.;
    zPrimary/=100.;        // m
    Float_t zAbs = TMath::Abs(zPrimary);
    zPrimary = TMath::Abs(zPrimary);
    
    if (zPrimary < 0.) 
    {
	if (zAbs > fZ1[19]) {
	    weight = 2.e4;
	} else {
	    for (Int_t i = 1; i < 20; i++) {
		if (zAbs < fZ1[i]) {
		    weight = fG1[i][fRunPeriod];
		    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][fRunPeriod];
		    break;
		}
	    }
	}
    }
    
    return weight;
}

/*
# 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
	16	/*
	17	$Log$
5bfdba97	18	Revision 1.6 2002/02/01 15:28:22 morsch
	19	Fastidious print statements removed.
	20
69818e67	21	Revision 1.5 2002/02/01 14:12:28 morsch
	22	Include new gas pressure estimates (LHC Pproject Note 273)
	23
eac2c8a3	24	Revision 1.4 2001/12/19 16:30:24 morsch
	25	Some bugs corrected and skip method added. (Rachid Guernane)
	26
05df024b	27	Revision 1.3 2001/07/27 17:09:36 morsch
	28	Use local SetTrack, KeepTrack and SetHighWaterMark methods
	29	to delegate either to local stack or to stack owned by AliRun.
	30	(Piotr Skowronski, A.M.)
	31
a99cf51f	32	Revision 1.2 2001/06/14 12:15:27 morsch
	33	Bugs corrected. SetSide() method added.
	34
154d524b	35	Revision 1.1 2001/01/23 15:04:33 morsch
	36	Generator to read beam halo file from Protvino group.
	37
14ab5373	38	*/
	39
	40	// Read background particles from a boundary source
	41	// Very specialized generator to simulate background from beam halo.
	42	// The input file is a text file specially prepared
	43	// for this purpose.
	44	// Author: andreas.morsch@cern.ch
	45
	46	#include "AliGenHaloProtvino.h"
	47	#include "AliRun.h"
	48	#include "AliMC.h"
	49	#include "AliPDG.h"
	50
	51	#include <TDatabasePDG.h>
eac2c8a3	52	#include <TSystem.h>
14ab5373	53	#include <stdlib.h>
	54
	55	ClassImp(AliGenHaloProtvino)
	56	AliGenHaloProtvino::AliGenHaloProtvino()
	57	:AliGenerator(-1)
	58	{
	59	// Constructor
154d524b	60
	61	fName = "HaloProtvino";
	62	fTitle = "Halo from LHC Tunnel";
14ab5373	63	//
14ab5373	64	// Read all particles
154d524b	65	fNpart = -1;
	66	fFile = 0;
	67	fSide = 1;
eac2c8a3	68	//
eac2c8a3	69	SetRunPeriod();
14ab5373	70	}
	71
	72	AliGenHaloProtvino::AliGenHaloProtvino(Int_t npart)
	73	:AliGenerator(npart)
	74	{
	75	// Constructor
154d524b	76	fName = "Halo";
154d524b	77	fTitle= "Halo from LHC Tunnel";
14ab5373	78	//
154d524b	79	fNpart = npart;
	80	fFile = 0;
	81	fSide = 1;
eac2c8a3	82	//
eac2c8a3	83	SetRunPeriod();
14ab5373	84	}
	85
	86	AliGenHaloProtvino::AliGenHaloProtvino(const AliGenHaloProtvino & HaloProtvino)
	87	{
	88	// copy constructor
	89	}
	90
	91
	92	//____________________________________________________________
	93	AliGenHaloProtvino::~AliGenHaloProtvino()
	94	{
	95	// Destructor
	96	}
	97
	98	//____________________________________________________________
	99	void AliGenHaloProtvino::Init()
	100	{
	101	// Initialisation
154d524b	102	fFile = fopen(fFileName,"r");
	103	if (fFile) {
	104	printf("\n File %s opened for reading, %p ! \n ", fFileName.Data(), fFile);
	105	} else {
	106	printf("\n Opening of file %s failed, %p ! \n ", fFileName.Data(), fFile);
	107	}
eac2c8a3	108	//
	109	//
	110	//
	111	// Read file with gas pressure values
	112	char *name;
	113
	114	name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
	115	FILE* file = fopen(name, "r");
	116	Float_t z;
	117	Int_t i, j;
	118
	119	//
	120	// Ring 1
	121	//
	122	for (i = 0; i < 20; i++)
	123	{
	124	fscanf(file, "%f %f %f %f %f %f", &z,
	125	&fG1[i][0], &fG1[i][1], &fG1[i][2], &fG1[i][3], &fG1[i][4]);
	126	if (i > 0) {
	127	fZ1[i] = fZ1[i-1] + z;
	128	} else {
	129	fZ1[i] = 20.;
	130	}
	131	}
	132	//
	133	// Ring 2
	134	//
	135	for (i = 0; i < 21; i++)
	136	{
	137	fscanf(file, "%f %f %f %f %f %f", &z,
	138	&fG2[i][0], &fG2[i][1], &fG2[i][2], &fG2[i][3], &fG2[i][4]);
	139	if (i > 0) {
	140	fZ2[i] = fZ2[i-1] + z;
	141	} else {
	142	fZ2[i] = 20.;
	143	}
	144	}
	145	//
	146	// Transform into interaction rates
	147	//
	148	const Float_t crossSection = 0.094e-28; // m^2
	149	const Float_t pFlux = 1.e11/25.e-9; // 1/s
	150
	151	for (j = 0; j < 5; j++) {
	152	for (i = 0; i < 21; i++)
	153	{
	154	if (i <20)
	155	fG1[i][j] = fG1[i][j] * crossSection * pFlux; // 1/m/s
	156	fG2[i][j] = fG2[i][j] * crossSection * pFlux; // 1/m/s
	157	}
	158	}
	159
	160
	161	Float_t sum;
	162
	163	for (Int_t i = 0; i < 250; i++) {
	164	Float_t z = 20.+i*1.;
	165	z*=100;
	166	Float_t wgt = GassPressureWeight(z);
69818e67	167	// printf("\n %f %f", z, wgt);
eac2c8a3	168	sum+=wgt;
	169	}
	170	sum/=250.;
	171	printf("\n %f \n \n", sum);
14ab5373	172	}
	173
	174	//____________________________________________________________
	175	void AliGenHaloProtvino::Generate()
	176	{
	177	// Generate from input file
14ab5373	178
	179	Float_t polar[3]= {0,0,0};
	180	Float_t origin[3];
	181	Float_t p[3], p0;
	182	Float_t ekin, wgt, tx, ty, tz, txy;
	183	Float_t zPrimary;
	184	Float_t amass;
	185	Int_t inuc;
	186	//
	187	Int_t ipart, ncols, nt;
05df024b	188	Int_t nskip = 0;
154d524b	189	Int_t nread = 0;
14ab5373	190	while(1) {
154d524b	191	ncols = fscanf(fFile,"%f %d %d %f %f %f %f %f %f",
14ab5373	192	&zPrimary, &inuc, &ipart, &wgt,
	193	&ekin, &origin[0], &origin[1],
	194	&tx, &ty);
154d524b	195	/*
14ab5373	196	printf(" \n %f %d %d %f %f %f %f %f %f",
05df024b	197	zPrimary, inuc, ipart, wgt,
	198	ekin, origin[0], origin[1],
	199	tx, ty);
154d524b	200	*/
05df024b	201
14ab5373	202	if (ncols < 0) break;
154d524b	203
05df024b	204	nskip++;
	205	if (fNpart !=-1 && nskip <= fNskip) continue;
	206
	207	nread++;
14ab5373	208	if (fNpart !=-1 && nread > fNpart) break;
14ab5373	209
14ab5373	210	amass = TDatabasePDG::Instance()->GetParticle(ipart)->Mass();
	211
	212	//
	213	// Momentum vector
	214	//
05df024b	215	p0=sqrt(ekinekin + 2.amass*ekin);
14ab5373	216
	217	txy=TMath::Sqrt(txtx+tyty);
	218	if (txy == 1.) {
	219	tz=0;
	220	} else {
	221	tz=-TMath::Sqrt(1.-txy);
	222	}
	223
	224	p[0]=p0*tx;
	225	p[1]=p0*ty;
154d524b	226	p[2]=-p0*tz;
05df024b	227
05df024b	228	origin[2] = -2196.5;
14ab5373	229
	230	//
	231	//
	232	// Particle weight
	233
154d524b	234	Float_t originP[3] = {0., 0., 0.};
	235	originP[2] = zPrimary;
	236
	237	Float_t pP[3] = {0., 0., 0.};
	238	Int_t ntP;
	239
	240	if (fSide == -1) {
	241	originP[2] = -zPrimary;
	242	origin[2] = -origin[2];
	243	p[2] = -p[2];
	244	}
05df024b	245
a99cf51f	246	SetTrack(0,-1,kProton,pP,originP,polar,0,kPNoProcess,ntP);
a99cf51f	247	KeepTrack(ntP);
14ab5373	248	fParentWeight=wgt*GassPressureWeight(zPrimary);
a99cf51f	249	SetTrack(fTrackIt,ntP,ipart,p,origin,polar,0,kPNoProcess,nt,fParentWeight);
14ab5373	250	//
14ab5373	251	// Assume particles come from two directions with same probability
154d524b	252
05df024b	253	// Both Side are considered
	254	if (fSide > 1) {
	255	origin[2]=-origin[2];
	256	p[2]=-p[2];
	257	fParentWeight=wgt*GassPressureWeight(-zPrimary);
	258	SetTrack(fTrackIt,ntP,ipart,p,origin,polar,0,kPNoProcess,nt,fParentWeight);
	259	}
5bfdba97	260	SetHighWaterMark(nt);
14ab5373	261	}
	262	}
	263
	264
	265	AliGenHaloProtvino& AliGenHaloProtvino::operator=(const AliGenHaloProtvino& rhs)
	266	{
	267	// Assignment operator
	268	return *this;
	269	}
	270
	271
eac2c8a3	272
14ab5373	273	Float_t AliGenHaloProtvino::GassPressureWeight(Float_t zPrimary)
14ab5373	274	{
eac2c8a3	275	//
	276	// Return z-dependent gasspressure weight = interaction rate [1/m/s].
	277	//
	278	Float_t weight = 0.;
	279	zPrimary/=100.; // m
	280	Float_t zAbs = TMath::Abs(zPrimary);
	281	zPrimary = TMath::Abs(zPrimary);
154d524b	282
eac2c8a3	283	if (zPrimary < 0.)
	284	{
	285	if (zAbs > fZ1[19]) {
	286	weight = 2.e4;
	287	} else {
	288	for (Int_t i = 1; i < 20; i++) {
	289	if (zAbs < fZ1[i]) {
	290	weight = fG1[i][fRunPeriod];
	291	break;
	292	}
	293	}
	294	}
	295	} else {
	296	if (zAbs > fZ2[20]) {
	297	weight = 2.e4;
	298	} else {
	299	for (Int_t i = 1; i < 21; i++) {
	300	if (zAbs < fZ2[i]) {
	301	weight = fG2[i][fRunPeriod];
	302	break;
	303	}
	304	}
	305	}
	306	}
154d524b	307
eac2c8a3	308	return weight;
14ab5373	309	}
	310
	311	/*
	312	# Title: README file for the sources of IR8 machine induced background
	313	# Author: Vadim Talanov <Vadim.Talanov@cern.ch>
	314	# Modified: 12-12-2000
	315
	316	0. Overview
	317
	318	There are three files, named ring.one.beta.[01,10,50].m, which
	319	contain the lists of background particles, induced by proton losses
	320	upstream of IP8 in the LHC ring one, for the beta* values of 1, 10
	321	and 50 m, respectively.
	322
	323	1. File contents
	324
	325	Each line in the files contains the coordinates of particle track
	326	crossing with the infinite plane, positioned at z=-1m, together with
	327	the physical properties of corresponding particle, namely:
	328
	329	S - S coordinate of the primary interaction vertex, cm;
	330	N - type of the gas nuclei at interaction, 1 is H, 2 - C and 3 - O;
	331	I - particle ID in PDG particle numbering scheme;
	332	W - particle weight;
	333	E - particle kinetic energy, GeV;
	334	X - x coordinate of the crossing point, cm;
	335	Y - y coordinate of the crossing point, cm;
	336	Dx - x direction cosine;
	337	Dy - y direction cosine.
	338
	339	2. Normalisation
	340
	341	Each file is given per unity of linear density of proton inelastic
	342	interactions with the gas nuclei, [1 inelastic interaction/m].
	343
	344	# ~/vtalanov/public/README.mib: the end.
	345
	346	*/
	347
	348
	349
	350