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