[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.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);
      SetHighWaterMark(nt);
      
      //
      // 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);
      }
      
  }
}
 

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