{
// ================================================================ //
// This is an example macro on usage of the fast simulation for the
// luminosity spectrometer.
//
// Author: Miro Helbich
//
// To run it: first download the .h and .cxx files for the classes
// GenGamma and LumiFast.
//
// then type:
//
// > root
// > .x lumi_example.C
//
// Good luck!
// ================================================================ //

// set the root ploting parameters
gROOT->SetStyle("Plain");
gStyle->SetLabelSize(0.05,"X");
gStyle->SetLabelSize(0.05,"Y");
gStyle->SetTitleXSize(0.05);
gStyle->SetStatH(0.35);
gStyle->SetStatW(0.3);
gStyle->SetOptStat();

// load the classes the files GenGamma.h,Gengamma.html,LumiFast.h,LumiFast.html
// should be in your home directory
gROOT->LoadMacro("GenGamma.cxx");
gROOT->LoadMacro("LumiFast.cxx");

// instantiate the generator and set the minimum energy cut and the tilt
GenGamma *gen = new GenGamma();
gen->SetEnergyMin(5.);
gen->SetYTilt(.5);

// instantiate the fast generator - switch on the physics and set Bdl
lf = new LumiFast();
lf->SetGenerator(gen);
lf->SetScattering(1);
lf->SetEnergyLoss(1);
lf->SetBdl(0.275);

//lf->SetNyStripsDn(13); // switch this to have antisymmetric spectrometer

Float_t Lumi = 1e7*0.0867; //luminosity * barn->m2*conversion_probability
Float_t xsec= lf->GetGenerator()->GetXsec(); // total cross-section

// book some histograms

TH1F *hgamma = new TH1F("hgamma","Generated spectrum",100,5.,30.);
TH1F *hcoinc = (TH1F*) hgamma->Clone("hcoinc");
hcoinc->SetTitle("Measured energy spectrum");
hcoinc->SetXTitle("E_{#gamma} [GeV]");
haccept = (TH1F*) hgamma->Clone("haccept");
haccept->SetTitle("The acceptance");
TH1F *hbdl = new TH1F("hbdl","Measured Bdl",100,lf->GetBdl()*0.8,lf->GetBdl()*1.2); 
hbdl->SetXTitle("Bdl [Tm]");
TH1F *hyprof = new TH1F("hyprof","y-profile",100,-5.,5.);
TH1F *hyup = new TH1F("hyup","y in LDS UP",11,4.7,4.7+11*0.8);


// start the run
Float_t bdl_norm = lf->GetBdl();

Int_t Nevents=100000; // number of events to be generated
Int_t NSingleUp=0,NSingleDn=0,NCoinc=0; // counters

gBenchmark->Start("anal"); // time the simmulation

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

  lf->NewEvent();

  // put all statements under condition like this to see if the photon
  // went through the slit
  //
  //if (lf->GetYPosGamTrue()<0.1&&lf->GetYPosGamTrue()>-0.1) {
  
  hgamma->Fill(lf->GetEnergyTrueUp()+lf->GetEnergyTrueDn());

  // chcek for the fixed collimator
  if (lf->GetYPosGamTrue()<3.5&&lf->GetYPosGamTrue()>-3.5) {
    hyprof->Fill(lf->GetYPosGam());
    hyup->Fill(lf->GetYPosUp());
  }

  // check if there was a coincidence
  if (lf->IsHitUp()&&lf->IsHitDn()) {
    hcoinc->Fill(lf->GetEnergyUp()+lf->GetEnergyDn());
    hbdl->Fill(lf->GetEnergyUp()*(lf->GetYPosUp()-lf->GetYPosGam())/0.3/(lf->GetZSpectro()-lf->GetZMagnet()));
  }

  // this includes the energy window
  if (lf->IsCoincidence()) haccept->Fill(gen->GetEnergy());

  // increment the counters
  if (lf->IsHitUp()) NSingleUp++;
  if (lf->IsHitDn()) NSingleDn++;
  if (lf->IsCoincidence()) NCoinc++;

}

gBenchmark->Show("anal");
{

  // show the results
cout << "singles rates dn= " << (Float_t) NSingleDn/Nevents * xsec *Lumi;
cout <<" +- " << (Float_t) sqrt(NSingleDn)/Nevents * xsec *Lumi <<endl;
cout << "singles rates up= " << (Float_t) NSingleUp/Nevents * xsec *Lumi;
cout <<" +- " << (Float_t) sqrt(NSingleUp)/Nevents * xsec *Lumi <<endl;
cout << "coincidence rate= " << (Float_t) NCoinc/Nevents * xsec *Lumi;
cout <<" +- " << (Float_t) sqrt(NCoinc)/Nevents * xsec *Lumi <<endl;

}

// plot some histograms
TCanvas *c1 = new TCanvas("c1");
c1->Divide(2,2); // same as zon 2 2

// show the generated energy spectrum
c1->cd(1); 
hgamma->Draw();

// show the accepetance: measured spectrum has to be divided by the
// generated one
c1->cd(2);
hcoinc->Divide(hgamma);

// show the measured y-profile
c1->cd(3);
hyprof->Draw();

// show the measure Bdl
c1->cd(4);
hbdl->Draw();

}