// ----------------------------------------------------------------------------
// $Id: root_GoThrust.C,v 1.7 2001/06/19 19:45:44 masako Exp $
// ----------------------------------------------------------------------------
// $Log: root_GoThrust.C,v $
// Revision 1.7  2001/06/19 19:45:44  masako
// Remove cross section.
//
// Revision 1.6  2001/06/19 16:52:06  toshi
// Use gSystem->Getenv() to get Parfiles directory.
//
// Revision 1.5  2001/05/16 17:46:56  masako
// Show the unit of cross section (fb)
//
// Revision 1.4  2001/05/11 21:52:41  toshi
// Remove .so suffix from gSystem->Load("...") in order to call macros from
// Windows system.
//
// Revision 1.3  2001/05/09 22:03:18  toshi
// Add $Id: root_GoThrust.C,v 1.7 2001/06/19 19:45:44 masako Exp $.
//
//
// Run Fast MC interactively with Root
// Sample program using EventShape       04/26/2001 Toshi
//
//
// :: How to run ::
// 
// 1) to call the root, type   
//    
//     root
//
// 2)
//     root [0]  .x root_GoThrust.C
// 
// .. that's it.
//
// For example, if you like to run 300 events, type like
//     root [0]  .x root_GoThrust.C(300)
//

#include <stdlib.h>

Int_t root_GoThrust(Int_t nEvent=100)
{
  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.
  gSystem->Load("libLCDEvent");
  gSystem->Load("libLCDRootAppsUtil");
  gSystem->Load("libLCDFastMC");
  gSystem->Load("libLCDPhUtil");

  LCDEvent event;

  // Input generator level root data
  Char_t* rootFile = "test.root";
  LCDreadRootFile source(rootFile,&event);

  // Open detector parameter file
  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);

  // If you want to change some detector parameter(s) for FastMC 
  // do here
  // For example..
  //gp.SetMagneticField(2.);  // Change Magnetif field
  //gp.SetCoilinHAD();        // Set HAD inside the Coil (This case
                              // you need to change some detector size also..)

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

  // Initialize FastMC
  LCDFastMC fmc(&gp,s_smearFile.Data(),&event); 
  // If you don't want the Cluster merging
  //fmc.SetNoClusterMerge();

  // 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;
  Int_t itrk;
  LCDTrack* trk;
  for (iEvent = 0; iEvent < nEvent; iEvent++) {
    if (!source.GetEvent()) break; // Read an event from the ROOT file.
    fmc.Doit();                    // Through FastMC

    if (iEvent == 0){
      // List Event Header information, for example
      LCDEventHeader* header = (LCDEventHeader*)event.EventHeader();
      printf("\n");
      printf("Process:: %s\n",header->GetProcessName());

      // List beam parameters 
      LCDBeam* beam = (LCDBeam*)event.BEAM();
      printf("Ecm = %f\n",beam->GetEcm());
      printf("electron polarization : %f  positron polarization : %f\n",
	     beam->GetPolElec(), beam->GetPolPosi());
      printf("==> effective polarization for electron : %f \n\n",
	     beam->GetEffectivePol());
    }

    // 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 (itrk=0; itrk < nTracks; itrk++) {
      trk = (LCDTrack*)event.Tracks()->At(itrk);
      new(veclist[itrk]) TVector3();
      TVector3* ptrk=(TVector3*)veclist[itrk];
      *ptrk=trk->GetMomentumVector(0.0);
      //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());
  }

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