#include "HadronSelection.h"
#include <iostream>
#include <math.h>

#include "EVENT/LCIO.h"
#include "EVENT/LCRunHeader.h"
#include "EVENT/LCParameters.h"
#include "EVENT/LCCollection.h"
#include "EVENT/SimCalorimeterHit.h"
#include "EVENT/SimTrackerHit.h"
#include "EVENT/CalorimeterHit.h"
#include "EVENT/TrackerHit.h"
#include "EVENT/Track.h"
#include "EVENT/MCParticle.h"
#include "IMPL/LCFlagImpl.h"
#include "EventHeader.hh"
#include "collection_names.hh"
#include "TriggerBits.hh"
#include "ErrorBits.hh"
#include "CellIndex.hh"
#include <sstream>
#include <assert.h>
#include <time.h>


using namespace std ;
using namespace lcio ;
using namespace marlin ;
using namespace CALICE ;


HadronSelection aHadronSelection ;


HadronSelection::HadronSelection() : Processor("HadronSelection") {
  
  // modify processor description
  _description = "HadronSelection makes histograms and writes to Root file ..." ;
  

  // register steering parameters: name, description, class-variable, default value

  registerProcessorParameter( "EcalCollectionName" , 
			      "Name of the Ecal hits collection"  ,
			      _EcalHitsColName ,
			      std::string("EmcCalorimeter_Hits") ) ;
  registerProcessorParameter( "HcalCollectionName" , 
			      "Name of the Hcal hits collection"  ,
			      _HcalHitsColName ,
			      std::string("AhcCalorimeter_Hits") ) ;
  registerProcessorParameter( "TcmtCollectionName" , 
			      "Name of the Tcmt hits collection"  ,
			      _TcmtHitsColName ,
			      std::string("TcmtCalorimeter_Hits") ) ;
  registerProcessorParameter( "Ebeam" , 
			      "Beam energy"  ,
			      Ebeam , float(1.) ) ;
  registerProcessorParameter( "EcalMIP" , 
			      "MIP value for ECAL"  ,
			      EcalMIP , float(1.) ) ;
  registerProcessorParameter( "EcalThresh" , 
			      "Threshold value for ECAL /MIP"  ,
			      EcalThresh , float(0.6) ) ;
  registerProcessorParameter( "HcalMIP" , 
			      "MIP value for HCAL"  ,
			      HcalMIP , float(1.) ) ;
  registerProcessorParameter( "HcalThresh" , 
			      "Threshold value for HCAL /MIP"  ,
			      HcalThresh , float(0.6) ) ;
  registerProcessorParameter( "TcmtMIP" , 
			      "MIP value for TCMT"  ,
			      TcmtMIP , float(1.) ) ;
  registerProcessorParameter( "TcmtThresh" , 
			      "Threshold value for TCMT /MIP"  ,
			      TcmtThresh , float(0.6) ) ;
  registerProcessorParameter( "Cerenkov" , 
			      "Cerenkov 1=on -1=off 0=not used"  ,
			      Cerenkov , int(0) ) ;
}


void HadronSelection::init() { 

  // usually a good idea to
  printParameters() ;

  _nRun = 0 ;
  _nEvt = 0 ;
  _nSel = 0 ;

   cout << " Hadron selection cuts" << endl;
   cout << " =======================" << endl;
   cout << " Ebeam = " << Ebeam << " GeV " << endl;
   EEcalMin = 300+100*(Ebeam-10)/70; 
   EHcalMin = 100+50*(Ebeam-10)/70; 
   ETcmtMin = 50+50*(Ebeam-10)/70; 
   ETotMax = 1.5*Ebeam; 
   E12Max = 50.;
   cout << " ERaw < " << ETotMax << " GeV" << endl;
   cout << " Muon veto: " << " EEcal > " << EEcalMin << " || EHcal > " << EHcalMin << " || ETcmt > " << ETcmtMin << " MIPs" << endl;
   cout << " Preshower veto: " << " E(layers 1, 2) > " << E12Max << " MIPs" << endl;
   if(Cerenkov==0) cout << " No Cerenkov requirement" << endl;
   if(Cerenkov==1) cout << " Cerenkov on required" << endl;
   if(Cerenkov==-1) cout << " Cerenkov off required" << endl;
   cout << " =======================" << endl;
   
}

void HadronSelection::processRunHeader( LCRunHeader* run) { 

  _nRun++ ;

  //  Ebeam=run->getParameters().getIntVal("beamEnergyMeV")*0.001;
 
  // cout << " Ebeam = " << Ebeam << endl;  

} 

void HadronSelection::processEvent( LCEvent * evt ) { 

  static bool firstEvent = true ;
    
  // this gets called for every event 
  // usually the working horse ...

  /*
  if(firstEvent) {
    char *file1 = new char[80];
    sprintf(file1,"%s%d%s","/r08/calice/drw/temp/Run",evt->getRunNumber(),".dat");
    pFile = fopen (file1,"w");
  }
  */

  loadHits(evt);

  Etotal=Eecal/250.+Ehcal/28.5;

  int state=evt->getParameters().getIntVal(PAR_RECO_STATE);
  TriggerBits trigger(evt->getParameters().getIntVal(PAR_TRIGGER_EVENT));

  int Class=0;
  Class=evt->getParameters().getIntVal("EventClass");
  if(Eecal<EEcalMin && Ehcal<EHcalMin && Etcmt<ETcmtMin) Class+=4;

  if((isMC || (state==kRecoStateBeam && trigger.isBeamTrigger()) )
     && Etotal<ETotMax// Total energy cut against double beams
     && (Eecal>EEcalMin || Ehcal>EHcalMin || Etcmt>ETcmtMin) // cut against muons
     && EecalLayer1<E12Max && EecalLayer2<E12Max // cut against preshowering 
     && (Cerenkov==0 || (Cerenkov>0 && trigger.isCherenkovTrigger()) || (Cerenkov<0 && !trigger.isCherenkovTrigger()) ) 
     ) 
    {
      _nSel++;
      Class+=2;
      setReturnValue(true);  
    }
 
  else {
    //      fprintf (pFile, "%d%s%d\n",evt->getRunNumber()," ",evt->getEventNumber());  
    setReturnValue(false);
  }
  
  evt->parameters().setValue("EventClass",Class);
  firstEvent = false ;
    _nEvt ++ ;
}



void HadronSelection::loadHits( LCEvent * evt ) { 

  // Zero some arrays

  Eecal=0;
  EecalLayer1=0;
  EecalLayer2=0;
  Ehcal=0;
  Etcmt=0;
  Nhits=0;
  isMC=false;

  typedef const std::vector<std::string> StringVec ;
  typedef std::vector<const std::string*> StringPVec ;
  StringVec* strVec = evt->getCollectionNames() ;

  // Get the calorimeter hits collections

  for( StringVec::const_iterator name = strVec->begin() ; name != strVec->end() ; name++){
    LCCollection* col = evt->getCollection( *name ) ;
    std::string sss;
    sss=name->c_str();
    int nHits =  col->getNumberOfElements() ;

    if(sss=="MCParticle") isMC = true;

    if(sss==_EcalHitsColName) {
      LCFlagImpl flag( col->getFlag() ) ;
      
      // Get total energy
      for( int i=0 ; i< nHits ; i++ ){        

	if(  col->getTypeName() == LCIO::SIMCALORIMETERHIT ) {
	  SimCalorimeterHit* hit = dynamic_cast<SimCalorimeterHit*>( col->getElementAt( i ) ) ;
	  Ehit=hit->getEnergy()/EcalMIP;
	  cellID = hit->getCellID0();
	}
	else{
	  CalorimeterHit* hit = dynamic_cast<CalorimeterHit*>( col->getElementAt( i ) ) ;
	  Ehit=hit->getEnergy()/EcalMIP;
	  cellID = hit->getCellID0();
	}
	CALICE::CellIndex CellID(cellID,0);
	K=CellID.getLayerIndex();

	if(Ehit>EcalThresh) {
	    Eecal+=Ehit;
	    if(K>10) Eecal+=Ehit;
	    if(K>20) Eecal+=Ehit;
	    if(K==1) EecalLayer1+=Ehit;
	    if(K==2) EecalLayer2+=Ehit;
	    Nhits++;
	}
      }
    }

    if(sss==_HcalHitsColName) {
      LCFlagImpl flag( col->getFlag() ) ;
      
      // Get total energy
      for( int i=0 ; i< nHits ; i++ ){        

	if(  col->getTypeName() == LCIO::SIMCALORIMETERHIT ) {
	  SimCalorimeterHit* hit = dynamic_cast<SimCalorimeterHit*>( col->getElementAt( i ) ) ;
	  Ehit=hit->getEnergy()/HcalMIP;
	  cellID = hit->getCellID0();
	}
	else{
	  CalorimeterHit* hit = dynamic_cast<CalorimeterHit*>( col->getElementAt( i ) ) ;
	  Ehit=hit->getEnergy()/HcalMIP;
	  cellID = hit->getCellID0();
	}
	CALICE::CellIndex CellID(cellID,0);
	K=CellID.getLayerIndex();

	if(Ehit>HcalThresh) {
	    Ehcal+=Ehit;
	}
      }
    }

    if(sss==_TcmtHitsColName) {
      LCFlagImpl flag( col->getFlag() ) ;
      
      // Get total energy
      for( int i=0 ; i< nHits ; i++ ){        

	if(  col->getTypeName() == LCIO::SIMCALORIMETERHIT ) {
	  SimCalorimeterHit* hit = dynamic_cast<SimCalorimeterHit*>( col->getElementAt( i ) ) ;
	  Ehit=hit->getEnergy()/TcmtMIP;
	  cellID = hit->getCellID0();
	}
	else{
	  CalorimeterHit* hit = dynamic_cast<CalorimeterHit*>( col->getElementAt( i ) ) ;
	  Ehit=hit->getEnergy()/TcmtMIP;
	  cellID = hit->getCellID0();
	}
	CALICE::CellIndex CellID(cellID,0);
	K=CellID.getLayerIndex();

	if(Ehit>TcmtThresh) {
	    Etcmt+=Ehit;
	}
      }
    }
  }
}


void HadronSelection::check( LCEvent * evt ) { 

  // nothing to check here - could be used to fill checkplots in reconstruction processor

}


void HadronSelection::end(){ 
  
  std::cout << "HadronSelection::end()  " << name() 
	    << " processed " << _nEvt << " events in " << _nRun << " runs.\n "
	    << "  *** Selected " << _nSel << " hadron candidate events. " 
	    << std::endl ;
  //  fclose (pFile);

}