// Read the output root file from FullRecon.  
// The root file can be produced from ASCII file, output from FullRecon,
// through asc2root. 
// (Obtained from GISMO Output::
//    GISMO output --> asc2root --> FullRecon --> asc2root)  
// Sample program using JetFinder           May/12/2000 M. Iwasaki
//
// :: How to run ::
// 
// 1) to call the root, type   
//    
//     root
//
// 2)
//     root [0]  .x Jets_full.C
// 
// .. that's it.
//
// For example, if you like to run 300 events, type like
//     root [0]  .x Jets_full.C(300)
//

#include <stdlib.h>

int Jets_full(int nEvent=10)
{
  gROOT->Reset();

  // In order to refer to shareable libraries in this simple form,
  // include their paths in a .rootrc file.
  gSystem->Load("libLCDEvent.so");
  gSystem->Load("libLCDRootAppsUtil.so");
  gSystem->Load("libLCDFastMC.so");
  gSystem->Load("libLCDFullRecon.so");
  gSystem->Load("libLCDPhUtil.so");

  // Detector parameter file
  char* geomFile = 
    "/afs/slac/g/nld/RootApps/V2.1/ParFiles/Small.par";

  FILE* m_geomfile =  fopen(geomFile,"r");
  LCDGetParameters* gp = new LCDGetParameters(m_geomfile);

  TChain chain("T");
  // Please input your root (after Full Recon) file here 
  // For example, using chain, so as to add one or more root files
  chain.Add("test_recon.root");
//chain.Add("Second file"); //if you like.

  LCDEvent* event = new LCDEvent();
  chain->SetBranchAddress("LCDEvent", &event);

  // Define histograms
  TH1F *nJets   = new TH1F("nJets","# of Jets per Event (y=0.004)",20,0.,20);
  TH1F *nJets_1 = new TH1F("nJets_1","# of Jets per Event(y=0.006)",20,0.,20);
  TH1F *nJets_2 = new TH1F("nJets_2","# of Jets per Event(y=0.006, JADE)",20,0.,20);
  TH1F *JetP    = new TH1F("JetP","Jet momentum",100,0.,200);
  TH1F *JetPhi  = new TH1F("JetPhi","JetPhi Angle",100,0.,3.14);
  TH1F *JetCosTh= new TH1F("JetCosTh","Jet Cos Theta",50,-1.0,1.0);

  // Setup of JetFinder
  LCDJetFinder* jetfind = new LCDJetFinder(0.004);  // Default.. DURHAM algorithm

  TObjArray* veclist = new TObjArray();

  // Event loop
  for (int iEvent = 0; iEvent < nEvent; iEvent++) {

    // get the event      
    T->GetEvent(iEvent);

    // Look at tracks. Can access them by the [] operator.
    int nTracks = event->Tracks()->GetEntries();

    //  Event cut # track >= 2
    if (nTracks < 2) continue;

    veclist -> Delete(); // Clean up veclist

    //  Track loop    
    for (Int_t ntrk=0; ntrk<nTracks; ntrk++) {
	LCDTrack* trk = (LCDTrack*)event->Tracks()->At(ntrk);

	Double_t* tP3 = trk->GetMomentum();
	Double_t  E_trk = 0.1396 * 0.1396; // Assume pion mass
	for (int i=0; i<3; i++) { E_trk += (*(tP3+i))*(*(tP3+i)); }
	E_trk = sqrt(E_trk);

	TVector3 vec3(tP3);
	TLorentzVector* vec4 = new TLorentzVector(vec3,E_trk);

	veclist->Add(vec4);  // Make particle 4-vector list
    }
    

//  Do Jet Finding
    jetfind->SetPartList(veclist); // Input particle 4-vector list
    jetfind->SetDURHAM();          // Use DURHAM algorithm
    jetfind->SetYCut(0.004);       // Set YCut value
    jetfind->doFindJets();         // JetFindig
    
    Int_t numjets = jetfind->njets();
    nJets->Fill(numjets);

    for(Int_t j = 0; j < numjets; j++){

	TLorentzVector* Jet_vec = jetfind->jet4vec(j);
	TVector3        Pjet    = Jet_vec ->Vect();   // Jet 3 vector

	Double_t  Jet_mom = Pjet->Mag();
	Double_t  Jet_cos = (Pjet->Z())/Jet_mom;
	Double_t  Jet_phi = Pjet->Phi();

	JetP    ->Fill(Jet_mom);
	JetPhi  ->Fill(Jet_phi);
	JetCosTh->Fill(Jet_cos);
    }

//  Example 1  ... Do JetFindong with different Y-cut
    jetfind->SetYCut(0.0060);  // Set different YCut value
    jetfind->doFindJets();     // JetFinding

    Int_t numjets_1 = jetfind->njets();
    nJets_1->Fill(numjets_1);


//  Example 2  ... Do JetFindong with different Algorithm
    jetfind ->SetJADE();      // Use JADE algorithm
    jetfind->doFindJets();    // JetFinding

    Int_t numjets_2 = jetfind->njets();
    nJets_2->Fill(numjets_2);

    event->Clear();

  }
  veclist->Delete(); delete veclist;    

  nJets   ->SetFillColor(4);
  nJets_1 ->SetFillColor(4);
  nJets_2 ->SetFillColor(4);
  JetP    ->SetFillColor(4);
  JetPhi  ->SetFillColor(4);
  JetCosTh->SetFillColor(4);

  TCanvas *c1 = new TCanvas("c1", "c1", 700,900);
  c1->Divide(2,3);
  c1->cd(1); nJets->Draw("B");
  c1->cd(2); nJets_1->Draw("B");
  c1->cd(3); nJets_2->Draw("B");
  c1->cd(4); JetP->Draw("B");
  c1->cd(5); JetPhi->Draw("B");
  c1->cd(6); JetCosTh->Draw("B");

  return 1;
}