#ifndef h1algo_TH1Algo_hh
#define h1algo_TH1Algo_hh
//_____________________________________________________________________________
// H1Algo algorithm to combine tracks and clusters for more precise
// determination of the hadronic final state
//
//Build and load
// libJetUser_forRoot.so (from JES group instructions)
// libStntuple_geom.so (gmake Stntuple._geom USESHLIBS=1)
// libStntuple_jet.so (gmake Stntuple._jet USESHLIBS=1)
//
//This class needs the following blocks 
//  THeaderBlock
//  TCalDataBlock
//  TStnJetBlock (whichever one you use)
//  TStnTrackBlock
//  TStnVertexBlock (ZVertexBlock)
//These need to be registered in your BeginJob method
// 
// To call the code you initialise the class only once eg in BeginJob() or
// constructor and, for example, make fH1Algo a data member
//
//   fH1Algo= new TH1Algo();
//
// Then in the event loop:
//
//  Int_t rc = fH1Algo->MakeH1Algo(TStnEvent* event, const char* jetName,
//		  Int_t level, Float_t zvtx=-999.0, Int_t nvtx=-99);
//
// event is pointer to the event, which you can get:
//  TStnEvent* event = fHeaderBlock->GetEvent();
//  jetName is the name of the block of jets to operate on, allowed:
//    "JetBlock"                     JetClu 0.4
//    "PROD@JetCluModule-cone0.7"    JetClu 0.7
//    "PROD@JetCluModule-cone1.0"    JetClu 1.0
//    "PROD@MidPointModule-cone0.4"  MidPoint 0.4
//    "PROD@MidPointModule-cone0.7"  MidPoint 0.7
//    "PROD@KtClusModule-cone0.7"    Kt 0.7
//  
//Level specifies the level of jet corrections to be applied
//  if(_level==0) no corrections
//  if(_level>=1) relativeEnergyCorrections
//  if(_level>=4) multipleInteractionEnergyCorrections
//  if(_level>=5) absoluteEnergyCorrections
//  if(_level>=6) underlyingEnergyCorrections
//  if(_level>=7) outOfConeEnergyCorrections
//
// zvtx, nvtx
//   the primary vertex and the number of class 12 ZVertices
//    these are optional, if you leave them off, then the code
//    will use the highest-sumpt zvertex as the vertex
//
//
// You can access the 4 vector of the final state by:
//
//  TLorentzVector* vec=fH1Algo->GetHad4Vector();
//
// Or loop over the jets (parallel to input jets):
//  for(Int_t ijet=0; ijet< fJetBlock->NJets(); ijet++) {
//    TLorentzVector* h1jet= fH1Algo->GetJet4Vector(ijet);
//   Float_t jetEta = jet->Eta();
//   Float_t jetPt  =jet->Pt();
//
//   Jet pt after upward systematic shift
//   Float_t jetPtUp  =jetPt*fH1Algo->GetJetSysUp(ijet);
//   Jet pt after downward systematic shift
//   Float_t jetPtDo  =jetPt*fH1Algo->GetJetSysDown(ijet);
//
//
//   Float_t trackfraction= fH1Algo->GetTrackFraction(ijet);
// }
//
// You can mask off tracks and/or towers, eg the electron or muon
// in W events, by:
//
//  LockTower(...);
//  LockTrack(...);      
//
// if you use these functions you must call
// 
//  ResetLockedTowers();
//  ResetLockedTracks(); 
//
//  at the begining of each event
//
// Any problems, bugs or suggestions please contact Andrew Mehta
//  (mehta@fnal.gov)
//
//
// Authors: Andrew Mehta
//  Latest:03/05/07 Bug fix locks now reset in constructor
//         19/10/06 All calibrations +systmatics implemented
//         26/09/06 Do not recalculate theta and phi
//         27/07/07 Systematic errors available
//         12/07/07 Apply L5-L7 calibration 
//         14/06/06 Stop using TClonesArray due to memory leak.
//                    Access to jets is now only via this class
//         07/05/07 Adapted to Stntuple style by Ray Culbertson
//_____________________________________________________________________________

#include <iostream>

#include <Stntuple/obj/TStnHeaderBlock.hh>
#include <Stntuple/obj/TStnTrackBlock.hh>
#include <Stntuple/obj/TStnVertexBlock.hh>
#include <Stntuple/obj/TCalDataBlock.hh>
//#include <Stntuple/obj/TStnElectronBlock.hh>
//#include <Stntuple/obj/TCesDataBlock.hh>
//#include <Stntuple/obj/TCprDataBlock.hh>
//#include <Stntuple/obj/TGenpBlock.hh>
//#include <Stntuple/obj/TStnClusterBlock.hh>
//#include <Stntuple/obj/TStnPhotonBlock.hh>
//#include <Stntuple/obj/TStnTriggerBlock.hh>
#include <Stntuple/obj/TStnVertexBlock.hh>
#include <Stntuple/obj/TStnJetBlock.hh>
#include "Stntuple/geom/TSimpleExtrapolator.hh"
#include "Stntuple/geom/TStnCalGeometry.hh"
#include "JetUser/JetEnergyCorrections.hh"
#include "TArrayS.h"

class TTrajectoryPoint;
class TSimpleExtrapolator;
class TStnCalGeometry;
class TStnEvent;

class TH1Algo: public TObject {
  
protected:

  Int_t fJetType;  // Jet type (0-2, cone, 3-4 midpoint, 5 kt)
  Int_t fCone;     // Cone size as used by JetEnergyCorrections 

  Float_t fConeSize[6];    //Cone size as used in h1 algo

  Double_t fSysMultiInter; //Mulitple interaction systematic
  Double_t fSysTrack;    //Track systematic

  Int_t fLevel;

  TSimpleExtrapolator* fExtrap;
  TStnCalGeometry* fGeom;

  //Standard  Ntuple blocks
  TStnHeaderBlock*    fHeaderBlock;
  TCalDataBlock*      fCalDataBlock;
  TStnJetBlock*       fJetBlock;
  TStnTrackBlock*     fTrackBlock;
  TStnVertexBlock*    fZVertexBlock;

  TLorentzVector fJet4Vector[500]; //H1 jet 4 vector
  Float_t fJetESysUp[500];   //Energy of jet after upwards systematic shift
  Float_t fJetESysDown[500]; //Energy of jet after downwards systematic shift

  Bool_t fSysUseMultiInter;  //Set false to turn off multiple interaction syst
  Bool_t fSysUseRel;         //Set false to turn off relative correction syst
  Bool_t fSysUseAbs;         //Set false to turn off absolute correction syst
  Bool_t fSysUseTrack;         //Set to false to turn off track systematic
  Bool_t fSysUseLevel6;         //Set to false to turn off level 6 systematic
  Bool_t fSysUseLevel7;         //Set to false to turn off level 7 systematic


  Float_t fH1trackfraction[500]; //Track fration of each jet
  
  Int_t NearJet(Double_t PhiTrack,Double_t EtaTrack);
  TArrayS fLockedTrack; //Array of locked tracks
  
  // 4 vector of final state after applying H1Algo algorithm
  TLorentzVector fHad4Vector; 
  // 4 vector of final state after applying H1Algo algorithm 
  // with upward energy shift
  TLorentzVector fHad4VectorESysUp; 
  // 4 vector of final state after applying H1Algo algorithm  
  // with downward energy shift
  TLorentzVector fHad4VectorESysDown; 
  Double_t fTrackFraction; // Track Fraction over all hadronic final state
  
  Int_t fNtow; //Internzl variable to keep track of #towers per track
  Float_t fTowerDistance[81];
  Int_t fNearTowersIndex[81];
  Int_t fNearTowersIeta[81];
  Int_t fNearTowersIphi[81];
  
  Int_t fNCombTowers; //Number of remaining calo towers after combination
  Int_t fNCombTracks; //Number of remaining tracks after combination
  
  void NearestTowers(Float_t etatrack,Float_t phitrack,
		     Int_t ietatrack,Int_t iphitrack);
  //Set if you want to include the mini plug in the algorithm
  Bool_t fUseMiniPlug; 
  //Unset if you do not want to include tracks in the algorithm
  Bool_t fUseTracks;  
  
  //Track quality cuts
  Float_t fTrMinPt; //minimum Track Pt
  Float_t fTrMaxPt; //maximum Track Pt
  Float_t fTrMinCotHits; //minimum COT hits in central region
  Float_t fTrEtaCotBoundary; // Boundary between COT and Si stand alone regions
  Float_t fTrZ0Max; // maximum Z0
  Float_t fTrZ0MinusZVetexMax; ; // maximum |Z0-zvertex|
  //Calo quality cuts
  Float_t fCaloMinEt; //minimum Et
  
  Float_t fExtraCaloFactor; //Extra factor applied to calo energy after H1 algo
  Float_t fExtraTrackFactor;//Extra factor applied to track energy after H1 algo
  Bool_t fOldCombinationMethod; //True if you wish to use the old combination
  // method. Warning this method has energy double counting
  Float_t fMaxTrackTowerEtaDist;//Maximum distance in eta for track tower match
  Float_t fMaxTrackTowerPhiDist;//Maximum distance in phi for track tower match
  
  Int_t fNvtx; //number of vertices
  Float_t fZv; //zvertex

  //Stores for quatities stored towerwise

  Short_t fTowerjetdo[52][48];   //error for calo energy up
  Short_t fTowerjetup[52][48];   //error for calo energy down
  Short_t fTowerjettrdo[52][48]; //error for track energy up
  Short_t fTowerjettrup[52][48]; //error for track energy down


  Float_t fETot[52][48];            //Energy in tower
  Float_t fETotdo[52][48];          //Energy in tower for systematics
  Float_t fETotup[52][48];          //    "  
  Float_t fETottrdo[52][48];        //    "  
  Float_t fETottrup[52][48];        //    "  
  Short_t fLockedTower[52][48]; //Stores which which towers are locked
  Short_t fTowerJet[52][48]; //Stores which towers have any energy
                       //This flag tells which
                       // jet the energy
                       //    in a tower belongs
                       // 0 no energy, 1 calo energy but no jet associated,
                       // >=2 jet index+2
                       //Same convention for fTowerTrack+asscociated errors

  Float_t fPtJet[52][48];        //pt of jet that tower belongs to
  
  Int_t fNtra[52][48]; //Number of tracks


  Short_t  fTowerTrack[52][48];  //=1 if there is a track in this tower
  Float_t fTrpt[52][48];         //total track 4 vector
  Float_t fTrpx[52][48];
  Float_t fTrpy[52][48];  
  Float_t fTrpz[52][48];
  Float_t fTre[52][48];

  Int_t fNjets;  //Number of jets

  JetEnergyCorrections* fJetcor; //Standard Jet Energy corrections

  Int_t fPrintLevel;

  void ResetArrarys();
  void JetLoop();
  void TowerLoop();
  void TrackLoop();
  void MakeH1Jets();
  void Level6Corrections(Int_t ijet);
  void Level7Corrections(Int_t ijet);
  Int_t UnpackBlocks(TStnEvent* event, const char* jetName);


 public:

  TH1Algo();
  ~TH1Algo();
  
  Float_t AbsoluteEnergyCorrections(const double pt, const double eta, 
				    int conesize);

  // run the algorithm
  Int_t MakeH1Algo(TStnEvent* event, const char* jetName,
		  Int_t level, Float_t zvtx=-999.0, Int_t nvtx=-99);
  

  TTrajectoryPoint* GetPointAtCalo(TStnTrack* track);
  // 4 vector of final state after applying H1Algo algorithm
  TLorentzVector* GetHad4Vector() {return &fHad4Vector;} 
  
  // 4 vector of jet after applying H1 algorithm
  TLorentzVector* GetJet4Vector(Int_t ijet) {return &fJet4Vector[ijet];} 

  Double_t GetTrackFraction() const {return fTrackFraction;} 
  Double_t GetTrackFraction(Int_t i) const {return fH1trackfraction[i];} 
  
  Double_t GetNCombTracks() const {
    //Number of remaining tracks after combination
    return fNCombTracks;
  } 
  Double_t GetNCombTowers() const {
    //Number of remaining calo towers after combination
    return fNCombTowers;
  } 

  //These functions allow user to exclude towers from total sums + jets
  // eg electron clusters or muon deposits
  // If any of these functions is used ResetLockedTowers() and/or  
  // ResetLockedTracks() must be called at begining of event loop
  void LockTowersInJet(Int_t ijet) ;

  inline void LockTower(Int_t ieta, Int_t iphi) {
    //Lock tower
    fLockedTower[ieta][iphi]=1;}
  
  inline void LockTower(TCalTower* to) {
    //Lock tower
    if(to) fLockedTower[to->IEta()][to->IPhi()]=1;
    else LockTowerError();
  }
  
  inline void LockTrack(Int_t TrackIndex) {
    //Tracks locked in similar way to towers
    if(TrackIndex>=fLockedTrack.GetSize()) fLockedTrack.Set(TrackIndex+1);
    fLockedTrack[TrackIndex]=1;
  }
  
  //Tracks locked in similar way to towers (slower than index method)
  void LockTrack(TStnTrack* track);     
      
  void LockTowerError() {
    std::cout <<"H1Algo, LockTower: Error in tower pointer. Tower not locked"<<std::endl;
  }

  void ResetLockedTowers() {memset(fLockedTower,0,52*48*sizeof(Short_t));}
  void ResetLockedTracks() {fLockedTrack.Reset();}
  Bool_t TrackSelection( TStnTrack* tr); //True if passes track selection
  Double_t CorrectedTheta(Double_t thetauncor,Double_t Zv);
  void MakeZVertexQuantities();
  //If Miniplug is to be used 
  void UseMiniPlug(Bool_t set=kTRUE) {fUseMiniPlug=set;} 
  //Set to false to ignore tracks
  void UseTracks(Bool_t set=kTRUE) {fUseTracks=set;} 

  void SetTrMinPt(Float_t TrMinPt) {fTrMinPt=TrMinPt;} 
  void SetTrMaxPt(Float_t TrMaxPt) {fTrMaxPt=TrMaxPt;} 
  void SetTrMinCotHits(Float_t TrMinCotHits)
    {fTrMinCotHits=TrMinCotHits;}
  void SetTrEtaCotBoundary(Float_t TrEtaCotBoundary)
    {fTrEtaCotBoundary=TrEtaCotBoundary;}
  void SetTrZ0Max(Float_t TrZ0Max){fTrZ0Max=TrZ0Max;}
  void SetCaloMinEt(Float_t CaloMinEt){fCaloMinEt=CaloMinEt;}
  void SetExtraCaloFactor(Float_t ExtraCaloFactor)
    {fExtraCaloFactor=ExtraCaloFactor;}
  void SetExtraTrackFactor(Float_t ExtraTrackFactor)
    {fExtraTrackFactor=ExtraTrackFactor;}
  void UseOldCombinationMethod (Bool_t in=kTRUE) 
    {fOldCombinationMethod=in;}

  void SetMaxTrackTowerEtaDist(Float_t MaxTrackTowerEtaDist)
  {fMaxTrackTowerEtaDist=MaxTrackTowerEtaDist;}
  void SetMaxTrackTowerPhiDist(Float_t MaxTrackTowerPhiDist)
  {fMaxTrackTowerPhiDist=MaxTrackTowerPhiDist;}

  //Energy of jet after upwards systematic shift
  Float_t GetJetEnSysUp(Int_t ijet) {return fJetESysUp[ijet];}   
  //Energy of jet after downwards systematic shift
  Float_t GetJetEnSysDown(Int_t ijet) {return fJetESysDown[ijet];}   

  //Fractional energy of jet after upwards systematic shift
  Float_t GetJetSysUp(Int_t ijet) {return ( fJet4Vector[ijet].E()>0)?  
      fJetESysUp[ijet]/fJet4Vector[ijet].E() : 0;}   

  //Fractional energy of jet after downwards systematic shift
  Float_t GetJetSysDown(Int_t ijet) {return ( fJet4Vector[ijet].E()>0)?  
      fJetESysDown[ijet]/fJet4Vector[ijet].E() : 0;}   


  //Set mulitple interaction systematic
  void SetSysMultiInter(Double_t sys) {fSysMultiInter=sys;} 
  //Set track systematic
  void SetSysTrack(Double_t sys) {fSysTrack=sys;}    

  //Set to false to turn off multiple interaction systematic
  void UseSysMultiInter(Bool_t in=kTRUE) {fSysUseMultiInter=in;}
  //Set to false to turn off relative correction systematic
  void UseSysRel(Bool_t in=kTRUE) {fSysUseRel=in;}
  //Set to false to turn off absolute correction systematic
  void UseSysAbs(Bool_t in=kTRUE) {fSysUseAbs=in;}
  //Set to false to turn off track systematic
  void UseSysTack(Bool_t in=kTRUE) {fSysUseTrack=in;}
  //Set to false to turn off level 6 systematic
  void UseSysLevel6(Bool_t in=kTRUE) {fSysUseLevel6=in;}
  //Set to false to turn off level 7 systematic
  void UseSysLevel7(Bool_t in=kTRUE) {fSysUseLevel7=in;}

  void SwitchOffAllSystematics();
  void SwitchOnAllSystematics();
  void SetPrintLevel(Int_t i) { fPrintLevel = i; }
  void Print();

  ClassDef(TH1Algo,1)

};
#endif