[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"
#include "AliLog.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");
    Int_t ir = 0;
    
    
    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);
	return;
    }

    if (fNskip > 0) {
      // Skip the first fNskip events
      SkipEvents();
    }
//
//
//
//    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 = 0;
    if (name) file = fopen(name, "r");
    if (!file) {
	AliError("No gas pressure file");
	return;
    }
    
    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++)
	{
	    ir = fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
	    if (ir == 0) break;
	    
	    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++)
	{
	    ir = fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
	    if (ir == 0) break;

	    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++) 
	{
	    ir = fscanf(file, "%f %e %e %e %e %e", &z, &p[0], &p[1], &p[2], &p[3], &p[4]);
	    if (ir == 0) break;
	    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);
    delete file;
}

//____________________________________________________________
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 && (fNskip == 0)) 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 - 1950./2.9979e10, 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);
}

void AliGenHalo::SkipEvents()
{
  //
  // Skip the first fNskip events
  Int_t skip = fNskip;
  Int_t oldID = -1;

  while (skip >= 0)
    {
      ReadNextParticle();
      if (oldID != fLossID) {
	oldID = fLossID;
	skip--;
      }
    } 
}
void AliGenHalo::CountEvents()
{
    // Count total number of events
    Int_t nev = 0;
    Int_t oldID = -1;
    Int_t nc = 1;
    while (nc != -1)
    {
	nc = ReadNextParticle();
	if (oldID != fLossID) {
	    oldID = fLossID;
	    nev++;
	    printf("Number of events %10d %10d \n", nev, oldID);
	}
    }


    rewind(fFile);
}
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"
418e091f	34	#include "AliLog.h"
f87cfe57	35
198bb1c7	36	ClassImp(AliGenHalo)
	37
	38	AliGenHalo::AliGenHalo()
d12ebe74	39	:AliGenerator(-1),
	40	fFile(0),
	41	fFileName(0),
	42	fSide(1),
	43	fRunPeriod(kY3D90),
	44	fTimePerEvent(1.e-4),
	45	fNskip(0),
	46	fZ1(0),
	47	fZ2(0),
	48	fG1(0),
	49	fG2(0),
	50	fGPASize(0),
	51	fLossID(0),
	52	fLossA(0),
	53	fPdg(0),
	54	fLossT0(0),
	55	fLossZ(0),
	56	fLossW(0),
	57	fXS(0),
	58	fYS(0),
	59	fZS(0),
	60	fDX(0),
	61	fDY(0),
	62	fEkin(0),
	63	fTS(0),
	64	fWS(0)
fe4da5cc	65	{
f87cfe57	66	// Constructor
d12ebe74	67
	68	fName = "Halo";
	69	fTitle = "Halo from LHC Beam";
fe4da5cc	70	//
fe4da5cc	71	// Read all particles
d12ebe74	72	fNpart = -1;
d12ebe74	73	SetAnalog(0);
fe4da5cc	74	}
	75
	76	AliGenHalo::AliGenHalo(Int_t npart)
1c56e311	77	:AliGenerator(npart),
d12ebe74	78	fFile(0),
	79	fFileName(0),
	80	fSide(1),
	81	fRunPeriod(kY3D90),
	82	fTimePerEvent(1.e-4),
	83	fNskip(0),
	84	fZ1(0),
	85	fZ2(0),
	86	fG1(0),
	87	fG2(0),
	88	fGPASize(0),
	89	fLossID(0),
	90	fLossA(0),
	91	fPdg(0),
	92	fLossT0(0),
	93	fLossZ(0),
	94	fLossW(0),
	95	fXS(0),
	96	fYS(0),
	97	fZS(0),
	98	fDX(0),
	99	fDY(0),
	100	fEkin(0),
	101	fTS(0),
	102	fWS(0)
fe4da5cc	103	{
f87cfe57	104	// Constructor
d12ebe74	105	fName = "Halo";
d12ebe74	106	fTitle= "Halo from LHC Beam";
fe4da5cc	107	//
d12ebe74	108	fNpart = npart;
	109	//
	110	SetAnalog(0);
fe4da5cc	111	}
	112
	113	//____________________________________________________________
	114	AliGenHalo::~AliGenHalo()
	115	{
f87cfe57	116	// Destructor
fe4da5cc	117	}
	118
	119	//____________________________________________________________
	120	void AliGenHalo::Init()
f87cfe57	121	{
f87cfe57	122	// Initialisation
d12ebe74	123	fFile = fopen(fFileName,"r");
1701cff6	124	Int_t ir = 0;
	125
	126
d12ebe74	127	if (fFile) {
	128	printf("\n File %s opened for reading, %p ! \n ", fFileName.Data(), (void*)fFile);
	129	} else {
	130	printf("\n Opening of file %s failed, %p ! \n ", fFileName.Data(), (void*)fFile);
418e091f	131	return;
d12ebe74	132	}
c321951c	133
	134	if (fNskip > 0) {
	135	// Skip the first fNskip events
	136	SkipEvents();
	137	}
d12ebe74	138	//
	139	//
	140	//
	141	// Read file with gas pressure values
	142	char *name = 0;
	143	if (fRunPeriod < 5) {
	144	name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/gasPressure.dat" );
	145	fGPASize = 21;
	146	fG1 = new Float_t[fGPASize];
	147	fG2 = new Float_t[fGPASize];
	148	fZ1 = new Float_t[fGPASize];
	149	fZ2 = new Float_t[fGPASize];
	150	} else if (fRunPeriod == 5) {
	151	name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_startup.dat");
	152	fGPASize = 18853;
	153	fG1 = new Float_t[fGPASize];
	154	fZ1 = new Float_t[fGPASize];
	155	} else if (fRunPeriod ==6) {
	156	name = gSystem->ExpandPathName("$(ALICE_ROOT)/LHC/pressure_2003_conditioned.dat");
	157	fGPASize = 12719;
	158	fG1 = new Float_t[fGPASize];
	159	fZ1 = new Float_t[fGPASize];
	160	} else {
	161	Fatal("Init()", "No gas pressure file for given run period !");
	162	}
	163
418e091f	164
	165	FILE* file = 0;
	166	if (name) file = fopen(name, "r");
	167	if (!file) {
	168	AliError("No gas pressure file");
	169	return;
	170	}
	171
d12ebe74	172	Float_t z;
	173	Int_t i;
	174	Float_t p[5];
	175
	176	const Float_t kCrossSection = 0.094e-28; // m^2
	177	const Float_t kFlux = 1.e11 / 25.e-9; // protons/s
	178	Float_t pFlux[5] = {0.2, 0.2, 0.3, 0.3, 1.0};
	179
	180	if (fRunPeriod < 5) {
	181	//
	182	// Ring 1
	183	//
	184
	185	for (i = 0; i < fGPASize; i++)
	186	{
1701cff6	187	ir = fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
	188	if (ir == 0) break;
	189
d12ebe74	190	fG1[i] = p[fRunPeriod];
	191
	192	if (i > 0) {
	193	fZ1[i] = fZ1[i-1] + z;
	194	} else {
	195	fZ1[i] = 20.;
	196	}
	197	}
	198	//
	199	// Ring 2
	200	//
	201	for (i = 0; i < fGPASize; i++)
	202	{
1701cff6	203	ir = fscanf(file, "%f %f %f %f %f %f", &z, &p[0], &p[1], &p[2] , &p[3], &p[4]);
	204	if (ir == 0) break;
	205
d12ebe74	206	fG2[i] = p[fRunPeriod];
	207	if (i > 0) {
	208	fZ2[i] = fZ2[i-1] + z;
	209	} else {
	210	fZ2[i] = 20.;
	211	}
	212	}
	213	//
	214	// Interaction rates
	215	//
	216	for (i = 0; i < fGPASize; i++)
	217	{
	218	fG1[i] = fG1[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
	219	fG2[i] = fG2[i] * kCrossSection * pFlux[fRunPeriod] * kFlux; // 1/m/s
	220	}
	221
	222	} else {
	223	for (i = 0; i < fGPASize; i++)
	224	{
1701cff6	225	ir = fscanf(file, "%f %e %e %e %e %e", &z, &p[0], &p[1], &p[2], &p[3], &p[4]);
1701cff6	226	if (ir == 0) break;
d12ebe74	227	z /= 1000.;
	228	fG1[i] = p[4] * kCrossSection * kFlux; // 1/m/s
	229	// 1/3 of nominal intensity at startup
	230	if (fRunPeriod == kLHCPR674Startup) fG1[i] /= 3.;
	231	fZ1[i] = z;
	232	}
	233	}
	234
	235
	236
	237
	238	//
	239	// Transform into interaction rates
	240	//
	241
	242
	243
	244
	245	Float_t sum1 = 0.;
	246	Float_t sum2 = 0.;
	247
	248	for (Int_t iz = 0; iz < 300; iz++) {
	249	Float_t zpos = 20. + iz * 1.;
	250	zpos *= 100;
	251	Float_t wgt1 = GasPressureWeight( zpos);
	252	Float_t wgt2 = GasPressureWeight(-zpos);
	253	sum1 += wgt1;
	254	sum2 += wgt2;
	255	}
	256	sum1/=250.;
	257	sum2/=250.;
	258	printf("\n %f %f \n \n", sum1, sum2);
de4ab12f	259	delete file;
f87cfe57	260	}
fe4da5cc	261
	262	//____________________________________________________________
	263	void AliGenHalo::Generate()
	264	{
d12ebe74	265	// Generate by reading particles from input file
fe4da5cc	266
d12ebe74	267	Float_t polar[3]= {0., 0., 0.};
fe4da5cc	268	Float_t origin[3];
fe4da5cc	269	Float_t p[3], p0;
d12ebe74	270	Float_t tz, txy;
d12ebe74	271	Float_t mass;
fe4da5cc	272	//
d12ebe74	273	Int_t nt;
	274	static Bool_t first = kTRUE;
	275	static Int_t oldID = -1;
	276	//
c321951c	277
	278	if (first && (fNskip == 0)) ReadNextParticle();
	279	first = kFALSE;
	280	oldID = fLossID;
fe4da5cc	281
fe4da5cc	282	while(1) {
d12ebe74	283	// Push particle to stack
d12ebe74	284	mass = TDatabasePDG::Instance()->GetParticle(fPdg)->Mass();
c321951c	285	p0 = TMath::Sqrt(fEkin * fEkin + 2. * mass * fEkin);
d12ebe74	286	txy = TMath::Sqrt(fDX * fDX + fDY * fDY);
73d273a4	287	if (txy > 1.) {
73d273a4	288	tz = 0.;
fe4da5cc	289	} else {
d12ebe74	290	tz = - TMath::Sqrt(1. - txy);
fe4da5cc	291	}
c321951c	292
d12ebe74	293	p[0] = p0 * fDX;
	294	p[1] = p0 * fDY;
	295	p[2] = p0 * tz;
fe4da5cc	296
d12ebe74	297	origin[0] = fXS;
	298	origin[1] = fYS;
	299	origin[2] = 1950.;
fe4da5cc	300
dc3747df	301	PushTrack(fTrackIt , -1, fPdg , p, origin, polar, fTS - 1950./2.9979e10, kPNoProcess, nt, fWS);
c321951c	302
d12ebe74	303	Int_t nc = ReadNextParticle();
	304
	305	if (fLossID != oldID \|\| nc == 0) {
	306	oldID = fLossID;
	307	break;
	308	}
	309	}
	310	SetHighWaterMark(nt);
	311	}
	312
fe4da5cc	313
d12ebe74	314	Float_t AliGenHalo::GasPressureWeight(Float_t zPrimary)
	315	{
	316	//
	317	// Return z-dependent gasspressure weight = interaction rate [1/m/s].
	318	//
	319	Float_t weight = 0.;
	320	zPrimary /= 100.; // m
	321	if (fRunPeriod < 5) {
	322	Float_t zAbs = TMath::Abs(zPrimary);
	323	if (zPrimary > 0.)
	324	{
	325	if (zAbs > fZ1[20]) {
	326	weight = 2.e4;
	327	} else {
	328	for (Int_t i = 1; i < 21; i++) {
	329	if (zAbs < fZ1[i]) {
	330	weight = fG1[i];
	331	break;
	332	}
	333	}
	334	}
	335	} else {
	336	if (zAbs > fZ2[20]) {
	337	weight = 2.e4;
	338	} else {
	339	for (Int_t i = 1; i < 21; i++) {
	340	if (zAbs < fZ2[i]) {
	341	weight = fG2[i];
	342	break;
	343	}
	344	}
	345	}
	346	}
	347	} else {
	348	Int_t index = TMath::BinarySearch(fGPASize, fZ1, zPrimary);
	349	weight = fG1[index];
	350	}
	351	return weight;
	352	}
fe4da5cc	353
d12ebe74	354	void AliGenHalo::Draw(Option_t *)
	355	{
	356	// Draws the gas pressure distribution
	357	Float_t z[400];
	358	Float_t p[400];
	359
	360	for (Int_t i = 0; i < 400; i++)
	361	{
	362	z[i] = -20000. + Float_t(i) * 100;
	363	p[i] = GasPressureWeight(z[i]);
	364	}
	365
	366	TGraph* gr = new TGraph(400, z, p);
	367	TCanvas* c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
	368	c1->cd();
	369	gr->Draw("AL");
	370	}
fe4da5cc	371
d12ebe74	372	Int_t AliGenHalo::ReadNextParticle()
	373	{
	374	// Read the next particle from the file
	375	Int_t ncols = fscanf(fFile,"%d %f %f %d %f %d %f %f %f %f %f %f %f %f",
	376	&fLossID, &fLossT0, &fLossZ, &fLossA, &fLossW, &fPdg, &fXS, &fYS, &fZS, &fDX, &fDY, &fEkin, &fTS, &fWS);
	377	fLossZ /= 10.;
	378	fXS /= 10.;
	379	fYS /= 10.;
	380	fZS /= 10.;
	381	fTS *= 1.e-9;
	382	return (ncols);
	383	}
c321951c	384
	385	void AliGenHalo::SkipEvents()
	386	{
	387	//
	388	// Skip the first fNskip events
	389	Int_t skip = fNskip;
	390	Int_t oldID = -1;
	391
	392	while (skip >= 0)
	393	{
	394	ReadNextParticle();
	395	if (oldID != fLossID) {
	396	oldID = fLossID;
	397	skip--;
	398	}
	399	}
	400	}
73d273a4	401	void AliGenHalo::CountEvents()
73d273a4	402	{
4a33c50d	403	// Count total number of events
73d273a4	404	Int_t nev = 0;
	405	Int_t oldID = -1;
	406	Int_t nc = 1;
	407	while (nc != -1)
	408	{
	409	nc = ReadNextParticle();
	410	if (oldID != fLossID) {
	411	oldID = fLossID;
	412	nev++;
	413	printf("Number of events %10d %10d \n", nev, oldID);
	414	}
	415	}
	416
	417
	418	rewind(fFile);
	419	}
	420
	421