// 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                      Aug/30/2000 M. Iwasaki
//
// :: How to run ::
// 
// 1) to call the root, type   
//    
//     root
//
// 2)
//     root [0]  .x GoFullRecon.C
// 
// .. that's it.
//
// For example, if you like to run 300 events, type like
//     root [0]  .x GoFullRecon.C(300)
//

#include <stdlib.h>

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

  // Detectro 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 a histogram or two
  TH1F *ClsMult = new TH1F("ClsMult","Cluster Multiplicity",50, 0.0, 200.0);
  TH1F *ClsE    = new TH1F("ClsE","Cluster Energy",50, 0.0, 50.0);
  TH1F *TrkMult = new TH1F("TrkMult","Track Multiplicity",50, 0.0, 100.0);
  TH1F *TrkE    = new TH1F("TrkE","Track Energy",50, 0.0, 50.0);

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

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

    // Look at EM Clusters. Can access them by the [] operator.
    int nClusters = event->EMClusterLst()->GetEntries();
    ClsMult->Fill(nClusters);


    for (Int_t iClus = 0; iClus < nClusters; iClus++) {
      LCDCluster* cls = (LCDCluster*)event->EMClusterLst()->At(iClus);

      if (cls) {
	Double_t     eClus     = cls->GetEnergy();	
	// Histogram energies
	ClsE->Fill(eClus);
      }

    }
    

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

    // Track loop    
    for (Int_t ntrk=0; ntrk<nTracks; ntrk++) {
      LCDTrack* trk = (LCDTrack*)event->Tracks()->At(ntrk);
      
      if(trk) {
	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);

	// Histogram energies
	TrkE->Fill(E_trk);
      }

    }

    event->Clear();

  }

  ClsMult->SetFillColor(4);
  ClsE->SetFillColor(4);
  TrkMult->SetFillColor(4);
  TrkE->SetFillColor(4);

  TCanvas *c1 = new TCanvas("c1", "c1", 400, 600);
  c1->Divide(2,2);
  c1->cd(1); ClsMult->Draw("B");
  c1->cd(2); ClsE->Draw("B");
  c1->cd(3); TrkMult->Draw("B");
  c1->cd(4); TrkE->Draw("B");

  return 1;
}