// ----------------------------------------------------------------------------
// $Id: sio_SeeCalHit.C,v 1.1 2001/06/21 22:20:41 masako Exp $
// ----------------------------------------------------------------------------
// $Log: sio_SeeCalHit.C,v $
// Revision 1.1  2001/06/21 22:20:41  masako
// Example file to handle CalHit.
//
//
//  Example file to handle CalHit.
//  
//  In this file, 1) read sio event (output from GISMO)
//                2) access CalHit
//                3) get E Theta Phi and R
//
//                                   Jun. 2001   M. Iwasaki
//
// :: How to run ::
// 
// 1) to call the root, type   
//    
//     root
//
// 2)
//     root [0]  .x sio_SeeCalHit.C
// 
// .. that's it.
//
// For example, if you like to run 3 events, type like
//     root [0]  .x sio_SeeCalHit(3)
//

#include <stdlib.h>

int sio_SeeCalHit(int nEvent=5) {
  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("libLCDFullRecon");
  //  gSystem->Load("libLCDPhUtil");
  gSystem->Load("libSIO2root");
  gSystem->Load("libLCDSIOUtil");

  LCDEvent event;

  // Input generated stdhep file
  Char_t* sioFile  = "test_L2_pi0.sio";
  //Char_t* sioFile  = "test_top_SD.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);

  // Setup CalHitUtil
  LCDCalHitUtil chutil(&gp);

  // Define a histogram or two
  TH1F *EM_E    = new TH1F("EM_E"    ,"EM energy",50, 0.0, 2.5);
  TH1F *EM_phi  = new TH1F("EM_phi"  ,"EM phi"   ,50, 0.0, 6.3);
  TH1F *EM_theta= new TH1F("EM_theta","EM theta" ,50, -1., 1.0);
  TH1F *EM_r    = new TH1F("EM_r"    ,"EM r"     ,50, 0.0, 300.0);
  TH1F *HAD_E    = new TH1F("HAD_E"    ,"HAD energy",50, 0.0, 2.5);
  TH1F *HAD_phi  = new TH1F("HAD_phi"  ,"HAD phi"   ,50, 0.0, 6.3);
  TH1F *HAD_theta= new TH1F("HAD_theta","HAD theta" ,50, -1., 1.0);
  TH1F *HAD_r    = new TH1F("HAD_r"    ,"HAD r"     ,50, 0.0, 300.0);

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

    if (!source.GetEvent()) break;

    LCDCalHit*  cE;
    Int_t       sysid;
    Double_t     Ecal, phi, theta, Rcyl;

    Int_t nCal = event->CAL()->GetEntries();
    if (nCal > 0) {
      for (int i_em=0; i_em<nCal; i_em++) {
	cE = (LCDCalHit*)event->CAL()->At(i_em);

	Ecal  = chutil.GetEtot(cE);
	phi   = chutil.GetPhi(cE);
	theta = TMath::Cos(chutil.GetTheta(cE));	
	Rcyl  = chutil.GetRcyl(cE);

	sysid=cE->GetTower()->GetSystem();

	switch(sysid) {
	case 10: // EM Barrel
	  EM_E    ->Fill(Ecal);
	  EM_phi  ->Fill(phi);
	  EM_theta->Fill(theta);
	  EM_r    ->Fill(Rcyl);
	  break;
	case 11: // EM EndCap
	  EM_E    ->Fill(Ecal);
	  EM_phi  ->Fill(phi);
	  EM_theta->Fill(theta);
	  EM_r    ->Fill(Rcyl);
	  break;
	case 12: // HAD Barrel
	  HAD_E    ->Fill(Ecal);
	  HAD_phi  ->Fill(phi);
	  HAD_theta->Fill(theta);
	  HAD_r    ->Fill(Rcyl);
	  break;
	case 13: // HAD EndCap
	  HAD_E    ->Fill(Ecal);
	  HAD_phi  ->Fill(phi);
	  HAD_theta->Fill(theta);
	  HAD_r    ->Fill(Rcyl);
	  break;
	default:
	  break;
	}
      }
    }    
  }

  EM_E->SetFillColor(4);
  EM_phi->SetFillColor(4);
  EM_theta->SetFillColor(4);
  EM_r->SetFillColor(4);
  HAD_E->SetFillColor(4);
  HAD_phi->SetFillColor(4);
  HAD_theta->SetFillColor(4);
  HAD_r->SetFillColor(4);

  TCanvas *c1 = new TCanvas("c1", "c1", 600, 600);
  c1->Divide(2,2);
  c1->cd(1); EM_E->Draw("B");
  c1->cd(2); EM_phi->Draw("B");
  c1->cd(3); EM_theta->Draw("B");
  c1->cd(4); EM_r->Draw("B");

  TCanvas *c2 = new TCanvas("c2", "c2", 600, 600);
  c2->Divide(2,2);
  c2->cd(1); HAD_E->Draw("B");
  c2->cd(2); HAD_phi->Draw("B");
  c2->cd(3); HAD_theta->Draw("B");
  c2->cd(4); HAD_r->Draw("B");

  return iEvent;
}