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