[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),
  fHmpMipTrkDistX(0x0),fHmpMipCharge3cm(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++)
    fHmpMipTrkDistPosY[i] = fHmpMipTrkDistNegY[i] = 0x0;
}

//___________________________________________________________________________
AliHMPIDTaskQA::AliHMPIDTaskQA(const Char_t* name) :
  AliAnalysisTaskSE(name),
  fESD(0x0), fMC(0x0), fUseMC(kTRUE),
  fHmpHistList(0x0),
  fHmpPesdPhmp(0x0),fHmpCkovPesd(0x0),fHmpCkovPhmp(0x0),
  fHmpMipTrkDistX(0x0),fHmpMipCharge3cm(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++)
    fHmpMipTrkDistPosY[i] = fHmpMipTrkDistNegY[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;
    fHmpMipTrkDistX  = c.fHmpMipTrkDistX;
    fHmpMipCharge3cm = c.fHmpMipCharge3cm;
    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++){
      fHmpMipTrkDistPosY[i] = c.fHmpMipTrkDistPosY[i];
      fHmpMipTrkDistNegY[i] = c.fHmpMipTrkDistNegY[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),
  fHmpMipTrkDistX(c.fHmpMipTrkDistX),
  fHmpMipCharge3cm(c.fHmpMipCharge3cm),
  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++){
    fHmpMipTrkDistPosY[i] = c.fHmpMipTrkDistPosY[i];
    fHmpMipTrkDistNegY[i] = c.fHmpMipTrkDistNegY[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[5]={1.,1.,1.,1.,1.}; //{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;
  if (fUseMC){
    pStack = fMC->Stack();
  }

  //
  // Main loop function, executed on Event basis
  //
  for (Int_t iTrack = 0; iTrack<fESD->GetNumberOfTracks(); iTrack++) {

    track = fESD->GetTrack(iTrack);
    if(!track) continue;
    Double_t rin[3], rout[3];
    track->GetInnerXYZ(rin);
    track->GetOuterXYZ(rout);
    Double_t ktol = 0.001;

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

    if(Equal(track->GetHMPIDsignal(),-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;
    Double_t dist = TMath::Sqrt( (xpc-x)*(xpc-x) + (ypc - y)*(ypc - y));
    fHmpMipTrkDistX->Fill(xpc - x);
    if (sign > 0.) fHmpMipTrkDistPosY[ch]->Fill(ypc - y);
    if (sign < 0.) fHmpMipTrkDistNegY[ch]->Fill(ypc - y);
    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 (dist <= 3.0) fHmpMipCharge3cm->Fill(q);

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

    if(track->GetHMPIDsignal() > 0 ){
      if (pHmp3) fHmpPesdPhmp->Fill(track->P(),pHmp3);
      fHmpCkovPesd->Fill(track->P(),track->GetHMPIDsignal());
      if (pHmp3) fHmpCkovPhmp->Fill(pHmp3,track->GetHMPIDsignal());

      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);

   fHmpMipTrkDistX = new TH1F("fHmpMipTrkDistX","HMPID MIP-Track distance in local X;distance (cm);Entries",800,-20,20);
   fHmpHistList->Add(fHmpMipTrkDistX);

   for (Int_t i=0; i<7; i++){
     TString 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]);
   }

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

   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),
15939b4c	46	fHmpMipTrkDistX(0x0),fHmpMipCharge3cm(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	//
15939b4c	66	for (Int_t i=0; i<7; i++)
15939b4c	67	fHmpMipTrkDistPosY[i] = fHmpMipTrkDistNegY[i] = 0x0;
3882540a	68	}
	69
	70	//___________________________________________________________________________
	71	AliHMPIDTaskQA::AliHMPIDTaskQA(const Char_t* name) :
	72	AliAnalysisTaskSE(name),
	73	fESD(0x0), fMC(0x0), fUseMC(kTRUE),
	74	fHmpHistList(0x0),
15939b4c	75	fHmpPesdPhmp(0x0),fHmpCkovPesd(0x0),fHmpCkovPhmp(0x0),
15939b4c	76	fHmpMipTrkDistX(0x0),fHmpMipCharge3cm(0x0),
3882540a	77	fHmpTrkFlags(0x0),
	78	fN1(6),
	79	fN2(8),
	80	fPionEff(0x0),
	81	fKaonEff(0x0),
	82	fProtEff(0x0),
	83	fPionTot(0x0),
	84	fKaonTot(0x0),
	85	fProtTot(0x0),
	86	fPionNot(0x0),
	87	fKaonNot(0x0),
	88	fProtNot(0x0),
	89	fPionCon(0x0),
	90	fKaonCon(0x0),
15939b4c	91	fProtCon(0x0)
3882540a	92	{
	93	//
	94	// Constructor. Initialization of Inputs and Outputs
	95	//
15939b4c	96	for (Int_t i=0; i<7; i++)
15939b4c	97	fHmpMipTrkDistPosY[i] = fHmpMipTrkDistNegY[i] = 0x0;
3882540a	98
3882540a	99	DefineOutput(1,TList::Class());
3882540a	100	}
	101
	102	//___________________________________________________________________________
	103	AliHMPIDTaskQA& AliHMPIDTaskQA::operator=(const AliHMPIDTaskQA& c)
	104	{
	105	//
	106	// Assignment operator
	107	//
	108	if (this!=&c) {
	109	AliAnalysisTaskSE::operator=(c);
	110	fESD = c.fESD;
	111	fMC = c.fMC;
	112	fUseMC = c.fUseMC;
	113	fHmpHistList = c.fHmpHistList;
15939b4c	114	fHmpPesdPhmp = c.fHmpPesdPhmp;
	115	fHmpCkovPesd = c.fHmpCkovPesd;
	116	fHmpCkovPhmp = c.fHmpCkovPhmp;
3882540a	117	fHmpMipTrkDistX = c.fHmpMipTrkDistX;
3882540a	118	fHmpMipCharge3cm = c.fHmpMipCharge3cm;
	119	fHmpTrkFlags = c.fHmpTrkFlags;
	120	fN1 = c.fN1;
	121	fN2 = c.fN2;
	122	fPionEff = c.fPionEff;
	123	fKaonEff = c.fKaonEff;
	124	fProtEff = c.fProtEff;
	125	fPionTot = c.fPionTot;
	126	fKaonTot = c.fKaonTot;
	127	fProtTot = c.fProtTot;
	128	fPionNot = c.fPionNot;
	129	fKaonNot = c.fKaonNot;
	130	fProtNot = c.fProtNot;
	131	fPionCon = c.fPionCon;
	132	fKaonCon = c.fKaonCon;
	133	fProtCon = c.fProtCon;
15939b4c	134	for (Int_t i=0; i<7; i++){
	135	fHmpMipTrkDistPosY[i] = c.fHmpMipTrkDistPosY[i];
	136	fHmpMipTrkDistNegY[i] = c.fHmpMipTrkDistNegY[i];
	137	}
3882540a	138	}
	139	return *this;
	140	}
	141
	142	//___________________________________________________________________________
	143	AliHMPIDTaskQA::AliHMPIDTaskQA(const AliHMPIDTaskQA& c) :
	144	AliAnalysisTaskSE(c),
	145	fESD(c.fESD),fMC(c.fMC),fUseMC(c.fUseMC),
	146	fHmpHistList(c.fHmpHistList),
15939b4c	147	fHmpPesdPhmp(c.fHmpPesdPhmp),fHmpCkovPesd(c.fHmpCkovPesd),fHmpCkovPhmp(c.fHmpCkovPhmp),
15939b4c	148	fHmpMipTrkDistX(c.fHmpMipTrkDistX),
3882540a	149	fHmpMipCharge3cm(c.fHmpMipCharge3cm),
	150	fHmpTrkFlags(c.fHmpTrkFlags),
	151	fN1(c.fN1),
	152	fN2(c.fN2),
	153	fPionEff(c.fPionEff),
	154	fKaonEff(c.fKaonEff),
	155	fProtEff(c.fProtEff),
	156	fPionTot(c.fPionTot),
	157	fKaonTot(c.fKaonTot),
	158	fProtTot(c.fProtTot),
	159	fPionNot(c.fPionNot),
	160	fKaonNot(c.fKaonNot),
	161	fProtNot(c.fProtNot),
	162	fPionCon(c.fPionCon),
	163	fKaonCon(c.fKaonCon),
15939b4c	164	fProtCon(c.fProtCon)
3882540a	165	{
	166	//
	167	// Copy Constructor
	168	//
15939b4c	169	for (Int_t i=0; i<7; i++){
	170	fHmpMipTrkDistPosY[i] = c.fHmpMipTrkDistPosY[i];
	171	fHmpMipTrkDistNegY[i] = c.fHmpMipTrkDistNegY[i];
	172	}
3882540a	173	}
	174
	175	//___________________________________________________________________________
	176	AliHMPIDTaskQA::~AliHMPIDTaskQA() {
	177	//
	178	//destructor
	179	//
	180	Info("~AliHMPIDTaskQA","Calling Destructor");
	181	if (fHmpHistList /&& !AliAnalysisManager::GetAnalysisManager()->IsProofMode()/) delete fHmpHistList;
	182	}
	183
	184	//___________________________________________________________________________
15939b4c	185	void AliHMPIDTaskQA::ConnectInputData(Option_t *option)
3882540a	186	{
15939b4c	187	AliAnalysisTaskSE::ConnectInputData(option);
3882540a	188
15939b4c	189	// Connect ESD here
3882540a	190	AliESDInputHandler esdH = dynamic_cast<AliESDInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
	191	if (!esdH) {
	192	AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
	193	} else
	194	fESD = esdH->GetEvent();
	195
	196	if (fUseMC){
	197	// Connect MC
	198	AliMCEventHandler mcH = dynamic_cast<AliMCEventHandler> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
	199	if (!mcH) {
	200	AliDebug(2,Form("ERROR: Could not get MCEventHandler"));
	201	fUseMC = kFALSE;
	202	} else
	203	fMC = mcH->MCEvent();
	204	if (!fMC) AliDebug(2,Form("ERROR: Could not get MCEvent"));
	205	}
	206	}
	207
	208	//___________________________________________________________________________
	209	void AliHMPIDTaskQA::UserExec(Option_t *)
	210	{
	211	Double_t priors[5]={1.,1.,1.,1.,1.}; //{0.01,0.01,0.83,0.10,0.5};
	212	Double_t probs[5];
	213	AliPID *pPid = new AliPID();
	214	pPid->SetPriors(priors);
3882540a	215	AliESDtrack *track=0;
	216	TParticle *pPart=0;
	217	AliStack* pStack = 0;
	218	Int_t label;
	219	if (fUseMC){
	220	pStack = fMC->Stack();
	221	}
	222
	223	//
	224	// Main loop function, executed on Event basis
	225	//
	226	for (Int_t iTrack = 0; iTrack<fESD->GetNumberOfTracks(); iTrack++) {
	227
	228	track = fESD->GetTrack(iTrack);
	229	if(!track) continue;
3882540a	230	Double_t rin[3], rout[3];
	231	track->GetInnerXYZ(rin);
	232	track->GetOuterXYZ(rout);
	233	Double_t ktol = 0.001;
	234
	235	if(Equal(track->GetHMPIDsignal(),-20.,ktol)) fHmpTrkFlags->Fill(0);
	236	else if(Equal(track->GetHMPIDsignal(),-9.,ktol)) fHmpTrkFlags->Fill(1);
	237	else if(Equal(track->GetHMPIDsignal(),-5.,ktol)) fHmpTrkFlags->Fill(2);
	238	else if(Equal(track->GetHMPIDsignal(),-11.,ktol)) fHmpTrkFlags->Fill(3);
	239	else if(Equal(track->GetHMPIDsignal(),-22.,ktol)) fHmpTrkFlags->Fill(4);
	240	else fHmpTrkFlags->Fill(4);
	241
	242	if(Equal(track->GetHMPIDsignal(),-20.,ktol)) continue;
	243	if(track->GetHMPIDcluIdx() < 0) continue;
	244
	245	Int_t q, nph;
	246	Float_t x, y;
	247	Float_t xpc, ypc, th, ph;
	248	track->GetHMPIDmip(x,y,q,nph);
	249	track->GetHMPIDtrk(xpc,ypc,th,ph);
	250
	251	if(Equal(x,0.,ktol) && Equal(y,0.,ktol) && Equal(xpc,0.,ktol) && Equal(ypc,0.,ktol)) continue;
	252
15939b4c	253	Float_t sign = track->GetSign();
15939b4c	254	Int_t ch = track->GetHMPIDcluIdx()/1000000;
3882540a	255	Double_t dist = TMath::Sqrt( (xpc-x)(xpc-x) + (ypc - y)(ypc - y));
15939b4c	256	fHmpMipTrkDistX->Fill(xpc - x);
	257	if (sign > 0.) fHmpMipTrkDistPosY[ch]->Fill(ypc - y);
	258	if (sign < 0.) fHmpMipTrkDistNegY[ch]->Fill(ypc - y);
3882540a	259	Double_t pHmp[3] = {0}, pHmp3 = 0;
	260	if (track->GetOuterHmpPxPyPz(pHmp)) pHmp3 = TMath::Sqrt(pHmp[0]pHmp[0]+pHmp[1]pHmp[1]+pHmp[2]*pHmp[2]);
	261	if (dist <= 3.0) fHmpMipCharge3cm->Fill(q);
	262
3882540a	263	track->GetHMPIDpid(probs);
	264	pPid->SetProbabilities(probs);
	265	if (fUseMC){
	266	if ((label = track->GetLabel()) < 0) continue;
	267	pPart = pStack->Particle(label);
	268	}
	269
	270	if(track->GetHMPIDsignal() > 0 ){
15939b4c	271	if (pHmp3) fHmpPesdPhmp->Fill(track->P(),pHmp3);
	272	fHmpCkovPesd->Fill(track->P(),track->GetHMPIDsignal());
	273	if (pHmp3) fHmpCkovPhmp->Fill(pHmp3,track->GetHMPIDsignal());
3882540a	274
	275	if (fUseMC && dist<0.5 && TMath::Abs(th)<0.13){
	276	if (!pStack->IsPhysicalPrimary(label)) continue;
	277	if (track->Pt()<1. \|\| track->Pt()>5.) continue;
	278	Int_t pdgCode = TMath::Abs(pPart->GetPdgCode());
	279	Int_t ptBin=(Int_t) (2*(track->Pt()-1));
	280	if (pdgCode!=2212) fProtCon->Fill(ptBin);
	281	if (pdgCode==2212){
	282	fProtTot->Fill(ptBin);
	283	fProtEff->Fill(ptBin,pPid->GetProbability(AliPID::kProton));
	284	fPionNot->Fill(ptBin,pPid->GetProbability(AliPID::kPion));
	285	fKaonNot->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
	286	}
	287	if (pdgCode!=211) fPionCon->Fill(ptBin);
	288	if (pdgCode!=321) fKaonCon->Fill(ptBin);
	289	if (pdgCode==211){
	290	if (ptBin < 6){
	291	Float_t weight=pPid->GetProbability(AliPID::kElectron)+
	292	pPid->GetProbability(AliPID::kMuon)+
	293	pPid->GetProbability(AliPID::kPion);
	294	fPionTot->Fill(ptBin);
	295	fPionEff->Fill(ptBin,weight);
	296	fKaonNot->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
	297	}
	298	fProtNot->Fill(ptBin,pPid->GetProbability(AliPID::kProton));
	299	}
	300	if (pdgCode==321){
	301	if (ptBin < 6){
	302	fKaonTot->Fill(ptBin);
	303	fKaonEff->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
	304	fPionNot->Fill(ptBin,(pPid->GetProbability(AliPID::kPion)));
	305	}
	306	fProtNot->Fill(ptBin,(pPid->GetProbability(AliPID::kProton)));
	307	}
	308	}
	309	}//there is signal
3882540a	310	}//track loop
	311	delete pPid;
	312
	313	/* PostData(0) is taken care of by AliAnalysisTaskSE */
	314	PostData(1,fHmpHistList);
3882540a	315	}
	316
	317	//___________________________________________________________________________
	318	void AliHMPIDTaskQA::Terminate(Option_t*)
	319	{
	320	// The Terminate() function is the last function to be called during
	321	// a query. It always runs on the client, it can be used to present
	322	// the results graphically or save the results to file.
	323
	324	Info("Terminate"," ");
	325	AliAnalysisTaskSE::Terminate();
	326
	327	}
	328
	329	//___________________________________________________________________________
	330	void AliHMPIDTaskQA::UserCreateOutputObjects() {
	331	//
	332	//HERE ONE CAN CREATE OUTPUT OBJECTS
	333	//TO BE SET BEFORE THE EXECUTION OF THE TASK
	334	//
	335
	336	//slot #1
15939b4c	337	// OpenFile(1);
3882540a	338	fHmpHistList = new TList();
15939b4c	339	fHmpHistList->SetOwner();
	340
	341	fHmpPesdPhmp = new TH2F("fHmpPesdPhmp","HMPID: ESD p - running p;HMP p (GeV/c);ESD p (Gev/c)",100,0,10,100,0,10);
	342	fHmpHistList->Add(fHmpPesdPhmp);
	343
	344	fHmpCkovPesd = new TH2F("fHmpCkovPesd","HMPID: ThetaCherenkov vs P;p_esd (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
	345	fHmpHistList->Add(fHmpCkovPesd);
	346
	347	fHmpCkovPhmp = new TH2F("fHmpCkovPhmp","HMPID: ThetaCherenkov vs P;p_hmp (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
	348	fHmpHistList->Add(fHmpCkovPhmp);
3882540a	349
	350	fHmpMipTrkDistX = new TH1F("fHmpMipTrkDistX","HMPID MIP-Track distance in local X;distance (cm);Entries",800,-20,20);
	351	fHmpHistList->Add(fHmpMipTrkDistX);
	352
15939b4c	353	for (Int_t i=0; i<7; i++){
	354	TString title=Form("MIP-Track distance in local Y, Chamber %d;distance (cm);Entries",i);
	355	fHmpMipTrkDistPosY[i] = new TH1F(Form("fHmpMipTrkDistPosY%d",i),title.Data(),800,-20,20);
	356	fHmpHistList->Add(fHmpMipTrkDistPosY[i]);
	357	fHmpMipTrkDistNegY[i] = new TH1F(Form("fHmpMipTrkDistNegY%d",i),title.Data(),800,-20,20);
	358	fHmpHistList->Add(fHmpMipTrkDistNegY[i]);
	359	}
3882540a	360
	361	fHmpMipCharge3cm = new TH1F("fHmpMipCharge3cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
	362	fHmpHistList->Add(fHmpMipCharge3cm);
	363
	364	fHmpTrkFlags = new TH1F("fHmpTrkFlags","HMPID track flags",6,0,6);
	365	TString summary[6] = {"NotPerformed","MipDistCut", "MipQdcCut", "NoPhotAccept", "kNoRad", "other"};
	366	for(Int_t ibin = 0; ibin < 6; ibin++) fHmpTrkFlags->GetXaxis()->SetBinLabel(ibin+1,Form("%i %s",ibin+1,summary[ibin].Data()));
	367	fHmpHistList->Add(fHmpTrkFlags);
	368
	369	fPionEff = new TH1F("PionEff","Identified pions",fN1,0,fN1);
	370	fKaonEff = new TH1F("KaonEff","Identified kaons",fN1,0,fN1);
	371	fProtEff = new TH1F("ProtEff","Identified protons",fN2,0,fN2);
	372	fPionTot = new TH1I("PionTot","Total MC pions",fN1,0,fN1);
	373	fKaonTot = new TH1I("KaonTot","Total MC kaons",fN1,0,fN1);
	374	fProtTot = new TH1I("ProtTot","Total MC protons",fN2,0,fN2);
	375	fPionNot = new TH1F("PionNot","Misidentified pions",fN1,0,fN1);
	376	fKaonNot = new TH1F("KaonNot","Misidentified kaons",fN1,0,fN1);
	377	fProtNot = new TH1F("ProtNot","Misidentified protons",fN2,0,fN2);
	378	fPionCon = new TH1I("PionCon","Total not MC pions",fN1,0,fN1);
	379	fKaonCon = new TH1I("KaonCon","Total not MC kaons",fN1,0,fN1);
	380	fProtCon = new TH1I("ProtCon","Total not MC protons",fN2,0,fN2);
	381
	382	fHmpHistList->Add(fPionEff); fHmpHistList->Add(fKaonEff); fHmpHistList->Add(fProtEff);
	383	fHmpHistList->Add(fPionTot); fHmpHistList->Add(fKaonTot); fHmpHistList->Add(fProtTot);
	384	fHmpHistList->Add(fPionNot); fHmpHistList->Add(fKaonNot); fHmpHistList->Add(fProtNot);
	385	fHmpHistList->Add(fPionCon); fHmpHistList->Add(fKaonCon); fHmpHistList->Add(fProtCon);
	386
15939b4c	387	PostData(1,fHmpHistList);
3882540a	388	}
	389
	390	//____________________________________________________________________________________________________________________________________
	391	Bool_t AliHMPIDTaskQA::Equal(Double_t x, Double_t y, Double_t tolerance)
	392	{
	393	return abs(x - y) <= tolerance ;
	394	}
	395
	396	#endif