#if !defined(__CINT__) || defined(__MAKECINT__)

#include "fstream.h"

#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TFile.h"
#include "TRandom.h"
#include "TDatime.h"
#include "TSystem.h"
#include "TMath.h"

#endif

void generate_kstar(Int_t nOut = 1, Int_t nevent = 10000)
{
  //
  // Macro to generate phasespace distributions needed as GEANT input
  //

  //________________________________________________________________________
  //
  // Definitions
  //________________________________________________________________________
  //

  // PID
  const Int_t    kPid    = 33;
  // Mass
  const Double_t kMass   = 0.8961;
  // Slope of the mt-distribution (GeV)  
  const Double_t kSlope  = 0.28;
  // Width of the y-distribution (Gaussian)
  const Double_t kWidth  = 1.0;
  // Center of the y-distribution (Gaussian)
  const Double_t kCenter = 2.9;
  // Range of the y-distribution (Flat)
  const Double_t kYlow   = 0.0;
  const Double_t kYhigh  = 2.0 * kCenter;
  // Number of tracks per event
  const Int_t    kNtrack = 1;

  const Double_t kPi     = 3.141592654;

  TH1F *hY   = new TH1F("hY"  ,"Rapidity", 12,0.0,  6.0);
  TH1F *hPt  = new TH1F("hPt" ,"p_{t}"   , 12,0.0,  3.0);
  TH1F *hPhi = new TH1F("hPhi","#phi"    ,100,0.0,360.0);
  TH1F *hPx  = new TH1F("hPx" ,"p_{x}"   ,100,0.0,  5.0);
  TH1F *hPy  = new TH1F("hPy" ,"p_{y}"   ,100,0.0,  5.0);
  TH1F *hPz  = new TH1F("hPz" ,"p_{z}"   ,100,0.0,100.0);

  TH2F *hYPt = new TH2F("hYPt","p_{t} vs rapidity",12,0.0,6.0,12,0.0,3.0);

  // Gaussian or flat rapidity distribution
  Bool_t gaussian = kTRUE;

  // The output filenames
  Char_t cFileDat[150];
  Char_t cFileRoot[150];
  sprintf(cFileDat ,"gen_kstar_%02d.d"   ,nOut);
  sprintf(cFileRoot,"gen_kstar_%02d.root",nOut);

  //________________________________________________________________________
  //
  // Execution part
  //________________________________________________________________________
  //

  // Set a new random generator seed
  TDatime *datime = new TDatime();
  Int_t time = datime->GetTime();
  Int_t date = datime->GetDate();
  Int_t prid = gSystem->GetPid();
  delete datime;
  Int_t seed = TMath::Abs(10000 * prid + time - date);
  gRandom->SetSeed(seed);

  printf(" seed   = %d\n",seed);
  printf(" nevent = %d\n",nevent);

  // The mt-distribution
  TF1 *funMt = new TF1("funMt","x*exp(-sqrt([0]*[0]+x*x)/[1])",0.0,5.0);
  funMt->SetParameter(0,kMass);
  funMt->SetParameter(1,kSlope);

  // Open the output file
  FILE *fileData = fopen(cFileDat,"w+");  

  for (Int_t ievent = 0; ievent < nevent; ievent++) {

    fprintf(fileData," %4d       %6d\n",kNtrack,ievent);       

    for (Int_t itrack = 0; itrack < kNtrack; itrack++) {

      // Rapidity
      Double_t y;
      if (gaussian) {
        y = gRandom->Gaus(kCenter,kWidth);
      }
      else {
        y = gRandom->Rndm() * (kYhigh - kYlow) + kYlow;
      }
      // Pt
      Double_t pt  = funMt->GetRandom();   
      // Phi
      Double_t phi = gRandom->Rndm() * 2.0 * kPi;

      Double_t px  = pt * TMath::Cos(phi);   
      Double_t py  = pt * TMath::Sin(phi);   
      Double_t pz  = TMath::Sqrt(kMass*kMass + pt*pt) * TMath::SinH(y);   

      hY->Fill(y);
      hPt->Fill(pt);
      hYPt->Fill(y,pt);
      hPhi->Fill(phi*180.0/kPi);
      hPx->Fill(px);
      hPy->Fill(py);
      hPz->Fill(pz);

      fprintf(fileData,"                  %2d  %8.3f %8.3f %8.3f\n"
                      ,kPid,px,py,pz);   

    }

  }

  fclose(fileData);

  TFile *fileHist = new TFile(cFileRoot,"RECREATE");
  hY->Write();
  hPt->Write();
  hYPt->Write();
  hPhi->Write();
  hPx->Write();
  hPy->Write();
  hPz->Write();
  fileHist->Close();

  printf("generate done\n");

}