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

#include <stdlib.h>

int Thrust_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");

  // Please input your GISMO output root file here 
  char* rootFile = "test.root";

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

  // Tracker smearing parameter file (used for Tracker reconstruction)
  char* smearFile = 
    "/afs/slac/g/nld/RootApps/V2.1/ParFiles/lookup_small.nobmcon";

  LCDreadRootFile *source = new LCDreadRootFile(rootFile);
  LCDFullRecon* frec = 0;
  LCDEvent*  event = 0;
  frec = new LCDFullRecon(source, geomFile, smearFile, 0);

  // Define histograms
  TH1F *THRUST = new TH1F("THRUST","Thrust", 50,0.,1.);
  TH1F *MAJOR  = new TH1F("MAJOR ","Major ", 50,0.,1.);
  TH1F *MINOR  = new TH1F("MINOR ","Minor ", 50,0.,1.);
  TH1F *THRUSTDIRZ = new TH1F("THRUSTDIRZ","Thrust Direction Z", 50,-1.,1.);  

  // Set up Thrust Finder
  LCDEventShape* eshape = new LCDEventShape();

  TObjArray* veclist = new TObjArray();

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

    // Do FullRecon
    frec->DoEventNewstyle();
    event = source->GetEventPointer();
    if (!event) break;

    // 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();
	TVector3* vec3 = new TVector3(tP3);

	veclist->Add(vec3);  // Make particle 3-vector list
    }
    
    // Do Thrust Finding
    eshape->SetPartList(veclist); // Input particle 3-vector list
    
    Double_t thrust = eshape->GetThrust();
    Double_t major  = eshape->GetMajor();
    Double_t minor  = eshape->GetMinor();
    
    THRUST->Fill(thrust);
    MAJOR ->Fill(major);
    MINOR ->Fill(minor);
    THRUSTDIRZ->Fill(eshape->GetThrustAxis()->Z());
    
    
    event->Clear();
  }

  veclist->Delete(); delete veclist;    

  THRUST    ->SetFillColor(4);
  MAJOR     ->SetFillColor(4);
  MINOR     ->SetFillColor(4);
  THRUSTDIRZ->SetFillColor(4);

  TCanvas *c1 = new TCanvas("c1", "c1", 700,900);
  c1->Divide(2,2);
  c1->cd(1); THRUST->Draw("B");
  c1->cd(2); MAJOR->Draw("B");
  c1->cd(3); MINOR->Draw("B");
  c1->cd(4); THRUSTDIRZ->Draw("B");

  return 1;
}