[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.7  2002/02/21 16:09:45  morsch
Move SetHighwaterMark() after last possible SetTrack()

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