[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 "AliGenEventHeader.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;
  Int_t np = 0;
  
  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);
      np++;
      Int_t nc = ReadNextParticle();
      
      if (fLossID != oldID || nc == 0) {
	  oldID = fLossID;
	  break;
      }
  }

  SetHighWaterMark(nt);
  AliGenEventHeader* header = new AliGenEventHeader("HALO");
  header->SetNProduced(np);
  // Passes header either to the container or to gAlice
  if (fContainer) {
      header->SetName(fName);
      fContainer->AddHeader(header);
  } else {
      gAlice->SetGenEventHeader(header);	
  }
}


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