--- /dev/null
+# -*- mode: CMake -*-
+#-----------------------------------------------------------------------#
+# Package File for PWG2spectra #
+# Author : Johny Jose (johny.jose@cern.ch) #
+# Variables Defined : #
+# #
+# SRCS - C++ source files #
+# HDRS - C++ header files #
+# DHDR - ROOT Dictionary Linkdef header file #
+# CSRCS - C source files #
+# CHDRS - C header files #
+# EINCLUDE - Include directories #
+# EDEFINE - Compiler definitions #
+# ELIBS - Extra libraries to link #
+# ELIBSDIR - Extra library directories #
+# PACKFFLAGS - Fortran compiler flags for package #
+# PACKCXXFLAGS - C++ compiler flags for package #
+# PACKCFLAGS - C compiler flags for package #
+# PACKSOFLAGS - Shared library linking flags #
+# PACKLDFLAGS - Module linker flags #
+# PACKBLIBS - Libraries to link (Executables only) #
+# EXPORT - Header files to be exported #
+# CINTHDRS - Dictionary header files #
+# CINTAUTOLINK - Set automatic dictionary generation #
+# ARLIBS - Archive Libraries and objects for linking (Executables only) #
+# SHLIBS - Shared Libraries and objects for linking (Executables only) #
+#-----------------------------------------------------------------------#
+
+set ( SRCS
+ QATasks/AliAnalysisTaskQAHighPtDeDx.cxx
+ QATasks/AliQAProdMultistrange.cxx
+)
+
+string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
+
+set ( DHDR PWGLFQATasksLinkDef.h)
+
+set ( EXPORT )
+
+set ( EINCLUDE ANALYSIS PWGLF/QATasks STEER/ESD STEER/STEERBase)
--- /dev/null
+#ifdef __CINT__
+
+#pragma link off all glols;
+#pragma link off all classes;
+#pragma link off all functions;
+
+#pragma link C++ class AliAnalysisTaskQAHighPtDeDx+;
+#pragma link C++ class AliQAProdMultistrange+;
+
+#endif
--- /dev/null
+ /*************************************************************************
+* Copyright(c) 1998-2008, 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. * **************************************************************************/
+
+
+#include "AliAnalysisTaskQAHighPtDeDx.h"
+
+// ROOT includes
+#include <TList.h>
+#include <TTree.h>
+#include <TMath.h>
+#include <TH1.h>
+#include <TF1.h>
+#include <TProfile.h>
+#include <TParticle.h>
+#include <TFile.h>
+
+// AliRoot includes
+#include <AliAnalysisManager.h>
+#include <AliAnalysisFilter.h>
+#include <AliESDInputHandler.h>
+#include <AliESDEvent.h>
+#include <AliESDVertex.h>
+#include <AliLog.h>
+#include <AliExternalTrackParam.h>
+#include <AliESDtrackCuts.h>
+#include <AliESDVZERO.h>
+#include <AliAODVZERO.h>
+
+#include <AliMCEventHandler.h>
+#include <AliMCEvent.h>
+#include <AliStack.h>
+
+#include <TTreeStream.h>
+
+#include <AliHeader.h>
+#include <AliGenPythiaEventHeader.h>
+#include <AliGenDPMjetEventHeader.h>
+
+#include "AliCentrality.h"
+#include <AliESDv0.h>
+#include <AliKFVertex.h>
+#include <AliAODVertex.h>
+
+#include <AliAODTrack.h>
+#include <AliAODPid.h>
+#include <AliAODMCHeader.h>
+
+
+// STL includes
+#include <iostream>
+using namespace std;
+
+
+//
+// Responsible:
+// Antonio Ortiz (Lund)
+// Peter Christiansen (Lund)
+//
+
+
+
+
+
+
+
+
+const Double_t AliAnalysisTaskQAHighPtDeDx::fgkClight = 2.99792458e-2;
+Float_t magf = -1;
+TF1* cutLow = new TF1("StandardPhiCutLow", "0.1/x/x+pi/18.0-0.025", 0, 50);
+TF1* cutHigh = new TF1("StandardPhiCutHigh", "0.12/x+pi/18.0+0.035", 0, 50);
+Double_t DeDxMIPMin = 30;
+Double_t DeDxMIPMax = 65;
+const Int_t nHists = 9;
+
+Int_t etaLow[nHists] = {-8, -8, -6, -4, -2, 0, 2, 4, 6};
+Int_t etaHigh[nHists] = { 8, -6, -4, -2, 0, 2, 4, 6, 8};
+
+Int_t nDeltaPiBins = 60;
+Double_t deltaPiLow = 40;
+Double_t deltaPiHigh = 100;
+const Char_t *Pid[7]={"Ch","Pion","Kaon","Proton","Electron","Muon","Oher"};
+ClassImp(AliAnalysisTaskQAHighPtDeDx)
+//_____________________________________________________________________________
+AliAnalysisTaskQAHighPtDeDx::AliAnalysisTaskQAHighPtDeDx(const char *name):
+ AliAnalysisTaskSE(name),
+ fESD(0x0),
+ fAOD(0x0),
+ fMC(0x0),
+ fMCStack(0x0),
+ fMCArray(0x0),
+ fTrackFilterGolden(0x0),
+ fTrackFilterTPC(0x0),
+ fCentEst("V0M"),
+ fAnalysisType("ESD"),
+ fAnalysisMC(kFALSE),
+ fAnalysisPbPb(kFALSE),
+ ftrigBit(0x0),
+ fRandom(0x0),
+ fPileUpRej(kFALSE),
+ fVtxCut(10.0),
+ fEtaCut(0.9),
+ fMinCent(0.0),
+ fMaxCent(100.0),
+ fStoreMcIn(kFALSE),//
+ fMcProcessType(-999),
+ fTriggeredEventMB(-999),
+ fVtxStatus(-999),
+ fZvtx(-999),
+ fZvtxMC(-999),
+ fRun(-999),
+ fEventId(-999),
+ fListOfObjects(0),
+ fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
+ fn1(0x0),
+ fcent(0x0),
+ hMIPVsEta(0x0),
+ pMIPVsEta(0x0),
+ hMIPVsEtaV0s(0x0),
+ pMIPVsEtaV0s(0x0),
+ hPlateauVsEta(0x0),
+ pPlateauVsEta(0x0),
+ hPhi(0x0)
+
+
+{
+ // Default constructor (should not be used)
+ for(Int_t i=0;i<9;++i){
+
+ hMIPVsNch[i]=0;//TH2D, MIP vs Nch for different eta intervals
+ pMIPVsNch[i]=0;//TProfile, MIP vs Nch for different eta intervals
+ hMIPVsPhi[i]=0;//TH2D, MIP vs phi for different eta intervals
+ pMIPVsPhi[i]=0;//TProfile, MIP vs phi for different eta intervals
+ hPlateauVsPhi[i]=0;//TH2D, dE/dx vs Phi, electrons 0.4<p<0.6 GeV/c
+ pPlateauVsPhi[i]=0;//TProfile, dE/dx vs Phi, electrons 0.4<p<0.6 GeV/c
+ histPiV0[i]=0;//TH2D, dE/dx vs p, pi id by V0s
+ histPV0[i]=0;// TH2D, dE/dx vs p, p id by V0s
+ histAllCh[i]=0;//TH2D, dE/dx vs p for all charged particles
+ histPiTof[i]=0;//TH2D, dE/dx vs p for a "clean" sample of pions, beta>1
+ histEV0[i]=0;
+
+ for(Int_t pid=0;pid<7;++pid){
+ hMcIn[pid][i]=0;
+ hMcOut[pid][i]=0;
+ }
+
+ }
+ DefineOutput(1, TList::Class());//esto es nuevo
+}
+
+//______________________________________________________________________________
+void AliAnalysisTaskQAHighPtDeDx::UserCreateOutputObjects()
+{
+ // This method is called once per worker node
+ // Here we define the output: histograms and debug tree if requested
+ // We also create the random generator here so it might get different seeds...
+ fRandom = new TRandom(0); // 0 means random seed
+
+
+
+ //OpenFile(1);
+ fListOfObjects = new TList();
+ fListOfObjects->SetOwner();
+
+
+
+
+ //
+ // Histograms
+ //
+ fEvents = new TH1I("fEvents","Number of analyzed events; Events; Counts", 3, 0, 3);
+ fListOfObjects->Add(fEvents);
+
+ fn1=new TH1F("fn1","fn1",11,-1,10);
+ fListOfObjects->Add(fn1);
+
+ fcent=new TH1F("fcent","fcent",104,-2,102);
+ fListOfObjects->Add(fcent);
+
+ fVtx = new TH1I("fVtx","Vtx info (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
+ fListOfObjects->Add(fVtx);
+
+ fVtxBeforeCuts = new TH1F("fVtxBeforeCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
+ fListOfObjects->Add(fVtxBeforeCuts);
+
+ fVtxAfterCuts = new TH1F("fVtxAfterCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
+ fListOfObjects->Add(fVtxAfterCuts);
+
+
+ const Int_t nPtBinsV0s = 25;
+ Double_t ptBinsV0s[nPtBinsV0s+1] = { 0.0 , 0.1 , 0.2 , 0.3 , 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 1.0 ,
+ 1.2 , 1.4 , 1.6 , 1.8 , 2.0 , 2.5 , 3.0 , 3.5 , 4.0 , 5.0 , 7.0 ,
+ 9.0 , 12.0, 15.0, 20.0 };
+
+
+
+
+ const Char_t* ending[nHists] = {"", "86", "64", "42", "20", "02", "24", "46", "68"};
+
+ const Char_t* LatexEta[nHists] = {
+ "|#eta|<0.8", "-0.8<#eta<-0.6", "-0.6<#eta<-0.4", "-0.4<#eta<-0.2", "-0.2<#eta<0",
+ "0<#eta<0.2", "0.2<#eta<0.4", "0.4<#eta<0.6", "0.6<#eta<0.8"
+ };
+ hPhi = new TH2D("histPhi", "pt; #phi'", nPtBinsV0s, ptBinsV0s, 90, -0.05, 0.4);
+ //dE/dx vs phi, pions at the MIP
+ fListOfObjects->Add(hPhi);
+
+
+
+
+ Int_t nPhiBins = 36;
+
+ for(Int_t i = 0; i < nHists; i++) {
+
+
+ hMIPVsPhi[i] = new TH2D(Form("hMIPVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx MIP",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(),
+ DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
+ hMIPVsPhi[i]->Sumw2();
+
+ pMIPVsPhi[i] = new TProfile(Form("pMIPVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx MIP",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(),
+ DeDxMIPMin, DeDxMIPMax);
+ pMIPVsPhi[i]->SetMarkerStyle(20);
+ pMIPVsPhi[i]->SetMarkerColor(1);
+ pMIPVsPhi[i]->SetLineColor(1);
+ pMIPVsPhi[i]->Sumw2();
+
+ fListOfObjects->Add(hMIPVsPhi[i]);
+ fListOfObjects->Add(pMIPVsPhi[i]);
+
+ hPlateauVsPhi[i] = new TH2D(Form("hPlateauVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx Plateau",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(),
+ 95-DeDxMIPMax, DeDxMIPMax, 95);
+ hPlateauVsPhi[i]->Sumw2();
+
+ pPlateauVsPhi[i] = new TProfile(Form("pPlateauVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx Plateau",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(),
+ DeDxMIPMax, 95);
+ pPlateauVsPhi[i]->SetMarkerStyle(20);
+ pPlateauVsPhi[i]->SetMarkerColor(1);
+ pPlateauVsPhi[i]->SetLineColor(1);
+ pPlateauVsPhi[i]->Sumw2();
+
+ fListOfObjects->Add(hPlateauVsPhi[i]);
+ fListOfObjects->Add(pPlateauVsPhi[i]);
+
+
+ hMIPVsNch[i] = new TH2D(Form("hMIPVsNch%s", ending[i]), Form("%s; TPC track mult. |#eta|<0.8; dE/dx MIP",LatexEta[i]), 400, 1, 2001,
+ DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
+ hMIPVsNch[i]->Sumw2();
+
+ pMIPVsNch[i] = new TProfile(Form("pMIPVsNch%s", ending[i]), Form("%s; TPC track mult. |#eta|<0.8; dE/dx MIP",LatexEta[i]), 400, 1, 2001, DeDxMIPMin, DeDxMIPMax);
+ pMIPVsNch[i]->SetMarkerStyle(20);
+ pMIPVsNch[i]->SetMarkerColor(1);
+ pMIPVsNch[i]->SetLineColor(1);
+ pMIPVsNch[i]->Sumw2();
+
+ fListOfObjects->Add(hMIPVsNch[i]);
+ fListOfObjects->Add(pMIPVsNch[i]);
+
+
+ histPiV0[i] = new TH2D(Form("histPiV0%s", ending[i]), "Pions id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
+ histPiV0[i]->Sumw2();
+ histPV0[i] = new TH2D(Form("histPV0%s", ending[i]), "Protons id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
+
+ fListOfObjects->Add(histPiV0[i]);
+ fListOfObjects->Add(histPV0[i]);
+
+ histPiTof[i] = new TH2D(Form("histPiTof%s", ending[i]), "all charged", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
+ histPiTof[i]->Sumw2();
+
+ histAllCh[i] = new TH2D(Form("histAllCh%s", ending[i]), "Pions id by TOF", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
+ histAllCh[i]->Sumw2();
+
+ fListOfObjects->Add(histPiTof[i]);
+ fListOfObjects->Add(histAllCh[i]);
+
+ histEV0[i] = new TH2D(Form("histEV0%s", ending[i]), "Electrons id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
+ histEV0[i]->Sumw2();
+ fListOfObjects->Add(histEV0[i]);
+
+
+
+ }
+
+
+ hMIPVsEta = new TH2D("hMIPVsEta","; #eta; dE/dx_{MIP, primary tracks}",16,-0.8,0.8,DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
+ pMIPVsEta = new TProfile("pMIPVsEta","; #eta; #LT dE/dx #GT_{MIP, primary tracks}",16,-0.8,0.8, DeDxMIPMin, DeDxMIPMax);
+ hMIPVsEtaV0s = new TH2D("hMIPVsEtaV0s","; #eta; dE/dx_{MIP, secondary tracks}",16,-0.8,0.8,DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
+ pMIPVsEtaV0s = new TProfile("pMIPVsEtaV0s","; #eta; #LT dE/dx #GT_{MIP, secondary tracks}",16,-0.8,0.8,DeDxMIPMin, DeDxMIPMax);
+
+ hPlateauVsEta = new TH2D("hPlateauVsEta","; #eta; dE/dx_{Plateau, primary tracks}",16,-0.8,0.8,95-DeDxMIPMax, DeDxMIPMax, 95);
+ pPlateauVsEta = new TProfile("pPlateauVsEta","; #eta; #LT dE/dx #GT_{Plateau, primary tracks}",16,-0.8,0.8, DeDxMIPMax, 95);
+
+ fListOfObjects->Add(hMIPVsEta);
+ fListOfObjects->Add(pMIPVsEta);
+ fListOfObjects->Add(hMIPVsEtaV0s);
+ fListOfObjects->Add(pMIPVsEtaV0s);
+ fListOfObjects->Add(hPlateauVsEta);
+ fListOfObjects->Add(pPlateauVsEta);
+
+
+
+
+
+
+ if (fAnalysisMC) {
+ for(Int_t i = 0; i < nHists; i++) {
+ for(Int_t pid = 0; pid < 7; pid++) {
+
+ hMcIn[pid][i] = new TH1D(Form("hIn%s%s", Pid[pid],ending[i]), Form("MC in (pid %s)", Pid[pid]),
+ nPtBinsV0s, ptBinsV0s);
+ fListOfObjects->Add(hMcIn[pid][i]);
+
+ hMcOut[pid][i] = new TH1D(Form("hMcOut%s%s", Pid[pid],ending[i]), Form("MC out (pid %s)", Pid[pid]),
+ nPtBinsV0s, ptBinsV0s);
+ fListOfObjects->Add(hMcOut[pid][i]);
+
+
+ }
+ }
+
+ fVtxMC = new TH1F("fVtxMC","mc vtx", 120, -30, 30);
+ fListOfObjects->Add(fVtxMC);
+
+ }
+
+ // Post output data.
+ PostData(1, fListOfObjects);
+}
+
+//______________________________________________________________________________
+void AliAnalysisTaskQAHighPtDeDx::UserExec(Option_t *)
+{
+ // Main loop
+
+ //
+ // First we make sure that we have valid input(s)!
+ //
+
+
+
+ AliVEvent *event = InputEvent();
+ if (!event) {
+ Error("UserExec", "Could not retrieve event");
+ return;
+ }
+
+
+
+ if (fAnalysisType == "ESD"){
+ fESD = dynamic_cast<AliESDEvent*>(event);
+ if(!fESD){
+ Printf("%s:%d ESDEvent not found in Input Manager",(char*)__FILE__,__LINE__);
+ this->Dump();
+ return;
+ }
+ } else {
+ fAOD = dynamic_cast<AliAODEvent*>(event);
+ if(!fAOD){
+ Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
+ this->Dump();
+ return;
+ }
+ }
+
+
+
+ if (fAnalysisMC) {
+
+ if (fAnalysisType == "ESD"){
+ fMC = dynamic_cast<AliMCEvent*>(MCEvent());
+ if(!fMC){
+ Printf("%s:%d MCEvent not found in Input Manager",(char*)__FILE__,__LINE__);
+ this->Dump();
+ return;
+ }
+
+ fMCStack = fMC->Stack();
+
+ if(!fMCStack){
+ Printf("%s:%d MCStack not found in Input Manager",(char*)__FILE__,__LINE__);
+ this->Dump();
+ return;
+ }
+ } else { // AOD
+
+ fMC = dynamic_cast<AliMCEvent*>(MCEvent());
+ if(fMC)
+ fMC->Dump();
+
+ fMCArray = (TClonesArray*)fAOD->FindListObject("mcparticles");
+ if(!fMCArray){
+ Printf("%s:%d AOD MC array not found in Input Manager",(char*)__FILE__,__LINE__);
+ this->Dump();
+ return;
+ }
+ }
+ }
+
+
+ // Get trigger decision
+ fTriggeredEventMB = 0; //init
+
+ fn1->Fill(0);
+
+ if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
+ ->IsEventSelected() & ftrigBit ){
+ fTriggeredEventMB = 1; //event triggered as minimum bias
+ }
+
+ // Get process type for MC
+ fMcProcessType = 0; // -1=invalid, 0=data, 1=ND, 2=SD, 3=DD
+
+ // real data that are not triggered we skip
+ if(!fAnalysisMC && !fTriggeredEventMB)
+ return;
+
+ fn1->Fill(1);
+
+
+ if (fAnalysisMC) {
+
+
+
+ if (fAnalysisType == "ESD"){
+
+
+
+ AliHeader* headerMC = fMC->Header();
+ if (headerMC) {
+
+ AliGenEventHeader* genHeader = headerMC->GenEventHeader();
+ TArrayF vtxMC(3); // primary vertex MC
+ vtxMC[0]=9999; vtxMC[1]=9999; vtxMC[2]=9999; //initialize with dummy
+ if (genHeader) {
+ genHeader->PrimaryVertex(vtxMC);
+ }
+ fZvtxMC = vtxMC[2];
+
+ // PYTHIA:
+ AliGenPythiaEventHeader* pythiaGenHeader =
+ dynamic_cast<AliGenPythiaEventHeader*>(headerMC->GenEventHeader());
+ if (pythiaGenHeader) { //works only for pythia
+ fMcProcessType = GetPythiaEventProcessType(pythiaGenHeader->ProcessType());
+ }
+ // PHOJET:
+ AliGenDPMjetEventHeader* dpmJetGenHeader =
+ dynamic_cast<AliGenDPMjetEventHeader*>(headerMC->GenEventHeader());
+ if (dpmJetGenHeader) {
+ fMcProcessType = GetDPMjetEventProcessType(dpmJetGenHeader->ProcessType());
+ }
+ }
+ } else { // AOD
+
+
+
+ AliAODMCHeader* mcHeader = dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject("mcHeader"));
+
+
+ if(mcHeader) {
+ fZvtxMC = mcHeader->GetVtxZ();
+
+
+
+ if(strstr(mcHeader->GetGeneratorName(), "Pythia")) {
+ fMcProcessType = GetPythiaEventProcessType(mcHeader->GetEventType());
+ } else {
+ fMcProcessType = GetDPMjetEventProcessType(mcHeader->GetEventType());
+ }
+ }
+ }
+
+
+ }
+
+
+
+ if (fAnalysisType == "ESD"){
+
+ const AliESDVertex *vtxESD = fESD->GetPrimaryVertexTracks();
+ if(vtxESD->GetNContributors()<1) {
+ // SPD vertex
+ vtxESD = fESD->GetPrimaryVertexSPD();
+ /* quality checks on SPD-vertex */
+ TString vertexType = vtxESD->GetTitle();
+ if (vertexType.Contains("vertexer: Z") && (vtxESD->GetDispersion() > 0.04 || vtxESD->GetZRes() > 0.25))
+ fZvtx = -1599; //vertex = 0x0; //
+ else if (vtxESD->GetNContributors()<1)
+ fZvtx = -999; //vertex = 0x0; //
+ else
+ fZvtx = vtxESD->GetZ();
+ }
+ else
+ fZvtx = vtxESD->GetZ();
+
+ }
+ else // AOD
+ fZvtx = GetVertex(fAOD);
+
+ fVtxBeforeCuts->Fill(fZvtx);
+
+ //cut on the z position of vertex
+ if (TMath::Abs(fZvtx) > fVtxCut) {
+ return;
+ }
+ fn1->Fill(2);
+
+
+
+
+
+
+ Float_t centrality = -10;
+
+ // only analyze triggered events
+ if(fTriggeredEventMB) {
+
+ if (fAnalysisType == "ESD"){
+ if(fAnalysisPbPb){
+ AliCentrality *centObject = fESD->GetCentrality();
+ centrality = centObject->GetCentralityPercentile(fCentEst);
+
+ if((centrality>fMaxCent)||(centrality<fMinCent))return;
+ }
+ fcent->Fill(centrality);
+ fn1->Fill(3);
+ if(fAnalysisMC){
+ if(TMath::Abs(fZvtxMC)<fVtxCut){
+ ProcessMCTruthESD();
+ fVtxMC->Fill(fZvtxMC);
+ }
+ }
+ AnalyzeESD(fESD);
+ } else { // AOD
+ if(fAnalysisPbPb){
+ AliCentrality *centObject = fAOD->GetCentrality();
+ if(centObject){
+ centrality = centObject->GetCentralityPercentile(fCentEst);
+ }
+ //cout<<"centrality="<<centrality<<endl;
+ if((centrality>fMaxCent)||(centrality<fMinCent))return;
+ }
+ fcent->Fill(centrality);
+ fn1->Fill(3);
+ if(fAnalysisMC){
+ if(TMath::Abs(fZvtxMC)<fVtxCut){
+
+ ProcessMCTruthAOD();
+ fVtxMC->Fill(fZvtxMC);
+ }
+ }
+ AnalyzeAOD(fAOD);
+ }
+ }
+
+ fVtxAfterCuts->Fill(fZvtx);
+
+
+
+
+ // Post output data.
+ PostData(1, fListOfObjects);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskQAHighPtDeDx::AnalyzeESD(AliESDEvent* esdEvent)
+{
+ fRun = esdEvent->GetRunNumber();
+ fEventId = 0;
+ if(esdEvent->GetHeader())
+ fEventId = GetEventIdAsLong(esdEvent->GetHeader());
+
+ Bool_t isPileup = esdEvent->IsPileupFromSPD();
+ if(fPileUpRej)
+ if(isPileup)
+ return;
+ fn1->Fill(4);
+
+
+ // Int_t event = esdEvent->GetEventNumberInFile();
+ //UInt_t time = esdEvent->GetTimeStamp();
+ // ULong64_t trigger = esdEvent->GetTriggerMask();
+ magf = esdEvent->GetMagneticField();
+
+
+
+
+
+ if(fTriggeredEventMB) {// Only MC case can we have not triggered events
+
+ // accepted event
+ fEvents->Fill(0);
+
+
+ //Change, 10/04/13. Now accept all events to do a correct normalization
+ //if(fVtxStatus!=1) return; // accepted vertex
+ //Int_t nESDTracks = esdEvent->GetNumberOfTracks();
+
+ ProduceArrayTrksESD( esdEvent );//produce array with global track parameters
+ ProduceArrayV0ESD( esdEvent );//v0's
+
+
+ fEvents->Fill(1);
+
+
+
+
+ } // end if triggered
+
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskQAHighPtDeDx::AnalyzeAOD(AliAODEvent* aodEvent)
+{
+ fRun = aodEvent->GetRunNumber();
+ fEventId = 0;
+ if(aodEvent->GetHeader())
+ fEventId = GetEventIdAsLong(aodEvent->GetHeader());
+
+ //UInt_t time = 0; // Missing AOD info? aodEvent->GetTimeStamp();
+ magf = aodEvent->GetMagneticField();
+
+ //Int_t trackmult = 0; // no pt cuts
+ //Int_t nadded = 0;
+
+ Bool_t isPileup = aodEvent->IsPileupFromSPD();
+ if(fPileUpRej)
+ if(isPileup)
+ return;
+ fn1->Fill(4);
+
+
+
+ if(fTriggeredEventMB) {// Only MC case can we have not triggered events
+
+ // accepted event
+ fEvents->Fill(0);
+
+ //if(fVtxStatus!=1) return; // accepted vertex
+ //Int_t nAODTracks = aodEvent->GetNumberOfTracks();
+
+ ProduceArrayTrksAOD( aodEvent );
+ ProduceArrayV0AOD( aodEvent );
+
+ fEvents->Fill(1);
+
+
+
+
+ } // end if triggered
+
+}
+
+//_____________________________________________________________________________
+Float_t AliAnalysisTaskQAHighPtDeDx::GetVertex(const AliVEvent* event) const
+{
+ Float_t zvtx = -999;
+
+ const AliVVertex* primaryVertex = event->GetPrimaryVertex();
+
+ if(primaryVertex->GetNContributors()>0)
+ zvtx = primaryVertex->GetZ();
+
+ return zvtx;
+}
+
+//_____________________________________________________________________________
+Short_t AliAnalysisTaskQAHighPtDeDx::GetPidCode(Int_t pdgCode) const
+{
+ // return our internal code for pions, kaons, and protons
+
+ Short_t pidCode = 6;
+
+ switch (TMath::Abs(pdgCode)) {
+ case 211:
+ pidCode = 1; // pion
+ break;
+ case 321:
+ pidCode = 2; // kaon
+ break;
+ case 2212:
+ pidCode = 3; // proton
+ break;
+ case 11:
+ pidCode = 4; // electron
+ break;
+ case 13:
+ pidCode = 5; // muon
+ break;
+ default:
+ pidCode = 6; // something else?
+ };
+
+ return pidCode;
+}
+
+//_____________________________________________________________________________
+void AliAnalysisTaskQAHighPtDeDx::ProcessMCTruthESD()
+{
+ // Fill the special MC histogram with the MC truth info
+
+ const Int_t nTracksMC = fMCStack->GetNtrack();
+
+ for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
+
+ //Cuts
+ if(!(fMCStack->IsPhysicalPrimary(iTracks)))
+ continue;
+
+ TParticle* trackMC = fMCStack->Particle(iTracks);
+
+ TParticlePDG* pdgPart = trackMC ->GetPDG();
+ Double_t chargeMC = pdgPart->Charge();
+
+ if(chargeMC==0)
+ continue;
+
+ if (TMath::Abs(trackMC->Eta()) > fEtaCut )
+ continue;
+
+ Double_t etaMC = trackMC->Eta();
+ Int_t pdgCode = trackMC->GetPdgCode();
+ Short_t pidCodeMC = 0;
+ pidCodeMC = GetPidCode(pdgCode);
+
+
+ for(Int_t nh = 0; nh < 9; nh++) {
+
+ if( etaMC > etaHigh[nh]/10.0 || etaMC < etaLow[nh]/10.0 )
+ continue;
+
+ hMcIn[0][nh]->Fill(trackMC->Pt());
+ hMcIn[pidCodeMC][nh]->Fill(trackMC->Pt());
+
+
+ }
+
+ }//MC track loop
+
+
+
+}
+
+//_____________________________________________________________________________
+void AliAnalysisTaskQAHighPtDeDx::ProcessMCTruthAOD()
+{
+ // Fill the special MC histogram with the MC truth info
+
+
+ const Int_t nTracksMC = fMCArray->GetEntriesFast();
+
+ for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
+
+ AliAODMCParticle* trackMC = dynamic_cast<AliAODMCParticle*>(fMCArray->At(iTracks));
+
+ //Cuts
+ if(!(trackMC->IsPhysicalPrimary()))
+ continue;
+
+
+ Double_t chargeMC = trackMC->Charge();
+ if(chargeMC==0)
+ continue;
+
+
+ if (TMath::Abs(trackMC->Eta()) > fEtaCut )
+ continue;
+
+ Double_t etaMC = trackMC->Eta();
+ Int_t pdgCode = trackMC->GetPdgCode();
+ Short_t pidCodeMC = 0;
+ pidCodeMC = GetPidCode(pdgCode);
+
+ //cout<<"pidcode="<<pidCodeMC<<endl;
+ for(Int_t nh = 0; nh < 9; nh++) {
+
+ if( etaMC > etaHigh[nh]/10.0 || etaMC < etaLow[nh]/10.0 )
+ continue;
+
+ hMcIn[0][nh]->Fill(trackMC->Pt());
+ hMcIn[pidCodeMC][nh]->Fill(trackMC->Pt());
+
+
+ }
+
+ }//MC track loop
+
+
+}
+
+
+//_____________________________________________________________________________
+Short_t AliAnalysisTaskQAHighPtDeDx::GetPythiaEventProcessType(Int_t pythiaType) {
+ //
+ // Get the process type of the event. PYTHIA
+ //
+ // source PWG0 dNdpt
+
+ Short_t globalType = -1; //init
+
+ if(pythiaType==92||pythiaType==93){
+ globalType = 2; //single diffractive
+ }
+ else if(pythiaType==94){
+ globalType = 3; //double diffractive
+ }
+ //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??
+ else {
+ globalType = 1; //non diffractive
+ }
+ return globalType;
+}
+
+//_____________________________________________________________________________
+Short_t AliAnalysisTaskQAHighPtDeDx::GetDPMjetEventProcessType(Int_t dpmJetType) {
+ //
+ // get the process type of the event. PHOJET
+ //
+ //source PWG0 dNdpt
+ // can only read pythia headers, either directly or from cocktalil header
+ Short_t globalType = -1;
+
+ if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
+ globalType = 1;
+ }
+ else if (dpmJetType==5 || dpmJetType==6) {
+ globalType = 2;
+ }
+ else if (dpmJetType==7) {
+ globalType = 3;
+ }
+ return globalType;
+}
+
+//_____________________________________________________________________________
+ULong64_t AliAnalysisTaskQAHighPtDeDx::GetEventIdAsLong(AliVHeader* header) const
+{
+ // To have a unique id for each event in a run!
+ // Modified from AliRawReader.h
+ return ((ULong64_t)header->GetBunchCrossNumber()+
+ (ULong64_t)header->GetOrbitNumber()*3564+
+ (ULong64_t)header->GetPeriodNumber()*16777215*3564);
+}
+
+
+//____________________________________________________________________
+TParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMother(AliStack* stack, Int_t label)
+{
+ //
+ // Finds the first mother among the primary particles of the particle identified by <label>,
+ // i.e. the primary that "caused" this particle
+ //
+ // Taken from AliPWG0Helper class
+ //
+
+ Int_t motherLabel = FindPrimaryMotherLabel(stack, label);
+ if (motherLabel < 0)
+ return 0;
+
+ return stack->Particle(motherLabel);
+}
+
+//____________________________________________________________________
+Int_t AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherLabel(AliStack* stack, Int_t label)
+{
+ //
+ // Finds the first mother among the primary particles of the particle identified by <label>,
+ // i.e. the primary that "caused" this particle
+ //
+ // returns its label
+ //
+ // Taken from AliPWG0Helper class
+ //
+ const Int_t nPrim = stack->GetNprimary();
+
+ while (label >= nPrim) {
+
+ //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
+
+ TParticle* particle = stack->Particle(label);
+ if (!particle) {
+
+ AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
+ return -1;
+ }
+
+ // find mother
+ if (particle->GetMother(0) < 0) {
+
+ AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
+ return -1;
+ }
+
+ label = particle->GetMother(0);
+ }
+
+ return label;
+}
+
+//____________________________________________________________________
+AliAODMCParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherAOD(AliAODMCParticle* startParticle)
+{
+ //
+ // Finds the first mother among the primary particles of the particle identified by <label>,
+ // i.e. the primary that "caused" this particle
+ //
+ // Taken from AliPWG0Helper class
+ //
+
+ AliAODMCParticle* mcPart = startParticle;
+
+ while (mcPart)
+ {
+
+ if(mcPart->IsPrimary())
+ return mcPart;
+
+ Int_t mother = mcPart->GetMother();
+
+ mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
+ }
+
+ return 0;
+}
+
+
+//V0______________________________________
+//____________________________________________________________________
+TParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherV0(AliStack* stack, Int_t label)
+{
+ //
+ // Finds the first mother among the primary particles of the particle identified by <label>,
+ // i.e. the primary that "caused" this particle
+ //
+ // Taken from AliPWG0Helper class
+ //
+
+ Int_t nSteps = 0;
+
+ Int_t motherLabel = FindPrimaryMotherLabelV0(stack, label, nSteps);
+ if (motherLabel < 0)
+ return 0;
+
+ return stack->Particle(motherLabel);
+}
+
+//____________________________________________________________________
+Int_t AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherLabelV0(AliStack* stack, Int_t label, Int_t& nSteps)
+{
+ //
+ // Finds the first mother among the primary particles of the particle identified by <label>,
+ // i.e. the primary that "caused" this particle
+ //
+ // returns its label
+ //
+ // Taken from AliPWG0Helper class
+ //
+ nSteps = 0;
+ const Int_t nPrim = stack->GetNprimary();
+
+ while (label >= nPrim) {
+
+ //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
+
+ nSteps++; // 1 level down
+
+ TParticle* particle = stack->Particle(label);
+ if (!particle) {
+
+ AliDebugGeneral("FindPrimaryMotherLabelV0", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
+ return -1;
+ }
+
+ // find mother
+ if (particle->GetMother(0) < 0) {
+
+ AliDebugGeneral("FindPrimaryMotherLabelV0", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
+ return -1;
+ }
+
+ label = particle->GetMother(0);
+ }
+
+ return label;
+}
+
+//____________________________________________________________________
+AliAODMCParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherAODV0(AliAODMCParticle* startParticle, Int_t& nSteps)
+{
+ //
+ // Finds the first mother among the primary particles of the particle identified by <label>,
+ // i.e. the primary that "caused" this particle
+ //
+ // Taken from AliPWG0Helper class
+ //
+
+ nSteps = 0;
+
+ AliAODMCParticle* mcPart = startParticle;
+
+ while (mcPart)
+ {
+
+ if(mcPart->IsPrimary())
+ return mcPart;
+
+ Int_t mother = mcPart->GetMother();
+
+ mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
+ nSteps++; // 1 level down
+ }
+
+ return 0;
+}
+
+
+
+//__________________________________________________________________
+void AliAnalysisTaskQAHighPtDeDx::ProduceArrayTrksESD( AliESDEvent *ESDevent ){
+
+ const Int_t nESDTracks = ESDevent->GetNumberOfTracks();
+ //Int_t trackmult=0;
+
+
+ Int_t multTPC = 0;
+
+ //get multiplicity tpc only track cuts
+ for(Int_t iT = 0; iT < nESDTracks; iT++) {
+
+ AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
+
+
+ if(TMath::Abs(esdTrack->Eta()) > fEtaCut)
+ continue;
+
+ //only golden track cuts
+ UInt_t selectDebug = 0;
+ if (fTrackFilterTPC) {
+ selectDebug = fTrackFilterTPC->IsSelected(esdTrack);
+ if (!selectDebug) {
+ //cout<<"this is not a golden track"<<endl;
+ continue;
+ }
+ }
+
+ multTPC++;
+
+ }
+
+
+
+ for(Int_t iT = 0; iT < nESDTracks; iT++) {
+
+ AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
+
+
+ if(TMath::Abs(esdTrack->Eta()) > fEtaCut)
+ continue;
+
+ //only golden track cuts
+ UInt_t selectDebug = 0;
+ if (fTrackFilterGolden) {
+ selectDebug = fTrackFilterGolden->IsSelected(esdTrack);
+ if (!selectDebug) {
+ //cout<<"this is not a golden track"<<endl;
+ continue;
+ }
+ }
+
+
+ //Int_t sharedtpcclusters = esdTrack->GetTPCnclsS();
+ Short_t ncl = esdTrack->GetTPCsignalN();
+
+
+ if(ncl<70)
+ continue;
+ Double_t eta = esdTrack->Eta();
+ Double_t phi = esdTrack->Phi();
+ Double_t momentum = esdTrack->P();
+ Float_t dedx = esdTrack->GetTPCsignal();
+
+ //cout<<"magf="<<magf<<endl;
+
+ if(!PhiCut(esdTrack->Pt(), phi, esdTrack->Charge(), magf, cutLow, cutHigh))
+ continue;
+
+
+
+
+ //TOF
+
+ Bool_t IsTOFout=kFALSE;
+ if ((esdTrack->GetStatus()&AliESDtrack::kTOFout)==0)
+ IsTOFout=kTRUE;
+ Float_t lengthtrack=esdTrack->GetIntegratedLength();
+ Float_t timeTOF=esdTrack->GetTOFsignal();
+ Double_t inttime[5]={0,0,0,0,0};
+ esdTrack->GetIntegratedTimes(inttime);// Returns the array with integrated times for each particle hypothesis
+ Float_t beta = -99;
+ if ( !IsTOFout ){
+ if ( ( lengthtrack != 0 ) && ( timeTOF != 0) )
+ beta = inttime[0] / timeTOF;
+ }
+
+ Short_t pidCode = 0;
+
+ if(fAnalysisMC) {
+
+ const Int_t label = TMath::Abs(esdTrack->GetLabel());
+ TParticle* mcTrack = fMCStack->Particle(label);
+ if (mcTrack){
+
+ Int_t pdgCode = mcTrack->GetPdgCode();
+ pidCode = GetPidCode(pdgCode);
+
+ }
+
+ }
+
+
+ if( momentum <= 0.6 && momentum >= 0.4 ){//only p:0.4-0.6 GeV, pion MIP
+
+ if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
+ if(momentum<0.6&&momentum>0.4){
+ hMIPVsEta->Fill(eta,dedx);
+ pMIPVsEta->Fill(eta,dedx);
+ }
+ }
+ if( dedx > DeDxMIPMax+5 && dedx < 95 ){
+ if(TMath::Abs(beta-1)<0.1){
+ hPlateauVsEta->Fill(eta,dedx);
+ pPlateauVsEta->Fill(eta,dedx);
+ }
+ }
+ }
+
+
+ for(Int_t nh = 0; nh < 9; nh++) {
+
+ if( eta > etaHigh[nh]/10.0 || eta < etaLow[nh]/10.0 )
+ continue;
+
+
+ if(fAnalysisMC){
+ hMcOut[0][nh]->Fill(esdTrack->Pt());
+ hMcOut[pidCode][nh]->Fill(esdTrack->Pt());
+ }
+
+ histAllCh[nh]->Fill(momentum, dedx);
+ if(beta>1)
+ histPiTof[nh]->Fill(momentum, dedx);
+
+ if( momentum <= 0.6 && momentum >= 0.4 ){//only p:0.4-0.6 GeV, pion MIP
+ //Fill pion MIP, before calibration
+ if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
+ hMIPVsPhi[nh]->Fill(phi,dedx);
+ pMIPVsPhi[nh]->Fill(phi,dedx);
+
+
+ hMIPVsNch[nh]->Fill(multTPC,dedx);
+ pMIPVsNch[nh]->Fill(multTPC,dedx);
+
+ }
+
+ //Fill electrons, before calibration
+ if( dedx > DeDxMIPMax+5 && dedx < 95 ){
+ if(TMath::Abs(beta-1)<0.1){
+ hPlateauVsPhi[nh]->Fill(phi,dedx);
+ pPlateauVsPhi[nh]->Fill(phi,dedx);
+ }
+ }
+
+ }
+
+
+ }//end loop over eta intervals
+
+
+
+
+
+ }//end of track loop
+
+
+
+
+}
+//__________________________________________________________________
+void AliAnalysisTaskQAHighPtDeDx::ProduceArrayTrksAOD( AliAODEvent *AODevent ){
+
+
+ Int_t nAODTracks = AODevent->GetNumberOfTracks();
+ Int_t multTPC = 0;
+
+ //get multiplicity tpc only track cuts
+ for(Int_t iT = 0; iT < nAODTracks; iT++) {
+
+ AliAODTrack* aodTrack = AODevent->GetTrack(iT);
+
+ if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
+ continue;
+
+
+ if (fTrackFilterTPC) {
+ // TPC only cuts is bit 1
+ if(!aodTrack->TestFilterBit(1))
+ continue;
+ }
+
+ multTPC++;
+
+ }
+
+
+ for(Int_t iT = 0; iT < nAODTracks; iT++) {
+
+ AliAODTrack* aodTrack = AODevent->GetTrack(iT);
+
+ if (fTrackFilterGolden) {
+ // "Global track RAA analysis QM2011 + Chi2ITS<36"; bit 1024
+ if(!aodTrack->TestFilterBit(1024))
+ continue;
+ }
+
+
+ if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
+ continue;
+
+
+ Double_t eta = aodTrack->Eta();
+ Double_t phi = aodTrack->Phi();
+ Double_t momentum = aodTrack->P();
+
+
+ if(!PhiCut(aodTrack->Pt(), phi, aodTrack->Charge(), magf, cutLow, cutHigh))
+ continue;
+
+
+
+ AliAODPid* aodPid = aodTrack->GetDetPid();
+ Short_t ncl = -10;
+ Float_t dedx = -10;
+
+ //TOF
+ Float_t beta = -99;
+
+
+ if(aodPid) {
+ ncl = aodPid->GetTPCsignalN();
+ dedx = aodPid->GetTPCsignal();
+ //TOF
+ Bool_t IsTOFout=kFALSE;
+ Float_t lengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF
+ Float_t timeTOF=-999;
+
+ if ((aodTrack->GetStatus()&AliESDtrack::kTOFout)==0)
+ IsTOFout=kTRUE;
+
+ lengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF
+
+ timeTOF=aodPid->GetTOFsignal();
+
+ Double_t inttime[5]={0,0,0,0,0};
+ aodPid->GetIntegratedTimes(inttime);// Returns the array with integrated times for each particle hypothesis
+
+
+ if ( !IsTOFout ){
+ if ( ( lengthtrack != 0 ) && ( timeTOF != 0) )
+ beta = inttime[0] / timeTOF;
+ }
+
+ }
+
+
+ if(ncl<70)
+ continue;
+
+
+ Short_t pidCode = 0;
+
+ if(fAnalysisMC) {
+
+ const Int_t label = TMath::Abs(aodTrack->GetLabel());
+ AliAODMCParticle* mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(label));
+
+ if (mcTrack){
+
+ Int_t pdgCode = mcTrack->GetPdgCode();
+ pidCode = GetPidCode(pdgCode);
+
+ }
+
+ }
+
+
+
+
+ //if(momentum>0.6||momentum<0.4)//only p:0.4-0.6 GeV
+ //continue;
+
+ //etaLow
+ //etaHigh
+ if( momentum <= 0.6 && momentum >= 0.4 ){//only p:0.4-0.6 GeV, pion MIP
+ if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
+ if(momentum<0.6&&momentum>0.4){
+ hMIPVsEta->Fill(eta,dedx);
+ pMIPVsEta->Fill(eta,dedx);
+ }
+ }
+ if( dedx > DeDxMIPMax+5 && dedx < 95 ){
+ if(TMath::Abs(beta-1)<0.1){
+ hPlateauVsEta->Fill(eta,dedx);
+ pPlateauVsEta->Fill(eta,dedx);
+ }
+ }
+ }
+
+
+ for(Int_t nh = 0; nh < 9; nh++) {
+
+ if( eta > etaHigh[nh]/10.0 || eta < etaLow[nh]/10.0 )
+ continue;
+
+
+ if(fAnalysisMC){
+ hMcOut[0][nh]->Fill(aodTrack->Pt());
+ hMcOut[pidCode][nh]->Fill(aodTrack->Pt());
+ }
+
+ histAllCh[nh]->Fill(momentum, dedx);
+ if(beta>1)
+ histPiTof[nh]->Fill(momentum, dedx);
+
+ if( momentum <= 0.6 && momentum >= 0.4 ){//only p:0.4-0.6 GeV, pion MIP
+ //Fill pion MIP, before calibration
+ if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
+ hMIPVsPhi[nh]->Fill(phi,dedx);
+ pMIPVsPhi[nh]->Fill(phi,dedx);
+
+
+ hMIPVsNch[nh]->Fill(multTPC,dedx);
+ pMIPVsNch[nh]->Fill(multTPC,dedx);
+
+ }
+
+ //Fill electrons, before calibration
+ if( dedx > DeDxMIPMax+5 && dedx < 95 ){
+ if(TMath::Abs(beta-1)<0.1){
+ hPlateauVsPhi[nh]->Fill(phi,dedx);
+ pPlateauVsPhi[nh]->Fill(phi,dedx);
+ }
+ }
+
+ }
+
+ }//end loop over eta intervals
+
+
+
+
+
+ }//end of track loop
+
+
+
+
+
+
+
+}
+//__________________________________________________
+void AliAnalysisTaskQAHighPtDeDx::ProduceArrayV0ESD( AliESDEvent *ESDevent ){
+ Int_t nv0s = ESDevent->GetNumberOfV0s();
+ /*
+ if(nv0s<1){
+ return;
+ }*/
+
+ const AliESDVertex *myBestPrimaryVertex = ESDevent->GetPrimaryVertex();
+ if (!myBestPrimaryVertex) return;
+ if (!(myBestPrimaryVertex->GetStatus())) return;
+ Double_t lPrimaryVtxPosition[3];
+ myBestPrimaryVertex->GetXYZ(lPrimaryVtxPosition);
+
+ Double_t lPrimaryVtxCov[6];
+ myBestPrimaryVertex->GetCovMatrix(lPrimaryVtxCov);
+ Double_t lPrimaryVtxChi2 = myBestPrimaryVertex->GetChi2toNDF();
+
+ AliAODVertex* myPrimaryVertex = new AliAODVertex(lPrimaryVtxPosition, lPrimaryVtxCov, lPrimaryVtxChi2, NULL, -1, AliAODVertex::kPrimary);
+
+
+
+ //
+ // LOOP OVER V0s, K0s, L, AL
+ //
+
+
+ for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
+
+ // This is the begining of the V0 loop
+ AliESDv0 *esdV0 = ESDevent->GetV0(iV0);
+ if (!esdV0) continue;
+
+ //check onfly status
+ if(esdV0->GetOnFlyStatus()!=0)
+ continue;
+
+
+ // AliESDTrack (V0 Daughters)
+ UInt_t lKeyPos = (UInt_t)TMath::Abs(esdV0->GetPindex());
+ UInt_t lKeyNeg = (UInt_t)TMath::Abs(esdV0->GetNindex());
+
+ AliESDtrack *pTrack = ESDevent->GetTrack(lKeyPos);
+ AliESDtrack *nTrack = ESDevent->GetTrack(lKeyNeg);
+ if (!pTrack || !nTrack) {
+ Printf("ERROR: Could not retreive one of the daughter track");
+ continue;
+ }
+
+ // Remove like-sign
+ if (pTrack->GetSign() == nTrack->GetSign()) {
+ //cout<< "like sign, continue"<< endl;
+ continue;
+ }
+
+ // Eta cut on decay products
+ if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
+ continue;
+
+
+ // Check if switch does anything!
+ Bool_t isSwitched = kFALSE;
+ if (pTrack->GetSign() < 0) { // switch
+
+ isSwitched = kTRUE;
+ AliESDtrack* helpTrack = nTrack;
+ nTrack = pTrack;
+ pTrack = helpTrack;
+ }
+
+ //Double_t lV0Position[3];
+ //esdV0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[2]);
+
+
+ AliKFVertex primaryVtxKF( *myPrimaryVertex );
+ AliKFParticle::SetField(ESDevent->GetMagneticField());
+
+ // Also implement switch here!!!!!!
+ AliKFParticle* negEKF = 0; // e-
+ AliKFParticle* posEKF = 0; // e+
+ AliKFParticle* negPiKF = 0; // pi -
+ AliKFParticle* posPiKF = 0; // pi +
+ AliKFParticle* posPKF = 0; // p
+ AliKFParticle* negAPKF = 0; // p-bar
+
+ if(!isSwitched) {
+ negEKF = new AliKFParticle( *(esdV0->GetParamN()) , 11);
+ posEKF = new AliKFParticle( *(esdV0->GetParamP()) ,-11);
+ negPiKF = new AliKFParticle( *(esdV0->GetParamN()) ,-211);
+ posPiKF = new AliKFParticle( *(esdV0->GetParamP()) , 211);
+ posPKF = new AliKFParticle( *(esdV0->GetParamP()) , 2212);
+ negAPKF = new AliKFParticle( *(esdV0->GetParamN()) ,-2212);
+ } else { // switch + and -
+ negEKF = new AliKFParticle( *(esdV0->GetParamP()) , 11);
+ posEKF = new AliKFParticle( *(esdV0->GetParamN()) ,-11);
+ negPiKF = new AliKFParticle( *(esdV0->GetParamP()) ,-211);
+ posPiKF = new AliKFParticle( *(esdV0->GetParamN()) , 211);
+ posPKF = new AliKFParticle( *(esdV0->GetParamN()) , 2212);
+ negAPKF = new AliKFParticle( *(esdV0->GetParamP()) ,-2212);
+ }
+
+ AliKFParticle v0GKF; // Gamma e.g. from pi0
+ v0GKF+=(*negEKF);
+ v0GKF+=(*posEKF);
+ v0GKF.SetProductionVertex(primaryVtxKF);
+
+ AliKFParticle v0K0sKF; // K0 short
+ v0K0sKF+=(*negPiKF);
+ v0K0sKF+=(*posPiKF);
+ v0K0sKF.SetProductionVertex(primaryVtxKF);
+
+ AliKFParticle v0LambdaKF; // Lambda
+ v0LambdaKF+=(*negPiKF);
+ v0LambdaKF+=(*posPKF);
+ v0LambdaKF.SetProductionVertex(primaryVtxKF);
+
+ AliKFParticle v0AntiLambdaKF; // Lambda-bar
+ v0AntiLambdaKF+=(*posPiKF);
+ v0AntiLambdaKF+=(*negAPKF);
+ v0AntiLambdaKF.SetProductionVertex(primaryVtxKF);
+
+ Double_t dmassG = v0GKF.GetMass();
+ Double_t dmassK = v0K0sKF.GetMass()-0.498;
+ Double_t dmassL = v0LambdaKF.GetMass()-1.116;
+ Double_t dmassAL = v0AntiLambdaKF.GetMass()-1.116;
+
+
+ for( Int_t case_v0 = 0; case_v0 < 2; ++case_v0 ){
+
+
+ switch(case_v0){
+ case 0:{
+
+ Bool_t fillPos = kFALSE;
+ Bool_t fillNeg = kFALSE;
+
+ if(dmassG < 0.1)
+ continue;
+
+ if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01){
+ continue;
+ }
+
+ if(dmassL<0.01){
+ fillPos = kTRUE;
+ }
+ if(dmassAL<0.01) {
+ if(fillPos)
+ continue;
+ fillNeg = kTRUE;
+ }
+
+ if(dmassK<0.01) {
+ if(fillPos||fillNeg)
+ continue;
+ fillPos = kTRUE;
+ fillNeg = kTRUE;
+ }
+
+
+ for(Int_t j = 0; j < 2; j++) {
+
+ AliESDtrack* track = 0;
+
+ if(j==0) {
+
+ if(fillNeg)
+ track = nTrack;
+ else
+ continue;
+ } else {
+
+ if(fillPos)
+ track = pTrack;
+ else
+ continue;
+ }
+
+ if(track->GetTPCsignalN()<=70)continue;
+ Double_t phi = track->Phi();
+
+ if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
+ continue;
+
+
+ //if(!PhiCut(pt, phi, charge, magf, cutLow, cutHigh, hPhi))
+ // continue;
+
+ Double_t eta = track->Eta();
+ Double_t momentum = track->Pt();
+ Double_t dedx = track->GetTPCsignal();
+
+ if(fillPos&&fillNeg){
+
+
+ if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
+ if(momentum<0.6&&momentum>0.4){
+ hMIPVsEtaV0s->Fill(eta,dedx);
+ pMIPVsEtaV0s->Fill(eta,dedx);
+ }
+ }
+
+
+ }
+
+ for(Int_t nh = 0; nh < nHists; nh++) {
+
+
+
+ if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
+ continue;
+
+ if(fillPos&&fillNeg){
+
+ histPiV0[nh]->Fill(momentum, dedx);
+
+ }
+ else{
+
+ histPV0[nh]->Fill(momentum, dedx);
+
+ }
+
+ }
+
+ }//end loop over two tracks
+
+ };
+ break;
+
+ case 1:{//gammas
+
+ Bool_t fillPos = kFALSE;
+ Bool_t fillNeg = kFALSE;
+
+
+
+ if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01) {
+ if(dmassG<0.01 && dmassG>0.0001) {
+
+
+ if( TMath::Abs(nTrack->GetTPCsignal() - 85.0) < 5)
+ fillPos = kTRUE;
+ if( TMath::Abs(pTrack->GetTPCsignal() - 85.0) < 5)
+ fillNeg = kTRUE;
+
+ } else {
+ continue;
+ }
+ }
+
+ //cout<<"fillPos="<<fillPos<<endl;
+
+ if(fillPos == kTRUE && fillNeg == kTRUE)
+ continue;
+
+
+ AliESDtrack* track = 0;
+ if(fillNeg)
+ track = nTrack;
+ else if(fillPos)
+ track = pTrack;
+ else
+ continue;
+
+ Double_t dedx = track->GetTPCsignal();
+ Double_t eta = track->Eta();
+ Double_t phi = track->Phi();
+ Double_t momentum = track->P();
+
+ if(track->GetTPCsignalN()<=70)continue;
+
+ if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
+ continue;
+
+ for(Int_t nh = 0; nh < nHists; nh++) {
+
+ if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
+ continue;
+
+ histEV0[nh]->Fill(momentum, dedx);
+
+ }
+
+ };
+ break;
+
+
+ }//end switch
+
+ }//end loop over case V0
+
+
+ // clean up loop over v0
+
+ delete negPiKF;
+ delete posPiKF;
+ delete posPKF;
+ delete negAPKF;
+
+
+
+ }
+
+
+ delete myPrimaryVertex;
+
+
+}
+//__________________________________________________________________________
+void AliAnalysisTaskQAHighPtDeDx::ProduceArrayV0AOD( AliAODEvent *AODevent ){
+ Int_t nv0s = AODevent->GetNumberOfV0s();
+ /*
+ if(nv0s<1){
+ return;
+ }*/
+
+ AliAODVertex *myBestPrimaryVertex = AODevent->GetPrimaryVertex();
+ if (!myBestPrimaryVertex) return;
+
+
+
+ // ################################
+ // #### BEGINNING OF V0 CODE ######
+ // ################################
+ // This is the begining of the V0 loop
+ for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
+ AliAODv0 *aodV0 = AODevent->GetV0(iV0);
+ if (!aodV0) continue;
+
+
+ //check onfly status
+ if(aodV0->GetOnFlyStatus()!=0)
+ continue;
+
+ // AliAODTrack (V0 Daughters)
+ AliAODVertex* vertex = aodV0->GetSecondaryVtx();
+ if (!vertex) {
+ Printf("ERROR: Could not retrieve vertex");
+ continue;
+ }
+
+ AliAODTrack *pTrack = (AliAODTrack*)vertex->GetDaughter(0);
+ AliAODTrack *nTrack = (AliAODTrack*)vertex->GetDaughter(1);
+ if (!pTrack || !nTrack) {
+ Printf("ERROR: Could not retrieve one of the daughter track");
+ continue;
+ }
+
+ // Remove like-sign
+ if (pTrack->Charge() == nTrack->Charge()) {
+ //cout<< "like sign, continue"<< endl;
+ continue;
+ }
+
+ // Make sure charge ordering is ok
+ if (pTrack->Charge() < 0) {
+ AliAODTrack* helpTrack = pTrack;
+ pTrack = nTrack;
+ nTrack = helpTrack;
+ }
+
+ // Eta cut on decay products
+ if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
+ continue;
+
+
+ Double_t dmassG = aodV0->InvMass2Prongs(0,1,11,11);
+ Double_t dmassK = aodV0->MassK0Short()-0.498;
+ Double_t dmassL = aodV0->MassLambda()-1.116;
+ Double_t dmassAL = aodV0->MassAntiLambda()-1.116;
+
+ for( Int_t case_v0 = 0; case_v0 < 2; ++case_v0 ){
+
+
+ switch(case_v0){
+ case 0:{
+ Bool_t fillPos = kFALSE;
+ Bool_t fillNeg = kFALSE;
+
+
+ if(dmassG < 0.1)
+ continue;
+
+ if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01){
+ continue;
+ }
+
+ if(dmassL<0.01){
+ fillPos = kTRUE;
+ }
+ if(dmassAL<0.01) {
+ if(fillPos)
+ continue;
+ fillNeg = kTRUE;
+ }
+
+ if(dmassK<0.01) {
+ if(fillPos||fillNeg)
+ continue;
+ fillPos = kTRUE;
+ fillNeg = kTRUE;
+ }
+
+
+
+ for(Int_t j = 0; j < 2; j++) {
+
+ AliAODTrack* track = 0;
+
+ if(j==0) {
+
+ if(fillNeg)
+ track = nTrack;
+ else
+ continue;
+ } else {
+
+ if(fillPos)
+ track = pTrack;
+ else
+ continue;
+ }
+
+ if(track->GetTPCsignalN()<=70)continue;
+
+ Double_t phi = track->Phi();
+
+ if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
+ continue;
+
+
+ //if(!PhiCut(pt, phi, charge, magf, cutLow, cutHigh, hPhi))
+ // continue;
+
+ Double_t eta = track->Eta();
+ Double_t momentum = track->Pt();
+ Double_t dedx = track->GetTPCsignal();
+
+ if(fillPos&&fillNeg){
+
+
+ if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
+ if(momentum<0.6&&momentum>0.4){
+ hMIPVsEtaV0s->Fill(eta,dedx);
+ pMIPVsEtaV0s->Fill(eta,dedx);
+ }
+ }
+
+
+ }
+
+ for(Int_t nh = 0; nh < nHists; nh++) {
+
+
+
+ if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
+ continue;
+
+ if(fillPos&&fillNeg){
+
+ histPiV0[nh]->Fill(momentum, dedx);
+
+ }
+ else{
+
+ histPV0[nh]->Fill(momentum, dedx);
+
+ }
+
+ }
+
+
+ }//end loop over two tracks
+ };
+ break;
+
+ case 1:{//gammas
+
+ Bool_t fillPos = kFALSE;
+ Bool_t fillNeg = kFALSE;
+
+ if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01) {
+ if(dmassG<0.01 && dmassG>0.0001) {
+
+ if( TMath::Abs(nTrack->GetTPCsignal() - 85.0) < 5)
+ fillPos = kTRUE;
+ if( TMath::Abs(pTrack->GetTPCsignal() - 85.0) < 5)
+ fillNeg = kTRUE;
+
+ } else {
+ continue;
+ }
+ }
+
+
+ if(fillPos == kTRUE && fillNeg == kTRUE)
+ continue;
+
+
+ AliAODTrack* track = 0;
+ if(fillNeg)
+ track = nTrack;
+ else if(fillPos)
+ track = pTrack;
+ else
+ continue;
+
+ Double_t dedx = track->GetTPCsignal();
+ Double_t eta = track->Eta();
+ Double_t phi = track->Phi();
+ Double_t momentum = track->P();
+
+ if(track->GetTPCsignalN()<=70)continue;
+
+ if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
+ continue;
+
+ for(Int_t nh = 0; nh < nHists; nh++) {
+
+ if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
+ continue;
+
+ histEV0[nh]->Fill(momentum, dedx);
+
+ }
+
+ };
+ break;
+
+
+ }//end switch
+ }//end loop over V0s cases
+
+ }//end loop over v0's
+
+
+
+
+}
+Bool_t AliAnalysisTaskQAHighPtDeDx::PhiCut(Double_t pt, Double_t phi, Double_t q, Float_t mag, TF1* phiCutLow, TF1* phiCutHigh)
+{
+ if(pt < 2.0)
+ return kTRUE;
+
+ //Double_t phi = track->Phi();
+ if(mag < 0) // for negatve polarity field
+ phi = TMath::TwoPi() - phi;
+ if(q < 0) // for negatve charge
+ phi = TMath::TwoPi()-phi;
+
+ phi += TMath::Pi()/18.0; // to center gap in the middle
+ phi = fmod(phi, TMath::Pi()/9.0);
+
+ if(phi<phiCutHigh->Eval(pt)
+ && phi>phiCutLow->Eval(pt))
+ return kFALSE; // reject track
+
+ hPhi->Fill(pt, phi);
+
+ return kTRUE;
+}
--- /dev/null
+#ifndef ALIANALYSISTASKQAHIGHPTDEDX_H
+#define ALIANALYSISTASKQAHIGHPTDEDX_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+* See cxx source for full Copyright notice */
+/* $Id$ */
+
+
+// ROOT includes
+#include <TList.h>
+#include <TH1.h>
+#include <TProfile.h>
+#include <TTreeStream.h>
+#include <TRandom.h>
+#include <TObject.h>
+
+// AliRoot includes
+#include <AliAnalysisTaskSE.h>
+#include <AliESDEvent.h>
+#include <AliAODEvent.h>
+#include <AliMCEvent.h>
+#include <AliAnalysisFilter.h>
+#include <AliStack.h>
+#include <AliGenEventHeader.h>
+#include <AliVHeader.h>
+#include <AliAODMCParticle.h>
+#include <AliESDtrackCuts.h>
+
+
+
+class AliAnalysisTaskQAHighPtDeDx : public AliAnalysisTaskSE {
+ public:
+
+ //AliAnalysisTaskQAHighPtDeDx();
+ //AliAnalysisTaskQAHighPtDeDx(const char *name);
+ //virtual ~AliAnalysisTaskQAHighPtDeDx();
+ AliAnalysisTaskQAHighPtDeDx(const char *name="<default name>");
+ virtual ~AliAnalysisTaskQAHighPtDeDx() { /*if (fOutputList) delete fOutputList;*/}//;
+
+ virtual void UserCreateOutputObjects();
+ virtual void UserExec(Option_t *option);
+
+ Bool_t GetAnalysisMC() { return fAnalysisMC; }
+ Double_t GetVtxCut() { return fVtxCut; }
+ Double_t GetEtaCut() { return fEtaCut; }
+ //Double_t GetMinPt() { return fMinPt; }
+ //Int_t GetTreeOption() { return fTreeOption; }
+
+ virtual void SetTrigger(UInt_t ktriggerInt) {ftrigBit = ktriggerInt;}
+ virtual void SetTrackFilterGolden(AliAnalysisFilter* trackF) {fTrackFilterGolden = trackF;}
+ virtual void SetTrackFilterTPC(AliAnalysisFilter* trackF) {fTrackFilterTPC = trackF;}
+ virtual void SetCentralityEstimator(const char * centEst) {fCentEst = centEst;}
+ virtual void SetAnalysisType(const char* analysisType) {fAnalysisType = analysisType;}
+ virtual void SetAnalysisMC(Bool_t isMC) {fAnalysisMC = isMC;}
+ virtual void SetVtxCut(Double_t vtxCut){fVtxCut = vtxCut;}
+ virtual void SetEtaCut(Double_t etaCut){fEtaCut = etaCut;}
+ virtual void SetPileUpRej(Bool_t isrej) {fPileUpRej = isrej;}
+ virtual void SetMinCent(Float_t minvalc) {fMinCent = minvalc;}
+ virtual void SetMaxCent(Float_t maxvalc) {fMaxCent = maxvalc;}
+ virtual void SetStoreMcIn(Bool_t value) {fStoreMcIn = value;}
+ virtual void SetAnalysisPbPb(Bool_t isanaPbPb) {fAnalysisPbPb = isanaPbPb;}
+
+ private:
+ virtual Float_t GetVertex(const AliVEvent* event) const;
+ virtual void AnalyzeESD(AliESDEvent* esd);
+ virtual void AnalyzeAOD(AliAODEvent* aod);
+ virtual void ProduceArrayTrksESD(AliESDEvent* event);
+ virtual void ProduceArrayV0ESD(AliESDEvent* event);
+ virtual void ProduceArrayTrksAOD(AliAODEvent* event);
+ virtual void ProduceArrayV0AOD(AliAODEvent* event);
+ Short_t GetPidCode(Int_t pdgCode) const;
+ void ProcessMCTruthESD();
+ void ProcessMCTruthAOD();
+
+ Short_t GetPythiaEventProcessType(Int_t pythiaType);
+ Short_t GetDPMjetEventProcessType(Int_t dpmJetType);
+ ULong64_t GetEventIdAsLong(AliVHeader* header) const;
+
+ TParticle* FindPrimaryMother(AliStack* stack, Int_t label);
+ Int_t FindPrimaryMotherLabel(AliStack* stack, Int_t label);
+
+ AliAODMCParticle* FindPrimaryMotherAOD(AliAODMCParticle* startParticle);
+
+ TParticle* FindPrimaryMotherV0(AliStack* stack, Int_t label);
+ Int_t FindPrimaryMotherLabelV0(AliStack* stack, Int_t label, Int_t& nSteps);
+ Bool_t PhiCut(Double_t pt, Double_t phi, Double_t q, Float_t mag, TF1* phiCutLow, TF1* phiCutHigh);
+
+
+ AliAODMCParticle* FindPrimaryMotherAODV0(AliAODMCParticle* startParticle, Int_t& nSteps);
+
+
+
+ static const Double_t fgkClight; // Speed of light (cm/ps)
+
+ AliESDEvent* fESD; //! ESD object
+ AliAODEvent* fAOD; //! AOD object
+ AliMCEvent* fMC; //! MC object
+ AliStack* fMCStack; //! MC ESD stack
+ TClonesArray* fMCArray; //! MC array for AOD
+ AliAnalysisFilter* fTrackFilterGolden; // Track Filter, set 2010 with golden cuts
+ AliAnalysisFilter* fTrackFilterTPC; // track filter for TPC only tracks
+ TString fCentEst; // V0A , V0M,
+ TString fAnalysisType; // "ESD" or "AOD"
+ Bool_t fAnalysisMC; // Real(kFALSE) or MC(kTRUE) flag
+ Bool_t fAnalysisPbPb; // true you want to analyze PbPb data, false for pp
+ UInt_t ftrigBit;
+ TRandom* fRandom; //! random number generator
+ Bool_t fPileUpRej; // kTRUE is pile-up is rejected
+
+
+
+ //
+ // Cuts and options
+ //
+
+ Double_t fVtxCut; // Vtx cut on z position in cm
+ Double_t fEtaCut; // Eta cut used to select particles
+ Float_t fMinCent; //minimum centrality
+ Float_t fMaxCent; //maximum centrality
+ Bool_t fStoreMcIn; // Store MC input tracks
+ //
+ // Help variables
+ //
+ Short_t fMcProcessType; // -1=invalid, 0=data, 1=ND, 2=SD, 3=DD
+ Short_t fTriggeredEventMB; // 1 = triggered, 0 = not trigged (MC only)
+ Short_t fVtxStatus; // -1 = no vtx, 0 = outside cut, 1 = inside cut
+ Float_t fZvtx; // z vertex
+ Float_t fZvtxMC; // z vertex MC (truth)
+ Int_t fRun; // run no
+ ULong64_t fEventId; // unique event id
+
+ //
+ // Output objects
+ //
+ TList* fListOfObjects; //! Output list of objects
+ TH1I* fEvents; //! No of accepted events
+ TH1I* fVtx; //! Event vertex info
+ TH1F* fVtxMC; //! Event vertex info for ALL MC events
+ TH1F* fVtxBeforeCuts; //! Vertex z dist before cuts
+ TH1F* fVtxAfterCuts; //! Vertex z dist after cuts
+ TH1F* fn1;
+ TH1F* fcent;
+
+ TH2D *hMIPVsEta;
+ TProfile *pMIPVsEta;
+ TH2D *hMIPVsEtaV0s;
+ TProfile *pMIPVsEtaV0s;
+ TH2D *hPlateauVsEta;
+ TProfile *pPlateauVsEta;
+ TH2D *hPhi;
+
+ TH2D *hMIPVsNch[9];
+ TProfile *pMIPVsNch[9];
+
+ TH2D *hMIPVsPhi[9];
+ TProfile *pMIPVsPhi[9];
+ TH2D *hPlateauVsPhi[9];
+ TProfile *pPlateauVsPhi[9];
+
+ TH2D* histPiV0[9];
+ TH2D* histPV0[9];
+
+ TH2D* histAllCh[9];
+ TH2D* histPiTof[9];
+ TH2D* histEV0[9];
+ TH1D* hMcIn[7][9];
+ TH1D* hMcOut[7][9];
+
+
+
+ //TTree* fTree; //! Debug tree
+
+ ClassDef(AliAnalysisTaskQAHighPtDeDx, 1); //Analysis task for high pt analysis
+};
+
+#endif
--- /dev/null
+/**************************************************************************
+ * Authors : Antonin Maire, Boris Hippolyte *
+ * 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. *
+ **************************************************************************/
+
+//-----------------------------------------------------------------
+// AliQAProdMultistrange class
+//
+// Origin AliAnalysisTaskCheckCascade which has four roles :
+// 1. QAing the Cascades from ESD and AOD
+// Origin: AliAnalysisTaskESDCheckV0 by Boris Hippolyte Nov2007, hippolyt@in2p3.fr
+// 2. Prepare the plots which stand as raw material for yield extraction (wi/wo PID)
+// 3. Supply an AliCFContainer meant to define the optimised topological selections
+// 4. Rough azimuthal correlation study (Eta, Phi)
+// Adapted to Cascade : A.Maire Mar2008, antonin.maire@ires.in2p3.fr
+// Modified : A.Maire Mar2010
+//
+// Adapted to PbPb analysis: M. Nicassio, maria.nicassio@ba.infn.it
+// Feb-August2011
+// - Physics selection moved to the run.C macro
+// - Centrality selection added (+ setters)
+// - setters added (vertex range)
+// - histo added and histo/container binning changed
+// - protection in the destructor for CAF usage
+// - AliWarning disabled
+// - automatic settings for PID
+// September2011
+// - proper time histos/container added (V0 and Cascades)
+// November2011
+// - re-run V0's and cascade's vertexers (SetCuts instead of SetDefaultCuts!!)
+// Genuary2012
+// - AOD analysis part completed
+// March2012
+// - min number of TPC clusters for track selection as a parameter
+// July2012
+// - cut on min pt for daughter tracks as a parameter (+control histos)
+// August2012
+// - cut on pseudorapidity for daughter tracks as a parameter (+control histos for Xi-)
+//-----------------------------------------------------------------
+
+class TTree;
+class TParticle;
+class TVector3;
+
+class AliESDVertex;
+class AliAODVertex;
+class AliESDv0;
+class AliAODv0;
+
+#include <Riostream.h>
+#include "THnSparse.h"
+#include "TVector3.h"
+#include "TMath.h"
+
+#include "AliLog.h"
+#include "AliCentrality.h"
+#include "AliESDEvent.h"
+#include "AliAODEvent.h"
+#include "AliESDtrackCuts.h"
+#include "AliPIDResponse.h"
+
+#include "AliInputEventHandler.h"
+#include "AliAnalysisManager.h"
+#include "AliESDInputHandler.h"
+#include "AliAODInputHandler.h"
+#include "AliCFContainer.h"
+#include "AliMultiplicity.h"
+
+#include "AliESDcascade.h"
+#include "AliAODcascade.h"
+#include "AliAODTrack.h"
+
+#include "AliQAProdMultistrange.h"
+
+ClassImp(AliQAProdMultistrange)
+
+
+
+//________________________________________________________________________
+AliQAProdMultistrange::AliQAProdMultistrange()
+ : AliAnalysisTaskSE(),
+ fAnalysisType ("ESD"),
+ fESDtrackCuts (0),
+ fCollidingSystem ("PbPb"),
+ fPIDResponse (0),
+ fkSDDSelectionOn (kTRUE),
+ fkQualityCutZprimVtxPos (kTRUE),
+ fkQualityCutNoTPConlyPrimVtx(kTRUE),
+ fkQualityCutTPCrefit (kTRUE),
+ fkQualityCutnTPCcls (kTRUE),
+ fkQualityCutPileup (kTRUE),
+ fwithSDD (kTRUE),
+ fMinnTPCcls (0),
+ fCentrLowLim (0),
+ fCentrUpLim (0),
+ fCentrEstimator (0),
+ fkUseCleaning (0),
+ fVtxRange (0),
+ fMinPtCutOnDaughterTracks (0),
+ fEtaCutOnDaughterTracks (0),
+
+
+ fCFContCascadeCuts(0)
+
+
+{
+ // Dummy Constructor
+}
+
+
+//________________________________________________________________________
+AliQAProdMultistrange::AliQAProdMultistrange(const char *name)
+ : AliAnalysisTaskSE(name),
+ fAnalysisType ("ESD"),
+ fESDtrackCuts (0),
+ fCollidingSystem ("PbPb"),
+ fPIDResponse (0),
+ fkSDDSelectionOn (kTRUE),
+ fkQualityCutZprimVtxPos (kTRUE),
+ fkQualityCutNoTPConlyPrimVtx(kTRUE),
+ fkQualityCutTPCrefit (kTRUE),
+ fkQualityCutnTPCcls (kTRUE),
+ fkQualityCutPileup (kTRUE),
+ fwithSDD (kTRUE),
+ fMinnTPCcls (0),
+ fCentrLowLim (0),
+ fCentrUpLim (0),
+ fCentrEstimator (0),
+ fkUseCleaning (0),
+ fVtxRange (0),
+ fMinPtCutOnDaughterTracks (0),
+ fEtaCutOnDaughterTracks (0),
+
+
+ fCFContCascadeCuts(0)
+
+
+{
+ // Constructor
+ // Output slot #0 writes into a TList container (Cascade)
+ DefineOutput(1, AliCFContainer::Class());
+
+ AliLog::SetClassDebugLevel("AliQAProdMultistrange",1); // MN this should (?) enable only AliFatal
+}
+
+
+AliQAProdMultistrange::~AliQAProdMultistrange() {
+ //
+ // Destructor
+ //
+ // For all TH1, 2, 3 HnSparse and CFContainer are in the fListCascade TList.
+ // They will be deleted when fListCascade is deleted by the TSelector dtor
+ // Because of TList::SetOwner() ...
+ if (fCFContCascadeCuts && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { delete fCFContCascadeCuts; fCFContCascadeCuts = 0x0; }
+ if (fESDtrackCuts) { delete fESDtrackCuts; fESDtrackCuts = 0x0; }
+}
+
+
+
+//________________________________________________________________________
+void AliQAProdMultistrange::UserCreateOutputObjects() {
+ // Create histograms
+ // Called once
+
+ //-----------------------------------------------
+ // Particle Identification Setup (new PID object)
+ //-----------------------------------------------
+ AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
+ AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
+ fPIDResponse = inputHandler->GetPIDResponse();
+
+
+ // Only used to get the number of primary reconstructed tracks
+ if (fAnalysisType == "ESD" && (! fESDtrackCuts )){
+ fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
+ //Printf("CheckCascade - ESDtrackCuts set up to 2010 std ITS-TPC cuts...");
+ }
+
+
+ //---------------------------------------------------
+ // Define the container for the topological variables
+ //---------------------------------------------------
+ if(! fCFContCascadeCuts) {
+ // Container meant to store all the relevant distributions corresponding to the cut variables.
+ // NB: overflow/underflow of variables on which we want to cut later should be 0!!!
+ const Int_t lNbSteps = 4 ;
+ const Int_t lNbVariables = 21 ;
+ //Array for the number of bins in each dimension :
+ Int_t lNbBinsPerVar[lNbVariables] = {0};
+ lNbBinsPerVar[0] = 25; //DcaCascDaughters : [0.0,2.4,3.0] -> Rec.Cut = 2.0;
+ lNbBinsPerVar[1] = 25; //DcaBachToPrimVertex : [0.0,0.24,100.0] -> Rec.Cut = 0.01;
+ lNbBinsPerVar[2] = 30; //CascCosineOfPointingAngle : [0.97,1.0] -> Rec.Cut = 0.98;
+ lNbBinsPerVar[3] = 40; //CascRadius : [0.0,3.9,1000.0] -> Rec.Cut = 0.2;
+ lNbBinsPerVar[4] = 30; //InvMassLambdaAsCascDghter : [1.1,1.3] -> Rec.Cut = 0.008;
+ lNbBinsPerVar[5] = 20; //DcaV0Daughters : [0.0,2.0] -> Rec.Cut = 1.5;
+ lNbBinsPerVar[6] = 201; //V0CosineOfPointingAngleToPV : [0.89,1.0] -> Rec.Cut = 0.9;
+ lNbBinsPerVar[7] = 40; //V0Radius : [0.0,3.9,1000.0] -> Rec.Cut = 0.2;
+ lNbBinsPerVar[8] = 40; //DcaV0ToPrimVertex : [0.0,0.39,110.0] -> Rec.Cut = 0.01;
+ lNbBinsPerVar[9] = 25; //DcaPosToPrimVertex : [0.0,0.24,100.0] -> Rec.Cut = 0.05;
+ lNbBinsPerVar[10] = 25; //DcaNegToPrimVertex : [0.0,0.24,100.0] -> Rec.Cut = 0.05
+ lNbBinsPerVar[11] = 150; //InvMassXi : 2-MeV/c2 bins
+ lNbBinsPerVar[12] = 120; //InvMassOmega : 2-MeV/c2 bins
+ lNbBinsPerVar[13] = 100; //XiTransvMom : [0.0,10.0]
+ lNbBinsPerVar[14] = 110; //Y(Xi) : 0.02 in rapidity units
+ lNbBinsPerVar[15] = 110; //Y(Omega) : 0.02 in rapidity units
+ lNbBinsPerVar[16] = 112; //Proper lenght of cascade
+ lNbBinsPerVar[17] = 112; //Proper lenght of V0
+ lNbBinsPerVar[18] = 201; //V0CosineOfPointingAngleToXiV
+ lNbBinsPerVar[19] = 11; //Centrality
+ lNbBinsPerVar[20] = 100; //ESD track multiplicity
+ //define the container
+ fCFContCascadeCuts = new AliCFContainer("fCFContCascadeCuts","Container for Cascade cuts", lNbSteps, lNbVariables, lNbBinsPerVar );
+ //Setting the bin limits
+ //0 - DcaXiDaughters
+ Double_t *lBinLim0 = new Double_t[ lNbBinsPerVar[0] + 1 ];
+ for(Int_t i=0; i< lNbBinsPerVar[0]; i++) lBinLim0[i] = (Double_t)0.0 + (2.4 - 0.0)/(lNbBinsPerVar[0] - 1) * (Double_t)i;
+ lBinLim0[ lNbBinsPerVar[0] ] = 3.0;
+ fCFContCascadeCuts -> SetBinLimits(0, lBinLim0);
+ delete [] lBinLim0;
+ //1 - DcaToPrimVertexXi
+ Double_t *lBinLim1 = new Double_t[ lNbBinsPerVar[1] + 1 ];
+ for(Int_t i=0; i<lNbBinsPerVar[1]; i++) lBinLim1[i] = (Double_t)0.0 + (0.24 - 0.0)/(lNbBinsPerVar[1] - 1) * (Double_t)i;
+ lBinLim1[ lNbBinsPerVar[1] ] = 100.0;
+ fCFContCascadeCuts -> SetBinLimits(1, lBinLim1);
+ delete [] lBinLim1;
+ //2 - CascCosineOfPointingAngle
+ fCFContCascadeCuts->SetBinLimits(2, 0.97, 1.);
+ //3 - CascRadius
+ Double_t *lBinLim3 = new Double_t[ lNbBinsPerVar[3]+1 ];
+ for(Int_t i=0; i< lNbBinsPerVar[3]; i++) lBinLim3[i] = (Double_t)0.0 + (3.9 - 0.0 )/(lNbBinsPerVar[3] - 1) * (Double_t)i ;
+ lBinLim3[ lNbBinsPerVar[3] ] = 1000.0;
+ fCFContCascadeCuts -> SetBinLimits(3, lBinLim3 );
+ delete [] lBinLim3;
+ //4 - InvMassLambdaAsCascDghter
+ fCFContCascadeCuts->SetBinLimits(4, 1.1, 1.13);
+ //5 - DcaV0Daughters
+ fCFContCascadeCuts -> SetBinLimits(5, 0., 2.);
+ //6 - V0CosineOfPointingAngleToPV
+ fCFContCascadeCuts -> SetBinLimits(6, 0.8, 1.001);
+ //7 - V0Radius
+ Double_t *lBinLim7 = new Double_t[ lNbBinsPerVar[7] + 1];
+ for(Int_t i=0; i< lNbBinsPerVar[7];i++) lBinLim7[i] = (Double_t)0.0 + (3.9 - 0.0)/(lNbBinsPerVar[7] - 1) * (Double_t)i;
+ lBinLim7[ lNbBinsPerVar[7] ] = 1000.0;
+ fCFContCascadeCuts -> SetBinLimits(7, lBinLim7);
+ delete [] lBinLim7;
+ //8 - DcaV0ToPrimVertex
+ Double_t *lBinLim8 = new Double_t[ lNbBinsPerVar[8]+1 ];
+ for(Int_t i=0; i< lNbBinsPerVar[8];i++) lBinLim8[i] = (Double_t)0.0 + (0.39 - 0.0 )/(lNbBinsPerVar[8]-1) * (Double_t)i ;
+ lBinLim8[ lNbBinsPerVar[8] ] = 100.0;
+ fCFContCascadeCuts -> SetBinLimits(8, lBinLim8 );
+ delete [] lBinLim8;
+ //9 - DcaPosToPrimVertex
+ Double_t *lBinLim9 = new Double_t[ lNbBinsPerVar[9]+1 ];
+ for(Int_t i=0; i< lNbBinsPerVar[9];i++) lBinLim9[i] = (Double_t)0.0 + (0.24 - 0.0 )/(lNbBinsPerVar[9]-1) * (Double_t)i ;
+ lBinLim9[ lNbBinsPerVar[9] ] = 100.0;
+ fCFContCascadeCuts -> SetBinLimits(9, lBinLim9 );
+ delete [] lBinLim9;
+ //10 - DcaNegToPrimVertex
+ Double_t *lBinLim10 = new Double_t[ lNbBinsPerVar[10]+1 ];
+ for(Int_t i=0; i< lNbBinsPerVar[10];i++) lBinLim10[i] = (Double_t)0.0 + (0.24 - 0.0 )/(lNbBinsPerVar[10]-1) * (Double_t)i ;
+ lBinLim10[ lNbBinsPerVar[10] ] = 100.0;
+ fCFContCascadeCuts -> SetBinLimits(10, lBinLim10 );
+ delete [] lBinLim10;
+ //11 - InvMassXi
+ fCFContCascadeCuts->SetBinLimits(11, 1.25, 1.40);
+ //12 - InvMassOmega
+ fCFContCascadeCuts->SetBinLimits(12, 1.62, 1.74);
+ //13 - XiTransvMom
+ fCFContCascadeCuts->SetBinLimits(13, 0.0, 10.0);
+ //14 - Y(Xi)
+ fCFContCascadeCuts->SetBinLimits(14, -1.1, 1.1);
+ //15 - Y(Omega)
+ fCFContCascadeCuts->SetBinLimits(15, -1.1, 1.1);
+ //16 - Proper time of cascade
+ Double_t *lBinLim16 = new Double_t[ lNbBinsPerVar[16]+1 ];
+ for(Int_t i=0; i< lNbBinsPerVar[16];i++) lBinLim16[i] = (Double_t) -1. + (110. + 1.0 ) / (lNbBinsPerVar[16] - 1) * (Double_t) i;
+ lBinLim16[ lNbBinsPerVar[16] ] = 2000.0;
+ fCFContCascadeCuts->SetBinLimits(16, lBinLim16);
+ //17 - Proper time of V0
+ fCFContCascadeCuts->SetBinLimits(17, lBinLim16);
+ //18 - V0CosineOfPointingAngleToXiV
+ fCFContCascadeCuts -> SetBinLimits(18, 0.8, 1.001);
+ //19
+ Double_t *lBinLim19 = new Double_t[ lNbBinsPerVar[19]+1 ];
+ for(Int_t i=3; i< lNbBinsPerVar[19]+1;i++) lBinLim19[i] = (Double_t)(i-1)*10.;
+ lBinLim19[0] = 0.0;
+ lBinLim19[1] = 5.0;
+ lBinLim19[2] = 10.0;
+ fCFContCascadeCuts->SetBinLimits(19, lBinLim19 );
+ delete [] lBinLim19;
+ //20
+ fCFContCascadeCuts->SetBinLimits(20, 0.0, 6000.0);
+ // Setting the number of steps : one for each cascade species (Xi-, Xi+ and Omega-, Omega+)
+ fCFContCascadeCuts->SetStepTitle(0, "#Xi^{-} candidates");
+ fCFContCascadeCuts->SetStepTitle(1, "#bar{#Xi}^{+} candidates");
+ fCFContCascadeCuts->SetStepTitle(2, "#Omega^{-} candidates");
+ fCFContCascadeCuts->SetStepTitle(3, "#bar{#Omega}^{+} candidates");
+ // Setting the variable title, per axis
+ fCFContCascadeCuts->SetVarTitle(0, "Dca(cascade daughters) (cm)");
+ fCFContCascadeCuts->SetVarTitle(1, "ImpactParamToPV(bachelor) (cm)");
+ fCFContCascadeCuts->SetVarTitle(2, "cos(cascade PA)");
+ fCFContCascadeCuts->SetVarTitle(3, "R_{2d}(cascade decay) (cm)");
+ fCFContCascadeCuts->SetVarTitle(4, "M_{#Lambda}(as casc dghter) (GeV/c^{2})");
+ fCFContCascadeCuts->SetVarTitle(5, "Dca(V0 daughters) in Xi (cm)");
+ fCFContCascadeCuts->SetVarTitle(6, "cos(V0 PA) in cascade to PV");
+ fCFContCascadeCuts->SetVarTitle(7, "R_{2d}(V0 decay) (cm)");
+ fCFContCascadeCuts->SetVarTitle(8, "ImpactParamToPV(V0) (cm)");
+ fCFContCascadeCuts->SetVarTitle(9, "ImpactParamToPV(Pos) (cm)");
+ fCFContCascadeCuts->SetVarTitle(10, "ImpactParamToPV(Neg) (cm)");
+ fCFContCascadeCuts->SetVarTitle(11, "Inv. Mass(Xi) (GeV/c^{2})");
+ fCFContCascadeCuts->SetVarTitle(12, "Inv. Mass(Omega) (GeV/c^{2})");
+ fCFContCascadeCuts->SetVarTitle(13, "pt(cascade) (GeV/c)");
+ fCFContCascadeCuts->SetVarTitle(14, "Y(Xi)");
+ fCFContCascadeCuts->SetVarTitle(15, "Y(Omega)");
+ fCFContCascadeCuts->SetVarTitle(16, "mL/p (cascade) (cm)");
+ fCFContCascadeCuts->SetVarTitle(17, "mL/p (V0) (cm)");
+ fCFContCascadeCuts->SetVarTitle(18, "cos(V0 PA) in cascade to XiV");
+ fCFContCascadeCuts->SetVarTitle(19, "Centrality");
+ fCFContCascadeCuts->SetVarTitle(20, "ESD track multiplicity");
+ }
+
+PostData(1, fCFContCascadeCuts);
+
+}// end UserCreateOutputObjects
+
+
+//________________________________________________________________________
+void AliQAProdMultistrange::UserExec(Option_t *) {
+
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ // Main loop (called for each event)
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ //-----------------------
+ //Define ESD/AOD handlers
+ AliESDEvent *lESDevent = 0x0;
+ AliAODEvent *lAODevent = 0x0;
+
+ //---------------------
+ //Check the PIDresponse
+ if(!fPIDResponse) {
+ AliError("Cannot get pid response");
+ return;
+ }
+
+ //__________________________________________________
+ // After these lines we should have an ESD/AOD event
+
+ //---------------------------------------------------------
+ //Load the InputEvent and check, before any event selection
+ //---------------------------------------------------------
+ Float_t lPrimaryTrackMultiplicity = -1.;
+ AliCentrality* centrality;
+ if (fAnalysisType == "ESD") {
+ lESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
+ if (!lESDevent) {
+ AliWarning("ERROR: lESDevent not available \n");
+ return;
+ }
+ if (fCollidingSystem == "PbPb") lPrimaryTrackMultiplicity = fESDtrackCuts->CountAcceptedTracks(lESDevent);
+ if (fCollidingSystem == "PbPb") centrality = lESDevent->GetCentrality();
+
+ } else if (fAnalysisType == "AOD") {
+ lAODevent = dynamic_cast<AliAODEvent*>( InputEvent() );
+ if (!lAODevent) {
+ AliWarning("ERROR: lAODevent not available \n");
+ return;
+ }
+ if (fCollidingSystem == "PbPb") {
+ lPrimaryTrackMultiplicity = 0;
+ Int_t nTrackMultiplicity = (InputEvent())->GetNumberOfTracks();
+ for (Int_t itrack = 0; itrack < nTrackMultiplicity; itrack++) {
+ AliAODTrack* track = lAODevent->GetTrack(itrack);
+ if (track->TestFilterBit(AliAODTrack::kTrkGlobalNoDCA)) lPrimaryTrackMultiplicity++;
+ }
+ }
+ if (fCollidingSystem == "PbPb") centrality = lAODevent->GetCentrality();
+ } else {
+ Printf("Analysis type (ESD or AOD) not specified \n");
+ return;
+ }
+
+ //-----------------------------------------
+ // Centrality selection for PbPb collisions
+ //-----------------------------------------
+ Float_t lcentrality = 0.;
+ if (fCollidingSystem == "PbPb") {
+ if (fkUseCleaning) lcentrality = centrality->GetCentralityPercentile(fCentrEstimator.Data());
+ else {
+ lcentrality = centrality->GetCentralityPercentileUnchecked(fCentrEstimator.Data());
+ if (centrality->GetQuality()>1) {
+ PostData(1, fCFContCascadeCuts);
+ return;
+ }
+ }
+ if (lcentrality<fCentrLowLim||lcentrality>=fCentrUpLim) {
+ PostData(1, fCFContCascadeCuts);
+ return;
+ }
+ } else if (fCollidingSystem == "pp") lcentrality = 0.;
+
+
+ //----------------------------------------
+ // SDD selection for pp@2.76TeV collisions
+ //----------------------------------------
+ if (fCollidingSystem == "pp") {
+ if (fAnalysisType == "ESD") {
+ if (fkSDDSelectionOn) {
+ TString trcl = lESDevent->GetFiredTriggerClasses();
+ if (fwithSDD) { if(!(trcl.Contains("ALLNOTRD"))) { PostData(1, fCFContCascadeCuts); return; } }
+ else if (!fwithSDD){ if((trcl.Contains("ALLNOTRD"))) { PostData(1, fCFContCascadeCuts); return; } }
+ }
+ } else if (fAnalysisType == "AOD") {
+ if (fkSDDSelectionOn) {
+ TString trcl = lAODevent->GetFiredTriggerClasses();
+ if (fwithSDD) { if(!(trcl.Contains("ALLNOTRD"))) { PostData(1, fCFContCascadeCuts); return; } }
+ else if (!fwithSDD) { if((trcl.Contains("ALLNOTRD"))) { PostData(1, fCFContCascadeCuts); return; } }
+ }
+ }
+ }
+
+ //--------------------------------------------
+ // Physics selection for pp@2.76TeV collisions
+ //--------------------------------------------
+ // - moved to the runGrid for the PbPb collisions
+ if (fCollidingSystem == "pp") {
+ if (fAnalysisType == "ESD") {
+ UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
+ Bool_t isSelected = 0;
+ isSelected = (maskIsSelected & AliVEvent::kMB) == AliVEvent::kMB;
+ if(! isSelected){ PostData(1, fCFContCascadeCuts); return; }
+ } else if (fAnalysisType == "AOD") {
+ UInt_t maskIsSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
+ Bool_t isSelected = 0;
+ isSelected = (maskIsSelected & AliVEvent::kMB) == AliVEvent::kMB;
+ if(! isSelected){ PostData(1, fCFContCascadeCuts); return; }
+ }
+ }
+
+ //------------------------------
+ // Well-established PV selection
+ //------------------------------
+ Double_t lBestPrimaryVtxPos[3] = {-100.0, -100.0, -100.0};
+ Double_t lMagneticField = -10.;
+ if (fAnalysisType == "ESD") {
+ const AliESDVertex *lPrimaryTrackingESDVtx = lESDevent->GetPrimaryVertexTracks();
+ const AliESDVertex *lPrimaryBestESDVtx = lESDevent->GetPrimaryVertex();
+ if (fkQualityCutNoTPConlyPrimVtx) {
+ const AliESDVertex *lPrimarySPDVtx = lESDevent->GetPrimaryVertexSPD();
+ if (!lPrimarySPDVtx->GetStatus() && !lPrimaryTrackingESDVtx->GetStatus() ){
+ AliWarning(" No SPD prim. vertex nor prim. Tracking vertex ... return !");
+ PostData(1, fCFContCascadeCuts);
+ return;
+ }
+ }
+ lPrimaryBestESDVtx->GetXYZ( lBestPrimaryVtxPos );
+ lMagneticField = lESDevent->GetMagneticField( );
+ } else if (fAnalysisType == "AOD") {
+ const AliAODVertex *lPrimaryBestAODVtx = lAODevent->GetPrimaryVertex();
+ if (!lPrimaryBestAODVtx){
+ AliWarning("No prim. vertex in AOD... return!");
+ PostData(1, fCFContCascadeCuts);
+ return;
+ }
+ lPrimaryBestAODVtx->GetXYZ( lBestPrimaryVtxPos );
+ lMagneticField = lAODevent->GetMagneticField();
+ }
+
+ //------------------------------------------
+ // Pilup selection for pp@2.76TeV collisions
+ //------------------------------------------
+ if (fCollidingSystem == "pp") {
+ if (fAnalysisType == "ESD") {
+ if (fkQualityCutPileup) { if(lESDevent->IsPileupFromSPD()){ PostData(1, fCFContCascadeCuts); return; } }
+ } else if (fAnalysisType == "AOD") {
+ if (fkQualityCutPileup) { if(lAODevent->IsPileupFromSPD()){ PostData(1, fCFContCascadeCuts); return; } }
+ }
+ }
+
+ //----------------------------
+ // Vertex Z position selection
+ //----------------------------
+ if (fkQualityCutZprimVtxPos) {
+ if (TMath::Abs(lBestPrimaryVtxPos[2]) > fVtxRange ) {
+ PostData(1, fCFContCascadeCuts);
+ return;
+ }
+ }
+
+
+
+ //////////////////////////////
+ // CASCADE RECONSTRUCTION PART
+ //////////////////////////////
+
+ //%%%%%%%%%%%%%
+ // Cascade loop
+ Int_t ncascades = 0;
+ if (fAnalysisType == "ESD") ncascades = lESDevent->GetNumberOfCascades();
+ else if (fAnalysisType == "AOD") ncascades = lAODevent->GetNumberOfCascades();
+
+ for (Int_t iXi = 0; iXi < ncascades; iXi++) {// This is the begining of the Cascade loop (ESD or AOD)
+
+ // -------------------------------------
+ // - Initialisation of the local variables that will be needed for ESD/AOD
+ // -- Container variables (1st round)
+ Double_t lDcaXiDaughters = -1. ; //[Container]
+ Double_t lXiCosineOfPointingAngle = -1. ; //[Container]
+ Double_t lPosXi[3] = { -1000.0, -1000.0, -1000.0 }; //Useful to define other variables: radius fid. vol., ctau, etc. for cascade
+ Double_t lXiRadius = -1000. ; //[Container]
+ UShort_t lPosTPCClusters = -1; //To check the quality of the tracks. For ESD only ...
+ UShort_t lNegTPCClusters = -1; //To check the quality of the tracks. For ESD only ...
+ UShort_t lBachTPCClusters = -1; //To check the quality of the tracks. For ESD only ...
+ Double_t lInvMassLambdaAsCascDghter = 0.; //[Container]
+ Double_t lDcaV0DaughtersXi = -1.; //[Container]
+ Double_t lDcaBachToPrimVertexXi = -1.; //[Container]
+ Double_t lDcaV0ToPrimVertexXi = -1.; //[Container]
+ Double_t lDcaPosToPrimVertexXi = -1.; //[Container]
+ Double_t lDcaNegToPrimVertexXi = -1.; //[Container]
+ Double_t lV0CosineOfPointingAngle = -1.; //[Container]
+ Double_t lV0toXiCosineOfPointingAngle = -1.; //[Container]
+ Double_t lPosV0Xi[3] = { -1000. , -1000., -1000. }; //Useful to define other variables: radius fid. vol., ctau, etc. for VO
+ Double_t lV0RadiusXi = -1000.0; //[Container]
+ Double_t lV0quality = 0.; // ??
+ Double_t lInvMassXiMinus = 0.; //[Container]
+ Double_t lInvMassXiPlus = 0.; //[Container]
+ Double_t lInvMassOmegaMinus = 0.; //[Container]
+ Double_t lInvMassOmegaPlus = 0.; //[Container]
+ // -- PID treatment
+ Bool_t lIsBachelorKaonForTPC = kFALSE;
+ Bool_t lIsBachelorPionForTPC = kFALSE;
+ Bool_t lIsNegPionForTPC = kFALSE;
+ Bool_t lIsPosPionForTPC = kFALSE;
+ Bool_t lIsNegProtonForTPC = kFALSE;
+ Bool_t lIsPosProtonForTPC = kFALSE;
+ // -- More container variables and quality checks
+ Double_t lXiMomX = 0.; //Useful to define other variables: lXiTransvMom, lXiTotMom
+ Double_t lXiMomY = 0.; //Useful to define other variables: lXiTransvMom, lXiTotMom
+ Double_t lXiMomZ = 0.; //Useful to define other variables: lXiTransvMom, lXiTotMom
+ Double_t lXiTransvMom = 0.; //[Container]
+ Double_t lXiTotMom = 0.; //Useful to define other variables: cTau
+ Double_t lV0PMomX = 0.; //Useful to define other variables: lV0TotMom, lpTrackTransvMom
+ Double_t lV0PMomY = 0.; //Useful to define other variables: lV0TotMom, lpTrackTransvMom
+ Double_t lV0PMomZ = 0.; //Useful to define other variables: lV0TotMom, lpTrackTransvMom
+ Double_t lV0NMomX = 0.; //Useful to define other variables: lV0TotMom, lnTrackTransvMom
+ Double_t lV0NMomY = 0.; //Useful to define other variables: lV0TotMom, lnTrackTransvMom
+ Double_t lV0NMomZ = 0.; //Useful to define other variables: lV0TotMom, lnTrackTransvMom
+ Double_t lV0TotMom = 0.; //Useful to define other variables: lctauV0
+ Double_t lBachMomX = 0.; //Useful to define other variables: lBachTransvMom
+ Double_t lBachMomY = 0.; //Useful to define other variables: lBachTransvMom
+ Double_t lBachMomZ = 0.; //Useful to define other variables: lBachTransvMom
+ Double_t lBachTransvMom = 0.; //Selection on the min bachelor pT
+ Double_t lpTrackTransvMom = 0.; //Selection on the min bachelor pT
+ Double_t lnTrackTransvMom = 0.; //Selection on the min bachelor pT
+ Short_t lChargeXi = -2; //Useful to select the particles based on the charge
+ Double_t lRapXi = -20.0; //[Container]
+ Double_t lRapOmega = -20.0; //[Container]
+ Float_t etaBach = 0.; //Selection on the eta range
+ Float_t etaPos = 0.; //Selection on the eta range
+ Float_t etaNeg = 0.; //Selection on the eta range
+ // -- variables for the AliCFContainer dedicated to cascade cut optmisiation: ESD and AOD
+ if (fAnalysisType == "ESD") {
+
+ // -------------------------------------
+ // - Load the cascades from the handler
+ AliESDcascade *xi = lESDevent->GetCascade(iXi);
+ if (!xi) continue;
+
+ // ---------------------------------------------------------------------------
+ // - Assigning the necessary variables for specific AliESDcascade data members
+ lV0quality = 0.;
+ xi->ChangeMassHypothesis(lV0quality , 3312); // default working hypothesis : cascade = Xi- decay
+ lDcaXiDaughters = xi->GetDcaXiDaughters();
+ lXiCosineOfPointingAngle = xi->GetCascadeCosineOfPointingAngle( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1], lBestPrimaryVtxPos[2] );
+ // Take care : the best available vertex should be used (like in AliCascadeVertexer)
+ xi->GetXYZcascade( lPosXi[0], lPosXi[1], lPosXi[2] );
+ lXiRadius = TMath::Sqrt( lPosXi[0]*lPosXi[0] + lPosXi[1]*lPosXi[1] );
+
+ // -------------------------------------------------------------------------------------------------------------------------------
+ // - Around the tracks : Bach + V0 (ESD). Necessary variables for ESDcascade data members coming from the ESDv0 part (inheritance)
+ UInt_t lIdxPosXi = (UInt_t) TMath::Abs( xi->GetPindex() );
+ UInt_t lIdxNegXi = (UInt_t) TMath::Abs( xi->GetNindex() );
+ UInt_t lBachIdx = (UInt_t) TMath::Abs( xi->GetBindex() );
+ // Care track label can be negative in MC production (linked with the track quality)
+ // However = normally, not the case for track index ...
+ // - Rejection of a double use of a daughter track (nothing but just a crosscheck of what is done in the cascade vertexer)
+ if (lBachIdx == lIdxNegXi) continue;
+ if (lBachIdx == lIdxPosXi) continue;
+ // - Get the track for the daughters
+ AliESDtrack *pTrackXi = lESDevent->GetTrack( lIdxPosXi );
+ AliESDtrack *nTrackXi = lESDevent->GetTrack( lIdxNegXi );
+ AliESDtrack *bachTrackXi = lESDevent->GetTrack( lBachIdx );
+ if (!pTrackXi || !nTrackXi || !bachTrackXi ) continue;
+ // - Get the TPCnumber of cluster for the daughters
+ lPosTPCClusters = pTrackXi->GetTPCNcls();
+ lNegTPCClusters = nTrackXi->GetTPCNcls();
+ lBachTPCClusters = bachTrackXi->GetTPCNcls();
+
+ // ------------------------------------
+ // - Rejection of a poor quality tracks
+ if (fkQualityCutTPCrefit) {
+ // 1 - Poor quality related to TPCrefit
+ ULong_t pStatus = pTrackXi->GetStatus();
+ ULong_t nStatus = nTrackXi->GetStatus();
+ ULong_t bachStatus = bachTrackXi->GetStatus();
+ if ((pStatus&AliESDtrack::kTPCrefit) == 0) { AliWarning(" V0 Pos. track has no TPCrefit ... continue!"); continue; }
+ if ((nStatus&AliESDtrack::kTPCrefit) == 0) { AliWarning(" V0 Neg. track has no TPCrefit ... continue!"); continue; }
+ if ((bachStatus&AliESDtrack::kTPCrefit) == 0) { AliWarning(" Bach. track has no TPCrefit ... continue!"); continue; }
+ }
+ if (fkQualityCutnTPCcls) {
+ // 2 - Poor quality related to TPC clusters
+ if (lPosTPCClusters < fMinnTPCcls) { AliWarning(" V0 Pos. track has less than minn TPC clusters ... continue!"); continue; }
+ if (lNegTPCClusters < fMinnTPCcls) { AliWarning(" V0 Neg. track has less than minn TPC clusters ... continue!"); continue; }
+ if (lBachTPCClusters < fMinnTPCcls) { AliWarning(" Bach. track has less than minn TPC clusters ... continue!"); continue; }
+ }
+
+ // ------------------------------
+ etaPos = pTrackXi->Eta();
+ etaNeg = nTrackXi->Eta();
+ etaBach = bachTrackXi->Eta();
+ lInvMassLambdaAsCascDghter = xi->GetEffMass();
+ lDcaV0DaughtersXi = xi->GetDcaV0Daughters();
+ lV0CosineOfPointingAngle = xi->GetV0CosineOfPointingAngle( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1], lBestPrimaryVtxPos[2] );
+ lDcaV0ToPrimVertexXi = xi->GetD( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1], lBestPrimaryVtxPos[2] );
+ lDcaBachToPrimVertexXi = TMath::Abs( bachTrackXi->GetD( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1], lMagneticField) );
+ xi->GetXYZ( lPosV0Xi[0], lPosV0Xi[1], lPosV0Xi[2] );
+ lV0RadiusXi = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0] + lPosV0Xi[1]*lPosV0Xi[1] );
+ lDcaPosToPrimVertexXi = TMath::Abs( pTrackXi->GetD( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1], lMagneticField ) );
+ lDcaNegToPrimVertexXi = TMath::Abs( nTrackXi->GetD( lBestPrimaryVtxPos[0], lBestPrimaryVtxPos[1], lMagneticField ) );
+
+ //----------------------------------------------------------------------------------------------------
+ // - Around effective masses. Change mass hypotheses to cover all the possibilities: Xi-/+, Omega -/+
+ if (bachTrackXi->Charge() < 0) {
+ // Calculate the effective mass of the Xi- candidate. pdg code 3312 = Xi-
+ lV0quality = 0.;
+ xi->ChangeMassHypothesis(lV0quality , 3312);
+ lInvMassXiMinus = xi->GetEffMassXi();
+ // Calculate the effective mass of the Xi- candidate. pdg code 3334 = Omega-
+ lV0quality = 0.;
+ xi->ChangeMassHypothesis(lV0quality , 3334);
+ lInvMassOmegaMinus = xi->GetEffMassXi();
+ // Back to default hyp.
+ lV0quality = 0.;
+ xi->ChangeMassHypothesis(lV0quality , 3312);
+ }// end if negative bachelor
+ if ( bachTrackXi->Charge() > 0 ) {
+ // Calculate the effective mass of the Xi+ candidate. pdg code -3312 = Xi+
+ lV0quality = 0.;
+ xi->ChangeMassHypothesis(lV0quality , -3312);
+ lInvMassXiPlus = xi->GetEffMassXi();
+ // Calculate the effective mass of the Xi+ candidate. pdg code -3334 = Omega+
+ lV0quality = 0.;
+ xi->ChangeMassHypothesis(lV0quality , -3334);
+ lInvMassOmegaPlus = xi->GetEffMassXi();
+ // Back to "default" hyp.
+ lV0quality = 0.;
+ xi->ChangeMassHypothesis(lV0quality , -3312);
+ }// end if positive bachelor
+
+ // ----------------------------------------------
+ // - TPC PID : 3-sigma bands on Bethe-Bloch curve
+ // Bachelor
+ if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 4) lIsBachelorKaonForTPC = kTRUE;
+ if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 4) lIsBachelorPionForTPC = kTRUE;
+ // Negative V0 daughter
+ if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kPion )) < 4) lIsNegPionForTPC = kTRUE;
+ if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kProton )) < 4) lIsNegProtonForTPC = kTRUE;
+ // Positive V0 daughter
+ if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kPion )) < 4) lIsPosPionForTPC = kTRUE;
+ if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 4) lIsPosProtonForTPC = kTRUE;
+
+ // ------------------------------
+ // - Miscellaneous pieces of info that may help regarding data quality assessment.
+ xi->GetPxPyPz( lXiMomX, lXiMomY, lXiMomZ );
+ lXiTransvMom = TMath::Sqrt( lXiMomX*lXiMomX + lXiMomY*lXiMomY );
+ lXiTotMom = TMath::Sqrt( lXiMomX*lXiMomX + lXiMomY*lXiMomY + lXiMomZ*lXiMomZ );
+ xi->GetNPxPyPz(lV0NMomX,lV0NMomY,lV0NMomZ);
+ xi->GetPPxPyPz(lV0PMomX,lV0PMomY,lV0PMomZ);
+ lV0TotMom = TMath::Sqrt(TMath::Power(lV0NMomX+lV0PMomX,2)+TMath::Power(lV0NMomY+lV0PMomY,2)+TMath::Power(lV0NMomZ+lV0PMomZ,2));
+ xi->GetBPxPyPz( lBachMomX, lBachMomY, lBachMomZ );
+ lBachTransvMom = TMath::Sqrt( lBachMomX*lBachMomX + lBachMomY*lBachMomY );
+ lnTrackTransvMom = TMath::Sqrt( lV0NMomX*lV0NMomX + lV0NMomY*lV0NMomY );
+ lpTrackTransvMom = TMath::Sqrt( lV0PMomX*lV0PMomX + lV0PMomY*lV0PMomY );
+ lChargeXi = xi->Charge();
+ lV0toXiCosineOfPointingAngle = xi->GetV0CosineOfPointingAngle( lPosXi[0], lPosXi[1], lPosXi[2] );
+ lRapXi = xi->RapXi();
+ lRapOmega = xi->RapOmega();
+
+ } else if (fAnalysisType == "AOD") {
+
+ // -------------------------------------
+ // - Load the cascades from the handler
+ const AliAODcascade *xi = lAODevent->GetCascade(iXi);
+ if (!xi) continue;
+
+ //----------------------------------------------------------------------------
+ // - Assigning the necessary variables for specific AliESDcascade data members
+ lDcaXiDaughters = xi->DcaXiDaughters();
+ lXiCosineOfPointingAngle = xi->CosPointingAngleXi( lBestPrimaryVtxPos[0],
+ lBestPrimaryVtxPos[1],
+ lBestPrimaryVtxPos[2] );
+ lPosXi[0] = xi->DecayVertexXiX();
+ lPosXi[1] = xi->DecayVertexXiY();
+ lPosXi[2] = xi->DecayVertexXiZ();
+ lXiRadius = TMath::Sqrt( lPosXi[0]*lPosXi[0] + lPosXi[1]*lPosXi[1] );
+
+ //-------------------------------------------------------------------------------------------------------------------------------
+ // - Around the tracks: Bach + V0 (ESD). Necessary variables for ESDcascade data members coming from the ESDv0 part (inheritance)
+ AliAODTrack *pTrackXi = dynamic_cast<AliAODTrack*>( xi->GetDaughter(0) );
+ AliAODTrack *nTrackXi = dynamic_cast<AliAODTrack*>( xi->GetDaughter(1) );
+ AliAODTrack *bachTrackXi = dynamic_cast<AliAODTrack*>( xi->GetDecayVertexXi()->GetDaughter(0) );
+ if (!pTrackXi || !nTrackXi || !bachTrackXi ) continue;
+ UInt_t lIdxPosXi = (UInt_t) TMath::Abs( pTrackXi->GetID() );
+ UInt_t lIdxNegXi = (UInt_t) TMath::Abs( nTrackXi->GetID() );
+ UInt_t lBachIdx = (UInt_t) TMath::Abs( bachTrackXi->GetID() );
+ // Care track label can be negative in MC production (linked with the track quality)
+ // However = normally, not the case for track index ...
+
+ // - Rejection of a double use of a daughter track (nothing but just a crosscheck of what is done in the cascade vertexer)
+ if (lBachIdx == lIdxNegXi) continue;
+ if (lBachIdx == lIdxPosXi) continue;
+ // - Get the TPCnumber of cluster for the daughters
+ lPosTPCClusters = pTrackXi->GetTPCNcls(); // FIXME: Is this ok? or something like in LambdaK0PbPb task AOD?
+ lNegTPCClusters = nTrackXi->GetTPCNcls();
+ lBachTPCClusters = bachTrackXi->GetTPCNcls();
+
+ // ------------------------------------
+ // - Rejection of a poor quality tracks
+ if (fkQualityCutTPCrefit) {
+ // - Poor quality related to TPCrefit
+ if (!(pTrackXi->IsOn(AliAODTrack::kTPCrefit))) { AliWarning(" V0 Pos. track has no TPCrefit ... continue!"); continue; }
+ if (!(nTrackXi->IsOn(AliAODTrack::kTPCrefit))) { AliWarning(" V0 Neg. track has no TPCrefit ... continue!"); continue; }
+ if (!(bachTrackXi->IsOn(AliAODTrack::kTPCrefit))) { AliWarning(" Bach. track has no TPCrefit ... continue!"); continue; }
+ }
+ if (fkQualityCutnTPCcls) {
+ // - Poor quality related to TPC clusters
+ if (lPosTPCClusters < fMinnTPCcls) continue;
+ if (lNegTPCClusters < fMinnTPCcls) continue;
+ if (lBachTPCClusters < fMinnTPCcls) continue;
+ }
+
+ // ------------------------------------------------------------------------------------------------------------------------------
+ // - Around the tracks: Bach + V0 (AOD). Necessary variables for AODcascade data members coming from the AODv0 part (inheritance)
+ etaPos = pTrackXi->Eta();
+ etaNeg = nTrackXi->Eta();
+ etaBach = bachTrackXi->Eta();
+ lChargeXi = xi->ChargeXi();
+ if ( lChargeXi < 0) lInvMassLambdaAsCascDghter = xi->MassLambda();
+ else lInvMassLambdaAsCascDghter = xi->MassAntiLambda();
+ lDcaV0DaughtersXi = xi->DcaV0Daughters();
+ lDcaV0ToPrimVertexXi = xi->DcaV0ToPrimVertex();
+ lDcaBachToPrimVertexXi = xi->DcaBachToPrimVertex();
+ lPosV0Xi[0] = xi->DecayVertexV0X();
+ lPosV0Xi[1] = xi->DecayVertexV0Y();
+ lPosV0Xi[2] = xi->DecayVertexV0Z();
+ lV0RadiusXi = TMath::Sqrt( lPosV0Xi[0]*lPosV0Xi[0] + lPosV0Xi[1]*lPosV0Xi[1] );
+ lV0CosineOfPointingAngle = xi->CosPointingAngle( lBestPrimaryVtxPos );
+ lDcaPosToPrimVertexXi = xi->DcaPosToPrimVertex();
+ lDcaNegToPrimVertexXi = xi->DcaNegToPrimVertex();
+
+ // ---------------------------------------------------------------------------------------------------
+ // - Around effective masses. Change mass hypotheses to cover all the possibilities: Xi-/+, Omega -/+
+ if ( lChargeXi < 0 ) lInvMassXiMinus = xi->MassXi();
+ if ( lChargeXi > 0 ) lInvMassXiPlus = xi->MassXi();
+ if ( lChargeXi < 0 ) lInvMassOmegaMinus = xi->MassOmega();
+ if ( lChargeXi > 0 ) lInvMassOmegaPlus = xi->MassOmega();
+
+ // ----------------------------------------------
+ // - TPC PID : 3-sigma bands on Bethe-Bloch curve
+ // Bachelor
+ if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kKaon)) < 4) lIsBachelorKaonForTPC = kTRUE;
+ if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( bachTrackXi,AliPID::kPion)) < 4) lIsBachelorPionForTPC = kTRUE;
+ // Negative V0 daughter
+ if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kPion )) < 4) lIsNegPionForTPC = kTRUE;
+ if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( nTrackXi,AliPID::kProton )) < 4) lIsNegProtonForTPC = kTRUE;
+ // Positive V0 daughter
+ if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kPion )) < 4) lIsPosPionForTPC = kTRUE;
+ if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC( pTrackXi,AliPID::kProton )) < 4) lIsPosProtonForTPC = kTRUE;
+
+ //---------------------------------
+ // - Extra info for QA (AOD)
+ // Miscellaneous pieces of info that may help regarding data quality assessment.
+ // Cascade transverse and total momentum
+ lXiMomX = xi->MomXiX();
+ lXiMomY = xi->MomXiY();
+ lXiMomZ = xi->MomXiZ();
+ lXiTransvMom = TMath::Sqrt( lXiMomX*lXiMomX + lXiMomY*lXiMomY );
+ lXiTotMom = TMath::Sqrt( lXiMomX*lXiMomX + lXiMomY*lXiMomY + lXiMomZ*lXiMomZ );
+ Double_t lV0MomX = xi->MomV0X();
+ Double_t lV0MomY = xi->MomV0Y();
+ Double_t lV0MomZ = xi->MomV0Z();
+ lV0TotMom = TMath::Sqrt(TMath::Power(lV0MomX,2)+TMath::Power(lV0MomY,2)+TMath::Power(lV0MomZ,2));
+ lBachMomX = xi->MomBachX();
+ lBachMomY = xi->MomBachY();
+ lBachMomZ = xi->MomBachZ();
+ lBachTransvMom = TMath::Sqrt( lBachMomX*lBachMomX + lBachMomY*lBachMomY );
+ lV0NMomX = xi->MomNegX();
+ lV0NMomY = xi->MomNegY();
+ lV0PMomX = xi->MomPosX();
+ lV0PMomY = xi->MomPosY();
+ lnTrackTransvMom = TMath::Sqrt( lV0NMomX*lV0NMomX + lV0NMomY*lV0NMomY );
+ lpTrackTransvMom = TMath::Sqrt( lV0PMomX*lV0PMomX + lV0PMomY*lV0PMomY );
+ lV0toXiCosineOfPointingAngle = xi->CosPointingAngle( xi->GetDecayVertexXi() );
+ lRapXi = xi->RapXi();
+ lRapOmega = xi->RapOmega();
+
+ }// end of AOD treatment
+
+ //---------------------------------------
+ // Cut on pt of the three daughter tracks
+ if (lBachTransvMom<fMinPtCutOnDaughterTracks) continue;
+ if (lpTrackTransvMom<fMinPtCutOnDaughterTracks) continue;
+ if (lnTrackTransvMom<fMinPtCutOnDaughterTracks) continue;
+
+ //---------------------------------------------------
+ // Cut on pseudorapidity of the three daughter tracks
+ if (TMath::Abs(etaBach)>fEtaCutOnDaughterTracks) continue;
+ if (TMath::Abs(etaPos)>fEtaCutOnDaughterTracks) continue;
+ if (TMath::Abs(etaNeg)>fEtaCutOnDaughterTracks) continue;
+
+ //----------------------------------
+ // Calculate proper time for cascade
+ Double_t cascadeMass = 0.;
+ if ( ( (lChargeXi<0) && lIsBachelorPionForTPC && lIsPosProtonForTPC && lIsNegPionForTPC ) ||
+ ( (lChargeXi>0) && lIsBachelorPionForTPC && lIsNegProtonForTPC && lIsPosPionForTPC ) ) cascadeMass = 1.321;
+ if ( ( (lChargeXi<0) && lIsBachelorKaonForTPC && lIsPosProtonForTPC && lIsNegPionForTPC ) ||
+ ( (lChargeXi>0) && lIsBachelorKaonForTPC && lIsNegProtonForTPC && lIsPosPionForTPC ) ) cascadeMass = 1.672;
+ Double_t lctau = TMath::Sqrt(TMath::Power((lPosXi[0]-lBestPrimaryVtxPos[0]),2)+TMath::Power((lPosXi[1]-lBestPrimaryVtxPos[1]),2)+TMath::Power(( lPosXi[2]-lBestPrimaryVtxPos[2]),2));
+ if (lXiTotMom!=0) lctau = lctau*cascadeMass/lXiTotMom;
+ else lctau = -1.;
+ // Calculate proper time for Lambda (reconstructed)
+ Float_t lambdaMass = 1.115683; // PDG mass
+ Float_t distV0Xi = TMath::Sqrt(TMath::Power((lPosV0Xi[0]-lPosXi[0]),2)+TMath::Power((lPosV0Xi[1]-lPosXi[1]),2)+TMath::Power((lPosV0Xi[2]-lPosXi[2]),2));
+ Float_t lctauV0 = -1.;
+ if (lV0TotMom!=0) lctauV0 = distV0Xi*lambdaMass/lV0TotMom;
+
+
+ // Fill the AliCFContainer (optimisation of topological selections)
+ Double_t lContainerCutVars[21] = {0.0};
+ lContainerCutVars[0] = lDcaXiDaughters;
+ lContainerCutVars[1] = lDcaBachToPrimVertexXi;
+ lContainerCutVars[2] = lXiCosineOfPointingAngle;
+ lContainerCutVars[3] = lXiRadius;
+ lContainerCutVars[4] = lInvMassLambdaAsCascDghter;
+ lContainerCutVars[5] = lDcaV0DaughtersXi;
+ lContainerCutVars[6] = lV0CosineOfPointingAngle;
+ lContainerCutVars[7] = lV0RadiusXi;
+ lContainerCutVars[8] = lDcaV0ToPrimVertexXi;
+ lContainerCutVars[9] = lDcaPosToPrimVertexXi;
+ lContainerCutVars[10] = lDcaNegToPrimVertexXi;
+ lContainerCutVars[13] = lXiTransvMom;
+ lContainerCutVars[16] = lctau;
+ lContainerCutVars[17] = lctauV0;
+ lContainerCutVars[18] = lV0toXiCosineOfPointingAngle;
+ lContainerCutVars[19] = lcentrality;
+ lContainerCutVars[20] = lPrimaryTrackMultiplicity;
+ if ( lChargeXi < 0 ) {
+ lContainerCutVars[11] = lInvMassXiMinus;
+ lContainerCutVars[12] = lInvMassOmegaMinus;
+ lContainerCutVars[14] = lRapXi;
+ lContainerCutVars[15] = -1.;
+ if (lIsBachelorPionForTPC && lIsPosProtonForTPC && lIsNegPionForTPC) fCFContCascadeCuts->Fill(lContainerCutVars,0); // for Xi-
+ lContainerCutVars[11] = lInvMassXiMinus;
+ lContainerCutVars[12] = lInvMassOmegaMinus;
+ lContainerCutVars[14] = -1.;
+ lContainerCutVars[15] = lRapOmega;
+ if (lIsBachelorKaonForTPC && lIsPosProtonForTPC && lIsNegPionForTPC) fCFContCascadeCuts->Fill(lContainerCutVars,2); // for Omega-
+ } else {
+ lContainerCutVars[11] = lInvMassXiPlus;
+ lContainerCutVars[12] = lInvMassOmegaPlus;
+ lContainerCutVars[14] = lRapXi;
+ lContainerCutVars[15] = -1.;
+ if (lIsBachelorPionForTPC && lIsNegProtonForTPC && lIsPosPionForTPC) fCFContCascadeCuts->Fill(lContainerCutVars,1); // for Xi+
+ lContainerCutVars[11] = lInvMassXiPlus;
+ lContainerCutVars[12] = lInvMassOmegaPlus;
+ lContainerCutVars[14] = -1.;
+ lContainerCutVars[15] = lRapOmega;
+ if (lIsBachelorKaonForTPC && lIsNegProtonForTPC && lIsPosPionForTPC) fCFContCascadeCuts->Fill(lContainerCutVars,3); // for Omega+
+ }
+
+
+ }// end of the Cascade loop (ESD or AOD)
+
+
+ // Post output data.
+ PostData(1, fCFContCascadeCuts);
+}
+
+//________________________________________________________________________
+void AliQAProdMultistrange::Terminate(Option_t *)
+{
+
+}
--- /dev/null
+#ifndef ALIANALYSISTASKCHECKCASCADEPBPB_H
+#define ALIANALYSISTASKCHECKCASCADEPBPB_H
+
+/* See cxx source for full Copyright notice */
+
+//-----------------------------------------------------------------
+// AliQAProdMultistrange class
+// Origin AliAnalysisTaskCheckCascade
+// This task has four roles :
+// 1. QAing the Cascades from ESD and AOD
+// Origin: AliAnalysisTaskESDCheckV0 by Boris Hippolyte Nov2007, hippolyt@in2p3.fr
+// 2. Prepare the plots which stand as raw material for yield extraction (wi/wo PID)
+// 3. Supply an AliCFContainer meant to define the optimised topological selections
+// 4. Rough azimuthal correlation study (Eta, Phi)
+// Adapted to Cascade : A.Maire Mar2008, antonin.maire@ires.in2p3.fr
+// Modified : A.Maire Mar2010, antonin.maire@ires.in2p3.fr
+// Modified for PbPb analysis: M. Nicassio Feb 2011, maria.nicassio@ba.infn.it
+//-----------------------------------------------------------------
+
+class TList;
+class TH1F;
+class TH2F;
+class TH3F;
+class TVector3;
+class THnSparse;
+
+class AliESDEvent;
+class AliPhysicsSelection;
+class AliCFContainer;
+class AliPIDResponse;
+
+#include "TString.h"
+
+#include "AliAnalysisTaskSE.h"
+
+class AliQAProdMultistrange : public AliAnalysisTaskSE {
+ public:
+ AliQAProdMultistrange();
+ AliQAProdMultistrange(const char *name);
+ virtual ~AliQAProdMultistrange();
+
+ virtual void UserCreateOutputObjects();
+ virtual void UserExec(Option_t *option);
+ virtual void Terminate(Option_t *);
+
+ void SetAnalysisType (const char* analysisType = "ESD" ) { fAnalysisType = analysisType; }
+ void SetCollidingSystem (const char* collidingSystem = "PbPb") { fCollidingSystem = collidingSystem; }
+ void SetSDDselection (Bool_t SDDSelection = kFALSE) { fkSDDSelectionOn = SDDSelection; }
+ void SetQualityCutZprimVtxPos (Bool_t qualityCutZprimVtxPos = kTRUE ) { fkQualityCutZprimVtxPos = qualityCutZprimVtxPos; }
+ void SetQualityCutNoTPConlyPrimVtx (Bool_t qualityCutNoTPConlyPrimVtx = kTRUE ) { fkQualityCutNoTPConlyPrimVtx = qualityCutNoTPConlyPrimVtx; }
+ void SetQualityCutTPCrefit (Bool_t qualityCutTPCrefit = kTRUE ) { fkQualityCutTPCrefit = qualityCutTPCrefit; }
+ void SetQualityCutnTPCcls (Bool_t qualityCutnTPCcls = kTRUE ) { fkQualityCutnTPCcls = qualityCutnTPCcls; }
+ void SetQualityCutMinnTPCcls (Int_t minnTPCcls = 70 ) { fMinnTPCcls = minnTPCcls; }
+ void SetQualityCutPileup (Bool_t qualitycutPileup = kFALSE) { fkQualityCutPileup = qualitycutPileup; }
+ void SetwithSDD (Bool_t withSDD = kTRUE ) { fwithSDD = withSDD; }
+ void SetCentralityLowLim (Float_t centrlowlim = 0. ) { fCentrLowLim = centrlowlim; }
+ void SetCentralityUpLim (Float_t centruplim = 100. ) { fCentrUpLim = centruplim; }
+ void SetCentralityEst (TString centrest = "V0M" ) { fCentrEstimator = centrest; }
+ void SetUseCleaning (Bool_t usecleaning = kTRUE ) { fkUseCleaning = usecleaning; }
+ void SetVertexRange (Float_t vtxrange = 0. ) { fVtxRange = vtxrange; }
+ void SetMinptCutOnDaughterTracks (Float_t minptdaughtrks = 0. ) { fMinPtCutOnDaughterTracks = minptdaughtrks; }
+ void SetEtaCutOnDaughterTracks (Float_t etadaughtrks = 0. ) { fEtaCutOnDaughterTracks = etadaughtrks; }
+
+ private:
+ // Note : In ROOT, "//!" means "do not stream the data from Master node to Worker node" ...
+ // your data member object is created on the worker nodes and streaming is not needed.
+ // http://root.cern.ch/download/doc/11InputOutput.pdf, page 14
+
+
+ TString fAnalysisType; // "ESD" or "AOD" analysis type
+ AliESDtrackCuts *fESDtrackCuts; // ESD track cuts used for primary track definition
+ TString fCollidingSystem; // "PbPb" or "pp" colliding system
+ AliPIDResponse *fPIDResponse; //! PID response object
+ Bool_t fkSDDSelectionOn; // Boolean : kTRUE = apply the selection on SDD status
+ Bool_t fkQualityCutZprimVtxPos; // Boolean : kTRUE = cut on the prim.vtx z-position
+ Bool_t fkQualityCutNoTPConlyPrimVtx; // Boolean : kTRUE = prim vtx should be SPD or Tracking vertex
+ Bool_t fkQualityCutTPCrefit; // Boolean : kTRUE = ask for TPCrefit for the 3 daughter tracks
+ Bool_t fkQualityCutnTPCcls; // Boolean : kTRUE = ask for at least n TPC clusters for each daughter track
+ Bool_t fkQualityCutPileup; // Boolean : kTRUE = ask for no pileup events
+ Bool_t fwithSDD; // Boolean : kTRUE = Select the events that has and use the info from the SDD
+ Int_t fMinnTPCcls; // minimum number of TPC cluster for daughter tracks
+ Float_t fCentrLowLim; // Lower limit for centrality percentile selection
+ Float_t fCentrUpLim; // Upper limit for centrality percentile selection
+ TString fCentrEstimator; // string for the centrality estimator
+ Bool_t fkUseCleaning; // Boolean : kTRUE = uses all the cleaning criteria of centrality selections (vertex cut + outliers) otherwise only outliers
+ Float_t fVtxRange; // to select events with |zvtx|<fVtxRange cm
+ Float_t fMinPtCutOnDaughterTracks; // minimum pt cut on daughter tracks
+ Float_t fEtaCutOnDaughterTracks; // pseudorapidity cut on daughter tracks
+
+
+
+ AliCFContainer *fCFContCascadeCuts; //! Container meant to store all the relevant distributions corresponding to the cut variables
+
+
+
+ AliQAProdMultistrange(const AliQAProdMultistrange&); // not implemented
+ AliQAProdMultistrange& operator=(const AliQAProdMultistrange&); // not implemented
+
+ ClassDef(AliQAProdMultistrange, 7);
+};
+
+#endif
git://git.uio.no / u/mrichter/AliRoot.git / commitdiff