#include "TChain.h"
#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliESDPmdTrack.h"
#include "AliESDVertex.h"
#include "AliAnalysisTaskPMDSim.h"
#include <TParticle.h>
#include <AliLog.h>
#include <AliStack.h>
#include <AliHeader.h>
#include <AliGenEventHeader.h>
#include <AliHeader.h>
#include <AliGenPythiaEventHeader.h>
#include <AliMCEventHandler.h>
#include <AliMCEvent.h>
#include "AliTriggerAnalysis.h"
#include "AliPhysicsSelection.h"

// AnalysisTask For PMD
// Authors: Sudipan De, Subhash Singha

ClassImp(AliAnalysisTaskPMDSim)

//________________________________________________________________________
AliAnalysisTaskPMDSim::AliAnalysisTaskPMDSim(const char *name) 
  : AliAnalysisTaskSE(name), 
  fESD(0), 
  fOutputList(0), 
  fHistTotEvent(0),
  fHistTotEventAfterPhySel(0),
  fHistTotEventAfterVtx(0),
  fVtxZ(0),
  fHistXYPre(0),
  fHistEtaPhM(0),
  fHistEtaPhM1(0),
  fHistEtaT(0),
  fMultMeasured(0),
  fMultMeasured1(0),
  fMultTrue(0),
  fMultCorr(0),
  fMultCorr1(0)
{
  for(Int_t i=0; i<10; i++){
    fHistMultMeasEtaBinA[i] = 0;
    fHistMultMeasEtaBinA1[i] = 0;
    fHistMultTrueEtaBinA[i] = 0;
    fHistMultCorrEtaBinA[i] = 0;
    fHistMultCorrEtaBinA1[i] = 0;
  }
  
  // Constructor
  
  // Define input and output slots here
  // Input slot #0 works with a TChain
  DefineInput(0, TChain::Class());
  // Output slot #0 id reserved by the base class for AOD
  // Output slot #1 writes into a TH1 container
  DefineOutput(1, TList::Class());
}

//________________________________________________________________________
void AliAnalysisTaskPMDSim::UserCreateOutputObjects()
{
  // Create histograms
  // Called once
  
  fOutputList = new TList();
  
  Int_t kNbinsMult = 100;  Float_t XminMult = -0.5; 
  Float_t XmaxMult = 99.5;//mult
  Int_t kNbinsMultA = 50;  Float_t XminMultA = -0.5; 
  Float_t XmaxMultA = 49.5;//mult
  Int_t kNbinEtaPMDCov = 10; Float_t XminEtaPMDCov = 2.1; Float_t XmaxEtaPMDCov = 4.1;//etaPMD
  Int_t kNBinsEvent = 10; Float_t XminEvent = 0; Float_t XmaxEvent = 10;
  Int_t kNbinsXY = 200; Float_t XminXY = -100.0; Float_t XmaxXY  = 100.0;
  fHistTotEvent = new TH1F("TotEvent","TotEvent",
			   kNBinsEvent,XminEvent,XmaxEvent); 
  fOutputList->Add(fHistTotEvent);
  fHistTotEventAfterPhySel = new TH1F("TotEventAfterPhySel","TotEventAfterPhySel",
				      kNBinsEvent,XminEvent,XmaxEvent); 
  fOutputList->Add(fHistTotEventAfterPhySel);
  fHistTotEventAfterVtx = new TH1F("TotEventAfterVtx","TotEventAfterVtx",
				  kNBinsEvent,XminEvent,XmaxEvent); 
  fOutputList->Add(fHistTotEventAfterVtx);
  fVtxZ     = new TH1F("VtxZ","VtxZ",200,-50.0,50.0);
  fOutputList->Add(fVtxZ);
  fHistXYPre = new TH2F("XYPre","XYPre",kNbinsXY,XminXY,XmaxXY,
                        kNbinsXY,XminXY,XmaxXY);
  fOutputList->Add(fHistXYPre);
  fHistEtaPhM = new TH1F("fHistEtaPhM","fHistEtaPhM",kNbinEtaPMDCov,XminEtaPMDCov,XmaxEtaPMDCov);
  fOutputList->Add(fHistEtaPhM);
  fHistEtaPhM1 = new TH1F("fHistEtaPhM1","fHistEtaPhM1",kNbinEtaPMDCov,XminEtaPMDCov,XmaxEtaPMDCov);
  fOutputList->Add(fHistEtaPhM1);
  fHistEtaT = new TH1F("fHistEtaT","fHistEtaT",kNbinEtaPMDCov,XminEtaPMDCov,XmaxEtaPMDCov);
  fOutputList->Add(fHistEtaT);
  fMultMeasured = new TH1F("MultM","MultM",kNbinsMult,XminMult,XmaxMult);
  fOutputList->Add(fMultMeasured);
  fMultMeasured1 = new TH1F("MultM1","MultM1",kNbinsMult,XminMult,XmaxMult);
  fOutputList->Add(fMultMeasured1);
  fMultTrue = new TH1F("MultT","MultT",kNbinsMult,XminMult,XmaxMult);
  fOutputList->Add(fMultTrue);
  fMultCorr = new TH2F("MultCorr","MultCorr",kNbinsMult,XminMult,XmaxMult,
  		       kNbinsMult,XminMult,XmaxMult);
  fOutputList->Add(fMultCorr);
  fMultCorr1 = new TH2F("MultCorr1","MultCorr1",kNbinsMult,XminMult,XmaxMult,
  		       kNbinsMult,XminMult,XmaxMult);
  fOutputList->Add(fMultCorr1);
   
  
  Char_t nameT[256], nameM[256], nameCorr[256],nameM1[256], nameCorr1[256];
  for(Int_t i=0; i<10; i++){
    sprintf(nameM,"MultM_EtaBin%d",i+1);
    sprintf(nameM1,"MultM1_EtaBin%d",i+1);
    sprintf(nameT,"MultT_EtaBin%d",i+1);
    sprintf(nameCorr,"MultCorr_EtaBin%d",i+1);
    sprintf(nameCorr1,"MultCorr1_EtaBin%d",i+1);
    fHistMultMeasEtaBinA[i] = 
      new TH1F(nameM,nameM,kNbinsMultA,XminMultA,XmaxMultA);
    fHistMultMeasEtaBinA1[i] = 
      new TH1F(nameM1,nameM1,kNbinsMultA,XminMultA,XmaxMultA);
    fHistMultTrueEtaBinA[i] = 
      new TH1F(nameT,nameT,kNbinsMultA,XminMultA,XmaxMultA);
    fHistMultCorrEtaBinA[i] = 
      new TH2F(nameCorr,nameCorr,kNbinsMultA,XminMultA,XmaxMultA,
	       kNbinsMultA,XminMultA,XmaxMultA);
    fHistMultCorrEtaBinA1[i] = 
      new TH2F(nameCorr1,nameCorr1,kNbinsMultA,XminMultA,XmaxMultA,
	       kNbinsMultA,XminMultA,XmaxMultA);
    fOutputList->Add(fHistMultMeasEtaBinA[i]);
    fOutputList->Add(fHistMultMeasEtaBinA1[i]);
    fOutputList->Add(fHistMultTrueEtaBinA[i]);
    fOutputList->Add(fHistMultCorrEtaBinA[i]);
    fOutputList->Add(fHistMultCorrEtaBinA1[i]);
  }//i loop
  //fOutputList->Add(fHistPt);
}

//________________________________________________________________________
void AliAnalysisTaskPMDSim::UserExec(Option_t *) 
{
  // Main loop
  // Called for each event
  Float_t MipCut1 = 432;//MPV=432=> 6*MPV=432
  Float_t MipCut2 = 648;//MPV=432=> 6*MPV=648 
  Float_t etacls=0., theta=0., rdist=0.;
  Int_t PhotonCls = 0;
  Int_t PhotonCls1 = 0;
  Int_t PhotonClsAEtaBin[10] = {0};//#of photon measured with diff Eta bin
  Int_t PhotonClsAEtaBin1[10] = {0};//#of photon measured with diff Eta bin
  Int_t PhotonTrueAEtaBin[10] = {0};//#of photon Incident with diff Eta bin
  // Post output data.
  fESD = dynamic_cast<AliESDEvent*>(InputEvent());
  if (!fESD) {
    printf("ERROR: fESD not available\n");
    return;
  }
  fHistTotEvent->Fill(5);
  //event selection
  Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
  if (! isSelected) return;
  fHistTotEventAfterPhySel->Fill(5);
  //Vertex selection
  const AliESDVertex *vertex = fESD->GetPrimaryVertex();
  Float_t Vz = vertex->GetZ();    
  // Float_t Vx = vertex->GetX();    
  //Float_t Vy = vertex->GetY();
  Bool_t zVerStatus = vertex->GetStatus();
  if(zVerStatus){
    fVtxZ->Fill(Vz);
    if(TMath::Abs(Vz)<10){   
      fHistTotEventAfterVtx->Fill(5);   
      Int_t ptracks = fESD->GetNumberOfPmdTracks();
      for(Int_t kk=0;kk<ptracks;kk++){
	AliESDPmdTrack *pmdtr = fESD->GetPmdTrack(kk);
	Int_t   det   = pmdtr->GetDetector();
	Float_t clsX  = pmdtr->GetClusterX();
	Float_t clsY  = pmdtr->GetClusterY();
	Float_t clsZ  = pmdtr->GetClusterZ();
	clsZ -= Vz;
	Float_t ncell = pmdtr->GetClusterCells();
	Float_t adc   = pmdtr->GetClusterADC();
	//Float_t pid   = pmdtr->GetClusterPID();
	//Float_t isotag = pmdtr->GetClusterSigmaX();
	//Int_t trno = pmdtr->GetClusterTrackNo();
	//calculation of #eta
	rdist = TMath::Sqrt(clsX*clsX + clsY*clsY);
	if(clsZ!=0) theta = TMath::ATan2(rdist,clsZ);
	etacls   = -TMath::Log(TMath::Tan(0.5*theta));
	
	if(det==0 && adc>0)fHistXYPre->Fill(clsX,clsY);
	if(det==0 && adc>MipCut1 && ncell>2){
	  fHistEtaPhM->Fill(etacls);
	  if(etacls > 2.3 && etacls <= 3.9)  PhotonCls++;
	  if(etacls > 2.1 && etacls <= 2.3)  PhotonClsAEtaBin[0]++;
	  if(etacls > 2.3 && etacls <= 2.5)  PhotonClsAEtaBin[1]++;
	  if(etacls > 2.5 && etacls <= 2.7)  PhotonClsAEtaBin[2]++;
	  if(etacls > 2.7 && etacls <= 2.9)  PhotonClsAEtaBin[3]++;
	  if(etacls > 2.9 && etacls <= 3.1)  PhotonClsAEtaBin[4]++;
	  if(etacls > 3.1 && etacls <= 3.3)  PhotonClsAEtaBin[5]++;
	  if(etacls > 3.3 && etacls <= 3.5)  PhotonClsAEtaBin[6]++;
	  if(etacls > 3.5 && etacls <= 3.7)  PhotonClsAEtaBin[7]++;
	  if(etacls > 3.7 && etacls <= 3.9)  PhotonClsAEtaBin[8]++;
	  if(etacls > 3.9 && etacls <= 4.1)  PhotonClsAEtaBin[9]++;
	}//if MipCut1
	if( det==0 && adc>MipCut2 && ncell>2){
	  fHistEtaPhM1->Fill(etacls);
	  if(etacls > 2.3 && etacls <= 3.9)  PhotonCls1++;
	  if(etacls > 2.1 && etacls <= 2.3)  PhotonClsAEtaBin1[0]++;
	  if(etacls > 2.3 && etacls <= 2.5)  PhotonClsAEtaBin1[1]++;
	  if(etacls > 2.5 && etacls <= 2.7)  PhotonClsAEtaBin1[2]++;
	  if(etacls > 2.7 && etacls <= 2.9)  PhotonClsAEtaBin1[3]++;
	  if(etacls > 2.9 && etacls <= 3.1)  PhotonClsAEtaBin1[4]++;
	  if(etacls > 3.1 && etacls <= 3.3)  PhotonClsAEtaBin1[5]++;
	  if(etacls > 3.3 && etacls <= 3.5)  PhotonClsAEtaBin1[6]++;
	  if(etacls > 3.5 && etacls <= 3.7)  PhotonClsAEtaBin1[7]++;
	  if(etacls > 3.7 && etacls <= 3.9)  PhotonClsAEtaBin1[8]++;
	  if(etacls > 3.9 && etacls <= 4.1)  PhotonClsAEtaBin1[9]++;
	}//if MipCut2
      } //track loop 
      //reading MC info.
      AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*>
	(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
      if (!eventHandler) {
	Printf("ERROR: Could not retrieve MC event handler");
	return;
      }
      AliMCEvent* mcEvent = eventHandler->MCEvent();
      if (!mcEvent) {
	Printf("ERROR: Could not retrieve MC event");
	return;
      }
      AliStack* stack = mcEvent->Stack();
      if (!stack)
	{
	  AliDebug(AliLog::kError, "Stack not available");
	  return;
	}
      if(stack){
	Int_t nPrim  = stack->GetNprimary();
	Int_t PhotonTrue = 0;
	for (Int_t iMc = 0; iMc < nPrim; ++iMc)
	  {
	    TParticle *MPart = stack->Particle(iMc);
	    Int_t mpart  = MPart->GetPdgCode();
	    Float_t eta   = MPart->Eta();
	    if(mpart == 22){  
	      fHistEtaT->Fill(eta);
	      if(eta > 2.3 && eta <= 3.9)PhotonTrue++;
	      if(eta > 2.1 && eta <= 2.3)PhotonTrueAEtaBin[0]++;
	      if(eta > 2.3 && eta <= 2.5)PhotonTrueAEtaBin[1]++;
	      if(eta > 2.5 && eta <= 2.7)PhotonTrueAEtaBin[2]++;
	      if(eta > 2.7 && eta <= 2.9)PhotonTrueAEtaBin[3]++;
	      if(eta > 2.9 && eta <= 3.1)PhotonTrueAEtaBin[4]++;
	      if(eta > 3.1 && eta <= 3.3)PhotonTrueAEtaBin[5]++;
	      if(eta > 3.3 && eta <= 3.5)PhotonTrueAEtaBin[6]++;
	      if(eta > 3.5 && eta <= 3.7)PhotonTrueAEtaBin[7]++;
	      if(eta > 3.7 && eta <= 3.9)PhotonTrueAEtaBin[8]++;
	      if(eta > 3.9 && eta <= 4.1)PhotonTrueAEtaBin[9]++;
	    }//mpart
	  }//track loop
	for(Int_t i=0; i<10; i++){
	  fHistMultMeasEtaBinA[i]->Fill(PhotonClsAEtaBin[i]);
	  fHistMultMeasEtaBinA1[i]->Fill(PhotonClsAEtaBin1[i]);
	  fHistMultTrueEtaBinA[i]->Fill(PhotonTrueAEtaBin[i]);
	  fHistMultCorrEtaBinA[i]->Fill(PhotonTrueAEtaBin[i],PhotonClsAEtaBin[i]);
	  fHistMultCorrEtaBinA1[i]->Fill(PhotonTrueAEtaBin[i],PhotonClsAEtaBin1[i]);
	}//i loop
	fMultMeasured->Fill(PhotonCls);
	fMultMeasured1->Fill(PhotonCls1);
	fMultTrue->Fill(PhotonTrue);
	fMultCorr->Fill(PhotonTrue,PhotonCls);
	fMultCorr1->Fill(PhotonTrue,PhotonCls1);
      }//if stack
    }//if Vz!= 0
  }//Vz loop
  PostData(1, fOutputList);
}//UserExec()     

//_______________________________________________________________________
void AliAnalysisTaskPMDSim::Terminate(Option_t *) 
{
  // Draw result to the screen
  // Called once at the end of the query
  fOutputList = dynamic_cast<TList*> (GetOutputData(1));
  if (!fOutputList) {
    printf("ERROR: Output list not available\n");
    return;
  }
}//Terminate