// ----------------------------------------------------------------------------
// $Id: sio_GoThrust.C,v 1.4 2001/06/19 16:52:14 toshi Exp $
// ----------------------------------------------------------------------------
// $Log: sio_GoThrust.C,v $
// Revision 1.4  2001/06/19 16:52:14  toshi
// Use gSystem->Getenv() to get Parfiles directory.
//
// Revision 1.3  2001/05/11 21:52:52  toshi
// Remove .so suffix from gSystem->Load("...") in order to call macros from
// Windows system.
//
// Revision 1.2  2001/05/09 22:07:19  toshi
// Add libSIO2root.so.
// Add $Id: sio_GoThrust.C,v 1.4 2001/06/19 16:52:14 toshi Exp $.
//
//
// Read the sio file from GISMO, and run FullRecon interactively.
// Sample program using EventShape         04/26/2001 Toshi
//
// :: How to run ::
// 
// 1) to call the root, type   
//    
//     root
//
// 2)
//     root [0]  .x sio_GoThrust.C
// 
// .. that's it.
//
// For example, if you like to run 20 events, type like
//     root [0]  .x sio_GoThrust.C(20)


Int_t sio_GoThrust(Int_t nEvent=1000) {
  
  gROOT->Reset();

  TString parfile_dir;
  parfile_dir += gSystem->Getenv("LCDROOT")   ;parfile_dir += "/";
  parfile_dir += gSystem->Getenv("LCDVERSION");parfile_dir += "/";
  parfile_dir += "ParFiles/";

  // In order to refer to shareable libraries in this simple form,
  // include their paths in a .rootrc file.

  // Using dev version Event Classes
  gSystem->Load("libLCDEvent");
  gSystem->Load("libLCDRootAppsUtil");
  gSystem->Load("libLCDFastMC");
  gSystem->Load("libLCDFullRecon");
  gSystem->Load("libLCDPhUtil");
  gSystem->Load("libSIO2root");
  gSystem->Load("libLCDSIOUtil");

  LCDEvent event;

  // Please input your GISMO output root file here 
  Char_t* sioFile = "test_L2_pi0.sio";
  LCDreadSIOFile source(sioFile,&event);

  // Must use the same detector name as Full simulation
  Char_t* geomFile = 
    "Silicon.par"; //Silicon
  //"Large.par";   // Large
  //"Small.par";   // Small
  TString s_geomFile(parfile_dir);
  s_geomFile += geomFile;
  FILE* m_geomfile =  fopen(s_geomFile.Data(),"r");
  LCDGetParameters gp(m_geomfile);
  fclose(m_geomfile);

  Char_t* smearFile = 
    "lookup_silicon.nobmcon";
  //"lookup_large.nobmcon";
  //"lookup_small.nobmcon";
  TString s_smearFile(parfile_dir);
  s_smearFile += smearFile;

  // Setup FullRecon
  LCDFullRecon frec(&gp, s_smearFile.Data(), &event);

  // 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;

  TClonesArray veclist("TVector3",100);

  // Event loop
  Int_t iEvent;
  for (iEvent = 0; iEvent < nEvent; iEvent++) {
    if (!source.GetEvent()) break;//read SIO data
    frec.Doit();//do full recon.

    // Look at tracks. Can access them by the [] operator.
    Int_t 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);
      TVector3  tP3 = trk->GetMomentumVector(0.0);
      TVector3* vec3= new(veclist[ntrk]) TVector3();
      *vec3 = tP3;
    }
    
    // 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());
    
  }

  veclist.Delete();

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

  TCanvas *c1 = new TCanvas("c1", "c1", 600,600);
  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 iEvent;
}