[u/mrichter/AliRoot.git] / EVGEN / AliGenHalo.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 beam halo background particles from a boundary source
// Boundary source is in the LHCb format http://www.hep.manchester.ac.uk/u/robert/LHC_backgrounds/Note-MIBStudies.pdf
// and has been provided by Robert Appleby
// Author: andreas.morsch@cern.ch


#include <stdlib.h>

#include <TDatabasePDG.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TPDGCode.h>
#include <TSystem.h>

#include "AliGenHalo.h"
#include "AliRun.h"

ClassImp(AliGenHalo)

AliGenHalo::AliGenHalo()
    :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),
     fLossID(0),   
     fLossA(0),   
     fPdg(0),     
     fLossT0(0),
     fLossZ(0), 
     fLossW(0), 
     fXS(0),    
     fYS(0),    
     fZS(0),    
     fDX(0),    
     fDY(0),    
     fEkin(0),  
     fTS(0),    
     fWS(0)     
{
// Constructor
    
    fName  = "Halo";
    fTitle = "Halo from LHC Beam";
//
//  Read all particles
    fNpart = -1;
    SetAnalog(0);
}

AliGenHalo::AliGenHalo(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),
     fLossID(0),   
     fLossA(0),   
     fPdg(0),     
     fLossT0(0),
     fLossZ(0), 
     fLossW(0), 
     fXS(0),    
     fYS(0),    
     fZS(0),    
     fDX(0),    
     fDY(0),    
     fEkin(0),  
     fTS(0),    
     fWS(0)     
{
// Constructor
    fName = "Halo";
    fTitle= "Halo from LHC Beam";
//
    fNpart   = npart;
//
    SetAnalog(0);
}

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

//____________________________________________________________
void AliGenHalo::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 iz = 0; iz < 300; iz++) {
	Float_t zpos = 20. + iz * 1.;
	zpos *= 100;
	Float_t wgt1 = GasPressureWeight( zpos);
	Float_t wgt2 = GasPressureWeight(-zpos);
	sum1 += wgt1;
	sum2 += wgt2;
    }
    sum1/=250.;
    sum2/=250.;
    printf("\n %f %f \n \n", sum1, sum2);
}

//____________________________________________________________
void AliGenHalo::Generate()
{
// Generate by reading particles 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 mass;
  //
  Int_t nt;
  static Bool_t first = kTRUE;
  static Int_t  oldID = -1;
//
  if (first) {
      ReadNextParticle();
      first = kFALSE;
      oldID = fLossID;
  }
  
  
  while(1) {
      // Push particle to stack
      mass = TDatabasePDG::Instance()->GetParticle(fPdg)->Mass();
      p0  = TMath::Sqrt(fEkin * fEkin + 2.* mass * fEkin);
      txy = TMath::Sqrt(fDX * fDX + fDY * fDY);
      if (txy == 1.) {
	  tz = 0;
      } else {
	  tz = - TMath::Sqrt(1. - txy);
      }
      p[0] =  p0 * fDX;
      p[1] =  p0 * fDY;
      p[2] =  p0 * tz;

      origin[0] = fXS;
      origin[1] = fYS;
      origin[2] = 1950.;

      PushTrack(fTrackIt , -1, fPdg , p, origin, polar, fTS, kPNoProcess, nt, fWS);

      Int_t nc = ReadNextParticle();
      
      if (fLossID != oldID || nc == 0) {
	  oldID = fLossID;
	  break;
      }
  }
  SetHighWaterMark(nt);
}
 

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

Int_t AliGenHalo::ReadNextParticle()
{
    // Read the next particle from the file
    Int_t ncols = fscanf(fFile,"%d %f %f %d %f %d %f %f %f %f %f %f %f %f",
		   &fLossID, &fLossT0, &fLossZ, &fLossA, &fLossW, &fPdg, &fXS, &fYS, &fZS, &fDX, &fDY, &fEkin, &fTS, &fWS);
    fLossZ /= 10.;
    fXS    /= 10.;
    fYS    /= 10.; 
    fZS    /= 10.;   
    fTS    *= 1.e-9;
    return (ncols);
}
Commit	Line	Data
4c039060	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
803d1ab0	16	/* $Id$ */
4c039060	17
d12ebe74	18	// Read beam halo background particles from a boundary source
	19	// Boundary source is in the LHCb format http://www.hep.manchester.ac.uk/u/robert/LHC_backgrounds/Note-MIBStudies.pdf
	20	// and has been provided by Robert Appleby
675e9664	21	// Author: andreas.morsch@cern.ch
675e9664	22
d12ebe74	23
116cbefd	24	#include <stdlib.h>
c5970876	25
1578254f	26	#include <TDatabasePDG.h>
d12ebe74	27	#include <TCanvas.h>
d12ebe74	28	#include <TGraph.h>
116cbefd	29	#include <TPDGCode.h>
d12ebe74	30	#include <TSystem.h>
116cbefd	31
	32	#include "AliGenHalo.h"
	33	#include "AliRun.h"
f87cfe57	34
198bb1c7	35	ClassImp(AliGenHalo)
	36
	37	AliGenHalo::AliGenHalo()
d12ebe74	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),
	50	fLossID(0),
	51	fLossA(0),
	52	fPdg(0),
	53	fLossT0(0),
	54	fLossZ(0),
	55	fLossW(0),
	56	fXS(0),
	57	fYS(0),
	58	fZS(0),
	59	fDX(0),
	60	fDY(0),
	61	fEkin(0),
	62	fTS(0),
	63	fWS(0)
fe4da5cc	64	{
f87cfe57	65	// Constructor
d12ebe74	66
	67	fName = "Halo";
	68	fTitle = "Halo from LHC Beam";
fe4da5cc	69	//
fe4da5cc	70	// Read all particles
d12ebe74	71	fNpart = -1;
d12ebe74	72	SetAnalog(0);
fe4da5cc	73	}
	74
	75	AliGenHalo::AliGenHalo(Int_t npart)
1c56e311	76	:AliGenerator(npart),
d12ebe74	77	fFile(0),
	78	fFileName(0),
	79	fSide(1),
	80	fRunPeriod(kY3D90),
	81	fTimePerEvent(1.e-4),
	82	fNskip(0),
	83	fZ1(0),
	84	fZ2(0),
	85	fG1(0),
	86	fG2(0),
	87	fGPASize(0),
	88	fLossID(0),
	89	fLossA(0),
	90	fPdg(0),
	91	fLossT0(0),
	92	fLossZ(0),
	93	fLossW(0),
	94	fXS(0),
	95	fYS(0),
	96	fZS(0),
	97	fDX(0),
	98	fDY(0),
	99	fEkin(0),
	100	fTS(0),
	101	fWS(0)
fe4da5cc	102	{
f87cfe57	103	// Constructor
d12ebe74	104	fName = "Halo";
d12ebe74	105	fTitle= "Halo from LHC Beam";
fe4da5cc	106	//
d12ebe74	107	fNpart = npart;
	108	//
	109	SetAnalog(0);
fe4da5cc	110	}
	111
	112	//____________________________________________________________
	113	AliGenHalo::~AliGenHalo()
	114	{
f87cfe57	115	// Destructor
fe4da5cc	116	}
	117
	118	//____________________________________________________________
	119	void AliGenHalo::Init()
f87cfe57	120	{
f87cfe57	121	// Initialisation
d12ebe74	122	fFile = fopen(fFileName,"r");
	123	if (fFile) {
	124	printf("\n File %s opened for reading, %p ! \n ", fFileName.Data(), (void*)fFile);
	125	} else {
	126	printf("\n Opening of file %s failed, %p ! \n ", fFileName.Data(), (void*)fFile);
	127	}
	128	//
	129	//
	130	//
	131	// Read file with gas pressure values
	132	char *name = 0;
	133	if (fRunPeriod < 5) {
	134	name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
	135	fGPASize = 21;
	136	fG1 = new Float_t[fGPASize];
	137	fG2 = new Float_t[fGPASize];
	138	fZ1 = new Float_t[fGPASize];
	139	fZ2 = new Float_t[fGPASize];
	140	} else if (fRunPeriod == 5) {
	141	name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_startup.dat");
	142	fGPASize = 18853;
	143	fG1 = new Float_t[fGPASize];
	144	fZ1 = new Float_t[fGPASize];
	145	} else if (fRunPeriod ==6) {
	146	name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_conditioned.dat");
	147	fGPASize = 12719;
	148	fG1 = new Float_t[fGPASize];
	149	fZ1 = new Float_t[fGPASize];
	150	} else {
	151	Fatal("Init()", "No gas pressure file for given run period !");
	152	}
	153
	154
	155	FILE* file = fopen(name, "r");
	156	Float_t z;
	157	Int_t i;
	158	Float_t p[5];
	159
	160	const Float_t kCrossSection = 0.094e-28; // m^2
	161	const Float_t kFlux = 1.e11 / 25.e-9; // protons/s
	162	Float_t pFlux[5] = {0.2, 0.2, 0.3, 0.3, 1.0};
	163
	164	if (fRunPeriod < 5) {
	165	//
	166	// Ring 1
	167	//
	168
	169	for (i = 0; i < fGPASize; i++)
	170	{
	171	fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
	172	fG1[i] = p[fRunPeriod];
	173
	174	if (i > 0) {
	175	fZ1[i] = fZ1[i-1] + z;
	176	} else {
	177	fZ1[i] = 20.;
	178	}
	179	}
	180	//
	181	// Ring 2
	182	//
	183	for (i = 0; i < fGPASize; i++)
	184	{
	185	fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
186	fG2[i] = p[fRunPeriod];
187	if (i > 0) {
188	fZ2[i] = fZ2[i-1] + z;
189	} else {
190	fZ2[i] = 20.;
191	}
192	}
193	//
194	// Interaction rates
195	//
196	for (i = 0; i < fGPASize; i++)
197	{
198	fG1[i] = fG1[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
199	fG2[i] = fG2[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
200	}
201
202	} else {
203	for (i = 0; i < fGPASize; i++)
204	{
205	fscanf(file, "%f %e %e %e %e %e", &z, &p[0], &p[1], &p[2], &p[3], &p[4]);
206	z /= 1000.;
207	fG1[i] = p[4] * kCrossSection * kFlux; // 1/m/s
208	// 1/3 of nominal intensity at startup
209	if (fRunPeriod == kLHCPR674Startup) fG1[i] /= 3.;
210	fZ1[i] = z;
211	}
212	}
213
214
215
216
217	//
218	// Transform into interaction rates
219	//
220
221
222
223
224	Float_t sum1 = 0.;
225	Float_t sum2 = 0.;
226
227	for (Int_t iz = 0; iz < 300; iz++) {
228	Float_t zpos = 20. + iz * 1.;
229	zpos *= 100;
230	Float_t wgt1 = GasPressureWeight( zpos);
231	Float_t wgt2 = GasPressureWeight(-zpos);
232	sum1 += wgt1;
233	sum2 += wgt2;
234	}
235	sum1/=250.;
236	sum2/=250.;
237	printf("\n %f %f \n \n", sum1, sum2);
f87cfe57	238	}
fe4da5cc	239
	240	//____________________________________________________________
	241	void AliGenHalo::Generate()
	242	{
d12ebe74	243	// Generate by reading particles from input file
fe4da5cc	244
d12ebe74	245	Float_t polar[3]= {0., 0., 0.};
fe4da5cc	246	Float_t origin[3];
fe4da5cc	247	Float_t p[3], p0;
d12ebe74	248	Float_t tz, txy;
d12ebe74	249	Float_t mass;
fe4da5cc	250	//
d12ebe74	251	Int_t nt;
	252	static Bool_t first = kTRUE;
	253	static Int_t oldID = -1;
	254	//
	255	if (first) {
	256	ReadNextParticle();
	257	first = kFALSE;
	258	oldID = fLossID;
	259	}
	260
fe4da5cc	261
fe4da5cc	262
fe4da5cc	263	while(1) {
d12ebe74	264	// Push particle to stack
	265	mass = TDatabasePDG::Instance()->GetParticle(fPdg)->Mass();
	266	p0 = TMath::Sqrt(fEkin * fEkin + 2.* mass * fEkin);
	267	txy = TMath::Sqrt(fDX * fDX + fDY * fDY);
fe4da5cc	268	if (txy == 1.) {
d12ebe74	269	tz = 0;
fe4da5cc	270	} else {
d12ebe74	271	tz = - TMath::Sqrt(1. - txy);
fe4da5cc	272	}
d12ebe74	273	p[0] = p0 * fDX;
	274	p[1] = p0 * fDY;
	275	p[2] = p0 * tz;
fe4da5cc	276
d12ebe74	277	origin[0] = fXS;
	278	origin[1] = fYS;
	279	origin[2] = 1950.;
fe4da5cc	280
d12ebe74	281	PushTrack(fTrackIt , -1, fPdg , p, origin, polar, fTS, kPNoProcess, nt, fWS);
fe4da5cc	282
d12ebe74	283	Int_t nc = ReadNextParticle();
	284
	285	if (fLossID != oldID \|\| nc == 0) {
	286	oldID = fLossID;
	287	break;
	288	}
	289	}
	290	SetHighWaterMark(nt);
	291	}
	292
fe4da5cc	293
d12ebe74	294	Float_t AliGenHalo::GasPressureWeight(Float_t zPrimary)
	295	{
	296	//
	297	// Return z-dependent gasspressure weight = interaction rate [1/m/s].
	298	//
	299	Float_t weight = 0.;
	300	zPrimary /= 100.; // m
	301	if (fRunPeriod < 5) {
	302	Float_t zAbs = TMath::Abs(zPrimary);
	303	if (zPrimary > 0.)
	304	{
	305	if (zAbs > fZ1[20]) {
	306	weight = 2.e4;
	307	} else {
	308	for (Int_t i = 1; i < 21; i++) {
	309	if (zAbs < fZ1[i]) {
	310	weight = fG1[i];
	311	break;
	312	}
	313	}
	314	}
	315	} else {
	316	if (zAbs > fZ2[20]) {
	317	weight = 2.e4;
	318	} else {
	319	for (Int_t i = 1; i < 21; i++) {
	320	if (zAbs < fZ2[i]) {
	321	weight = fG2[i];
	322	break;
	323	}
	324	}
	325	}
	326	}
	327	} else {
	328	Int_t index = TMath::BinarySearch(fGPASize, fZ1, zPrimary);
	329	weight = fG1[index];
	330	}
	331	return weight;
	332	}
fe4da5cc	333
d12ebe74	334	void AliGenHalo::Draw(Option_t *)
	335	{
	336	// Draws the gas pressure distribution
	337	Float_t z[400];
	338	Float_t p[400];
	339
	340	for (Int_t i = 0; i < 400; i++)
	341	{
	342	z[i] = -20000. + Float_t(i) * 100;
	343	p[i] = GasPressureWeight(z[i]);
	344	}
	345
	346	TGraph* gr = new TGraph(400, z, p);
	347	TCanvas* c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
	348	c1->cd();
	349	gr->Draw("AL");
	350	}
fe4da5cc	351
d12ebe74	352	Int_t AliGenHalo::ReadNextParticle()
	353	{
	354	// Read the next particle from the file
	355	Int_t ncols = fscanf(fFile,"%d %f %f %d %f %d %f %f %f %f %f %f %f %f",
	356	&fLossID, &fLossT0, &fLossZ, &fLossA, &fLossW, &fPdg, &fXS, &fYS, &fZS, &fDX, &fDY, &fEkin, &fTS, &fWS);
	357	fLossZ /= 10.;
	358	fXS /= 10.;
	359	fYS /= 10.;
	360	fZS /= 10.;
	361	fTS *= 1.e-9;
	362	return (ncols);
	363	}