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