// ----------------------------------------------------------------------------
// $Id: rootg4_GoThrust.C,v 1.2 2002/01/17 23:49:05 masako Exp $
// ----------------------------------------------------------------------------
//
// Read the root file from GEANT4, 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 rootg4_GoThrust.C
// 
// .. that's it.
//
// For example, if you like to run 20 events, type like
//     root [0]  .x rootg4_GoThrust.C(20)


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

  TString parfile_dir;
  parfile_dir += gSystem->Getenv("LCDROOT")   ;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");

  LCDEvent* event=new LCDEvent();

  // Input generated root file from GEANT4
  Char_t* rootFile  = "test_geant4_SD.root";
  LCDreadRootFile source(rootFile,event);

  // Must use the same detector name as Full simulation
  Char_t* geomFile = 
    "SDMar01.par"; // SD
  //"LDMar01.par"; // LD
  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";
  TString s_smearFile(parfile_dir);
  s_smearFile += smearFile;

  // Setup FullRecon
  LCDFullRecon *frec=
    new LCDFullRecon(&gp,s_smearFile.Data(), event);
  //new LCDFullRecon(&gp,s_smearFile.Data(), event, 1); // For GISMO

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