fix for coverity (O. Busch)
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetChem.cxx
1 /*************************************************************************
2  *                                                                       *
3  *                                                                       *
4  *      Task for Jet Chemistry Analysis in PWG4 Jet Task Force Train     *
5  *                                                                       *
6  *                                                                       *
7  *************************************************************************/
8
9 /**************************************************************************
10  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
11  *                                                                        *
12  * Author: The ALICE Off-line Project.                                    *
13  * Contributors are mentioned in the code where appropriate.              *
14  *                                                                        *
15  * Permission to use, copy, modify and distribute this software and its   *
16  * documentation strictly for non-commercial purposes is hereby granted   *
17  * without fee, provided that the above copyright notice appears in all   *
18  * copies and that both the copyright notice and this permission notice   *
19  * appear in the supporting documentation. The authors make no claims     *
20  * about the suitability of this software for any purpose. It is          *
21  * provided "as is" without express or implied warranty.                  *
22  **************************************************************************/
23
24 /* $Id: */
25
26 #include "TH2F.h"
27 #include "TH3F.h"
28 #include "TProfile.h"
29 #include "THnSparse.h"
30
31 #include "AliAnalysisHelperJetTasks.h"
32 #include "AliAnalysisManager.h"
33 #include "AliAODHandler.h" 
34 #include "AliAODInputHandler.h" 
35 #include "AliESDEvent.h"
36 #include "AliGenPythiaEventHeader.h"
37 #include "AliGenHijingEventHeader.h"
38
39 #include "AliAODEvent.h"
40 #include "AliAODJet.h"
41 #include "AliAODv0.h"
42 #include "AliAODTrack.h"
43
44
45 #include "AliPID.h" 
46 #include "AliExternalTrackParam.h"
47
48 #include "AliAnalysisTaskJetChem.h"
49
50 ClassImp(AliAnalysisTaskJetChem)
51
52 //____________________________________________________________________________
53 AliAnalysisTaskJetChem::AliAnalysisTaskJetChem()
54    : AliAnalysisTaskFragmentationFunction()
55    ,fK0Type(0)
56    ,fFilterMaskK0(0)
57    ,fListK0s(0)
58    ,fV0QAK0(0)
59    ,fFFHistosRecCutsK0Evt(0)      
60    ,fFFHistosIMK0AllEvt(0)        
61    ,fFFHistosIMK0Jet(0)           
62    ,fFFHistosIMK0Cone(0)
63    ,fFFHistosPhiCorrIMK0(0)
64    ,fFFIMNBinsJetPt(0)    
65    ,fFFIMJetPtMin(0) 
66    ,fFFIMJetPtMax(0)
67    ,fFFIMNBinsInvM(0) 
68    ,fFFIMInvMMin(0)   
69    ,fFFIMInvMMax(0)   
70    ,fFFIMNBinsPt(0)      
71    ,fFFIMPtMin(0)        
72    ,fFFIMPtMax(0)        
73    ,fFFIMNBinsXi(0)      
74    ,fFFIMXiMin(0)        
75    ,fFFIMXiMax(0)        
76    ,fFFIMNBinsZ(0)       
77    ,fFFIMZMin(0)         
78    ,fFFIMZMax(0)
79    ,fPhiCorrIMNBinsPt(0)
80    ,fPhiCorrIMPtMin(0)
81    ,fPhiCorrIMPtMax(0)
82    ,fPhiCorrIMNBinsPhi(0)
83    ,fPhiCorrIMPhiMin(0)
84    ,fPhiCorrIMPhiMax(0)
85    ,fPhiCorrIMNBinsInvM(0)
86    ,fPhiCorrIMInvMMin(0)
87    ,fPhiCorrIMInvMMax(0)
88    ,fh1K0Mult(0)
89    ,fh1dPhiJetK0(0)
90 {
91    // default constructor
92 }
93
94 //__________________________________________________________________________________________
95 AliAnalysisTaskJetChem::AliAnalysisTaskJetChem(const char *name) 
96   : AliAnalysisTaskFragmentationFunction(name)
97   ,fK0Type(0)  
98   ,fFilterMaskK0(0)
99   ,fListK0s(0)
100   ,fV0QAK0(0)
101   ,fFFHistosRecCutsK0Evt(0)      
102   ,fFFHistosIMK0AllEvt(0)        
103   ,fFFHistosIMK0Jet(0)           
104   ,fFFHistosIMK0Cone(0)
105   ,fFFHistosPhiCorrIMK0(0)
106   ,fFFIMNBinsJetPt(0)    
107   ,fFFIMJetPtMin(0) 
108   ,fFFIMJetPtMax(0)
109   ,fFFIMNBinsInvM(0) 
110   ,fFFIMInvMMin(0)   
111   ,fFFIMInvMMax(0)   
112   ,fFFIMNBinsPt(0)      
113   ,fFFIMPtMin(0)        
114   ,fFFIMPtMax(0)        
115   ,fFFIMNBinsXi(0)      
116   ,fFFIMXiMin(0)        
117   ,fFFIMXiMax(0)        
118   ,fFFIMNBinsZ(0)       
119   ,fFFIMZMin(0)         
120   ,fFFIMZMax(0)
121   ,fPhiCorrIMNBinsPt(0)
122   ,fPhiCorrIMPtMin(0)
123   ,fPhiCorrIMPtMax(0)
124   ,fPhiCorrIMNBinsPhi(0)
125   ,fPhiCorrIMPhiMin(0)
126   ,fPhiCorrIMPhiMax(0)
127   ,fPhiCorrIMNBinsInvM(0)
128   ,fPhiCorrIMInvMMin(0)
129   ,fPhiCorrIMInvMMax(0)
130   ,fh1K0Mult(0)
131   ,fh1dPhiJetK0(0)     
132 {
133   // constructor
134   
135   DefineOutput(1,TList::Class());  
136 }
137
138 //__________________________________________________________________________________________________________________________
139 AliAnalysisTaskJetChem::AliAnalysisTaskJetChem(const  AliAnalysisTaskJetChem &copy)
140   : AliAnalysisTaskFragmentationFunction()
141   ,fK0Type(copy.fK0Type)
142   ,fFilterMaskK0(copy.fFilterMaskK0)
143   ,fListK0s(copy.fListK0s)
144   ,fV0QAK0(copy.fV0QAK0)
145   ,fFFHistosRecCutsK0Evt(copy.fFFHistosRecCutsK0Evt)      
146   ,fFFHistosIMK0AllEvt(copy.fFFHistosIMK0AllEvt)        
147   ,fFFHistosIMK0Jet(copy.fFFHistosIMK0Jet)           
148   ,fFFHistosIMK0Cone(copy.fFFHistosIMK0Cone)          
149   ,fFFHistosPhiCorrIMK0(copy.fFFHistosPhiCorrIMK0)          
150   ,fFFIMNBinsJetPt(copy.fFFIMNBinsJetPt)   
151   ,fFFIMJetPtMin(copy.fFFIMJetPtMin)     
152   ,fFFIMJetPtMax(copy.fFFIMJetPtMax)     
153   ,fFFIMNBinsInvM(copy.fFFIMNBinsInvM)  
154   ,fFFIMInvMMin(copy.fFFIMInvMMin)    
155   ,fFFIMInvMMax(copy.fFFIMInvMMax)    
156   ,fFFIMNBinsPt(copy.fFFIMNBinsPt)      
157   ,fFFIMPtMin(copy.fFFIMPtMin)        
158   ,fFFIMPtMax(copy.fFFIMPtMax)        
159   ,fFFIMNBinsXi(copy.fFFIMNBinsXi)      
160   ,fFFIMXiMin(copy.fFFIMXiMin)        
161   ,fFFIMXiMax(copy.fFFIMXiMax)        
162   ,fFFIMNBinsZ(copy.fFFIMNBinsZ)       
163   ,fFFIMZMin(copy.fFFIMZMin)         
164   ,fFFIMZMax(copy.fFFIMZMax) 
165   ,fPhiCorrIMNBinsPt(copy.fPhiCorrIMNBinsPt)
166   ,fPhiCorrIMPtMin(copy.fPhiCorrIMPtMin)
167   ,fPhiCorrIMPtMax(copy.fPhiCorrIMPtMax)
168   ,fPhiCorrIMNBinsPhi(copy.fPhiCorrIMNBinsPhi)
169   ,fPhiCorrIMPhiMin(copy.fPhiCorrIMPhiMin)
170   ,fPhiCorrIMPhiMax(copy.fPhiCorrIMPhiMax)
171   ,fPhiCorrIMNBinsInvM(copy.fPhiCorrIMNBinsInvM)
172   ,fPhiCorrIMInvMMin(copy.fPhiCorrIMInvMMin)
173   ,fPhiCorrIMInvMMax(copy.fPhiCorrIMInvMMax)
174   ,fh1K0Mult(copy.fh1K0Mult)
175   ,fh1dPhiJetK0(copy.fh1dPhiJetK0)
176 {
177   // copy constructor
178   
179 }
180
181 // _________________________________________________________________________________________________________________________________
182 AliAnalysisTaskJetChem& AliAnalysisTaskJetChem::operator=(const AliAnalysisTaskJetChem& o)
183 {
184   // assignment
185   
186   if(this!=&o){
187     AliAnalysisTaskFragmentationFunction::operator=(o);
188
189     fK0Type                         = o.fK0Type;
190     fFilterMaskK0                   = o.fFilterMaskK0;
191     fListK0s                        = o.fListK0s;
192     fV0QAK0                         = o.fV0QAK0;
193     fFFHistosRecCutsK0Evt           = o.fFFHistosRecCutsK0Evt;      
194     fFFHistosIMK0AllEvt             = o.fFFHistosIMK0AllEvt;        
195     fFFHistosIMK0Jet                = o.fFFHistosIMK0Jet;           
196     fFFHistosIMK0Cone               = o.fFFHistosIMK0Cone;          
197     fFFHistosPhiCorrIMK0            = o.fFFHistosPhiCorrIMK0;
198     fFFIMNBinsJetPt                 = o.fFFIMNBinsJetPt;    
199     fFFIMJetPtMin                   = o.fFFIMJetPtMin; 
200     fFFIMJetPtMax                   = o.fFFIMJetPtMax;
201     fFFIMNBinsPt                    = o.fFFIMNBinsPt;      
202     fFFIMPtMin                      = o.fFFIMPtMin;        
203     fFFIMPtMax                      = o.fFFIMPtMax;        
204     fFFIMNBinsXi                    = o.fFFIMNBinsXi;      
205     fFFIMXiMin                      = o.fFFIMXiMin;        
206     fFFIMXiMax                      = o.fFFIMXiMax;        
207     fFFIMNBinsZ                     = o.fFFIMNBinsZ;       
208     fFFIMZMin                       = o.fFFIMZMin;         
209     fFFIMZMax                       = o.fFFIMZMax;         
210     fPhiCorrIMNBinsPt               = o.fPhiCorrIMNBinsPt;
211     fPhiCorrIMPtMin                 = o.fPhiCorrIMPtMin;
212     fPhiCorrIMPtMax                 = o.fPhiCorrIMPtMax;
213     fPhiCorrIMNBinsPhi              = o.fPhiCorrIMNBinsPhi;
214     fPhiCorrIMPhiMin                = o.fPhiCorrIMPhiMin;
215     fPhiCorrIMPhiMax                = o.fPhiCorrIMPhiMax;
216     fPhiCorrIMNBinsInvM             = o.fPhiCorrIMNBinsInvM;
217     fPhiCorrIMInvMMin               = o.fPhiCorrIMInvMMin;
218     fPhiCorrIMInvMMax               = o.fPhiCorrIMInvMMax;
219     fh1K0Mult                       = o.fh1K0Mult;
220     fh1dPhiJetK0                    = o.fh1dPhiJetK0;
221   }
222     
223   return *this;
224 }
225
226 //_______________________________________________
227 AliAnalysisTaskJetChem::~AliAnalysisTaskJetChem()
228 {
229   // destructor  
230
231   
232   if(fTracksRecCuts) delete fTracksRecCuts;
233   if(fJetsRecCuts) delete fJetsRecCuts;
234   if(fListK0s) delete fListK0s;
235   if(fCommonHistList) delete fCommonHistList;
236
237 }
238
239 //________________________________________________________________________________________________________________________________
240 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AliFragFuncHistosInvMass(const char* name, 
241                                                                            Int_t nJetPt, Float_t jetPtMin, Float_t jetPtMax,  
242                                                                            Int_t nInvMass, Float_t invMassMin, Float_t invMassMax,
243                                                                            Int_t nPt, Float_t ptMin, Float_t ptMax,
244                                                                            Int_t nXi, Float_t xiMin, Float_t xiMax,
245                                                                            Int_t nZ , Float_t zMin , Float_t zMax )
246   : TObject()
247   ,fNBinsJetPt(nJetPt)
248   ,fJetPtMin(jetPtMin)
249   ,fJetPtMax(jetPtMax)
250   ,fNBinsInvMass(nInvMass)
251   ,fInvMassMin(invMassMin)  
252   ,fInvMassMax(invMassMax)
253   ,fNBinsPt(nPt) 
254   ,fPtMin(ptMin)   
255   ,fPtMax(ptMax)   
256   ,fNBinsXi(nXi) 
257   ,fXiMin(xiMin)   
258   ,fXiMax(xiMax)   
259   ,fNBinsZ(nZ)  
260   ,fZMin(zMin)    
261   ,fZMax(zMax)    
262   ,fh3TrackPt(0)
263   ,fh3Xi(0)
264   ,fh3Z(0)
265   ,fh1JetPt(0)
266   ,fNameFF(name)
267 {
268   // default constructor
269
270 }
271
272 //______________________________________________________________________________________________________________
273 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AliFragFuncHistosInvMass(const AliFragFuncHistosInvMass& copy)
274   : TObject()
275   ,fNBinsJetPt(copy.fNBinsJetPt)
276   ,fJetPtMin(copy.fJetPtMin)
277   ,fJetPtMax(copy.fJetPtMax)
278   ,fNBinsInvMass(copy.fNBinsInvMass)
279   ,fInvMassMin(copy.fInvMassMin)  
280   ,fInvMassMax(copy.fInvMassMax)
281   ,fNBinsPt(copy.fNBinsPt) 
282   ,fPtMin(copy.fPtMin)   
283   ,fPtMax(copy.fPtMax)   
284   ,fNBinsXi(copy.fNBinsXi) 
285   ,fXiMin(copy.fXiMin)   
286   ,fXiMax(copy.fXiMax)   
287   ,fNBinsZ(copy.fNBinsZ)  
288   ,fZMin(copy.fZMin)    
289   ,fZMax(copy.fZMax)    
290   ,fh3TrackPt(copy.fh3TrackPt)
291   ,fh3Xi(copy.fh3Xi)
292   ,fh3Z(copy.fh3Z)
293   ,fh1JetPt(copy.fh1JetPt)
294   ,fNameFF(copy.fNameFF)
295 {
296   // copy constructor
297 }
298
299 //______________________________________________________________________________________________________________________________________________________________________
300 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass& AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::operator=(const AliAnalysisTaskJetChem::AliFragFuncHistosInvMass& o)
301 {
302   // assignment
303   
304   if(this!=&o){
305     TObject::operator=(o);
306     fNBinsJetPt   = o.fNBinsJetPt;
307     fJetPtMin     = o.fJetPtMin;
308     fJetPtMax     = o.fJetPtMax;
309     fNBinsInvMass = o.fNBinsInvMass;
310     fInvMassMin   = o.fInvMassMin;  
311     fInvMassMax   = o.fInvMassMax;
312     fNBinsPt      = o.fNBinsPt; 
313     fPtMin        = o.fPtMin;   
314     fPtMax        = o.fPtMax;   
315     fNBinsXi      = o.fNBinsXi; 
316     fXiMin        = o.fXiMin;   
317     fXiMax        = o.fXiMax;   
318     fNBinsZ       = o.fNBinsZ;  
319     fZMin         = o.fZMin;    
320     fZMax         = o.fZMax;    
321     fh3TrackPt    = o.fh3TrackPt;
322     fh3Xi         = o.fh3Xi;
323     fh3Z          = o.fh3Z;
324     fh1JetPt      = o.fh1JetPt;
325     fNameFF       = o.fNameFF;
326   }
327     
328   return *this;
329 }
330
331 //___________________________________________________________________________
332 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::~AliFragFuncHistosInvMass()
333
334   // destructor 
335
336   if(fh1JetPt)   delete fh1JetPt;
337   if(fh3TrackPt) delete fh3TrackPt;
338   if(fh3Xi)      delete fh3Xi;
339   if(fh3Z)       delete fh3Z;
340 }
341
342 //_________________________________________________________________
343 void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::DefineHistos()
344 {
345   // book FF histos
346
347   fh1JetPt   = new TH1F(Form("fh1FFJetPtIM%s", fNameFF.Data()),"",fNBinsJetPt,fJetPtMin,fJetPtMax);
348   fh3TrackPt = new TH3F(Form("fh3FFTrackPtIM%s",fNameFF.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsPt, fPtMin, fPtMax);
349   fh3Xi      = new TH3F(Form("fh3FFXiIM%s", fNameFF.Data()),"", fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsXi, fXiMin, fXiMax);
350   fh3Z       = new TH3F(Form("fh3FFZIM%s", fNameFF.Data()),"", fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsZ, fZMin, fZMax);
351
352   AliAnalysisTaskFragmentationFunction::SetProperties(fh1JetPt, "p_{t} (GeV/c)", "entries"); 
353   AliAnalysisTaskJetChem::SetProperties(fh3TrackPt,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","p_{t} (GeV/c)");
354   AliAnalysisTaskJetChem::SetProperties(fh3Xi,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","#xi");
355   AliAnalysisTaskJetChem::SetProperties(fh3Z,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","z");
356 }
357
358 //________________________________________________________________________________________________________________________________
359 void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::FillFF(Float_t trackPt, Float_t invM, Float_t jetPt, Bool_t incrementJetPt)
360 {
361   // fill FF
362  
363   if(incrementJetPt) fh1JetPt->Fill(jetPt);    
364   fh3TrackPt->Fill(jetPt,invM,trackPt);
365   
366   Double_t z = 0.;
367   if(jetPt>0) z = trackPt / jetPt;
368   Double_t xi = 0;
369   if(z>0) xi = TMath::Log(1/z);
370   
371   fh3Xi->Fill(jetPt,invM,xi);
372   fh3Z->Fill(jetPt,invM,z);
373 }
374
375 //___________________________________________________________________________________
376 void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AddToOutput(TList* list) const
377 {
378   // add histos to list
379
380   list->Add(fh1JetPt);
381   
382   list->Add(fh3TrackPt);
383   list->Add(fh3Xi);
384   list->Add(fh3Z);
385 }
386
387 // ---
388
389
390 //_______________________________________________________________________________________________________
391 AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::AliFragFuncHistosPhiCorrInvMass(const char* name,
392                                                                                          Int_t nPt,   Float_t ptMin,   Float_t ptMax,
393                                                                                          Int_t nPhi,  Float_t phiMin,  Float_t phiMax,
394                                                                                          Int_t nInvMass,  Float_t invMassMin,  Float_t invMassMax)
395   : TObject()
396   ,fNBinsPt(nPt)
397   ,fPtMin(ptMin)
398   ,fPtMax(ptMax)
399   ,fNBinsPhi(nPhi)
400   ,fPhiMin(phiMin)
401   ,fPhiMax(phiMax)
402   ,fNBinsInvMass(nInvMass)
403   ,fInvMassMin(invMassMin)  
404   ,fInvMassMax(invMassMax)
405   ,fh3PhiCorr(0) 
406   ,fNamePhiCorr(name) 
407 {
408   // default constructor
409 }
410
411 //____________________________________________________________________________________________________________________________________
412 AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::AliFragFuncHistosPhiCorrInvMass(const AliFragFuncHistosPhiCorrInvMass& copy)
413   : TObject()
414   ,fNBinsPt(copy.fNBinsPt)
415   ,fPtMin(copy.fPtMin)
416   ,fPtMax(copy.fPtMax)
417   ,fNBinsPhi(copy.fNBinsPhi)
418   ,fPhiMin(copy.fPhiMin)
419   ,fPhiMax(copy.fPhiMax)
420   ,fNBinsInvMass(copy.fNBinsInvMass)
421   ,fInvMassMin(copy.fInvMassMin)  
422   ,fInvMassMax(copy.fInvMassMax)
423   ,fh3PhiCorr(copy.fh3PhiCorr)
424   ,fNamePhiCorr(copy.fNamePhiCorr)
425 {
426   // copy constructor
427 }
428
429 //________________________________________________________________________________________________________________________________________________________________________
430 AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass& AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::operator=(const AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass& o)
431 {
432   // assignment
433   
434   if(this!=&o){
435     TObject::operator=(o);
436     fNBinsPt      = o.fNBinsPt;
437     fPtMin        = o.fPtMin;
438     fPtMax        = o.fPtMax;
439     fNBinsPhi     = o.fNBinsPhi;
440     fPhiMin       = o.fPhiMin;
441     fPhiMax       = o.fPhiMax;
442     fNBinsInvMass = o.fNBinsInvMass;
443     fInvMassMin   = o.fInvMassMin;  
444     fInvMassMax   = o.fInvMassMax;
445     
446     fh3PhiCorr    = o.fh3PhiCorr;
447     fNamePhiCorr  = o.fNamePhiCorr;
448   }
449   
450   return *this;
451 }
452
453 //_________________________________________________________________________________________
454 AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::~AliFragFuncHistosPhiCorrInvMass()
455 {
456   // destructor 
457   
458   if(fh3PhiCorr) delete fh3PhiCorr;
459 }
460
461 //__________________________________________________________________________
462 void AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::DefineHistos()
463 {
464   // book jet QA histos
465
466   fh3PhiCorr  = new TH3F(Form("fh3PhiCorrIM%s", fNamePhiCorr.Data()), 
467                          Form("%s: p_{t} - #phi - m_{inv} distribution",fNamePhiCorr.Data()), 
468                          fNBinsPt, fPtMin, fPtMax, 
469                          fNBinsPhi, fPhiMin, fPhiMax,
470                          fNBinsInvMass, fInvMassMin, fInvMassMax);
471   
472   AliAnalysisTaskJetChem::SetProperties(fh3PhiCorr, "p_{t} (GeV/c)", "#phi", "m_{inv} (GeV/c^2)"); 
473 }
474
475 //___________________________________________________________________________________________________________
476 void AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::FillPhiCorr(Float_t pt, Float_t phi, Float_t invM)
477 {
478   // fill jet QA histos 
479
480   fh3PhiCorr->Fill(pt, phi, invM);
481 }
482
483 //______________________________________________________________________________________________
484 void AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::AddToOutput(TList* list) const 
485 {
486   // add histos to list
487
488   list->Add(fh3PhiCorr);
489 }
490
491 //____________________________________________________
492 void AliAnalysisTaskJetChem::UserCreateOutputObjects()
493 {
494   // create output objects
495
496   if(fDebug > 1) Printf("AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()");
497  
498   // create list of tracks and jets 
499   
500   fTracksRecCuts = new TList();
501   fTracksRecCuts->SetOwner(kFALSE);  
502
503   fJetsRecCuts = new TList();
504   fJetsRecCuts->SetOwner(kFALSE);
505
506   fListK0s = new TList(); 
507   fListK0s->SetOwner(kFALSE);
508
509   //
510   // Create histograms / output container
511   //
512
513   OpenFile(1);
514   fCommonHistList = new TList();
515   
516   Bool_t oldStatus = TH1::AddDirectoryStatus();
517   TH1::AddDirectory(kFALSE);
518         
519
520   // histograms inherited from AliAnalysisTaskFragmentationFunction
521
522   fh1EvtSelection            = new TH1F("fh1EvtSelection", "Event Selection", 6, -0.5, 5.5);
523   fh1EvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
524   fh1EvtSelection->GetXaxis()->SetBinLabel(2,"event selection: rejected");
525   fh1EvtSelection->GetXaxis()->SetBinLabel(3,"event class: rejected");
526   fh1EvtSelection->GetXaxis()->SetBinLabel(4,"vertex Ncontr: rejected");
527   fh1EvtSelection->GetXaxis()->SetBinLabel(5,"vertex z: rejected");
528   fh1EvtSelection->GetXaxis()->SetBinLabel(6,"vertex type: rejected");
529
530   fh1EvtCent                 = new TH1F("fh1EvtCent","centrality",100,0.,100.);
531   fh1VertexNContributors     = new TH1F("fh1VertexNContributors", "Vertex N contributors", 11,-.5, 10.5);
532   fh1VertexZ                 = new TH1F("fh1VertexZ", "Vertex z distribution", 30, -15., 15.);
533   fh1Xsec                    = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
534   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
535   fh1Trials                  = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
536   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
537   fh1PtHard                  = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
538   fh1PtHardTrials            = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
539
540   fh1nRecJetsCuts            = new TH1F("fh1nRecJetsCuts","reconstructed jets per event",10,-0.5,9.5);
541
542
543   // histograms jetChem proper 
544
545   fh1EvtMult                 = new TH1F("fh1EvtMult","multiplicity",1200,0.,12000.);
546   fh1K0Mult                  = new TH1F("fh1K0Mult","K0 multiplicity",500,0.,500.);
547   fh1dPhiJetK0               = new TH1F("fh1dPhiJetK0","",640,-1,5.4);
548
549
550
551   fFFHistosRecCuts           = new AliFragFuncHistos("RecCuts", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
552                                                      fFFNBinsPt, fFFPtMin, fFFPtMax, 
553                                                      fFFNBinsXi, fFFXiMin, fFFXiMax,  
554                                                      fFFNBinsZ , fFFZMin , fFFZMax);
555   
556   fV0QAK0                   = new AliFragFuncQATrackHistos("V0QAK0",fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax, 
557                                                             fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
558                                                             fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax, 
559                                                             fQATrackHighPtThreshold);
560   
561   fFFHistosRecCutsK0Evt      = new AliFragFuncHistos("RecCutsK0Evt", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax, 
562                                                      fFFNBinsPt, fFFPtMin, fFFPtMax, 
563                                                      fFFNBinsXi, fFFXiMin, fFFXiMax,  
564                                                      fFFNBinsZ , fFFZMin , fFFZMax);
565   
566   
567   fFFHistosIMK0AllEvt        = new AliFragFuncHistosInvMass("K0AllEvt", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax, 
568                                                             fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
569                                                             fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax, 
570                                                             fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,  
571                                                             fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
572   
573   fFFHistosIMK0Jet           = new AliFragFuncHistosInvMass("K0Jet", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax, 
574                                                             fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
575                                                             fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax, 
576                                                             fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,  
577                                                             fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
578   
579   
580   fFFHistosIMK0Cone          = new AliFragFuncHistosInvMass("K0Cone", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax, 
581                                                             fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
582                                                             fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax, 
583                                                             fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,  
584                                                             fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
585   
586   fFFHistosPhiCorrIMK0       = new AliFragFuncHistosPhiCorrInvMass("K0",fPhiCorrIMNBinsPt, fPhiCorrIMPtMin, fPhiCorrIMPtMax, 
587                                                                    fPhiCorrIMNBinsPhi, fPhiCorrIMPhiMin, fPhiCorrIMPhiMax,  
588                                                                    fPhiCorrIMNBinsInvM , fPhiCorrIMInvMMin , fPhiCorrIMInvMMax);
589   
590   
591   fV0QAK0->DefineHistos();
592   fFFHistosRecCuts->DefineHistos();
593   fFFHistosRecCutsK0Evt->DefineHistos();
594   fFFHistosIMK0AllEvt->DefineHistos();
595   fFFHistosIMK0Jet->DefineHistos();
596   fFFHistosIMK0Cone->DefineHistos();
597   fFFHistosPhiCorrIMK0->DefineHistos();
598
599   const Int_t saveLevel = 5;
600   if(saveLevel>0){
601     
602     fCommonHistList->Add(fh1EvtSelection);
603     fCommonHistList->Add(fh1EvtCent);
604     fCommonHistList->Add(fh1VertexNContributors);
605     fCommonHistList->Add(fh1VertexZ);
606     fCommonHistList->Add(fh1Xsec);
607     fCommonHistList->Add(fh1Trials);
608     fCommonHistList->Add(fh1PtHard);
609     fCommonHistList->Add(fh1PtHardTrials);
610     fCommonHistList->Add(fh1nRecJetsCuts);
611     fCommonHistList->Add(fh1EvtMult);
612     fCommonHistList->Add(fh1K0Mult);
613     fCommonHistList->Add(fh1dPhiJetK0);
614
615     fV0QAK0->AddToOutput(fCommonHistList);
616     fFFHistosRecCuts->AddToOutput(fCommonHistList);
617     fFFHistosRecCutsK0Evt->AddToOutput(fCommonHistList);
618
619     fFFHistosIMK0AllEvt->AddToOutput(fCommonHistList);
620     fFFHistosIMK0Jet->AddToOutput(fCommonHistList);
621     fFFHistosIMK0Cone->AddToOutput(fCommonHistList);
622     fFFHistosPhiCorrIMK0->AddToOutput(fCommonHistList);
623   }
624   
625   // =========== Switch on Sumw2 for all histos ===========
626   for (Int_t i=0; i<fCommonHistList->GetEntries(); ++i){
627     TH1 *h1 = dynamic_cast<TH1*>(fCommonHistList->At(i));
628     if (h1) h1->Sumw2();
629     else{
630       THnSparse *hnSparse = dynamic_cast<THnSparse*>(fCommonHistList->At(i));
631       if(hnSparse) hnSparse->Sumw2();
632     }
633   }
634   
635   TH1::AddDirectory(oldStatus);
636 }
637
638 //_______________________________________________
639 void AliAnalysisTaskJetChem::UserExec(Option_t *) 
640 {
641   // Main loop
642   // Called for each event
643
644   if(fDebug > 1) Printf("Analysis event #%5d", (Int_t) fEntry);
645
646   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
647     ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
648   if(inputHandler->IsEventSelected() & AliVEvent::kMB){
649     if(fDebug > 1)  Printf(" Trigger Selection: event ACCEPTED ... ");
650     fh1EvtSelection->Fill(1.);
651   } else {
652     fh1EvtSelection->Fill(0.);
653     if(inputHandler->InheritsFrom("AliESDInputHandler") && fUsePhysicsSelection){ // PhysicsSelection only with ESD input
654       if (fDebug > 1 ) Printf(" Trigger Selection: event REJECTED ... ");
655       PostData(1, fCommonHistList);
656       return;
657     }
658   }
659
660   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
661   if(!fESD){
662     if(fDebug>3) Printf("%s:%d ESDEvent not found in the input", (char*)__FILE__,__LINE__);
663   }
664   
665   fMCEvent = MCEvent();
666   if(!fMCEvent){
667     if(fDebug>3) Printf("%s:%d MCEvent not found in the input", (char*)__FILE__,__LINE__);
668   }
669   
670   // get AOD event from input/ouput
671   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
672   if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
673     fAOD  =  ((AliAODInputHandler*)handler)->GetEvent();
674     if(fUseAODInputJets) fAODJets = fAOD;
675     if (fDebug > 1)  Printf("%s:%d AOD event from input", (char*)__FILE__,__LINE__);
676   }
677   else {
678     handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
679     if( handler && handler->InheritsFrom("AliAODHandler") ) {
680       fAOD = ((AliAODHandler*)handler)->GetAOD();
681       fAODJets = fAOD;
682       if (fDebug > 1)  Printf("%s:%d AOD event from output", (char*)__FILE__,__LINE__);
683     }
684   }
685   
686   if(!fAODJets && !fUseAODInputJets){ // case we have AOD in input & output and want jets from output
687     TObject* outHandler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
688     if( outHandler && outHandler->InheritsFrom("AliAODHandler") ) {
689       fAODJets = ((AliAODHandler*)outHandler)->GetAOD();
690       if (fDebug > 1)  Printf("%s:%d jets from output AOD", (char*)__FILE__,__LINE__);
691     }
692   }
693   
694   if(fNonStdFile.Length()!=0){
695     // case we have an AOD extension - fetch the jets from the extended output
696     
697     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
698     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);    
699     if(!fAODExtension){
700       if(fDebug>1)Printf("AODExtension not found for %s",fNonStdFile.Data());
701     }
702   }
703   
704   if(!fAOD){
705     Printf("%s:%d AODEvent not found", (char*)__FILE__,__LINE__);
706     return;
707   }
708   if(!fAODJets){
709     Printf("%s:%d AODEvent with jet branch not found", (char*)__FILE__,__LINE__);
710     return;
711   }
712   
713
714   // event selection  *****************************************
715   
716   Double_t centPercent = -1;
717   if(fEventClass>0){
718     Int_t cl = 0;
719     if(handler && handler->InheritsFrom("AliAODInputHandler")){ 
720       // since it is not supported by the helper task define own classes
721       centPercent = fAOD->GetHeader()->GetCentrality();
722       cl = 1;
723       if(centPercent>10) cl = 2;
724       if(centPercent>30) cl = 3;
725       if(centPercent>50) cl = 4;
726     }
727     else {
728       cl = AliAnalysisHelperJetTasks::EventClass();
729       if(fESD) centPercent = fESD->GetCentrality()->GetCentralityPercentile("V0M"); // OB added
730     }
731     
732     if(cl!=fEventClass){
733       // event not in selected event class, reject event
734       if (fDebug > 1) Printf("%s:%d event not in selected event class: event REJECTED ...",(char*)__FILE__,__LINE__);
735       fh1EvtSelection->Fill(2.);
736       PostData(1, fCommonHistList);
737       return;
738     }
739   }
740   
741   // *** vertex cut ***
742   AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
743   Int_t nTracksPrim = primVtx->GetNContributors();
744   fh1VertexNContributors->Fill(nTracksPrim);
745   
746   if (fDebug > 1) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim);
747   if(!nTracksPrim){
748     if (fDebug > 1) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__); 
749     fh1EvtSelection->Fill(3.);
750     PostData(1, fCommonHistList);
751     return;
752   }
753   
754   fh1VertexZ->Fill(primVtx->GetZ());
755   
756   if(TMath::Abs(primVtx->GetZ())>fMaxVertexZ){
757     if (fDebug > 1) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ()); 
758     fh1EvtSelection->Fill(4.);
759     PostData(1, fCommonHistList);
760     return; 
761   }
762   
763   TString primVtxName(primVtx->GetName());
764   
765   if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){
766     if (fDebug > 1) Printf("%s:%d primary vertex selection: TPC vertex, event REJECTED...",(char*)__FILE__,__LINE__);
767     fh1EvtSelection->Fill(5.);
768     PostData(1, fCommonHistList);
769     return;
770   }
771   
772   if (fDebug > 1) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__); 
773   fh1EvtSelection->Fill(0.);
774   fh1EvtCent->Fill(centPercent);
775     
776
777   //___ get MC information __________________________________________________________________
778
779   Double_t ptHard = 0.;
780   Double_t nTrials = 1; // trials for MC trigger weight for real data
781   
782   if(fMCEvent){
783      AliGenEventHeader* genHeader = fMCEvent->GenEventHeader();
784      AliGenPythiaEventHeader*  pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
785      AliGenHijingEventHeader*  hijingGenHeader = 0x0;
786
787      if(pythiaGenHeader){
788          if(fDebug>3) Printf("%s:%d pythiaGenHeader found", (char*)__FILE__,__LINE__);
789          nTrials = pythiaGenHeader->Trials();
790          ptHard  = pythiaGenHeader->GetPtHard();
791
792          fh1PtHard->Fill(ptHard);
793          fh1PtHardTrials->Fill(ptHard,nTrials);
794
795
796      } else { // no pythia, hijing?
797
798          if(fDebug>3) Printf("%s:%d no pythiaGenHeader found", (char*)__FILE__,__LINE__);
799
800          hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
801          if(!hijingGenHeader){
802             Printf("%s:%d no pythiaGenHeader or hjingGenHeader found", (char*)__FILE__,__LINE__);
803          } else {
804             if(fDebug>3) Printf("%s:%d hijingGenHeader found", (char*)__FILE__,__LINE__);
805          }
806      }
807
808      fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
809   }
810
811
812   //____ fetch jets ______________________________________________________________
813
814   Int_t nJCuts = GetListOfJets(fJetsRecCuts, kJetsRecAcceptance);
815   Int_t nRecJetsCuts = 0;
816   if(nJCuts>=0) nRecJetsCuts = fJetsRecCuts->GetEntries();
817   if(fDebug>2)Printf("%s:%d Selected Rec jets after cuts: %d %d",(char*)__FILE__,__LINE__,nJCuts,nRecJetsCuts);
818   if(nRecJetsCuts != nJCuts) Printf("%s:%d Mismatch selected Rec jets after cuts: %d %d",(char*)__FILE__,__LINE__,nJCuts,nRecJetsCuts);
819   fh1nRecJetsCuts->Fill(nRecJetsCuts);
820   
821
822   //____ fetch particles __________________________________________________________
823  
824   Int_t nTCuts = GetListOfTracks(fTracksRecCuts, kTrackAODCuts);
825   if(fDebug>2)Printf("%s:%d Selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,fTracksRecCuts->GetEntries());
826   if(fTracksRecCuts->GetEntries() != nTCuts) 
827     Printf("%s:%d Mismatch selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,fTracksRecCuts->GetEntries());
828   fh1EvtMult->Fill(fTracksRecCuts->GetEntries());
829
830   Int_t nK0s = GetListOfK0s(fListK0s,fK0Type);
831   if(fDebug>2)Printf("%s:%d Selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nK0s,fListK0s->GetEntries());
832   if(nK0s != fListK0s->GetEntries()) Printf("%s:%d Mismatch selected K0s: %d %d",(char*)__FILE__,__LINE__,nK0s,fListK0s->GetEntries());
833   fh1K0Mult->Fill(fListK0s->GetEntries());
834
835    // ___ V0 QA + K0 pt spectra all events _______________________________________________
836   
837   if(fListK0s->GetEntries()>0){
838     
839     for(Int_t it=0; it<fListK0s->GetSize(); ++it){ // loop all K0s 
840       
841       AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListK0s->At(it));
842       
843       Float_t trackPt       = v0->Pt();
844       Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
845       Double_t invM         = v0->MassK0Short();
846       Double_t jetPt        = fFFIMJetPtMin; // assign pro forma jet energy
847
848       fV0QAK0->FillTrackQA(v0->Eta(), TVector2::Phi_0_2pi(v0->Phi()), v0->Pt()); 
849       fFFHistosIMK0AllEvt->FillFF(trackPt, invM, jetPt, incrementJetPt);
850     }
851   }
852
853
854   //____ fill FF histos  __________________________________________________________
855
856   for(Int_t ij=0; ij<nRecJetsCuts; ++ij){
857
858     AliAODJet* jet = (AliAODJet*) (fJetsRecCuts->At(ij));
859     
860     if(ij==0){ // leading jet
861       
862       TList* jettracklist = new TList();
863       Double_t sumPt      = 0.;
864       
865       if(GetFFRadius()<=0){
866         GetJetTracksTrackrefs(jettracklist, jet);
867       } else {
868         GetJetTracksPointing(fTracksRecCuts, jettracklist, jet, GetFFRadius(), sumPt);
869       }
870       
871       for(Int_t it=0; it<jettracklist->GetSize(); ++it){
872
873         AliVParticle* trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));      
874         if(!trackVP)continue;
875
876         Float_t jetPt   = jet->Pt();
877         Float_t trackPt = trackVP->Pt();
878         
879         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
880         
881         fFFHistosRecCuts->FillFF(trackPt, jetPt, incrementJetPt);
882         if(nK0s>0) fFFHistosRecCutsK0Evt->FillFF(trackPt, jetPt, incrementJetPt);
883       }
884       
885       delete jettracklist;
886     }
887   }
888
889
890   for(Int_t ij=0; ij<nRecJetsCuts; ++ij){ // jet loop
891
892     AliAODJet* jet = dynamic_cast<AliAODJet*>(fJetsRecCuts->At(ij));
893     
894     if(ij==0){ // leading jet
895       
896       //fQAJetHistosRecCutsLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt() );
897       
898       Float_t jetPt   = jet->Pt();
899
900       for(Int_t it=0; it<fListK0s->GetSize(); ++it){ // loop all K0s 
901
902         AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListK0s->At(it));
903
904         Double_t v0Mom[3];
905         v0->PxPyPz(v0Mom);
906         TVector3 v0MomVect(v0Mom);
907
908         Double_t dPhiJetK0 = (jet->MomentumVector()->Vect()).DeltaPhi(v0MomVect);
909         
910         Float_t trackPt       = v0->Pt();
911         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
912         Double_t invM         = v0->MassK0Short();
913         
914         fFFHistosIMK0Jet->FillFF(trackPt, invM, jetPt, incrementJetPt);
915         fFFHistosPhiCorrIMK0->FillPhiCorr(trackPt,TVector2::Phi_0_2pi(dPhiJetK0),invM);
916
917         if(dPhiJetK0<fh1dPhiJetK0->GetXaxis()->GetXmin()) dPhiJetK0 += 2*TMath::Pi();
918         fh1dPhiJetK0->Fill(dPhiJetK0);
919
920       }
921       if(fListK0s->GetSize() == 0){ // no K0: increment jet pt spectrum 
922
923         Bool_t incrementJetPt = kTRUE;
924         fFFHistosIMK0Jet->FillFF(-1, -1, jetPt, incrementJetPt);
925       }
926       
927
928       TList* jetConeK0list = new TList();
929       Double_t sumPt   = 0.;
930       
931       GetJetTracksPointing(fListK0s, jetConeK0list, jet, GetFFRadius(), sumPt);
932  
933       if(fDebug>2)Printf("%s:%d nK0s total: %d, in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nK0s,jetConeK0list->GetEntries(),GetFFRadius());
934       
935       for(Int_t it=0; it<jetConeK0list->GetSize(); ++it){ // loop K0s in jet cone
936         
937         AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetConeK0list->At(it));
938         if(!v0) continue;
939
940         Double_t invM           = v0->MassK0Short();
941         Float_t  trackPt        = v0->Pt();
942         Bool_t   incrementJetPt = (it==0) ? kTRUE : kFALSE;
943         
944         //std::cout<<" trackPt "<<trackPt<<" invM "<<invM<<std::endl;
945         fFFHistosIMK0Cone->FillFF(trackPt, invM, jetPt, incrementJetPt);
946       }
947      if(jetConeK0list->GetSize() == 0){ // no K0: increment jet pt spectrum 
948
949         Bool_t incrementJetPt = kTRUE;
950         fFFHistosIMK0Cone->FillFF(-1, -1, jetPt, incrementJetPt);
951       }
952
953       delete jetConeK0list;
954     }
955   }
956   
957   fTracksRecCuts->Clear();
958   fJetsRecCuts->Clear();
959   fListK0s->Clear();
960
961   //Post output data.
962   PostData(1, fCommonHistList);    
963 }
964
965 // ____________________________________________________________________________________________
966 void AliAnalysisTaskJetChem::SetProperties(TH3F* h,const char* x, const char* y, const char* z)
967 {
968   //Set properties of histos (x,y and z title)
969
970   h->SetXTitle(x);
971   h->SetYTitle(y);
972   h->SetZTitle(z);
973   h->GetXaxis()->SetTitleColor(1);
974   h->GetYaxis()->SetTitleColor(1);
975   h->GetZaxis()->SetTitleColor(1);
976 }
977
978 // ____________________________________________________________________________________________________________________________________________
979 Bool_t AliAnalysisTaskJetChem::IsAccepteddEdx(const Double_t mom,const Double_t signal, AliPID::EParticleType n, const Double_t cutnSig) const{
980
981   // apply TPC dE/dx cut similar as in AliTPCpidESD 
982   // note: AliTPCidESD uses inner track param for momentum - not avaiable on AOD,  
983   //       so we use global track momentum 
984   // should use separate parametrisation for MC and data, but probably ALEPH param & 7% resolution used here anyway not the last word 
985  
986  
987   const Double_t kBBMIP(50.);
988   const Double_t kBBRes(0.07);
989   //const Double_t kBBRange(5.);
990   const Double_t kBBp1(0.76176e-1);
991   const Double_t kBBp2(10.632);
992   const Double_t kBBp3(0.13279e-4);
993   const Double_t kBBp4(1.8631);
994   const Double_t kBBp5(1.9479);
995
996   Double_t mass=AliPID::ParticleMass(n); 
997   Double_t betaGamma = mom/mass;
998
999   const Float_t kmeanCorrection =0.1;
1000   Double_t bb = AliExternalTrackParam::BetheBlochAleph(betaGamma,kBBp1,kBBp2,kBBp3,kBBp4,kBBp5);
1001   Double_t meanCorrection =(1+(bb-1)*kmeanCorrection);
1002   Double_t bethe = bb * meanCorrection; // expected
1003   Double_t sigma = bethe * kBBRes;
1004         
1005
1006   Double_t dedx = signal/kBBMIP; // measured
1007
1008   Double_t nSig = (TMath::Abs(dedx - bethe))/sigma;
1009   
1010   if(nSig > cutnSig) return kFALSE; 
1011
1012   return kTRUE;
1013 }
1014
1015 //___________________________________________________________________
1016 Bool_t AliAnalysisTaskJetChem::IsK0InvMass(const Double_t mass) const
1017 {
1018   // K0 mass ? Use FF histo limits
1019   
1020   if(fFFIMInvMMin <= mass && mass <= fFFIMInvMMax) return kTRUE;
1021
1022   return kFALSE;
1023 }
1024
1025 //_____________________________________________________________________________________
1026 Int_t AliAnalysisTaskJetChem::GetListOfK0s(TList *list, const Int_t type)
1027 {
1028   // fill list of V0s selected according to type
1029
1030   if(!list){
1031     if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
1032     return -1;
1033   }
1034
1035   if(type==kTrackUndef) return 0;
1036
1037   for(int i=0; i<fAOD->GetNumberOfV0s(); i++){ // loop over V0s
1038     
1039     AliAODv0* v0 = fAOD->GetV0(i);
1040
1041     Bool_t isOnFly = v0->GetOnFlyStatus();
1042     
1043     if(!isOnFly &&  (type == kOnFly || type == kOnFlyPID || type == kOnFlydEdx || type == kOnFlyPrim)) continue; 
1044     if( isOnFly &&  (type == kOffl  || type == kOfflPID  || type == kOffldEdx  || type == kOfflPrim))  continue; 
1045
1046     Double_t massK0 = v0->MassK0Short();
1047
1048     if(!(IsK0InvMass(massK0))) continue; // moved invMass cut for HI - otherwise too slow
1049
1050     if(type == kOnFlyPID || type == kOfflPID){
1051       
1052       AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0)); // slow 
1053       AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1)); // slow  
1054      
1055       // AOD pid - cuts strongly into signal
1056
1057       AliAODTrack::AODTrkPID_t mpPIDNeg = trackNeg->GetMostProbablePID();
1058       AliAODTrack::AODTrkPID_t mpPIDPos = trackPos->GetMostProbablePID();
1059         
1060       if(!( (mpPIDNeg == AliAODTrack::kPion) && (mpPIDPos == AliAODTrack::kPion) ) ) continue;
1061     }
1062    
1063     if(type == kOnFlydEdx || type == kOffldEdx){
1064
1065       AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0)); // slow 
1066       AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1)); // slow  
1067
1068       AliAODPid*  aodPidPos = trackPos->GetDetPid();
1069       AliAODPid*  aodPidNeg = trackNeg->GetDetPid();
1070
1071       Double_t  dEdxPos = aodPidPos->GetTPCsignal();
1072       Double_t  dEdxNeg = aodPidNeg->GetTPCsignal();
1073
1074       Double_t momPos  = trackPos->P();
1075       Double_t momNeg  = trackNeg->P();
1076
1077       Int_t cutnSigdEdx = 2;
1078       if(! (IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,cutnSigdEdx)) && (IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,cutnSigdEdx)) ) continue;
1079       
1080     }   
1081     
1082     if(type == kOnFlyPrim || type == kOfflPrim){
1083       
1084       AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0)); // slow 
1085       AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1)); // slow  
1086       
1087       //std::cout<<"  filer map trackPos "<<trackPos->GetFilterMap()<<" trackNeg "<<trackNeg->GetFilterMap()<<std::endl;
1088       
1089       if((fFilterMaskK0>0) && !(trackPos->TestFilterBit(fFilterMaskK0)))   continue;
1090       if((fFilterMaskK0>0) && !(trackNeg->TestFilterBit(fFilterMaskK0)))   continue;
1091     }
1092     
1093     list->Add(v0);
1094   }
1095   
1096   Int_t nK0s = list->GetSize();
1097   
1098   return nK0s;
1099 }
1100
1101