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
142 AliAnalysisTaskQAHighPtDeDx::AliAnalysisTaskQAHighPtDeDx(const char *name):
143 AliAnalysisTaskSE(name),
149 fTrackFilterGolden(0x0),
150 fTrackFilterTPC(0x0),
152 fAnalysisType("ESD"),
154 fAnalysisPbPb(kFALSE),
162 fStoreMcIn(kFALSE),//
163 fMcProcessType(-999),
164 fTriggeredEventMB(-999),
171 fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
184 // Default constructor (should not be used)
185 for(Int_t i=0;i<9;++i){
187 hMIPVsNch[i]=0;//TH2D, MIP vs Nch for different eta intervals
188 pMIPVsNch[i]=0;//TProfile, MIP vs Nch for different eta intervals
189 hMIPVsPhi[i]=0;//TH2D, MIP vs phi for different eta intervals
190 pMIPVsPhi[i]=0;//TProfile, MIP vs phi for different eta intervals
191 hPlateauVsPhi[i]=0;//TH2D, dE/dx vs Phi, electrons 0.4<p<0.6 GeV/c
192 pPlateauVsPhi[i]=0;//TProfile, dE/dx vs Phi, electrons 0.4<p<0.6 GeV/c
193 histPiV0[i]=0;//TH2D, dE/dx vs p, pi id by V0s
194 histpPiV0[i]=0;//TH1D, pi id by V0s
195 histPV0[i]=0;// TH2D, dE/dx vs p, p id by V0s
196 histpPV0[i]=0;// TH1D, p id by V0s
197 histAllCh[i]=0;//TH2D, dE/dx vs p for all charged particles
198 histPiTof[i]=0;//TH2D, dE/dx vs p for a "clean" sample of pions, beta>1
199 histpPiTof[i]=0;//TH1D, for a "clean" sample of pions, beta>1
202 for(Int_t pid=0;pid<7;++pid){
208 DefineOutput(1, TList::Class());//esto es nuevo
214 AliAnalysisTaskQAHighPtDeDx::~AliAnalysisTaskQAHighPtDeDx() {
228 //______________________________________________________________________________
229 void AliAnalysisTaskQAHighPtDeDx::UserCreateOutputObjects()
231 // This method is called once per worker node
232 // Here we define the output: histograms and debug tree if requested
233 // We also create the random generator here so it might get different seeds...
234 fRandom = new TRandom(0); // 0 means random seed
239 fListOfObjects = new TList();
240 fListOfObjects->SetOwner();
248 fEvents = new TH1I("fEvents","Number of analyzed events; Events; Counts", 3, 0, 3);
249 fListOfObjects->Add(fEvents);
251 fn1=new TH1F("fn1","fn1",11,-1,10);
252 fListOfObjects->Add(fn1);
254 fcent=new TH1F("fcent","fcent",104,-2,102);
255 fListOfObjects->Add(fcent);
257 fVtx = new TH1I("fVtx","Vtx info (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
258 fListOfObjects->Add(fVtx);
260 fVtxBeforeCuts = new TH1F("fVtxBeforeCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
261 fListOfObjects->Add(fVtxBeforeCuts);
263 fVtxAfterCuts = new TH1F("fVtxAfterCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
264 fListOfObjects->Add(fVtxAfterCuts);
267 const Int_t nPtBinsV0s = 25;
268 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 ,
269 1.2 , 1.4 , 1.6 , 1.8 , 2.0 , 2.5 , 3.0 , 3.5 , 4.0 , 5.0 , 7.0 ,
270 9.0 , 12.0, 15.0, 20.0 };
275 const Char_t* ending[nHists] = {"", "86", "64", "42", "20", "02", "24", "46", "68"};
277 const Char_t* LatexEta[nHists] = {
278 "|#eta|<0.8", "-0.8<#eta<-0.6", "-0.6<#eta<-0.4", "-0.4<#eta<-0.2", "-0.2<#eta<0",
279 "0<#eta<0.2", "0.2<#eta<0.4", "0.4<#eta<0.6", "0.6<#eta<0.8"
281 hPhi = new TH2D("histPhi", "pt; #phi'", nPtBinsV0s, ptBinsV0s, 90, -0.05, 0.4);
282 //dE/dx vs phi, pions at the MIP
283 fListOfObjects->Add(hPhi);
290 for(Int_t i = 0; i < nHists; i++) {
293 hMIPVsPhi[i] = new TH2D(Form("hMIPVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx MIP",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(),
294 DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
295 hMIPVsPhi[i]->Sumw2();
297 pMIPVsPhi[i] = new TProfile(Form("pMIPVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx MIP",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(),
298 DeDxMIPMin, DeDxMIPMax);
299 pMIPVsPhi[i]->SetMarkerStyle(20);
300 pMIPVsPhi[i]->SetMarkerColor(1);
301 pMIPVsPhi[i]->SetLineColor(1);
302 pMIPVsPhi[i]->Sumw2();
304 fListOfObjects->Add(hMIPVsPhi[i]);
305 fListOfObjects->Add(pMIPVsPhi[i]);
307 hPlateauVsPhi[i] = new TH2D(Form("hPlateauVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx Plateau",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(),
308 95-DeDxMIPMax, DeDxMIPMax, 95);
309 hPlateauVsPhi[i]->Sumw2();
311 pPlateauVsPhi[i] = new TProfile(Form("pPlateauVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx Plateau",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(),
313 pPlateauVsPhi[i]->SetMarkerStyle(20);
314 pPlateauVsPhi[i]->SetMarkerColor(1);
315 pPlateauVsPhi[i]->SetLineColor(1);
316 pPlateauVsPhi[i]->Sumw2();
318 fListOfObjects->Add(hPlateauVsPhi[i]);
319 fListOfObjects->Add(pPlateauVsPhi[i]);
322 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,
323 DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
324 hMIPVsNch[i]->Sumw2();
326 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);
327 pMIPVsNch[i]->SetMarkerStyle(20);
328 pMIPVsNch[i]->SetMarkerColor(1);
329 pMIPVsNch[i]->SetLineColor(1);
330 pMIPVsNch[i]->Sumw2();
332 fListOfObjects->Add(hMIPVsNch[i]);
333 fListOfObjects->Add(pMIPVsNch[i]);
335 //two dimmesional distributions dE/dx vs p for secondary pions
336 histPiV0[i] = new TH2D(Form("histPiV0%s", ending[i]), "Pions id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
337 histPiV0[i]->Sumw2();
338 fListOfObjects->Add(histPiV0[i]);
340 histpPiV0[i] = new TH1D(Form("histPiV01D%s", ending[i]), "Pions id by V0; #it{p} (GeV/#it{c}); counts", 200, 0, 20);
341 histpPiV0[i]->Sumw2();
342 fListOfObjects->Add(histpPiV0[i]);
344 //two dimmesional distributions dE/dx vs p for secondary protons
345 histPV0[i] = new TH2D(Form("histPV0%s", ending[i]), "Protons id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
347 fListOfObjects->Add(histPV0[i]);
349 histpPV0[i] = new TH1D(Form("histPV01D%s", ending[i]), "Protons id by V0; #it{p} (GeV/#it{c}); counts", 200, 0, 20);
350 histpPV0[i]->Sumw2();
351 fListOfObjects->Add(histpPV0[i]);
353 //two dimmesional distributions dE/dx vs p for primary pions
354 histPiTof[i] = new TH2D(Form("histPiTof%s", ending[i]), "all charged", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
355 histPiTof[i]->Sumw2();
357 histpPiTof[i] = new TH1D(Form("histPiTof1D%s", ending[i]), "Protons id by V0; #it{p} (GeV/#it{c}); counts", 200, 0, 20);
358 histpPiTof[i]->Sumw2();
359 fListOfObjects->Add(histpPiTof[i]);
362 histAllCh[i] = new TH2D(Form("histAllCh%s", ending[i]), "Pions id by TOF", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
363 histAllCh[i]->Sumw2();
365 fListOfObjects->Add(histPiTof[i]);
366 fListOfObjects->Add(histAllCh[i]);
368 histEV0[i] = new TH2D(Form("histEV0%s", ending[i]), "Electrons id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
370 fListOfObjects->Add(histEV0[i]);
377 hMIPVsEta = new TH2D("hMIPVsEta","; #eta; dE/dx_{MIP, primary tracks}",16,-0.8,0.8,DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
378 pMIPVsEta = new TProfile("pMIPVsEta","; #eta; #LT dE/dx #GT_{MIP, primary tracks}",16,-0.8,0.8, DeDxMIPMin, DeDxMIPMax);
379 hMIPVsEtaV0s = new TH2D("hMIPVsEtaV0s","; #eta; dE/dx_{MIP, secondary tracks}",16,-0.8,0.8,DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
380 pMIPVsEtaV0s = new TProfile("pMIPVsEtaV0s","; #eta; #LT dE/dx #GT_{MIP, secondary tracks}",16,-0.8,0.8,DeDxMIPMin, DeDxMIPMax);
382 hPlateauVsEta = new TH2D("hPlateauVsEta","; #eta; dE/dx_{Plateau, primary tracks}",16,-0.8,0.8,95-DeDxMIPMax, DeDxMIPMax, 95);
383 pPlateauVsEta = new TProfile("pPlateauVsEta","; #eta; #LT dE/dx #GT_{Plateau, primary tracks}",16,-0.8,0.8, DeDxMIPMax, 95);
385 fListOfObjects->Add(hMIPVsEta);
386 fListOfObjects->Add(pMIPVsEta);
387 fListOfObjects->Add(hMIPVsEtaV0s);
388 fListOfObjects->Add(pMIPVsEtaV0s);
389 fListOfObjects->Add(hPlateauVsEta);
390 fListOfObjects->Add(pPlateauVsEta);
398 for(Int_t i = 0; i < nHists; i++) {
399 for(Int_t pid = 0; pid < 7; pid++) {
401 hMcIn[pid][i] = new TH1D(Form("hIn%s%s", Pid[pid],ending[i]), Form("MC in (pid %s)", Pid[pid]),
402 nPtBinsV0s, ptBinsV0s);
403 fListOfObjects->Add(hMcIn[pid][i]);
405 hMcOut[pid][i] = new TH1D(Form("hMcOut%s%s", Pid[pid],ending[i]), Form("MC out (pid %s)", Pid[pid]),
406 nPtBinsV0s, ptBinsV0s);
407 fListOfObjects->Add(hMcOut[pid][i]);
413 fVtxMC = new TH1F("fVtxMC","mc vtx", 120, -30, 30);
414 fListOfObjects->Add(fVtxMC);
419 PostData(1, fListOfObjects);
422 //______________________________________________________________________________
423 void AliAnalysisTaskQAHighPtDeDx::UserExec(Option_t *)
428 // First we make sure that we have valid input(s)!
433 AliVEvent *event = InputEvent();
435 Error("UserExec", "Could not retrieve event");
441 if (fAnalysisType == "ESD"){
442 fESD = dynamic_cast<AliESDEvent*>(event);
444 Printf("%s:%d ESDEvent not found in Input Manager",(char*)__FILE__,__LINE__);
449 fAOD = dynamic_cast<AliAODEvent*>(event);
451 Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
461 if (fAnalysisType == "ESD"){
462 fMC = dynamic_cast<AliMCEvent*>(MCEvent());
464 Printf("%s:%d MCEvent not found in Input Manager",(char*)__FILE__,__LINE__);
469 fMCStack = fMC->Stack();
472 Printf("%s:%d MCStack not found in Input Manager",(char*)__FILE__,__LINE__);
478 fMC = dynamic_cast<AliMCEvent*>(MCEvent());
482 fMCArray = (TClonesArray*)fAOD->FindListObject("mcparticles");
484 Printf("%s:%d AOD MC array not found in Input Manager",(char*)__FILE__,__LINE__);
492 // Get trigger decision
493 fTriggeredEventMB = 0; //init
497 if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
498 ->IsEventSelected() & ftrigBit ){
499 fTriggeredEventMB = 1; //event triggered as minimum bias
502 // Get process type for MC
503 fMcProcessType = 0; // -1=invalid, 0=data, 1=ND, 2=SD, 3=DD
505 // real data that are not triggered we skip
506 if(!fAnalysisMC && !fTriggeredEventMB)
516 if (fAnalysisType == "ESD"){
520 AliHeader* headerMC = fMC->Header();
523 AliGenEventHeader* genHeader = headerMC->GenEventHeader();
524 TArrayF vtxMC(3); // primary vertex MC
525 vtxMC[0]=9999; vtxMC[1]=9999; vtxMC[2]=9999; //initialize with dummy
527 genHeader->PrimaryVertex(vtxMC);
532 AliGenPythiaEventHeader* pythiaGenHeader =
533 dynamic_cast<AliGenPythiaEventHeader*>(headerMC->GenEventHeader());
534 if (pythiaGenHeader) { //works only for pythia
535 fMcProcessType = GetPythiaEventProcessType(pythiaGenHeader->ProcessType());
538 AliGenDPMjetEventHeader* dpmJetGenHeader =
539 dynamic_cast<AliGenDPMjetEventHeader*>(headerMC->GenEventHeader());
540 if (dpmJetGenHeader) {
541 fMcProcessType = GetDPMjetEventProcessType(dpmJetGenHeader->ProcessType());
548 AliAODMCHeader* mcHeader = dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject("mcHeader"));
552 fZvtxMC = mcHeader->GetVtxZ();
556 if(strstr(mcHeader->GetGeneratorName(), "Pythia")) {
557 fMcProcessType = GetPythiaEventProcessType(mcHeader->GetEventType());
559 fMcProcessType = GetDPMjetEventProcessType(mcHeader->GetEventType());
569 if (fAnalysisType == "ESD"){
571 const AliESDVertex *vtxESD = fESD->GetPrimaryVertexTracks();
572 if(vtxESD->GetNContributors()<1) {
574 vtxESD = fESD->GetPrimaryVertexSPD();
575 /* quality checks on SPD-vertex */
576 TString vertexType = vtxESD->GetTitle();
577 if (vertexType.Contains("vertexer: Z") && (vtxESD->GetDispersion() > 0.04 || vtxESD->GetZRes() > 0.25))
578 fZvtx = -1599; //vertex = 0x0; //
579 else if (vtxESD->GetNContributors()<1)
580 fZvtx = -999; //vertex = 0x0; //
582 fZvtx = vtxESD->GetZ();
585 fZvtx = vtxESD->GetZ();
589 fZvtx = GetVertex(fAOD);
591 fVtxBeforeCuts->Fill(fZvtx);
593 //cut on the z position of vertex
594 if (TMath::Abs(fZvtx) > fVtxCut) {
604 Float_t centrality = -10;
606 // only analyze triggered events
607 if(fTriggeredEventMB) {
609 if (fAnalysisType == "ESD"){
611 AliCentrality *centObject = fESD->GetCentrality();
612 centrality = centObject->GetCentralityPercentile(fCentEst);
613 centralityGlobal = centrality;
614 if((centrality>fMaxCent)||(centrality<fMinCent))return;
616 fcent->Fill(centrality);
619 if(TMath::Abs(fZvtxMC)<fVtxCut){
621 fVtxMC->Fill(fZvtxMC);
627 AliCentrality *centObject = fAOD->GetCentrality();
629 centrality = centObject->GetCentralityPercentile(fCentEst);
631 //cout<<"centrality="<<centrality<<endl;
632 if((centrality>fMaxCent)||(centrality<fMinCent))return;
634 fcent->Fill(centrality);
637 if(TMath::Abs(fZvtxMC)<fVtxCut){
640 fVtxMC->Fill(fZvtxMC);
647 fVtxAfterCuts->Fill(fZvtx);
653 PostData(1, fListOfObjects);
656 //________________________________________________________________________
657 void AliAnalysisTaskQAHighPtDeDx::AnalyzeESD(AliESDEvent* esdEvent)
659 fRun = esdEvent->GetRunNumber();
661 if(esdEvent->GetHeader())
662 fEventId = GetEventIdAsLong(esdEvent->GetHeader());
664 cout << "centrality=" << centralityGlobal << endl;
667 Bool_t isPileup = esdEvent->IsPileupFromSPD();
674 // Int_t event = esdEvent->GetEventNumberInFile();
675 //UInt_t time = esdEvent->GetTimeStamp();
676 // ULong64_t trigger = esdEvent->GetTriggerMask();
677 magf = esdEvent->GetMagneticField();
683 if(fTriggeredEventMB) {// Only MC case can we have not triggered events
689 //Change, 10/04/13. Now accept all events to do a correct normalization
690 //if(fVtxStatus!=1) return; // accepted vertex
691 //Int_t nESDTracks = esdEvent->GetNumberOfTracks();
693 ProduceArrayTrksESD( esdEvent );//produce array with global track parameters
694 ProduceArrayV0ESD( esdEvent );//v0's
702 } // end if triggered
707 //________________________________________________________________________
708 void AliAnalysisTaskQAHighPtDeDx::AnalyzeAOD(AliAODEvent* aodEvent)
710 fRun = aodEvent->GetRunNumber();
712 if(aodEvent->GetHeader())
713 fEventId = GetEventIdAsLong(aodEvent->GetHeader());
715 //UInt_t time = 0; // Missing AOD info? aodEvent->GetTimeStamp();
716 magf = aodEvent->GetMagneticField();
718 //Int_t trackmult = 0; // no pt cuts
721 Bool_t isPileup = aodEvent->IsPileupFromSPD();
729 if(fTriggeredEventMB) {// Only MC case can we have not triggered events
734 //if(fVtxStatus!=1) return; // accepted vertex
735 //Int_t nAODTracks = aodEvent->GetNumberOfTracks();
737 ProduceArrayTrksAOD( aodEvent );
738 ProduceArrayV0AOD( aodEvent );
745 } // end if triggered
749 //_____________________________________________________________________________
750 Float_t AliAnalysisTaskQAHighPtDeDx::GetVertex(const AliVEvent* event) const
754 const AliVVertex* primaryVertex = event->GetPrimaryVertex();
756 if(primaryVertex->GetNContributors()>0)
757 zvtx = primaryVertex->GetZ();
762 //_____________________________________________________________________________
763 Short_t AliAnalysisTaskQAHighPtDeDx::GetPidCode(Int_t pdgCode) const
765 // return our internal code for pions, kaons, and protons
769 switch (TMath::Abs(pdgCode)) {
777 pidCode = 3; // proton
780 pidCode = 4; // electron
786 pidCode = 6; // something else?
792 //_____________________________________________________________________________
793 void AliAnalysisTaskQAHighPtDeDx::ProcessMCTruthESD()
795 // Fill the special MC histogram with the MC truth info
797 const Int_t nTracksMC = fMCStack->GetNtrack();
799 for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
802 if(!(fMCStack->IsPhysicalPrimary(iTracks)))
805 TParticle* trackMC = fMCStack->Particle(iTracks);
807 TParticlePDG* pdgPart = trackMC ->GetPDG();
808 Double_t chargeMC = pdgPart->Charge();
813 if (TMath::Abs(trackMC->Eta()) > fEtaCut )
816 Double_t etaMC = trackMC->Eta();
817 Int_t pdgCode = trackMC->GetPdgCode();
818 Short_t pidCodeMC = 0;
819 pidCodeMC = GetPidCode(pdgCode);
822 for(Int_t nh = 0; nh < 9; nh++) {
824 if( etaMC > etaHigh[nh]/10.0 || etaMC < etaLow[nh]/10.0 )
827 hMcIn[0][nh]->Fill(trackMC->Pt());
828 hMcIn[pidCodeMC][nh]->Fill(trackMC->Pt());
839 //_____________________________________________________________________________
840 void AliAnalysisTaskQAHighPtDeDx::ProcessMCTruthAOD()
842 // Fill the special MC histogram with the MC truth info
845 const Int_t nTracksMC = fMCArray->GetEntriesFast();
847 for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
849 AliAODMCParticle* trackMC = dynamic_cast<AliAODMCParticle*>(fMCArray->At(iTracks));
852 AliError("Cannot get MC particle");
857 if(!(trackMC->IsPhysicalPrimary()))
861 Double_t chargeMC = trackMC->Charge();
866 if (TMath::Abs(trackMC->Eta()) > fEtaCut )
869 Double_t etaMC = trackMC->Eta();
870 Int_t pdgCode = trackMC->GetPdgCode();
871 Short_t pidCodeMC = 0;
872 pidCodeMC = GetPidCode(pdgCode);
874 //cout<<"pidcode="<<pidCodeMC<<endl;
875 for(Int_t nh = 0; nh < 9; nh++) {
877 if( etaMC > etaHigh[nh]/10.0 || etaMC < etaLow[nh]/10.0 )
880 hMcIn[0][nh]->Fill(trackMC->Pt());
881 hMcIn[pidCodeMC][nh]->Fill(trackMC->Pt());
892 //_____________________________________________________________________________
893 Short_t AliAnalysisTaskQAHighPtDeDx::GetPythiaEventProcessType(Int_t pythiaType) {
895 // Get the process type of the event. PYTHIA
899 Short_t globalType = -1; //init
901 if(pythiaType==92||pythiaType==93){
902 globalType = 2; //single diffractive
904 else if(pythiaType==94){
905 globalType = 3; //double diffractive
907 //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??
909 globalType = 1; //non diffractive
914 //_____________________________________________________________________________
915 Short_t AliAnalysisTaskQAHighPtDeDx::GetDPMjetEventProcessType(Int_t dpmJetType) {
917 // get the process type of the event. PHOJET
920 // can only read pythia headers, either directly or from cocktalil header
921 Short_t globalType = -1;
923 if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
926 else if (dpmJetType==5 || dpmJetType==6) {
929 else if (dpmJetType==7) {
935 //_____________________________________________________________________________
936 ULong64_t AliAnalysisTaskQAHighPtDeDx::GetEventIdAsLong(AliVHeader* header) const
938 // To have a unique id for each event in a run!
939 // Modified from AliRawReader.h
940 return ((ULong64_t)header->GetBunchCrossNumber()+
941 (ULong64_t)header->GetOrbitNumber()*3564+
942 (ULong64_t)header->GetPeriodNumber()*16777215*3564);
946 //____________________________________________________________________
947 TParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMother(AliStack* stack, Int_t label)
950 // Finds the first mother among the primary particles of the particle identified by <label>,
951 // i.e. the primary that "caused" this particle
953 // Taken from AliPWG0Helper class
956 Int_t motherLabel = FindPrimaryMotherLabel(stack, label);
960 return stack->Particle(motherLabel);
963 //____________________________________________________________________
964 Int_t AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherLabel(AliStack* stack, Int_t label)
967 // Finds the first mother among the primary particles of the particle identified by <label>,
968 // i.e. the primary that "caused" this particle
972 // Taken from AliPWG0Helper class
974 const Int_t nPrim = stack->GetNprimary();
976 while (label >= nPrim) {
978 //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
980 TParticle* particle = stack->Particle(label);
983 AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
988 if (particle->GetMother(0) < 0) {
990 AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
994 label = particle->GetMother(0);
1000 //____________________________________________________________________
1001 AliAODMCParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherAOD(AliAODMCParticle* startParticle)
1004 // Finds the first mother among the primary particles of the particle identified by <label>,
1005 // i.e. the primary that "caused" this particle
1007 // Taken from AliPWG0Helper class
1010 AliAODMCParticle* mcPart = startParticle;
1015 if(mcPart->IsPrimary())
1018 Int_t mother = mcPart->GetMother();
1020 mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
1027 //V0______________________________________
1028 //____________________________________________________________________
1029 TParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherV0(AliStack* stack, Int_t label)
1032 // Finds the first mother among the primary particles of the particle identified by <label>,
1033 // i.e. the primary that "caused" this particle
1035 // Taken from AliPWG0Helper class
1040 Int_t motherLabel = FindPrimaryMotherLabelV0(stack, label, nSteps);
1041 if (motherLabel < 0)
1044 return stack->Particle(motherLabel);
1047 //____________________________________________________________________
1048 Int_t AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherLabelV0(AliStack* stack, Int_t label, Int_t& nSteps)
1051 // Finds the first mother among the primary particles of the particle identified by <label>,
1052 // i.e. the primary that "caused" this particle
1054 // returns its label
1056 // Taken from AliPWG0Helper class
1059 const Int_t nPrim = stack->GetNprimary();
1061 while (label >= nPrim) {
1063 //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
1065 nSteps++; // 1 level down
1067 TParticle* particle = stack->Particle(label);
1070 AliDebugGeneral("FindPrimaryMotherLabelV0", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
1075 if (particle->GetMother(0) < 0) {
1077 AliDebugGeneral("FindPrimaryMotherLabelV0", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
1081 label = particle->GetMother(0);
1087 //____________________________________________________________________
1088 AliAODMCParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherAODV0(AliAODMCParticle* startParticle, Int_t& nSteps)
1091 // Finds the first mother among the primary particles of the particle identified by <label>,
1092 // i.e. the primary that "caused" this particle
1094 // Taken from AliPWG0Helper class
1099 AliAODMCParticle* mcPart = startParticle;
1104 if(mcPart->IsPrimary())
1107 Int_t mother = mcPart->GetMother();
1109 mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
1110 nSteps++; // 1 level down
1118 //__________________________________________________________________
1119 void AliAnalysisTaskQAHighPtDeDx::ProduceArrayTrksESD( AliESDEvent *ESDevent ){
1121 const Int_t nESDTracks = ESDevent->GetNumberOfTracks();
1122 //Int_t trackmult=0;
1127 //get multiplicity tpc only track cuts
1128 for(Int_t iT = 0; iT < nESDTracks; iT++) {
1130 AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
1133 if(TMath::Abs(esdTrack->Eta()) > fEtaCut)
1136 //only golden track cuts
1137 UInt_t selectDebug = 0;
1138 if (fTrackFilterTPC) {
1139 selectDebug = fTrackFilterTPC->IsSelected(esdTrack);
1141 //cout<<"this is not a golden track"<<endl;
1152 for(Int_t iT = 0; iT < nESDTracks; iT++) {
1154 AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
1157 if(TMath::Abs(esdTrack->Eta()) > fEtaCut)
1160 //only golden track cuts
1161 UInt_t selectDebug = 0;
1162 if (fTrackFilterGolden) {
1163 selectDebug = fTrackFilterGolden->IsSelected(esdTrack);
1165 //cout<<"this is not a golden track"<<endl;
1171 //Int_t sharedtpcclusters = esdTrack->GetTPCnclsS();
1172 Short_t ncl = esdTrack->GetTPCsignalN();
1177 Double_t eta = esdTrack->Eta();
1178 Double_t phi = esdTrack->Phi();
1179 Double_t momentum = esdTrack->P();
1180 Float_t dedx = esdTrack->GetTPCsignal();
1182 //cout<<"magf="<<magf<<endl;
1184 if(!PhiCut(esdTrack->Pt(), phi, esdTrack->Charge(), magf, cutLow, cutHigh))
1192 Bool_t IsTOFout=kFALSE;
1193 if ((esdTrack->GetStatus()&AliESDtrack::kTOFout)==0)
1195 Float_t lengthtrack=esdTrack->GetIntegratedLength();
1196 Float_t timeTOF=esdTrack->GetTOFsignal();
1197 Double_t inttime[5]={0,0,0,0,0};
1198 esdTrack->GetIntegratedTimes(inttime);// Returns the array with integrated times for each particle hypothesis
1201 if ( ( lengthtrack != 0 ) && ( timeTOF != 0) )
1202 beta = inttime[0] / timeTOF;
1205 Short_t pidCode = 0;
1209 const Int_t label = TMath::Abs(esdTrack->GetLabel());
1210 TParticle* mcTrack = fMCStack->Particle(label);
1213 Int_t pdgCode = mcTrack->GetPdgCode();
1214 pidCode = GetPidCode(pdgCode);
1221 if( momentum <= 0.6 && momentum >= 0.4 ){//only p:0.4-0.6 GeV, pion MIP
1223 if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
1224 if(momentum<0.6&&momentum>0.4){
1225 hMIPVsEta->Fill(eta,dedx);
1226 pMIPVsEta->Fill(eta,dedx);
1229 if( dedx > DeDxMIPMax+1 && dedx < 95 ){
1230 if(TMath::Abs(beta-1)<0.1){
1231 hPlateauVsEta->Fill(eta,dedx);
1232 pPlateauVsEta->Fill(eta,dedx);
1238 for(Int_t nh = 0; nh < 9; nh++) {
1240 if( eta > etaHigh[nh]/10.0 || eta < etaLow[nh]/10.0 )
1245 hMcOut[0][nh]->Fill(esdTrack->Pt());
1246 hMcOut[pidCode][nh]->Fill(esdTrack->Pt());
1249 histAllCh[nh]->Fill(momentum, dedx);
1251 histPiTof[nh]->Fill(momentum, dedx);
1252 histpPiTof[nh]->Fill(momentum);
1255 if( momentum <= 0.6 && momentum >= 0.4 ){//only p:0.4-0.6 GeV, pion MIP
1256 //Fill pion MIP, before calibration
1257 if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
1258 hMIPVsPhi[nh]->Fill(phi,dedx);
1259 pMIPVsPhi[nh]->Fill(phi,dedx);
1262 hMIPVsNch[nh]->Fill(multTPC,dedx);
1263 pMIPVsNch[nh]->Fill(multTPC,dedx);
1267 //Fill electrons, before calibration
1268 if( dedx > DeDxMIPMax+1 && dedx < 95 ){
1269 if(TMath::Abs(beta-1)<0.1){
1270 hPlateauVsPhi[nh]->Fill(phi,dedx);
1271 pPlateauVsPhi[nh]->Fill(phi,dedx);
1278 }//end loop over eta intervals
1284 }//end of track loop
1290 //__________________________________________________________________
1291 void AliAnalysisTaskQAHighPtDeDx::ProduceArrayTrksAOD( AliAODEvent *AODevent ){
1294 Int_t nAODTracks = AODevent->GetNumberOfTracks();
1297 //get multiplicity tpc only track cuts
1298 for(Int_t iT = 0; iT < nAODTracks; iT++) {
1300 AliAODTrack* aodTrack = AODevent->GetTrack(iT);
1302 if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
1306 if (fTrackFilterTPC) {
1307 // TPC only cuts is bit 1
1308 if(!aodTrack->TestFilterBit(1))
1317 for(Int_t iT = 0; iT < nAODTracks; iT++) {
1319 AliAODTrack* aodTrack = AODevent->GetTrack(iT);
1321 if (fTrackFilterGolden) {
1322 // "Global track RAA analysis QM2011 + Chi2ITS<36"; bit 1024
1323 if(!aodTrack->TestFilterBit(1024))
1328 if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
1332 Double_t eta = aodTrack->Eta();
1333 Double_t phi = aodTrack->Phi();
1334 Double_t momentum = aodTrack->P();
1337 if(!PhiCut(aodTrack->Pt(), phi, aodTrack->Charge(), magf, cutLow, cutHigh))
1342 AliAODPid* aodPid = aodTrack->GetDetPid();
1351 ncl = aodPid->GetTPCsignalN();
1352 dedx = aodPid->GetTPCsignal();
1354 Bool_t IsTOFout=kFALSE;
1355 Float_t lengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF
1356 Float_t timeTOF=-999;
1358 if ((aodTrack->GetStatus()&AliESDtrack::kTOFout)==0)
1361 lengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF
1363 timeTOF=aodPid->GetTOFsignal();
1365 Double_t inttime[5]={0,0,0,0,0};
1366 aodPid->GetIntegratedTimes(inttime);// Returns the array with integrated times for each particle hypothesis
1370 if ( ( lengthtrack != 0 ) && ( timeTOF != 0) )
1371 beta = inttime[0] / timeTOF;
1381 Short_t pidCode = 0;
1385 const Int_t label = TMath::Abs(aodTrack->GetLabel());
1386 AliAODMCParticle* mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(label));
1390 Int_t pdgCode = mcTrack->GetPdgCode();
1391 pidCode = GetPidCode(pdgCode);
1400 //if(momentum>0.6||momentum<0.4)//only p:0.4-0.6 GeV
1405 if( momentum <= 0.6 && momentum >= 0.4 ){//only p:0.4-0.6 GeV, pion MIP
1406 if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
1407 if(momentum<0.6&&momentum>0.4){
1408 hMIPVsEta->Fill(eta,dedx);
1409 pMIPVsEta->Fill(eta,dedx);
1412 if( dedx > DeDxMIPMax+1 && dedx < 95 ){
1413 if(TMath::Abs(beta-1)<0.1){
1414 hPlateauVsEta->Fill(eta,dedx);
1415 pPlateauVsEta->Fill(eta,dedx);
1421 for(Int_t nh = 0; nh < 9; nh++) {
1423 if( eta > etaHigh[nh]/10.0 || eta < etaLow[nh]/10.0 )
1428 hMcOut[0][nh]->Fill(aodTrack->Pt());
1429 hMcOut[pidCode][nh]->Fill(aodTrack->Pt());
1432 histAllCh[nh]->Fill(momentum, dedx);
1434 histPiTof[nh]->Fill(momentum, dedx);
1435 histpPiTof[nh]->Fill(momentum);
1438 if( momentum <= 0.6 && momentum >= 0.4 ){//only p:0.4-0.6 GeV, pion MIP
1439 //Fill pion MIP, before calibration
1440 if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
1441 hMIPVsPhi[nh]->Fill(phi,dedx);
1442 pMIPVsPhi[nh]->Fill(phi,dedx);
1445 hMIPVsNch[nh]->Fill(multTPC,dedx);
1446 pMIPVsNch[nh]->Fill(multTPC,dedx);
1450 //Fill electrons, before calibration
1451 if( dedx > DeDxMIPMax+1 && dedx < 95 ){
1452 if(TMath::Abs(beta-1)<0.1){
1453 hPlateauVsPhi[nh]->Fill(phi,dedx);
1454 pPlateauVsPhi[nh]->Fill(phi,dedx);
1460 }//end loop over eta intervals
1466 }//end of track loop
1475 //__________________________________________________
1476 void AliAnalysisTaskQAHighPtDeDx::ProduceArrayV0ESD( AliESDEvent *ESDevent ){
1477 Int_t nv0s = ESDevent->GetNumberOfV0s();
1483 const AliESDVertex *myBestPrimaryVertex = ESDevent->GetPrimaryVertex();
1484 if (!myBestPrimaryVertex) return;
1485 if (!(myBestPrimaryVertex->GetStatus())) return;
1486 Double_t lPrimaryVtxPosition[3];
1487 myBestPrimaryVertex->GetXYZ(lPrimaryVtxPosition);
1489 Double_t lPrimaryVtxCov[6];
1490 myBestPrimaryVertex->GetCovMatrix(lPrimaryVtxCov);
1491 Double_t lPrimaryVtxChi2 = myBestPrimaryVertex->GetChi2toNDF();
1493 AliAODVertex* myPrimaryVertex = new AliAODVertex(lPrimaryVtxPosition, lPrimaryVtxCov, lPrimaryVtxChi2, NULL, -1, AliAODVertex::kPrimary);
1498 // LOOP OVER V0s, K0s, L, AL
1502 for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
1504 // This is the begining of the V0 loop
1505 AliESDv0 *esdV0 = ESDevent->GetV0(iV0);
1506 if (!esdV0) continue;
1508 //check onfly status
1509 if(esdV0->GetOnFlyStatus()!=0)
1513 // AliESDTrack (V0 Daughters)
1514 UInt_t lKeyPos = (UInt_t)TMath::Abs(esdV0->GetPindex());
1515 UInt_t lKeyNeg = (UInt_t)TMath::Abs(esdV0->GetNindex());
1517 AliESDtrack *pTrack = ESDevent->GetTrack(lKeyPos);
1518 AliESDtrack *nTrack = ESDevent->GetTrack(lKeyNeg);
1519 if (!pTrack || !nTrack) {
1520 Printf("ERROR: Could not retreive one of the daughter track");
1525 if (pTrack->GetSign() == nTrack->GetSign()) {
1526 //cout<< "like sign, continue"<< endl;
1530 // Eta cut on decay products
1531 if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
1535 // Check if switch does anything!
1536 Bool_t isSwitched = kFALSE;
1537 if (pTrack->GetSign() < 0) { // switch
1540 AliESDtrack* helpTrack = nTrack;
1545 //Double_t lV0Position[3];
1546 //esdV0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[2]);
1549 AliKFVertex primaryVtxKF( *myPrimaryVertex );
1550 AliKFParticle::SetField(ESDevent->GetMagneticField());
1552 // Also implement switch here!!!!!!
1553 AliKFParticle* negEKF = 0; // e-
1554 AliKFParticle* posEKF = 0; // e+
1555 AliKFParticle* negPiKF = 0; // pi -
1556 AliKFParticle* posPiKF = 0; // pi +
1557 AliKFParticle* posPKF = 0; // p
1558 AliKFParticle* negAPKF = 0; // p-bar
1561 negEKF = new AliKFParticle( *(esdV0->GetParamN()) , 11);
1562 posEKF = new AliKFParticle( *(esdV0->GetParamP()) ,-11);
1563 negPiKF = new AliKFParticle( *(esdV0->GetParamN()) ,-211);
1564 posPiKF = new AliKFParticle( *(esdV0->GetParamP()) , 211);
1565 posPKF = new AliKFParticle( *(esdV0->GetParamP()) , 2212);
1566 negAPKF = new AliKFParticle( *(esdV0->GetParamN()) ,-2212);
1567 } else { // switch + and -
1568 negEKF = new AliKFParticle( *(esdV0->GetParamP()) , 11);
1569 posEKF = new AliKFParticle( *(esdV0->GetParamN()) ,-11);
1570 negPiKF = new AliKFParticle( *(esdV0->GetParamP()) ,-211);
1571 posPiKF = new AliKFParticle( *(esdV0->GetParamN()) , 211);
1572 posPKF = new AliKFParticle( *(esdV0->GetParamN()) , 2212);
1573 negAPKF = new AliKFParticle( *(esdV0->GetParamP()) ,-2212);
1576 AliKFParticle v0GKF; // Gamma e.g. from pi0
1579 v0GKF.SetProductionVertex(primaryVtxKF);
1581 AliKFParticle v0K0sKF; // K0 short
1582 v0K0sKF+=(*negPiKF);
1583 v0K0sKF+=(*posPiKF);
1584 v0K0sKF.SetProductionVertex(primaryVtxKF);
1586 AliKFParticle v0LambdaKF; // Lambda
1587 v0LambdaKF+=(*negPiKF);
1588 v0LambdaKF+=(*posPKF);
1589 v0LambdaKF.SetProductionVertex(primaryVtxKF);
1591 AliKFParticle v0AntiLambdaKF; // Lambda-bar
1592 v0AntiLambdaKF+=(*posPiKF);
1593 v0AntiLambdaKF+=(*negAPKF);
1594 v0AntiLambdaKF.SetProductionVertex(primaryVtxKF);
1596 Double_t dmassG = v0GKF.GetMass();
1597 Double_t dmassK = v0K0sKF.GetMass()-0.498;
1598 Double_t dmassL = v0LambdaKF.GetMass()-1.116;
1599 Double_t dmassAL = v0AntiLambdaKF.GetMass()-1.116;
1602 for( Int_t case_v0 = 0; case_v0 < 2; ++case_v0 ){
1608 Bool_t fillPos = kFALSE;
1609 Bool_t fillNeg = kFALSE;
1614 if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01){
1628 if(fillPos||fillNeg)
1635 for(Int_t j = 0; j < 2; j++) {
1637 AliESDtrack* track = 0;
1653 if(track->GetTPCsignalN()<=70)continue;
1654 Double_t phi = track->Phi();
1656 if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
1660 //if(!PhiCut(pt, phi, charge, magf, cutLow, cutHigh, hPhi))
1663 Double_t eta = track->Eta();
1664 Double_t momentum = track->Pt();
1665 Double_t dedx = track->GetTPCsignal();
1667 if(fillPos&&fillNeg){
1670 if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
1671 if(momentum<0.6&&momentum>0.4){
1672 hMIPVsEtaV0s->Fill(eta,dedx);
1673 pMIPVsEtaV0s->Fill(eta,dedx);
1680 for(Int_t nh = 0; nh < nHists; nh++) {
1684 if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
1687 if(fillPos&&fillNeg){
1689 histPiV0[nh]->Fill(momentum, dedx);
1690 histpPiV0[nh]->Fill(momentum);
1695 histPV0[nh]->Fill(momentum, dedx);
1696 histpPV0[nh]->Fill(momentum);
1702 }//end loop over two tracks
1709 Bool_t fillPos = kFALSE;
1710 Bool_t fillNeg = kFALSE;
1714 if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01) {
1715 if(dmassG<0.01 && dmassG>0.0001) {
1718 if( TMath::Abs(nTrack->GetTPCsignal() - 85.0) < 5)
1720 if( TMath::Abs(pTrack->GetTPCsignal() - 85.0) < 5)
1728 //cout<<"fillPos="<<fillPos<<endl;
1730 if(fillPos == kTRUE && fillNeg == kTRUE)
1734 AliESDtrack* track = 0;
1742 Double_t dedx = track->GetTPCsignal();
1743 Double_t eta = track->Eta();
1744 Double_t phi = track->Phi();
1745 Double_t momentum = track->P();
1747 if(track->GetTPCsignalN()<=70)continue;
1749 if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
1752 for(Int_t nh = 0; nh < nHists; nh++) {
1754 if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
1757 histEV0[nh]->Fill(momentum, dedx);
1767 }//end loop over case V0
1770 // clean up loop over v0
1782 delete myPrimaryVertex;
1786 //__________________________________________________________________________
1787 void AliAnalysisTaskQAHighPtDeDx::ProduceArrayV0AOD( AliAODEvent *AODevent ){
1788 Int_t nv0s = AODevent->GetNumberOfV0s();
1794 AliAODVertex *myBestPrimaryVertex = AODevent->GetPrimaryVertex();
1795 if (!myBestPrimaryVertex) return;
1799 // ################################
1800 // #### BEGINNING OF V0 CODE ######
1801 // ################################
1802 // This is the begining of the V0 loop
1803 for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
1804 AliAODv0 *aodV0 = AODevent->GetV0(iV0);
1805 if (!aodV0) continue;
1808 //check onfly status
1809 if(aodV0->GetOnFlyStatus()!=0)
1812 // AliAODTrack (V0 Daughters)
1813 AliAODVertex* vertex = aodV0->GetSecondaryVtx();
1815 Printf("ERROR: Could not retrieve vertex");
1819 AliAODTrack *pTrack = (AliAODTrack*)vertex->GetDaughter(0);
1820 AliAODTrack *nTrack = (AliAODTrack*)vertex->GetDaughter(1);
1821 if (!pTrack || !nTrack) {
1822 Printf("ERROR: Could not retrieve one of the daughter track");
1827 if (pTrack->Charge() == nTrack->Charge()) {
1828 //cout<< "like sign, continue"<< endl;
1832 // Make sure charge ordering is ok
1833 if (pTrack->Charge() < 0) {
1834 AliAODTrack* helpTrack = pTrack;
1839 // Eta cut on decay products
1840 if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
1844 Double_t dmassG = aodV0->InvMass2Prongs(0,1,11,11);
1845 Double_t dmassK = aodV0->MassK0Short()-0.498;
1846 Double_t dmassL = aodV0->MassLambda()-1.116;
1847 Double_t dmassAL = aodV0->MassAntiLambda()-1.116;
1849 for( Int_t case_v0 = 0; case_v0 < 2; ++case_v0 ){
1854 Bool_t fillPos = kFALSE;
1855 Bool_t fillNeg = kFALSE;
1861 if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01){
1875 if(fillPos||fillNeg)
1883 for(Int_t j = 0; j < 2; j++) {
1885 AliAODTrack* track = 0;
1901 if(track->GetTPCsignalN()<=70)continue;
1903 Double_t phi = track->Phi();
1905 if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
1909 //if(!PhiCut(pt, phi, charge, magf, cutLow, cutHigh, hPhi))
1912 Double_t eta = track->Eta();
1913 Double_t momentum = track->Pt();
1914 Double_t dedx = track->GetTPCsignal();
1916 if(fillPos&&fillNeg){
1919 if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){
1920 if(momentum<0.6&&momentum>0.4){
1921 hMIPVsEtaV0s->Fill(eta,dedx);
1922 pMIPVsEtaV0s->Fill(eta,dedx);
1929 for(Int_t nh = 0; nh < nHists; nh++) {
1933 if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
1936 if(fillPos&&fillNeg){
1938 histPiV0[nh]->Fill(momentum, dedx);
1939 histpPiV0[nh]->Fill(momentum);
1944 histPV0[nh]->Fill(momentum, dedx);
1945 histpPV0[nh]->Fill(momentum);
1952 }//end loop over two tracks
1958 Bool_t fillPos = kFALSE;
1959 Bool_t fillNeg = kFALSE;
1961 if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01) {
1962 if(dmassG<0.01 && dmassG>0.0001) {
1964 if( TMath::Abs(nTrack->GetTPCsignal() - 85.0) < 5)
1966 if( TMath::Abs(pTrack->GetTPCsignal() - 85.0) < 5)
1975 if(fillPos == kTRUE && fillNeg == kTRUE)
1979 AliAODTrack* track = 0;
1987 Double_t dedx = track->GetTPCsignal();
1988 Double_t eta = track->Eta();
1989 Double_t phi = track->Phi();
1990 Double_t momentum = track->P();
1992 if(track->GetTPCsignalN()<=70)continue;
1994 if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
1997 for(Int_t nh = 0; nh < nHists; nh++) {
1999 if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
2002 histEV0[nh]->Fill(momentum, dedx);
2011 }//end loop over V0s cases
2013 }//end loop over v0's
2019 Bool_t AliAnalysisTaskQAHighPtDeDx::PhiCut(Double_t pt, Double_t phi, Double_t q, Float_t mag, TF1* phiCutLow, TF1* phiCutHigh)
2024 //Double_t phi = track->Phi();
2025 if(mag < 0) // for negatve polarity field
2026 phi = TMath::TwoPi() - phi;
2027 if(q < 0) // for negatve charge
2028 phi = TMath::TwoPi()-phi;
2030 phi += TMath::Pi()/18.0; // to center gap in the middle
2031 phi = fmod(phi, TMath::Pi()/9.0);
2033 if(phi<phiCutHigh->Eval(pt)
2034 && phi>phiCutLow->Eval(pt))
2035 return kFALSE; // reject track
2037 hPhi->Fill(pt, phi);