/***************************************************************************
 *
 * $Id: StHiSpectra.cxx,v 1.2 2002/04/03 00:23:27 jklay Exp $                                    
 *
 * Author: Bum Choi, UT Austin, Apr 2002
 *
 ***************************************************************************
 *
 * Description:  Class for making highpt inclusive spectra from highpt
 *               uDST's
 *
 ***************************************************************************
 * 
 * $Log: StHiSpectra.cxx,v $
 * Revision 1.2  2002/04/03 00:23:27  jklay
 * Fixed private member access bugs in analysis code
 *
 * Revision 1.1  2002/04/02 20:05:18  jklay
 * Bums analysis tools for highpt uDSTs
 *
 *
 **************************************************************************/
#include "StHiSpectra.h"
#include <iostream>
using namespace std;

const char* border = "*************************************************";

//__________________

StHiSpectra::StHiSpectra(const char* inputDir,const char* match,
			 const char* outRootName,
			 const char* f)
  : StHiBaseAnalysis(inputDir,match,outRootName)
     
{

  mEfficiencyFileName = f;
  mEfficiencyMap = 0;
}

//__________________

StHiSpectra::~StHiSpectra()
{
}

//__________________
//
// Get the efficiency file...
//

Int_t
StHiSpectra::initMore()
{

  //  char name[100];
  Int_t stat=0;

  TString fName = "fEfficiency";
  TFile* file=0; 
  if(mEfficiencyFileName==""){
    cout << "no efficiency file set " << endl;
  }
  else{
    file=new TFile(mEfficiencyFileName.Data());
    if(!file ||!file->IsOpen()){
      cout << "Cannot find efficiency file : " << mEfficiencyFileName << endl;
    }
    else{
        mEfficiencyMap = (TF2*)file->Get(fName.Data());
	delete file;
    }
  }
  if(!mEfficiencyMap){
    cout << "Cannot find efficency function : " << fName.Data() <<endl;
    cout << "WARNING-USING DUMMY EFFIENCY" << endl;
    mEfficiencyMap = new TF2("fEfficiency","1");
  }
  
  
  cout << "\tefficiency file : " << mEfficiencyFileName << endl;
  cout << "\tefficiency fcn title : " << mEfficiencyMap->GetTitle() << endl;

  return stat;
}


void
StHiSpectra::initHistograms()
{
  cout << "StHiSpectra::initHistograms()" << endl;

  //*********************

  using namespace Bin;

  //*********************

  //  gStyle->SetPalette(1,0);

  char title[500],name[500];

  h1CentralityCut = new TH1D("h1CentralityCut","centrality cut",
			     nFlowCentBin,flowCentMin,flowCentMax);

  h1Centrality = new TH1D("h1Centrality","centrality",
			     nFlowCentBin,flowCentMin,flowCentMax);

  h1VertexZ = new TH1D("h1VertexZ","vertex z",
		       nVertexZEvtBin,vertexZEvtMin,vertexZEvtMax);

  cout << "using bins" << endl;
  for(int iBin=0; iBin<mNVarBin; iBin++){
    varBin[iBin].ptBinsAry = new TArrayD;
    initPtAry(varBin[iBin].ptBinsAry,iBin);
    int nBin = varBin[iBin].ptBinsAry->GetSize()-1;
    varBin[iBin].nPtBinAry = nBin;;
    cout << "Bin " << iBin << " : " << endl;
    for(int i=0; i<nBin; i++){
      cout << varBin[iBin].ptBinsAry->At(i) << ", ";
    }
    cout << varBin[iBin].ptBinsAry->At(nBin) << endl;
  }
  //------------------------------------------------------------

  TString sPM[3] = { "Plus","Minus","Charged"};

  //********************************************************
  // number of events used to scale the spectra

  h1NEvent = new TH1D("h1NEvent","h1NEvent",1,1,2);

  // eta cut

  h1EtaCut = new TH1D("h1EtaCut","h1EtaCut",2,0,2);

  //********************************************************
  // yield for both signs

  h2YieldPr
    = new TH2D("h2YieldPr","h2YieldPr",
	       nEtaSmallBin,etaSmallMin,etaSmallMax,
	       nPtBin,ptMin,ptMax);
  h2YieldPr->Sumw2();
  h2YieldGl
    = new TH2D("h2YieldGl","h2YieldGl",
	       nEtaSmallBin,etaSmallMin,etaSmallMax,
	       nPtBin,ptMin,ptMax);
  h2YieldGl->Sumw2();


  //********************************************************
  // var bin
  
  for(Int_t iBin=0; iBin<mNVarBin; iBin++){

    //**** spectra plus, minus, charged

    for(Int_t iCharge=0; iCharge<3; iCharge++){

      // mean stuff
      //

      setName(name,"h1OneOverPtGl",iBin,sPM[iCharge].Data());
      varBin[iBin].meanPM[iCharge].h1OneOverPtGl =
	new TH1D(name,name,varBin[iBin].nPtBinAry,
		 varBin[iBin].ptBinsAry->GetArray());
      varBin[iBin].meanPM[iCharge].h1OneOverPtGl->Sumw2();
      setName(name,"h1OneOverPtPr",iBin,sPM[iCharge].Data());
      varBin[iBin].meanPM[iCharge].h1OneOverPtPr =
	new TH1D(name,name,varBin[iBin].nPtBinAry,
		 varBin[iBin].ptBinsAry->GetArray());
      varBin[iBin].meanPM[iCharge].h1OneOverPtPr->Sumw2();

      setName(name,"h1WeightedMeanGl",iBin,sPM[iCharge].Data());
      varBin[iBin].meanPM[iCharge].h1WeightedMeanGl =
	new TH1D(name,name,varBin[iBin].nPtBinAry,
		 varBin[iBin].ptBinsAry->GetArray());
      varBin[iBin].meanPM[iCharge].h1WeightedMeanGl->Sumw2();
      setName(name,"h1WeightedMeanPr",iBin,sPM[iCharge].Data());
      varBin[iBin].meanPM[iCharge].h1WeightedMeanPr =
	new TH1D(name,name,varBin[iBin].nPtBinAry,
		 varBin[iBin].ptBinsAry->GetArray());
      varBin[iBin].meanPM[iCharge].h1WeightedMeanPr->Sumw2();

      //
      // raw
      //
      setName(name,"h1RawGl",iBin,sPM[iCharge].Data());
      varBin[iBin].specPM[iCharge].h1RawGl =
	new TH1D(name,name,varBin[iBin].nPtBinAry,
		 varBin[iBin].ptBinsAry->GetArray());
      varBin[iBin].specPM[iCharge].h1RawGl->Sumw2();
      setName(name,"h1RawPr",iBin,sPM[iCharge].Data());
      varBin[iBin].specPM[iCharge].h1RawPr =
	new TH1D(name,name,varBin[iBin].nPtBinAry,
		 varBin[iBin].ptBinsAry->GetArray());
      varBin[iBin].specPM[iCharge].h1RawPr->Sumw2();

      //
      // efficiency corrected
      //
      setName(name,"h1EffCorrectedGl",iBin,sPM[iCharge].Data());
      varBin[iBin].specPM[iCharge].h1EffCorrectedGl =
	new TH1D(name,name,varBin[iBin].nPtBinAry,
		 varBin[iBin].ptBinsAry->GetArray());
      varBin[iBin].specPM[iCharge].h1EffCorrectedGl->Sumw2();
      setName(name,"h1EffCorrectedPr",iBin,sPM[iCharge].Data());
      varBin[iBin].specPM[iCharge].h1EffCorrectedPr =
	new TH1D(name,name,varBin[iBin].nPtBinAry,
		 varBin[iBin].ptBinsAry->GetArray());
      varBin[iBin].specPM[iCharge].h1EffCorrectedPr->Sumw2();
      
      //
      // corrected (background + momentum resolution)
      //
      
      setName(name,"h1CorrectedGl",iBin,sPM[iCharge].Data());
      varBin[iBin].specPM[iCharge].h1CorrectedGl =
	new TH1D(name,name,varBin[iBin].nPtBinAry,
		 varBin[iBin].ptBinsAry->GetArray());
      varBin[iBin].specPM[iCharge].h1CorrectedGl->Sumw2();
      setName(name,"h1CorrectedPr",iBin,sPM[iCharge].Data());
      varBin[iBin].specPM[iCharge].h1CorrectedPr =
	new TH1D(name,name,varBin[iBin].nPtBinAry,
		 varBin[iBin].ptBinsAry->GetArray());
      varBin[iBin].specPM[iCharge].h1CorrectedPr->Sumw2();
      
      
    }

    //***** ratios

    setName(name,"h1RawRatioGl",iBin);
    varBin[iBin].h1RawRatioGl
      = new TH1D(name,name,varBin[iBin].nPtBinAry,
		 varBin[iBin].ptBinsAry->GetArray());
    varBin[iBin].h1RawRatioGl->Sumw2();
    setName(name,"h1RawRatioPr",iBin);
    varBin[iBin].h1RawRatioPr
      = new TH1D(name,name,varBin[iBin].nPtBinAry,
		 varBin[iBin].ptBinsAry->GetArray());
    varBin[iBin].h1RawRatioPr->Sumw2();

    setName(name,"h1CorrectedRatioGl",iBin);
    varBin[iBin].h1CorrectedRatioGl
      = new TH1D(name,name,varBin[iBin].nPtBinAry,
		 varBin[iBin].ptBinsAry->GetArray());
    varBin[iBin].h1CorrectedRatioGl->Sumw2();
    setName(name,"h1CorrectedRatioPr",iBin);
    varBin[iBin].h1CorrectedRatioPr
      = new TH1D(name,name,varBin[iBin].nPtBinAry,
		 varBin[iBin].ptBinsAry->GetArray());
    varBin[iBin].h1CorrectedRatioPr->Sumw2();

  }  // varBin

}

//_____________________

void 
StHiSpectra::trackLoop()
{
  if(mDebug)
    cout << "StHiSpectra::spectraTrackLoop()" << endl;

//  cout << "<WARNING>::SPECTRA TRACK LOOP IS ASSUMING MUDST INPUT!!!!" << endl;
  Int_t nTrack = mHiMicroEvent->NTrack();
  StHiMicroTrack* track;
  
  Float_t vertexZ = mHiMicroEvent->VertexZ();

  if(mDebug)
    cout << nTrack << " tracks in this event" << endl;

  for(Int_t i=0; i<nTrack; i++){
    track =(StHiMicroTrack*) mHiMicroEvent->tracks()->At(i);

    //Need to set the sector variable for the tracks if we are analysing
    //MuDst files...
    track->SetFirstSector( findSector(track->PhiGl(),track->FirstZ()) );
    track->SetLastSector( findSector(track->PhiGl(),track->LastZ()) );
    if(mDebug)
      cout << "FirstSector : " << track->FirstSector() << ", LastSector: " << track->LastSector() << endl;
    

    //***************** spectra  ***********************************
    
    if(!CutRc::AcceptTrackHalf(track)) {
        continue; }
    if(!CutRc::Accept(track)) continue; 

    if(mDebug)
      cout << "Track accepted" << endl;

    //**************************************************************
    Float_t ptPr = track->PtPr();
    Float_t etaPr = track->EtaPr();
    Float_t ptGl = track->PtGl();
    Float_t etaGl = track->EtaGl();
    //
    // efficiency correction here
    //
    Float_t effPr = mEfficiencyMap->Eval(etaPr,ptPr);
    Float_t effGl = mEfficiencyMap->Eval(etaGl,ptGl);
    
    Float_t weightPr = (1./effPr);
    Float_t weightGl = (1./effGl);

    //
    // yield for both charge signs

    h2YieldPr->Fill(etaPr,ptPr);
    h2YieldGl->Fill(etaGl,ptGl);

    //
    // counts
    //
    if(ptPr>1.5&&ptPr<6){
      mNHiPtTrack++;
      if(mNHiPtTrack%10000==0){
//	cout << "ptPr: " << ptPr << " etaPr: " << etaPr << " effPr : " << effPr << endl;
      }
    }
    //
    // spectra stuff here
    //

    Int_t iCharge = (track->Charge()>0) ? 0 : 1; //plus is 0

    for(Int_t iBin=0; iBin<mNVarBin; iBin++){
      
      //
      // both
      //
      if (iBin == 2) {
	ptGl = 0.01*ptGl; //Scale to xT
        ptPr = 0.01*ptPr; //Scale to xT
      }

      varBin[iBin].specPM[2].h1RawGl->Fill(ptGl);
      varBin[iBin].meanPM[2].h1OneOverPtGl->Fill(ptGl,1./ptGl);
      varBin[iBin].specPM[2].h1EffCorrectedGl->Fill(ptGl,weightGl);
      varBin[iBin].specPM[2].h1CorrectedGl->Fill(ptGl,weightGl);

      varBin[iBin].specPM[2].h1RawPr->Fill(ptPr);
      varBin[iBin].meanPM[2].h1OneOverPtPr->Fill(ptPr,1./ptPr);
      varBin[iBin].specPM[2].h1EffCorrectedPr->Fill(ptPr,weightPr);
      varBin[iBin].specPM[2].h1CorrectedPr->Fill(ptPr,weightPr);
     

      //
      // plus/minus
      //

      varBin[iBin].specPM[iCharge].h1RawGl->Fill(ptGl);
      varBin[iBin].meanPM[iCharge].h1OneOverPtGl->Fill(ptGl,1./ptGl);
      varBin[iBin].specPM[iCharge].h1EffCorrectedGl->Fill(ptGl,weightGl);
      varBin[iBin].specPM[iCharge].h1CorrectedGl->Fill(ptGl,weightGl);

      varBin[iBin].specPM[iCharge].h1RawPr->Fill(ptPr);
      varBin[iBin].meanPM[iCharge].h1OneOverPtPr->Fill(ptPr,1./ptPr);
      varBin[iBin].specPM[iCharge].h1EffCorrectedPr->Fill(ptPr,weightPr);
      varBin[iBin].specPM[iCharge].h1CorrectedPr->Fill(ptPr,weightPr);


    } // varbin

  }
}

//_____________________

void
StHiSpectra::fillEventHistograms()
{
  Int_t cent = flowCentrality(mHiMicroEvent->NUncorrectedPrimaries());
  h1Centrality->Fill(cent);

  h1VertexZ->Fill(mHiMicroEvent->VertexZ());
		  
  if(CutRc::AcceptVertexZ(mHiMicroEvent))
    h1CentralityCut->Fill(cent);

}


//______________________

Bool_t
StHiSpectra::acceptEvent(StHiMicroEvent* event)
{
  return CutRc::Accept(event);

}


//______________________

void
StHiSpectra::finishHistograms()
{
  h1NEvent->SetBinContent(1,mNEventAccepted);
  h1EtaCut->SetBinContent(1,CutRc::mEta[0]);
  h1EtaCut->SetBinContent(2,CutRc::mEta[1]);

  for(Int_t iBin=0; iBin<mNVarBin; iBin++){

    for(Int_t iCharge=0; iCharge<3; iCharge++){

      varBin[iBin].meanPM[iCharge].h1WeightedMeanPr->Divide(varBin[iBin].specPM[iCharge].h1RawPr,varBin[iBin].meanPM[iCharge].h1OneOverPtPr);
      varBin[iBin].meanPM[iCharge].h1WeightedMeanGl->Divide(varBin[iBin].specPM[iCharge].h1RawGl,varBin[iBin].meanPM[iCharge].h1OneOverPtGl);

    }
  }

}

//______________________________________________
//This method added to compute the first and last sector on a track...
//JLK 19-Nov-2002

Int_t StHiSpectra::findSector(float phi, float z) {  

  Int_t clock = 0;
  Float_t pi = (Float_t) TMath::Pi();
  Float_t sectorWidth = pi/6.;

  Int_t sector[] = {8,7,6,5,4,3,2,1,12,11,10};

  if(phi>=(pi-sectorWidth/2.) || phi < (-pi+sectorWidth/2.)){
      clock = 9;
  }
  for(int i=0; i<11; i++){
    if(phi>=(-pi +sectorWidth*(i+1) - sectorWidth/2.) &&
       phi<(-pi +sectorWidth*(i+1) + sectorWidth/2.)){

        clock = sector[i];
        break;
    }
  }
  if(!clock){
   cerr << "Huh? phi = " << phi << endl; exit(1);
  }
  //Need to verify that this gets the right sector on the other side...
  //JLK 19-Nov-2002
  if (z < 0) clock+=12;   

  return clock;
   
}

//_____________________________________

ClassImp(StHiSpectra)