//
// Runs StMuEEmcPreCluster maker, performs a (trivial) clustering
//   around the high-tower in each event, and generates some histograms
//   of tower, preshower and postshower response, and SMD response...
//
// Binning on the histograms may not be optimal...
//


// forward declarations
class StChain;
class StIOMaker;
class StMuDstMaker;
class StMuEEmcPreAnalysisMaker;
class StMuEEmcPreClusterMaker;
class StMuDbReader;
class StEEmcDbMaker;
class St_db_Maker;

StChain                  *chain;
StIOMaker                *ioMaker;
StMuDstMaker             *muDstMaker;
StMuEEmcPreAnalysisMaker *muEEmcAnal;
StMuEEmcPreClusterMaker  *muEEmcClust;
StMuDbReader             *muDb;
StEEmcDbMaker            *eemdDb;
St_db_Maker              *starDb;

#include <iomanip>
#include <vector>

void runEEmcPreCluster ( Int_t nevents = 1000,  // Number of events
			 Char_t *muDst = "MuDst/mcPi0n_field_ongamma_10000_06TC05_10.MuDst.root",
			 Char_t *output = "output.root" // Output file
			 );

/////////////////////////////////////////////////////////////////////////////

void runEEmcPreCluster ( Int_t nevents, 
			 Char_t *muDst, 
			 Char_t *output ) {


  if ( output != "" ) 
    TFile *f = new TFile( output, "RECREATE" );

  ////////////////////////////////////////////
  //
  // Preshower, postshower and tower 2D histos
  //
  TH2F *hPre1VsTow  = new TH2F("hPre1VsTow", "Preshower(1) energy deposit vs reco tower E", 125, 0., 25., 500, 0., 80.);
  TH2F *hPre2VsTow  = new TH2F("hPre2VsTow", "Preshower(2) energy deposit vs reco tower E", 125, 0., 25., 500, 0., 80.);
  TH2F *hPostVsTow  = new TH2F("hPostVsTow", "Postshower   energy deposit vs reco tower E", 125, 0., 25., 500, 0., 80.);
  TH2F *hPre1VsPre2 = new TH2F("hPre1VsPre2","Preshower(1) energy deposit vs preshower(2) energy deposit", 500, 0., 80., 500, 0., 80.);
  //--
  TH2F *hPre1VsTowR  = new TH2F("hPre1VsTowR", "Preshower(1) energy deposit / tower E vs tower E", 125, 0., 25., 500, 0., 10.0);
  TH2F *hPre2VsTowR  = new TH2F("hPre2VsTowR", "Preshower(2) energy deposit / tower E vs tower E", 125, 0., 25., 500, 0., 10.0);
  TH2F *hPostVsTowR  = new TH2F("hPostVsTowR", "Postshower   energy deposit / tower E vs tower E", 125, 0., 25., 500, 0., 10.0);
  //--

  TH1F *hClusterEnergy = new TH1F("hClusterEnergy","Cluster Energy [GeV]",125,0.,25.);
  TH1F *hSeedEnergy    = new TH1F("hSeedEnergy",   "Seed energy [GeV]",125,0.,25.);
  TH2F *hClustVsSeed   = new TH2F("hClustVsSeed",  "Cluster Energy [GeV] vs Seed Energy [GeV]",125,0.,25.,125,0.,25.);
  


  ////////////////////////////////////////////
  //
  // Load shared libraries
  //
  gROOT -> LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
  loadSharedLibraries();
  gSystem -> Load("StEEmcUtil.so");
  gSystem -> Load("StMuEEmcPreAnalysisMaker.so");
  gSystem->Load("StDbLib");
  gSystem->Load("StDbBroker");
  gSystem->Load("St_db_Maker");
  gSystem -> Load("StEEmcDbMaker.so");
 
  // Load in and compile pre/postshower pid cuts
  gROOT -> LoadMacro("particleId.C+");


  ////////////////////////////////////////////
  //
  // Create the chain and add the muDst maker
  //
  chain      = new StChain("StChain");
  muDstMaker = new StMuDstMaker(0,0,"",muDst,"MuDst.root",1);
  muDb       = StMuDbReader::instance();

  eemcDb = new StEEmcDbMaker("eemcDb");
  starDb = new St_db_Maker("StarDb", "MySQL:StarDb");

  assert(eemcDb);
  //  eemcDb -> setTimeStampDay(20030516);
  eemcDb -> setTimeStampDay(20040101);
  eemcDb -> setPreferedFlavor( "set492", "eemcPMTcal" );

  muEEmcAnaly = new StMuEEmcPreAnalysisMaker();
  muEEmcAnaly -> setScaleFactor(0.8); // approx to get MC correct

  muEEmcClust = new StMuEEmcPreClusterMaker();
  muEEmcClust -> setSeedEnergy(0.3); // 0.3 GeV seeds (default 0.7)

  ////////////////////////////////////////////
  //
  // Initialize the chain in preparation to loop over events
  //
  chain -> Init();
  chain -> ls(3);

  ////////////////////////////////////////////
  // 
  // Event Loop
  //
  
  Int_t status = 0;
  Int_t ievent = 0;
  Int_t npass  = 0;
  Int_t ncluster = 0;
  while ( !status && ievent < nevents ) {

    //    if ( !ievent % 100 ) 
    std::cout << "Processing event number " << ievent << std::endl;

    // Clear all makers on the chain
    chain -> Clear();

    // Call all of the makers on the chain
    status = chain -> Make();


    // Get high-tower (aka seed for this application) energy
    Float_t energy_seed = muEEmcAnaly -> getHighTower() -> getEnergy(0);
    hSeedEnergy -> Fill(energy_seed);


    for ( Int_t i = 0; i < muEEmcClust -> getNumClusters(); i++ ) {

      // Get the cluster
      StMuEEmcCluster cluster = muEEmcClust -> getCluster(i);
      ncluster++;

      // Abort if it fails event selection.
      //if ( !select_pi0( &cluster ) ) continue;
      npass++;

      Float_t energy_tow  = cluster.getEnergy(0) ;
      Float_t energy_pre1 = cluster.getEnergy(1) * 1E3;
      Float_t energy_pre2 = cluster.getEnergy(2) * 1E3;
      Float_t energy_post = cluster.getEnergy(3) * 1E3;

      hClusterEnergy -> Fill(energy_tow);
      hClustVsSeed -> Fill(energy_seed,energy_tow);

      hPre1VsTow -> Fill( energy_tow, energy_pre1 );
      hPre2VsTow -> Fill( energy_tow, energy_pre2 );
      hPostVsTow -> Fill( energy_tow, energy_post );
      hPre1VsPre2 -> Fill ( energy_pre2, energy_pre1 );

      if ( energy_tow == 0 ) continue;
      hPre1VsTowR -> Fill( energy_tow, energy_pre1/energy_tow );
      hPre2VsTowR -> Fill( energy_tow, energy_pre2/energy_tow );
      hPostVsTowR -> Fill( energy_tow, energy_post/energy_tow );



      std::cout << " Filling " 
		<< energy_tow << " "
		<< energy_pre1 / energy_tow << " "
		<< energy_pre2 / energy_tow << " "
		<< energy_post / energy_tow << " "
		<< std::endl;
    
    }

    // Increment the event counter
    ievent++; 

  }

  std::cout << std::endl << std::endl;
  std::cout << "Found " << npass << " pi0 candidates out of " << ncluster << " clusters in " << ievent << " events" << std::endl;



  f->cd();
  muEEmcClust -> GetHistList() -> Write();

  f -> Write();

}