void runGenStarlight(Int_t nEvents=100, Int_t ifile = 0)
{
  gSystem->Load("libStarLight.so");
  gSystem->Load("libAliStarLight.so");

  TString HistName = "galice_";
  HistName += ifile;
  HistName += ".root";

  AliRunLoader* rl = AliRunLoader::Open(HistName.Data(),AliConfig::GetDefaultEventFolderName(),"RECREATE");
  rl->MakeTree("E");
  rl->MakeStack();
  AliStack* stack = rl->Stack();
  gAlice->SetRunLoader(rl);

  Int_t projA = 208; 
  Int_t projZ = 82;
  Int_t targA = 208; 
  Int_t targZ = 82;
  Double_t energyConfig = 5020.; // GeV/c2
  Float_t beam1energy = TMath::Sqrt(Double_t(projZ)/projA*targA/targZ)*energyConfig/2;
  Float_t beam2energy = TMath::Sqrt(Double_t(projA)/projZ*targZ/targA)*energyConfig/2;
  Float_t gamma1  = beam1energy/0.938272;
  Float_t gamma2  = beam2energy/0.938272;

  Int_t nwbins = 5000;
  Double_t wmin = 0.001022; // GeV/c2
  Double_t wmax = 0.03;
  Int_t nybins = 20;
  Double_t ymax = 7.05;

  Int_t seedConfig = ifile+1;

  AliGenStarLight* genStarLight = new AliGenStarLight(1000);
  genStarLight->SetParameter(Form("BEAM_1_Z     =    %3i    #Z of target",targZ));
  genStarLight->SetParameter(Form("BEAM_1_A     =    %3i    #A of target",targA));
  genStarLight->SetParameter(Form("BEAM_2_Z     =    %3i    #Z of projectile",projZ));
  genStarLight->SetParameter(Form("BEAM_2_A     =    %3i    #A of projectile",projA));
  genStarLight->SetParameter(Form("BEAM_1_GAMMA =    %6.1f    #Gamma of the target",gamma1));
  genStarLight->SetParameter(Form("BEAM_2_GAMMA =    %6.1f    #Gamma of the projectile",gamma2));
  genStarLight->SetParameter(Form("W_MAX        =    %.6f    #Max value of w",wmax));
  genStarLight->SetParameter(Form("W_MIN        =    %.6f    #Min value of w",wmin));
  genStarLight->SetParameter(Form("W_N_BINS     =     %3i    #Bins i w",nwbins));
  genStarLight->SetParameter(Form("RAP_MAX      =    %.2f    #max y",ymax));
  genStarLight->SetParameter(Form("RAP_N_BINS   =    %i      #Bins i y",nybins));
  genStarLight->SetParameter("CUT_PT       =    0    #Cut in pT? 0 = (no, 1 = yes)");
  genStarLight->SetParameter("PT_MIN       =    0    #Minimum pT in GeV");
  genStarLight->SetParameter("PT_MAX       =   10    #Maximum pT in GeV");
  genStarLight->SetParameter("CUT_ETA      =    0    #Cut in pseudorapidity? (0 = no, 1 = yes)");
  genStarLight->SetParameter("ETA_MIN      =   -5    #Minimum pseudorapidity");
  genStarLight->SetParameter("ETA_MAX      =    5    #Maximum pseudorapidity");
  genStarLight->SetParameter(Form("PROD_MODE    =    %i    #Production mode",1));
  genStarLight->SetParameter(Form("PROD_PID     =    %6i    #Channel of interest (not relevant for photonuclear processes)",11));
  genStarLight->SetParameter(Form("RND_SEED     =    %i    #Random number seed", seedConfig));
  genStarLight->SetParameter("BREAKUP_MODE  =   5    #Controls the nuclear breakup");
  genStarLight->SetParameter("INTERFERENCE  =   0    #Interference (0 = off, 1 = on)");
  genStarLight->SetParameter("IF_STRENGTH   =   1.   #% of interfernce (0.0 - 0.1)");
  genStarLight->SetParameter("INT_PT_MAX    =   0.24 #Maximum pt considered, when interference is turned on");
  genStarLight->SetParameter("INT_PT_N_BINS = 120    #Number of pt bins when interference is turned on");
  genStarLight->SetParameter("XSEC_METHOD   = 1      # Set to 0 to use old method for calculating gamma-gamma luminosity");

  genStarLight->SetStack(stack);
  genStarLight->Init();

  //TH1D* hM = new TH1D("hM",";m_{ee} (GeV/c^{2})",100,0.,1.); // pair mass
  TH2D* hPtY  = new TH2D("hPtY","ptY", 1000, 0. ,100., 200, -10, 10); // electron pt, y

  for (Int_t iev=0; iev<nEvents; iev++){
    if (iev%100==0) printf("Event: %i\n",iev);
    stack->Reset();
    genStarLight->Generate();
    TLorentzVector p;
    for (Int_t i=0;i<stack->GetNtrack();i++){
      TParticle* part = stack->Particle(i);
      part->Momentum(p);
      hPtY->Fill(p.Pt()*1000., p.PseudoRapidity());
    }
  }
  hPtY->Write();
}