/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//==============================================================================
// AliHMPIDTaskQA - Class representing a quality check tool of HMPID 
// A set of histograms is created.
//==============================================================================

#ifndef AliHMPIDTaskQA_CXX
#define AliHMPIDTaskQA_CXX


#include "TH1.h"
#include "TH2.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TGraphErrors.h"
#include "AliAnalysisManager.h"
#include "AliESDInputHandler.h"
#include "AliMCEventHandler.h"
#include "AliMCEvent.h"
#include "AliESDtrack.h"
#include "AliPID.h"
#include "AliLog.h"
#include "AliHMPIDTaskQA.h"

ClassImp(AliHMPIDTaskQA)

//__________________________________________________________________________
AliHMPIDTaskQA::AliHMPIDTaskQA() :
  fESD(0x0),fMC(0x0),fUseMC(kTRUE),
  fHmpHistList(0x0),
  fHmpPesdPhmp(0x0),fHmpCkovPesd(0x0),fHmpCkovPhmp(0x0),
  fHmpMipCharge3cm(0x0),fHmpMipTrkDist(0x0),
  fHmpTrkFlags(0x0),
  fN1(6),
  fN2(8),
  fPionEff(0x0),
  fKaonEff(0x0),
  fProtEff(0x0),
  fPionTot(0x0),
  fKaonTot(0x0),
  fProtTot(0x0),
  fPionNot(0x0),
  fKaonNot(0x0),
  fProtNot(0x0),
  fPionCon(0x0),
  fKaonCon(0x0),
  fProtCon(0x0)
{
  //
  //Default ctor
  //
  for (Int_t i=0; i<7; i++){
    fHmpPhotons[i]    = 0x0;
    fHmpPhotP[i]      = 0x0;
    fHmpPhotSin2th[i] = 0x0;
    fHmpMipTrkDistPosX[i] = fHmpMipTrkDistNegX[i] = 0x0;
    fHmpMipTrkDistPosY[i] = fHmpMipTrkDistNegY[i] = 0x0;
    fHmpMipCharge[i] = 0x0;
  }
}

//___________________________________________________________________________
AliHMPIDTaskQA::AliHMPIDTaskQA(const Char_t* name) :
  AliAnalysisTaskSE(name),
  fESD(0x0), fMC(0x0), fUseMC(kTRUE),
  fHmpHistList(0x0),
  fHmpPesdPhmp(0x0),fHmpCkovPesd(0x0),fHmpCkovPhmp(0x0),
  fHmpMipCharge3cm(0x0),fHmpMipTrkDist(0x0),
  fHmpTrkFlags(0x0),
  fN1(6),
  fN2(8),
  fPionEff(0x0),
  fKaonEff(0x0),
  fProtEff(0x0),
  fPionTot(0x0),
  fKaonTot(0x0),
  fProtTot(0x0),
  fPionNot(0x0),
  fKaonNot(0x0),
  fProtNot(0x0),
  fPionCon(0x0),
  fKaonCon(0x0),
  fProtCon(0x0)
{
  //
  // Constructor. Initialization of Inputs and Outputs
  //
  for (Int_t i=0; i<7; i++){
    fHmpPhotons[i] = 0x0;
    fHmpPhotP[i]   = 0x0;
    fHmpPhotSin2th[i] = 0x0;
    fHmpMipTrkDistPosX[i] = fHmpMipTrkDistNegX[i] = 0x0;
    fHmpMipTrkDistPosY[i] = fHmpMipTrkDistNegY[i] = 0x0;
    fHmpMipCharge[i] = 0x0;
  }

  DefineOutput(1,TList::Class());
}

//___________________________________________________________________________
AliHMPIDTaskQA& AliHMPIDTaskQA::operator=(const AliHMPIDTaskQA& c)
{
  //
  // Assignment operator
  //
  if (this!=&c) {
    AliAnalysisTaskSE::operator=(c);
    fESD             = c.fESD;
    fMC              = c.fMC;
    fUseMC           = c.fUseMC;
    fHmpHistList     = c.fHmpHistList;
    fHmpPesdPhmp     = c.fHmpPesdPhmp;
    fHmpCkovPesd     = c.fHmpCkovPesd;
    fHmpCkovPhmp     = c.fHmpCkovPhmp;
    fHmpMipCharge3cm = c.fHmpMipCharge3cm;
    fHmpMipTrkDist   = c.fHmpMipTrkDist;
    fHmpTrkFlags     = c.fHmpTrkFlags;
    fN1              = c.fN1;
    fN2              = c.fN2;
    fPionEff         = c.fPionEff;
    fKaonEff         = c.fKaonEff;
    fProtEff         = c.fProtEff;
    fPionTot         = c.fPionTot;
    fKaonTot         = c.fKaonTot;
    fProtTot         = c.fProtTot;
    fPionNot         = c.fPionNot;
    fKaonNot         = c.fKaonNot;
    fProtNot         = c.fProtNot;
    fPionCon         = c.fPionCon;
    fKaonCon         = c.fKaonCon;
    fProtCon         = c.fProtCon;
    for (Int_t i=0; i<7; i++){
      fHmpPhotons[i] = c.fHmpPhotons[i];
      fHmpPhotP[i]   = c.fHmpPhotP[i];
      fHmpPhotSin2th[i] = c.fHmpPhotSin2th[i];
      fHmpMipTrkDistPosX[i] = c.fHmpMipTrkDistPosX[i];
      fHmpMipTrkDistNegX[i] = c.fHmpMipTrkDistNegX[i];
      fHmpMipTrkDistPosY[i] = c.fHmpMipTrkDistPosY[i];
      fHmpMipTrkDistNegY[i] = c.fHmpMipTrkDistNegY[i];
      fHmpMipCharge[i] = c.fHmpMipCharge[i];
    }
  }
  return *this;
}

//___________________________________________________________________________
AliHMPIDTaskQA::AliHMPIDTaskQA(const AliHMPIDTaskQA& c) :
  AliAnalysisTaskSE(c),
  fESD(c.fESD),fMC(c.fMC),fUseMC(c.fUseMC),
  fHmpHistList(c.fHmpHistList),
  fHmpPesdPhmp(c.fHmpPesdPhmp),fHmpCkovPesd(c.fHmpCkovPesd),fHmpCkovPhmp(c.fHmpCkovPhmp),
  fHmpMipCharge3cm(c.fHmpMipCharge3cm),fHmpMipTrkDist(c.fHmpMipTrkDist),
  fHmpTrkFlags(c.fHmpTrkFlags),
  fN1(c.fN1),
  fN2(c.fN2),
  fPionEff(c.fPionEff),
  fKaonEff(c.fKaonEff),
  fProtEff(c.fProtEff),
  fPionTot(c.fPionTot),
  fKaonTot(c.fKaonTot),
  fProtTot(c.fProtTot),
  fPionNot(c.fPionNot),
  fKaonNot(c.fKaonNot),
  fProtNot(c.fProtNot),
  fPionCon(c.fPionCon),
  fKaonCon(c.fKaonCon),
  fProtCon(c.fProtCon)
{
  //
  // Copy Constructor
  //
  for (Int_t i=0; i<7; i++){
    fHmpPhotons[i] = c.fHmpPhotons[i];
    fHmpPhotP[i]   = c.fHmpPhotP[i];
    fHmpPhotSin2th[i] = c.fHmpPhotSin2th[i];
    fHmpMipTrkDistPosX[i] = c.fHmpMipTrkDistPosX[i];
    fHmpMipTrkDistNegX[i] = c.fHmpMipTrkDistNegX[i];
    fHmpMipTrkDistPosY[i] = c.fHmpMipTrkDistPosY[i];
    fHmpMipTrkDistNegY[i] = c.fHmpMipTrkDistNegY[i];
    fHmpMipCharge[i] = c.fHmpMipCharge[i];
  }
}
 
//___________________________________________________________________________
AliHMPIDTaskQA::~AliHMPIDTaskQA() {
  //
  //destructor
  //
  Info("~AliHMPIDTaskQA","Calling Destructor");
  if (fHmpHistList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fHmpHistList;
}

//___________________________________________________________________________
void AliHMPIDTaskQA::ConnectInputData(Option_t *option)
{
  AliAnalysisTaskSE::ConnectInputData(option);

  // Connect ESD here
  AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
  if (!esdH) {
    AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
  } else
    fESD = esdH->GetEvent();

  if (fUseMC){
    // Connect MC
    AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
    if (!mcH) {
      AliDebug(2,Form("ERROR: Could not get MCEventHandler"));
      fUseMC = kFALSE;
    }
  }
}

//___________________________________________________________________________
void AliHMPIDTaskQA::UserExec(Option_t *)
{
  Double_t priors[10]={1.,1.,1.,1.,1.,0.,0.,0.,0.,0.}; //{0.01,0.01,0.83,0.10,0.5};
  Double_t probs[5];
  AliPID *pPid = new AliPID();
  pPid->SetPriors(priors);
  AliESDtrack *track=0;
  TParticle *pPart=0;
  AliStack* pStack = 0;
  Int_t label = -1;
  if (fUseMC){
    AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
    if (mcH) {
      fMC = mcH->MCEvent();
      pStack = fMC->Stack();
    } else {
      AliFatal("fMC flag set but MC handler not available");
      return;
    }  
  }

  //
  // Main loop function, executed on Event basis
  //
  
  if( !fESD->GetPrimaryVertex()) return;
  if( fESD->GetPrimaryVertex()->GetNContributors() < 1 ) return;
  if( TMath::Abs(fESD->GetPrimaryVertex()->GetZ()) > 10.0 /* cm */) return;
  
  
  for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTracks(); iTrack++) {

    track = fESD->GetTrack(iTrack);
    if(!track) continue;
    if(! track->IsOn(AliESDtrack::kITSrefit) || !track->IsOn(AliESDtrack::kITSrefit) ) continue;
     
    Double_t rin[3], rout[3];
    track->GetInnerXYZ(rin);
    track->GetOuterXYZ(rout);
    Double_t ktol = 0.001;
    Double_t thetaCh = track->GetHMPIDsignal();

    if(Equal(thetaCh,-20.,ktol))      fHmpTrkFlags->Fill(0);
    else if(Equal(thetaCh,-9.,ktol))  fHmpTrkFlags->Fill(1);
    else if(Equal(thetaCh,-5.,ktol))  fHmpTrkFlags->Fill(2);
    else if(Equal(thetaCh,-11.,ktol)) fHmpTrkFlags->Fill(3);
    else if(Equal(thetaCh,-22.,ktol)) fHmpTrkFlags->Fill(4);
    else fHmpTrkFlags->Fill(4);

    if(Equal(thetaCh,-20.,ktol)) continue;
    if(track->GetHMPIDcluIdx() < 0) continue;

    Int_t q, nph;
    Float_t x, y;
    Float_t xpc, ypc, th, ph;
    track->GetHMPIDmip(x,y,q,nph);
    track->GetHMPIDtrk(xpc,ypc,th,ph);

    if(Equal(x,0.,ktol) && Equal(y,0.,ktol) && Equal(xpc,0.,ktol) && Equal(ypc,0.,ktol)) continue;

    Float_t sign = track->GetSign();
    Int_t ch = track->GetHMPIDcluIdx()/1000000; 
    if( ch < 0 || ch > 6 ) continue; //Chamber numbering is [0,6]
    
    if (sign > 0.) fHmpMipTrkDistPosX[ch]->Fill(xpc - x), fHmpMipTrkDistPosY[ch]->Fill(ypc - y);
    if (sign < 0.) fHmpMipTrkDistNegX[ch]->Fill(xpc - x), fHmpMipTrkDistNegY[ch]->Fill(ypc - y);
    fHmpMipCharge[ch]->Fill(q);
    Double_t dist = TMath::Sqrt((xpc - x)*(xpc - x) + (ypc - y)*(ypc - y));
    if (dist <= 3.0) fHmpMipCharge3cm->Fill(q);
    fHmpMipTrkDist->Fill(dist);

    track->GetHMPIDpid(probs);
    pPid->SetProbabilities(probs);
    if (fUseMC){
      if ((label = track->GetLabel()) < 0) continue;
      pPart = pStack->Particle(label);
    }

    if(thetaCh > 0 ){
      Double_t pHmp[3] = {0}, pHmp3 = 0;
      if (track->GetOuterHmpPxPyPz(pHmp)) pHmp3 = TMath::Sqrt(pHmp[0]*pHmp[0]+pHmp[1]*pHmp[1]+pHmp[2]*pHmp[2]);
      if (pHmp3) fHmpPesdPhmp->Fill(track->P(),pHmp3);
      fHmpCkovPesd->Fill(track->P(),thetaCh);
      if (pHmp3) fHmpCkovPhmp->Fill(pHmp3,thetaCh);
      fHmpPhotons[ch]->Fill(nph);
      fHmpPhotP[ch]->Fill(pHmp3,nph);
      Double_t sin2 = TMath::Power(TMath::Sin(th),2);
      fHmpPhotSin2th[ch]->Fill(sin2,nph);

      if (fUseMC && dist<0.5 && TMath::Abs(th)<0.13){
        if (!pStack->IsPhysicalPrimary(label)) continue;
        if (track->Pt()<1. || track->Pt()>5.) continue;
        Int_t pdgCode = TMath::Abs(pPart->GetPdgCode());
        Int_t ptBin=(Int_t) (2*(track->Pt()-1));
        if (pdgCode!=2212) fProtCon->Fill(ptBin);
        if (pdgCode==2212){
          fProtTot->Fill(ptBin);
          fProtEff->Fill(ptBin,pPid->GetProbability(AliPID::kProton));
          fPionNot->Fill(ptBin,pPid->GetProbability(AliPID::kPion));
          fKaonNot->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
        }
        if (pdgCode!=211) fPionCon->Fill(ptBin);
        if (pdgCode!=321) fKaonCon->Fill(ptBin);
        if (pdgCode==211){
          if (ptBin < 6){
            Float_t weight=pPid->GetProbability(AliPID::kElectron)+
                           pPid->GetProbability(AliPID::kMuon)+
                           pPid->GetProbability(AliPID::kPion);
            fPionTot->Fill(ptBin);
            fPionEff->Fill(ptBin,weight);
            fKaonNot->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
          }
          fProtNot->Fill(ptBin,pPid->GetProbability(AliPID::kProton));
        }
        if (pdgCode==321){
          if (ptBin < 6){
            fKaonTot->Fill(ptBin);
            fKaonEff->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
            fPionNot->Fill(ptBin,(pPid->GetProbability(AliPID::kPion)));
          }
          fProtNot->Fill(ptBin,(pPid->GetProbability(AliPID::kProton)));
        }
      }
    }//there is signal
  }//track loop
  delete pPid;

  /* PostData(0) is taken care of by AliAnalysisTaskSE */
  PostData(1,fHmpHistList);
}

//___________________________________________________________________________
void AliHMPIDTaskQA::Terminate(Option_t*)
{
  // The Terminate() function is the last function to be called during
  // a query. It always runs on the client, it can be used to present
  // the results graphically or save the results to file.

  Info("Terminate"," ");
  AliAnalysisTaskSE::Terminate();

}

//___________________________________________________________________________
void AliHMPIDTaskQA::UserCreateOutputObjects() {
  //
  //HERE ONE CAN CREATE OUTPUT OBJECTS
  //TO BE SET BEFORE THE EXECUTION OF THE TASK
  //

  //slot #1
//   OpenFile(1);
   fHmpHistList = new TList();
   fHmpHistList->SetOwner();

   fHmpPesdPhmp = new TH2F("fHmpPesdPhmp","HMPID: ESD p - running p;HMP p (GeV/c);ESD p (Gev/c)",100,0,10,100,0,10);
   fHmpHistList->Add(fHmpPesdPhmp);

   fHmpCkovPesd = new TH2F("fHmpCkovPesd","HMPID: ThetaCherenkov vs P;p_esd (GeV/c);#Theta_C;Entries",500,0,10,500,0,1.0);
   fHmpHistList->Add(fHmpCkovPesd);

   fHmpCkovPhmp = new TH2F("fHmpCkovPhmp","HMPID: ThetaCherenkov vs P;p_hmp (GeV/c);#Theta_C;Entries",500,0,10,500,0,1.0);
   fHmpHistList->Add(fHmpCkovPhmp);

   for (Int_t i=0; i<7; i++){
     TString title=Form("MIP-Track distance in local X, Chamber %d;distance (cm);Entries",i);
     fHmpMipTrkDistPosX[i] = new TH1F(Form("fHmpMipTrkDistPosX%d",i),title.Data(),800,-20,20);
     fHmpHistList->Add(fHmpMipTrkDistPosX[i]);
     fHmpMipTrkDistNegX[i] = new TH1F(Form("fHmpMipTrkDistNegX%d",i),title.Data(),800,-20,20);
     fHmpHistList->Add(fHmpMipTrkDistNegX[i]);

     title=Form("MIP-Track distance in local Y, Chamber %d;distance (cm);Entries",i);
     fHmpMipTrkDistPosY[i] = new TH1F(Form("fHmpMipTrkDistPosY%d",i),title.Data(),800,-20,20);
     fHmpHistList->Add(fHmpMipTrkDistPosY[i]);
     fHmpMipTrkDistNegY[i] = new TH1F(Form("fHmpMipTrkDistNegY%d",i),title.Data(),800,-20,20);
     fHmpHistList->Add(fHmpMipTrkDistNegY[i]);

     title=Form("Mip charge distribution, Chamber %d;Charge;Entries",i);
     fHmpMipCharge[i] = new TH1F(Form("fHmpMipCharge%d",i),title.Data(),5001,-0.5,5000.5);
     fHmpHistList->Add(fHmpMipCharge[i]);

     title=Form("Photons, Chamber %d;N of Photons;Entries",i);
     fHmpPhotons[i] = new TH1F(Form("fHmpPhotons%d",i),title.Data(),100,0,100);
     fHmpHistList->Add(fHmpPhotons[i]);

     title=Form("Photons versus HMP momentum, Chamber %d;P (GeV);N of Photons",i);
     fHmpPhotP[i] = new TH2F(Form("fHmpPhotP%d",i),title.Data(),200,0.,20.,100,0,100);
     fHmpHistList->Add(fHmpPhotP[i]);

     title=Form("Photons versus sin(th)^2 (uncorrected), Chamber %d;P (GeV);N of Photons",i);
     fHmpPhotSin2th[i] = new TH2F(Form("fHmpPhotSin2th%d",i),title.Data(),100,0.,1.,100,0,100);
     fHmpHistList->Add(fHmpPhotSin2th[i]);
   }

   fHmpMipCharge3cm = new TH1F("fHmpMipCharge3cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
   fHmpHistList->Add(fHmpMipCharge3cm);

   fHmpMipTrkDist = new TH1F("fHmpMipTrkDist","Mip-track distance in all the chambers",100,0.,10.);
   fHmpHistList->Add(fHmpMipTrkDist);

   fHmpTrkFlags = new TH1F("fHmpTrkFlags","HMPID track flags",6,0,6);
   TString summary[6] =  {"NotPerformed","MipDistCut", "MipQdcCut", "NoPhotAccept", "kNoRad", "other"};
   for(Int_t ibin = 0; ibin < 6; ibin++) fHmpTrkFlags->GetXaxis()->SetBinLabel(ibin+1,Form("%i  %s",ibin+1,summary[ibin].Data()));
   fHmpHistList->Add(fHmpTrkFlags);

   fPionEff = new TH1F("PionEff","Identified pions",fN1,0,fN1);
   fKaonEff = new TH1F("KaonEff","Identified kaons",fN1,0,fN1);
   fProtEff = new TH1F("ProtEff","Identified protons",fN2,0,fN2);
   fPionTot = new TH1I("PionTot","Total MC pions",fN1,0,fN1);
   fKaonTot = new TH1I("KaonTot","Total MC kaons",fN1,0,fN1);
   fProtTot = new TH1I("ProtTot","Total MC protons",fN2,0,fN2);
   fPionNot = new TH1F("PionNot","Misidentified pions",fN1,0,fN1);
   fKaonNot = new TH1F("KaonNot","Misidentified kaons",fN1,0,fN1);
   fProtNot = new TH1F("ProtNot","Misidentified protons",fN2,0,fN2);
   fPionCon = new TH1I("PionCon","Total not MC pions",fN1,0,fN1);
   fKaonCon = new TH1I("KaonCon","Total not MC kaons",fN1,0,fN1);
   fProtCon = new TH1I("ProtCon","Total not MC protons",fN2,0,fN2);

   fHmpHistList->Add(fPionEff); fHmpHistList->Add(fKaonEff); fHmpHistList->Add(fProtEff);
   fHmpHistList->Add(fPionTot); fHmpHistList->Add(fKaonTot); fHmpHistList->Add(fProtTot);
   fHmpHistList->Add(fPionNot); fHmpHistList->Add(fKaonNot); fHmpHistList->Add(fProtNot);
   fHmpHistList->Add(fPionCon); fHmpHistList->Add(fKaonCon); fHmpHistList->Add(fProtCon);

   PostData(1,fHmpHistList);
}

//____________________________________________________________________________________________________________________________________
Bool_t AliHMPIDTaskQA::Equal(Double_t x, Double_t y, Double_t tolerance)
{
 return abs(x - y) <= tolerance ;
}
   
#endif