//
// 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;
class StMuEEmcCluster;

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

#include <iomanip>
#include <vector>

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

void runEEmcPreClusterPid ( Int_t nevents = 1000,  // Number of events
			 Char_t *muDst = "MuDst/mcPi0n_field_onpi0_10000_06TC05_20.MuDst.root"
			 );


// PID cuts

Int_t reject_electron ( StMuEEmcCluster *cluster );
Int_t reject_hadron   ( StMuEEmcCluster *cluster );

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

void runEEmcPreClusterPid ( Int_t nevents, 
			    Char_t *muDst
			    ){


  ////////////////////////////////////////////
  //
  // 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");
 

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


  Int_t nRejectElectron = 0;
  Int_t nRejectHadron   = 0;
  Int_t nRejectNet      = 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();

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

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

      if ( reject_electron( &cluster ) ) nRejectElectron++;
      if ( reject_hadron  ( &cluster ) ) nRejectHadron++;

      if ( reject_electron( &cluster ) ||
	   reject_hadron  ( &cluster ) ) nRejectNet++;
    
    }

    // Increment the event counter
    ievent++; 

  }

  std::cout << "In " << nevents << " events:" << std::endl;
  std::cout << "Found " << ncluster << " clusters" << std::endl;
  std::cout << "Rejected " << nRejectElectron << " e+/- candidates" << std::endl;
  std::cout << "Rejected " << nRejectHadron << " h+/- candiates" << std::endl;
  std::cout << "Rejected " << nRejectNet << " in total" << std::endl;

}

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

Int_t reject_electron ( StMuEEmcCluster *cluster ) {

  if ( cluster -> getEnergy(1) * 1E3 > 0.8 &&
       cluster -> getEnergy(1) * 1E3 < 1.4 &&
       cluster -> getEnergy(2) * 1E3 > 0.8 ) return 1;

  return 0;

}

Int_t reject_hadron ( StMuEEmcCluster *cluster ) { 

  // Reject if postshower shows energy
  if ( cluster -> getEnergy(3) * 1E3 > 0.075 ) return 1;

  Float_t tower = cluster -> getEnergy(0);
  assert ( tower > 0. ); // Need to handle this case

  Float_t pre1 = cluster -> getEnergy(1) * 1E3 / tower;
  Float_t pre2 = cluster -> getEnergy(2) * 1E3 / tower;

  if ( ( pre1 > 0.75 / tower && pre1 < 2.5 / tower ) &&
       ( pre2 > 0.75 / tower && pre2 < 2.5 / tower ) ) return 1;
  
  return 0;

}