1 /*************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. * **************************************************************************/
16 #include "AliAnalysisTaskQAHighPtDeDx.h"
25 #include <TParticle.h>
29 #include <AliAnalysisManager.h>
30 #include <AliAnalysisFilter.h>
31 #include <AliESDInputHandler.h>
32 #include <AliESDEvent.h>
33 #include <AliESDVertex.h>
35 #include <AliExternalTrackParam.h>
36 #include <AliESDtrackCuts.h>
37 #include <AliESDVZERO.h>
38 #include <AliAODVZERO.h>
40 #include <AliMCEventHandler.h>
41 #include <AliMCEvent.h>
44 #include <TTreeStream.h>
46 #include <AliHeader.h>
47 #include <AliGenPythiaEventHeader.h>
48 #include <AliGenDPMjetEventHeader.h>
50 #include "AliCentrality.h"
52 #include <AliKFVertex.h>
53 #include <AliAODVertex.h>
55 #include <AliAODTrack.h>
56 #include <AliAODPid.h>
57 #include <AliAODMCHeader.h>
67 // Antonio Ortiz (Lund)
68 // Peter Christiansen (Lund)
78 const Double_t AliAnalysisTaskQAHighPtDeDx::fgkClight = 2.99792458e-2;
80 TF1* cutLow = new TF1("StandardPhiCutLow", "0.1/x/x+pi/18.0-0.025", 0, 50);
81 TF1* cutHigh = new TF1("StandardPhiCutHigh", "0.12/x+pi/18.0+0.035", 0, 50);
82 Double_t DeDxMIPMin = 30;
83 Double_t DeDxMIPMax = 65;
84 const Int_t nHists = 9;
85 Float_t centralityGlobal = -10;
86 Int_t etaLow[nHists] = {-8, -8, -6, -4, -2, 0, 2, 4, 6};
87 Int_t etaHigh[nHists] = { 8, -6, -4, -2, 0, 2, 4, 6, 8};
89 Int_t nDeltaPiBins = 80;
90 Double_t deltaPiLow = 20;
91 Double_t deltaPiHigh = 100;
92 const Char_t *Pid[7]={"Ch","Pion","Kaon","Proton","Electron","Muon","Oher"};
93 ClassImp(AliAnalysisTaskQAHighPtDeDx)
94 //_____________________________________________________________________________
95 //AliAnalysisTaskQAHighPtDeDx::AliAnalysisTaskQAHighPtDeDx(const char *name):
96 AliAnalysisTaskQAHighPtDeDx::AliAnalysisTaskQAHighPtDeDx():
103 fTrackFilterGolden(0x0),
104 fTrackFilterTPC(0x0),
106 fAnalysisType("ESD"),
108 fAnalysisPbPb(kFALSE),
116 fStoreMcIn(kFALSE),//
117 fMcProcessType(-999),
118 fTriggeredEventMB(-999),
125 fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
138 //default constructor
139 for(Int_t i=0;i<9;++i){
141 hMIPVsNch[i]=0;//TH2D, MIP vs Nch for different eta intervals
142 pMIPVsNch[i]=0;//TProfile, MIP vs Nch for different eta intervals
143 hMIPVsPhi[i]=0;//TH2D, MIP vs phi for different eta intervals
144 pMIPVsPhi[i]=0;//TProfile, MIP vs phi for different eta intervals
145 hPlateauVsPhi[i]=0;//TH2D, dE/dx vs Phi, electrons 0.4<p<0.6 GeV/c
146 pPlateauVsPhi[i]=0;//TProfile, dE/dx vs Phi, electrons 0.4<p<0.6 GeV/c
147 histPiV0[i]=0;//TH2D, dE/dx vs p, pi id by V0s
148 histpPiV0[i]=0;//TH1D, pi id by V0s
149 histPV0[i]=0;// TH2D, dE/dx vs p, p id by V0s
150 histpPV0[i]=0;// TH1D, p id by V0s
151 histAllCh[i]=0;//TH2D, dE/dx vs p for all charged particles
152 histPiTof[i]=0;//TH2D, dE/dx vs p for a "clean" sample of pions, beta>1
153 histpPiTof[i]=0;//TH1D, for a "clean" sample of pions, beta>1
156 for(Int_t pid=0;pid<7;++pid){
168 AliAnalysisTaskQAHighPtDeDx::AliAnalysisTaskQAHighPtDeDx(const char *name):
169 AliAnalysisTaskSE(name),
175 fTrackFilterGolden(0x0),
176 fTrackFilterTPC(0x0),
178 fAnalysisType("ESD"),
180 fAnalysisPbPb(kFALSE),
188 fStoreMcIn(kFALSE),//
189 fMcProcessType(-999),
190 fTriggeredEventMB(-999),
197 fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
210 // Default constructor (should not be used)
211 for(Int_t i=0;i<9;++i){
213 hMIPVsNch[i]=0;//TH2D, MIP vs Nch for different eta intervals
214 pMIPVsNch[i]=0;//TProfile, MIP vs Nch for different eta intervals
215 hMIPVsPhi[i]=0;//TH2D, MIP vs phi for different eta intervals
216 pMIPVsPhi[i]=0;//TProfile, MIP vs phi for different eta intervals
217 hPlateauVsPhi[i]=0;//TH2D, dE/dx vs Phi, electrons 0.4<p<0.6 GeV/c
218 pPlateauVsPhi[i]=0;//TProfile, dE/dx vs Phi, electrons 0.4<p<0.6 GeV/c
219 histPiV0[i]=0;//TH2D, dE/dx vs p, pi id by V0s
220 histpPiV0[i]=0;//TH1D, pi id by V0s
221 histPV0[i]=0;// TH2D, dE/dx vs p, p id by V0s
222 histpPV0[i]=0;// TH1D, p id by V0s
223 histAllCh[i]=0;//TH2D, dE/dx vs p for all charged particles
224 histPiTof[i]=0;//TH2D, dE/dx vs p for a "clean" sample of pions, beta>1
225 histpPiTof[i]=0;//TH1D, for a "clean" sample of pions, beta>1
228 for(Int_t pid=0;pid<7;++pid){
234 DefineOutput(1, TList::Class());//esto es nuevo
240 AliAnalysisTaskQAHighPtDeDx::~AliAnalysisTaskQAHighPtDeDx() {
254 //______________________________________________________________________________
255 void AliAnalysisTaskQAHighPtDeDx::UserCreateOutputObjects()
257 // This method is called once per worker node
258 // Here we define the output: histograms and debug tree if requested
259 // We also create the random generator here so it might get different seeds...
260 fRandom = new TRandom(0); // 0 means random seed
265 fListOfObjects = new TList();
266 fListOfObjects->SetOwner();
274 fEvents = new TH1I("fEvents","Number of analyzed events; Events; Counts", 3, 0, 3);
275 fListOfObjects->Add(fEvents);
277 fn1=new TH1F("fn1","fn1",11,-1,10);
278 fListOfObjects->Add(fn1);
280 fcent=new TH1F("fcent","fcent",104,-2,102);
281 fListOfObjects->Add(fcent);
283 fVtx = new TH1I("fVtx","Vtx info (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
284 fListOfObjects->Add(fVtx);
286 fVtxBeforeCuts = new TH1F("fVtxBeforeCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
287 fListOfObjects->Add(fVtxBeforeCuts);
289 fVtxAfterCuts = new TH1F("fVtxAfterCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
290 fListOfObjects->Add(fVtxAfterCuts);
293 const Int_t nPtBinsV0s = 25;
294 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 ,
295 1.2 , 1.4 , 1.6 , 1.8 , 2.0 , 2.5 , 3.0 , 3.5 , 4.0 , 5.0 , 7.0 ,
296 9.0 , 12.0, 15.0, 20.0 };
301 const Char_t* ending[nHists] = {"", "86", "64", "42", "20", "02", "24", "46", "68"};
303 const Char_t* LatexEta[nHists] = {
304 "|#eta|<0.8", "-0.8<#eta<-0.6", "-0.6<#eta<-0.4", "-0.4<#eta<-0.2", "-0.2<#eta<0",
305 "0<#eta<0.2", "0.2<#eta<0.4", "0.4<#eta<0.6", "0.6<#eta<0.8"
307 hPhi = new TH2D("histPhi", "pt; #phi'", nPtBinsV0s, ptBinsV0s, 90, -0.05, 0.4);
308 //dE/dx vs phi, pions at the MIP
309 fListOfObjects->Add(hPhi);
316 for(Int_t i = 0; i < nHists; i++) {
319 hMIPVsPhi[i] = new TH2D(Form("hMIPVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx MIP",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(),
320 DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
321 hMIPVsPhi[i]->Sumw2();
323 pMIPVsPhi[i] = new TProfile(Form("pMIPVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx MIP",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(),
324 DeDxMIPMin, DeDxMIPMax);
325 pMIPVsPhi[i]->SetMarkerStyle(20);
326 pMIPVsPhi[i]->SetMarkerColor(1);
327 pMIPVsPhi[i]->SetLineColor(1);
328 pMIPVsPhi[i]->Sumw2();
330 fListOfObjects->Add(hMIPVsPhi[i]);
331 fListOfObjects->Add(pMIPVsPhi[i]);
333 hPlateauVsPhi[i] = new TH2D(Form("hPlateauVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx Plateau",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(),
334 95-DeDxMIPMax, DeDxMIPMax, 95);
335 hPlateauVsPhi[i]->Sumw2();
337 pPlateauVsPhi[i] = new TProfile(Form("pPlateauVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx Plateau",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(),
339 pPlateauVsPhi[i]->SetMarkerStyle(20);
340 pPlateauVsPhi[i]->SetMarkerColor(1);
341 pPlateauVsPhi[i]->SetLineColor(1);
342 pPlateauVsPhi[i]->Sumw2();
344 fListOfObjects->Add(hPlateauVsPhi[i]);
345 fListOfObjects->Add(pPlateauVsPhi[i]);
348 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,
349 DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
350 hMIPVsNch[i]->Sumw2();
352 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);
353 pMIPVsNch[i]->SetMarkerStyle(20);
354 pMIPVsNch[i]->SetMarkerColor(1);
355 pMIPVsNch[i]->SetLineColor(1);
356 pMIPVsNch[i]->Sumw2();
358 fListOfObjects->Add(hMIPVsNch[i]);
359 fListOfObjects->Add(pMIPVsNch[i]);
361 //two dimmesional distributions dE/dx vs p for secondary pions
362 histPiV0[i] = new TH2D(Form("histPiV0%s", ending[i]), "Pions id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
363 histPiV0[i]->Sumw2();
364 fListOfObjects->Add(histPiV0[i]);
366 histpPiV0[i] = new TH1D(Form("histPiV01D%s", ending[i]), "Pions id by V0; #it{p} (GeV/#it{c}); counts", 200, 0, 20);
367 histpPiV0[i]->Sumw2();
368 fListOfObjects->Add(histpPiV0[i]);
370 //two dimmesional distributions dE/dx vs p for secondary protons
371 histPV0[i] = new TH2D(Form("histPV0%s", ending[i]), "Protons id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
373 fListOfObjects->Add(histPV0[i]);
375 histpPV0[i] = new TH1D(Form("histPV01D%s", ending[i]), "Protons id by V0; #it{p} (GeV/#it{c}); counts", 200, 0, 20);
376 histpPV0[i]->Sumw2();
377 fListOfObjects->Add(histpPV0[i]);
379 //two dimmesional distributions dE/dx vs p for primary pions
380 histPiTof[i] = new TH2D(Form("histPiTof%s", ending[i]), "all charged", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
381 histPiTof[i]->Sumw2();
383 histpPiTof[i] = new TH1D(Form("histPiTof1D%s", ending[i]), "Protons id by V0; #it{p} (GeV/#it{c}); counts", 200, 0, 20);
384 histpPiTof[i]->Sumw2();
385 fListOfObjects->Add(histpPiTof[i]);
388 histAllCh[i] = new TH2D(Form("histAllCh%s", ending[i]), "Pions id by TOF", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
389 histAllCh[i]->Sumw2();
391 fListOfObjects->Add(histPiTof[i]);
392 fListOfObjects->Add(histAllCh[i]);
394 histEV0[i] = new TH2D(Form("histEV0%s", ending[i]), "Electrons id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
396 fListOfObjects->Add(histEV0[i]);
403 hMIPVsEta = new TH2D("hMIPVsEta","; #eta; dE/dx_{MIP, primary tracks}",16,-0.8,0.8,DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
404 pMIPVsEta = new TProfile("pMIPVsEta","; #eta; #LT dE/dx #GT_{MIP, primary tracks}",16,-0.8,0.8, DeDxMIPMin, DeDxMIPMax);
405 hMIPVsEtaV0s = new TH2D("hMIPVsEtaV0s","; #eta; dE/dx_{MIP, secondary tracks}",16,-0.8,0.8,DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
406 pMIPVsEtaV0s = new TProfile("pMIPVsEtaV0s","; #eta; #LT dE/dx #GT_{MIP, secondary tracks}",16,-0.8,0.8,DeDxMIPMin, DeDxMIPMax);
408 hPlateauVsEta = new TH2D("hPlateauVsEta","; #eta; dE/dx_{Plateau, primary tracks}",16,-0.8,0.8,95-DeDxMIPMax, DeDxMIPMax, 95);
409 pPlateauVsEta = new TProfile("pPlateauVsEta","; #eta; #LT dE/dx #GT_{Plateau, primary tracks}",16,-0.8,0.8, DeDxMIPMax, 95);
411 fListOfObjects->Add(hMIPVsEta);
412 fListOfObjects->Add(pMIPVsEta);
413 fListOfObjects->Add(hMIPVsEtaV0s);
414 fListOfObjects->Add(pMIPVsEtaV0s);
415 fListOfObjects->Add(hPlateauVsEta);
416 fListOfObjects->Add(pPlateauVsEta);
424 for(Int_t i = 0; i < nHists; i++) {
425 for(Int_t pid = 0; pid < 7; pid++) {
427 hMcIn[pid][i] = new TH1D(Form("hIn%s%s", Pid[pid],ending[i]), Form("MC in (pid %s)", Pid[pid]),
428 nPtBinsV0s, ptBinsV0s);
429 fListOfObjects->Add(hMcIn[pid][i]);
431 hMcOut[pid][i] = new TH1D(Form("hMcOut%s%s", Pid[pid],ending[i]), Form("MC out (pid %s)", Pid[pid]),
432 nPtBinsV0s, ptBinsV0s);
433 fListOfObjects->Add(hMcOut[pid][i]);
439 fVtxMC = new TH1F("fVtxMC","mc vtx", 120, -30, 30);
440 fListOfObjects->Add(fVtxMC);
445 PostData(1, fListOfObjects);
448 //______________________________________________________________________________
449 void AliAnalysisTaskQAHighPtDeDx::UserExec(Option_t *)
454 // First we make sure that we have valid input(s)!
459 AliVEvent *event = InputEvent();
461 Error("UserExec", "Could not retrieve event");
467 if (fAnalysisType == "ESD"){
468 fESD = dynamic_cast<AliESDEvent*>(event);
470 Printf("%s:%d ESDEvent not found in Input Manager",(char*)__FILE__,__LINE__);
475 fAOD = dynamic_cast<AliAODEvent*>(event);
477 Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
487 if (fAnalysisType == "ESD"){
488 fMC = dynamic_cast<AliMCEvent*>(MCEvent());
490 Printf("%s:%d MCEvent not found in Input Manager",(char*)__FILE__,__LINE__);
495 fMCStack = fMC->Stack();
498 Printf("%s:%d MCStack not found in Input Manager",(char*)__FILE__,__LINE__);
504 fMC = dynamic_cast<AliMCEvent*>(MCEvent());
508 fMCArray = (TClonesArray*)fAOD->FindListObject("mcparticles");
510 Printf("%s:%d AOD MC array not found in Input Manager",(char*)__FILE__,__LINE__);
518 // Get trigger decision
519 fTriggeredEventMB = 0; //init
523 if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
524 ->IsEventSelected() & ftrigBit ){
525 fTriggeredEventMB = 1; //event triggered as minimum bias
528 // Get process type for MC
529 fMcProcessType = 0; // -1=invalid, 0=data, 1=ND, 2=SD, 3=DD
531 // real data that are not triggered we skip
532 if(!fAnalysisMC && !fTriggeredEventMB)
542 if (fAnalysisType == "ESD"){
546 AliHeader* headerMC = fMC->Header();
549 AliGenEventHeader* genHeader = headerMC->GenEventHeader();
550 TArrayF vtxMC(3); // primary vertex MC
551 vtxMC[0]=9999; vtxMC[1]=9999; vtxMC[2]=9999; //initialize with dummy
553 genHeader->PrimaryVertex(vtxMC);
558 AliGenPythiaEventHeader* pythiaGenHeader =
559 dynamic_cast<AliGenPythiaEventHeader*>(headerMC->GenEventHeader());
560 if (pythiaGenHeader) { //works only for pythia
561 fMcProcessType = GetPythiaEventProcessType(pythiaGenHeader->ProcessType());
564 AliGenDPMjetEventHeader* dpmJetGenHeader =
565 dynamic_cast<AliGenDPMjetEventHeader*>(headerMC->GenEventHeader());
566 if (dpmJetGenHeader) {
567 fMcProcessType = GetDPMjetEventProcessType(dpmJetGenHeader->ProcessType());
574 AliAODMCHeader* mcHeader = dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject("mcHeader"));
578 fZvtxMC = mcHeader->GetVtxZ();
582 if(strstr(mcHeader->GetGeneratorName(), "Pythia")) {
583 fMcProcessType = GetPythiaEventProcessType(mcHeader->GetEventType());
585 fMcProcessType = GetDPMjetEventProcessType(mcHeader->GetEventType());
595 if (fAnalysisType == "ESD"){
597 const AliESDVertex *vtxESD = fESD->GetPrimaryVertexTracks();
598 if(vtxESD->GetNContributors()<1) {
600 vtxESD = fESD->GetPrimaryVertexSPD();
601 /* quality checks on SPD-vertex */
602 TString vertexType = vtxESD->GetTitle();
603 if (vertexType.Contains("vertexer: Z") && (vtxESD->GetDispersion() > 0.04 || vtxESD->GetZRes() > 0.25))
604 fZvtx = -1599; //vertex = 0x0; //
605 else if (vtxESD->GetNContributors()<1)
606 fZvtx = -999; //vertex = 0x0; //
608 fZvtx = vtxESD->GetZ();
611 fZvtx = vtxESD->GetZ();
615 fZvtx = GetVertex(fAOD);
617 fVtxBeforeCuts->Fill(fZvtx);
619 //cut on the z position of vertex
620 if (TMath::Abs(fZvtx) > fVtxCut) {
630 Float_t centrality = -10;
632 // only analyze triggered events
633 if(fTriggeredEventMB) {
635 if (fAnalysisType == "ESD"){
637 AliCentrality *centObject = fESD->GetCentrality();
638 centrality = centObject->GetCentralityPercentile(fCentEst);
639 centralityGlobal = centrality;
640 if((centrality>fMaxCent)||(centrality<fMinCent))return;
642 fcent->Fill(centrality);
645 if(TMath::Abs(fZvtxMC)<fVtxCut){
647 fVtxMC->Fill(fZvtxMC);
653 AliCentrality *centObject = fAOD->GetCentrality();
655 centrality = centObject->GetCentralityPercentile(fCentEst);
657 //cout<<"centrality="<<centrality<<endl;
658 if((centrality>fMaxCent)||(centrality<fMinCent))return;
660 fcent->Fill(centrality);
663 if(TMath::Abs(fZvtxMC)<fVtxCut){
666 fVtxMC->Fill(fZvtxMC);
673 fVtxAfterCuts->Fill(fZvtx);
679 PostData(1, fListOfObjects);
682 //________________________________________________________________________
683 void AliAnalysisTaskQAHighPtDeDx::AnalyzeESD(AliESDEvent* esdEvent)
685 fRun = esdEvent->GetRunNumber();
687 if(esdEvent->GetHeader())
688 fEventId = GetEventIdAsLong(esdEvent->GetHeader());
690 cout << "centrality=" << centralityGlobal << endl;
693 Bool_t isPileup = esdEvent->IsPileupFromSPD();
700 // Int_t event = esdEvent->GetEventNumberInFile();
701 //UInt_t time = esdEvent->GetTimeStamp();
702 // ULong64_t trigger = esdEvent->GetTriggerMask();
703 magf = esdEvent->GetMagneticField();
709 if(fTriggeredEventMB) {// Only MC case can we have not triggered events
715 //Change, 10/04/13. Now accept all events to do a correct normalization
716 //if(fVtxStatus!=1) return; // accepted vertex
717 //Int_t nESDTracks = esdEvent->GetNumberOfTracks();
719 ProduceArrayTrksESD( esdEvent );//produce array with global track parameters
720 ProduceArrayV0ESD( esdEvent );//v0's
728 } // end if triggered
733 //________________________________________________________________________
734 void AliAnalysisTaskQAHighPtDeDx::AnalyzeAOD(AliAODEvent* aodEvent)
736 fRun = aodEvent->GetRunNumber();
738 if(aodEvent->GetHeader())
739 fEventId = GetEventIdAsLong(aodEvent->GetHeader());
741 //UInt_t time = 0; // Missing AOD info? aodEvent->GetTimeStamp();
742 magf = aodEvent->GetMagneticField();
744 //Int_t trackmult = 0; // no pt cuts
747 Bool_t isPileup = aodEvent->IsPileupFromSPD();
755 if(fTriggeredEventMB) {// Only MC case can we have not triggered events
760 //if(fVtxStatus!=1) return; // accepted vertex
761 //Int_t nAODTracks = aodEvent->GetNumberOfTracks();
763 ProduceArrayTrksAOD( aodEvent );
764 ProduceArrayV0AOD( aodEvent );
771 } // end if triggered
775 //_____________________________________________________________________________
776 Float_t AliAnalysisTaskQAHighPtDeDx::GetVertex(const AliVEvent* event) const
780 const AliVVertex* primaryVertex = event->GetPrimaryVertex();
782 if(primaryVertex->GetNContributors()>0)
783 zvtx = primaryVertex->GetZ();
788 //_____________________________________________________________________________
789 Short_t AliAnalysisTaskQAHighPtDeDx::GetPidCode(Int_t pdgCode) const
791 // return our internal code for pions, kaons, and protons
795 switch (TMath::Abs(pdgCode)) {
803 pidCode = 3; // proton
806 pidCode = 4; // electron
812 pidCode = 6; // something else?
818 //_____________________________________________________________________________
819 void AliAnalysisTaskQAHighPtDeDx::ProcessMCTruthESD()
821 // Fill the special MC histogram with the MC truth info
823 const Int_t nTracksMC = fMCStack->GetNtrack();
825 for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
828 if(!(fMCStack->IsPhysicalPrimary(iTracks)))
831 TParticle* trackMC = fMCStack->Particle(iTracks);
833 TParticlePDG* pdgPart = trackMC ->GetPDG();
834 Double_t chargeMC = pdgPart->Charge();
839 if (TMath::Abs(trackMC->Eta()) > fEtaCut )
842 Double_t etaMC = trackMC->Eta();
843 Int_t pdgCode = trackMC->GetPdgCode();
844 Short_t pidCodeMC = 0;
845 pidCodeMC = GetPidCode(pdgCode);
848 for(Int_t nh = 0; nh < 9; nh++) {
850 if( etaMC > etaHigh[nh]/10.0 || etaMC < etaLow[nh]/10.0 )
853 hMcIn[0][nh]->Fill(trackMC->Pt());
854 hMcIn[pidCodeMC][nh]->Fill(trackMC->Pt());
865 //_____________________________________________________________________________
866 void AliAnalysisTaskQAHighPtDeDx::ProcessMCTruthAOD()
868 // Fill the special MC histogram with the MC truth info
871 const Int_t nTracksMC = fMCArray->GetEntriesFast();
873 for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
875 AliAODMCParticle* trackMC = dynamic_cast<AliAODMCParticle*>(fMCArray->At(iTracks));
878 AliError("Cannot get MC particle");
883 if(!(trackMC->IsPhysicalPrimary()))
887 Double_t chargeMC = trackMC->Charge();
892 if (TMath::Abs(trackMC->Eta()) > fEtaCut )
895 Double_t etaMC = trackMC->Eta();
896 Int_t pdgCode = trackMC->GetPdgCode();
897 Short_t pidCodeMC = 0;
898 pidCodeMC = GetPidCode(pdgCode);
900 //cout<<"pidcode="<<pidCodeMC<<endl;
901 for(Int_t nh = 0; nh < 9; nh++) {
903 if( etaMC > etaHigh[nh]/10.0 || etaMC < etaLow[nh]/10.0 )
906 hMcIn[0][nh]->Fill(trackMC->Pt());
907 hMcIn[pidCodeMC][nh]->Fill(trackMC->Pt());
918 //_____________________________________________________________________________
919 Short_t AliAnalysisTaskQAHighPtDeDx::GetPythiaEventProcessType(Int_t pythiaType) {
921 // Get the process type of the event. PYTHIA
925 Short_t globalType = -1; //init
927 if(pythiaType==92||pythiaType==93){
928 globalType = 2; //single diffractive
930 else if(pythiaType==94){
931 globalType = 3; //double diffractive
933 //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??
935 globalType = 1; //non diffractive
940 //_____________________________________________________________________________
941 Short_t AliAnalysisTaskQAHighPtDeDx::GetDPMjetEventProcessType(Int_t dpmJetType) {
943 // get the process type of the event. PHOJET
946 // can only read pythia headers, either directly or from cocktalil header
947 Short_t globalType = -1;
949 if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
952 else if (dpmJetType==5 || dpmJetType==6) {
955 else if (dpmJetType==7) {
961 //_____________________________________________________________________________
962 ULong64_t AliAnalysisTaskQAHighPtDeDx::GetEventIdAsLong(AliVHeader* header) const
964 // To have a unique id for each event in a run!
965 // Modified from AliRawReader.h
966 return ((ULong64_t)header->GetBunchCrossNumber()+
967 (ULong64_t)header->GetOrbitNumber()*3564+
968 (ULong64_t)header->GetPeriodNumber()*16777215*3564);
972 //____________________________________________________________________
973 TParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMother(AliStack* stack, Int_t label)
976 // Finds the first mother among the primary particles of the particle identified by <label>,
977 // i.e. the primary that "caused" this particle
979 // Taken from AliPWG0Helper class
982 Int_t motherLabel = FindPrimaryMotherLabel(stack, label);
986 return stack->Particle(motherLabel);
989 //____________________________________________________________________
990 Int_t AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherLabel(AliStack* stack, Int_t label)
993 // Finds the first mother among the primary particles of the particle identified by <label>,
994 // i.e. the primary that "caused" this particle
998 // Taken from AliPWG0Helper class
1000 const Int_t nPrim = stack->GetNprimary();
1002 while (label >= nPrim) {
1004 //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
1006 TParticle* particle = stack->Particle(label);
1009 AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
1014 if (particle->GetMother(0) < 0) {
1016 AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
1020 label = particle->GetMother(0);
1026 //____________________________________________________________________
1027 AliAODMCParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherAOD(AliAODMCParticle* startParticle)
1030 // Finds the first mother among the primary particles of the particle identified by <label>,
1031 // i.e. the primary that "caused" this particle
1033 // Taken from AliPWG0Helper class
1036 AliAODMCParticle* mcPart = startParticle;
1041 if(mcPart->IsPrimary())
1044 Int_t mother = mcPart->GetMother();
1046 mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
1053 //V0______________________________________
1054 //____________________________________________________________________
1055 TParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherV0(AliStack* stack, Int_t label)
1058 // Finds the first mother among the primary particles of the particle identified by <label>,
1059 // i.e. the primary that "caused" this particle
1061 // Taken from AliPWG0Helper class
1066 Int_t motherLabel = FindPrimaryMotherLabelV0(stack, label, nSteps);
1067 if (motherLabel < 0)
1070 return stack->Particle(motherLabel);
1073 //____________________________________________________________________
1074 Int_t AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherLabelV0(AliStack* stack, Int_t label, Int_t& nSteps)
1077 // Finds the first mother among the primary particles of the particle identified by <label>,
1078 // i.e. the primary that "caused" this particle
1080 // returns its label
1082 // Taken from AliPWG0Helper class
1085 const Int_t nPrim = stack->GetNprimary();
1087 while (label >= nPrim) {
1089 //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
1091 nSteps++; // 1 level down
1093 TParticle* particle = stack->Particle(label);
1096 AliDebugGeneral("FindPrimaryMotherLabelV0", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
1101 if (particle->GetMother(0) < 0) {
1103 AliDebugGeneral("FindPrimaryMotherLabelV0", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
1107 label = particle->GetMother(0);
1113 //____________________________________________________________________
1114 AliAODMCParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherAODV0(AliAODMCParticle* startParticle, Int_t& nSteps)
1117 // Finds the first mother among the primary particles of the particle identified by <label>,
1118 // i.e. the primary that "caused" this particle
1120 // Taken from AliPWG0Helper class
1125 AliAODMCParticle* mcPart = startParticle;
1130 if(mcPart->IsPrimary())
1133 Int_t mother = mcPart->GetMother();
1135 mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
1136 nSteps++; // 1 level down
1144 //__________________________________________________________________
1145 void AliAnalysisTaskQAHighPtDeDx::ProduceArrayTrksESD( AliESDEvent *ESDevent ){
1147 const Int_t nESDTracks = ESDevent->GetNumberOfTracks();
1148 //Int_t trackmult=0;
1153 //get multiplicity tpc only track cuts
1154 for(Int_t iT = 0; iT < nESDTracks; iT++) {
1156 AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
1159 if(TMath::Abs(esdTrack->Eta()) > fEtaCut)
1162 //only golden track cuts
1163 UInt_t selectDebug = 0;
1164 if (fTrackFilterTPC) {
1165 selectDebug = fTrackFilterTPC->IsSelected(esdTrack);
1167 //cout<<"this is not a golden track"<<endl;
1178 for(Int_t iT = 0; iT < nESDTracks; iT++) {
1180 AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
1183 if(TMath::Abs(esdTrack->Eta()) > fEtaCut)
1186 //only golden track cuts
1187 UInt_t selectDebug = 0;
1188 if (fTrackFilterGolden) {
1189 selectDebug = fTrackFilterGolden->IsSelected(esdTrack);
1191 //cout<<"this is not a golden track"<<endl;
1197 //Int_t sharedtpcclusters = esdTrack->GetTPCnclsS();
1198 Short_t ncl = esdTrack->GetTPCsignalN();
1203 Double_t eta = esdTrack->Eta();
1204 Double_t phi = esdTrack->Phi();
1205 Double_t momentum = esdTrack->P();
1206 Float_t dedx = esdTrack->GetTPCsignal();
1208 //cout<<"magf="<<magf<<endl;
1210 if(!PhiCut(esdTrack->Pt(), phi, esdTrack->Charge(), magf, cutLow, cutHigh))
1218 Bool_t IsTOFout=kFALSE;
1219 if ((esdTrack->GetStatus()&AliESDtrack::kTOFout)==0)
1221 Float_t lengthtrack=esdTrack->GetIntegratedLength();
1222 Float_t timeTOF=esdTrack->GetTOFsignal();
1223 Double_t inttime[5]={0,0,0,0,0};
1224 esdTrack->GetIntegratedTimes(inttime);// Returns the array with integrated times for each particle hypothesis
1227 if ( ( lengthtrack != 0 ) && ( timeTOF != 0) )
1228 beta = inttime[0] / timeTOF;
1231 Short_t pidCode = 0;
1235 const Int_t label = TMath::Abs(esdTrack->GetLabel());
1236 TParticle* mcTrack = fMCStack->Particle(label);
1239 Int_t pdgCode = mcTrack->GetPdgCode();
1240 pidCode = GetPidCode(pdgCode);
1247 if( momentum <= 0.6 && momentum >= 0.4 ){//only p:0.4-0.6 GeV, pion MIP
1249 if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
1250 if(momentum<0.6&&momentum>0.4){
1251 hMIPVsEta->Fill(eta,dedx);
1252 pMIPVsEta->Fill(eta,dedx);
1255 if( dedx > DeDxMIPMax+1 && dedx < 95 ){
1256 if(TMath::Abs(beta-1)<0.1){
1257 hPlateauVsEta->Fill(eta,dedx);
1258 pPlateauVsEta->Fill(eta,dedx);
1264 for(Int_t nh = 0; nh < 9; nh++) {
1266 if( eta > etaHigh[nh]/10.0 || eta < etaLow[nh]/10.0 )
1271 hMcOut[0][nh]->Fill(esdTrack->Pt());
1272 hMcOut[pidCode][nh]->Fill(esdTrack->Pt());
1275 histAllCh[nh]->Fill(momentum, dedx);
1277 histPiTof[nh]->Fill(momentum, dedx);
1278 histpPiTof[nh]->Fill(momentum);
1281 if( momentum <= 0.6 && momentum >= 0.4 ){//only p:0.4-0.6 GeV, pion MIP
1282 //Fill pion MIP, before calibration
1283 if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
1284 hMIPVsPhi[nh]->Fill(phi,dedx);
1285 pMIPVsPhi[nh]->Fill(phi,dedx);
1288 hMIPVsNch[nh]->Fill(multTPC,dedx);
1289 pMIPVsNch[nh]->Fill(multTPC,dedx);
1293 //Fill electrons, before calibration
1294 if( dedx > DeDxMIPMax+1 && dedx < 95 ){
1295 if(TMath::Abs(beta-1)<0.1){
1296 hPlateauVsPhi[nh]->Fill(phi,dedx);
1297 pPlateauVsPhi[nh]->Fill(phi,dedx);
1304 }//end loop over eta intervals
1310 }//end of track loop
1316 //__________________________________________________________________
1317 void AliAnalysisTaskQAHighPtDeDx::ProduceArrayTrksAOD( AliAODEvent *AODevent ){
1320 Int_t nAODTracks = AODevent->GetNumberOfTracks();
1323 //get multiplicity tpc only track cuts
1324 for(Int_t iT = 0; iT < nAODTracks; iT++) {
1326 AliAODTrack* aodTrack = AODevent->GetTrack(iT);
1328 if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
1332 if (fTrackFilterTPC) {
1333 // TPC only cuts is bit 1
1334 if(!aodTrack->TestFilterBit(1))
1343 for(Int_t iT = 0; iT < nAODTracks; iT++) {
1345 AliAODTrack* aodTrack = AODevent->GetTrack(iT);
1347 if (fTrackFilterGolden) {
1348 // "Global track RAA analysis QM2011 + Chi2ITS<36"; bit 1024
1349 if(!aodTrack->TestFilterBit(1024))
1354 if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
1358 Double_t eta = aodTrack->Eta();
1359 Double_t phi = aodTrack->Phi();
1360 Double_t momentum = aodTrack->P();
1363 if(!PhiCut(aodTrack->Pt(), phi, aodTrack->Charge(), magf, cutLow, cutHigh))
1368 AliAODPid* aodPid = aodTrack->GetDetPid();
1377 ncl = aodPid->GetTPCsignalN();
1378 dedx = aodPid->GetTPCsignal();
1380 Bool_t IsTOFout=kFALSE;
1381 Float_t lengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF
1382 Float_t timeTOF=-999;
1384 if ((aodTrack->GetStatus()&AliESDtrack::kTOFout)==0)
1387 lengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF
1389 timeTOF=aodPid->GetTOFsignal();
1391 Double_t inttime[5]={0,0,0,0,0};
1392 aodPid->GetIntegratedTimes(inttime);// Returns the array with integrated times for each particle hypothesis
1396 if ( ( lengthtrack != 0 ) && ( timeTOF != 0) )
1397 beta = inttime[0] / timeTOF;
1407 Short_t pidCode = 0;
1411 const Int_t label = TMath::Abs(aodTrack->GetLabel());
1412 AliAODMCParticle* mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(label));
1416 Int_t pdgCode = mcTrack->GetPdgCode();
1417 pidCode = GetPidCode(pdgCode);
1426 //if(momentum>0.6||momentum<0.4)//only p:0.4-0.6 GeV
1431 if( momentum <= 0.6 && momentum >= 0.4 ){//only p:0.4-0.6 GeV, pion MIP
1432 if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
1433 if(momentum<0.6&&momentum>0.4){
1434 hMIPVsEta->Fill(eta,dedx);
1435 pMIPVsEta->Fill(eta,dedx);
1438 if( dedx > DeDxMIPMax+1 && dedx < 95 ){
1439 if(TMath::Abs(beta-1)<0.1){
1440 hPlateauVsEta->Fill(eta,dedx);
1441 pPlateauVsEta->Fill(eta,dedx);
1447 for(Int_t nh = 0; nh < 9; nh++) {
1449 if( eta > etaHigh[nh]/10.0 || eta < etaLow[nh]/10.0 )
1454 hMcOut[0][nh]->Fill(aodTrack->Pt());
1455 hMcOut[pidCode][nh]->Fill(aodTrack->Pt());
1458 histAllCh[nh]->Fill(momentum, dedx);
1460 histPiTof[nh]->Fill(momentum, dedx);
1461 histpPiTof[nh]->Fill(momentum);
1464 if( momentum <= 0.6 && momentum >= 0.4 ){//only p:0.4-0.6 GeV, pion MIP
1465 //Fill pion MIP, before calibration
1466 if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
1467 hMIPVsPhi[nh]->Fill(phi,dedx);
1468 pMIPVsPhi[nh]->Fill(phi,dedx);
1471 hMIPVsNch[nh]->Fill(multTPC,dedx);
1472 pMIPVsNch[nh]->Fill(multTPC,dedx);
1476 //Fill electrons, before calibration
1477 if( dedx > DeDxMIPMax+1 && dedx < 95 ){
1478 if(TMath::Abs(beta-1)<0.1){
1479 hPlateauVsPhi[nh]->Fill(phi,dedx);
1480 pPlateauVsPhi[nh]->Fill(phi,dedx);
1486 }//end loop over eta intervals
1492 }//end of track loop
1501 //__________________________________________________
1502 void AliAnalysisTaskQAHighPtDeDx::ProduceArrayV0ESD( AliESDEvent *ESDevent ){
1503 Int_t nv0s = ESDevent->GetNumberOfV0s();
1509 const AliESDVertex *myBestPrimaryVertex = ESDevent->GetPrimaryVertex();
1510 if (!myBestPrimaryVertex) return;
1511 if (!(myBestPrimaryVertex->GetStatus())) return;
1512 Double_t lPrimaryVtxPosition[3];
1513 myBestPrimaryVertex->GetXYZ(lPrimaryVtxPosition);
1515 Double_t lPrimaryVtxCov[6];
1516 myBestPrimaryVertex->GetCovMatrix(lPrimaryVtxCov);
1517 Double_t lPrimaryVtxChi2 = myBestPrimaryVertex->GetChi2toNDF();
1519 AliAODVertex* myPrimaryVertex = new AliAODVertex(lPrimaryVtxPosition, lPrimaryVtxCov, lPrimaryVtxChi2, NULL, -1, AliAODVertex::kPrimary);
1524 // LOOP OVER V0s, K0s, L, AL
1528 for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
1530 // This is the begining of the V0 loop
1531 AliESDv0 *esdV0 = ESDevent->GetV0(iV0);
1532 if (!esdV0) continue;
1534 //check onfly status
1535 if(esdV0->GetOnFlyStatus()!=0)
1539 // AliESDTrack (V0 Daughters)
1540 UInt_t lKeyPos = (UInt_t)TMath::Abs(esdV0->GetPindex());
1541 UInt_t lKeyNeg = (UInt_t)TMath::Abs(esdV0->GetNindex());
1543 AliESDtrack *pTrack = ESDevent->GetTrack(lKeyPos);
1544 AliESDtrack *nTrack = ESDevent->GetTrack(lKeyNeg);
1545 if (!pTrack || !nTrack) {
1546 Printf("ERROR: Could not retreive one of the daughter track");
1551 if (pTrack->GetSign() == nTrack->GetSign()) {
1552 //cout<< "like sign, continue"<< endl;
1556 // Eta cut on decay products
1557 if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
1561 // Check if switch does anything!
1562 Bool_t isSwitched = kFALSE;
1563 if (pTrack->GetSign() < 0) { // switch
1566 AliESDtrack* helpTrack = nTrack;
1571 //Double_t lV0Position[3];
1572 //esdV0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[2]);
1575 AliKFVertex primaryVtxKF( *myPrimaryVertex );
1576 AliKFParticle::SetField(ESDevent->GetMagneticField());
1578 // Also implement switch here!!!!!!
1579 AliKFParticle* negEKF = 0; // e-
1580 AliKFParticle* posEKF = 0; // e+
1581 AliKFParticle* negPiKF = 0; // pi -
1582 AliKFParticle* posPiKF = 0; // pi +
1583 AliKFParticle* posPKF = 0; // p
1584 AliKFParticle* negAPKF = 0; // p-bar
1587 negEKF = new AliKFParticle( *(esdV0->GetParamN()) , 11);
1588 posEKF = new AliKFParticle( *(esdV0->GetParamP()) ,-11);
1589 negPiKF = new AliKFParticle( *(esdV0->GetParamN()) ,-211);
1590 posPiKF = new AliKFParticle( *(esdV0->GetParamP()) , 211);
1591 posPKF = new AliKFParticle( *(esdV0->GetParamP()) , 2212);
1592 negAPKF = new AliKFParticle( *(esdV0->GetParamN()) ,-2212);
1593 } else { // switch + and -
1594 negEKF = new AliKFParticle( *(esdV0->GetParamP()) , 11);
1595 posEKF = new AliKFParticle( *(esdV0->GetParamN()) ,-11);
1596 negPiKF = new AliKFParticle( *(esdV0->GetParamP()) ,-211);
1597 posPiKF = new AliKFParticle( *(esdV0->GetParamN()) , 211);
1598 posPKF = new AliKFParticle( *(esdV0->GetParamN()) , 2212);
1599 negAPKF = new AliKFParticle( *(esdV0->GetParamP()) ,-2212);
1602 AliKFParticle v0GKF; // Gamma e.g. from pi0
1605 v0GKF.SetProductionVertex(primaryVtxKF);
1607 AliKFParticle v0K0sKF; // K0 short
1608 v0K0sKF+=(*negPiKF);
1609 v0K0sKF+=(*posPiKF);
1610 v0K0sKF.SetProductionVertex(primaryVtxKF);
1612 AliKFParticle v0LambdaKF; // Lambda
1613 v0LambdaKF+=(*negPiKF);
1614 v0LambdaKF+=(*posPKF);
1615 v0LambdaKF.SetProductionVertex(primaryVtxKF);
1617 AliKFParticle v0AntiLambdaKF; // Lambda-bar
1618 v0AntiLambdaKF+=(*posPiKF);
1619 v0AntiLambdaKF+=(*negAPKF);
1620 v0AntiLambdaKF.SetProductionVertex(primaryVtxKF);
1622 Double_t dmassG = v0GKF.GetMass();
1623 Double_t dmassK = v0K0sKF.GetMass()-0.498;
1624 Double_t dmassL = v0LambdaKF.GetMass()-1.116;
1625 Double_t dmassAL = v0AntiLambdaKF.GetMass()-1.116;
1628 for( Int_t case_v0 = 0; case_v0 < 2; ++case_v0 ){
1634 Bool_t fillPos = kFALSE;
1635 Bool_t fillNeg = kFALSE;
1640 if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01){
1654 if(fillPos||fillNeg)
1661 for(Int_t j = 0; j < 2; j++) {
1663 AliESDtrack* track = 0;
1679 if(track->GetTPCsignalN()<=70)continue;
1680 Double_t phi = track->Phi();
1682 if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
1686 //if(!PhiCut(pt, phi, charge, magf, cutLow, cutHigh, hPhi))
1689 Double_t eta = track->Eta();
1690 Double_t momentum = track->Pt();
1691 Double_t dedx = track->GetTPCsignal();
1693 if(fillPos&&fillNeg){
1696 if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
1697 if(momentum<0.6&&momentum>0.4){
1698 hMIPVsEtaV0s->Fill(eta,dedx);
1699 pMIPVsEtaV0s->Fill(eta,dedx);
1706 for(Int_t nh = 0; nh < nHists; nh++) {
1710 if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
1713 if(fillPos&&fillNeg){
1715 histPiV0[nh]->Fill(momentum, dedx);
1716 histpPiV0[nh]->Fill(momentum);
1721 histPV0[nh]->Fill(momentum, dedx);
1722 histpPV0[nh]->Fill(momentum);
1728 }//end loop over two tracks
1735 Bool_t fillPos = kFALSE;
1736 Bool_t fillNeg = kFALSE;
1740 if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01) {
1741 if(dmassG<0.01 && dmassG>0.0001) {
1744 if( TMath::Abs(nTrack->GetTPCsignal() - 85.0) < 5)
1746 if( TMath::Abs(pTrack->GetTPCsignal() - 85.0) < 5)
1754 //cout<<"fillPos="<<fillPos<<endl;
1756 if(fillPos == kTRUE && fillNeg == kTRUE)
1760 AliESDtrack* track = 0;
1768 Double_t dedx = track->GetTPCsignal();
1769 Double_t eta = track->Eta();
1770 Double_t phi = track->Phi();
1771 Double_t momentum = track->P();
1773 if(track->GetTPCsignalN()<=70)continue;
1775 if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
1778 for(Int_t nh = 0; nh < nHists; nh++) {
1780 if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
1783 histEV0[nh]->Fill(momentum, dedx);
1793 }//end loop over case V0
1796 // clean up loop over v0
1808 delete myPrimaryVertex;
1812 //__________________________________________________________________________
1813 void AliAnalysisTaskQAHighPtDeDx::ProduceArrayV0AOD( AliAODEvent *AODevent ){
1814 Int_t nv0s = AODevent->GetNumberOfV0s();
1820 AliAODVertex *myBestPrimaryVertex = AODevent->GetPrimaryVertex();
1821 if (!myBestPrimaryVertex) return;
1825 // ################################
1826 // #### BEGINNING OF V0 CODE ######
1827 // ################################
1828 // This is the begining of the V0 loop
1829 for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
1830 AliAODv0 *aodV0 = AODevent->GetV0(iV0);
1831 if (!aodV0) continue;
1834 //check onfly status
1835 if(aodV0->GetOnFlyStatus()!=0)
1838 // AliAODTrack (V0 Daughters)
1839 AliAODVertex* vertex = aodV0->GetSecondaryVtx();
1841 Printf("ERROR: Could not retrieve vertex");
1845 AliAODTrack *pTrack = (AliAODTrack*)vertex->GetDaughter(0);
1846 AliAODTrack *nTrack = (AliAODTrack*)vertex->GetDaughter(1);
1847 if (!pTrack || !nTrack) {
1848 Printf("ERROR: Could not retrieve one of the daughter track");
1853 if (pTrack->Charge() == nTrack->Charge()) {
1854 //cout<< "like sign, continue"<< endl;
1858 // Make sure charge ordering is ok
1859 if (pTrack->Charge() < 0) {
1860 AliAODTrack* helpTrack = pTrack;
1865 // Eta cut on decay products
1866 if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
1870 Double_t dmassG = aodV0->InvMass2Prongs(0,1,11,11);
1871 Double_t dmassK = aodV0->MassK0Short()-0.498;
1872 Double_t dmassL = aodV0->MassLambda()-1.116;
1873 Double_t dmassAL = aodV0->MassAntiLambda()-1.116;
1875 for( Int_t case_v0 = 0; case_v0 < 2; ++case_v0 ){
1880 Bool_t fillPos = kFALSE;
1881 Bool_t fillNeg = kFALSE;
1887 if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01){
1901 if(fillPos||fillNeg)
1909 for(Int_t j = 0; j < 2; j++) {
1911 AliAODTrack* track = 0;
1927 if(track->GetTPCsignalN()<=70)continue;
1929 Double_t phi = track->Phi();
1931 if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
1935 //if(!PhiCut(pt, phi, charge, magf, cutLow, cutHigh, hPhi))
1938 Double_t eta = track->Eta();
1939 Double_t momentum = track->Pt();
1940 Double_t dedx = track->GetTPCsignal();
1942 if(fillPos&&fillNeg){
1945 if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
1946 if(momentum<0.6&&momentum>0.4){
1947 hMIPVsEtaV0s->Fill(eta,dedx);
1948 pMIPVsEtaV0s->Fill(eta,dedx);
1955 for(Int_t nh = 0; nh < nHists; nh++) {
1959 if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
1962 if(fillPos&&fillNeg){
1964 histPiV0[nh]->Fill(momentum, dedx);
1965 histpPiV0[nh]->Fill(momentum);
1970 histPV0[nh]->Fill(momentum, dedx);
1971 histpPV0[nh]->Fill(momentum);
1978 }//end loop over two tracks
1984 Bool_t fillPos = kFALSE;
1985 Bool_t fillNeg = kFALSE;
1987 if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01) {
1988 if(dmassG<0.01 && dmassG>0.0001) {
1990 if( TMath::Abs(nTrack->GetTPCsignal() - 85.0) < 5)
1992 if( TMath::Abs(pTrack->GetTPCsignal() - 85.0) < 5)
2001 if(fillPos == kTRUE && fillNeg == kTRUE)
2005 AliAODTrack* track = 0;
2013 Double_t dedx = track->GetTPCsignal();
2014 Double_t eta = track->Eta();
2015 Double_t phi = track->Phi();
2016 Double_t momentum = track->P();
2018 if(track->GetTPCsignalN()<=70)continue;
2020 if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
2023 for(Int_t nh = 0; nh < nHists; nh++) {
2025 if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
2028 histEV0[nh]->Fill(momentum, dedx);
2037 }//end loop over V0s cases
2039 }//end loop over v0's
2045 Bool_t AliAnalysisTaskQAHighPtDeDx::PhiCut(Double_t pt, Double_t phi, Double_t q, Float_t mag, TF1* phiCutLow, TF1* phiCutHigh)
2050 //Double_t phi = track->Phi();
2051 if(mag < 0) // for negatve polarity field
2052 phi = TMath::TwoPi() - phi;
2053 if(q < 0) // for negatve charge
2054 phi = TMath::TwoPi()-phi;
2056 phi += TMath::Pi()/18.0; // to center gap in the middle
2057 phi = fmod(phi, TMath::Pi()/9.0);
2059 if(phi<phiCutHigh->Eval(pt)
2060 && phi>phiCutLow->Eval(pt))
2061 return kFALSE; // reject track
2063 hPhi->Fill(pt, phi);