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       if(!v0) continue;
 843
 844       Float_t trackPt       = v0->Pt();
 845       Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
 846       Double_t invM         = v0->MassK0Short();
 847       Double_t jetPt        = fFFIMJetPtMin; // assign pro forma jet energy
 848
 849       fV0QAK0->FillTrackQA(v0->Eta(), TVector2::Phi_0_2pi(v0->Phi()), v0->Pt());
 850       fFFHistosIMK0AllEvt->FillFF(trackPt, invM, jetPt, incrementJetPt);
 851     }
 852   }
 853
 854
 855   //____ fill FF histos  __________________________________________________________
 856
 857   for(Int_t ij=0; ij<nRecJetsCuts; ++ij){
 858
 859     AliAODJet* jet = (AliAODJet*) (fJetsRecCuts->At(ij));
 860
 861     if(ij==0){ // leading jet
 862
 863       TList* jettracklist = new TList();
 864       Double_t sumPt      = 0.;
 865
 866       if(GetFFRadius()<=0){
 867         GetJetTracksTrackrefs(jettracklist, jet);
 868       } else {
 869         GetJetTracksPointing(fTracksRecCuts, jettracklist, jet, GetFFRadius(), sumPt);
 870       }
 871
 872       for(Int_t it=0; it<jettracklist->GetSize(); ++it){
 873
 874         AliVParticle* trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));
 875         if(!trackVP)continue;
 876
 877         Float_t jetPt   = jet->Pt();
 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   }
 889
 890
 891   for(Int_t ij=0; ij<nRecJetsCuts; ++ij){ // jet loop
 892
 893     AliAODJet* jet = dynamic_cast<AliAODJet*>(fJetsRecCuts->At(ij));
 894     if(!jet) continue;
 895
 896     if(ij==0){ // leading jet
 897
 898       //fQAJetHistosRecCutsLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt() );
 899
 900       Float_t jetPt   = jet->Pt();
 901
 902       for(Int_t it=0; it<fListK0s->GetSize(); ++it){ // loop all K0s
 903
 904         AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListK0s->At(it));
 905         if(!v0) continue;
 906
 907         Double_t v0Mom[3];
 908         v0->PxPyPz(v0Mom);
 909         TVector3 v0MomVect(v0Mom);
 910
 911         Double_t dPhiJetK0 = (jet->MomentumVector()->Vect()).DeltaPhi(v0MomVect);
 912
 913         Float_t trackPt       = v0->Pt();
 914         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
 915         Double_t invM         = v0->MassK0Short();
 916
 917         fFFHistosIMK0Jet->FillFF(trackPt, invM, jetPt, incrementJetPt);
 918         fFFHistosPhiCorrIMK0->FillPhiCorr(trackPt,TVector2::Phi_0_2pi(dPhiJetK0),invM);
 919
 920         if(dPhiJetK0<fh1dPhiJetK0->GetXaxis()->GetXmin()) dPhiJetK0 += 2*TMath::Pi();
 921         fh1dPhiJetK0->Fill(dPhiJetK0);
 922
 923       }
 924       if(fListK0s->GetSize() == 0){ // no K0: increment jet pt spectrum
 925
 926         Bool_t incrementJetPt = kTRUE;
 927         fFFHistosIMK0Jet->FillFF(-1, -1, jetPt, incrementJetPt);
 928       }
 929
 930
 931       TList* jetConeK0list = new TList();
 932       Double_t sumPt   = 0.;
 933
 934       GetJetTracksPointing(fListK0s, jetConeK0list, jet, GetFFRadius(), sumPt);
 935
 936       if(fDebug>2)Printf("%s:%d nK0s total: %d, in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nK0s,jetConeK0list->GetEntries(),GetFFRadius());
 937
 938       for(Int_t it=0; it<jetConeK0list->GetSize(); ++it){ // loop K0s in jet cone
 939
 940         AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetConeK0list->At(it));
 941         if(!v0) continue;
 942
 943         Double_t invM           = v0->MassK0Short();
 944         Float_t  trackPt        = v0->Pt();
 945         Bool_t   incrementJetPt = (it==0) ? kTRUE : kFALSE;
 946
 947         //std::cout<<" trackPt "<<trackPt<<" invM "<<invM<<std::endl;
 948         fFFHistosIMK0Cone->FillFF(trackPt, invM, jetPt, incrementJetPt);
 949       }
 950      if(jetConeK0list->GetSize() == 0){ // no K0: increment jet pt spectrum
 951
 952         Bool_t incrementJetPt = kTRUE;
 953         fFFHistosIMK0Cone->FillFF(-1, -1, jetPt, incrementJetPt);
 954       }
 955
 956       delete jetConeK0list;
 957     }
 958   }
 959
 960   fTracksRecCuts->Clear();
 961   fJetsRecCuts->Clear();
 962   fListK0s->Clear();
 963
 964   //Post output data.
 965   PostData(1, fCommonHistList);
 966 }
 967
 968 // ____________________________________________________________________________________________
 969 void AliAnalysisTaskJetChem::SetProperties(TH3F* h,const char* x, const char* y, const char* z)
 970 {
 971   //Set properties of histos (x,y and z title)
 972
 973   h->SetXTitle(x);
 974   h->SetYTitle(y);
 975   h->SetZTitle(z);
 976   h->GetXaxis()->SetTitleColor(1);
 977   h->GetYaxis()->SetTitleColor(1);
 978   h->GetZaxis()->SetTitleColor(1);
 979 }
 980
 981 // ____________________________________________________________________________________________________________________________________________
 982 Bool_t AliAnalysisTaskJetChem::IsAccepteddEdx(const Double_t mom,const Double_t signal, AliPID::EParticleType n, const Double_t cutnSig) const{
 983
 984   // apply TPC dE/dx cut similar as in AliTPCpidESD
 985   // note: AliTPCidESD uses inner track param for momentum - not avaiable on AOD,
 986   //       so we use global track momentum
 987   // should use separate parametrisation for MC and data, but probably ALEPH param & 7% resolution used here anyway not the last word
 988
 989
 990   const Double_t kBBMIP(50.);
 991   const Double_t kBBRes(0.07);
 992   //const Double_t kBBRange(5.);
 993   const Double_t kBBp1(0.76176e-1);
 994   const Double_t kBBp2(10.632);
 995   const Double_t kBBp3(0.13279e-4);
 996   const Double_t kBBp4(1.8631);
 997   const Double_t kBBp5(1.9479);
 998
 999   Double_t mass=AliPID::ParticleMass(n);
1000   Double_t betaGamma = mom/mass;
1001
1002   const Float_t kmeanCorrection =0.1;
1003   Double_t bb = AliExternalTrackParam::BetheBlochAleph(betaGamma,kBBp1,kBBp2,kBBp3,kBBp4,kBBp5);
1004   Double_t meanCorrection =(1+(bb-1)*kmeanCorrection);
1005   Double_t bethe = bb * meanCorrection; // expected
1006   Double_t sigma = bethe * kBBRes;
1007
1008
1009   Double_t dedx = signal/kBBMIP; // measured
1010
1011   Double_t nSig = (TMath::Abs(dedx - bethe))/sigma;
1012
1013   if(nSig > cutnSig) return kFALSE;
1014
1015   return kTRUE;
1016 }
1017
1018 //___________________________________________________________________
1019 Bool_t AliAnalysisTaskJetChem::IsK0InvMass(const Double_t mass) const
1020 {
1021   // K0 mass ? Use FF histo limits
1022
1023   if(fFFIMInvMMin <= mass && mass <= fFFIMInvMMax) return kTRUE;
1024
1025   return kFALSE;
1026 }
1027
1028 //_____________________________________________________________________________________
1029 Int_t AliAnalysisTaskJetChem::GetListOfK0s(TList *list, const Int_t type)
1030 {
1031   // fill list of V0s selected according to type
1032
1033   if(!list){
1034     if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
1035     return -1;
1036   }
1037
1038   if(type==kTrackUndef) return 0;
1039
1040   for(int i=0; i<fAOD->GetNumberOfV0s(); i++){ // loop over V0s
1041
1042     AliAODv0* v0 = fAOD->GetV0(i);
1043
1044     Bool_t isOnFly = v0->GetOnFlyStatus();
1045
1046     if(!isOnFly &&  (type == kOnFly || type == kOnFlyPID || type == kOnFlydEdx || type == kOnFlyPrim)) continue;
1047     if( isOnFly &&  (type == kOffl  || type == kOfflPID  || type == kOffldEdx  || type == kOfflPrim))  continue;
1048
1049     Double_t massK0 = v0->MassK0Short();
1050
1051     if(!(IsK0InvMass(massK0))) continue; // moved invMass cut for HI - otherwise too slow
1052
1053     if(type == kOnFlyPID || type == kOfflPID){
1054
1055       AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0)); // slow
1056       AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1)); // slow
1057
1058       // AOD pid - cuts strongly into signal
1059
1060       AliAODTrack::AODTrkPID_t mpPIDNeg = trackNeg->GetMostProbablePID();
1061       AliAODTrack::AODTrkPID_t mpPIDPos = trackPos->GetMostProbablePID();
1062
1063       if(!( (mpPIDNeg == AliAODTrack::kPion) && (mpPIDPos == AliAODTrack::kPion) ) ) continue;
1064     }
1065
1066     if(type == kOnFlydEdx || type == kOffldEdx){
1067
1068       AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0)); // slow
1069       AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1)); // slow
1070
1071       AliAODPid*  aodPidPos = trackPos->GetDetPid();
1072       AliAODPid*  aodPidNeg = trackNeg->GetDetPid();
1073
1074       Double_t  dEdxPos = aodPidPos->GetTPCsignal();
1075       Double_t  dEdxNeg = aodPidNeg->GetTPCsignal();
1076
1077       Double_t momPos  = trackPos->P();
1078       Double_t momNeg  = trackNeg->P();
1079
1080       Int_t cutnSigdEdx = 2;
1081       if(! (IsAccepteddEdx(momPos,dEdxPos,AliPID::kPion,cutnSigdEdx)) && (IsAccepteddEdx(momNeg,dEdxNeg,AliPID::kPion,cutnSigdEdx)) ) continue;
1082
1083     }
1084
1085     if(type == kOnFlyPrim || type == kOfflPrim){
1086
1087       AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0)); // slow
1088       AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1)); // slow
1089
1090       //std::cout<<"  filer map trackPos "<<trackPos->GetFilterMap()<<" trackNeg "<<trackNeg->GetFilterMap()<<std::endl;
1091
1092       if((fFilterMaskK0>0) && !(trackPos->TestFilterBit(fFilterMaskK0)))   continue;
1093       if((fFilterMaskK0>0) && !(trackNeg->TestFilterBit(fFilterMaskK0)))   continue;
1094     }
1095
1096     list->Add(v0);
1097   }
1098
1099   Int_t nK0s = list->GetSize();
1100
1101   return nK0s;
1102 }
1103
1104