//   Author:      Gero Flucke
//   Date:        June 1st, 2002
//   last update: $Date: 2006/01/11 12:23:12 $
//   by:          $Author: flucke $

#ifndef __GFDSTARANALYSEMC_H
#define __GFDSTARANALYSEMC_H

 // ROOT includes
#include "GFDstarAnalyse.h"


class TF1;
class TH2;
class H1PartDstar;
class GFHistArray;
class H1PartMC;
class H1ElectronTagger;

class GFDstarAnalyseMc : public GFDstarAnalyse {
public:
  GFDstarAnalyseMc();
  virtual ~GFDstarAnalyseMc();
  virtual void StartJob();
  virtual void StartRun(Int_t run);
  virtual void EndJob();
  virtual void PrepareEvent();

  Bool_t SetWeight();
  void SwitchOffExtraW();// {delete fWeightFunc; fWeightFunc = NULL;}

  virtual void OpenTree();
  
//   void FillHistsEffGen(const H1GtrRow* genDstar);
  void FillHistsGen(const H1PartMC* genDstar);
  void FillHistsHadronisation(const H1PartMC *dstar);
  void FillHistsHadronisationDsDiJet(const H1PartMC *dstar);
  Bool_t FillHistsEffTrigAccept(const H1PartDstar* trigDstar);
//   void FillHistsTrackLength(const H1PartDstar* dstar);
  void FillHistsYjbRes(const H1PartDstar* dstar);
  virtual void FillHistsTest(const H1PartDstar* dstar);
  void FillHistsWeight(TObjArray* dstars, TObjArray* genParts);
  void FillHistsEflow(const H1PartDstar*ds);
  void FillHistsJetQuark(const H1PartDstar* ds);
  void FillHistsDsJetQuark(const H1PartDstar *ds);
//   void FillHistsJetPurStab(const H1PartDstar* ds);
  void FillHistsDstarJetsGen(const H1PartMC *);
  void FillHistsDstarDiJetGen(const H1PartMC *genDstar);
  void FillHistsDsJetPurStab(const H1PartDstar* ds, Bool_t all);
  void FillHistsDsDiJetPurStab(const H1PartDstar* ds);
  void FillHistsDsJetTrig(const H1PartDstar* trigDstar);
  void FillHistsDsDiJetTrig(const H1PartDstar* trigDstar);

  virtual Bool_t IsS83(const H1PartDstar* dstar = NULL) const;
  virtual Bool_t IsS84(const H1PartDstar* dstar = NULL) const;
  virtual Bool_t IsS83NoL1Check(const H1PartDstar* dstar) const;
  virtual Bool_t IsS84NoL1Check(const H1PartDstar* dstar) const;
  virtual Bool_t IsS83NoL4Check(const H1PartDstar* dstar) const;
  virtual Bool_t IsS84NoL4Check(const H1PartDstar* dstar) const;
  virtual Bool_t IsL4Ref(Int_t trigger) const;
  virtual Bool_t IsL4Found(Int_t hqsel45bit, const H1PartDstar *ds, Int_t hqselMode = -2) const;
  virtual Bool_t IsS83Gen(const H1PartMC* dstar = NULL) const;
  virtual Bool_t IsS84Gen(const H1PartMC* dstar = NULL) const;

  Double_t GetYgen() const;
  Double_t GetYjbGen() const;
  Double_t GetYgenVis() const;
  Double_t GetQ2gen() const;
  Float_t GetWgammaPgen() const;
  Double_t GetZgen(const H1PartMC *genDs) const;
  Double_t GetXgammaMass(const H1PartMC* dstar,const H1PartJet* jet,Double_t yHfs) const;
  virtual Double_t GetXgammaMassHelp(const H1Part* dstar, const H1PartJet* jet) const;
  Double_t GetXgammaMassDiJet(const H1PartMC *dstar, const H1PartJet *dsJet, 
			      const H1PartJet *othJet, Double_t yHfs) const;
  virtual Double_t GetXgammaMassDiJetHelp(const H1Part *dstar, const H1PartJet *dsJet, 
					  const H1PartJet *othJet) const;
  Double_t GetXpMass(const H1PartMC* dstar,const H1PartJet* jet) const;
  virtual Double_t GetXpMassHelp(const H1Part* dstar, const H1PartJet* jet) const;
  Double_t GetGenEnergies(Double_t& forwBeam, Double_t& forward, Double_t& forwardBar, 
			  Double_t& barrel, Double_t& spacal) const;
  virtual Bool_t DstarCutLoose(const H1PartDstar* ds) const;
//   virtual Bool_t AreTracksOk(const H1PartDstar* dstar) const;

  Bool_t DstarCutLoose(const H1PartMC* ds) const;

  Bool_t DstarIsGenerated(const H1PartDstar* ds) const;
  Bool_t GetHeavyQuarks(H1PartMC*& quark, H1PartMC*& antiQuark) const;//, Int_t quarkPDG = 4) const;
  static Double_t Associate(const H1Part* ref1, const H1Part* ref2, 
			    H1PartMC*& outCloseTo1, H1PartMC*& outCloseTo2);
  TObjArray GetGenJets() const;
  Int_t     GetNonDstarGenJets(const H1PartMC* dstar, TObjArray& jets) const;
  Bool_t GetGenDiJets(const H1PartMC* dstar, const H1PartJet *&dsJet,
		      const H1PartJet *&oJet) const;
  TObjArray GetPartonJets() const;
  Bool_t GetPartDiJets(const H1PartMC* dstar, const H1PartJet *&dsJet,
		       const H1PartJet *&oJet, const TObjArray *allPartMc) const;

  const H1PartMC* GetCloseGenDsKpipi(const H1Part* ds, Option_t *opt) const;
  Bool_t IsAncestor(const H1PartMC *const ancest, const H1PartMC *part,
		    const TObjArray *allMc)const;
  void FillHistsPurStab(Double_t varRec, Double_t varGen, GFHistArray* hists, 
			Double_t weight = 1., const Bool_t* fill = NULL) const;

  void FillHistsBranching();

  Bool_t HasOffShellParton(const TObjArray *allPartMc) const;
  H1PartJet* GetDstarPartonJet(const TObjArray *partJets, const H1PartMC *ds,
			       const TObjArray *allPartMc) const;
  //  virtual void EndJob();  
protected:
  virtual Double_t GetWeight(Int_t trigger, const H1PartDstar *ds) const;
  Double_t GetMeanAccept44(Double_t y);
//   virtual void CreateHistsDm(); // overwritten/extended!
  virtual void CreateHistsDm(GFHistArray*& dmHist, GFHistArray*& dmHist83, 
			     GFHistArray*& dmHist84, const TString& nameAdd);
  virtual void CreateHistsD0();
  virtual void CreateHistsDstarJets();
  virtual void CreateHistsDstarDiJet();
  virtual void CreateHistsTracks(); 
  virtual void CreateHistsSumEt(); 
  void CreateHistsTrig();
  void CreateHistsAccept();
  void CreateHistsGen();
  void CreateHistsHadronisation();
  void CreateHistsHadronisationDsDiJet();

  void CreateHistsYjbRes(const char * namAdd = NULL);
  virtual void CreateHistsL4Eff2();
  void CreateHistsWeight();
  void CreateHistsEflow();
  void CreateHistsJetQuark();
  void CreateHistsDsJetQuark();

//   void CreateHistsJetPurStab();
  void CreateHistsDsJetPurStab();
  void CreateHistsDsDiJetPurStab();
  void CreateHistsDstarJetsGen();
  void CreateHistsDstarDiJetGen();
  void CreateHistsDsJetTrig();
  void CreateHistsDsDiJetTrig();

  void CreateHistsBranching();
  virtual void CreateHistsJetMult();


  void FillHistsGen(const H1PartMC* genDstar, const TObjArray* jets, 
		    GFHistArray* hists, Double_t weight) const;
//   void FillHistsRunDepAccept84(const H1PartDstar* trigDstar, Double_t y);
//   void CreateHistsRunDepAccept84(GFHistArray* notRunDepPtEtaPhiWgammaP, TDirectory* dir);
//   void CreateHistsRunDepTrig84(GFHistArray* notRunDepPtEtaPhiWgammaP, TDirectory* dir);

  void FillHistsJetQuark(const GFHistArray* jetQuarkHists, 
			 const H1PartJet* jet1, const H1PartJet* jet2, 
			 const H1PartMC* quark1, const H1PartMC* quark2);

//   void OpenDirCharmA();
//   void OpenDirCharmANoForwTr();
//   void OpenDirCharmANoTr();
//   void OpenDirCharmB();
//   void OpenDirCharmC();
//   void OpenDirCharmD();
//   void OpenResCharmA();
//   void OpenResCharmB();
//   void OpenCharmExcitationA();
//   void OpenCharmExcitationB();
//   void OpenCascadeA();
//   void OpenCascadeB();
//   void OpenCascadeC();
//   void OpenDirBeauty();
//   void OpenResBeauty();
//   void OpenTest();
  void HatSelectYear(const TString& year);
  virtual Float_t GetWgammaP33() const;
  virtual Float_t GetWgammaP44() const;
  virtual Float_t GetY33() const;
  virtual Float_t GetY44() const;

  Double_t GetEffCorrFactor(const H1PartDstar *ds) const;
  Double_t GetEffCorrFactor(Double_t pt, Double_t theta) const;

  Double_t     fEventCutQ2Gen;

  GFHistArray* fTrigHistsS83; // rec and trig (TE 31 && 19)
  GFHistArray* fTrigHistsS84;
  GFHistArray* fNonTrigHistsS83; // rec but maybe not trig
  GFHistArray* fNonTrigHistsS84;
  GFHistArray* fTrigHists19S83; //rec and TE 19
  GFHistArray* fTrigHists19S84;
  GFHistArray* fTrigHists31S83; // rec and TE 31
  GFHistArray* fTrigHists31S84;
//   GFHistArray* fTrigHistsS84Pt;     // histos for each run ranges of ETAG44
//   GFHistArray* fTrigHistsS84Eta;    // dito
//   GFHistArray* fTrigHistsS84Phi;    // dito
//   GFHistArray* fTrigHistsS84WgammaP;// dito

  GFHistArray* fGenHistsS83;
  GFHistArray* fGenHistsS84;
  GFHistArray* fDetAccHistsS83;

  GFHistArray* fAcceptHistsS83;
  GFHistArray* fAcceptHistsS84;
//   GFHistArray* fAcceptHistsS84Pt;       // histos for each run ranges of ETAG44
//   GFHistArray* fAcceptHistsS84Eta;      // dito
//   GFHistArray* fAcceptHistsS84Phi;      // dito
//   GFHistArray* fAcceptHistsS84WgammaP;  // dito

  GFHistArray* fGenAcceptHistsS83;
  GFHistArray* fGenAcceptHistsS84;
  GFHistArray* fGenNoAccWeightHistsS83;
  GFHistArray* fGenNoAccWeightHistsS84;
  GFHistArray* fGenWeightHistsS83;
  GFHistArray* fGenWeightHistsS84;


  GFHistArray* fWeightHistsS83;
  GFHistArray* fWeightHistsS84;

  GFHistArray* fYjbResHists;    // hists for y_jb resolution hists
  GFHistArray* fYjbResHistsS83; // hists for y_jb resolution hists
  GFHistArray* fYjbResHistsS84; // hists for y_jb resolution hists
  GFHistArray* fEflowHistsS83; // hists for energy flow in different detector regions
  GFHistArray* fEflowHistsS84; // hists for energy flow in different detector regions

  GFHistArray* fJetQuarkHists;
  GFHistArray* fJetQuarkHistsS83;
  GFHistArray* fJetQuarkHistsS84;
  Int_t        fNumJetQuarkHistsPerJet;
  GFHistArray* fDsJetQuarkHistsS83;

//   GFHistArray* fJetPurStabHists;  
//   GFHistArray* fJetPurStabHistsS83;
//   GFHistArray* fJetPurStabHistsS84;
//   GFHistArray* fJetPurStabPtInEtaHists;  
//   GFHistArray* fJetPurStabPtInEtaHistsS83;
//   GFHistArray* fJetPurStabPtInEtaHistsS84;

  // generator level D* + 1 jet hists:
  GFHistArray* fDstar1JetHistsGS83;        // just count
  GFHistArray* fDstar1JetHistsGAS83;       // acceptance weighted
  GFHistArray* fDstar1JetHistsGWS83;       // weighted
  GFHistArray* fDstar1JetHistsGNAWS83;     // weighted w/o acceptance part
  GFHistArray* fDstarBack1JetHistsGS83;    // now for eta(D*) < eta(jet)
  GFHistArray* fDstarBack1JetHistsGAS83;
  GFHistArray* fDstarBack1JetHistsGWS83;
  GFHistArray* fDstarBack1JetHistsGNAWS83;
  GFHistArray* fDstarForw1JetHistsGS83;    // now for eta(D*) > eta(jet)
  GFHistArray* fDstarForw1JetHistsGAS83;
  GFHistArray* fDstarForw1JetHistsGWS83;
  GFHistArray* fDstarForw1JetHistsGNAWS83;

  GFHistArray* fDstar1JetHistsGS84;        // now for S84
  GFHistArray* fDstar1JetHistsGAS84;
  GFHistArray* fDstar1JetHistsGWS84;
  GFHistArray* fDstar1JetHistsGNAWS84;
  GFHistArray* fDstarBack1JetHistsGS84;
  GFHistArray* fDstarBack1JetHistsGAS84;
  GFHistArray* fDstarBack1JetHistsGWS84;
  GFHistArray* fDstarBack1JetHistsGNAWS84;
  GFHistArray* fDstarForw1JetHistsGS84;
  GFHistArray* fDstarForw1JetHistsGAS84;
  GFHistArray* fDstarForw1JetHistsGWS84;
  GFHistArray* fDstarForw1JetHistsGNAWS84;

  GFHistArray* fDstarDiJetHistsGS83;
  GFHistArray* fDstarDiJetHistsGAS83;
  GFHistArray* fDstarDiJetHistsGWS83;
  GFHistArray* fDstarDiJetHistsGNAWS83;

  GFHistArray* fDstarDiJetHistsGS84;
  GFHistArray* fDstarDiJetHistsGAS84;
  GFHistArray* fDstarDiJetHistsGWS84;
  GFHistArray* fDstarDiJetHistsGNAWS84;

  // trig, non-trig and acceptance reco-level hists for D* + 1 jet
  GFHistArray* fNonTrigDs1JetHistsS83; // rec but maybe not trig
  GFHistArray* fNonTrigDs1JetHistsS84;
  GFHistArray* fTrigDs1JetHistsS83;    // rec and trig (TE 31 && 19)
  GFHistArray* fTrigDs1JetHistsS84;
  GFHistArray* fTrig19Ds1JetHistsS83;  //rec and TE 19
  GFHistArray* fTrig19Ds1JetHistsS84;
  GFHistArray* fTrig31Ds1JetHistsS83;  // rec and TE 31
  GFHistArray* fTrig31Ds1JetHistsS84;
  // now for eta(D*) < eta(jet)
  GFHistArray* fNonTrigDs1JetBackHistsS83;
  GFHistArray* fNonTrigDs1JetBackHistsS84;
  GFHistArray* fTrigDs1JetBackHistsS83;
  GFHistArray* fTrigDs1JetBackHistsS84;
  GFHistArray* fTrig19Ds1JetBackHistsS83;
  GFHistArray* fTrig19Ds1JetBackHistsS84;
  GFHistArray* fTrig31Ds1JetBackHistsS83;
  GFHistArray* fTrig31Ds1JetBackHistsS84;
  // now for eta(D*) > eta(jet)
  GFHistArray* fNonTrigDs1JetForwHistsS83;
  GFHistArray* fNonTrigDs1JetForwHistsS84;
  GFHistArray* fTrigDs1JetForwHistsS83;
  GFHistArray* fTrigDs1JetForwHistsS84;
  GFHistArray* fTrig19Ds1JetForwHistsS83;
  GFHistArray* fTrig19Ds1JetForwHistsS84;
  GFHistArray* fTrig31Ds1JetForwHistsS83;
  GFHistArray* fTrig31Ds1JetForwHistsS84;

  GFHistArray* fNonTrigDsDiJetHistsS83; // rec but maybe not trig
  GFHistArray* fNonTrigDsDiJetHistsS84;
  GFHistArray* fTrigDsDiJetHistsS83;    // rec and trig (TE 31 && 19)
  GFHistArray* fTrigDsDiJetHistsS84;
  GFHistArray* fTrig19DsDiJetHistsS83;  //rec and TE 19
  GFHistArray* fTrig19DsDiJetHistsS84;
  GFHistArray* fTrig31DsDiJetHistsS83;  // rec and TE 31
  GFHistArray* fTrig31DsDiJetHistsS84;

  GFHistArray* fDsJetPurStabHists;
  GFHistArray* fDsJetPurStabHistsDphi;
  GFHistArray* fDsJetPurStabHistsDeta;
  GFHistArray* fDsJetPurStabHistsDeltaR;
  GFHistArray* fDsJetPurStabHistsXgam;
  GFHistArray* fDsJetPurStabHistsDsPt;
  GFHistArray* fDsJetPurStabHistsJetPt;
  GFHistArray* fDsJetPurStabHistsJetEta;
  GFHistArray* fDsJetPurStabHistsJetPhi;
  GFHistArray* fDsJetPurStabHistsJetDaugh;
  GFHistArray* fDsJetPurStabHistsChFrac;
  GFHistArray* fDsJetPurStabHistsDsJetPt;
//   GFHistArray* fDsJetPurStabHistsDsJetPt2;
  GFHistArray* fDsJetPurStabHistsPtDsJetOrth;
  GFHistArray* fDsJetPurStabHistsDsJetM;
  GFHistArray* fDsJetPurStabHistsCosTh;

  GFHistArray* fDsJetPurStabHistsS83;
  GFHistArray* fDsJetPurStabHistsDphiS83;
  GFHistArray* fDsJetPurStabHistsDetaS83;
  GFHistArray* fDsJetPurStabHistsDeltaRS83;
  GFHistArray* fDsJetPurStabHistsXgamS83;
  GFHistArray* fDsJetPurStabHistsDsPtS83;
  GFHistArray* fDsJetPurStabHistsJetPtS83;
  GFHistArray* fDsJetPurStabHistsJetEtaS83;
  GFHistArray* fDsJetPurStabHistsJetPhiS83;
  GFHistArray* fDsJetPurStabHistsJetDaughS83;
  GFHistArray* fDsJetPurStabHistsChFracS83;
  GFHistArray* fDsJetPurStabHistsDsJetPtS83;
  GFHistArray* fDsJetPurStabHistsDsJetPt2S83;
  GFHistArray* fDsJetPurStabHistsPtDsJetOrthS83;
  GFHistArray* fDsJetPurStabHistsDsJetMS83;
  GFHistArray* fDsJetPurStabHistsCosThS83;

  GFHistArray *fDsPurStabHistsZDsS83;

  TObjArray*   fDsJetPurStabEtaInPtHistArrays;
  TObjArray*   fDsJetPurStabPtInEtaHistArrays;
  TObjArray*   fDsJetPurStabPtInPtHistArrays;
  TObjArray*   fDsJetPurStabEtaInPtHistArraysS83;
  TObjArray*   fDsJetPurStabPtInEtaHistArraysS83;
  TObjArray*   fDsJetPurStabPtInPtHistArraysS83;

  GFHistArray *fDsDiJetPurStabHistsS83; 
  GFHistArray *fDsDiJetPurStabHistsDetaS83;
  GFHistArray *fDsDiJetPurStabHistsDphiS83;
  GFHistArray *fDsDiJetPurStabHistsXgamS83;
  GFHistArray *fDsDiJetPurStabHistsCosThS83;
  GFHistArray *fDsDiJetPurStabHistsDijetPtS83;
  GFHistArray *fDsDiJetPurStabHistsDijetPt2S83;
  GFHistArray *fDsDiJetPurStabHistsDJetPtS83;
  GFHistArray *fDsDiJetPurStabHistsOthJetPtS83;
  GFHistArray *fDsDiJetPurStabHistsDJetEtaS83;
  GFHistArray *fDsDiJetPurStabHistsOthJetEtaS83;

  GFHistArray* fBranchingHists;

  GFHistArray *fHadronisationHistsS83;
  GFHistArray *fHadronisationGenHistsS83;
  GFHistArray *fHadPurStabHists;
  GFHistArray *fHadPurStabHistsJetPt;
  GFHistArray *fHadPurStabHistsJetEta;
  GFHistArray *fHadPurStabHistsDphi;
  GFHistArray *fHadPurStabHistsDeta;
  GFHistArray *fHadPurStabHistsXgam;
  GFHistArray *fHadPurStabHistsPtDsJet;
  GFHistArray *fHadPurStabHistsMDsJet;

  GFHistArray *fHadronisationDsDiJetHistsS83;
  GFHistArray *fHadronisationDsDiJetGenHistsS83;
  GFHistArray *fHadPurStabDijetHists;
  GFHistArray *fHadPurStabDijetHistsDphi;
  GFHistArray *fHadPurStabDijetHistsXgam;


  H1ElectronTagger** fTag84OfPeriod; // taggers for each 84 period

  Double_t fWeight83;           // storage for weights to apply
  Double_t fWeight84;           // dito
  Double_t fWeightPt;           // in case of pt-reweighting store that separately
  Double_t fInversePrescale33;  // to be taken into account in weight
  Double_t fInversePrescale44;  // to be taken into account in weight
  Double_t fY33Rec;             // store y calculated from smearing of yGen
  TF1*     fWeightFunc;         // weighting function
  TH1*     fMeanAccept84;
  TH2 *fHistEffCorrHighPt;
  TH2 *fHistEffCorrLowPt;

  ClassDef(GFDstarAnalyseMc,0)  // histogram book keeping, filling etc. for MC
};
#endif // __GFDSTARANALYSEMC_H