//-------------------------------------------------------------------------------------
// skeleton_aleph.cc
//
// This is a skeleton program to read/analyse Aleph files via root or standalnoe c++
//
// * You can run this macro using
//     $ root skeleton_aleph.cc
//
// * or if you can compile this file as follows (you may require root-shell):
//     $ g++ skeleton_aleph.cc -o skeleton_aleph `root-config --cflags --glibs`
//   and run:
//     $ ./skeleton_aleph
//
// Apr 2026
//-------------------------------------------------------------------------------------

// C++ headers
#include <iostream>
#include <fstream>
#include <cmath>
#include <vector>
#include <iomanip>
using namespace std;

// ROOT headers
#include "TFile.h"
#include "TTree.h"
#include "TChain.h"
#include "TH1F.h"
#include "TF1.h"
#include "TLegend.h"
#include "TLorentzVector.h"
#include "TVector3.h"
#include "TNtuple.h"

//makefile headers
#include "Track.C"
#include "Truth.C"
#include "Gamma.C"
#include "Vertex.C"
#include "Conversion.C"
#include "EventInfo.C"

// --- particle class derived from TLorentzVector ---
class Particle : public TLorentzVector{
 public:
  int ch, id, m1id, m1brc, m2id, m2brc;
  // default values
  Particle(){
    ch = 0;
    id = -1;
    m1brc = m1id = -1;
    m2brc = m2id = -1;
  }
  // Add two particle
  Particle operator+ (const Particle& rhs){
     Particle temp;
     temp.SetPxPyPzE((this->Px()+rhs.Px()),
                     (this->Py()+rhs.Py()),
                     (this->Pz()+rhs.Pz()),
                     (this->E() +rhs.E()) );
     return temp;
  }
  // copy all properties from TLorenzVector to Particle
  Particle& operator=(TLorentzVector rhs){
     this->SetPxPyPzE(rhs.Px(),
                      rhs.Py(),
                      rhs.Pz(),
                      rhs.E() );
     return *this;
  }
};

// --- analysis function ---
void skeleton_aleph(){

  // Output root file
  TFile *fileout = new TFile("alephoutput.root","recreate");

  // chains to read tree in data files
  TChain *ch1 = new TChain("Gamma","Gamma");
  TChain *ch2 = new TChain("Truth","Truth");
  TChain *ch3 = new TChain("EventInfo","EventInfo");

  // add all MC files for each chain
  // (MC: simulation, DA: recorded data)
  ch1->Add("/hepdata/Aleph/AlephMC9*.root");
  ch2->Add("/hepdata/Aleph/AlephMC9*.root");
  ch3->Add("/hepdata/Aleph/AlephMC9*.root");

  Gamma     *Gam = new Gamma(ch1);
  Truth     *Tru = new Truth(ch2);
  EventInfo *Evt = new EventInfo(ch3);

  // total events to process
  double fraction = 0.01;
  Long64_t nEvents = fraction * ch1->GetEntries();
  cout << "Total events desired = " << nEvents
       << " ( " << int(100*fraction) << "% ) "<< endl;

  //histograms for signal and background
  TH1F *hpt = new TH1F("hpt","hpt;pT [GeV];Entries",100,0,5);

  cout << "Looping over events ... " << endl;
  for(int ev = 0; ev<nEvents; ev++)
  {
    //get entries in the event
    Gam->GetEntry(ev);
    Tru->GetEntry(ev);
    Evt->GetEntry(ev);

    // Select good hadronic events (e- e+ -> Z -> q qbar)
    bool isGoodEvent = Evt->evt_goodevent->at(0);
    if(isGoodEvent==false) continue;

    // list of photons in the event
    vector<Particle> pht;

    // loop over photons
    for(unsigned int i=0; i<Gam->gam_px->size(); i++){
       // get the order in truth list
       unsigned int truthindex = Gam->gam_truthindex->at(i);
       if (truthindex>998) continue;

       //define temporary particle, set mass, and truth informations
       Particle g;
       g.SetXYZM(Gam->gam_px->at(i),
                 Gam->gam_py->at(i),
                 Gam->gam_pz->at(i),
                 0.0);

       // gamma momentum cut 1 GeV
       if(g.P() < 1.0) continue;

       hpt->Fill(g.Pt());

       g.m1id  = Tru->tru_m1id->at(truthindex);
       g.m1brc = Tru->tru_m1brc->at(truthindex);
       pht.push_back(g);
    }

    // combining two photons
    if(pht.size()>1){
      for(unsigned int i=0;   i<pht.size()-1; i++){
      for(unsigned int j=i+1; j<pht.size()  ; j++){
          Particle pi0 = pht[i] + pht[j];
          //cout << pi0.M() << endl;
      }
      }
    }

  } //end of event loop
  cout << "\nEvent loop completed." << endl;

  hpt->Write("hpt");
  fileout->Close();

} // end of anaysis

int main(int argc, char** argv){
  skeleton_aleph();
  return 0;
}