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