First preliminary tasks for PWGLF production QA (cont'd)
[u/mrichter/AliRoot.git] / PWGLF / QATasks / AliAnalysisTaskQAHighPtDeDx.cxx
1  /*************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3 *                                                                        *
4 * Author: The ALICE Off-line Project.                                    *
5 * Contributors are mentioned in the code where appropriate.              *
6 *                                                                        *
7 * Permission to use, copy, modify and distribute this software and its   *
8 * documentation strictly for non-commercial purposes is hereby granted   *
9 * without fee, provided that the above copyright notice appears in all   *
10 * copies and that both the copyright notice and this permission notice   *
11 * appear in the supporting documentation. The authors make no claims     *
12 * about the suitability of this software for any purpose. It is          *
13 * provided "as is" without express or implied warranty.                  * **************************************************************************/
14
15
16 #include "AliAnalysisTaskQAHighPtDeDx.h"
17
18 // ROOT includes
19 #include <TList.h>
20 #include <TTree.h>
21 #include <TMath.h>
22 #include <TH1.h>
23 #include <TF1.h>
24 #include <TProfile.h>
25 #include <TParticle.h>
26 #include <TFile.h>
27
28 // AliRoot includes
29 #include <AliAnalysisManager.h>
30 #include <AliAnalysisFilter.h>
31 #include <AliESDInputHandler.h>
32 #include <AliESDEvent.h>
33 #include <AliESDVertex.h>
34 #include <AliLog.h>
35 #include <AliExternalTrackParam.h>
36 #include <AliESDtrackCuts.h>
37 #include <AliESDVZERO.h>
38 #include <AliAODVZERO.h>
39
40 #include <AliMCEventHandler.h>
41 #include <AliMCEvent.h>
42 #include <AliStack.h>
43
44 #include <TTreeStream.h>
45
46 #include <AliHeader.h>
47 #include <AliGenPythiaEventHeader.h>
48 #include <AliGenDPMjetEventHeader.h>
49
50 #include "AliCentrality.h" 
51 #include <AliESDv0.h>
52 #include <AliKFVertex.h>
53 #include <AliAODVertex.h>
54
55 #include <AliAODTrack.h> 
56 #include <AliAODPid.h> 
57 #include <AliAODMCHeader.h> 
58
59
60 // STL includes
61 #include <iostream>
62 using namespace std;
63
64
65 //
66 // Responsible:
67 // Antonio Ortiz (Lund)
68 // Peter Christiansen (Lund)
69 //
70
71
72
73
74
75
76
77
78 const Double_t AliAnalysisTaskQAHighPtDeDx::fgkClight = 2.99792458e-2;
79 Float_t magf = -1;
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
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};
88
89 Int_t nDeltaPiBins   = 60;
90 Double_t deltaPiLow  = 40;
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   AliAnalysisTaskSE(name),
97   fESD(0x0),
98   fAOD(0x0),
99   fMC(0x0),
100   fMCStack(0x0),
101   fMCArray(0x0),
102   fTrackFilterGolden(0x0),
103   fTrackFilterTPC(0x0),
104   fCentEst("V0M"),
105   fAnalysisType("ESD"),
106   fAnalysisMC(kFALSE),
107   fAnalysisPbPb(kFALSE),
108   ftrigBit(0x0),
109   fRandom(0x0),
110   fPileUpRej(kFALSE),
111   fVtxCut(10.0),  
112   fEtaCut(0.9),  
113   fMinCent(0.0),
114   fMaxCent(100.0),
115   fStoreMcIn(kFALSE),//
116   fMcProcessType(-999),
117   fTriggeredEventMB(-999),
118   fVtxStatus(-999),
119   fZvtx(-999),
120   fZvtxMC(-999),
121   fRun(-999),
122   fEventId(-999),
123   fListOfObjects(0), 
124   fEvents(0x0), fVtx(0x0), fVtxMC(0x0), fVtxBeforeCuts(0x0), fVtxAfterCuts(0x0),
125   fn1(0x0),
126   fcent(0x0),
127   hMIPVsEta(0x0),
128   pMIPVsEta(0x0),
129   hMIPVsEtaV0s(0x0),
130   pMIPVsEtaV0s(0x0),
131   hPlateauVsEta(0x0),
132   pPlateauVsEta(0x0),
133   hPhi(0x0)
134
135
136 {
137   // Default constructor (should not be used)
138   for(Int_t i=0;i<9;++i){
139  
140     hMIPVsNch[i]=0;//TH2D, MIP vs Nch for different eta intervals
141     pMIPVsNch[i]=0;//TProfile, MIP vs Nch for different eta intervals
142     hMIPVsPhi[i]=0;//TH2D, MIP vs phi for different eta intervals
143     pMIPVsPhi[i]=0;//TProfile, MIP vs phi for different eta intervals
144     hPlateauVsPhi[i]=0;//TH2D, dE/dx vs Phi, electrons 0.4<p<0.6 GeV/c
145     pPlateauVsPhi[i]=0;//TProfile, dE/dx vs Phi, electrons 0.4<p<0.6 GeV/c
146     histPiV0[i]=0;//TH2D, dE/dx vs p, pi id by V0s
147     histPV0[i]=0;// TH2D, dE/dx vs p, p id by V0s
148     histAllCh[i]=0;//TH2D, dE/dx vs p for all charged particles
149     histPiTof[i]=0;//TH2D, dE/dx vs p for a "clean" sample of pions, beta>1
150     histEV0[i]=0;
151
152     for(Int_t pid=0;pid<7;++pid){
153       hMcIn[pid][i]=0;
154       hMcOut[pid][i]=0;
155     }
156
157   }
158   DefineOutput(1, TList::Class());//esto es nuevo
159 }
160
161 //______________________________________________________________________________
162 void AliAnalysisTaskQAHighPtDeDx::UserCreateOutputObjects()
163
164   // This method is called once per worker node
165   // Here we define the output: histograms and debug tree if requested 
166   // We also create the random generator here so it might get different seeds...
167   fRandom = new TRandom(0); // 0 means random seed
168  
169
170
171   //OpenFile(1);
172   fListOfObjects = new TList();
173   fListOfObjects->SetOwner();
174   
175
176
177
178   //
179   // Histograms
180   //  
181   fEvents = new TH1I("fEvents","Number of analyzed events; Events; Counts", 3, 0, 3);
182   fListOfObjects->Add(fEvents);
183
184   fn1=new TH1F("fn1","fn1",11,-1,10);
185   fListOfObjects->Add(fn1);
186
187   fcent=new TH1F("fcent","fcent",104,-2,102);
188   fListOfObjects->Add(fcent);
189
190   fVtx = new TH1I("fVtx","Vtx info (0=no, 1=yes); Vtx; Counts", 2, -0.5, 1.5);
191   fListOfObjects->Add(fVtx);
192
193   fVtxBeforeCuts = new TH1F("fVtxBeforeCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
194   fListOfObjects->Add(fVtxBeforeCuts);
195   
196   fVtxAfterCuts = new TH1F("fVtxAfterCuts", "Vtx distribution (before cuts); Vtx z [cm]; Counts", 120, -30, 30);
197   fListOfObjects->Add(fVtxAfterCuts);
198
199
200   const Int_t nPtBinsV0s = 25;
201   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 ,
202                                        1.2 , 1.4 , 1.6 , 1.8 , 2.0 , 2.5 , 3.0 , 3.5 , 4.0 , 5.0 , 7.0 ,
203                                        9.0 , 12.0, 15.0, 20.0 };
204   
205
206
207
208   const Char_t* ending[nHists] = {"", "86", "64", "42", "20", "02", "24", "46", "68"};
209   
210   const Char_t* LatexEta[nHists] = {
211     "|#eta|<0.8", "-0.8<#eta<-0.6", "-0.6<#eta<-0.4", "-0.4<#eta<-0.2", "-0.2<#eta<0", 
212     "0<#eta<0.2", "0.2<#eta<0.4", "0.4<#eta<0.6", "0.6<#eta<0.8" 
213   };
214   hPhi = new TH2D("histPhi", "pt; #phi'", nPtBinsV0s, ptBinsV0s, 90, -0.05, 0.4); 
215   //dE/dx vs phi, pions at the MIP
216   fListOfObjects->Add(hPhi);
217
218
219   
220
221   Int_t nPhiBins = 36;
222   
223   for(Int_t i = 0; i < nHists; i++) {
224     
225    
226     hMIPVsPhi[i]      = new TH2D(Form("hMIPVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx MIP",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(), 
227                                  DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
228     hMIPVsPhi[i]->Sumw2();
229     
230     pMIPVsPhi[i]     = new TProfile(Form("pMIPVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx MIP",LatexEta[i]),  nPhiBins, 0, 2*TMath::Pi(),
231                                     DeDxMIPMin, DeDxMIPMax);
232     pMIPVsPhi[i]->SetMarkerStyle(20);
233     pMIPVsPhi[i]->SetMarkerColor(1);
234     pMIPVsPhi[i]->SetLineColor(1);
235     pMIPVsPhi[i]->Sumw2();
236
237     fListOfObjects->Add(hMIPVsPhi[i]);
238     fListOfObjects->Add(pMIPVsPhi[i]);
239
240     hPlateauVsPhi[i]  = new TH2D(Form("hPlateauVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx Plateau",LatexEta[i]),  nPhiBins, 0, 2*TMath::Pi(),
241                                   95-DeDxMIPMax, DeDxMIPMax, 95);
242     hPlateauVsPhi[i]->Sumw2();
243     
244     pPlateauVsPhi[i] = new TProfile(Form("pPlateauVsPhi%s", ending[i]), Form("%s; #phi (rad); dE/dx Plateau",LatexEta[i]), nPhiBins, 0, 2*TMath::Pi(),
245                                      DeDxMIPMax, 95);
246     pPlateauVsPhi[i]->SetMarkerStyle(20);
247     pPlateauVsPhi[i]->SetMarkerColor(1);
248     pPlateauVsPhi[i]->SetLineColor(1);
249     pPlateauVsPhi[i]->Sumw2();
250
251     fListOfObjects->Add(hPlateauVsPhi[i]);
252     fListOfObjects->Add(pPlateauVsPhi[i]);
253
254
255     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, 
256                                   DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
257     hMIPVsNch[i]->Sumw2();
258     
259     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);
260     pMIPVsNch[i]->SetMarkerStyle(20);
261     pMIPVsNch[i]->SetMarkerColor(1);
262     pMIPVsNch[i]->SetLineColor(1);
263     pMIPVsNch[i]->Sumw2();
264
265     fListOfObjects->Add(hMIPVsNch[i]);
266     fListOfObjects->Add(pMIPVsNch[i]);
267
268
269     histPiV0[i]  = new TH2D(Form("histPiV0%s", ending[i]), "Pions id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
270     histPiV0[i]->Sumw2();
271     histPV0[i]   = new TH2D(Form("histPV0%s", ending[i]), "Protons id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
272
273     fListOfObjects->Add(histPiV0[i]);
274     fListOfObjects->Add(histPV0[i]);
275
276     histPiTof[i] = new TH2D(Form("histPiTof%s", ending[i]), "all charged", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
277     histPiTof[i]->Sumw2();
278
279     histAllCh[i] = new TH2D(Form("histAllCh%s", ending[i]), "Pions id by TOF", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
280     histAllCh[i]->Sumw2();
281
282     fListOfObjects->Add(histPiTof[i]);
283     fListOfObjects->Add(histAllCh[i]);
284
285     histEV0[i]   = new TH2D(Form("histEV0%s", ending[i]), "Electrons id by V0", nPtBinsV0s, ptBinsV0s, nDeltaPiBins, deltaPiLow, deltaPiHigh);
286     histEV0[i]->Sumw2();
287     fListOfObjects->Add(histEV0[i]);
288
289
290
291   }
292
293
294   hMIPVsEta = new TH2D("hMIPVsEta","; #eta; dE/dx_{MIP, primary tracks}",16,-0.8,0.8,DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
295   pMIPVsEta = new TProfile("pMIPVsEta","; #eta; #LT dE/dx #GT_{MIP, primary tracks}",16,-0.8,0.8, DeDxMIPMin, DeDxMIPMax);
296   hMIPVsEtaV0s = new TH2D("hMIPVsEtaV0s","; #eta; dE/dx_{MIP, secondary tracks}",16,-0.8,0.8,DeDxMIPMax-DeDxMIPMin, DeDxMIPMin, DeDxMIPMax);
297   pMIPVsEtaV0s = new TProfile("pMIPVsEtaV0s","; #eta; #LT dE/dx #GT_{MIP, secondary tracks}",16,-0.8,0.8,DeDxMIPMin, DeDxMIPMax);
298   
299   hPlateauVsEta = new TH2D("hPlateauVsEta","; #eta; dE/dx_{Plateau, primary tracks}",16,-0.8,0.8,95-DeDxMIPMax, DeDxMIPMax, 95);
300   pPlateauVsEta = new TProfile("pPlateauVsEta","; #eta; #LT dE/dx #GT_{Plateau, primary tracks}",16,-0.8,0.8, DeDxMIPMax, 95);
301
302   fListOfObjects->Add(hMIPVsEta);
303   fListOfObjects->Add(pMIPVsEta);
304   fListOfObjects->Add(hMIPVsEtaV0s);
305   fListOfObjects->Add(pMIPVsEtaV0s);
306   fListOfObjects->Add(hPlateauVsEta);
307   fListOfObjects->Add(pPlateauVsEta);
308
309
310
311
312
313
314   if (fAnalysisMC) {    
315     for(Int_t i = 0; i < nHists; i++) {
316       for(Int_t pid = 0; pid < 7; pid++) {
317         
318         hMcIn[pid][i] = new TH1D(Form("hIn%s%s", Pid[pid],ending[i]), Form("MC in (pid %s)", Pid[pid]), 
319                                  nPtBinsV0s, ptBinsV0s);
320         fListOfObjects->Add(hMcIn[pid][i]);
321         
322         hMcOut[pid][i] = new TH1D(Form("hMcOut%s%s", Pid[pid],ending[i]), Form("MC out (pid %s)", Pid[pid]), 
323                                   nPtBinsV0s, ptBinsV0s);
324         fListOfObjects->Add(hMcOut[pid][i]);
325         
326         
327       }
328     }
329
330     fVtxMC = new TH1F("fVtxMC","mc vtx", 120, -30, 30);
331     fListOfObjects->Add(fVtxMC);
332     
333   }
334  
335   // Post output data.
336   PostData(1, fListOfObjects);
337 }
338
339 //______________________________________________________________________________
340 void AliAnalysisTaskQAHighPtDeDx::UserExec(Option_t *) 
341 {
342   // Main loop
343
344   //
345   // First we make sure that we have valid input(s)!
346   //
347
348
349
350   AliVEvent *event = InputEvent();
351   if (!event) {
352      Error("UserExec", "Could not retrieve event");
353      return;
354   }
355
356
357
358   if (fAnalysisType == "ESD"){
359     fESD = dynamic_cast<AliESDEvent*>(event);
360     if(!fESD){
361       Printf("%s:%d ESDEvent not found in Input Manager",(char*)__FILE__,__LINE__);
362       this->Dump();
363       return;
364     }    
365   } else {
366     fAOD = dynamic_cast<AliAODEvent*>(event);
367     if(!fAOD){
368       Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
369       this->Dump();
370       return;
371     }    
372   }
373
374
375
376   if (fAnalysisMC) {
377
378     if (fAnalysisType == "ESD"){
379       fMC = dynamic_cast<AliMCEvent*>(MCEvent());
380       if(!fMC){
381         Printf("%s:%d MCEvent not found in Input Manager",(char*)__FILE__,__LINE__);
382         this->Dump();
383         return;
384       }    
385
386       fMCStack = fMC->Stack();
387       
388       if(!fMCStack){
389         Printf("%s:%d MCStack not found in Input Manager",(char*)__FILE__,__LINE__);
390         this->Dump();
391         return;
392       }    
393     } else { // AOD
394
395       fMC = dynamic_cast<AliMCEvent*>(MCEvent());
396       if(fMC)
397         fMC->Dump();
398
399       fMCArray = (TClonesArray*)fAOD->FindListObject("mcparticles");
400       if(!fMCArray){
401         Printf("%s:%d AOD MC array not found in Input Manager",(char*)__FILE__,__LINE__);
402         this->Dump();
403         return;
404       }    
405     }
406   }
407
408   
409   // Get trigger decision
410   fTriggeredEventMB = 0; //init
411   
412   fn1->Fill(0);
413   
414   if(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))
415      ->IsEventSelected() & ftrigBit ){
416     fTriggeredEventMB = 1;  //event triggered as minimum bias
417   }
418
419   // Get process type for MC
420   fMcProcessType = 0; // -1=invalid, 0=data, 1=ND, 2=SD, 3=DD
421
422   // real data that are not triggered we skip
423   if(!fAnalysisMC && !fTriggeredEventMB)
424     return; 
425   
426   fn1->Fill(1);
427
428
429   if (fAnalysisMC) {
430     
431
432
433     if (fAnalysisType == "ESD"){
434
435   
436
437       AliHeader* headerMC = fMC->Header();
438       if (headerMC) {
439         
440         AliGenEventHeader* genHeader = headerMC->GenEventHeader();
441         TArrayF vtxMC(3); // primary vertex  MC 
442         vtxMC[0]=9999; vtxMC[1]=9999;  vtxMC[2]=9999; //initialize with dummy
443         if (genHeader) {
444           genHeader->PrimaryVertex(vtxMC);
445         }
446         fZvtxMC = vtxMC[2];
447         
448         // PYTHIA:
449         AliGenPythiaEventHeader* pythiaGenHeader =
450           dynamic_cast<AliGenPythiaEventHeader*>(headerMC->GenEventHeader());
451         if (pythiaGenHeader) {  //works only for pythia
452           fMcProcessType =  GetPythiaEventProcessType(pythiaGenHeader->ProcessType());
453         }
454         // PHOJET:
455         AliGenDPMjetEventHeader* dpmJetGenHeader =
456           dynamic_cast<AliGenDPMjetEventHeader*>(headerMC->GenEventHeader());
457         if (dpmJetGenHeader) {
458           fMcProcessType = GetDPMjetEventProcessType(dpmJetGenHeader->ProcessType());
459         }
460       }
461     } else { // AOD
462       
463   
464
465       AliAODMCHeader* mcHeader = dynamic_cast<AliAODMCHeader*>(fAOD->FindListObject("mcHeader")); 
466
467
468       if(mcHeader) {
469         fZvtxMC = mcHeader->GetVtxZ();
470         
471
472
473         if(strstr(mcHeader->GetGeneratorName(), "Pythia")) {
474           fMcProcessType =  GetPythiaEventProcessType(mcHeader->GetEventType());
475         } else {
476           fMcProcessType =  GetDPMjetEventProcessType(mcHeader->GetEventType());
477         }
478       }
479     }
480
481
482   }
483   
484   
485
486   if (fAnalysisType == "ESD"){
487     
488     const AliESDVertex *vtxESD = fESD->GetPrimaryVertexTracks();
489     if(vtxESD->GetNContributors()<1) {
490       // SPD vertex
491       vtxESD = fESD->GetPrimaryVertexSPD();
492       /* quality checks on SPD-vertex */
493       TString vertexType = vtxESD->GetTitle();
494       if (vertexType.Contains("vertexer: Z") && (vtxESD->GetDispersion() > 0.04 || vtxESD->GetZRes() > 0.25))  
495         fZvtx  = -1599; //vertex = 0x0; //
496       else if (vtxESD->GetNContributors()<1) 
497         fZvtx  = -999; //vertex = 0x0; //
498       else
499         fZvtx = vtxESD->GetZ();
500     }  
501     else
502       fZvtx = vtxESD->GetZ();
503     
504   }
505   else // AOD
506     fZvtx = GetVertex(fAOD);
507     
508   fVtxBeforeCuts->Fill(fZvtx);
509   
510   //cut on the z position of vertex
511   if (TMath::Abs(fZvtx) > fVtxCut) {    
512     return;
513   }
514   fn1->Fill(2);
515   
516
517
518
519
520
521   Float_t centrality = -10;
522
523   // only analyze triggered events
524   if(fTriggeredEventMB) {
525     
526     if (fAnalysisType == "ESD"){
527       if(fAnalysisPbPb){
528         AliCentrality *centObject = fESD->GetCentrality();
529         centrality = centObject->GetCentralityPercentile(fCentEst);
530         
531         if((centrality>fMaxCent)||(centrality<fMinCent))return;
532       }
533       fcent->Fill(centrality);
534       fn1->Fill(3);
535       if(fAnalysisMC){
536         if(TMath::Abs(fZvtxMC)<fVtxCut){
537           ProcessMCTruthESD();
538           fVtxMC->Fill(fZvtxMC);
539         }
540       }
541       AnalyzeESD(fESD);
542     } else { // AOD
543       if(fAnalysisPbPb){
544         AliCentrality *centObject = fAOD->GetCentrality();
545         if(centObject){
546           centrality = centObject->GetCentralityPercentile(fCentEst); 
547         }
548         //cout<<"centrality="<<centrality<<endl;
549         if((centrality>fMaxCent)||(centrality<fMinCent))return;
550       }
551       fcent->Fill(centrality);
552       fn1->Fill(3);
553       if(fAnalysisMC){
554         if(TMath::Abs(fZvtxMC)<fVtxCut){
555
556           ProcessMCTruthAOD();
557           fVtxMC->Fill(fZvtxMC);
558         }
559       }
560       AnalyzeAOD(fAOD);
561     }
562   }
563
564   fVtxAfterCuts->Fill(fZvtx);
565
566   
567   
568
569   // Post output data.
570   PostData(1, fListOfObjects);
571 }
572
573 //________________________________________________________________________
574 void AliAnalysisTaskQAHighPtDeDx::AnalyzeESD(AliESDEvent* esdEvent)
575 {
576   fRun  = esdEvent->GetRunNumber();
577   fEventId = 0;
578   if(esdEvent->GetHeader())
579     fEventId = GetEventIdAsLong(esdEvent->GetHeader());
580   
581   Bool_t isPileup = esdEvent->IsPileupFromSPD();
582   if(fPileUpRej)
583     if(isPileup)
584       return;
585   fn1->Fill(4);
586
587
588   //  Int_t     event     = esdEvent->GetEventNumberInFile();
589   //UInt_t    time      = esdEvent->GetTimeStamp();
590   //  ULong64_t trigger   = esdEvent->GetTriggerMask();
591   magf      = esdEvent->GetMagneticField();
592
593
594
595
596
597   if(fTriggeredEventMB) {// Only MC case can we have not triggered events
598     
599     // accepted event
600     fEvents->Fill(0);
601     
602     
603     //Change, 10/04/13. Now accept all events to do a correct normalization
604     //if(fVtxStatus!=1) return; // accepted vertex
605     //Int_t nESDTracks = esdEvent->GetNumberOfTracks();
606     
607     ProduceArrayTrksESD( esdEvent );//produce array with global track parameters
608     ProduceArrayV0ESD( esdEvent );//v0's
609
610     
611     fEvents->Fill(1);
612
613
614
615
616   } // end if triggered
617
618
619 }
620
621 //________________________________________________________________________
622 void AliAnalysisTaskQAHighPtDeDx::AnalyzeAOD(AliAODEvent* aodEvent)
623 {
624   fRun  = aodEvent->GetRunNumber();
625   fEventId = 0;
626   if(aodEvent->GetHeader())
627     fEventId = GetEventIdAsLong(aodEvent->GetHeader());
628    
629   //UInt_t    time      = 0; // Missing AOD info? aodEvent->GetTimeStamp();
630   magf      = aodEvent->GetMagneticField();
631
632   //Int_t     trackmult = 0; // no pt cuts
633   //Int_t     nadded    = 0;
634
635   Bool_t isPileup = aodEvent->IsPileupFromSPD();
636   if(fPileUpRej)
637     if(isPileup)
638       return;
639   fn1->Fill(4);
640
641
642
643   if(fTriggeredEventMB) {// Only MC case can we have not triggered events
644     
645     // accepted event
646     fEvents->Fill(0);
647     
648     //if(fVtxStatus!=1) return; // accepted vertex
649     //Int_t nAODTracks = aodEvent->GetNumberOfTracks();  
650
651     ProduceArrayTrksAOD( aodEvent );
652     ProduceArrayV0AOD( aodEvent );
653
654     fEvents->Fill(1);
655  
656
657
658
659   } // end if triggered
660
661 }
662
663 //_____________________________________________________________________________
664 Float_t AliAnalysisTaskQAHighPtDeDx::GetVertex(const AliVEvent* event) const
665 {
666   Float_t zvtx = -999;
667   
668   const AliVVertex* primaryVertex = event->GetPrimaryVertex(); 
669   
670   if(primaryVertex->GetNContributors()>0)
671     zvtx = primaryVertex->GetZ();
672
673   return zvtx;
674 }
675
676 //_____________________________________________________________________________
677 Short_t AliAnalysisTaskQAHighPtDeDx::GetPidCode(Int_t pdgCode) const 
678 {
679   // return our internal code for pions, kaons, and protons
680   
681   Short_t pidCode = 6;
682   
683   switch (TMath::Abs(pdgCode)) {
684   case 211:
685     pidCode = 1; // pion
686     break;
687   case 321:
688     pidCode = 2; // kaon
689     break;
690   case 2212:
691     pidCode = 3; // proton
692     break;
693   case 11:
694     pidCode = 4; // electron
695     break;
696   case 13:
697     pidCode = 5; // muon
698     break;
699   default:
700     pidCode = 6;  // something else?
701   };
702   
703   return pidCode;
704 }
705
706 //_____________________________________________________________________________
707 void AliAnalysisTaskQAHighPtDeDx::ProcessMCTruthESD() 
708 {
709   // Fill the special MC histogram with the MC truth info
710   
711   const Int_t nTracksMC = fMCStack->GetNtrack();
712
713   for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
714     
715     //Cuts
716     if(!(fMCStack->IsPhysicalPrimary(iTracks)))
717       continue;
718     
719     TParticle* trackMC = fMCStack->Particle(iTracks);
720     
721     TParticlePDG* pdgPart = trackMC ->GetPDG();
722     Double_t chargeMC = pdgPart->Charge();
723
724     if(chargeMC==0)
725       continue;
726     
727     if (TMath::Abs(trackMC->Eta()) > fEtaCut )
728       continue;
729
730     Double_t etaMC = trackMC->Eta();
731     Int_t pdgCode = trackMC->GetPdgCode();
732     Short_t pidCodeMC = 0;
733     pidCodeMC = GetPidCode(pdgCode);
734     
735     
736     for(Int_t nh = 0; nh < 9; nh++) {
737       
738       if( etaMC > etaHigh[nh]/10.0 || etaMC < etaLow[nh]/10.0 )
739         continue;
740       
741       hMcIn[0][nh]->Fill(trackMC->Pt());
742       hMcIn[pidCodeMC][nh]->Fill(trackMC->Pt());
743       
744       
745     }
746
747   }//MC track loop
748   
749  
750   
751 }
752
753 //_____________________________________________________________________________
754 void AliAnalysisTaskQAHighPtDeDx::ProcessMCTruthAOD() 
755 {
756   // Fill the special MC histogram with the MC truth info
757   
758
759   const Int_t nTracksMC = fMCArray->GetEntriesFast();
760   
761   for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) {
762     
763     AliAODMCParticle* trackMC = dynamic_cast<AliAODMCParticle*>(fMCArray->At(iTracks));
764     
765     //Cuts
766     if(!(trackMC->IsPhysicalPrimary()))
767       continue;
768     
769    
770     Double_t chargeMC = trackMC->Charge();
771     if(chargeMC==0)
772       continue;
773
774     
775     if (TMath::Abs(trackMC->Eta()) > fEtaCut )
776       continue;
777     
778     Double_t etaMC = trackMC->Eta();
779     Int_t pdgCode = trackMC->GetPdgCode();
780     Short_t pidCodeMC = 0;
781     pidCodeMC = GetPidCode(pdgCode);
782     
783     //cout<<"pidcode="<<pidCodeMC<<endl;
784     for(Int_t nh = 0; nh < 9; nh++) {
785       
786       if( etaMC > etaHigh[nh]/10.0 || etaMC < etaLow[nh]/10.0 )
787         continue;
788       
789       hMcIn[0][nh]->Fill(trackMC->Pt());
790       hMcIn[pidCodeMC][nh]->Fill(trackMC->Pt());
791       
792       
793     }
794     
795   }//MC track loop
796   
797   
798 }
799
800
801 //_____________________________________________________________________________
802 Short_t AliAnalysisTaskQAHighPtDeDx::GetPythiaEventProcessType(Int_t pythiaType) {
803   //
804   // Get the process type of the event.  PYTHIA
805   //
806   // source PWG0   dNdpt 
807
808   Short_t globalType = -1; //init
809       
810   if(pythiaType==92||pythiaType==93){
811     globalType = 2; //single diffractive
812   }
813   else if(pythiaType==94){
814     globalType = 3; //double diffractive
815   }
816   //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB??
817   else {
818     globalType = 1;  //non diffractive
819   }
820   return globalType;
821 }
822
823 //_____________________________________________________________________________
824 Short_t AliAnalysisTaskQAHighPtDeDx::GetDPMjetEventProcessType(Int_t dpmJetType) {
825   //
826   // get the process type of the event.  PHOJET
827   //
828   //source PWG0   dNdpt 
829   // can only read pythia headers, either directly or from cocktalil header
830   Short_t globalType = -1;
831   
832   if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction
833     globalType = 1;
834   }
835   else if (dpmJetType==5 || dpmJetType==6) {
836     globalType = 2;
837   }
838   else if (dpmJetType==7) {
839     globalType = 3;
840   }
841   return globalType;
842 }
843
844 //_____________________________________________________________________________
845 ULong64_t AliAnalysisTaskQAHighPtDeDx::GetEventIdAsLong(AliVHeader* header) const
846 {
847   // To have a unique id for each event in a run!
848   // Modified from AliRawReader.h
849   return ((ULong64_t)header->GetBunchCrossNumber()+
850           (ULong64_t)header->GetOrbitNumber()*3564+
851           (ULong64_t)header->GetPeriodNumber()*16777215*3564);
852 }
853
854
855 //____________________________________________________________________
856 TParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMother(AliStack* stack, Int_t label)
857 {
858   //
859   // Finds the first mother among the primary particles of the particle identified by <label>,
860   // i.e. the primary that "caused" this particle
861   //
862   // Taken from AliPWG0Helper class
863   //
864
865   Int_t motherLabel = FindPrimaryMotherLabel(stack, label);
866   if (motherLabel < 0)
867     return 0;
868
869   return stack->Particle(motherLabel);
870 }
871
872 //____________________________________________________________________
873 Int_t AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherLabel(AliStack* stack, Int_t label)
874 {
875   //
876   // Finds the first mother among the primary particles of the particle identified by <label>,
877   // i.e. the primary that "caused" this particle
878   //
879   // returns its label
880   //
881   // Taken from AliPWG0Helper class
882   //
883   const Int_t nPrim  = stack->GetNprimary();
884   
885   while (label >= nPrim) {
886
887     //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
888
889     TParticle* particle = stack->Particle(label);
890     if (!particle) {
891       
892       AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
893       return -1;
894     }
895     
896     // find mother
897     if (particle->GetMother(0) < 0) {
898
899       AliDebugGeneral("FindPrimaryMotherLabel", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
900       return -1;
901     }
902     
903     label = particle->GetMother(0);
904   }
905   
906   return label;
907 }
908
909 //____________________________________________________________________
910 AliAODMCParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherAOD(AliAODMCParticle* startParticle)
911 {
912   //
913   // Finds the first mother among the primary particles of the particle identified by <label>,
914   // i.e. the primary that "caused" this particle
915   //
916   // Taken from AliPWG0Helper class
917   //
918
919   AliAODMCParticle* mcPart = startParticle;
920
921   while (mcPart)
922     {
923       
924       if(mcPart->IsPrimary())
925         return mcPart;
926       
927       Int_t mother = mcPart->GetMother();
928
929       mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
930     }
931
932   return 0;
933 }
934
935
936 //V0______________________________________
937 //____________________________________________________________________
938 TParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherV0(AliStack* stack, Int_t label)
939 {
940   //
941   // Finds the first mother among the primary particles of the particle identified by <label>,
942   // i.e. the primary that "caused" this particle
943   //
944   // Taken from AliPWG0Helper class
945   //
946
947   Int_t nSteps = 0;
948
949   Int_t motherLabel = FindPrimaryMotherLabelV0(stack, label, nSteps);
950   if (motherLabel < 0)
951     return 0;
952
953   return stack->Particle(motherLabel);
954 }
955
956 //____________________________________________________________________
957 Int_t AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherLabelV0(AliStack* stack, Int_t label, Int_t& nSteps)
958 {
959   //
960   // Finds the first mother among the primary particles of the particle identified by <label>,
961   // i.e. the primary that "caused" this particle
962   //
963   // returns its label
964   //
965   // Taken from AliPWG0Helper class
966   //
967   nSteps = 0;
968   const Int_t nPrim  = stack->GetNprimary();
969   
970   while (label >= nPrim) {
971
972     //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0));
973
974     nSteps++; // 1 level down
975     
976     TParticle* particle = stack->Particle(label);
977     if (!particle) {
978       
979       AliDebugGeneral("FindPrimaryMotherLabelV0", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label));
980       return -1;
981     }
982     
983     // find mother
984     if (particle->GetMother(0) < 0) {
985
986       AliDebugGeneral("FindPrimaryMotherLabelV0", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label));
987       return -1;
988     }
989     
990     label = particle->GetMother(0);
991   }
992   
993   return label;
994 }
995
996 //____________________________________________________________________
997 AliAODMCParticle* AliAnalysisTaskQAHighPtDeDx::FindPrimaryMotherAODV0(AliAODMCParticle* startParticle, Int_t& nSteps)
998 {
999   //
1000   // Finds the first mother among the primary particles of the particle identified by <label>,
1001   // i.e. the primary that "caused" this particle
1002   //
1003   // Taken from AliPWG0Helper class
1004   //
1005
1006   nSteps = 0;
1007
1008   AliAODMCParticle* mcPart = startParticle;
1009
1010   while (mcPart)
1011     {
1012       
1013       if(mcPart->IsPrimary())
1014         return mcPart;
1015       
1016       Int_t mother = mcPart->GetMother();
1017
1018       mcPart = dynamic_cast<AliAODMCParticle*>(fMCArray->At(mother));
1019       nSteps++; // 1 level down
1020     }
1021
1022   return 0;
1023 }
1024
1025
1026
1027 //__________________________________________________________________
1028 void AliAnalysisTaskQAHighPtDeDx::ProduceArrayTrksESD( AliESDEvent *ESDevent ){
1029   
1030   const Int_t nESDTracks = ESDevent->GetNumberOfTracks();
1031   //Int_t trackmult=0;
1032
1033
1034   Int_t multTPC = 0;
1035   
1036   //get multiplicity tpc only track cuts
1037   for(Int_t iT = 0; iT < nESDTracks; iT++) {
1038     
1039     AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
1040     
1041     
1042     if(TMath::Abs(esdTrack->Eta()) > fEtaCut)
1043       continue;
1044     
1045     //only golden track cuts
1046     UInt_t selectDebug = 0;
1047     if (fTrackFilterTPC) {
1048       selectDebug = fTrackFilterTPC->IsSelected(esdTrack);
1049       if (!selectDebug) {
1050         //cout<<"this is not a golden track"<<endl;
1051         continue;
1052       }
1053     }
1054    
1055     multTPC++;
1056     
1057   }
1058
1059
1060
1061   for(Int_t iT = 0; iT < nESDTracks; iT++) {
1062     
1063     AliESDtrack* esdTrack = ESDevent->GetTrack(iT);
1064     
1065       
1066     if(TMath::Abs(esdTrack->Eta()) > fEtaCut)
1067       continue;
1068     
1069     //only golden track cuts
1070     UInt_t selectDebug = 0;
1071     if (fTrackFilterGolden) {
1072       selectDebug = fTrackFilterGolden->IsSelected(esdTrack);
1073       if (!selectDebug) {
1074         //cout<<"this is not a golden track"<<endl;
1075         continue;
1076       }
1077     }
1078     
1079     
1080     //Int_t sharedtpcclusters = esdTrack->GetTPCnclsS();
1081     Short_t ncl     = esdTrack->GetTPCsignalN();
1082
1083
1084     if(ncl<70)
1085       continue;
1086     Double_t eta  = esdTrack->Eta();
1087     Double_t phi  = esdTrack->Phi();
1088     Double_t momentum = esdTrack->P();
1089     Float_t  dedx    = esdTrack->GetTPCsignal();
1090
1091     //cout<<"magf="<<magf<<endl;
1092
1093     if(!PhiCut(esdTrack->Pt(), phi, esdTrack->Charge(), magf, cutLow, cutHigh))
1094       continue;
1095
1096        
1097
1098
1099     //TOF
1100     
1101     Bool_t IsTOFout=kFALSE;
1102     if ((esdTrack->GetStatus()&AliESDtrack::kTOFout)==0)
1103       IsTOFout=kTRUE;
1104     Float_t lengthtrack=esdTrack->GetIntegratedLength();
1105     Float_t timeTOF=esdTrack->GetTOFsignal();
1106     Double_t inttime[5]={0,0,0,0,0}; 
1107     esdTrack->GetIntegratedTimes(inttime);// Returns the array with integrated times for each particle hypothesis
1108     Float_t beta = -99;
1109     if ( !IsTOFout ){
1110       if ( ( lengthtrack != 0 ) && ( timeTOF != 0) )
1111         beta = inttime[0] / timeTOF;
1112     }
1113
1114     Short_t pidCode     = 0; 
1115     
1116     if(fAnalysisMC) {
1117       
1118       const Int_t label = TMath::Abs(esdTrack->GetLabel());
1119       TParticle* mcTrack = fMCStack->Particle(label);       
1120       if (mcTrack){
1121         
1122         Int_t pdgCode = mcTrack->GetPdgCode();
1123         pidCode = GetPidCode(pdgCode);
1124         
1125       }
1126       
1127     }
1128
1129
1130     if( momentum <= 0.6 && momentum >= 0.4  ){//only p:0.4-0.6 GeV, pion MIP
1131       
1132       if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){       
1133         if(momentum<0.6&&momentum>0.4){
1134           hMIPVsEta->Fill(eta,dedx);
1135           pMIPVsEta->Fill(eta,dedx);
1136         }
1137       }
1138       if( dedx > DeDxMIPMax+5 && dedx < 95 ){
1139         if(TMath::Abs(beta-1)<0.1){
1140           hPlateauVsEta->Fill(eta,dedx);
1141           pPlateauVsEta->Fill(eta,dedx); 
1142         }
1143       }
1144     }
1145     
1146     
1147     for(Int_t nh = 0; nh < 9; nh++) {
1148       
1149       if( eta > etaHigh[nh]/10.0 || eta < etaLow[nh]/10.0 )
1150         continue;
1151       
1152   
1153       if(fAnalysisMC){
1154         hMcOut[0][nh]->Fill(esdTrack->Pt());
1155         hMcOut[pidCode][nh]->Fill(esdTrack->Pt());
1156       }
1157
1158       histAllCh[nh]->Fill(momentum, dedx);
1159       if(beta>1)
1160         histPiTof[nh]->Fill(momentum, dedx);
1161       
1162       if( momentum <= 0.6 && momentum >= 0.4  ){//only p:0.4-0.6 GeV, pion MIP
1163         //Fill  pion MIP, before calibration
1164         if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){     
1165           hMIPVsPhi[nh]->Fill(phi,dedx);
1166           pMIPVsPhi[nh]->Fill(phi,dedx);
1167           
1168           
1169           hMIPVsNch[nh]->Fill(multTPC,dedx);
1170           pMIPVsNch[nh]->Fill(multTPC,dedx);
1171           
1172         }
1173         
1174         //Fill electrons, before calibration
1175         if( dedx > DeDxMIPMax+5 && dedx < 95 ){
1176           if(TMath::Abs(beta-1)<0.1){
1177             hPlateauVsPhi[nh]->Fill(phi,dedx);
1178             pPlateauVsPhi[nh]->Fill(phi,dedx);
1179           }
1180         }
1181
1182       }
1183
1184
1185     }//end loop over eta intervals
1186
1187     
1188
1189     
1190     
1191   }//end of track loop
1192   
1193   
1194
1195   
1196 }
1197 //__________________________________________________________________
1198 void AliAnalysisTaskQAHighPtDeDx::ProduceArrayTrksAOD( AliAODEvent *AODevent ){
1199
1200
1201   Int_t nAODTracks = AODevent->GetNumberOfTracks();
1202   Int_t multTPC = 0;
1203   
1204   //get multiplicity tpc only track cuts
1205   for(Int_t iT = 0; iT < nAODTracks; iT++) {
1206     
1207     AliAODTrack* aodTrack = AODevent->GetTrack(iT);
1208     
1209     if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
1210       continue;
1211     
1212
1213     if (fTrackFilterTPC) {
1214       // TPC only cuts is bit 1
1215       if(!aodTrack->TestFilterBit(1))
1216         continue;
1217     }
1218    
1219     multTPC++;
1220     
1221   }
1222
1223
1224   for(Int_t iT = 0; iT < nAODTracks; iT++) {
1225     
1226     AliAODTrack* aodTrack = AODevent->GetTrack(iT);
1227     
1228     if (fTrackFilterGolden) {     
1229       // "Global track RAA analysis QM2011 + Chi2ITS<36"; bit 1024
1230       if(!aodTrack->TestFilterBit(1024))
1231         continue;
1232     }
1233     
1234     
1235     if(TMath::Abs(aodTrack->Eta()) > fEtaCut)
1236       continue;
1237     
1238     
1239     Double_t eta  = aodTrack->Eta();
1240     Double_t phi  = aodTrack->Phi();
1241     Double_t momentum = aodTrack->P();
1242     
1243  
1244     if(!PhiCut(aodTrack->Pt(), phi, aodTrack->Charge(), magf, cutLow, cutHigh))
1245       continue;
1246     
1247
1248     
1249     AliAODPid* aodPid = aodTrack->GetDetPid();
1250     Short_t ncl     = -10;
1251     Float_t dedx    = -10;
1252
1253     //TOF    
1254     Float_t beta = -99;
1255
1256
1257     if(aodPid) {
1258       ncl     = aodPid->GetTPCsignalN();
1259       dedx    = aodPid->GetTPCsignal();
1260       //TOF
1261       Bool_t IsTOFout=kFALSE;
1262       Float_t lengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF 
1263       Float_t timeTOF=-999;
1264
1265       if ((aodTrack->GetStatus()&AliESDtrack::kTOFout)==0)
1266         IsTOFout=kTRUE;
1267       
1268       lengthtrack=-999;//in aod we do not have that information, beta must be: beta=inttime/timeTOF 
1269       
1270       timeTOF=aodPid->GetTOFsignal();
1271       
1272       Double_t inttime[5]={0,0,0,0,0}; 
1273       aodPid->GetIntegratedTimes(inttime);// Returns the array with integrated times for each particle hypothesis
1274
1275       
1276       if ( !IsTOFout ){
1277         if ( ( lengthtrack != 0 ) && ( timeTOF != 0) )
1278           beta = inttime[0] / timeTOF;
1279       }
1280       
1281     }
1282
1283
1284     if(ncl<70)
1285       continue;
1286
1287
1288     Short_t pidCode     = 0; 
1289     
1290     if(fAnalysisMC) {
1291       
1292       const Int_t label = TMath::Abs(aodTrack->GetLabel());
1293       AliAODMCParticle* mcTrack = dynamic_cast<AliAODMCParticle*>(fMCArray->At(label));
1294       
1295       if (mcTrack){
1296         
1297         Int_t pdgCode = mcTrack->GetPdgCode();
1298         pidCode = GetPidCode(pdgCode);
1299         
1300       }
1301       
1302     }
1303
1304
1305
1306
1307     //if(momentum>0.6||momentum<0.4)//only p:0.4-0.6 GeV
1308         //continue;
1309     
1310     //etaLow
1311     //etaHigh
1312     if( momentum <= 0.6 && momentum >= 0.4  ){//only p:0.4-0.6 GeV, pion MIP
1313       if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){       
1314         if(momentum<0.6&&momentum>0.4){
1315           hMIPVsEta->Fill(eta,dedx);
1316           pMIPVsEta->Fill(eta,dedx);
1317         }
1318       }
1319       if( dedx > DeDxMIPMax+5 && dedx < 95 ){
1320         if(TMath::Abs(beta-1)<0.1){
1321           hPlateauVsEta->Fill(eta,dedx);
1322           pPlateauVsEta->Fill(eta,dedx); 
1323         }
1324       }
1325     }
1326     
1327     
1328     for(Int_t nh = 0; nh < 9; nh++) {
1329       
1330       if( eta > etaHigh[nh]/10.0 || eta < etaLow[nh]/10.0 )
1331         continue;
1332       
1333       
1334       if(fAnalysisMC){
1335         hMcOut[0][nh]->Fill(aodTrack->Pt());
1336         hMcOut[pidCode][nh]->Fill(aodTrack->Pt());
1337       }
1338       
1339       histAllCh[nh]->Fill(momentum, dedx);
1340       if(beta>1)
1341         histPiTof[nh]->Fill(momentum, dedx);
1342       
1343       if( momentum <= 0.6 && momentum >= 0.4  ){//only p:0.4-0.6 GeV, pion MIP
1344         //Fill  pion MIP, before calibration
1345         if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){     
1346           hMIPVsPhi[nh]->Fill(phi,dedx);
1347           pMIPVsPhi[nh]->Fill(phi,dedx);
1348           
1349           
1350           hMIPVsNch[nh]->Fill(multTPC,dedx);
1351           pMIPVsNch[nh]->Fill(multTPC,dedx);
1352           
1353         }
1354         
1355         //Fill electrons, before calibration
1356         if( dedx > DeDxMIPMax+5 && dedx < 95 ){
1357           if(TMath::Abs(beta-1)<0.1){
1358             hPlateauVsPhi[nh]->Fill(phi,dedx);
1359             pPlateauVsPhi[nh]->Fill(phi,dedx);
1360           }
1361         }
1362
1363       }
1364
1365     }//end loop over eta intervals
1366     
1367     
1368
1369     
1370     
1371   }//end of track loop
1372   
1373  
1374   
1375   
1376   
1377  
1378   
1379 }
1380 //__________________________________________________
1381 void AliAnalysisTaskQAHighPtDeDx::ProduceArrayV0ESD( AliESDEvent *ESDevent ){
1382   Int_t nv0s = ESDevent->GetNumberOfV0s();
1383   /*
1384   if(nv0s<1){
1385     return;
1386     }*/
1387   
1388   const AliESDVertex *myBestPrimaryVertex = ESDevent->GetPrimaryVertex();
1389   if (!myBestPrimaryVertex) return;
1390   if (!(myBestPrimaryVertex->GetStatus())) return;
1391   Double_t  lPrimaryVtxPosition[3];
1392   myBestPrimaryVertex->GetXYZ(lPrimaryVtxPosition);
1393   
1394   Double_t  lPrimaryVtxCov[6];
1395   myBestPrimaryVertex->GetCovMatrix(lPrimaryVtxCov);
1396   Double_t  lPrimaryVtxChi2 = myBestPrimaryVertex->GetChi2toNDF();
1397   
1398   AliAODVertex* myPrimaryVertex = new AliAODVertex(lPrimaryVtxPosition, lPrimaryVtxCov, lPrimaryVtxChi2, NULL, -1, AliAODVertex::kPrimary);
1399
1400
1401     
1402   //
1403   // LOOP OVER V0s, K0s, L, AL
1404   //
1405   
1406   
1407   for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
1408     
1409     // This is the begining of the V0 loop  
1410     AliESDv0 *esdV0 = ESDevent->GetV0(iV0);
1411     if (!esdV0) continue;
1412     
1413     //check onfly status
1414     if(esdV0->GetOnFlyStatus()!=0)
1415       continue;
1416
1417
1418     // AliESDTrack (V0 Daughters)
1419     UInt_t lKeyPos = (UInt_t)TMath::Abs(esdV0->GetPindex());
1420     UInt_t lKeyNeg = (UInt_t)TMath::Abs(esdV0->GetNindex());
1421     
1422     AliESDtrack *pTrack = ESDevent->GetTrack(lKeyPos);
1423     AliESDtrack *nTrack = ESDevent->GetTrack(lKeyNeg);
1424     if (!pTrack || !nTrack) {
1425       Printf("ERROR: Could not retreive one of the daughter track");
1426       continue;
1427     }
1428     
1429     // Remove like-sign
1430     if (pTrack->GetSign() == nTrack->GetSign()) {
1431       //cout<< "like sign, continue"<< endl;
1432       continue;
1433     } 
1434     
1435     // Eta cut on decay products
1436     if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
1437       continue;
1438     
1439     
1440     // Check if switch does anything!
1441     Bool_t isSwitched = kFALSE;
1442     if (pTrack->GetSign() < 0) { // switch
1443       
1444       isSwitched = kTRUE;
1445       AliESDtrack* helpTrack = nTrack;
1446       nTrack = pTrack;
1447       pTrack = helpTrack;
1448     }   
1449     
1450     //Double_t  lV0Position[3];
1451     //esdV0->GetXYZ(lV0Position[0], lV0Position[1], lV0Position[2]);
1452     
1453
1454     AliKFVertex primaryVtxKF( *myPrimaryVertex );
1455     AliKFParticle::SetField(ESDevent->GetMagneticField());
1456     
1457     // Also implement switch here!!!!!!
1458     AliKFParticle* negEKF  = 0; // e-
1459     AliKFParticle* posEKF  = 0; // e+
1460     AliKFParticle* negPiKF = 0; // pi -
1461     AliKFParticle* posPiKF = 0; // pi +
1462     AliKFParticle* posPKF  = 0; // p
1463     AliKFParticle* negAPKF = 0; // p-bar
1464     
1465     if(!isSwitched) {
1466       negEKF  = new AliKFParticle( *(esdV0->GetParamN()) , 11);
1467       posEKF  = new AliKFParticle( *(esdV0->GetParamP()) ,-11);
1468       negPiKF = new AliKFParticle( *(esdV0->GetParamN()) ,-211);
1469       posPiKF = new AliKFParticle( *(esdV0->GetParamP()) , 211);
1470       posPKF  = new AliKFParticle( *(esdV0->GetParamP()) , 2212);
1471       negAPKF = new AliKFParticle( *(esdV0->GetParamN()) ,-2212);
1472     } else { // switch + and - 
1473       negEKF  = new AliKFParticle( *(esdV0->GetParamP()) , 11);
1474       posEKF  = new AliKFParticle( *(esdV0->GetParamN()) ,-11);
1475       negPiKF = new AliKFParticle( *(esdV0->GetParamP()) ,-211);
1476       posPiKF = new AliKFParticle( *(esdV0->GetParamN()) , 211);
1477       posPKF  = new AliKFParticle( *(esdV0->GetParamN()) , 2212);
1478       negAPKF = new AliKFParticle( *(esdV0->GetParamP()) ,-2212);
1479     }
1480     
1481     AliKFParticle v0GKF;  // Gamma e.g. from pi0
1482     v0GKF+=(*negEKF);
1483     v0GKF+=(*posEKF);
1484     v0GKF.SetProductionVertex(primaryVtxKF);
1485     
1486     AliKFParticle v0K0sKF; // K0 short
1487     v0K0sKF+=(*negPiKF);
1488     v0K0sKF+=(*posPiKF);
1489     v0K0sKF.SetProductionVertex(primaryVtxKF);
1490     
1491     AliKFParticle v0LambdaKF; // Lambda
1492     v0LambdaKF+=(*negPiKF);
1493     v0LambdaKF+=(*posPKF);      
1494     v0LambdaKF.SetProductionVertex(primaryVtxKF);
1495     
1496     AliKFParticle v0AntiLambdaKF; // Lambda-bar
1497     v0AntiLambdaKF+=(*posPiKF);
1498     v0AntiLambdaKF+=(*negAPKF);
1499     v0AntiLambdaKF.SetProductionVertex(primaryVtxKF);
1500     
1501     Double_t dmassG     = v0GKF.GetMass();
1502     Double_t dmassK     = v0K0sKF.GetMass()-0.498;
1503     Double_t dmassL     = v0LambdaKF.GetMass()-1.116;
1504     Double_t dmassAL    = v0AntiLambdaKF.GetMass()-1.116;
1505
1506
1507     for( Int_t case_v0 = 0; case_v0 < 2; ++case_v0 ){
1508
1509
1510       switch(case_v0){
1511       case 0:{    
1512         
1513         Bool_t fillPos = kFALSE;
1514         Bool_t fillNeg = kFALSE;
1515         
1516         if(dmassG < 0.1)
1517           continue;
1518         
1519         if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01){
1520           continue;
1521         }
1522         
1523         if(dmassL<0.01){
1524           fillPos = kTRUE;
1525         }
1526         if(dmassAL<0.01) {
1527           if(fillPos)
1528             continue;
1529           fillNeg = kTRUE;
1530         }
1531         
1532         if(dmassK<0.01) {
1533           if(fillPos||fillNeg)
1534             continue;
1535           fillPos = kTRUE;
1536           fillNeg = kTRUE;
1537         }
1538         
1539         
1540         for(Int_t j = 0; j < 2; j++) {
1541           
1542           AliESDtrack* track = 0;
1543           
1544           if(j==0) {
1545             
1546             if(fillNeg)
1547               track = nTrack;
1548             else
1549               continue;
1550           } else {
1551             
1552             if(fillPos)
1553               track = pTrack;
1554             else
1555               continue;
1556           }
1557           
1558           if(track->GetTPCsignalN()<=70)continue;
1559           Double_t phi     = track->Phi();
1560           
1561           if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
1562             continue;
1563           
1564           
1565           //if(!PhiCut(pt, phi, charge, magf, cutLow, cutHigh, hPhi))
1566           //    continue;
1567           
1568           Double_t eta  = track->Eta();
1569           Double_t momentum = track->Pt();
1570           Double_t dedx = track->GetTPCsignal();
1571           
1572           if(fillPos&&fillNeg){
1573             
1574             
1575             if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){         
1576               if(momentum<0.6&&momentum>0.4){
1577                 hMIPVsEtaV0s->Fill(eta,dedx);
1578                 pMIPVsEtaV0s->Fill(eta,dedx);
1579               }
1580             }
1581             
1582             
1583           }
1584           
1585           for(Int_t nh = 0; nh < nHists; nh++) {
1586             
1587             
1588             
1589             if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
1590               continue;
1591             
1592             if(fillPos&&fillNeg){
1593               
1594               histPiV0[nh]->Fill(momentum, dedx);           
1595               
1596             }
1597             else{
1598               
1599               histPV0[nh]->Fill(momentum, dedx);
1600               
1601             }
1602             
1603           }
1604                   
1605         }//end loop over two tracks
1606
1607       };
1608         break;
1609
1610       case 1:{//gammas    
1611         
1612         Bool_t fillPos = kFALSE;
1613         Bool_t fillNeg = kFALSE;
1614         
1615
1616
1617         if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01) {
1618           if(dmassG<0.01 && dmassG>0.0001) {
1619             
1620
1621             if( TMath::Abs(nTrack->GetTPCsignal() - 85.0) < 5)
1622               fillPos = kTRUE;
1623             if( TMath::Abs(pTrack->GetTPCsignal() - 85.0) < 5)
1624               fillNeg = kTRUE;
1625             
1626           } else {
1627             continue;
1628           }
1629         }
1630
1631         //cout<<"fillPos="<<fillPos<<endl;
1632
1633         if(fillPos == kTRUE && fillNeg == kTRUE)      
1634           continue;
1635         
1636         
1637         AliESDtrack* track = 0;
1638         if(fillNeg)
1639           track = nTrack;
1640         else if(fillPos)
1641           track = pTrack;
1642         else
1643           continue;
1644
1645         Double_t dedx  = track->GetTPCsignal();
1646         Double_t eta  = track->Eta();
1647         Double_t phi  = track->Phi();
1648         Double_t momentum = track->P();
1649
1650         if(track->GetTPCsignalN()<=70)continue;
1651
1652         if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
1653           continue;
1654         
1655         for(Int_t nh = 0; nh < nHists; nh++) {
1656       
1657           if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
1658             continue;
1659
1660           histEV0[nh]->Fill(momentum, dedx);
1661
1662         }
1663         
1664       };
1665         break;
1666         
1667         
1668       }//end switch
1669       
1670     }//end loop over case V0
1671
1672     
1673     // clean up loop over v0
1674     
1675     delete negPiKF;
1676     delete posPiKF;
1677     delete posPKF;
1678     delete negAPKF;
1679     
1680     
1681     
1682   }
1683   
1684   
1685   delete myPrimaryVertex;
1686   
1687   
1688 }
1689 //__________________________________________________________________________
1690 void AliAnalysisTaskQAHighPtDeDx::ProduceArrayV0AOD( AliAODEvent *AODevent ){
1691   Int_t nv0s = AODevent->GetNumberOfV0s();
1692   /*
1693   if(nv0s<1){
1694     return;
1695     }*/
1696   
1697   AliAODVertex *myBestPrimaryVertex = AODevent->GetPrimaryVertex();
1698   if (!myBestPrimaryVertex) return;
1699   
1700
1701      
1702   // ################################
1703   // #### BEGINNING OF V0 CODE ######
1704   // ################################
1705   // This is the begining of the V0 loop  
1706   for (Int_t iV0 = 0; iV0 < nv0s; iV0++) {
1707     AliAODv0 *aodV0 = AODevent->GetV0(iV0);
1708     if (!aodV0) continue;
1709     
1710     
1711     //check onfly status
1712     if(aodV0->GetOnFlyStatus()!=0)
1713       continue;
1714
1715     // AliAODTrack (V0 Daughters)
1716     AliAODVertex* vertex = aodV0->GetSecondaryVtx();
1717     if (!vertex) {
1718       Printf("ERROR: Could not retrieve vertex");
1719       continue;
1720     }
1721     
1722     AliAODTrack *pTrack = (AliAODTrack*)vertex->GetDaughter(0);
1723     AliAODTrack *nTrack = (AliAODTrack*)vertex->GetDaughter(1);
1724     if (!pTrack || !nTrack) {
1725       Printf("ERROR: Could not retrieve one of the daughter track");
1726       continue;
1727     }
1728     
1729     // Remove like-sign
1730     if (pTrack->Charge() == nTrack->Charge()) {
1731       //cout<< "like sign, continue"<< endl;
1732       continue;
1733     } 
1734     
1735     // Make sure charge ordering is ok
1736     if (pTrack->Charge() < 0) {
1737       AliAODTrack* helpTrack = pTrack;
1738       pTrack = nTrack;
1739       nTrack = helpTrack;
1740     } 
1741     
1742     // Eta cut on decay products
1743     if(TMath::Abs(pTrack->Eta()) > fEtaCut || TMath::Abs(nTrack->Eta()) > fEtaCut)
1744       continue;
1745     
1746     
1747     Double_t dmassG  = aodV0->InvMass2Prongs(0,1,11,11);
1748     Double_t dmassK  = aodV0->MassK0Short()-0.498;
1749     Double_t dmassL  = aodV0->MassLambda()-1.116;
1750     Double_t dmassAL = aodV0->MassAntiLambda()-1.116;
1751     
1752     for( Int_t case_v0 = 0; case_v0 < 2; ++case_v0 ){
1753       
1754       
1755       switch(case_v0){
1756       case 0:{   
1757         Bool_t fillPos = kFALSE;
1758         Bool_t fillNeg = kFALSE;
1759         
1760         
1761         if(dmassG < 0.1)
1762           continue;
1763         
1764         if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01){
1765           continue;
1766         }
1767         
1768         if(dmassL<0.01){
1769           fillPos = kTRUE;
1770         }
1771         if(dmassAL<0.01) {
1772           if(fillPos)
1773             continue;
1774           fillNeg = kTRUE;
1775         }
1776         
1777         if(dmassK<0.01) {
1778           if(fillPos||fillNeg)
1779             continue;
1780           fillPos = kTRUE;
1781           fillNeg = kTRUE;
1782         }
1783         
1784         
1785         
1786         for(Int_t j = 0; j < 2; j++) {
1787           
1788           AliAODTrack* track = 0;
1789           
1790           if(j==0) {
1791             
1792             if(fillNeg)
1793               track = nTrack;
1794             else
1795               continue;
1796           } else {
1797             
1798             if(fillPos)
1799               track = pTrack;
1800             else
1801               continue;
1802           }
1803           
1804           if(track->GetTPCsignalN()<=70)continue;
1805           
1806           Double_t phi     = track->Phi();
1807           
1808           if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
1809             continue;
1810           
1811           
1812           //if(!PhiCut(pt, phi, charge, magf, cutLow, cutHigh, hPhi))
1813           //    continue;
1814           
1815           Double_t eta  = track->Eta();
1816           Double_t momentum = track->Pt();
1817           Double_t dedx = track->GetTPCsignal();
1818           
1819           if(fillPos&&fillNeg){
1820             
1821             
1822             if( dedx < DeDxMIPMax && dedx > DeDxMIPMin ){         
1823               if(momentum<0.6&&momentum>0.4){
1824                 hMIPVsEtaV0s->Fill(eta,dedx);
1825                 pMIPVsEtaV0s->Fill(eta,dedx);
1826               }
1827             }
1828             
1829             
1830           }
1831           
1832           for(Int_t nh = 0; nh < nHists; nh++) {
1833             
1834             
1835             
1836             if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
1837               continue;
1838             
1839             if(fillPos&&fillNeg){
1840               
1841               histPiV0[nh]->Fill(momentum, dedx);           
1842               
1843             }
1844             else{
1845               
1846               histPV0[nh]->Fill(momentum, dedx);
1847               
1848             }
1849             
1850           }
1851           
1852
1853         }//end loop over two tracks
1854       };
1855         break;
1856         
1857       case 1:{//gammas    
1858         
1859         Bool_t fillPos = kFALSE;
1860         Bool_t fillNeg = kFALSE;
1861         
1862         if(dmassK>0.01 && dmassL>0.01 && dmassAL>0.01) {
1863           if(dmassG<0.01 && dmassG>0.0001) {
1864             
1865             if( TMath::Abs(nTrack->GetTPCsignal() - 85.0) < 5)
1866               fillPos = kTRUE;
1867             if( TMath::Abs(pTrack->GetTPCsignal() - 85.0) < 5)
1868               fillNeg = kTRUE;
1869             
1870           } else {
1871             continue;
1872           }
1873         }
1874
1875         
1876         if(fillPos == kTRUE && fillNeg == kTRUE)      
1877           continue;
1878         
1879         
1880         AliAODTrack* track = 0;
1881         if(fillNeg)
1882           track = nTrack;
1883         else if(fillPos)
1884           track = pTrack;
1885         else
1886           continue;
1887
1888         Double_t dedx  = track->GetTPCsignal();
1889         Double_t eta  = track->Eta();
1890         Double_t phi  = track->Phi();
1891         Double_t momentum = track->P();
1892
1893         if(track->GetTPCsignalN()<=70)continue;
1894
1895         if(!PhiCut(track->Pt(), phi, track->Charge(), magf, cutLow, cutHigh))
1896           continue;
1897         
1898         for(Int_t nh = 0; nh < nHists; nh++) {
1899       
1900           if( eta < etaLow[nh]/10.0 || eta > etaHigh[nh]/10.0 )
1901             continue;
1902
1903           histEV0[nh]->Fill(momentum, dedx);
1904
1905         }
1906         
1907       };
1908         break;
1909
1910
1911       }//end switch
1912     }//end loop over V0s cases
1913     
1914   }//end loop over v0's
1915   
1916   
1917   
1918   
1919 }
1920 Bool_t AliAnalysisTaskQAHighPtDeDx::PhiCut(Double_t pt, Double_t phi, Double_t q, Float_t   mag, TF1* phiCutLow, TF1* phiCutHigh)
1921 {
1922   if(pt < 2.0)
1923     return kTRUE;
1924   
1925   //Double_t phi = track->Phi();
1926   if(mag < 0)    // for negatve polarity field
1927     phi = TMath::TwoPi() - phi;
1928   if(q < 0) // for negatve charge
1929     phi = TMath::TwoPi()-phi;
1930   
1931   phi += TMath::Pi()/18.0; // to center gap in the middle
1932   phi = fmod(phi, TMath::Pi()/9.0);
1933     
1934   if(phi<phiCutHigh->Eval(pt) 
1935      && phi>phiCutLow->Eval(pt))
1936   return kFALSE; // reject track
1937
1938   hPhi->Fill(pt, phi);
1939
1940   return kTRUE;
1941 }