// ----------------------------------------------------------------------------
// $Id: root_GoFullRecon.C,v 1.4 2001/06/23 16:05:29 toshi Exp $
// ----------------------------------------------------------------------------
// $Log: root_GoFullRecon.C,v $
// Revision 1.4  2001/06/23 16:05:29  toshi
// LCDFullRecon frec(...) -> LCDFullRecon *frec=new LCDFullRecon(...) .
//
// Revision 1.3  2001/06/23 05:47:28  toshi
// LCDEvent event; -> LCDEvent *event;
//
// Revision 1.2  2001/06/22 20:05:34  toshi
// Change input file name.
//
// Revision 1.1  2001/06/21 21:16:23  toshi
// This program is to display full simulation data written in ROOT file format.
//
//
// Read root output data, then do FullRecon
//
// Sample program                      06/21/2001 Toshi
//
// :: How to run ::
// 
// 1) to call the root, type   
//    
//     root
//
// 2)
//     root [0]  .x root_GoFullReconMC.C
// 
// .. that's it.
//
// For example, if you like to run 300 events, type like
//     root [0]  .x root_GoFullReconMC.C(300)
//

Int_t root_GoFullRecon(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  = "../";
  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("libLCDFullRecon");
  gSystem->Load("libLCDPhUtil");

  LCDEvent *event=new LCDEvent();

  // Input sio file (GISMO output)
  Char_t* rootFile  = "test_L2_pi0.root";
  LCDreadRootFile source(rootFile,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=new LCDFullRecon(&gp, s_smearFile.Data(), 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
  Int_t iEvent;
  for (iEvent = 0; iEvent < nEvent; iEvent++) {

    if (!source.GetEvent()) break; // Get an event from the input

    frec->Doit();                   // Do FullRecon

    LCDBeam* beam = event->BEAM();
    Double_t Ecm = beam->GetEcm();
    TVector3 IP  = beam->GetPos();
    if (iEvent == 0) printf(" Ecm = %f \n", Ecm);
    //printf(" beam position = %f %f %f \n",
    //	  IP.X(), IP.Y(), IP.Z());    

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

    for (Int_t iClus = 0; iClus < nClusters; iClus++) {
      LCDCluster* cls = (LCDCluster*)event->ClusterLst()->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", 600, 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 iEvent;
}