- /*************************************************************************\r
-* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *\r
-* *\r
-* Author: The ALICE Off-line Project. *\r
-* Contributors are mentioned in the code where appropriate. *\r
-* *\r
-* Permission to use, copy, modify and distribute this software and its *\r
-* documentation strictly for non-commercial purposes is hereby granted *\r
-* without fee, provided that the above copyright notice appears in all *\r
-* copies and that both the copyright notice and this permission notice *\r
-* appear in the supporting documentation. The authors make no claims *\r
-* about the suitability of this software for any purpose. It is *\r
-* provided "as is" without express or implied warranty. * **************************************************************************/\r
-\r
-\r
-#include "AliAnalysisTaskQAHighPtDeDx.h"\r
-\r
-// ROOT includes\r
-#include <TList.h>\r
-#include <TTree.h>\r
-#include <TMath.h>\r
-#include <TH1.h>\r
-#include <TF1.h>\r
-#include <TProfile.h>\r
-#include <TParticle.h>\r
-#include <TFile.h>\r
-\r
-// AliRoot includes\r
-#include <AliAnalysisManager.h>\r
-#include <AliAnalysisFilter.h>\r
-#include <AliESDInputHandler.h>\r
-#include <AliESDEvent.h>\r
-#include <AliESDVertex.h>\r
-#include <AliLog.h>\r
-#include <AliExternalTrackParam.h>\r
-#include <AliESDtrackCuts.h>\r
-#include <AliESDVZERO.h>\r
-#include <AliAODVZERO.h>\r
-\r
-#include <AliMCEventHandler.h>\r
-#include <AliMCEvent.h>\r
-#include <AliStack.h>\r
-\r
-#include <TTreeStream.h>\r
-\r
-#include <AliHeader.h>\r
-#include <AliGenPythiaEventHeader.h>\r
-#include <AliGenDPMjetEventHeader.h>\r
-\r
-#include "AliCentrality.h" \r
-#include <AliESDv0.h>\r
-#include <AliKFVertex.h>\r
-#include <AliAODVertex.h>\r
-\r
-#include <AliAODTrack.h> \r
-#include <AliAODPid.h> \r
-#include <AliAODMCHeader.h> \r
-\r
-\r
-// STL includes\r
-#include <iostream>\r
-using namespace std;\r
-\r
-\r
-//\r
-// Responsible:\r
-// Antonio Ortiz (Lund)\r
-// Peter Christiansen (Lund)\r
-//\r
-\r
-\r
-\r
-\r
-\r
-\r
-\r
-\r
-const Double_t AliAnalysisTaskQAHighPtDeDx::fgkClight = 2.99792458e-2;\r
-Float_t magf = -1;\r
-TF1* cutLow = new TF1("StandardPhiCutLow", "0.1/x/x+pi/18.0-0.025", 0, 50);\r
-TF1* cutHigh = new TF1("StandardPhiCutHigh", "0.12/x+pi/18.0+0.035", 0, 50);\r
-Double_t DeDxMIPMin = 30;\r
-Double_t DeDxMIPMax = 65;\r
-const Int_t nHists = 9;\r
-Float_t centralityGlobal = -10;\r
-Int_t etaLow[nHists] = {-8, -8, -6, -4, -2, 0, 2, 4, 6};\r
-Int_t etaHigh[nHists] = { 8, -6, -4, -2, 0, 2, 4, 6, 8};\r
-\r
-Int_t nDeltaPiBins = 80;\r
-Double_t deltaPiLow = 20;\r
-Double_t deltaPiHigh = 100;\r
-const Char_t *Pid[7]={"Ch","Pion","Kaon","Proton","Electron","Muon","Oher"};\r
-ClassImp(AliAnalysisTaskQAHighPtDeDx)\r
-//_____________________________________________________________________________\r
-//AliAnalysisTaskQAHighPtDeDx::AliAnalysisTaskQAHighPtDeDx(const char *name):\r
-AliAnalysisTaskQAHighPtDeDx::AliAnalysisTaskQAHighPtDeDx():\r
- AliAnalysisTaskSE(),\r
- fESD(0x0),\r
- fAOD(0x0),\r
- fMC(0x0),\r
- fMCStack(0x0),\r
- fMCArray(0x0),\r
- fTrackFilterGolden(0x0),\r
- fTrackFilterTPC(0x0),\r
- fCentEst("V0M"),\r
- fAnalysisType("ESD"),\r
- fAnalysisMC(kFALSE),\r
- fAnalysisPbPb(kFALSE),\r
- ftrigBit(0x0),\r
- fRandom(0x0),\r
- fPileUpRej(kFALSE),\r
- fVtxCut(10.0), \r
- fEtaCut(0.9), \r
- fMinCent(0.0),\r
- fMaxCent(100.0),\r
- fStoreMcIn(kFALSE),//\r
- fMcProcessType(-999),\r
- fTriggeredEventMB(-999),\r
- fVtxStatus(-999),\r
- fZvtx(-999),\r
- fZvtxMC(-999),\r
- fRun(-999),\r
- fEventId(-999),\r
- fListOfObjects(0), \r
- fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),\r
- fn1(0x0),\r
- fcent(0x0),\r
- hMIPVsEta(0x0),\r
- pMIPVsEta(0x0),\r
- hMIPVsEtaV0s(0x0),\r
- pMIPVsEtaV0s(0x0),\r
- hPlateauVsEta(0x0),\r
- pPlateauVsEta(0x0),\r
- hPhi(0x0)\r
-\r
-\r
-{\r
- //default constructor\r
-}\r
-\r
-\r
-AliAnalysisTaskQAHighPtDeDx::AliAnalysisTaskQAHighPtDeDx(const char *name):\r
- AliAnalysisTaskSE(name),\r
- fESD(0x0),\r
- fAOD(0x0),\r
- fMC(0x0),\r
- fMCStack(0x0),\r
- fMCArray(0x0),\r
- fTrackFilterGolden(0x0),\r
- fTrackFilterTPC(0x0),\r
- fCentEst("V0M"),\r
- fAnalysisType("ESD"),\r
- fAnalysisMC(kFALSE),\r
- fAnalysisPbPb(kFALSE),\r
- ftrigBit(0x0),\r
- fRandom(0x0),\r
- fPileUpRej(kFALSE),\r
- fVtxCut(10.0), \r
- fEtaCut(0.9), \r
- fMinCent(0.0),\r
- fMaxCent(100.0),\r
- fStoreMcIn(kFALSE),//\r
- fMcProcessType(-999),\r
- fTriggeredEventMB(-999),\r
- fVtxStatus(-999),\r
- fZvtx(-999),\r
- fZvtxMC(-999),\r
- fRun(-999),\r
- fEventId(-999),\r
- fListOfObjects(0), \r
- fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),\r
- fn1(0x0),\r
- fcent(0x0),\r
- hMIPVsEta(0x0),\r
- pMIPVsEta(0x0),\r
- hMIPVsEtaV0s(0x0),\r
- pMIPVsEtaV0s(0x0),\r
- hPlateauVsEta(0x0),\r
- pPlateauVsEta(0x0),\r
- hPhi(0x0)\r
-\r
-\r
-{\r
- // Default constructor (should not be used)\r
- for(Int_t i=0;i<9;++i){\r
- \r
- hMIPVsNch[i]=0;//TH2D, MIP vs Nch for different eta intervals\r
- pMIPVsNch[i]=0;//TProfile, MIP vs Nch for different eta intervals\r
- hMIPVsPhi[i]=0;//TH2D, MIP vs phi for different eta intervals\r
- pMIPVsPhi[i]=0;//TProfile, MIP vs phi for different eta intervals\r
- hPlateauVsPhi[i]=0;//TH2D, dE/dx vs Phi, electrons 0.4<p<0.6 GeV/c\r
- pPlateauVsPhi[i]=0;//TProfile, dE/dx vs Phi, electrons 0.4<p<0.6 GeV/c\r
- histPiV0[i]=0;//TH2D, dE/dx vs p, pi id by V0s\r
- histpPiV0[i]=0;//TH1D, pi id by V0s\r
- histPV0[i]=0;// TH2D, dE/dx vs p, p id by V0s\r
- histpPV0[i]=0;// TH1D, p id by V0s\r
- histAllCh[i]=0;//TH2D, dE/dx vs p for all charged particles\r
- histPiTof[i]=0;//TH2D, dE/dx vs p for a "clean" sample of pions, beta>1\r
- histpPiTof[i]=0;//TH1D, for a "clean" sample of pions, beta>1\r
- histEV0[i]=0;\r
-\r
- for(Int_t pid=0;pid<7;++pid){\r
- hMcIn[pid][i]=0;\r
- hMcOut[pid][i]=0;\r
- }\r
-\r
- }\r
- DefineOutput(1, TList::Class());//esto es nuevo\r
-}\r
-\r
-\r
-\r
-\r
-AliAnalysisTaskQAHighPtDeDx::~AliAnalysisTaskQAHighPtDeDx() {\r
- //\r
- // Destructor\r
- //\r
-\r
-}\r
-\r
-\r
-\r
-\r
-\r
-\r
-\r
-\r
-//______________________________________________________________________________\r
-void AliAnalysisTaskQAHighPtDeDx::UserCreateOutputObjects()\r
-{ \r
- // This method is called once per worker node\r
- // Here we define the output: histograms and debug tree if requested \r
- // We also create the random generator here so it might get different seeds...\r
- fRandom = new TRandom(0); // 0 means random seed\r
- \r
-\r
-\r
- //OpenFile(1);\r
- fListOfObjects = new TList();\r
- fListOfObjects->SetOwner();\r
- \r
-\r
-\r
-\r
- //\r
- // Histograms\r
- // \r
- fEvents = new TH1I("fEvents","Number of analyzed events; Events; Counts", 3, 0, 3);\r
- fListOfObjects->Add(fEvents);\r
-\r
- fn1=new TH1F("fn1","fn1",11,-1,10);\r
- fListOfObjects->Add(fn1);\r
-\r
- fcent=new TH1F("fcent","fcent",104,-2,102);\r
- fListOfObjects->Add(fcent);\r
-\r
- fVtx = new TH1I("fVtx","Vtx info (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);\r
- fListOfObjects->Add(fVtx);\r
-\r
- fVtxBeforeCuts = new TH1F("fVtxBeforeCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);\r
- fListOfObjects->Add(fVtxBeforeCuts);\r
- \r
- fVtxAfterCuts = new TH1F("fVtxAfterCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);\r
- fListOfObjects->Add(fVtxAfterCuts);\r
-\r
-\r
- const Int_t nPtBinsV0s = 25;\r
- 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 ,\r
- 1.2 , 1.4 , 1.6 , 1.8 , 2.0 , 2.5 , 3.0 , 3.5 , 4.0 , 5.0 , 7.0 ,\r
- 9.0 , 12.0, 15.0, 20.0 };\r
- \r
-\r
-\r
-\r
- const Char_t* ending[nHists] = {"", "86", "64", "42", "20", "02", "24", "46", "68"};\r
- \r
- const Char_t* LatexEta[nHists] = {\r
- "|#eta|<0.8", "-0.8<#eta<-0.6", "-0.6<#eta<-0.4", "-0.4<#eta<-0.2", "-0.2<#eta<0", \r
- "0<#eta<0.2", "0.2<#eta<0.4", "0.4<#eta<0.6", "0.6<#eta<0.8" \r
- };\r
- hPhi = new TH2D("histPhi", "pt; #phi'", nPtBinsV0s, ptBinsV0s, 90, -0.05, 0.4); \r
- //dE/dx vs phi, pions at the MIP\r
- fListOfObjects->Add(hPhi);\r
-\r
-\r
- \r
-\r
- Int_t nPhiBins = 36;\r
- \r
- for(Int_t i = 0; i < nHists; i++) {\r
- \r
- \r
- hMIPVsPhi[i] = new TH2D(Form("hMIPVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx MIP",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(), \r
- DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);\r
- hMIPVsPhi[i]->Sumw2();\r
- \r
- pMIPVsPhi[i] = new TProfile(Form("pMIPVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx MIP",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(),\r
- DeDxMIPMin, DeDxMIPMax);\r
- pMIPVsPhi[i]->SetMarkerStyle(20);\r
- pMIPVsPhi[i]->SetMarkerColor(1);\r
- pMIPVsPhi[i]->SetLineColor(1);\r
- pMIPVsPhi[i]->Sumw2();\r
-\r
- fListOfObjects->Add(hMIPVsPhi[i]);\r
- fListOfObjects->Add(pMIPVsPhi[i]);\r
-\r
- hPlateauVsPhi[i] = new TH2D(Form("hPlateauVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx Plateau",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(),\r
- 95-DeDxMIPMax, DeDxMIPMax, 95);\r
- hPlateauVsPhi[i]->Sumw2();\r
- \r
- pPlateauVsPhi[i] = new TProfile(Form("pPlateauVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx Plateau",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(),\r
- DeDxMIPMax, 95);\r
- pPlateauVsPhi[i]->SetMarkerStyle(20);\r
- pPlateauVsPhi[i]->SetMarkerColor(1);\r
- pPlateauVsPhi[i]->SetLineColor(1);\r
- pPlateauVsPhi[i]->Sumw2();\r
-\r
- fListOfObjects->Add(hPlateauVsPhi[i]);\r
- fListOfObjects->Add(pPlateauVsPhi[i]);\r
-\r
-\r
- 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, \r
- DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);\r
- hMIPVsNch[i]->Sumw2();\r
- \r
- 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);\r
- pMIPVsNch[i]->SetMarkerStyle(20);\r
- pMIPVsNch[i]->SetMarkerColor(1);\r
- pMIPVsNch[i]->SetLineColor(1);\r
- pMIPVsNch[i]->Sumw2();\r
-\r
- fListOfObjects->Add(hMIPVsNch[i]);\r
- fListOfObjects->Add(pMIPVsNch[i]);\r
-\r
- //two dimmesional distributions dE/dx vs p for secondary pions\r
- histPiV0[i] = new TH2D(Form("histPiV0%s", ending[i]), "Pions id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);\r
- histPiV0[i]->Sumw2();\r
- fListOfObjects->Add(histPiV0[i]);\r
-\r
- histpPiV0[i] = new TH1D(Form("histPiV01D%s", ending[i]), "Pions id by V0; #it{p} (GeV/#it{c}); counts", 200, 0, 20);\r
- histpPiV0[i]->Sumw2();\r
- fListOfObjects->Add(histpPiV0[i]);\r
-\r
- //two dimmesional distributions dE/dx vs p for secondary protons\r
- histPV0[i] = new TH2D(Form("histPV0%s", ending[i]), "Protons id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);\r
- histPV0[i]->Sumw2();\r
- fListOfObjects->Add(histPV0[i]);\r
-\r
- histpPV0[i] = new TH1D(Form("histPV01D%s", ending[i]), "Protons id by V0; #it{p} (GeV/#it{c}); counts", 200, 0, 20);\r
- histpPV0[i]->Sumw2();\r
- fListOfObjects->Add(histpPV0[i]);\r
-\r
- //two dimmesional distributions dE/dx vs p for primary pions\r
- histPiTof[i] = new TH2D(Form("histPiTof%s", ending[i]), "all charged", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);\r
- histPiTof[i]->Sumw2();\r
-\r
- histpPiTof[i] = new TH1D(Form("histPiTof1D%s", ending[i]), "Protons id by V0; #it{p} (GeV/#it{c}); counts", 200, 0, 20);\r
- histpPiTof[i]->Sumw2();\r
- fListOfObjects->Add(histpPiTof[i]);\r
-\r
-\r
- histAllCh[i] = new TH2D(Form("histAllCh%s", ending[i]), "Pions id by TOF", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);\r
- histAllCh[i]->Sumw2();\r
-\r
- fListOfObjects->Add(histPiTof[i]);\r
- fListOfObjects->Add(histAllCh[i]);\r
-\r
- histEV0[i] = new TH2D(Form("histEV0%s", ending[i]), "Electrons id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);\r
- histEV0[i]->Sumw2();\r
- fListOfObjects->Add(histEV0[i]);\r
-\r
-\r
-\r
- }\r
-\r
-\r
- hMIPVsEta = new TH2D("hMIPVsEta","; #eta; dE/dx_{MIP, primary tracks}",16,-0.8,0.8,DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);\r
- pMIPVsEta = new TProfile("pMIPVsEta","; #eta; #LT dE/dx #GT_{MIP, primary tracks}",16,-0.8,0.8, DeDxMIPMin, DeDxMIPMax);\r
- hMIPVsEtaV0s = new TH2D("hMIPVsEtaV0s","; #eta; dE/dx_{MIP, secondary tracks}",16,-0.8,0.8,DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);\r
- pMIPVsEtaV0s = new TProfile("pMIPVsEtaV0s","; #eta; #LT dE/dx #GT_{MIP, secondary tracks}",16,-0.8,0.8,DeDxMIPMin, DeDxMIPMax);\r
- \r
- hPlateauVsEta = new TH2D("hPlateauVsEta","; #eta; dE/dx_{Plateau, primary tracks}",16,-0.8,0.8,95-DeDxMIPMax, DeDxMIPMax, 95);\r
- pPlateauVsEta = new TProfile("pPlateauVsEta","; #eta; #LT dE/dx #GT_{Plateau, primary tracks}",16,-0.8,0.8, DeDxMIPMax, 95);\r
-\r
- fListOfObjects->Add(hMIPVsEta);\r
- fListOfObjects->Add(pMIPVsEta);\r
- fListOfObjects->Add(hMIPVsEtaV0s);\r
- fListOfObjects->Add(pMIPVsEtaV0s);\r
- fListOfObjects->Add(hPlateauVsEta);\r
- fListOfObjects->Add(pPlateauVsEta);\r
-\r
-\r
-\r
-\r
-\r
-\r
- if (fAnalysisMC) { \r
- for(Int_t i = 0; i < nHists; i++) {\r
- for(Int_t pid = 0; pid < 7; pid++) {\r
- \r
- hMcIn[pid][i] = new TH1D(Form("hIn%s%s", Pid[pid],ending[i]), Form("MC in (pid %s)", Pid[pid]), \r
- nPtBinsV0s, ptBinsV0s);\r
- fListOfObjects->Add(hMcIn[pid][i]);\r
- \r
- hMcOut[pid][i] = new TH1D(Form("hMcOut%s%s", Pid[pid],ending[i]), Form("MC out (pid %s)", Pid[pid]), \r
- nPtBinsV0s, ptBinsV0s);\r
- fListOfObjects->Add(hMcOut[pid][i]);\r
- \r
- \r
- }\r
- }\r
-\r
- fVtxMC = new TH1F("fVtxMC","mc vtx", 120, -30, 30);\r
- fListOfObjects->Add(fVtxMC);\r
- \r
- }\r
- \r
- // Post output data.\r
- PostData(1, fListOfObjects);\r
-}\r
-\r
-//______________________________________________________________________________\r
-void AliAnalysisTaskQAHighPtDeDx::UserExec(Option_t *) \r
-{\r
- // Main loop\r
-\r
- //\r
- // First we make sure that we have valid input(s)!\r
- //\r
-\r
-\r
-\r
- AliVEvent *event = InputEvent();\r
- if (!event) {\r
- Error("UserExec", "Could not retrieve event");\r
- return;\r
- }\r
-\r
-\r
-\r
- if (fAnalysisType == "ESD"){\r
- fESD = dynamic_cast<AliESDEvent*>(event);\r
- if(!fESD){\r
- Printf("%s:%d ESDEvent not found in Input Manager",(char*)__FILE__,__LINE__);\r
- this->Dump();\r
- return;\r
- } \r
- } else {\r
- fAOD = dynamic_cast<AliAODEvent*>(event);\r
- if(!fAOD){\r
- Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);\r
- this->Dump();\r
- return;\r
- } \r
- }\r
-\r
-\r
-\r
- if (fAnalysisMC) {\r
-\r
- if (fAnalysisType == "ESD"){\r
- fMC = dynamic_cast<AliMCEvent*>(MCEvent());\r
- if(!fMC){\r
- Printf("%s:%d MCEvent not found in Input Manager",(char*)__FILE__,__LINE__);\r
- this->Dump();\r
- return;\r
- } \r
-\r
- fMCStack = fMC->Stack();\r
- \r
- if(!fMCStack){\r
- Printf("%s:%d MCStack not found in Input Manager",(char*)__FILE__,__LINE__);\r
- this->Dump();\r
- return;\r
- } \r
- } else { // AOD\r
-\r
- fMC = dynamic_cast<AliMCEvent*>(MCEvent());\r
- if(fMC)\r
- fMC->Dump();\r
-\r
- fMCArray = (TClonesArray*)fAOD->FindListObject("mcparticles");\r
- if(!fMCArray){\r
- Printf("%s:%d AOD MC array not found in Input Manager",(char*)__FILE__,__LINE__);\r
- this->Dump();\r
- return;\r
- } \r
- }\r
- }\r
-\r
- \r
- // Get trigger decision\r
- fTriggeredEventMB = 0; //init\r
- \r
- fn1->Fill(0);\r
- \r
- if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))\r
- ->IsEventSelected() & ftrigBit ){\r
- fTriggeredEventMB = 1; //event triggered as minimum bias\r
- }\r
-\r
- // Get process type for MC\r
- fMcProcessType = 0; // -1=invalid, 0=data, 1=ND, 2=SD, 3=DD\r
-\r
- // real data that are not triggered we skip\r
- if(!fAnalysisMC && !fTriggeredEventMB)\r
- return; \r
- \r
- fn1->Fill(1);\r
-\r
-\r
- if (fAnalysisMC) {\r
- \r
-\r
-\r
- if (fAnalysisType == "ESD"){\r
-\r
- \r
-\r
- AliHeader* headerMC = fMC->Header();\r
- if (headerMC) {\r
- \r
- AliGenEventHeader* genHeader = headerMC->GenEventHeader();\r
- TArrayF vtxMC(3); // primary vertex MC \r
- vtxMC[0]=9999; vtxMC[1]=9999; vtxMC[2]=9999; //initialize with dummy\r
- if (genHeader) {\r
- genHeader->PrimaryVertex(vtxMC);\r
- }\r
- fZvtxMC = vtxMC[2];\r
- \r
- // PYTHIA:\r
- AliGenPythiaEventHeader* pythiaGenHeader =\r
- dynamic_cast<AliGenPythiaEventHeader*>(headerMC->GenEventHeader());\r
- if (pythiaGenHeader) { //works only for pythia\r
- fMcProcessType = GetPythiaEventProcessType(pythiaGenHeader->ProcessType());\r
- }\r
- // PHOJET:\r
- AliGenDPMjetEventHeader* dpmJetGenHeader =\r
- dynamic_cast<AliGenDPMjetEventHeader*>(headerMC->GenEventHeader());\r
- if (dpmJetGenHeader) {\r
- fMcProcessType = GetDPMjetEventProcessType(dpmJetGenHeader->ProcessType());\r
- }\r
- }\r
- } else { // AOD\r
- \r
- \r
-\r
- AliAODMCHeader* mcHeader = dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject("mcHeader")); \r
-\r
-\r
- if(mcHeader) {\r
- fZvtxMC = mcHeader->GetVtxZ();\r
- \r
-\r
-\r
- if(strstr(mcHeader->GetGeneratorName(), "Pythia")) {\r
- fMcProcessType = GetPythiaEventProcessType(mcHeader->GetEventType());\r
- } else {\r
- fMcProcessType = GetDPMjetEventProcessType(mcHeader->GetEventType());\r
- }\r
- }\r
- }\r
-\r
-\r
- }\r
- \r
- \r
-\r
- if (fAnalysisType == "ESD"){\r
- \r
- const AliESDVertex *vtxESD = fESD->GetPrimaryVertexTracks();\r
- if(vtxESD->GetNContributors()<1) {\r
- // SPD vertex\r
- vtxESD = fESD->GetPrimaryVertexSPD();\r
- /* quality checks on SPD-vertex */\r
- TString vertexType = vtxESD->GetTitle();\r
- if (vertexType.Contains("vertexer: Z") && (vtxESD->GetDispersion() > 0.04 || vtxESD->GetZRes() > 0.25)) \r
- fZvtx = -1599; //vertex = 0x0; //\r
- else if (vtxESD->GetNContributors()<1) \r
- fZvtx = -999; //vertex = 0x0; //\r
- else\r
- fZvtx = vtxESD->GetZ();\r
- } \r
- else\r
- fZvtx = vtxESD->GetZ();\r
- \r
- }\r
- else // AOD\r
- fZvtx = GetVertex(fAOD);\r
- \r
- fVtxBeforeCuts->Fill(fZvtx);\r
- \r
- //cut on the z position of vertex\r
- if (TMath::Abs(fZvtx) > fVtxCut) { \r
- return;\r
- }\r
- fn1->Fill(2);\r
- \r
-\r
-\r
-\r
-\r
-\r
- Float_t centrality = -10;\r
-\r
- // only analyze triggered events\r
- if(fTriggeredEventMB) {\r
- \r
- if (fAnalysisType == "ESD"){\r
- if(fAnalysisPbPb){\r
- AliCentrality *centObject = fESD->GetCentrality();\r
- centrality = centObject->GetCentralityPercentile(fCentEst);\r
- centralityGlobal = centrality;\r
- if((centrality>fMaxCent)||(centrality<fMinCent))return;\r
- }\r
- fcent->Fill(centrality);\r
- fn1->Fill(3);\r
- if(fAnalysisMC){\r
- if(TMath::Abs(fZvtxMC)<fVtxCut){\r
- ProcessMCTruthESD();\r
- fVtxMC->Fill(fZvtxMC);\r
- }\r
- }\r
- AnalyzeESD(fESD);\r
- } else { // AOD\r
- if(fAnalysisPbPb){\r
- AliCentrality *centObject = fAOD->GetCentrality();\r
- if(centObject){\r
- centrality = centObject->GetCentralityPercentile(fCentEst); \r
- }\r
- //cout<<"centrality="<<centrality<<endl;\r
- if((centrality>fMaxCent)||(centrality<fMinCent))return;\r
- }\r
- fcent->Fill(centrality);\r
- fn1->Fill(3);\r
- if(fAnalysisMC){\r
- if(TMath::Abs(fZvtxMC)<fVtxCut){\r
-\r
- ProcessMCTruthAOD();\r
- fVtxMC->Fill(fZvtxMC);\r
- }\r
- }\r
- AnalyzeAOD(fAOD);\r
- }\r
- }\r
-\r
- fVtxAfterCuts->Fill(fZvtx);\r
-\r
- \r
- \r
-\r
- // Post output data.\r
- PostData(1, fListOfObjects);\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskQAHighPtDeDx::AnalyzeESD(AliESDEvent* esdEvent)\r
-{\r
- fRun = esdEvent->GetRunNumber();\r
- fEventId = 0;\r
- if(esdEvent->GetHeader())\r
- fEventId = GetEventIdAsLong(esdEvent->GetHeader());\r
- \r
- cout << "centrality=" << centralityGlobal << endl;\r
-\r
-\r
- Bool_t isPileup = esdEvent->IsPileupFromSPD();\r
- if(fPileUpRej)\r
- if(isPileup)\r
- return;\r
- fn1->Fill(4);\r
-\r
-\r
- // Int_t event = esdEvent->GetEventNumberInFile();\r
- //UInt_t time = esdEvent->GetTimeStamp();\r
- // ULong64_t trigger = esdEvent->GetTriggerMask();\r
- magf = esdEvent->GetMagneticField();\r
-\r
-\r
-\r
-\r
-\r
- if(fTriggeredEventMB) {// Only MC case can we have not triggered events\r
- \r
- // accepted event\r
- fEvents->Fill(0);\r
- \r
- \r
- //Change, 10/04/13. Now accept all events to do a correct normalization\r
- //if(fVtxStatus!=1) return; // accepted vertex\r
- //Int_t nESDTracks = esdEvent->GetNumberOfTracks();\r
- \r
- ProduceArrayTrksESD( esdEvent );//produce array with global track parameters\r
- ProduceArrayV0ESD( esdEvent );//v0's\r
-\r
- \r
- fEvents->Fill(1);\r
-\r
-\r
-\r
-\r
- } // end if triggered\r
-\r
-\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskQAHighPtDeDx::AnalyzeAOD(AliAODEvent* aodEvent)\r
-{\r
- fRun = aodEvent->GetRunNumber();\r
- fEventId = 0;\r
- if(aodEvent->GetHeader())\r
- fEventId = GetEventIdAsLong(aodEvent->GetHeader());\r
- \r
- //UInt_t time = 0; // Missing AOD info? aodEvent->GetTimeStamp();\r
- magf = aodEvent->GetMagneticField();\r
-\r
- //Int_t trackmult = 0; // no pt cuts\r
- //Int_t nadded = 0;\r
-\r
- Bool_t isPileup = aodEvent->IsPileupFromSPD();\r
- if(fPileUpRej)\r
- if(isPileup)\r
- return;\r
- fn1->Fill(4);\r
-\r
-\r
-\r
- if(fTriggeredEventMB) {// Only MC case can we have not triggered events\r
- \r
- // accepted event\r
- fEvents->Fill(0);\r
- \r
- //if(fVtxStatus!=1) return; // accepted vertex\r
- //Int_t nAODTracks = aodEvent->GetNumberOfTracks(); \r
-\r
- ProduceArrayTrksAOD( aodEvent );\r
- ProduceArrayV0AOD( aodEvent );\r
-\r
- fEvents->Fill(1);\r
- \r
-\r
-\r
-\r
- } // end if triggered\r
-\r
-}\r
-\r
-//_____________________________________________________________________________\r
-Float_t AliAnalysisTaskQAHighPtDeDx::GetVertex(const AliVEvent* event) const\r
-{\r
- Float_t zvtx = -999;\r
- \r
- const AliVVertex* primaryVertex = event->GetPrimaryVertex(); \r
- \r
- if(primaryVertex->GetNContributors()>0)\r
- zvtx = primaryVertex->GetZ();\r
-\r
- return zvtx;\r
-}\r
-\r
-//_____________________________________________________________________________\r
-Short_t AliAnalysisTaskQAHighPtDeDx::GetPidCode(Int_t pdgCode) const \r
-{\r
- // return our internal code for pions, kaons, and protons\r
- \r
- Short_t pidCode = 6;\r
- \r
- switch (TMath::Abs(pdgCode)) {\r
- case 211:\r
- pidCode = 1; // pion\r
- break;\r
- case 321:\r
- pidCode = 2; // kaon\r
- break;\r
- case 2212:\r
- pidCode = 3; // proton\r
- break;\r
- case 11:\r
- pidCode = 4; // electron\r
- break;\r
- case 13:\r
- pidCode = 5; // muon\r
- break;\r
- default:\r
- pidCode = 6; // something else?\r
- };\r
- \r
- return pidCode;\r
-}\r
-\r
-//_____________________________________________________________________________\r
-void AliAnalysisTaskQAHighPtDeDx::ProcessMCTruthESD() \r
-{\r
- // Fill the special MC histogram with the MC truth info\r
- \r
- const Int_t nTracksMC = fMCStack->GetNtrack();\r
-\r
- for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {\r
- \r
- //Cuts\r
- if(!(fMCStack->IsPhysicalPrimary(iTracks)))\r
- continue;\r
- \r
- TParticle* trackMC = fMCStack->Particle(iTracks);\r
- \r
- TParticlePDG* pdgPart = trackMC ->GetPDG();\r
- Double_t chargeMC = pdgPart->Charge();\r
-\r
- if(chargeMC==0)\r
- continue;\r
- \r
- if (TMath::Abs(trackMC->Eta()) > fEtaCut )\r
- continue;\r
-\r
- Double_t etaMC = trackMC->Eta();\r
- Int_t pdgCode = trackMC->GetPdgCode();\r
- Short_t pidCodeMC = 0;\r
- pidCodeMC = GetPidCode(pdgCode);\r
- \r
- \r
- for(Int_t nh = 0; nh < 9; nh++) {\r
- \r
- if( etaMC > etaHigh[nh]/10.0 || etaMC < etaLow[nh]/10.0 )\r
- continue;\r
- \r
- hMcIn[0][nh]->Fill(trackMC->Pt());\r
- hMcIn[pidCodeMC][nh]->Fill(trackMC->Pt());\r
- \r
- \r
- }\r
-\r
- }//MC track loop\r
- \r
- \r
- \r
-}\r
-\r
-//_____________________________________________________________________________\r
-void AliAnalysisTaskQAHighPtDeDx::ProcessMCTruthAOD() \r
-{\r
- // Fill the special MC histogram with the MC truth info\r
- \r
-\r
- const Int_t nTracksMC = fMCArray->GetEntriesFast();\r
- \r
- for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {\r
- \r
- AliAODMCParticle* trackMC = dynamic_cast<AliAODMCParticle*>(fMCArray->At(iTracks));\r
-\r
- if(!trackMC){\r
- AliError("Cannot get MC particle");\r
- continue;\r
- } \r
- \r
- //Cuts\r
- if(!(trackMC->IsPhysicalPrimary()))\r
- continue;\r
- \r
- \r
- Double_t chargeMC = trackMC->Charge();\r
- if(chargeMC==0)\r
- continue;\r
-\r
- \r
- if (TMath::Abs(trackMC->Eta()) > fEtaCut )\r
- continue;\r
- \r
- Double_t etaMC = trackMC->Eta();\r
- Int_t pdgCode = trackMC->GetPdgCode();\r
- Short_t pidCodeMC = 0;\r
- pidCodeMC = GetPidCode(pdgCode);\r
- \r
- //cout<<"pidcode="<<pidCodeMC<<endl;\r
- for(Int_t nh = 0; nh < 9; nh++) {\r
- \r
- if( etaMC > etaHigh[nh]/10.0 || etaMC < etaLow[nh]/10.0 )\r
- continue;\r
- \r
- hMcIn[0][nh]->Fill(trackMC->Pt());\r
- hMcIn[pidCodeMC][nh]->Fill(trackMC->Pt());\r
- \r
- \r
- }\r
- \r
- }//MC track loop\r
- \r
- \r
-}\r
-\r
-\r
-//_____________________________________________________________________________\r
-Short_t AliAnalysisTaskQAHighPtDeDx::GetPythiaEventProcessType(Int_t pythiaType) {\r
- //\r
- // Get the process type of the event. PYTHIA\r
- //\r
- // source PWG0 dNdpt \r
-\r
- Short_t globalType = -1; //init\r
- \r
- if(pythiaType==92||pythiaType==93){\r
- globalType = 2; //single diffractive\r
- }\r
- else if(pythiaType==94){\r
- globalType = 3; //double diffractive\r
- }\r
- //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??\r
- else {\r
- globalType = 1; //non diffractive\r
- }\r
- return globalType;\r
-}\r
-\r
-//_____________________________________________________________________________\r
-Short_t AliAnalysisTaskQAHighPtDeDx::GetDPMjetEventProcessType(Int_t dpmJetType) {\r
- //\r
- // get the process type of the event. PHOJET\r
- //\r
- //source PWG0 dNdpt \r
- // can only read pythia headers, either directly or from cocktalil header\r
- Short_t globalType = -1;\r
- \r
- if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction\r
- globalType = 1;\r
- }\r
- else if (dpmJetType==5 || dpmJetType==6) {\r
- globalType = 2;\r
- }\r
- else if (dpmJetType==7) {\r
- globalType = 3;\r
- }\r
- return globalType;\r
-}\r
-\r
-//_____________________________________________________________________________\r
-ULong64_t AliAnalysisTaskQAHighPtDeDx::GetEventIdAsLong(AliVHeader* header) const\r
-{\r
- // To have a unique id for each event in a run!\r
- // Modified from AliRawReader.h\r
- return ((ULong64_t)header->GetBunchCrossNumber()+\r
- (ULong64_t)header->GetOrbitNumber()*3564+\r
- (ULong64_t)header->GetPeriodNumber()*16777215*3564);\r
-}\r
-\r
-\r
-//____________________________________________________________________\r
-TParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMother(AliStack* stack, Int_t label)\r
-{\r
- //\r
- // Finds the first mother among the primary particles of the particle identified by <label>,\r
- // i.e. the primary that "caused" this particle\r
- //\r
- // Taken from AliPWG0Helper class\r
- //\r
-\r
- Int_t motherLabel = FindPrimaryMotherLabel(stack, label);\r
- if (motherLabel < 0)\r
- return 0;\r
-\r
- return stack->Particle(motherLabel);\r
-}\r
-\r
-//____________________________________________________________________\r
-Int_t AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherLabel(AliStack* stack, Int_t label)\r
-{\r
- //\r
- // Finds the first mother among the primary particles of the particle identified by <label>,\r
- // i.e. the primary that "caused" this particle\r
- //\r
- // returns its label\r
- //\r
- // Taken from AliPWG0Helper class\r
- //\r
- const Int_t nPrim = stack->GetNprimary();\r
- \r
- while (label >= nPrim) {\r
-\r
- //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));\r
-\r
- TParticle* particle = stack->Particle(label);\r
- if (!particle) {\r
- \r
- AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));\r
- return -1;\r
- }\r
- \r
- // find mother\r
- if (particle->GetMother(0) < 0) {\r
-\r
- AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));\r
- return -1;\r
- }\r
- \r
- label = particle->GetMother(0);\r
- }\r
- \r
- return label;\r
-}\r
-\r
-//____________________________________________________________________\r
-AliAODMCParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherAOD(AliAODMCParticle* startParticle)\r
-{\r
- //\r
- // Finds the first mother among the primary particles of the particle identified by <label>,\r
- // i.e. the primary that "caused" this particle\r
- //\r
- // Taken from AliPWG0Helper class\r
- //\r
-\r
- AliAODMCParticle* mcPart = startParticle;\r
-\r
- while (mcPart)\r
- {\r
- \r
- if(mcPart->IsPrimary())\r
- return mcPart;\r
- \r
- Int_t mother = mcPart->GetMother();\r
-\r
- mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));\r
- }\r
-\r
- return 0;\r
-}\r
-\r
-\r
-//V0______________________________________\r
-//____________________________________________________________________\r
-TParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherV0(AliStack* stack, Int_t label)\r
-{\r
- //\r
- // Finds the first mother among the primary particles of the particle identified by <label>,\r
- // i.e. the primary that "caused" this particle\r
- //\r
- // Taken from AliPWG0Helper class\r
- //\r
-\r
- Int_t nSteps = 0;\r
-\r
- Int_t motherLabel = FindPrimaryMotherLabelV0(stack, label, nSteps);\r
- if (motherLabel < 0)\r
- return 0;\r
-\r
- return stack->Particle(motherLabel);\r
-}\r
-\r
-//____________________________________________________________________\r
-Int_t AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherLabelV0(AliStack* stack, Int_t label, Int_t& nSteps)\r
-{\r
- //\r
- // Finds the first mother among the primary particles of the particle identified by <label>,\r
- // i.e. the primary that "caused" this particle\r
- //\r
- // returns its label\r
- //\r
- // Taken from AliPWG0Helper class\r
- //\r
- nSteps = 0;\r
- const Int_t nPrim = stack->GetNprimary();\r
- \r
- while (label >= nPrim) {\r
-\r
- //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));\r
-\r
- nSteps++; // 1 level down\r
- \r
- TParticle* particle = stack->Particle(label);\r
- if (!particle) {\r
- \r
- AliDebugGeneral("FindPrimaryMotherLabelV0", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));\r
- return -1;\r
- }\r
- \r
- // find mother\r
- if (particle->GetMother(0) < 0) {\r
-\r
- AliDebugGeneral("FindPrimaryMotherLabelV0", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));\r
- return -1;\r
- }\r
- \r
- label = particle->GetMother(0);\r
- }\r
- \r
- return label;\r
-}\r
-\r
-//____________________________________________________________________\r
-AliAODMCParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherAODV0(AliAODMCParticle* startParticle, Int_t& nSteps)\r
-{\r
- //\r
- // Finds the first mother among the primary particles of the particle identified by <label>,\r
- // i.e. the primary that "caused" this particle\r
- //\r
- // Taken from AliPWG0Helper class\r
- //\r
-\r
- nSteps = 0;\r
-\r
- AliAODMCParticle* mcPart = startParticle;\r
-\r
- while (mcPart)\r
- {\r
- \r
- if(mcPart->IsPrimary())\r
- return mcPart;\r
- \r
- Int_t mother = mcPart->GetMother();\r
-\r
- mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));\r
- nSteps++; // 1 level down\r
- }\r
-\r
- return 0;\r
-}\r
-\r
-\r
-\r
-//__________________________________________________________________\r
-void AliAnalysisTaskQAHighPtDeDx::ProduceArrayTrksESD( AliESDEvent *ESDevent ){\r
- \r
- const Int_t nESDTracks = ESDevent->GetNumberOfTracks();\r
- //Int_t trackmult=0;\r
-\r
-\r
- Int_t multTPC = 0;\r
- \r
- //get multiplicity tpc only track cuts\r
- for(Int_t iT = 0; iT < nESDTracks; iT++) {\r
- \r
- AliESDtrack* esdTrack = ESDevent->GetTrack(iT);\r
- \r
- \r
- if(TMath::Abs(esdTrack->Eta()) > fEtaCut)\r
- continue;\r
- \r
- //only golden track cuts\r
- UInt_t selectDebug = 0;\r
- if (fTrackFilterTPC) {\r
- selectDebug = fTrackFilterTPC->IsSelected(esdTrack);\r
- if (!selectDebug) {\r
- //cout<<"this is not a golden track"<<endl;\r
- continue;\r
- }\r
- }\r
- \r
- multTPC++;\r
- \r
- }\r
-\r
-\r
-\r
- for(Int_t iT = 0; iT < nESDTracks; iT++) {\r
- \r
- AliESDtrack* esdTrack = ESDevent->GetTrack(iT);\r
- \r
- \r
- if(TMath::Abs(esdTrack->Eta()) > fEtaCut)\r
- continue;\r
- \r
- //only golden track cuts\r
- UInt_t selectDebug = 0;\r
- if (fTrackFilterGolden) {\r
- selectDebug = fTrackFilterGolden->IsSelected(esdTrack);\r
- if (!selectDebug) {\r
- //cout<<"this is not a golden track"<<endl;\r
- continue;\r
- }\r
- }\r
- \r
- \r
- //Int_t sharedtpcclusters = esdTrack->GetTPCnclsS();\r
- Short_t ncl = esdTrack->GetTPCsignalN();\r
-\r
-\r
- if(ncl<70)\r
- continue;\r
- Double_t eta = esdTrack->Eta();\r
- Double_t phi = esdTrack->Phi();\r
- Double_t momentum = esdTrack->P();\r
- Float_t dedx = esdTrack->GetTPCsignal();\r
-\r
- //cout<<"magf="<<magf<<endl;\r
-\r
- if(!PhiCut(esdTrack->Pt(), phi, esdTrack->Charge(), magf, cutLow, cutHigh))\r
- continue;\r
-\r
- \r
-\r
-\r
- //TOF\r
- \r
- Bool_t IsTOFout=kFALSE;\r
- if ((esdTrack->GetStatus()&AliESDtrack::kTOFout)==0)\r
- IsTOFout=kTRUE;\r
- Float_t lengthtrack=esdTrack->GetIntegratedLength();\r
- Float_t timeTOF=esdTrack->GetTOFsignal();\r
- Double_t inttime[5]={0,0,0,0,0}; \r
- esdTrack->GetIntegratedTimes(inttime);// Returns the array with integrated times for each particle hypothesis\r
- Float_t beta = -99;\r
- if ( !IsTOFout ){\r
- if ( ( lengthtrack != 0 ) && ( timeTOF != 0) )\r
- beta = inttime[0] / timeTOF;\r
- }\r
-\r
- Short_t pidCode = 0; \r
- \r
- if(fAnalysisMC) {\r
- \r
- const Int_t label = TMath::Abs(esdTrack->GetLabel());\r
- TParticle* mcTrack = fMCStack->Particle(label); \r
- if (mcTrack){\r
- \r
- Int_t pdgCode = mcTrack->GetPdgCode();\r
- pidCode = GetPidCode(pdgCode);\r
- \r
- }\r
- \r
- }\r
-\r
-\r
- if( momentum <= 0.6 && momentum >= 0.4 ){//only p:0.4-0.6 GeV, pion MIP\r
- \r
- if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){ \r
- if(momentum<0.6&&momentum>0.4){\r
- hMIPVsEta->Fill(eta,dedx);\r
- pMIPVsEta->Fill(eta,dedx);\r
- }\r
- }\r
- if( dedx > DeDxMIPMax+1 && dedx < 95 ){\r
- if(TMath::Abs(beta-1)<0.1){\r
- hPlateauVsEta->Fill(eta,dedx);\r
- pPlateauVsEta->Fill(eta,dedx); \r
- }\r
- }\r
- }\r
- \r
- \r
- for(Int_t nh = 0; nh < 9; nh++) {\r
- \r
- if( eta > etaHigh[nh]/10.0 || eta < etaLow[nh]/10.0 )\r
- continue;\r
- \r
- \r
- if(fAnalysisMC){\r
- hMcOut[0][nh]->Fill(esdTrack->Pt());\r
- hMcOut[pidCode][nh]->Fill(esdTrack->Pt());\r
- }\r
-\r
- histAllCh[nh]->Fill(momentum, dedx);\r
- if(beta>1){\r
- histPiTof[nh]->Fill(momentum, dedx);\r
- histpPiTof[nh]->Fill(momentum);\r
- }\r
-\r
- if( momentum <= 0.6 && momentum >= 0.4 ){//only p:0.4-0.6 GeV, pion MIP\r
- //Fill pion MIP, before calibration\r
- if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){ \r
- hMIPVsPhi[nh]->Fill(phi,dedx);\r
- pMIPVsPhi[nh]->Fill(phi,dedx);\r
- \r
- \r
- hMIPVsNch[nh]->Fill(multTPC,dedx);\r
- pMIPVsNch[nh]->Fill(multTPC,dedx);\r
- \r
- }\r
- \r
- //Fill electrons, before calibration\r
- if( dedx > DeDxMIPMax+1 && dedx < 95 ){\r
- if(TMath::Abs(beta-1)<0.1){\r
- hPlateauVsPhi[nh]->Fill(phi,dedx);\r
- pPlateauVsPhi[nh]->Fill(phi,dedx);\r
- }\r
- }\r
-\r
- }\r
-\r
-\r
- }//end loop over eta intervals\r
-\r
- \r
-\r
- \r
- \r
- }//end of track loop\r
- \r
- \r
-\r
- \r
-}\r
-//__________________________________________________________________\r
-void AliAnalysisTaskQAHighPtDeDx::ProduceArrayTrksAOD( AliAODEvent *AODevent ){\r
-\r
-\r
- Int_t nAODTracks = AODevent->GetNumberOfTracks();\r
- Int_t multTPC = 0;\r
- \r
- //get multiplicity tpc only track cuts\r
- for(Int_t iT = 0; iT < nAODTracks; iT++) {\r
- \r
- AliAODTrack* aodTrack = AODevent->GetTrack(iT);\r
- \r
- if(TMath::Abs(aodTrack->Eta()) > fEtaCut)\r
- continue;\r
- \r
-\r
- if (fTrackFilterTPC) {\r
- // TPC only cuts is bit 1\r
- if(!aodTrack->TestFilterBit(1))\r
- continue;\r
- }\r
- \r
- multTPC++;\r
- \r
- }\r
-\r
-\r
- for(Int_t iT = 0; iT < nAODTracks; iT++) {\r
- \r
- AliAODTrack* aodTrack = AODevent->GetTrack(iT);\r
- \r
- if (fTrackFilterGolden) { \r
- // "Global track RAA analysis QM2011 + Chi2ITS<36"; bit 1024\r
- if(!aodTrack->TestFilterBit(1024))\r
- continue;\r
- }\r
- \r
- \r
- if(TMath::Abs(aodTrack->Eta()) > fEtaCut)\r
- continue;\r
- \r
- \r
- Double_t eta = aodTrack->Eta();\r
- Double_t phi = aodTrack->Phi();\r
- Double_t momentum = aodTrack->P();\r
- \r
- \r
- if(!PhiCut(aodTrack->Pt(), phi, aodTrack->Charge(), magf, cutLow, cutHigh))\r
- continue;\r
- \r
-\r
- \r
- AliAODPid* aodPid = aodTrack->GetDetPid();\r
- Short_t ncl = -10;\r
- Float_t dedx = -10;\r
-\r
- //TOF \r
- Float_t beta = -99;\r
-\r
-\r
- if(aodPid) {\r
- ncl = aodPid->GetTPCsignalN();\r
- dedx = aodPid->GetTPCsignal();\r
- //TOF\r
- Bool_t IsTOFout=kFALSE;\r
- Float_t lengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF \r
- Float_t timeTOF=-999;\r
-\r
- if ((aodTrack->GetStatus()&AliESDtrack::kTOFout)==0)\r
- IsTOFout=kTRUE;\r
- \r
- lengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF \r
- \r
- timeTOF=aodPid->GetTOFsignal();\r
- \r
- Double_t inttime[5]={0,0,0,0,0}; \r
- aodPid->GetIntegratedTimes(inttime);// Returns the array with integrated times for each particle hypothesis\r
-\r
- \r
- if ( !IsTOFout ){\r
- if ( ( lengthtrack != 0 ) && ( timeTOF != 0) )\r
- beta = inttime[0] / timeTOF;\r
- }\r
- \r
- }\r
-\r
-\r
- if(ncl<70)\r
- continue;\r
-\r
-\r
- Short_t pidCode = 0; \r
- \r
- if(fAnalysisMC) {\r
- \r
- const Int_t label = TMath::Abs(aodTrack->GetLabel());\r
- AliAODMCParticle* mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(label));\r
- \r
- if (mcTrack){\r
- \r
- Int_t pdgCode = mcTrack->GetPdgCode();\r
- pidCode = GetPidCode(pdgCode);\r
- \r
- }\r
- \r
- }\r
-\r
-\r
-\r
-\r
- //if(momentum>0.6||momentum<0.4)//only p:0.4-0.6 GeV\r
- //continue;\r
- \r
- //etaLow\r
- //etaHigh\r
- if( momentum <= 0.6 && momentum >= 0.4 ){//only p:0.4-0.6 GeV, pion MIP\r
- if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){ \r
- if(momentum<0.6&&momentum>0.4){\r
- hMIPVsEta->Fill(eta,dedx);\r
- pMIPVsEta->Fill(eta,dedx);\r
- }\r
- }\r
- if( dedx > DeDxMIPMax+1 && dedx < 95 ){\r
- if(TMath::Abs(beta-1)<0.1){\r
- hPlateauVsEta->Fill(eta,dedx);\r
- pPlateauVsEta->Fill(eta,dedx); \r
- }\r
- }\r
- }\r
- \r
- \r
- for(Int_t nh = 0; nh < 9; nh++) {\r
- \r
- if( eta > etaHigh[nh]/10.0 || eta < etaLow[nh]/10.0 )\r
- continue;\r
- \r
- \r
- if(fAnalysisMC){\r
- hMcOut[0][nh]->Fill(aodTrack->Pt());\r
- hMcOut[pidCode][nh]->Fill(aodTrack->Pt());\r
- }\r
- \r
- histAllCh[nh]->Fill(momentum, dedx);\r
- if(beta>1){\r
- histPiTof[nh]->Fill(momentum, dedx);\r
- histpPiTof[nh]->Fill(momentum);\r
- }\r
-\r
- if( momentum <= 0.6 && momentum >= 0.4 ){//only p:0.4-0.6 GeV, pion MIP\r
- //Fill pion MIP, before calibration\r
- if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){ \r
- hMIPVsPhi[nh]->Fill(phi,dedx);\r
- pMIPVsPhi[nh]->Fill(phi,dedx);\r
- \r
- \r
- hMIPVsNch[nh]->Fill(multTPC,dedx);\r
- pMIPVsNch[nh]->Fill(multTPC,dedx);\r
- \r
- }\r
- \r
- //Fill electrons, before calibration\r
- if( dedx > DeDxMIPMax+1 && dedx < 95 ){\r
- if(TMath::Abs(beta-1)<0.1){\r
- hPlateauVsPhi[nh]->Fill(phi,dedx);\r
- pPlateauVsPhi[nh]->Fill(phi,dedx);\r
- }\r
- }\r
-\r
- }\r
-\r
- }//end loop over eta intervals\r
- \r
- \r
-\r
- \r
- \r
- }//end of track loop\r
- \r
- \r
- \r
- \r
- \r
- \r
- \r
-}\r
-//__________________________________________________\r
-void AliAnalysisTaskQAHighPtDeDx::ProduceArrayV0ESD( AliESDEvent *ESDevent ){\r
- Int_t nv0s = ESDevent->GetNumberOfV0s();\r
- /*\r
- if(nv0s<1){\r
- return;\r
- }*/\r
- \r
- const AliESDVertex *myBestPrimaryVertex = ESDevent->GetPrimaryVertex();\r
- if (!myBestPrimaryVertex) return;\r
- if (!(myBestPrimaryVertex->GetStatus())) return;\r
- Double_t lPrimaryVtxPosition[3];\r
- myBestPrimaryVertex->GetXYZ(lPrimaryVtxPosition);\r
- \r
- Double_t lPrimaryVtxCov[6];\r
- myBestPrimaryVertex->GetCovMatrix(lPrimaryVtxCov);\r
- Double_t lPrimaryVtxChi2 = myBestPrimaryVertex->GetChi2toNDF();\r
- \r
- AliAODVertex* myPrimaryVertex = new AliAODVertex(lPrimaryVtxPosition, lPrimaryVtxCov, lPrimaryVtxChi2, NULL, -1, AliAODVertex::kPrimary);\r
-\r
-\r
- \r
- //\r
- // LOOP OVER V0s, K0s, L, AL\r
- //\r
- \r
- \r
- for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {\r
- \r
- // This is the begining of the V0 loop \r
- AliESDv0 *esdV0 = ESDevent->GetV0(iV0);\r
- if (!esdV0) continue;\r
- \r
- //check onfly status\r
- if(esdV0->GetOnFlyStatus()!=0)\r
- continue;\r
-\r
-\r
- // AliESDTrack (V0 Daughters)\r
- UInt_t lKeyPos = (UInt_t)TMath::Abs(esdV0->GetPindex());\r
- UInt_t lKeyNeg = (UInt_t)TMath::Abs(esdV0->GetNindex());\r
- \r
- AliESDtrack *pTrack = ESDevent->GetTrack(lKeyPos);\r
- AliESDtrack *nTrack = ESDevent->GetTrack(lKeyNeg);\r
- if (!pTrack || !nTrack) {\r
- Printf("ERROR: Could not retreive one of the daughter track");\r
- continue;\r
- }\r
- \r
- // Remove like-sign\r
- if (pTrack->GetSign() == nTrack->GetSign()) {\r
- //cout<< "like sign, continue"<< endl;\r
- continue;\r
- } \r
- \r
- // Eta cut on decay products\r
- if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)\r
- continue;\r
- \r
- \r
- // Check if switch does anything!\r
- Bool_t isSwitched = kFALSE;\r
- if (pTrack->GetSign() < 0) { // switch\r
- \r
- isSwitched = kTRUE;\r
- AliESDtrack* helpTrack = nTrack;\r
- nTrack = pTrack;\r
- pTrack = helpTrack;\r
- } \r
- \r
- //Double_t lV0Position[3];\r
- //esdV0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[2]);\r
- \r
-\r
- AliKFVertex primaryVtxKF( *myPrimaryVertex );\r
- AliKFParticle::SetField(ESDevent->GetMagneticField());\r
- \r
- // Also implement switch here!!!!!!\r
- AliKFParticle* negEKF = 0; // e-\r
- AliKFParticle* posEKF = 0; // e+\r
- AliKFParticle* negPiKF = 0; // pi -\r
- AliKFParticle* posPiKF = 0; // pi +\r
- AliKFParticle* posPKF = 0; // p\r
- AliKFParticle* negAPKF = 0; // p-bar\r
- \r
- if(!isSwitched) {\r
- negEKF = new AliKFParticle( *(esdV0->GetParamN()) , 11);\r
- posEKF = new AliKFParticle( *(esdV0->GetParamP()) ,-11);\r
- negPiKF = new AliKFParticle( *(esdV0->GetParamN()) ,-211);\r
- posPiKF = new AliKFParticle( *(esdV0->GetParamP()) , 211);\r
- posPKF = new AliKFParticle( *(esdV0->GetParamP()) , 2212);\r
- negAPKF = new AliKFParticle( *(esdV0->GetParamN()) ,-2212);\r
- } else { // switch + and - \r
- negEKF = new AliKFParticle( *(esdV0->GetParamP()) , 11);\r
- posEKF = new AliKFParticle( *(esdV0->GetParamN()) ,-11);\r
- negPiKF = new AliKFParticle( *(esdV0->GetParamP()) ,-211);\r
- posPiKF = new AliKFParticle( *(esdV0->GetParamN()) , 211);\r
- posPKF = new AliKFParticle( *(esdV0->GetParamN()) , 2212);\r
- negAPKF = new AliKFParticle( *(esdV0->GetParamP()) ,-2212);\r
- }\r
- \r
- AliKFParticle v0GKF; // Gamma e.g. from pi0\r
- v0GKF+=(*negEKF);\r
- v0GKF+=(*posEKF);\r
- v0GKF.SetProductionVertex(primaryVtxKF);\r
- \r
- AliKFParticle v0K0sKF; // K0 short\r
- v0K0sKF+=(*negPiKF);\r
- v0K0sKF+=(*posPiKF);\r
- v0K0sKF.SetProductionVertex(primaryVtxKF);\r
- \r
- AliKFParticle v0LambdaKF; // Lambda\r
- v0LambdaKF+=(*negPiKF);\r
- v0LambdaKF+=(*posPKF); \r
- v0LambdaKF.SetProductionVertex(primaryVtxKF);\r
- \r
- AliKFParticle v0AntiLambdaKF; // Lambda-bar\r
- v0AntiLambdaKF+=(*posPiKF);\r
- v0AntiLambdaKF+=(*negAPKF);\r
- v0AntiLambdaKF.SetProductionVertex(primaryVtxKF);\r
- \r
- Double_t dmassG = v0GKF.GetMass();\r
- Double_t dmassK = v0K0sKF.GetMass()-0.498;\r
- Double_t dmassL = v0LambdaKF.GetMass()-1.116;\r
- Double_t dmassAL = v0AntiLambdaKF.GetMass()-1.116;\r
-\r
-\r
- for( Int_t case_v0 = 0; case_v0 < 2; ++case_v0 ){\r
-\r
-\r
- switch(case_v0){\r
- case 0:{ \r
- \r
- Bool_t fillPos = kFALSE;\r
- Bool_t fillNeg = kFALSE;\r
- \r
- if(dmassG < 0.1)\r
- continue;\r
- \r
- if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01){\r
- continue;\r
- }\r
- \r
- if(dmassL<0.01){\r
- fillPos = kTRUE;\r
- }\r
- if(dmassAL<0.01) {\r
- if(fillPos)\r
- continue;\r
- fillNeg = kTRUE;\r
- }\r
- \r
- if(dmassK<0.01) {\r
- if(fillPos||fillNeg)\r
- continue;\r
- fillPos = kTRUE;\r
- fillNeg = kTRUE;\r
- }\r
- \r
- \r
- for(Int_t j = 0; j < 2; j++) {\r
- \r
- AliESDtrack* track = 0;\r
- \r
- if(j==0) {\r
- \r
- if(fillNeg)\r
- track = nTrack;\r
- else\r
- continue;\r
- } else {\r
- \r
- if(fillPos)\r
- track = pTrack;\r
- else\r
- continue;\r
- }\r
- \r
- if(track->GetTPCsignalN()<=70)continue;\r
- Double_t phi = track->Phi();\r
- \r
- if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))\r
- continue;\r
- \r
- \r
- //if(!PhiCut(pt, phi, charge, magf, cutLow, cutHigh, hPhi))\r
- // continue;\r
- \r
- Double_t eta = track->Eta();\r
- Double_t momentum = track->Pt();\r
- Double_t dedx = track->GetTPCsignal();\r
- \r
- if(fillPos&&fillNeg){\r
- \r
- \r
- if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){ \r
- if(momentum<0.6&&momentum>0.4){\r
- hMIPVsEtaV0s->Fill(eta,dedx);\r
- pMIPVsEtaV0s->Fill(eta,dedx);\r
- }\r
- }\r
- \r
- \r
- }\r
- \r
- for(Int_t nh = 0; nh < nHists; nh++) {\r
- \r
- \r
- \r
- if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )\r
- continue;\r
- \r
- if(fillPos&&fillNeg){\r
- \r
- histPiV0[nh]->Fill(momentum, dedx); \r
- histpPiV0[nh]->Fill(momentum);\r
-\r
- }\r
- else{\r
- \r
- histPV0[nh]->Fill(momentum, dedx);\r
- histpPV0[nh]->Fill(momentum);\r
-\r
- }\r
- \r
- }\r
- \r
- }//end loop over two tracks\r
-\r
- };\r
- break;\r
-\r
- case 1:{//gammas \r
- \r
- Bool_t fillPos = kFALSE;\r
- Bool_t fillNeg = kFALSE;\r
- \r
-\r
-\r
- if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01) {\r
- if(dmassG<0.01 && dmassG>0.0001) {\r
- \r
-\r
- if( TMath::Abs(nTrack->GetTPCsignal() - 85.0) < 5)\r
- fillPos = kTRUE;\r
- if( TMath::Abs(pTrack->GetTPCsignal() - 85.0) < 5)\r
- fillNeg = kTRUE;\r
- \r
- } else {\r
- continue;\r
- }\r
- }\r
-\r
- //cout<<"fillPos="<<fillPos<<endl;\r
-\r
- if(fillPos == kTRUE && fillNeg == kTRUE) \r
- continue;\r
- \r
- \r
- AliESDtrack* track = 0;\r
- if(fillNeg)\r
- track = nTrack;\r
- else if(fillPos)\r
- track = pTrack;\r
- else\r
- continue;\r
-\r
- Double_t dedx = track->GetTPCsignal();\r
- Double_t eta = track->Eta();\r
- Double_t phi = track->Phi();\r
- Double_t momentum = track->P();\r
-\r
- if(track->GetTPCsignalN()<=70)continue;\r
-\r
- if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))\r
- continue;\r
- \r
- for(Int_t nh = 0; nh < nHists; nh++) {\r
- \r
- if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )\r
- continue;\r
-\r
- histEV0[nh]->Fill(momentum, dedx);\r
-\r
- }\r
- \r
- };\r
- break;\r
- \r
- \r
- }//end switch\r
- \r
- }//end loop over case V0\r
-\r
- \r
- // clean up loop over v0\r
- \r
- delete negPiKF;\r
- delete posPiKF;\r
- delete posPKF;\r
- delete negAPKF;\r
- \r
- \r
- \r
- }\r
- \r
- \r
- delete myPrimaryVertex;\r
- \r
- \r
-}\r
-//__________________________________________________________________________\r
-void AliAnalysisTaskQAHighPtDeDx::ProduceArrayV0AOD( AliAODEvent *AODevent ){\r
- Int_t nv0s = AODevent->GetNumberOfV0s();\r
- /*\r
- if(nv0s<1){\r
- return;\r
- }*/\r
- \r
- AliAODVertex *myBestPrimaryVertex = AODevent->GetPrimaryVertex();\r
- if (!myBestPrimaryVertex) return;\r
- \r
-\r
- \r
- // ################################\r
- // #### BEGINNING OF V0 CODE ######\r
- // ################################\r
- // This is the begining of the V0 loop \r
- for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {\r
- AliAODv0 *aodV0 = AODevent->GetV0(iV0);\r
- if (!aodV0) continue;\r
- \r
- \r
- //check onfly status\r
- if(aodV0->GetOnFlyStatus()!=0)\r
- continue;\r
-\r
- // AliAODTrack (V0 Daughters)\r
- AliAODVertex* vertex = aodV0->GetSecondaryVtx();\r
- if (!vertex) {\r
- Printf("ERROR: Could not retrieve vertex");\r
- continue;\r
- }\r
- \r
- AliAODTrack *pTrack = (AliAODTrack*)vertex->GetDaughter(0);\r
- AliAODTrack *nTrack = (AliAODTrack*)vertex->GetDaughter(1);\r
- if (!pTrack || !nTrack) {\r
- Printf("ERROR: Could not retrieve one of the daughter track");\r
- continue;\r
- }\r
- \r
- // Remove like-sign\r
- if (pTrack->Charge() == nTrack->Charge()) {\r
- //cout<< "like sign, continue"<< endl;\r
- continue;\r
- } \r
- \r
- // Make sure charge ordering is ok\r
- if (pTrack->Charge() < 0) {\r
- AliAODTrack* helpTrack = pTrack;\r
- pTrack = nTrack;\r
- nTrack = helpTrack;\r
- } \r
- \r
- // Eta cut on decay products\r
- if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)\r
- continue;\r
- \r
- \r
- Double_t dmassG = aodV0->InvMass2Prongs(0,1,11,11);\r
- Double_t dmassK = aodV0->MassK0Short()-0.498;\r
- Double_t dmassL = aodV0->MassLambda()-1.116;\r
- Double_t dmassAL = aodV0->MassAntiLambda()-1.116;\r
- \r
- for( Int_t case_v0 = 0; case_v0 < 2; ++case_v0 ){\r
- \r
- \r
- switch(case_v0){\r
- case 0:{ \r
- Bool_t fillPos = kFALSE;\r
- Bool_t fillNeg = kFALSE;\r
- \r
- \r
- if(dmassG < 0.1)\r
- continue;\r
- \r
- if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01){\r
- continue;\r
- }\r
- \r
- if(dmassL<0.01){\r
- fillPos = kTRUE;\r
- }\r
- if(dmassAL<0.01) {\r
- if(fillPos)\r
- continue;\r
- fillNeg = kTRUE;\r
- }\r
- \r
- if(dmassK<0.01) {\r
- if(fillPos||fillNeg)\r
- continue;\r
- fillPos = kTRUE;\r
- fillNeg = kTRUE;\r
- }\r
- \r
- \r
- \r
- for(Int_t j = 0; j < 2; j++) {\r
- \r
- AliAODTrack* track = 0;\r
- \r
- if(j==0) {\r
- \r
- if(fillNeg)\r
- track = nTrack;\r
- else\r
- continue;\r
- } else {\r
- \r
- if(fillPos)\r
- track = pTrack;\r
- else\r
- continue;\r
- }\r
- \r
- if(track->GetTPCsignalN()<=70)continue;\r
- \r
- Double_t phi = track->Phi();\r
- \r
- if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))\r
- continue;\r
- \r
- \r
- //if(!PhiCut(pt, phi, charge, magf, cutLow, cutHigh, hPhi))\r
- // continue;\r
- \r
- Double_t eta = track->Eta();\r
- Double_t momentum = track->Pt();\r
- Double_t dedx = track->GetTPCsignal();\r
- \r
- if(fillPos&&fillNeg){\r
- \r
- \r
- if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){ \r
- if(momentum<0.6&&momentum>0.4){\r
- hMIPVsEtaV0s->Fill(eta,dedx);\r
- pMIPVsEtaV0s->Fill(eta,dedx);\r
- }\r
- }\r
- \r
- \r
- }\r
- \r
- for(Int_t nh = 0; nh < nHists; nh++) {\r
- \r
- \r
- \r
- if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )\r
- continue;\r
- \r
- if(fillPos&&fillNeg){\r
- \r
- histPiV0[nh]->Fill(momentum, dedx); \r
- histpPiV0[nh]->Fill(momentum); \r
-\r
- }\r
- else{\r
- \r
- histPV0[nh]->Fill(momentum, dedx);\r
- histpPV0[nh]->Fill(momentum);\r
-\r
- }\r
- \r
- }\r
- \r
-\r
- }//end loop over two tracks\r
- };\r
- break;\r
- \r
- case 1:{//gammas \r
- \r
- Bool_t fillPos = kFALSE;\r
- Bool_t fillNeg = kFALSE;\r
- \r
- if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01) {\r
- if(dmassG<0.01 && dmassG>0.0001) {\r
- \r
- if( TMath::Abs(nTrack->GetTPCsignal() - 85.0) < 5)\r
- fillPos = kTRUE;\r
- if( TMath::Abs(pTrack->GetTPCsignal() - 85.0) < 5)\r
- fillNeg = kTRUE;\r
- \r
- } else {\r
- continue;\r
- }\r
- }\r
-\r
- \r
- if(fillPos == kTRUE && fillNeg == kTRUE) \r
- continue;\r
- \r
- \r
- AliAODTrack* track = 0;\r
- if(fillNeg)\r
- track = nTrack;\r
- else if(fillPos)\r
- track = pTrack;\r
- else\r
- continue;\r
-\r
- Double_t dedx = track->GetTPCsignal();\r
- Double_t eta = track->Eta();\r
- Double_t phi = track->Phi();\r
- Double_t momentum = track->P();\r
-\r
- if(track->GetTPCsignalN()<=70)continue;\r
-\r
- if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))\r
- continue;\r
- \r
- for(Int_t nh = 0; nh < nHists; nh++) {\r
- \r
- if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )\r
- continue;\r
-\r
- histEV0[nh]->Fill(momentum, dedx);\r
-\r
- }\r
- \r
- };\r
- break;\r
-\r
-\r
- }//end switch\r
- }//end loop over V0s cases\r
- \r
- }//end loop over v0's\r
- \r
- \r
- \r
- \r
-}\r
-Bool_t AliAnalysisTaskQAHighPtDeDx::PhiCut(Double_t pt, Double_t phi, Double_t q, Float_t mag, TF1* phiCutLow, TF1* phiCutHigh)\r
-{\r
- if(pt < 2.0)\r
- return kTRUE;\r
- \r
- //Double_t phi = track->Phi();\r
- if(mag < 0) // for negatve polarity field\r
- phi = TMath::TwoPi() - phi;\r
- if(q < 0) // for negatve charge\r
- phi = TMath::TwoPi()-phi;\r
- \r
- phi += TMath::Pi()/18.0; // to center gap in the middle\r
- phi = fmod(phi, TMath::Pi()/9.0);\r
- \r
- if(phi<phiCutHigh->Eval(pt) \r
- && phi>phiCutLow->Eval(pt))\r
- return kFALSE; // reject track\r
-\r
- hPhi->Fill(pt, phi);\r
-\r
- return kTRUE;\r
-}\r
+ /*************************************************************************
+* 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;
+Float_t centralityGlobal = -10;
+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 = 80;
+Double_t deltaPiLow = 20;
+Double_t deltaPiHigh = 100;
+const Char_t *Pid[7]={"Ch","Pion","Kaon","Proton","Electron","Muon","Oher"};
+ClassImp(AliAnalysisTaskQAHighPtDeDx)
+//_____________________________________________________________________________
+//AliAnalysisTaskQAHighPtDeDx::AliAnalysisTaskQAHighPtDeDx(const char *name):
+AliAnalysisTaskQAHighPtDeDx::AliAnalysisTaskQAHighPtDeDx():
+ AliAnalysisTaskSE(),
+ 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
+}
+
+
+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
+ histpPiV0[i]=0;//TH1D, pi id by V0s
+ histPV0[i]=0;// TH2D, dE/dx vs p, p id by V0s
+ histpPV0[i]=0;// TH1D, 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
+ histpPiTof[i]=0;//TH1D, 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
+}
+
+
+
+
+AliAnalysisTaskQAHighPtDeDx::~AliAnalysisTaskQAHighPtDeDx() {
+ //
+ // Destructor
+ //
+
+}
+
+
+
+
+
+
+
+
+//______________________________________________________________________________
+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]);
+
+ //two dimmesional distributions dE/dx vs p for secondary pions
+ histPiV0[i] = new TH2D(Form("histPiV0%s", ending[i]), "Pions id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
+ histPiV0[i]->Sumw2();
+ fListOfObjects->Add(histPiV0[i]);
+
+ histpPiV0[i] = new TH1D(Form("histPiV01D%s", ending[i]), "Pions id by V0; #it{p} (GeV/#it{c}); counts", 200, 0, 20);
+ histpPiV0[i]->Sumw2();
+ fListOfObjects->Add(histpPiV0[i]);
+
+ //two dimmesional distributions dE/dx vs p for secondary protons
+ histPV0[i] = new TH2D(Form("histPV0%s", ending[i]), "Protons id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
+ histPV0[i]->Sumw2();
+ fListOfObjects->Add(histPV0[i]);
+
+ histpPV0[i] = new TH1D(Form("histPV01D%s", ending[i]), "Protons id by V0; #it{p} (GeV/#it{c}); counts", 200, 0, 20);
+ histpPV0[i]->Sumw2();
+ fListOfObjects->Add(histpPV0[i]);
+
+ //two dimmesional distributions dE/dx vs p for primary pions
+ histPiTof[i] = new TH2D(Form("histPiTof%s", ending[i]), "all charged", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
+ histPiTof[i]->Sumw2();
+
+ histpPiTof[i] = new TH1D(Form("histPiTof1D%s", ending[i]), "Protons id by V0; #it{p} (GeV/#it{c}); counts", 200, 0, 20);
+ histpPiTof[i]->Sumw2();
+ fListOfObjects->Add(histpPiTof[i]);
+
+
+ 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);
+ centralityGlobal = centrality;
+ 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());
+
+ cout << "centrality=" << centralityGlobal << endl;
+
+
+ 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));
+
+ if(!trackMC){
+ AliError("Cannot get MC particle");
+ continue;
+ }
+
+ //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+1 && 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);
+ histpPiTof[nh]->Fill(momentum);
+ }
+
+ 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+1 && 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+1 && 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);
+ histpPiTof[nh]->Fill(momentum);
+ }
+
+ 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+1 && 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);
+ histpPiV0[nh]->Fill(momentum);
+
+ }
+ else{
+
+ histPV0[nh]->Fill(momentum, dedx);
+ histpPV0[nh]->Fill(momentum);
+
+ }
+
+ }
+
+ }//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);
+ histpPiV0[nh]->Fill(momentum);
+
+ }
+ else{
+
+ histPV0[nh]->Fill(momentum, dedx);
+ histpPV0[nh]->Fill(momentum);
+
+ }
+
+ }
+
+
+ }//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;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff