[u/mrichter/AliRoot.git] / PWG1 / HMPID / AliHMPIDTaskQA.cxx

/**************************************************************************
 * 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;
    } else
      fMC = mcH->MCEvent();
      if (!fMC) AliDebug(2,Form("ERROR: Could not get MCEvent"));
  }
}

//___________________________________________________________________________
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){
    pStack = fMC->Stack();
  }

  //
  // 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",100,0,10,110,0,1.1);
   fHmpHistList->Add(fHmpCkovPesd);

   fHmpCkovPhmp = new TH2F("fHmpCkovPhmp","HMPID: ThetaCherenkov vs P;p_hmp (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
   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
Commit	Line	Data
3882540a	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//==============================================================================
	17	// AliHMPIDTaskQA - Class representing a quality check tool of HMPID
	18	// A set of histograms is created.
	19	//==============================================================================
	20
	21	#ifndef AliHMPIDTaskQA_CXX
	22	#define AliHMPIDTaskQA_CXX
	23
	24
	25	#include "TH1.h"
	26	#include "TH2.h"
	27	#include "TFile.h"
	28	#include "TCanvas.h"
	29	#include "TGraphErrors.h"
	30	#include "AliAnalysisManager.h"
	31	#include "AliESDInputHandler.h"
	32	#include "AliMCEventHandler.h"
	33	#include "AliMCEvent.h"
	34	#include "AliESDtrack.h"
	35	#include "AliPID.h"
	36	#include "AliLog.h"
	37	#include "AliHMPIDTaskQA.h"
	38
	39	ClassImp(AliHMPIDTaskQA)
	40
	41	//__________________________________________________________________________
	42	AliHMPIDTaskQA::AliHMPIDTaskQA() :
	43	fESD(0x0),fMC(0x0),fUseMC(kTRUE),
	44	fHmpHistList(0x0),
15939b4c	45	fHmpPesdPhmp(0x0),fHmpCkovPesd(0x0),fHmpCkovPhmp(0x0),
65896ecb	46	fHmpMipCharge3cm(0x0),fHmpMipTrkDist(0x0),
3882540a	47	fHmpTrkFlags(0x0),
	48	fN1(6),
	49	fN2(8),
	50	fPionEff(0x0),
	51	fKaonEff(0x0),
	52	fProtEff(0x0),
	53	fPionTot(0x0),
	54	fKaonTot(0x0),
	55	fProtTot(0x0),
	56	fPionNot(0x0),
	57	fKaonNot(0x0),
	58	fProtNot(0x0),
	59	fPionCon(0x0),
	60	fKaonCon(0x0),
15939b4c	61	fProtCon(0x0)
3882540a	62	{
	63	//
	64	//Default ctor
	65	//
65896ecb	66	for (Int_t i=0; i<7; i++){
	67	fHmpPhotons[i] = 0x0;
	68	fHmpPhotP[i] = 0x0;
	69	fHmpPhotSin2th[i] = 0x0;
	70	fHmpMipTrkDistPosX[i] = fHmpMipTrkDistNegX[i] = 0x0;
15939b4c	71	fHmpMipTrkDistPosY[i] = fHmpMipTrkDistNegY[i] = 0x0;
65896ecb	72	fHmpMipCharge[i] = 0x0;
65896ecb	73	}
3882540a	74	}
	75
	76	//___________________________________________________________________________
	77	AliHMPIDTaskQA::AliHMPIDTaskQA(const Char_t* name) :
	78	AliAnalysisTaskSE(name),
	79	fESD(0x0), fMC(0x0), fUseMC(kTRUE),
	80	fHmpHistList(0x0),
15939b4c	81	fHmpPesdPhmp(0x0),fHmpCkovPesd(0x0),fHmpCkovPhmp(0x0),
65896ecb	82	fHmpMipCharge3cm(0x0),fHmpMipTrkDist(0x0),
3882540a	83	fHmpTrkFlags(0x0),
	84	fN1(6),
	85	fN2(8),
	86	fPionEff(0x0),
	87	fKaonEff(0x0),
	88	fProtEff(0x0),
	89	fPionTot(0x0),
	90	fKaonTot(0x0),
	91	fProtTot(0x0),
	92	fPionNot(0x0),
	93	fKaonNot(0x0),
	94	fProtNot(0x0),
	95	fPionCon(0x0),
	96	fKaonCon(0x0),
15939b4c	97	fProtCon(0x0)
3882540a	98	{
	99	//
	100	// Constructor. Initialization of Inputs and Outputs
	101	//
65896ecb	102	for (Int_t i=0; i<7; i++){
	103	fHmpPhotons[i] = 0x0;
	104	fHmpPhotP[i] = 0x0;
	105	fHmpPhotSin2th[i] = 0x0;
	106	fHmpMipTrkDistPosX[i] = fHmpMipTrkDistNegX[i] = 0x0;
15939b4c	107	fHmpMipTrkDistPosY[i] = fHmpMipTrkDistNegY[i] = 0x0;
65896ecb	108	fHmpMipCharge[i] = 0x0;
65896ecb	109	}
3882540a	110
3882540a	111	DefineOutput(1,TList::Class());
3882540a	112	}
	113
	114	//___________________________________________________________________________
	115	AliHMPIDTaskQA& AliHMPIDTaskQA::operator=(const AliHMPIDTaskQA& c)
	116	{
	117	//
	118	// Assignment operator
	119	//
	120	if (this!=&c) {
	121	AliAnalysisTaskSE::operator=(c);
	122	fESD = c.fESD;
	123	fMC = c.fMC;
	124	fUseMC = c.fUseMC;
	125	fHmpHistList = c.fHmpHistList;
15939b4c	126	fHmpPesdPhmp = c.fHmpPesdPhmp;
	127	fHmpCkovPesd = c.fHmpCkovPesd;
	128	fHmpCkovPhmp = c.fHmpCkovPhmp;
3882540a	129	fHmpMipCharge3cm = c.fHmpMipCharge3cm;
65896ecb	130	fHmpMipTrkDist = c.fHmpMipTrkDist;
3882540a	131	fHmpTrkFlags = c.fHmpTrkFlags;
	132	fN1 = c.fN1;
	133	fN2 = c.fN2;
	134	fPionEff = c.fPionEff;
	135	fKaonEff = c.fKaonEff;
	136	fProtEff = c.fProtEff;
	137	fPionTot = c.fPionTot;
	138	fKaonTot = c.fKaonTot;
	139	fProtTot = c.fProtTot;
	140	fPionNot = c.fPionNot;
	141	fKaonNot = c.fKaonNot;
	142	fProtNot = c.fProtNot;
	143	fPionCon = c.fPionCon;
	144	fKaonCon = c.fKaonCon;
	145	fProtCon = c.fProtCon;
15939b4c	146	for (Int_t i=0; i<7; i++){
65896ecb	147	fHmpPhotons[i] = c.fHmpPhotons[i];
	148	fHmpPhotP[i] = c.fHmpPhotP[i];
	149	fHmpPhotSin2th[i] = c.fHmpPhotSin2th[i];
	150	fHmpMipTrkDistPosX[i] = c.fHmpMipTrkDistPosX[i];
	151	fHmpMipTrkDistNegX[i] = c.fHmpMipTrkDistNegX[i];
15939b4c	152	fHmpMipTrkDistPosY[i] = c.fHmpMipTrkDistPosY[i];
15939b4c	153	fHmpMipTrkDistNegY[i] = c.fHmpMipTrkDistNegY[i];
65896ecb	154	fHmpMipCharge[i] = c.fHmpMipCharge[i];
15939b4c	155	}
3882540a	156	}
	157	return *this;
	158	}
	159
	160	//___________________________________________________________________________
	161	AliHMPIDTaskQA::AliHMPIDTaskQA(const AliHMPIDTaskQA& c) :
	162	AliAnalysisTaskSE(c),
	163	fESD(c.fESD),fMC(c.fMC),fUseMC(c.fUseMC),
	164	fHmpHistList(c.fHmpHistList),
15939b4c	165	fHmpPesdPhmp(c.fHmpPesdPhmp),fHmpCkovPesd(c.fHmpCkovPesd),fHmpCkovPhmp(c.fHmpCkovPhmp),
65896ecb	166	fHmpMipCharge3cm(c.fHmpMipCharge3cm),fHmpMipTrkDist(c.fHmpMipTrkDist),
3882540a	167	fHmpTrkFlags(c.fHmpTrkFlags),
	168	fN1(c.fN1),
	169	fN2(c.fN2),
	170	fPionEff(c.fPionEff),
	171	fKaonEff(c.fKaonEff),
	172	fProtEff(c.fProtEff),
	173	fPionTot(c.fPionTot),
	174	fKaonTot(c.fKaonTot),
	175	fProtTot(c.fProtTot),
	176	fPionNot(c.fPionNot),
	177	fKaonNot(c.fKaonNot),
	178	fProtNot(c.fProtNot),
	179	fPionCon(c.fPionCon),
	180	fKaonCon(c.fKaonCon),
15939b4c	181	fProtCon(c.fProtCon)
3882540a	182	{
	183	//
	184	// Copy Constructor
	185	//
15939b4c	186	for (Int_t i=0; i<7; i++){
65896ecb	187	fHmpPhotons[i] = c.fHmpPhotons[i];
	188	fHmpPhotP[i] = c.fHmpPhotP[i];
	189	fHmpPhotSin2th[i] = c.fHmpPhotSin2th[i];
	190	fHmpMipTrkDistPosX[i] = c.fHmpMipTrkDistPosX[i];
	191	fHmpMipTrkDistNegX[i] = c.fHmpMipTrkDistNegX[i];
15939b4c	192	fHmpMipTrkDistPosY[i] = c.fHmpMipTrkDistPosY[i];
15939b4c	193	fHmpMipTrkDistNegY[i] = c.fHmpMipTrkDistNegY[i];
65896ecb	194	fHmpMipCharge[i] = c.fHmpMipCharge[i];
15939b4c	195	}
3882540a	196	}
	197
	198	//___________________________________________________________________________
	199	AliHMPIDTaskQA::~AliHMPIDTaskQA() {
	200	//
	201	//destructor
	202	//
	203	Info("~AliHMPIDTaskQA","Calling Destructor");
4860c9db	204	if (fHmpHistList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fHmpHistList;
3882540a	205	}
	206
	207	//___________________________________________________________________________
15939b4c	208	void AliHMPIDTaskQA::ConnectInputData(Option_t *option)
3882540a	209	{
15939b4c	210	AliAnalysisTaskSE::ConnectInputData(option);
3882540a	211
15939b4c	212	// Connect ESD here
3882540a	213	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
	214	if (!esdH) {
	215	AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
	216	} else
	217	fESD = esdH->GetEvent();
	218
	219	if (fUseMC){
	220	// Connect MC
	221	AliMCEventHandler mcH = dynamic_cast<AliMCEventHandler> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
	222	if (!mcH) {
	223	AliDebug(2,Form("ERROR: Could not get MCEventHandler"));
	224	fUseMC = kFALSE;
	225	} else
	226	fMC = mcH->MCEvent();
	227	if (!fMC) AliDebug(2,Form("ERROR: Could not get MCEvent"));
	228	}
	229	}
	230
	231	//___________________________________________________________________________
	232	void AliHMPIDTaskQA::UserExec(Option_t *)
	233	{
735bc477	234	Double_t priors[10]={1.,1.,1.,1.,1.,0.,0.,0.,0.,0.}; //{0.01,0.01,0.83,0.10,0.5};
3882540a	235	Double_t probs[5];
	236	AliPID *pPid = new AliPID();
	237	pPid->SetPriors(priors);
3882540a	238	AliESDtrack *track=0;
	239	TParticle *pPart=0;
	240	AliStack* pStack = 0;
001fa2ee	241	Int_t label = -1;
3882540a	242	if (fUseMC){
	243	pStack = fMC->Stack();
	244	}
	245
	246	//
	247	// Main loop function, executed on Event basis
	248	//
f5d9a192	249
	250	if( !fESD->GetPrimaryVertex()) return;
	251	if( fESD->GetPrimaryVertex()->GetNContributors() < 1 ) return;
	252	if( TMath::Abs(fESD->GetPrimaryVertex()->GetZ()) > 10.0 /* cm */) return;
	253
	254
	255	for (Int_t iTrack = 0; iTrack < fESD->GetNumberOfTracks(); iTrack++) {
3882540a	256
	257	track = fESD->GetTrack(iTrack);
	258	if(!track) continue;
f5d9a192	259	if(! track->IsOn(AliESDtrack::kITSrefit) \|\| !track->IsOn(AliESDtrack::kITSrefit) ) continue;
f5d9a192	260
3882540a	261	Double_t rin[3], rout[3];
	262	track->GetInnerXYZ(rin);
	263	track->GetOuterXYZ(rout);
	264	Double_t ktol = 0.001;
65896ecb	265	Double_t thetaCh = track->GetHMPIDsignal();
3882540a	266
65896ecb	267	if(Equal(thetaCh,-20.,ktol)) fHmpTrkFlags->Fill(0);
	268	else if(Equal(thetaCh,-9.,ktol)) fHmpTrkFlags->Fill(1);
	269	else if(Equal(thetaCh,-5.,ktol)) fHmpTrkFlags->Fill(2);
	270	else if(Equal(thetaCh,-11.,ktol)) fHmpTrkFlags->Fill(3);
	271	else if(Equal(thetaCh,-22.,ktol)) fHmpTrkFlags->Fill(4);
3882540a	272	else fHmpTrkFlags->Fill(4);
3882540a	273
65896ecb	274	if(Equal(thetaCh,-20.,ktol)) continue;
3882540a	275	if(track->GetHMPIDcluIdx() < 0) continue;
	276
	277	Int_t q, nph;
	278	Float_t x, y;
	279	Float_t xpc, ypc, th, ph;
	280	track->GetHMPIDmip(x,y,q,nph);
	281	track->GetHMPIDtrk(xpc,ypc,th,ph);
	282
	283	if(Equal(x,0.,ktol) && Equal(y,0.,ktol) && Equal(xpc,0.,ktol) && Equal(ypc,0.,ktol)) continue;
	284
15939b4c	285	Float_t sign = track->GetSign();
85c88881	286	Int_t ch = track->GetHMPIDcluIdx()/1000000;
	287	if( ch < 0 \|\| ch > 6 ) continue; //Chamber numbering is [0,6]
	288
65896ecb	289	if (sign > 0.) fHmpMipTrkDistPosX[ch]->Fill(xpc - x), fHmpMipTrkDistPosY[ch]->Fill(ypc - y);
	290	if (sign < 0.) fHmpMipTrkDistNegX[ch]->Fill(xpc - x), fHmpMipTrkDistNegY[ch]->Fill(ypc - y);
	291	fHmpMipCharge[ch]->Fill(q);
	292	Double_t dist = TMath::Sqrt((xpc - x)(xpc - x) + (ypc - y)(ypc - y));
3882540a	293	if (dist <= 3.0) fHmpMipCharge3cm->Fill(q);
65896ecb	294	fHmpMipTrkDist->Fill(dist);
3882540a	295
3882540a	296	track->GetHMPIDpid(probs);
	297	pPid->SetProbabilities(probs);
	298	if (fUseMC){
	299	if ((label = track->GetLabel()) < 0) continue;
	300	pPart = pStack->Particle(label);
	301	}
	302
65896ecb	303	if(thetaCh > 0 ){
	304	Double_t pHmp[3] = {0}, pHmp3 = 0;
	305	if (track->GetOuterHmpPxPyPz(pHmp)) pHmp3 = TMath::Sqrt(pHmp[0]pHmp[0]+pHmp[1]pHmp[1]+pHmp[2]*pHmp[2]);
15939b4c	306	if (pHmp3) fHmpPesdPhmp->Fill(track->P(),pHmp3);
65896ecb	307	fHmpCkovPesd->Fill(track->P(),thetaCh);
	308	if (pHmp3) fHmpCkovPhmp->Fill(pHmp3,thetaCh);
	309	fHmpPhotons[ch]->Fill(nph);
	310	fHmpPhotP[ch]->Fill(pHmp3,nph);
	311	Double_t sin2 = TMath::Power(TMath::Sin(th),2);
	312	fHmpPhotSin2th[ch]->Fill(sin2,nph);
3882540a	313
	314	if (fUseMC && dist<0.5 && TMath::Abs(th)<0.13){
	315	if (!pStack->IsPhysicalPrimary(label)) continue;
	316	if (track->Pt()<1. \|\| track->Pt()>5.) continue;
	317	Int_t pdgCode = TMath::Abs(pPart->GetPdgCode());
	318	Int_t ptBin=(Int_t) (2*(track->Pt()-1));
	319	if (pdgCode!=2212) fProtCon->Fill(ptBin);
	320	if (pdgCode==2212){
	321	fProtTot->Fill(ptBin);
	322	fProtEff->Fill(ptBin,pPid->GetProbability(AliPID::kProton));
	323	fPionNot->Fill(ptBin,pPid->GetProbability(AliPID::kPion));
	324	fKaonNot->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
	325	}
	326	if (pdgCode!=211) fPionCon->Fill(ptBin);
	327	if (pdgCode!=321) fKaonCon->Fill(ptBin);
	328	if (pdgCode==211){
	329	if (ptBin < 6){
	330	Float_t weight=pPid->GetProbability(AliPID::kElectron)+
	331	pPid->GetProbability(AliPID::kMuon)+
	332	pPid->GetProbability(AliPID::kPion);
	333	fPionTot->Fill(ptBin);
	334	fPionEff->Fill(ptBin,weight);
	335	fKaonNot->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
	336	}
	337	fProtNot->Fill(ptBin,pPid->GetProbability(AliPID::kProton));
	338	}
	339	if (pdgCode==321){
	340	if (ptBin < 6){
	341	fKaonTot->Fill(ptBin);
	342	fKaonEff->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
	343	fPionNot->Fill(ptBin,(pPid->GetProbability(AliPID::kPion)));
	344	}
	345	fProtNot->Fill(ptBin,(pPid->GetProbability(AliPID::kProton)));
	346	}
	347	}
	348	}//there is signal
3882540a	349	}//track loop
	350	delete pPid;
	351
	352	/* PostData(0) is taken care of by AliAnalysisTaskSE */
	353	PostData(1,fHmpHistList);
3882540a	354	}
	355
	356	//___________________________________________________________________________
	357	void AliHMPIDTaskQA::Terminate(Option_t*)
	358	{
	359	// The Terminate() function is the last function to be called during
	360	// a query. It always runs on the client, it can be used to present
	361	// the results graphically or save the results to file.
	362
	363	Info("Terminate"," ");
	364	AliAnalysisTaskSE::Terminate();
	365
	366	}
	367
	368	//___________________________________________________________________________
	369	void AliHMPIDTaskQA::UserCreateOutputObjects() {
	370	//
	371	//HERE ONE CAN CREATE OUTPUT OBJECTS
	372	//TO BE SET BEFORE THE EXECUTION OF THE TASK
	373	//
	374
	375	//slot #1
15939b4c	376	// OpenFile(1);
3882540a	377	fHmpHistList = new TList();
15939b4c	378	fHmpHistList->SetOwner();
	379
	380	fHmpPesdPhmp = new TH2F("fHmpPesdPhmp","HMPID: ESD p - running p;HMP p (GeV/c);ESD p (Gev/c)",100,0,10,100,0,10);
	381	fHmpHistList->Add(fHmpPesdPhmp);
	382
	383	fHmpCkovPesd = new TH2F("fHmpCkovPesd","HMPID: ThetaCherenkov vs P;p_esd (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
	384	fHmpHistList->Add(fHmpCkovPesd);
	385
	386	fHmpCkovPhmp = new TH2F("fHmpCkovPhmp","HMPID: ThetaCherenkov vs P;p_hmp (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
	387	fHmpHistList->Add(fHmpCkovPhmp);
3882540a	388
15939b4c	389	for (Int_t i=0; i<7; i++){
65896ecb	390	TString title=Form("MIP-Track distance in local X, Chamber %d;distance (cm);Entries",i);
	391	fHmpMipTrkDistPosX[i] = new TH1F(Form("fHmpMipTrkDistPosX%d",i),title.Data(),800,-20,20);
	392	fHmpHistList->Add(fHmpMipTrkDistPosX[i]);
	393	fHmpMipTrkDistNegX[i] = new TH1F(Form("fHmpMipTrkDistNegX%d",i),title.Data(),800,-20,20);
	394	fHmpHistList->Add(fHmpMipTrkDistNegX[i]);
	395
	396	title=Form("MIP-Track distance in local Y, Chamber %d;distance (cm);Entries",i);
15939b4c	397	fHmpMipTrkDistPosY[i] = new TH1F(Form("fHmpMipTrkDistPosY%d",i),title.Data(),800,-20,20);
	398	fHmpHistList->Add(fHmpMipTrkDistPosY[i]);
	399	fHmpMipTrkDistNegY[i] = new TH1F(Form("fHmpMipTrkDistNegY%d",i),title.Data(),800,-20,20);
	400	fHmpHistList->Add(fHmpMipTrkDistNegY[i]);
65896ecb	401
	402	title=Form("Mip charge distribution, Chamber %d;Charge;Entries",i);
	403	fHmpMipCharge[i] = new TH1F(Form("fHmpMipCharge%d",i),title.Data(),5001,-0.5,5000.5);
	404	fHmpHistList->Add(fHmpMipCharge[i]);
	405
	406	title=Form("Photons, Chamber %d;N of Photons;Entries",i);
	407	fHmpPhotons[i] = new TH1F(Form("fHmpPhotons%d",i),title.Data(),100,0,100);
	408	fHmpHistList->Add(fHmpPhotons[i]);
	409
	410	title=Form("Photons versus HMP momentum, Chamber %d;P (GeV);N of Photons",i);
	411	fHmpPhotP[i] = new TH2F(Form("fHmpPhotP%d",i),title.Data(),200,0.,20.,100,0,100);
	412	fHmpHistList->Add(fHmpPhotP[i]);
	413
	414	title=Form("Photons versus sin(th)^2 (uncorrected), Chamber %d;P (GeV);N of Photons",i);
	415	fHmpPhotSin2th[i] = new TH2F(Form("fHmpPhotSin2th%d",i),title.Data(),100,0.,1.,100,0,100);
	416	fHmpHistList->Add(fHmpPhotSin2th[i]);
15939b4c	417	}
3882540a	418
	419	fHmpMipCharge3cm = new TH1F("fHmpMipCharge3cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
	420	fHmpHistList->Add(fHmpMipCharge3cm);
	421
65896ecb	422	fHmpMipTrkDist = new TH1F("fHmpMipTrkDist","Mip-track distance in all the chambers",100,0.,10.);
	423	fHmpHistList->Add(fHmpMipTrkDist);
	424
3882540a	425	fHmpTrkFlags = new TH1F("fHmpTrkFlags","HMPID track flags",6,0,6);
	426	TString summary[6] = {"NotPerformed","MipDistCut", "MipQdcCut", "NoPhotAccept", "kNoRad", "other"};
	427	for(Int_t ibin = 0; ibin < 6; ibin++) fHmpTrkFlags->GetXaxis()->SetBinLabel(ibin+1,Form("%i %s",ibin+1,summary[ibin].Data()));
	428	fHmpHistList->Add(fHmpTrkFlags);
	429
	430	fPionEff = new TH1F("PionEff","Identified pions",fN1,0,fN1);
	431	fKaonEff = new TH1F("KaonEff","Identified kaons",fN1,0,fN1);
	432	fProtEff = new TH1F("ProtEff","Identified protons",fN2,0,fN2);
	433	fPionTot = new TH1I("PionTot","Total MC pions",fN1,0,fN1);
	434	fKaonTot = new TH1I("KaonTot","Total MC kaons",fN1,0,fN1);
	435	fProtTot = new TH1I("ProtTot","Total MC protons",fN2,0,fN2);
	436	fPionNot = new TH1F("PionNot","Misidentified pions",fN1,0,fN1);
	437	fKaonNot = new TH1F("KaonNot","Misidentified kaons",fN1,0,fN1);
	438	fProtNot = new TH1F("ProtNot","Misidentified protons",fN2,0,fN2);
	439	fPionCon = new TH1I("PionCon","Total not MC pions",fN1,0,fN1);
	440	fKaonCon = new TH1I("KaonCon","Total not MC kaons",fN1,0,fN1);
	441	fProtCon = new TH1I("ProtCon","Total not MC protons",fN2,0,fN2);
	442
	443	fHmpHistList->Add(fPionEff); fHmpHistList->Add(fKaonEff); fHmpHistList->Add(fProtEff);
	444	fHmpHistList->Add(fPionTot); fHmpHistList->Add(fKaonTot); fHmpHistList->Add(fProtTot);
	445	fHmpHistList->Add(fPionNot); fHmpHistList->Add(fKaonNot); fHmpHistList->Add(fProtNot);
	446	fHmpHistList->Add(fPionCon); fHmpHistList->Add(fKaonCon); fHmpHistList->Add(fProtCon);
	447
15939b4c	448	PostData(1,fHmpHistList);
3882540a	449	}
	450
	451	//____________________________________________________________________________________________________________________________________
	452	Bool_t AliHMPIDTaskQA::Equal(Double_t x, Double_t y, Double_t tolerance)
	453	{
	454	return abs(x - y) <= tolerance ;
	455	}
	456
	457	#endif