PWGJE/StrangenessInJets/AliAnalysisTaskJetChem.cxx

   1  /*************************************************************************
   2  *                                                                       *
   3  *                                                                       *
   4  *      Task for Jet Chemistry Analysis in PWG4 Jet Task Force Train     *
   5  *    Analysis of K0s, Lambda and Antilambda with and without Jetevents  *
   6  *                                                                       *
   7  *************************************************************************/
   8
   9 /**************************************************************************
  10  * Copyright(c) 1998-2012, 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 grante    *
  17  *                                                                        *
  18  * without fee, provided that the above copyright notice appears in all   *
  19  * copies and that both the copyright notice and this permission notice   *
  20  * appear in the supporting documentation. The authors make no claims     *
  21  * about the suitability of this software for any purpose. It is          *
  22  * provided "as is" without express or implied warranty.                  *
  23  **************************************************************************/
  24
  25
  26 /* $Id: */
  27
  28 #include <iostream>
  29 #include "TH2.h"
  30 #include "TH3.h"
  31 #include "TH2F.h"
  32 #include "TH3F.h"
  33 #include "TH2D.h"
  34 #include "TH3D.h"
  35 #include "TChain.h"
  36 #include "TTree.h"
  37 #include "TList.h"
  38 #include "TCanvas.h"
  39 #include "TPDGCode.h"
  40 #include "TProfile.h"
  41 #include "THnSparse.h"
  42 #include <algorithm>
  43 #include <string>
  44 #include "AliAnalysisHelperJetTasks.h"
  45 #include "TDatabasePDG.h"
  46 #include "TPDGCode.h"
  47 #include "AliAnalysisManager.h"
  48 #include "AliAODHandler.h"
  49 #include "AliAODInputHandler.h"
  50 #include "AliESDEvent.h"
  51 #include "AliGenPythiaEventHeader.h"
  52 #include "AliGenHijingEventHeader.h"
  53 #include "AliGenEventHeader.h"
  54 #include "TLorentzVector.h"
  55 #include "AliAODEvent.h"
  56 #include "AliAODJet.h"
  57 #include "AliAODv0.h"
  58 #include "AliAODTrack.h"
  59 #include "AliCentrality.h"
  60 #include "AliAnalysisTaskSE.h"
  61 #include "AliESDtrack.h"
  62 #include "AliESDtrackCuts.h"
  63 #include "AliESDEvent.h"
  64 #include "AliESDInputHandler.h"
  65 #include "AliPID.h"
  66 #include "AliPIDResponse.h"
  67 #include "AliAODPid.h"
  68 #include "AliExternalTrackParam.h"
  69 #include "AliAnalysisTaskJetChem.h"
  70 #include "AliPhysicsSelection.h"
  71 #include "AliBackgroundSelection.h"
  72 #include "AliInputEventHandler.h"
  73 #include "AliAODMCHeader.h"
  74 #include "AliAODPid.h"
  75 #include "AliVEvent.h"
  76 #include "AliAODMCParticle.h"
  77 #include "TVector3.h"
  78 #include "TRandom.h"
  79
  80 ClassImp(AliAnalysisTaskJetChem)
  81
  82 //____________________________________________________________________________
  83 AliAnalysisTaskJetChem::AliAnalysisTaskJetChem()
  84    : AliAnalysisTaskFragmentationFunction()
  85
  86    ,fAnalysisMC(0)
  87    ,fDeltaVertexZ(0)
  88    ,fCuttrackNegNcls(0)
  89    ,fCuttrackPosNcls(0)
  90    ,fCutPostrackRap(0)
  91    ,fCutNegtrackRap(0)
  92    ,fCutRap(0)
  93    ,fCutPostrackEta(0)
  94    ,fCutNegtrackEta(0)
  95    ,fCutEta(0)
  96    ,fCutV0cosPointAngle(0)
  97    ,fCutChi2PosDaughter(0)
  98    ,fCutChi2NegDaughter(0)
  99    ,fKinkDaughters(0)
 100    ,fRequireTPCRefit(0)
 101    ,fCutArmenteros(0)
 102    ,fCutV0DecayMin(0)
 103    ,fCutV0DecayMax(0)
 104    ,fCutV0totMom(0)
 105    ,fCutDcaV0Daughters(0)
 106    ,fCutDcaPosToPrimVertex(0)
 107    ,fCutDcaNegToPrimVertex(0)
 108    ,fCutV0RadiusMin(0)
 109    ,fCutV0RadiusMax(0)
 110    ,fCutBetheBloch(0)
 111    ,fCutRatio(0)
 112    ,fK0Type(0)
 113    ,fFilterMaskK0(0)
 114    ,fListK0s(0)
 115    ,fPIDResponse(0)
 116    ,fV0QAK0(0)
 117    ,fFFHistosRecCutsK0Evt(0)
 118    ,fFFHistosIMK0AllEvt(0)
 119    ,fFFHistosIMK0Jet(0)
 120    ,fFFHistosIMK0Cone(0)
 121    ,fFFHistosPhiCorrIMK0(0)
 122    ,fLaType(0)
 123    ,fFilterMaskLa(0)
 124    ,fListLa(0)
 125    ,fFFHistosIMLaAllEvt(0)
 126    ,fFFHistosIMLaJet(0)
 127    ,fFFHistosIMLaCone(0)
 128    ,fFFHistosPhiCorrIMLa(0)
 129    ,fALaType(0)
 130    ,fFilterMaskALa(0)
 131    ,fListALa(0)
 132    ,fListFeeddownLaCand(0)
 133    ,fListFeeddownALaCand(0)
 134    ,jetConeFDLalist(0)
 135    ,jetConeFDALalist(0)
 136    ,fListMCgenK0s(0)
 137    ,fListMCgenLa(0)
 138    ,fListMCgenALa(0)
 139    ,fListMCgenK0sCone(0)
 140    ,fListMCgenLaCone(0)
 141    ,fListMCgenALaCone(0)
 142    ,IsArmenterosSelected(0)
 143    ,fFFHistosIMALaAllEvt(0)
 144    ,fFFHistosIMALaJet(0)
 145    ,fFFHistosIMALaCone(0)
 146    ,fFFHistosPhiCorrIMALa(0)
 147    ,fFFIMNBinsJetPt(0)
 148    ,fFFIMJetPtMin(0)
 149    ,fFFIMJetPtMax(0)
 150    ,fFFIMNBinsInvM(0)
 151    ,fFFIMInvMMin(0)
 152    ,fFFIMInvMMax(0)
 153    ,fFFIMNBinsPt(0)
 154    ,fFFIMPtMin(0)
 155    ,fFFIMPtMax(0)
 156    ,fFFIMNBinsXi(0)
 157    ,fFFIMXiMin(0)
 158    ,fFFIMXiMax(0)
 159    ,fFFIMNBinsZ(0)
 160    ,fFFIMZMin(0)
 161    ,fFFIMZMax(0)
 162    ,fFFIMLaNBinsJetPt(0)
 163    ,fFFIMLaJetPtMin(0)
 164    ,fFFIMLaJetPtMax(0)
 165    ,fFFIMLaNBinsInvM(0)
 166    ,fFFIMLaInvMMin(0)
 167    ,fFFIMLaInvMMax(0)
 168    ,fFFIMLaNBinsPt(0)
 169    ,fFFIMLaPtMin(0)
 170    ,fFFIMLaPtMax(0)
 171    ,fFFIMLaNBinsXi(0)
 172    ,fFFIMLaXiMin(0)
 173    ,fFFIMLaXiMax(0)
 174    ,fFFIMLaNBinsZ(0)
 175    ,fFFIMLaZMin(0)
 176    ,fFFIMLaZMax(0)
 177    ,fPhiCorrIMNBinsPt(0)
 178    ,fPhiCorrIMPtMin(0)
 179    ,fPhiCorrIMPtMax(0)
 180    ,fPhiCorrIMNBinsPhi(0)
 181    ,fPhiCorrIMPhiMin(0)
 182    ,fPhiCorrIMPhiMax(0)
 183    ,fPhiCorrIMNBinsInvM(0)
 184    ,fPhiCorrIMInvMMin(0)
 185    ,fPhiCorrIMInvMMax(0)
 186    ,fPhiCorrIMLaNBinsPt(0)
 187    ,fPhiCorrIMLaPtMin(0)
 188    ,fPhiCorrIMLaPtMax(0)
 189    ,fPhiCorrIMLaNBinsPhi(0)
 190    ,fPhiCorrIMLaPhiMin(0)
 191    ,fPhiCorrIMLaPhiMax(0)
 192    ,fPhiCorrIMLaNBinsInvM(0)
 193    ,fPhiCorrIMLaInvMMin(0)
 194    ,fPhiCorrIMLaInvMMax(0)
 195    ,fh1EvtAllCent(0)
 196    ,fh1Evt(0)
 197    ,fh1K0Mult(0)
 198    ,fh1dPhiJetK0(0)
 199    ,fh1LaMult(0)
 200    ,fh1dPhiJetLa(0)
 201    ,fh1ALaMult(0)
 202    ,fh1dPhiJetALa(0)
 203    ,fh1JetEta(0)
 204    ,fh1JetPhi(0)
 205    ,fh2JetEtaPhi(0)
 206    ,fh1V0JetPt(0)
 207    ,fh2FFJetTrackEta(0)
 208    ,fh1trackPosNCls(0)
 209    ,fh1trackNegNCls(0)
 210    ,fh1trackPosRap(0)
 211    ,fh1trackNegRap(0)
 212    ,fh1V0Rap(0)
 213    ,fh1trackPosEta(0)
 214    ,fh1trackNegEta(0)
 215    ,fh1V0Eta(0)
 216    ,fh1V0totMom(0)
 217    ,fh1CosPointAngle(0)
 218    ,fh1Chi2Pos(0)
 219    ,fh1Chi2Neg(0)
 220    ,fh1DecayLengthV0(0)
 221    ,fh2ProperLifetimeK0sVsPtBeforeCut(0)
 222    ,fh2ProperLifetimeK0sVsPtAfterCut(0)
 223    ,fh1ProperLifetimeV0BeforeCut(0)
 224    ,fh1ProperLifetimeV0AfterCut(0)
 225    ,fh1V0Radius(0)
 226    ,fh1DcaV0Daughters(0)
 227    ,fh1DcaPosToPrimVertex(0)
 228    ,fh1DcaNegToPrimVertex(0)
 229    ,fh2ArmenterosBeforeCuts(0)
 230    ,fh2ArmenterosAfterCuts(0)
 231    ,fh2BBLaPos(0)
 232    ,fh2BBLaNeg(0)
 233    ,fh1CrossedRowsOverFindableNeg(0)
 234    ,fh1CrossedRowsOverFindablePos(0)
 235    ,fh1PosDaughterCharge(0)
 236    ,fh1NegDaughterCharge(0)
 237    ,fh1PtMCK0s(0)
 238    ,fh1PtMCLa(0)
 239    ,fh1PtMCALa(0)
 240    ,fh1EtaK0s(0)
 241    ,fh1EtaLa(0)
 242    ,fh1EtaALa(0)
 243    ,fh3InvMassEtaTrackPtK0s(0)
 244    ,fh3InvMassEtaTrackPtLa(0)
 245    ,fh3InvMassEtaTrackPtALa(0)
 246
 247    ,fh1TrackMultCone(0)
 248    ,fh2TrackMultCone(0)
 249    ,fh2MCgenK0Cone(0)
 250    ,fh2MCgenLaCone(0)
 251    ,fh2MCgenALaCone(0)
 252    ,fh2MCEtagenK0Cone(0)
 253    ,fh2MCEtagenLaCone(0)
 254    ,fh2MCEtagenALaCone(0)
 255    ,fh1FFIMK0ConeSmear(0)
 256    ,fh1FFIMLaConeSmear(0)
 257    ,fh1FFIMALaConeSmear(0)
 258    ,fh3MCrecK0Cone(0)
 259    ,fh3MCrecLaCone(0)
 260    ,fh3MCrecALaCone(0)
 261    ,fh3MCrecK0ConeSmear(0)
 262    ,fh3MCrecLaConeSmear(0)
 263    ,fh3MCrecALaConeSmear(0)
 264    ,fh3SecContinCone(0)
 265    ,fh3StrContinCone(0)
 266    ,fh3IMK0PerpCone(0)
 267    ,fh3IMLaPerpCone(0)
 268    ,fh3IMALaPerpCone(0)
 269    ,fh3IMK0MedianCone(0)
 270    ,fh3IMLaMedianCone(0)
 271    ,fh3IMALaMedianCone(0)
 272    ,fh1MCMultiplicityPrimary(0)
 273    ,fh1MCMultiplicityTracks(0)
 274    ,fh1MCmotherLa(0)
 275    ,fh1MCmotherALa(0)
 276    ,fh3FeedDownLa(0)
 277    ,fh3FeedDownALa(0)
 278    ,fh3FeedDownLaCone(0)
 279    ,fh3FeedDownALaCone(0)
 280    ,fh1MCProdRadiusK0s(0)
 281    ,fh1MCProdRadiusLambda(0)
 282    ,fh1MCProdRadiusAntiLambda(0)
 283    ,fh1MCPtV0s(0)
 284    ,fh1MCPtK0s(0)
 285    ,fh1MCPtLambda(0)
 286    ,fh1MCPtAntiLambda(0)
 287    ,fh1MCXiPt(0)
 288    ,fh1MCXibarPt(0)
 289    ,fh2MCEtaVsPtK0s(0)
 290    ,fh2MCEtaVsPtLa(0)
 291    ,fh2MCEtaVsPtALa(0)
 292    ,fh1MCRapK0s(0)
 293    ,fh1MCRapLambda(0)
 294    ,fh1MCRapAntiLambda(0)
 295    ,fh1MCEtaAllK0s(0)
 296    ,fh1MCEtaK0s(0)
 297    ,fh1MCEtaLambda(0)
 298    ,fh1MCEtaAntiLambda(0)
 299
 300 {
 301    // default constructor
 302 }
 303
 304 //__________________________________________________________________________________________
 305 AliAnalysisTaskJetChem::AliAnalysisTaskJetChem(const char *name)
 306   : AliAnalysisTaskFragmentationFunction(name)
 307
 308   ,fAnalysisMC(0)
 309   ,fDeltaVertexZ(0)
 310   ,fCuttrackNegNcls(0)
 311   ,fCuttrackPosNcls(0)
 312   ,fCutPostrackRap(0)
 313   ,fCutNegtrackRap(0)
 314   ,fCutRap(0)
 315   ,fCutPostrackEta(0)
 316   ,fCutNegtrackEta(0)
 317   ,fCutEta(0)
 318   ,fCutV0cosPointAngle(0)
 319   ,fCutChi2PosDaughter(0)
 320   ,fCutChi2NegDaughter(0)
 321   ,fKinkDaughters(0)
 322   ,fRequireTPCRefit(0)
 323   ,fCutArmenteros(0)
 324   ,fCutV0DecayMin(0)
 325   ,fCutV0DecayMax(0)
 326   ,fCutV0totMom(0)
 327   ,fCutDcaV0Daughters(0)
 328   ,fCutDcaPosToPrimVertex(0)
 329   ,fCutDcaNegToPrimVertex(0)
 330   ,fCutV0RadiusMin(0)
 331   ,fCutV0RadiusMax(0)
 332   ,fCutBetheBloch(0)
 333   ,fCutRatio(0)
 334   ,fK0Type(0)
 335   ,fFilterMaskK0(0)
 336   ,fListK0s(0)
 337   ,fPIDResponse(0)
 338   ,fV0QAK0(0)
 339   ,fFFHistosRecCutsK0Evt(0)
 340   ,fFFHistosIMK0AllEvt(0)
 341   ,fFFHistosIMK0Jet(0)
 342   ,fFFHistosIMK0Cone(0)
 343   ,fFFHistosPhiCorrIMK0(0)
 344   ,fLaType(0)
 345   ,fFilterMaskLa(0)
 346   ,fListLa(0)
 347   ,fFFHistosIMLaAllEvt(0)
 348   ,fFFHistosIMLaJet(0)
 349   ,fFFHistosIMLaCone(0)
 350   ,fFFHistosPhiCorrIMLa(0)
 351   ,fALaType(0)
 352   ,fFilterMaskALa(0)
 353   ,fListALa(0)
 354   ,fListFeeddownLaCand(0)
 355   ,fListFeeddownALaCand(0)
 356   ,jetConeFDLalist(0)
 357   ,jetConeFDALalist(0)
 358   ,fListMCgenK0s(0)
 359   ,fListMCgenLa(0)
 360   ,fListMCgenALa(0)
 361   ,fListMCgenK0sCone(0)
 362   ,fListMCgenLaCone(0)
 363   ,fListMCgenALaCone(0)
 364   ,IsArmenterosSelected(0)
 365   ,fFFHistosIMALaAllEvt(0)
 366   ,fFFHistosIMALaJet(0)
 367   ,fFFHistosIMALaCone(0)
 368   ,fFFHistosPhiCorrIMALa(0)
 369   ,fFFIMNBinsJetPt(0)
 370   ,fFFIMJetPtMin(0)
 371   ,fFFIMJetPtMax(0)
 372   ,fFFIMNBinsInvM(0)
 373   ,fFFIMInvMMin(0)
 374   ,fFFIMInvMMax(0)
 375   ,fFFIMNBinsPt(0)
 376   ,fFFIMPtMin(0)
 377   ,fFFIMPtMax(0)
 378   ,fFFIMNBinsXi(0)
 379   ,fFFIMXiMin(0)
 380   ,fFFIMXiMax(0)
 381   ,fFFIMNBinsZ(0)
 382   ,fFFIMZMin(0)
 383   ,fFFIMZMax(0)
 384   ,fFFIMLaNBinsJetPt(0)
 385   ,fFFIMLaJetPtMin(0)
 386   ,fFFIMLaJetPtMax(0)
 387   ,fFFIMLaNBinsInvM(0)
 388   ,fFFIMLaInvMMin(0)
 389   ,fFFIMLaInvMMax(0)
 390   ,fFFIMLaNBinsPt(0)
 391   ,fFFIMLaPtMin(0)
 392   ,fFFIMLaPtMax(0)
 393   ,fFFIMLaNBinsXi(0)
 394   ,fFFIMLaXiMin(0)
 395   ,fFFIMLaXiMax(0)
 396   ,fFFIMLaNBinsZ(0)
 397   ,fFFIMLaZMin(0)
 398   ,fFFIMLaZMax(0)
 399   ,fPhiCorrIMNBinsPt(0)
 400   ,fPhiCorrIMPtMin(0)
 401   ,fPhiCorrIMPtMax(0)
 402   ,fPhiCorrIMNBinsPhi(0)
 403   ,fPhiCorrIMPhiMin(0)
 404   ,fPhiCorrIMPhiMax(0)
 405   ,fPhiCorrIMNBinsInvM(0)
 406   ,fPhiCorrIMInvMMin(0)
 407   ,fPhiCorrIMInvMMax(0)
 408   ,fPhiCorrIMLaNBinsPt(0)
 409   ,fPhiCorrIMLaPtMin(0)
 410   ,fPhiCorrIMLaPtMax(0)
 411   ,fPhiCorrIMLaNBinsPhi(0)
 412   ,fPhiCorrIMLaPhiMin(0)
 413   ,fPhiCorrIMLaPhiMax(0)
 414   ,fPhiCorrIMLaNBinsInvM(0)
 415   ,fPhiCorrIMLaInvMMin(0)
 416   ,fPhiCorrIMLaInvMMax(0)
 417   ,fh1EvtAllCent(0)
 418   ,fh1Evt(0)
 419   ,fh1K0Mult(0)
 420   ,fh1dPhiJetK0(0)
 421   ,fh1LaMult(0)
 422   ,fh1dPhiJetLa(0)
 423   ,fh1ALaMult(0)
 424   ,fh1dPhiJetALa(0)
 425   ,fh1JetEta(0)
 426   ,fh1JetPhi(0)
 427   ,fh2JetEtaPhi(0)
 428   ,fh1V0JetPt(0)
 429   ,fh2FFJetTrackEta(0)
 430   ,fh1trackPosNCls(0)
 431   ,fh1trackNegNCls(0)
 432   ,fh1trackPosRap(0)
 433   ,fh1trackNegRap(0)
 434   ,fh1V0Rap(0)
 435   ,fh1trackPosEta(0)
 436   ,fh1trackNegEta(0)
 437   ,fh1V0Eta(0)
 438   ,fh1V0totMom(0)
 439   ,fh1CosPointAngle(0)
 440   ,fh1Chi2Pos(0)
 441   ,fh1Chi2Neg(0)
 442   ,fh1DecayLengthV0(0)
 443   ,fh2ProperLifetimeK0sVsPtBeforeCut(0)
 444   ,fh2ProperLifetimeK0sVsPtAfterCut(0)
 445   ,fh1ProperLifetimeV0BeforeCut(0)
 446   ,fh1ProperLifetimeV0AfterCut(0)
 447   ,fh1V0Radius(0)
 448   ,fh1DcaV0Daughters(0)
 449   ,fh1DcaPosToPrimVertex(0)
 450   ,fh1DcaNegToPrimVertex(0)
 451   ,fh2ArmenterosBeforeCuts(0)
 452   ,fh2ArmenterosAfterCuts(0)
 453   ,fh2BBLaPos(0)
 454   ,fh2BBLaNeg(0)
 455   ,fh1CrossedRowsOverFindableNeg(0)
 456   ,fh1CrossedRowsOverFindablePos(0)
 457   ,fh1PosDaughterCharge(0)
 458   ,fh1NegDaughterCharge(0)
 459   ,fh1PtMCK0s(0)
 460   ,fh1PtMCLa(0)
 461   ,fh1PtMCALa(0)
 462   ,fh1EtaK0s(0)
 463   ,fh1EtaLa(0)
 464   ,fh1EtaALa(0)
 465   ,fh3InvMassEtaTrackPtK0s(0)
 466   ,fh3InvMassEtaTrackPtLa(0)
 467   ,fh3InvMassEtaTrackPtALa(0)
 468   ,fh1TrackMultCone(0)
 469   ,fh2TrackMultCone(0)
 470   ,fh2MCgenK0Cone(0)
 471   ,fh2MCgenLaCone(0)
 472   ,fh2MCgenALaCone(0)
 473   ,fh2MCEtagenK0Cone(0)
 474   ,fh2MCEtagenLaCone(0)
 475   ,fh2MCEtagenALaCone(0)
 476   ,fh1FFIMK0ConeSmear(0)
 477   ,fh1FFIMLaConeSmear(0)
 478   ,fh1FFIMALaConeSmear(0)
 479   ,fh3MCrecK0Cone(0)
 480   ,fh3MCrecLaCone(0)
 481   ,fh3MCrecALaCone(0)
 482   ,fh3MCrecK0ConeSmear(0)
 483   ,fh3MCrecLaConeSmear(0)
 484   ,fh3MCrecALaConeSmear(0)
 485   ,fh3SecContinCone(0)
 486   ,fh3StrContinCone(0)
 487   ,fh3IMK0PerpCone(0)
 488   ,fh3IMLaPerpCone(0)
 489   ,fh3IMALaPerpCone(0)
 490   ,fh3IMK0MedianCone(0)
 491   ,fh3IMLaMedianCone(0)
 492   ,fh3IMALaMedianCone(0)
 493   ,fh1MCMultiplicityPrimary(0)
 494   ,fh1MCMultiplicityTracks(0)
 495   ,fh1MCmotherLa(0)
 496   ,fh1MCmotherALa(0)
 497   ,fh3FeedDownLa(0)
 498   ,fh3FeedDownALa(0)
 499   ,fh3FeedDownLaCone(0)
 500   ,fh3FeedDownALaCone(0)
 501   ,fh1MCProdRadiusK0s(0)
 502   ,fh1MCProdRadiusLambda(0)
 503   ,fh1MCProdRadiusAntiLambda(0)
 504   ,fh1MCPtV0s(0)
 505   ,fh1MCPtK0s(0)
 506   ,fh1MCPtLambda(0)
 507   ,fh1MCPtAntiLambda(0)
 508   ,fh1MCXiPt(0)
 509   ,fh1MCXibarPt(0)
 510   ,fh2MCEtaVsPtK0s(0)
 511   ,fh2MCEtaVsPtLa(0)
 512   ,fh2MCEtaVsPtALa(0)
 513   ,fh1MCRapK0s(0)
 514   ,fh1MCRapLambda(0)
 515   ,fh1MCRapAntiLambda(0)
 516   ,fh1MCEtaAllK0s(0)
 517   ,fh1MCEtaK0s(0)
 518   ,fh1MCEtaLambda(0)
 519   ,fh1MCEtaAntiLambda(0)
 520
 521
 522 {
 523   // constructor
 524
 525   DefineOutput(1,TList::Class());
 526 }
 527
 528 //__________________________________________________________________________________________________________________________
 529 AliAnalysisTaskJetChem::AliAnalysisTaskJetChem(const  AliAnalysisTaskJetChem &copy)
 530   : AliAnalysisTaskFragmentationFunction()
 531
 532   ,fAnalysisMC(copy.fAnalysisMC)
 533   ,fDeltaVertexZ(copy.fDeltaVertexZ)
 534   ,fCuttrackNegNcls(copy.fCuttrackNegNcls)
 535   ,fCuttrackPosNcls(copy.fCuttrackPosNcls)
 536   ,fCutPostrackRap(copy.fCutPostrackRap)
 537   ,fCutNegtrackRap(copy.fCutNegtrackRap)
 538   ,fCutRap(copy.fCutRap)
 539   ,fCutPostrackEta(copy.fCutPostrackEta)
 540   ,fCutNegtrackEta(copy.fCutNegtrackEta)
 541   ,fCutEta(copy.fCutEta)
 542   ,fCutV0cosPointAngle(copy.fCutV0cosPointAngle)
 543   ,fCutChi2PosDaughter(copy.fCutChi2PosDaughter)
 544   ,fCutChi2NegDaughter(copy.fCutChi2NegDaughter)
 545   ,fKinkDaughters(copy.fKinkDaughters)
 546   ,fRequireTPCRefit(copy.fRequireTPCRefit)
 547   ,fCutArmenteros(copy.fCutArmenteros)
 548   ,fCutV0DecayMin(copy.fCutV0DecayMin)
 549   ,fCutV0DecayMax(copy.fCutV0DecayMax)
 550   ,fCutV0totMom(copy.fCutV0totMom)
 551   ,fCutDcaV0Daughters(copy.fCutDcaV0Daughters)
 552   ,fCutDcaPosToPrimVertex(copy.fCutDcaPosToPrimVertex)
 553   ,fCutDcaNegToPrimVertex(copy.fCutDcaNegToPrimVertex)
 554   ,fCutV0RadiusMin(copy.fCutV0RadiusMin)
 555   ,fCutV0RadiusMax(copy.fCutV0RadiusMax)
 556   ,fCutBetheBloch(copy.fCutBetheBloch)
 557   ,fCutRatio(copy.fCutRatio)
 558   ,fK0Type(copy.fK0Type)
 559   ,fFilterMaskK0(copy.fFilterMaskK0)
 560   ,fListK0s(copy.fListK0s)
 561   ,fPIDResponse(copy.fPIDResponse)
 562   ,fV0QAK0(copy.fV0QAK0)
 563   ,fFFHistosRecCutsK0Evt(copy.fFFHistosRecCutsK0Evt)
 564   ,fFFHistosIMK0AllEvt(copy.fFFHistosIMK0AllEvt)
 565   ,fFFHistosIMK0Jet(copy.fFFHistosIMK0Jet)
 566   ,fFFHistosIMK0Cone(copy.fFFHistosIMK0Cone)
 567   ,fFFHistosPhiCorrIMK0(copy.fFFHistosPhiCorrIMK0)
 568   ,fLaType(copy.fLaType)
 569   ,fFilterMaskLa(copy.fFilterMaskLa)
 570   ,fListLa(copy.fListLa)
 571   ,fFFHistosIMLaAllEvt(copy.fFFHistosIMLaAllEvt)
 572   ,fFFHistosIMLaJet(copy.fFFHistosIMLaJet)
 573   ,fFFHistosIMLaCone(copy.fFFHistosIMLaCone)
 574   ,fFFHistosPhiCorrIMLa(copy.fFFHistosPhiCorrIMLa)
 575   ,fALaType(copy.fALaType)
 576   ,fFilterMaskALa(copy.fFilterMaskALa)
 577   ,fListALa(copy.fListALa)
 578   ,fListFeeddownLaCand(copy.fListFeeddownLaCand)
 579   ,fListFeeddownALaCand(copy.fListFeeddownALaCand)
 580   ,jetConeFDLalist(copy.jetConeFDLalist)
 581   ,jetConeFDALalist(copy.jetConeFDALalist)
 582   ,fListMCgenK0s(copy.fListMCgenK0s)
 583   ,fListMCgenLa(copy.fListMCgenLa)
 584   ,fListMCgenALa(copy.fListMCgenALa)
 585   ,fListMCgenK0sCone(copy.fListMCgenK0sCone)
 586   ,fListMCgenLaCone(copy.fListMCgenLaCone)
 587   ,fListMCgenALaCone(copy.fListMCgenALaCone)
 588   ,IsArmenterosSelected(copy.IsArmenterosSelected)
 589   ,fFFHistosIMALaAllEvt(copy.fFFHistosIMALaAllEvt)
 590   ,fFFHistosIMALaJet(copy.fFFHistosIMALaJet)
 591   ,fFFHistosIMALaCone(copy.fFFHistosIMALaCone)
 592   ,fFFHistosPhiCorrIMALa(copy.fFFHistosPhiCorrIMALa)
 593   ,fFFIMNBinsJetPt(copy.fFFIMNBinsJetPt)
 594   ,fFFIMJetPtMin(copy.fFFIMJetPtMin)
 595   ,fFFIMJetPtMax(copy.fFFIMJetPtMax)
 596   ,fFFIMNBinsInvM(copy.fFFIMNBinsInvM)
 597   ,fFFIMInvMMin(copy.fFFIMInvMMin)
 598   ,fFFIMInvMMax(copy.fFFIMInvMMax)
 599   ,fFFIMNBinsPt(copy.fFFIMNBinsPt)
 600   ,fFFIMPtMin(copy.fFFIMPtMin)
 601   ,fFFIMPtMax(copy.fFFIMPtMax)
 602   ,fFFIMNBinsXi(copy.fFFIMNBinsXi)
 603   ,fFFIMXiMin(copy.fFFIMXiMin)
 604   ,fFFIMXiMax(copy.fFFIMXiMax)
 605   ,fFFIMNBinsZ(copy.fFFIMNBinsZ)
 606   ,fFFIMZMin(copy.fFFIMZMin)
 607   ,fFFIMZMax(copy.fFFIMZMax)
 608   ,fFFIMLaNBinsJetPt(copy.fFFIMLaNBinsJetPt)
 609   ,fFFIMLaJetPtMin(copy.fFFIMLaJetPtMin)
 610   ,fFFIMLaJetPtMax(copy.fFFIMLaJetPtMax)
 611   ,fFFIMLaNBinsInvM(copy.fFFIMLaNBinsInvM)
 612   ,fFFIMLaInvMMin(copy.fFFIMLaInvMMin)
 613   ,fFFIMLaInvMMax(copy.fFFIMLaInvMMax)
 614   ,fFFIMLaNBinsPt(copy.fFFIMLaNBinsPt)
 615   ,fFFIMLaPtMin(copy.fFFIMLaPtMin)
 616   ,fFFIMLaPtMax(copy.fFFIMLaPtMax)
 617   ,fFFIMLaNBinsXi(copy.fFFIMLaNBinsXi)
 618   ,fFFIMLaXiMin(copy.fFFIMLaXiMin)
 619   ,fFFIMLaXiMax(copy.fFFIMLaXiMax)
 620   ,fFFIMLaNBinsZ(copy.fFFIMLaNBinsZ)
 621   ,fFFIMLaZMin(copy.fFFIMLaZMin)
 622   ,fFFIMLaZMax(copy.fFFIMLaZMax)
 623   ,fPhiCorrIMNBinsPt(copy.fPhiCorrIMNBinsPt)
 624   ,fPhiCorrIMPtMin(copy.fPhiCorrIMPtMin)
 625   ,fPhiCorrIMPtMax(copy.fPhiCorrIMPtMax)
 626   ,fPhiCorrIMNBinsPhi(copy.fPhiCorrIMNBinsPhi)
 627   ,fPhiCorrIMPhiMin(copy.fPhiCorrIMPhiMin)
 628   ,fPhiCorrIMPhiMax(copy.fPhiCorrIMPhiMax)
 629   ,fPhiCorrIMNBinsInvM(copy.fPhiCorrIMNBinsInvM)
 630   ,fPhiCorrIMInvMMin(copy.fPhiCorrIMInvMMin)
 631   ,fPhiCorrIMInvMMax(copy.fPhiCorrIMInvMMax)
 632   ,fPhiCorrIMLaNBinsPt(copy.fPhiCorrIMLaNBinsPt)
 633   ,fPhiCorrIMLaPtMin(copy.fPhiCorrIMLaPtMin)
 634   ,fPhiCorrIMLaPtMax(copy.fPhiCorrIMLaPtMax)
 635   ,fPhiCorrIMLaNBinsPhi(copy.fPhiCorrIMLaNBinsPhi)
 636   ,fPhiCorrIMLaPhiMin(copy.fPhiCorrIMLaPhiMin)
 637   ,fPhiCorrIMLaPhiMax(copy.fPhiCorrIMLaPhiMax)
 638   ,fPhiCorrIMLaNBinsInvM(copy.fPhiCorrIMLaNBinsInvM)
 639   ,fPhiCorrIMLaInvMMin(copy.fPhiCorrIMLaInvMMin)
 640   ,fPhiCorrIMLaInvMMax(copy.fPhiCorrIMLaInvMMax)
 641   ,fh1EvtAllCent(copy.fh1EvtAllCent)
 642   ,fh1Evt(copy.fh1Evt)
 643   ,fh1K0Mult(copy.fh1K0Mult)
 644   ,fh1dPhiJetK0(copy.fh1dPhiJetK0)
 645   ,fh1LaMult(copy.fh1LaMult)
 646   ,fh1dPhiJetLa(copy.fh1dPhiJetLa)
 647   ,fh1ALaMult(copy.fh1ALaMult)
 648   ,fh1dPhiJetALa(copy.fh1dPhiJetALa)
 649   ,fh1JetEta(copy.fh1JetEta)
 650   ,fh1JetPhi(copy.fh1JetPhi)
 651   ,fh2JetEtaPhi(copy.fh2JetEtaPhi)
 652   ,fh1V0JetPt(copy.fh1V0JetPt)
 653   ,fh2FFJetTrackEta(copy.fh2FFJetTrackEta)
 654   ,fh1trackPosNCls(copy.fh1trackPosNCls)
 655   ,fh1trackNegNCls(copy.fh1trackNegNCls)
 656   ,fh1trackPosRap(copy.fh1trackPosRap)
 657   ,fh1trackNegRap(copy.fh1trackNegRap)
 658   ,fh1V0Rap(copy.fh1V0Rap)
 659   ,fh1trackPosEta(copy.fh1trackPosEta)
 660   ,fh1trackNegEta(copy.fh1trackNegEta)
 661   ,fh1V0Eta(copy.fh1V0Eta)
 662   ,fh1V0totMom(copy.fh1V0totMom)
 663   ,fh1CosPointAngle(copy.fh1CosPointAngle)
 664   ,fh1Chi2Pos(copy.fh1Chi2Pos)
 665   ,fh1Chi2Neg(copy.fh1Chi2Neg)
 666   ,fh1DecayLengthV0(copy.fh1DecayLengthV0)
 667   ,fh2ProperLifetimeK0sVsPtBeforeCut(copy.fh2ProperLifetimeK0sVsPtBeforeCut)
 668   ,fh2ProperLifetimeK0sVsPtAfterCut(copy.fh2ProperLifetimeK0sVsPtAfterCut)
 669   ,fh1ProperLifetimeV0BeforeCut(copy.fh1ProperLifetimeV0BeforeCut)
 670   ,fh1ProperLifetimeV0AfterCut(copy.fh1ProperLifetimeV0AfterCut)
 671   ,fh1V0Radius(copy.fh1V0Radius)
 672   ,fh1DcaV0Daughters(copy.fh1DcaV0Daughters)
 673   ,fh1DcaPosToPrimVertex(copy.fh1DcaPosToPrimVertex)
 674   ,fh1DcaNegToPrimVertex(copy.fh1DcaNegToPrimVertex)
 675   ,fh2ArmenterosBeforeCuts(copy.fh2ArmenterosBeforeCuts)
 676   ,fh2ArmenterosAfterCuts(copy.fh2ArmenterosAfterCuts)
 677   ,fh2BBLaPos(copy.fh2BBLaPos)
 678   ,fh2BBLaNeg(copy.fh2BBLaPos)
 679   ,fh1CrossedRowsOverFindableNeg(copy.fh1CrossedRowsOverFindableNeg)
 680   ,fh1CrossedRowsOverFindablePos(copy.fh1CrossedRowsOverFindablePos)
 681   ,fh1PosDaughterCharge(copy.fh1PosDaughterCharge)
 682   ,fh1NegDaughterCharge(copy.fh1NegDaughterCharge)
 683   ,fh1PtMCK0s(copy.fh1PtMCK0s)
 684   ,fh1PtMCLa(copy.fh1PtMCLa)
 685   ,fh1PtMCALa(copy.fh1PtMCALa)
 686   ,fh1EtaK0s(copy.fh1EtaK0s)
 687   ,fh1EtaLa(copy.fh1EtaLa)
 688   ,fh1EtaALa(copy.fh1EtaALa)
 689   ,fh3InvMassEtaTrackPtK0s(copy.fh3InvMassEtaTrackPtK0s)
 690   ,fh3InvMassEtaTrackPtLa(copy.fh3InvMassEtaTrackPtLa)
 691   ,fh3InvMassEtaTrackPtALa(copy.fh3InvMassEtaTrackPtALa)
 692
 693   ,fh1TrackMultCone(copy.fh1TrackMultCone)
 694   ,fh2TrackMultCone(copy.fh2TrackMultCone)
 695   ,fh2MCgenK0Cone(copy.fh2MCgenK0Cone)
 696   ,fh2MCgenLaCone(copy.fh2MCgenLaCone)
 697   ,fh2MCgenALaCone(copy.fh2MCgenALaCone)
 698   ,fh2MCEtagenK0Cone(copy.fh2MCEtagenK0Cone)
 699   ,fh2MCEtagenLaCone(copy.fh2MCEtagenLaCone)
 700   ,fh2MCEtagenALaCone(copy.fh2MCEtagenALaCone)
 701   ,fh1FFIMK0ConeSmear(copy.fh1FFIMK0ConeSmear)
 702   ,fh1FFIMLaConeSmear(copy.fh1FFIMLaConeSmear)
 703   ,fh1FFIMALaConeSmear(copy.fh1FFIMALaConeSmear)
 704   ,fh3MCrecK0Cone(copy.fh3MCrecK0Cone)
 705   ,fh3MCrecLaCone(copy.fh3MCrecLaCone)
 706   ,fh3MCrecALaCone(copy.fh3MCrecALaCone)
 707   ,fh3MCrecK0ConeSmear(copy.fh3MCrecK0ConeSmear)
 708   ,fh3MCrecLaConeSmear(copy.fh3MCrecLaConeSmear)
 709   ,fh3MCrecALaConeSmear(copy.fh3MCrecALaConeSmear)
 710   ,fh3SecContinCone(copy.fh3SecContinCone)
 711   ,fh3StrContinCone(copy.fh3StrContinCone)
 712   ,fh3IMK0PerpCone(copy.fh3IMK0PerpCone)
 713   ,fh3IMLaPerpCone(copy.fh3IMLaPerpCone)
 714   ,fh3IMALaPerpCone(copy.fh3IMALaPerpCone)
 715   ,fh3IMK0MedianCone(copy.fh3IMK0MedianCone)
 716   ,fh3IMLaMedianCone(copy.fh3IMLaMedianCone)
 717   ,fh3IMALaMedianCone(copy.fh3IMALaMedianCone)
 718   ,fh1MCMultiplicityPrimary(copy.fh1MCMultiplicityPrimary)
 719   ,fh1MCMultiplicityTracks(copy.fh1MCMultiplicityTracks)
 720   ,fh1MCmotherLa(copy.fh1MCmotherLa)
 721   ,fh1MCmotherALa(copy.fh1MCmotherALa)
 722   ,fh3FeedDownLa(copy.fh3FeedDownLa)
 723   ,fh3FeedDownALa(copy.fh3FeedDownALa)
 724   ,fh3FeedDownLaCone(copy.fh3FeedDownLaCone)
 725   ,fh3FeedDownALaCone(copy.fh3FeedDownALaCone)
 726   ,fh1MCProdRadiusK0s(copy.fh1MCProdRadiusK0s)
 727   ,fh1MCProdRadiusLambda(copy.fh1MCProdRadiusLambda)
 728   ,fh1MCProdRadiusAntiLambda(copy.fh1MCProdRadiusAntiLambda)
 729   ,fh1MCPtV0s(copy.fh1MCPtV0s)
 730   ,fh1MCPtK0s(copy.fh1MCPtK0s)
 731   ,fh1MCPtLambda(copy.fh1MCPtLambda)
 732   ,fh1MCPtAntiLambda(copy.fh1MCPtAntiLambda)
 733   ,fh1MCXiPt(copy.fh1MCXiPt)
 734   ,fh1MCXibarPt(copy.fh1MCXibarPt)
 735   ,fh2MCEtaVsPtK0s(copy.fh2MCEtaVsPtK0s)
 736   ,fh2MCEtaVsPtLa(copy.fh2MCEtaVsPtLa)
 737   ,fh2MCEtaVsPtALa(copy.fh2MCEtaVsPtALa)
 738   ,fh1MCRapK0s(copy.fh1MCRapK0s)
 739   ,fh1MCRapLambda(copy.fh1MCRapLambda)
 740   ,fh1MCRapAntiLambda(copy.fh1MCRapAntiLambda)
 741   ,fh1MCEtaAllK0s(copy.fh1MCEtaAllK0s)
 742   ,fh1MCEtaK0s(copy.fh1MCEtaK0s)
 743   ,fh1MCEtaLambda(copy.fh1MCEtaLambda)
 744   ,fh1MCEtaAntiLambda(copy.fh1MCEtaAntiLambda)
 745
 746 {
 747   // copy constructor
 748
 749 }
 750
 751 // _________________________________________________________________________________________________________________________________
 752 AliAnalysisTaskJetChem& AliAnalysisTaskJetChem::operator=(const AliAnalysisTaskJetChem& o)
 753 {
 754   // assignment
 755
 756   if(this!=&o){
 757     AliAnalysisTaskFragmentationFunction::operator=(o);
 758
 759     fAnalysisMC                     = o.fAnalysisMC;
 760     fDeltaVertexZ                   = o.fDeltaVertexZ;
 761     fCuttrackNegNcls                = o.fCuttrackNegNcls;
 762     fCuttrackPosNcls                = o.fCuttrackPosNcls;
 763     fCutPostrackRap                 = o.fCutPostrackRap;
 764     fCutNegtrackRap                 = o.fCutNegtrackRap;
 765     fCutRap                         = o.fCutRap;
 766     fCutPostrackEta                 = o.fCutPostrackEta;
 767     fCutNegtrackEta                 = o.fCutNegtrackEta;
 768     fCutEta                         = o.fCutEta;
 769     fCutV0cosPointAngle             = o.fCutV0cosPointAngle;
 770     fCutChi2PosDaughter             = o.fCutChi2PosDaughter;
 771     fCutChi2NegDaughter             = o.fCutChi2NegDaughter;
 772     fKinkDaughters                  = o.fKinkDaughters;
 773     fRequireTPCRefit                = o.fRequireTPCRefit;
 774     fCutArmenteros                  = o.fCutArmenteros;
 775     fCutV0DecayMin                  = o.fCutV0DecayMin;
 776     fCutV0DecayMax                  = o.fCutV0DecayMax;
 777     fCutV0totMom                    = o.fCutV0totMom;
 778     fCutDcaV0Daughters              = o.fCutDcaV0Daughters;
 779     fCutDcaPosToPrimVertex          = o.fCutDcaPosToPrimVertex;
 780     fCutDcaNegToPrimVertex          = o.fCutDcaNegToPrimVertex;
 781     fCutV0RadiusMin                 = o.fCutV0RadiusMin;
 782     fCutV0RadiusMax                 = o.fCutV0RadiusMax;
 783     fCutBetheBloch                  = o.fCutBetheBloch;
 784     fCutRatio                       = o.fCutRatio;
 785     fK0Type                         = o.fK0Type;
 786     fFilterMaskK0                   = o.fFilterMaskK0;
 787     fListK0s                        = o.fListK0s;
 788     fPIDResponse                    = o.fPIDResponse;
 789     fV0QAK0                         = o.fV0QAK0;
 790     fFFHistosRecCutsK0Evt           = o.fFFHistosRecCutsK0Evt;
 791     fFFHistosIMK0AllEvt             = o.fFFHistosIMK0AllEvt;
 792     fFFHistosIMK0Jet                = o.fFFHistosIMK0Jet;
 793     fFFHistosIMK0Cone               = o.fFFHistosIMK0Cone;
 794     fFFHistosPhiCorrIMK0            = o.fFFHistosPhiCorrIMK0;
 795     fLaType                         = o.fLaType;
 796     fFilterMaskLa                   = o.fFilterMaskLa;
 797     fListLa                         = o.fListLa;
 798     fFFHistosIMLaAllEvt             = o.fFFHistosIMLaAllEvt;
 799     fFFHistosIMLaJet                = o.fFFHistosIMLaJet;
 800     fFFHistosIMLaCone               = o.fFFHistosIMLaCone;
 801     fFFHistosPhiCorrIMLa            = o.fFFHistosPhiCorrIMLa;
 802     fALaType                        = o.fALaType;
 803     fFilterMaskALa                  = o.fFilterMaskALa;
 804     fListFeeddownLaCand             = o.fListFeeddownLaCand;
 805     fListFeeddownALaCand            = o.fListFeeddownALaCand;
 806     jetConeFDLalist                 = o.jetConeFDLalist;
 807     jetConeFDALalist                = o.jetConeFDALalist;
 808     fListMCgenK0s                   = o.fListMCgenK0s;
 809     fListMCgenLa                    = o.fListMCgenLa;
 810     fListMCgenALa                   = o.fListMCgenALa;
 811     fListMCgenK0sCone               = o.fListMCgenK0sCone;
 812     fListMCgenLaCone                = o.fListMCgenLaCone;
 813     fListMCgenALaCone               = o.fListMCgenALaCone;
 814     IsArmenterosSelected            = o.IsArmenterosSelected;
 815     fFFHistosIMALaAllEvt            = o.fFFHistosIMALaAllEvt;
 816     fFFHistosIMALaJet               = o.fFFHistosIMALaJet;
 817     fFFHistosIMALaCone              = o.fFFHistosIMALaCone;
 818     fFFHistosPhiCorrIMALa           = o.fFFHistosPhiCorrIMALa;
 819     fFFIMNBinsJetPt                 = o.fFFIMNBinsJetPt;
 820     fFFIMJetPtMin                   = o.fFFIMJetPtMin;
 821     fFFIMJetPtMax                   = o.fFFIMJetPtMax;
 822     fFFIMNBinsPt                    = o.fFFIMNBinsPt;
 823     fFFIMPtMin                      = o.fFFIMPtMin;
 824     fFFIMPtMax                      = o.fFFIMPtMax;
 825     fFFIMNBinsXi                    = o.fFFIMNBinsXi;
 826     fFFIMXiMin                      = o.fFFIMXiMin;
 827     fFFIMXiMax                      = o.fFFIMXiMax;
 828     fFFIMNBinsZ                     = o.fFFIMNBinsZ;
 829     fFFIMZMin                       = o.fFFIMZMin;
 830     fFFIMZMax                       = o.fFFIMZMax;
 831     fFFIMLaNBinsJetPt               = o.fFFIMLaNBinsJetPt;
 832     fFFIMLaJetPtMin                 = o.fFFIMLaJetPtMin;
 833     fFFIMLaJetPtMax                 = o.fFFIMLaJetPtMax;
 834     fFFIMLaNBinsPt                  = o.fFFIMLaNBinsPt;
 835     fFFIMLaPtMin                    = o.fFFIMLaPtMin;
 836     fFFIMLaPtMax                    = o.fFFIMLaPtMax;
 837     fFFIMLaNBinsXi                  = o.fFFIMLaNBinsXi;
 838     fFFIMLaXiMin                    = o.fFFIMLaXiMin;
 839     fFFIMLaXiMax                    = o.fFFIMLaXiMax;
 840     fFFIMLaNBinsZ                   = o.fFFIMLaNBinsZ;
 841     fFFIMLaZMin                     = o.fFFIMLaZMin;
 842     fFFIMLaZMax                     = o.fFFIMLaZMax;
 843     fPhiCorrIMNBinsPt               = o.fPhiCorrIMNBinsPt;
 844     fPhiCorrIMPtMin                 = o.fPhiCorrIMPtMin;
 845     fPhiCorrIMPtMax                 = o.fPhiCorrIMPtMax;
 846     fPhiCorrIMNBinsPhi              = o.fPhiCorrIMNBinsPhi;
 847     fPhiCorrIMPhiMin                = o.fPhiCorrIMPhiMin;
 848     fPhiCorrIMPhiMax                = o.fPhiCorrIMPhiMax;
 849     fPhiCorrIMNBinsInvM             = o.fPhiCorrIMNBinsInvM;
 850     fPhiCorrIMInvMMin               = o.fPhiCorrIMInvMMin;
 851     fPhiCorrIMInvMMax               = o.fPhiCorrIMInvMMax;
 852     fPhiCorrIMLaNBinsPt             = o.fPhiCorrIMLaNBinsPt;
 853     fPhiCorrIMLaPtMin               = o.fPhiCorrIMLaPtMin;
 854     fPhiCorrIMLaPtMax               = o.fPhiCorrIMLaPtMax;
 855     fPhiCorrIMLaNBinsPhi            = o.fPhiCorrIMLaNBinsPhi;
 856     fPhiCorrIMLaPhiMin              = o.fPhiCorrIMLaPhiMin;
 857     fPhiCorrIMLaPhiMax              = o.fPhiCorrIMLaPhiMax;
 858     fPhiCorrIMLaNBinsInvM           = o.fPhiCorrIMLaNBinsInvM;
 859     fPhiCorrIMLaInvMMin             = o.fPhiCorrIMLaInvMMin;
 860     fPhiCorrIMLaInvMMax             = o.fPhiCorrIMLaInvMMax;
 861     fh1EvtAllCent                   = o.fh1EvtAllCent;
 862     fh1Evt                          = o.fh1Evt;
 863     fh1K0Mult                       = o.fh1K0Mult;
 864     fh1dPhiJetK0                    = o.fh1dPhiJetK0;
 865     fh1LaMult                       = o.fh1LaMult;
 866     fh1dPhiJetLa                    = o.fh1dPhiJetLa;
 867     fh1ALaMult                      = o.fh1ALaMult;
 868     fh1dPhiJetALa                   = o.fh1dPhiJetALa;
 869     fh1JetEta                       = o.fh1JetEta;
 870     fh1JetPhi                       = o.fh1JetPhi;
 871     fh2JetEtaPhi                    = o.fh2JetEtaPhi;
 872     fh1V0JetPt                     = o.fh1V0JetPt;
 873     fh2FFJetTrackEta                = o.fh2FFJetTrackEta;
 874     fh1trackPosNCls                 = o.fh1trackPosNCls;
 875     fh1trackNegNCls                 = o.fh1trackNegNCls;
 876     fh1trackPosRap                  = o.fh1trackPosRap;
 877     fh1trackNegRap                  = o.fh1trackNegRap;
 878     fh1V0Rap                        = o.fh1V0Rap;
 879     fh1trackPosEta                  = o.fh1trackPosEta;
 880     fh1trackNegEta                  = o.fh1trackNegEta;
 881     fh1V0Eta                        = o.fh1V0Eta;
 882     fh1V0totMom                     = o.fh1V0totMom;
 883     fh1CosPointAngle                = o.fh1CosPointAngle;
 884     fh1Chi2Pos                      = o.fh1Chi2Pos;
 885     fh1Chi2Neg                      = o.fh1Chi2Neg;
 886     fh1DecayLengthV0                = o.fh1DecayLengthV0;
 887     fh2ProperLifetimeK0sVsPtBeforeCut = o.fh2ProperLifetimeK0sVsPtBeforeCut;
 888     fh2ProperLifetimeK0sVsPtAfterCut= o.fh2ProperLifetimeK0sVsPtAfterCut;
 889     fh1ProperLifetimeV0BeforeCut    = o.fh1ProperLifetimeV0BeforeCut;
 890     fh1ProperLifetimeV0AfterCut     = o.fh1ProperLifetimeV0AfterCut;
 891     fh1V0Radius                     = o.fh1V0Radius;
 892     fh1DcaV0Daughters               = o.fh1DcaV0Daughters;
 893     fh1DcaPosToPrimVertex           = o.fh1DcaPosToPrimVertex;
 894     fh1DcaNegToPrimVertex           = o.fh1DcaNegToPrimVertex;
 895     fh2ArmenterosBeforeCuts         = o.fh2ArmenterosBeforeCuts;
 896     fh2ArmenterosAfterCuts          = o.fh2ArmenterosAfterCuts;
 897     fh2BBLaPos                      = o.fh2BBLaPos;
 898     fh2BBLaNeg                      = o.fh2BBLaPos;
 899     fh1CrossedRowsOverFindableNeg   = o.fh1CrossedRowsOverFindableNeg;
 900     fh1CrossedRowsOverFindablePos   = o.fh1CrossedRowsOverFindablePos;
 901     fh1PosDaughterCharge            = o.fh1PosDaughterCharge;
 902     fh1NegDaughterCharge            = o.fh1NegDaughterCharge;
 903     fh1PtMCK0s                      = o.fh1PtMCK0s;
 904     fh1PtMCLa                       = o.fh1PtMCLa;
 905     fh1PtMCALa                      = o.fh1PtMCALa;
 906     fh1EtaK0s                       = o.fh1EtaK0s;
 907     fh1EtaLa                        = o.fh1EtaLa;
 908     fh1EtaALa                       = o.fh1EtaALa;
 909     fh3InvMassEtaTrackPtK0s         = o.fh3InvMassEtaTrackPtK0s;
 910     fh3InvMassEtaTrackPtLa          = o.fh3InvMassEtaTrackPtLa;
 911     fh3InvMassEtaTrackPtALa         = o.fh3InvMassEtaTrackPtALa;
 912     fh1TrackMultCone                = o.fh1TrackMultCone;
 913     fh2TrackMultCone                = o.fh2TrackMultCone;
 914     fh2MCgenK0Cone                  = o.fh2MCgenK0Cone;
 915     fh2MCgenLaCone                  = o.fh2MCgenLaCone;
 916     fh2MCgenALaCone                 = o.fh2MCgenALaCone;
 917     fh2MCEtagenK0Cone               = o.fh2MCEtagenK0Cone;
 918     fh2MCEtagenLaCone               = o.fh2MCEtagenLaCone;
 919     fh2MCEtagenALaCone              = o.fh2MCEtagenALaCone;
 920     fh1FFIMK0ConeSmear              = o.fh1FFIMK0ConeSmear;
 921     fh1FFIMLaConeSmear              = o.fh1FFIMLaConeSmear;
 922     fh1FFIMALaConeSmear             = o.fh1FFIMALaConeSmear;
 923     fh3MCrecK0Cone                  = o.fh3MCrecK0Cone;
 924     fh3MCrecLaCone                  = o.fh3MCrecLaCone;
 925     fh3MCrecALaCone                 = o.fh3MCrecALaCone;
 926     fh3MCrecK0ConeSmear             = o.fh3MCrecK0ConeSmear;
 927     fh3MCrecLaConeSmear             = o.fh3MCrecLaConeSmear;
 928     fh3MCrecALaConeSmear            = o.fh3MCrecALaConeSmear;
 929     fh3SecContinCone                = o.fh3SecContinCone;
 930     fh3StrContinCone                = o.fh3StrContinCone;
 931     fh3IMK0PerpCone                 = o.fh3IMK0PerpCone;
 932     fh3IMLaPerpCone                 = o.fh3IMLaPerpCone;
 933     fh3IMALaPerpCone                = o.fh3IMALaPerpCone;
 934     fh3IMK0MedianCone               = o.fh3IMK0MedianCone;
 935     fh3IMLaMedianCone               = o.fh3IMLaMedianCone;
 936     fh3IMALaMedianCone              = o.fh3IMALaMedianCone;
 937     fh1MCMultiplicityPrimary        = o.fh1MCMultiplicityPrimary;
 938     fh1MCMultiplicityTracks         = o.fh1MCMultiplicityTracks;
 939     fh1MCmotherLa                   = o.fh1MCmotherLa;
 940     fh1MCmotherALa                  = o.fh1MCmotherALa;
 941     fh3FeedDownLa                   = o.fh3FeedDownLa;
 942     fh3FeedDownALa                  = o.fh3FeedDownALa;
 943     fh3FeedDownLaCone               = o.fh3FeedDownLaCone;
 944     fh3FeedDownALaCone              = o.fh3FeedDownALaCone;
 945     fh1MCProdRadiusK0s              = o.fh1MCProdRadiusK0s;
 946     fh1MCProdRadiusLambda           = o.fh1MCProdRadiusLambda;
 947     fh1MCProdRadiusAntiLambda       = o.fh1MCProdRadiusAntiLambda;
 948     fh1MCPtV0s                      = o.fh1MCPtV0s;
 949     fh1MCPtK0s                      = o.fh1MCPtK0s;
 950     fh1MCPtLambda                   = o.fh1MCPtLambda;
 951     fh1MCPtAntiLambda               = o.fh1MCPtAntiLambda;
 952     fh1MCXiPt                       = o.fh1MCXiPt;
 953     fh1MCXibarPt                    = o.fh1MCXibarPt;
 954     fh2MCEtaVsPtK0s                 = o.fh2MCEtaVsPtK0s;
 955     fh2MCEtaVsPtLa                  = o.fh2MCEtaVsPtLa;
 956     fh2MCEtaVsPtALa                 = o.fh2MCEtaVsPtALa;
 957     fh1MCRapK0s                     = o.fh1MCRapK0s;
 958     fh1MCRapLambda                  = o.fh1MCRapLambda;
 959     fh1MCRapAntiLambda              = o.fh1MCRapAntiLambda;
 960     fh1MCEtaAllK0s                  = o.fh1MCEtaAllK0s;
 961     fh1MCEtaK0s                     = o.fh1MCEtaK0s;
 962     fh1MCEtaLambda                  = o.fh1MCEtaLambda;
 963     fh1MCEtaAntiLambda              = o.fh1MCEtaAntiLambda;
 964 }
 965
 966   return *this;
 967 }
 968
 969 //_______________________________________________
 970 AliAnalysisTaskJetChem::~AliAnalysisTaskJetChem()
 971 {
 972   // destructor
 973
 974
 975   if(fListK0s) delete fListK0s;
 976   if(fListLa) delete fListLa;
 977   if(fListALa) delete fListALa;
 978   if(fListFeeddownLaCand) delete fListFeeddownLaCand;
 979   if(fListFeeddownALaCand) delete fListFeeddownALaCand;
 980   if(jetConeFDLalist) delete jetConeFDLalist;
 981   if(jetConeFDALalist) delete jetConeFDALalist;
 982   if(fListMCgenK0s) delete fListMCgenK0s;
 983   if(fListMCgenLa) delete fListMCgenLa;
 984   if(fListMCgenALa) delete fListMCgenALa;
 985
 986
 987
 988 }
 989
 990 //________________________________________________________________________________________________________________________________
 991 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AliFragFuncHistosInvMass(const char* name,
 992                                                                            Int_t nJetPt, Float_t jetPtMin, Float_t jetPtMax,
 993                                                                            Int_t nInvMass, Float_t invMassMin, Float_t invMassMax,
 994                                                                            Int_t nPt, Float_t ptMin, Float_t ptMax,
 995                                                                            Int_t nXi, Float_t xiMin, Float_t xiMax,
 996                                                                            Int_t nZ , Float_t zMin , Float_t zMax )
 997   : TObject()
 998   ,fNBinsJetPt(nJetPt)
 999   ,fJetPtMin(jetPtMin)
1000   ,fJetPtMax(jetPtMax)
1001   ,fNBinsInvMass(nInvMass)
1002   ,fInvMassMin(invMassMin)
1003   ,fInvMassMax(invMassMax)
1004   ,fNBinsPt(nPt)
1005   ,fPtMin(ptMin)
1006   ,fPtMax(ptMax)
1007   ,fNBinsXi(nXi)
1008   ,fXiMin(xiMin)
1009   ,fXiMax(xiMax)
1010   ,fNBinsZ(nZ)
1011   ,fZMin(zMin)
1012   ,fZMax(zMax)
1013   ,fh3TrackPt(0)
1014   ,fh3Xi(0)
1015   ,fh3Z(0)
1016   ,fh1JetPt(0)
1017   ,fNameFF(name)
1018 {
1019   // default constructor
1020
1021 }
1022
1023 //______________________________________________________________________________________________________________
1024 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AliFragFuncHistosInvMass(const AliFragFuncHistosInvMass& copy)
1025   : TObject()
1026   ,fNBinsJetPt(copy.fNBinsJetPt)
1027   ,fJetPtMin(copy.fJetPtMin)
1028   ,fJetPtMax(copy.fJetPtMax)
1029   ,fNBinsInvMass(copy.fNBinsInvMass)
1030   ,fInvMassMin(copy.fInvMassMin)
1031   ,fInvMassMax(copy.fInvMassMax)
1032   ,fNBinsPt(copy.fNBinsPt)
1033   ,fPtMin(copy.fPtMin)
1034
1035   ,fPtMax(copy.fPtMax)
1036   ,fNBinsXi(copy.fNBinsXi)
1037   ,fXiMin(copy.fXiMin)
1038   ,fXiMax(copy.fXiMax)
1039   ,fNBinsZ(copy.fNBinsZ)
1040   ,fZMin(copy.fZMin)
1041   ,fZMax(copy.fZMax)
1042   ,fh3TrackPt(copy.fh3TrackPt)
1043   ,fh3Xi(copy.fh3Xi)
1044   ,fh3Z(copy.fh3Z)
1045   ,fh1JetPt(copy.fh1JetPt)
1046   ,fNameFF(copy.fNameFF)
1047 {
1048   // copy constructor
1049 }
1050
1051 //______________________________________________________________________________________________________________________________________________________________________
1052 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass& AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::operator=(const AliAnalysisTaskJetChem::AliFragFuncHistosInvMass& o)
1053 {
1054   // assignment
1055
1056   if(this!=&o){
1057     TObject::operator=(o);
1058     fNBinsJetPt   = o.fNBinsJetPt;
1059     fJetPtMin     = o.fJetPtMin;
1060     fJetPtMax     = o.fJetPtMax;
1061     fNBinsInvMass = o.fNBinsInvMass;
1062     fInvMassMin   = o.fInvMassMin;
1063     fInvMassMax   = o.fInvMassMax;
1064     fNBinsPt      = o.fNBinsPt;
1065     fPtMin        = o.fPtMin;
1066     fPtMax        = o.fPtMax;
1067     fNBinsXi      = o.fNBinsXi;
1068     fXiMin        = o.fXiMin;
1069     fXiMax        = o.fXiMax;
1070     fNBinsZ       = o.fNBinsZ;
1071     fZMin         = o.fZMin;
1072     fZMax         = o.fZMax;
1073     fh3TrackPt    = o.fh3TrackPt;
1074     fh3Xi         = o.fh3Xi;
1075     fh3Z          = o.fh3Z;
1076     fh1JetPt      = o.fh1JetPt;
1077     fNameFF       = o.fNameFF;
1078   }
1079
1080   return *this;
1081 }
1082
1083 //___________________________________________________________________________
1084 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::~AliFragFuncHistosInvMass()
1085 {
1086   // destructor
1087
1088   if(fh1JetPt)   delete fh1JetPt;
1089   if(fh3TrackPt) delete fh3TrackPt;
1090   if(fh3Xi)      delete fh3Xi;
1091   if(fh3Z)       delete fh3Z;
1092 }
1093
1094 //_________________________________________________________________
1095 void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::DefineHistos()
1096 {
1097   // book FF histos
1098
1099   fh1JetPt   = new TH1F(Form("fh1FFJetPtIM%s", fNameFF.Data()),"",fNBinsJetPt,fJetPtMin,fJetPtMax);
1100   fh3TrackPt = new TH3F(Form("fh3FFTrackPtIM%s",fNameFF.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsPt, fPtMin, fPtMax);
1101   fh3Xi      = new TH3F(Form("fh3FFXiIM%s", fNameFF.Data()),"", fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsXi, fXiMin, fXiMax);
1102   fh3Z       = new TH3F(Form("fh3FFZIM%s", fNameFF.Data()),"", fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsZ, fZMin, fZMax);
1103
1104   AliAnalysisTaskFragmentationFunction::SetProperties(fh1JetPt, "p_{t} (GeV/c)", "entries");
1105   AliAnalysisTaskJetChem::SetProperties(fh3TrackPt,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","p_{t} (GeV/c)");
1106   AliAnalysisTaskJetChem::SetProperties(fh3Xi,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","#xi");
1107   AliAnalysisTaskJetChem::SetProperties(fh3Z,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","z");
1108 }
1109
1110 //________________________________________________________________________________________________________________________________
1111 void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::FillFF(Float_t trackPt, Float_t invM, Float_t jetPt, Bool_t incrementJetPt)
1112 {
1113   // fill FF
1114
1115   if(incrementJetPt) fh1JetPt->Fill(jetPt);
1116   fh3TrackPt->Fill(jetPt,invM,trackPt);//Fill(x,y,z)
1117
1118   Double_t z = 0.;
1119   if(jetPt>0) z = trackPt / jetPt;
1120   Double_t xi = 0;
1121   if(z>0) xi = TMath::Log(1/z);
1122
1123   fh3Xi->Fill(jetPt,invM,xi);
1124   fh3Z->Fill(jetPt,invM,z);
1125 }
1126
1127 //___________________________________________________________________________________
1128 void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AddToOutput(TList* list) const
1129 {
1130   // add histos to list
1131
1132   list->Add(fh1JetPt);
1133   list->Add(fh3TrackPt);
1134   list->Add(fh3Xi);
1135   list->Add(fh3Z);
1136 }
1137
1138 // ---
1139
1140
1141 //_______________________________________________________________________________________________________
1142 AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::AliFragFuncHistosPhiCorrInvMass(const char* name,
1143                                                                                          Int_t nPt,   Float_t ptMin,   Float_t ptMax,
1144                                                                                          Int_t nPhi,  Float_t phiMin,  Float_t phiMax,
1145                                                                                          Int_t nInvMass,  Float_t invMassMin,  Float_t invMassMax)
1146   : TObject()
1147   ,fNBinsPt(nPt)
1148   ,fPtMin(ptMin)
1149   ,fPtMax(ptMax)
1150   ,fNBinsPhi(nPhi)
1151   ,fPhiMin(phiMin)
1152   ,fPhiMax(phiMax)
1153   ,fNBinsInvMass(nInvMass)
1154   ,fInvMassMin(invMassMin)
1155   ,fInvMassMax(invMassMax)
1156   ,fh3PhiCorr(0)
1157   ,fNamePhiCorr(name)
1158 {
1159   // default constructor
1160 }
1161
1162 //____________________________________________________________________________________________________________________________________
1163 AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::AliFragFuncHistosPhiCorrInvMass(const AliFragFuncHistosPhiCorrInvMass& copy)
1164   : TObject()
1165   ,fNBinsPt(copy.fNBinsPt)
1166   ,fPtMin(copy.fPtMin)
1167   ,fPtMax(copy.fPtMax)
1168   ,fNBinsPhi(copy.fNBinsPhi)
1169   ,fPhiMin(copy.fPhiMin)
1170   ,fPhiMax(copy.fPhiMax)
1171   ,fNBinsInvMass(copy.fNBinsInvMass)
1172   ,fInvMassMin(copy.fInvMassMin)
1173   ,fInvMassMax(copy.fInvMassMax)
1174   ,fh3PhiCorr(copy.fh3PhiCorr)
1175   ,fNamePhiCorr(copy.fNamePhiCorr)
1176 {
1177   // copy constructor
1178 }
1179
1180 //________________________________________________________________________________________________________________________________________________________________________
1181 AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass& AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::operator=(const AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass& o)
1182 {
1183   // assignment
1184
1185   if(this!=&o){
1186     TObject::operator=(o);
1187     fNBinsPt      = o.fNBinsPt;
1188     fPtMin        = o.fPtMin;
1189     fPtMax        = o.fPtMax;
1190     fNBinsPhi     = o.fNBinsPhi;
1191     fPhiMin       = o.fPhiMin;
1192     fPhiMax       = o.fPhiMax;
1193     fNBinsInvMass = o.fNBinsInvMass;
1194     fInvMassMin   = o.fInvMassMin;
1195     fInvMassMax   = o.fInvMassMax;
1196     fh3PhiCorr    = o.fh3PhiCorr;
1197     fNamePhiCorr  = o.fNamePhiCorr;
1198   }
1199
1200   return *this;
1201 }
1202
1203 //_________________________________________________________________________________________
1204 AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::~AliFragFuncHistosPhiCorrInvMass()
1205 {
1206   // destructor
1207
1208   if(fh3PhiCorr) delete fh3PhiCorr;
1209 }
1210
1211 //__________________________________________________________________________
1212 void AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::DefineHistos()
1213 {
1214   // book jet QA histos
1215
1216   fh3PhiCorr  = new TH3F(Form("fh3PhiCorrIM%s", fNamePhiCorr.Data()),
1217                          Form("%s: p_{t} - #phi - m_{inv} distribution",fNamePhiCorr.Data()),
1218                          fNBinsPt, fPtMin, fPtMax,
1219                          fNBinsPhi, fPhiMin, fPhiMax,
1220                          fNBinsInvMass, fInvMassMin, fInvMassMax);
1221
1222   AliAnalysisTaskJetChem::SetProperties(fh3PhiCorr, "p_{t} (GeV/c)", "#phi", "m_{inv} (GeV/c^2)");
1223 }
1224
1225 //___________________________________________________________________________________________________________
1226 void AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::FillPhiCorr(Float_t pt, Float_t phi, Float_t invM)
1227 {
1228   // fill jet QA histos
1229
1230   fh3PhiCorr->Fill(pt, phi, invM);
1231 }
1232
1233 //______________________________________________________________________________________________
1234 void AliAnalysisTaskJetChem::AliFragFuncHistosPhiCorrInvMass::AddToOutput(TList* list) const
1235 {
1236   // add histos to list
1237
1238   list->Add(fh3PhiCorr);
1239 }
1240
1241 //____________________________________________________
1242 void AliAnalysisTaskJetChem::UserCreateOutputObjects()
1243 {
1244   // create output objects
1245
1246   if(fDebug > 1) Printf("AliAnalysisTaskJetChem::UserCreateOutputObjects()");
1247
1248   // create list of tracks and jets
1249
1250   fTracksRecCuts = new TList();
1251   fTracksRecCuts->SetOwner(kFALSE); //objects in TList wont be deleted when TList is deleted
1252   fJetsRecCuts = new TList();
1253   fJetsRecCuts->SetOwner(kFALSE);
1254   fBckgJetsRec = new TList();
1255   fBckgJetsRec->SetOwner(kFALSE);
1256   fListK0s = new TList();
1257   fListK0s->SetOwner(kFALSE);
1258   fListLa = new TList();
1259   fListLa->SetOwner(kFALSE);
1260   fListALa = new TList();
1261   fListALa->SetOwner(kFALSE);
1262   fListFeeddownLaCand = new TList();    //feeddown Lambda candidates
1263   fListFeeddownLaCand->SetOwner(kFALSE);
1264   fListFeeddownALaCand = new TList();   //feeddown Antilambda candidates
1265   fListFeeddownALaCand->SetOwner(kFALSE);
1266   jetConeFDLalist = new TList();
1267   jetConeFDLalist->SetOwner(kFALSE);  //feeddown Lambda candidates in jet cone
1268   jetConeFDALalist = new TList();
1269   jetConeFDALalist->SetOwner(kFALSE); //feeddown Antilambda candidates in jet cone
1270   fListMCgenK0s = new TList();          //MC generated K0s
1271   fListMCgenK0s->SetOwner(kFALSE);
1272   fListMCgenLa = new TList();           //MC generated Lambdas
1273   fListMCgenLa->SetOwner(kFALSE);
1274   fListMCgenALa = new TList();          //MC generated Antilambdas
1275   fListMCgenALa->SetOwner(kFALSE);
1276
1277
1278   // Create histograms / output container
1279
1280   //for AliPIDResponse:
1281   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
1282   AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1283   fPIDResponse = inputHandler->GetPIDResponse();
1284
1285   OpenFile(1);
1286   fCommonHistList = new TList();
1287
1288   Bool_t oldStatus = TH1::AddDirectoryStatus();
1289   TH1::AddDirectory(kFALSE);//By default (fAddDirectory = kTRUE), histograms are automatically added to the list of objects in memory
1290
1291   // histograms inherited from AliAnalysisTaskFragmentationFunction
1292
1293   fh1EvtSelection            = new TH1F("fh1EvtSelection", "Event Selection", 6, -0.5, 5.5);
1294   fh1EvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
1295   fh1EvtSelection->GetXaxis()->SetBinLabel(2,"event trigger selection: rejected");
1296   fh1EvtSelection->GetXaxis()->SetBinLabel(3,"event class: rejected");
1297   fh1EvtSelection->GetXaxis()->SetBinLabel(4,"vertex Ncontr: rejected");
1298   fh1EvtSelection->GetXaxis()->SetBinLabel(5,"vertex z: rejected");
1299   fh1EvtSelection->GetXaxis()->SetBinLabel(6,"vertex type: rejected");
1300
1301
1302   fh1EvtCent                 = new TH1F("fh1EvtCent","centrality",100,0.,100.);
1303   fh1VertexNContributors     = new TH1F("fh1VertexNContributors", "Vertex N contributors", 11,-.5, 10.5);
1304   fh1VertexZ                 = new TH1F("fh1VertexZ", "Vertex z distribution", 30, -15., 15.);
1305   fh1Xsec                    = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
1306   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
1307   fh1Trials                  = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
1308   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
1309   fh1PtHard                  = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
1310   fh1PtHardTrials            = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
1311   fh1nRecJetsCuts            = new TH1F("fh1nRecJetsCuts","reconstructed jets per event",100,-0.5,99.5);
1312
1313   // histograms JetChem task
1314
1315   fh1EvtAllCent                 = new TH1F("fh1EvtAllCent","before centrality selection",100,0.,100.);
1316   fh1Evt                        = new TH1F("fh1Evt", "All events runned over", 3, 0.,1.);
1317   fh1EvtMult                    = new TH1F("fh1EvtMult","multiplicity",1200,0.,12000.);
1318   fh1K0Mult                     = new TH1F("fh1K0Mult","K0 multiplicity",1000,0.,1000.);//500. all
1319   fh1dPhiJetK0                  = new TH1F("fh1dPhiJetK0","",640,-1,5.4);
1320   fh1LaMult                     = new TH1F("fh1LaMult","La multiplicity",1000,0.,1000.);
1321   fh1dPhiJetLa                  = new TH1F("fh1dPhiJetLa","",640,-1,5.4);
1322   fh1ALaMult                    = new TH1F("fh1ALaMult","ALa multiplicity",1000,0.,1000.);
1323   fh1dPhiJetALa                 = new TH1F("fh1dPhiJetALa","",640,-1,5.4);
1324   fh1JetEta                     = new TH1F("fh1JetEta","#eta distribution of all jets",400,-2.,2.);
1325   fh1JetPhi                     = new TH1F("fh1JetPhi","#phi distribution of all jets",630,0.,6.3);
1326   fh2JetEtaPhi                  = new TH2F("fh2JetEtaPhi","#eta and #phi distribution of all jets",400,-2.,2.,630,0.,6.3);
1327   fh1V0JetPt                    = new TH1F("fh1V0JetPt","#p_{T} distribution of all jets containing v0s",200,0.,200.);
1328   fh2FFJetTrackEta              = new TH2F("fh2FFJetTrackEta","charged track eta distr. in jet cone",200,-1.,1.,40,0.,200.);
1329   fh1trackPosNCls               = new TH1F("fh1trackPosNCls","NTPC clusters positive daughters",250,0.,250.);
1330   fh1trackNegNCls               = new TH1F("fh1trackNegNCls","NTPC clusters negative daughters",250,0.,250.);
1331   fh1trackPosEta                = new TH1F("fh1trackPosEta","eta positive daughters",100,-2.,2.);
1332   fh1trackNegEta                = new TH1F("fh1trackNegEta","eta negative daughters",100,-2.,2.);
1333   fh1V0Eta                      = new TH1F("fh1V0Eta","V0 eta",60,-1.5,1.5);
1334   fh1V0totMom                   = new TH1F("fh1V0totMom","V0 tot mom",240,0.,20.);
1335   fh1CosPointAngle              = new TH1F("fh1CosPointAngle", "Cosine of V0's pointing angle",100,0.99,1.0);
1336   fh1Chi2Pos                    = new TH1F("fh1Chi2Pos", "V0s chi2",100,0.,5.);
1337   fh1Chi2Neg                    = new TH1F("fh1Chi2Neg", "V0s chi2",100,0.,5.);
1338   fh1DecayLengthV0              = new TH1F("fh1DecayLengthV0", "V0s decay Length;decay length(cm)",1200,0.,120.);
1339   fh2ProperLifetimeK0sVsPtBeforeCut = new TH2F("fh2ProperLifetimeK0sVsPtBeforeCut"," K0s ProperLifetime vs Pt; p_{T} (GeV/#it{c})",150,0.,15.,250,0.,250.);
1340   fh2ProperLifetimeK0sVsPtAfterCut = new TH2F("fh2ProperLifetimeK0sVsPtAfterCut"," K0s ProperLifetime vs Pt; p_{T} (GeV/#it{c})",1500,0.,15.,250,0.,250.);
1341   fh1ProperLifetimeV0BeforeCut  = new TH1F("fh1ProperLifetimeV0BeforeCut", "V0s 2D distance over transerse mom.;(cm)",120,0.,120.);
1342   fh1ProperLifetimeV0AfterCut   = new TH1F("fh1ProperLifetimeV0AfterCut", "V0s 2D distance over transverse mom.;(cm)",120,0.,120.);
1343   fh1V0Radius                   = new TH1F("fh1V0Radius", "V0s Radius;Radius(cm)",200,0.,40.);
1344   fh1DcaV0Daughters             = new TH1F("fh1DcaV0Daughters", "DCA between daughters;dca(cm)",200,0.,2.);
1345   fh1DcaPosToPrimVertex         = new TH1F("fh1DcaPosToPrimVertex", "Positive V0 daughter;dca(cm)",100,0.,10.);
1346   fh1DcaNegToPrimVertex         = new TH1F("fh1DcaNegToPrimVertex", "Negative V0 daughter;dca(cm)",100,0.,10.);
1347   fh2ArmenterosBeforeCuts       = new TH2F("fh2ArmenterosBeforeCuts","Armenteros Podolanski Plot for K0s Candidates;#alpha;(p^{arm})_{T}/(GeV/#it{c})",200,-1.2,1.2,350,0.,0.35);
1348   fh2ArmenterosAfterCuts        = new TH2F("fh2ArmenterosAfterCuts","Armenteros Podolanski Plot for K0s Candidates;#alpha;(p^{arm})_{T}/(GeV/#it{c});",200,-1.2,1.2,350,0.,0.35);
1349   fh2BBLaPos                    = new TH2F("fh2BBLaPos","PID of the positive daughter of La candidates; P (GeV); -dE/dx (keV/cm ?)",100,0,10,200,0,200);
1350   fh2BBLaNeg                    = new TH2F("fh2BBLaNeg","PID of the negative daughter of La candidates; P (GeV); -dE/dx (keV/cm ?)",100,0,10,200,0,200);
1351   fh1CrossedRowsOverFindableNeg = new TH1F("fh1CrossedRowsOverFindableNeg","pos daughter crossed rows over findable in TPC;counts",20,0.,2.);
1352   fh1CrossedRowsOverFindablePos = new TH1F("fh1CrossedRowsOverFindablePos","neg daughter crossed rows over findable in TPC;counts",20,0.,2.);
1353   fh1PosDaughterCharge          = new TH1F("fh1PosDaughterCharge","charge of V0 positive daughters; V0 daughters",3,-2.,2.);
1354   fh1NegDaughterCharge          = new TH1F("fh1NegDaughterCharge","charge of V0 negative daughters; V0 daughters",3,-2.,2.);
1355   fh1PtMCK0s                    = new TH1F("fh1PtMCK0s","Pt of MC rec K0s; #it{p}_{T} (GeV/#it{c})",200,0.,20.);
1356   fh1PtMCLa                     = new TH1F("fh1PtMCLa","Pt of MC rec La; #it{p}_{T} (GeV/#it{c})",200,0.,20.);
1357   fh1PtMCALa                    = new TH1F("fh1PtMCALa","Pt of MC rec ALa; #it{p}_{T} (GeV/#it{c})",200,0.,20.);
1358   fh1EtaK0s                     = new TH1F("fh1EtaK0s","K^{0}_{s} entries ;#eta",200,-1.,1.);
1359   fh1EtaLa                      = new TH1F("fh1EtaLa","#Lambda entries ;#eta",200,-1.,1.);
1360   fh1EtaALa                     = new TH1F("fh1EtaALa","#bar{#Lambda} entries ;#eta",200,-1.,1.);
1361   fh3InvMassEtaTrackPtK0s       = new TH3F("fh3InvMassEtaTrackPtK0s","#eta; invMass (GeV/{#it{c}}^{2}); #it{p}_{T} (GeV/#it{c})", 200, -1., 1., 240, 0.4, 0.6, 140, 0., 14.);
1362   fh3InvMassEtaTrackPtLa        = new TH3F("fh3InvMassEtaTrackPtLa", "#eta; invMass (GeV/{#it{c}}^{2}; #it{p}_{T} (GeV/#it{c}))",  200, -1., 1., 140, 1.06, 1.2, 140, 0., 14.);
1363   fh3InvMassEtaTrackPtALa       = new TH3F("fh3InvMassEtaTrackPtALa","#eta; invMass (GeV/#it{c}^{2}); #it{p}_{T} (GeV/#it{c})",  200, -1., 1., 140, 1.06, 1.2, 140, 0., 14.);
1364   fh3IMK0PerpCone               = new TH3F("fh3IMK0PerpCone","{K_{0}}^{s} content in perpendicular cone",39,5.,200., 400,0.3,0.7, 200,0.,20.);
1365   fh3IMLaPerpCone               = new TH3F("fh3IMLaPerpCone","#Lambda content in perpendicular cone",39,5.,200., 140,1.06,1.2, 200,0.,20.);
1366   fh3IMALaPerpCone              = new TH3F("fh3IMALaPerpCone","#Antilambda content in perpendicular cone",39,5.,200., 140,1.06,1.2, 200,0.,20.);
1367   fh3IMK0MedianCone             = new TH3F("fh3IMK0MedianCone","{K_{0}}^{s} content in median cluster cone",39,5.,200., 400,0.3,0.7, 200,0.,20.);
1368   fh3IMLaMedianCone             = new TH3F("fh3IMLaMedianCone","#Lambda content in median cluster cone",39,5.,200., 140,1.06,1.2, 200,0.,20.);
1369   fh3IMALaMedianCone            = new TH3F("fh3IMALaMedianCone","#Antilambda content in median cluster cone",39,5.,200., 140,1.06,1.2, 200,0.,20.);
1370
1371
1372   fh1TrackMultCone          = new TH1F("fh1TrackMultCone","track multiplicity in jet cone; number of tracks",200,0.,1000.);
1373
1374   fh2TrackMultCone          = new TH2F("fh2TrackMultCone","track multiplicity in jet cone vs. jet momentum; number of tracks; jet it{p}_{T} (GeV/it{c})",200,0.,1000.,39,5.,200.);
1375
1376   fFFHistosRecCuts          = new AliFragFuncHistos("RecCuts", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1377                                                      fFFNBinsPt, fFFPtMin, fFFPtMax,
1378                                                      fFFNBinsXi, fFFXiMin, fFFXiMax,
1379                                                      fFFNBinsZ , fFFZMin , fFFZMax);
1380
1381   fV0QAK0                   = new AliFragFuncQATrackHistos("V0QAK0",fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
1382                                                             fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1383                                                             fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
1384                                                             fQATrackHighPtThreshold);
1385
1386   fFFHistosRecCutsK0Evt      = new AliFragFuncHistos("RecCutsK0Evt", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1387                                                      fFFNBinsPt, fFFPtMin, fFFPtMax,
1388                                                      fFFNBinsXi, fFFXiMin, fFFXiMax,
1389                                                      fFFNBinsZ , fFFZMin , fFFZMax);
1390
1391
1392   fFFHistosIMK0AllEvt        = new AliFragFuncHistosInvMass("K0AllEvt", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax,
1393                                                             fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
1394                                                             fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax,
1395                                                             fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,
1396                                                             fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
1397
1398   fFFHistosIMK0Jet           = new AliFragFuncHistosInvMass("K0Jet", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax,
1399                                                             fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
1400                                                             fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax,
1401                                                             fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,
1402                                                             fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
1403
1404   fFFHistosIMK0Cone          = new AliFragFuncHistosInvMass("K0Cone", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax,
1405                                                             fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
1406                                                             fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax,
1407                                                             fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,
1408                                                             fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
1409
1410   fFFHistosPhiCorrIMK0       = new AliFragFuncHistosPhiCorrInvMass("K0",fPhiCorrIMNBinsPt, fPhiCorrIMPtMin, fPhiCorrIMPtMax,
1411                                                                    fPhiCorrIMNBinsPhi, fPhiCorrIMPhiMin, fPhiCorrIMPhiMax,
1412                                                                    fPhiCorrIMNBinsInvM , fPhiCorrIMInvMMin , fPhiCorrIMInvMMax);
1413
1414   fFFHistosIMLaAllEvt        = new AliFragFuncHistosInvMass("LaAllEvt", fFFIMLaNBinsJetPt, fFFIMLaJetPtMin, fFFIMLaJetPtMax,
1415                                                             fFFIMLaNBinsInvM,fFFIMLaInvMMin,fFFIMLaInvMMax,
1416                                                             fFFIMLaNBinsPt, fFFIMLaPtMin, fFFIMLaPtMax,
1417                                                             fFFIMLaNBinsXi, fFFIMLaXiMin, fFFIMLaXiMax,
1418                                                             fFFIMLaNBinsZ , fFFIMLaZMin , fFFIMLaZMax);
1419
1420   fFFHistosIMLaJet           = new AliFragFuncHistosInvMass("LaJet", fFFIMLaNBinsJetPt, fFFIMLaJetPtMin, fFFIMLaJetPtMax,
1421                                                             fFFIMLaNBinsInvM,fFFIMLaInvMMin,fFFIMLaInvMMax,
1422                                                             fFFIMLaNBinsPt, fFFIMLaPtMin, fFFIMLaPtMax,
1423                                                             fFFIMLaNBinsXi, fFFIMLaXiMin, fFFIMLaXiMax,
1424                                                             fFFIMLaNBinsZ , fFFIMLaZMin , fFFIMLaZMax);
1425
1426
1427   fFFHistosIMLaCone          = new AliFragFuncHistosInvMass("LaCone", fFFIMLaNBinsJetPt, fFFIMLaJetPtMin, fFFIMLaJetPtMax,
1428                                                             fFFIMLaNBinsInvM,fFFIMLaInvMMin,fFFIMLaInvMMax,
1429                                                             fFFIMLaNBinsPt, fFFIMLaPtMin, fFFIMLaPtMax,
1430                                                             fFFIMLaNBinsXi, fFFIMLaXiMin, fFFIMLaXiMax,
1431                                                             fFFIMLaNBinsZ , fFFIMLaZMin , fFFIMLaZMax);
1432
1433   fFFHistosPhiCorrIMLa       = new AliFragFuncHistosPhiCorrInvMass("La",fPhiCorrIMLaNBinsPt, fPhiCorrIMLaPtMin, fPhiCorrIMLaPtMax,
1434                                                                    fPhiCorrIMLaNBinsPhi, fPhiCorrIMLaPhiMin, fPhiCorrIMLaPhiMax,
1435                                                                    fPhiCorrIMLaNBinsInvM , fPhiCorrIMLaInvMMin , fPhiCorrIMLaInvMMax);
1436
1437
1438   fFFHistosIMALaAllEvt        = new AliFragFuncHistosInvMass("ALaAllEvt", fFFIMLaNBinsJetPt, fFFIMLaJetPtMin, fFFIMLaJetPtMax,
1439                                                             fFFIMLaNBinsInvM,fFFIMLaInvMMin,fFFIMLaInvMMax,
1440                                                             fFFIMLaNBinsPt, fFFIMLaPtMin, fFFIMLaPtMax,
1441                                                             fFFIMLaNBinsXi, fFFIMLaXiMin, fFFIMLaXiMax,
1442                                                             fFFIMLaNBinsZ , fFFIMLaZMin , fFFIMLaZMax);
1443
1444   fFFHistosIMALaJet           = new AliFragFuncHistosInvMass("ALaJet", fFFIMLaNBinsJetPt, fFFIMLaJetPtMin, fFFIMLaJetPtMax,
1445                                                             fFFIMLaNBinsInvM,fFFIMLaInvMMin,fFFIMLaInvMMax,
1446                                                             fFFIMLaNBinsPt, fFFIMLaPtMin, fFFIMLaPtMax,
1447                                                             fFFIMLaNBinsXi, fFFIMLaXiMin, fFFIMLaXiMax,
1448                                                             fFFIMLaNBinsZ , fFFIMLaZMin , fFFIMLaZMax);
1449
1450   fFFHistosIMALaCone          = new AliFragFuncHistosInvMass("ALaCone", fFFIMLaNBinsJetPt, fFFIMLaJetPtMin, fFFIMLaJetPtMax,
1451                                                             fFFIMLaNBinsInvM,fFFIMLaInvMMin,fFFIMLaInvMMax,
1452                                                             fFFIMLaNBinsPt, fFFIMLaPtMin, fFFIMLaPtMax,
1453                                                             fFFIMLaNBinsXi, fFFIMLaXiMin, fFFIMLaXiMax,
1454                                                             fFFIMLaNBinsZ , fFFIMLaZMin , fFFIMLaZMax);
1455
1456   fFFHistosPhiCorrIMALa       = new AliFragFuncHistosPhiCorrInvMass("ALa",fPhiCorrIMLaNBinsPt, fPhiCorrIMLaPtMin, fPhiCorrIMLaPtMax,
1457                                                                    fPhiCorrIMLaNBinsPhi, fPhiCorrIMLaPhiMin, fPhiCorrIMLaPhiMax,
1458                                                                    fPhiCorrIMLaNBinsInvM , fPhiCorrIMLaInvMMin , fPhiCorrIMLaInvMMax);
1459
1460   //***************
1461   // MC histograms
1462   //***************
1463
1464   fh2MCgenK0Cone                = new TH2F("fh2MCgenK0Cone", "MC gen {K^{0}}^{s} #it{p}_{T}  in cone around rec jet axis versus jet #it{p}_{T}; jet #it{p}_{T}",39,5.,200.,200,0.,20.);
1465   fh2MCgenLaCone                = new TH2F("fh2MCgenLaCone", "MC gen #Lambda #it{p}_{T} in cone around rec jet axis versus jet #it{p}_{T} ; jet #it{p}_{T}",39,5.,200.,200,0.,20.);
1466   fh2MCgenALaCone               = new TH2F("fh2MCgenALaCone", "MC gen #Antilambda #it{p}_{T} in cone around rec jet axis versus jet #it{p}_{T}; jet #it{p}_{T}",39,5.,200.,200,0.,20.);
1467
1468   fh2MCgenK0Cone->GetYaxis()->SetTitle("MC gen K^{0}}^{s} #it{p}_{T}");
1469   fh2MCgenLaCone->GetYaxis()->SetTitle("MC gen #Lambda #it{p}_{T}");
1470   fh2MCgenALaCone->GetYaxis()->SetTitle("MC gen #Antilambda #it{p}_{T}");
1471
1472   fh2MCEtagenK0Cone             = new TH2F("fh2MCEtagenK0Cone","MC gen {K^{0}}^{s} #it{p}_{T} #eta distribution in jet cone;#eta",39,5.,200.,200,-1.,1.);
1473   fh2MCEtagenLaCone             = new TH2F("fh2MCEtagenLaCone","MC gen #Lambda #it{p}_{T} #eta distribution in jet cone;#eta",39,5.,200.,200,-1.,1.);
1474   fh2MCEtagenALaCone            = new TH2F("fh2MCEtagenALaCone","MC gen #Antilambda #it{p}_{T} #eta distribution in jet cone;#eta",39,5.,200.,200,-1.,1.);
1475   fh1FFIMK0ConeSmear            = new TH1F("fh1FFIMK0ConeSmear","Smeared jet pt study for K0s-in-cone-jets; smeared jet #it{p}_{T}", 39,5.,200.);
1476   fh1FFIMLaConeSmear            = new TH1F("fh1FFIMLaConeSmear","Smeared jet pt study for La-in-cone-jets; smeared jet #it{p}_{T}", 39,5.,200.);
1477   fh1FFIMALaConeSmear           = new TH1F("fh1FFIMALaConeSmear","Smeared jet pt study for ALa-in-cone-jets; smeared jet #it{p}_{T}", 39,5.,200.);
1478
1479   fh3MCrecK0Cone                = new TH3F("fh3MCrecK0Cone", "MC rec {K^{0}}^{s} #it{p}_{T}  in cone around jet axis matching MC gen particle; jet #it{p}_{T}; inv mass (GeV/#it{c}^{2};#it{p}_{T}",39,5.,200., 400,0.3,0.7, 200,0.,20.);
1480   fh3MCrecLaCone                = new TH3F("fh3MCrecLaCone", "MC rec {#Lambda #it{p}_{T}  in cone around jet axis matching MC gen particle; jet #it{p}_{T}; inv mass (GeV/#it{c}^{2});#it{p}_{T}",39,5.,200., 140,1.05,1.25, 200,0.,20.);
1481   fh3MCrecALaCone               = new TH3F("fh3MCrecALaCone", "MC rec {#Antilambda #it{p}_{T}  in cone around jet axis matching MC gen particle; jet #it{p}_{T}; inv mass (GeV/#it{c}^{2});#it{p}_{T}",39,5.,200.,140,1.05,1.25, 200,0.,20.);
1482   fh3MCrecK0ConeSmear           = new TH3F("fh3MCrecK0ConeSmear", "MC rec {K^{0}}^{s} #it{p}_{T}  in cone around jet axis matching MC gen particle, with jet p_{T} smeared; jet #it{p}_{T}; inv mass (GeV/#it{c}^{2};#it{p}_{T}",39,5.,200., 400,0.3,0.7, 200,0.,20.);
1483   fh3MCrecLaConeSmear           = new TH3F("fh3MCrecLaConeSmear", "MC rec {#Lambda #it{p}_{T}  in cone around jet axis matching MC gen particle, with jet p_{T} smeared; jet #it{p}_{T}; inv mass (GeV/#it{c}^{2});#it{p}_{T}",39,5.,200., 140,1.06,1.2, 200,0.,20.);
1484   fh3MCrecALaConeSmear          = new TH3F("fh3MCrecALaConeSmear", "MC rec {#Antilambda #it{p}_{T}  in cone around jet axis matching MC gen particle, with jet p_{T} smeared; jet #it{p}_{T}; inv mass (GeV/#it{c}^{2});#it{p}_{T}",39,5.,200.,140,1.06,1.2, 200,0.,20.);
1485   fh3SecContinCone              = new TH3F("fh3SecContinCone","secondary contamination of jet cones; jet #it{p}_{T}; track #it{p}_{T}, #eta",39,5.,200.,200,0.,20.,200,-1.,1.);
1486   fh3StrContinCone              = new TH3F("fh3StrContinCone","strange particle contamination of jet cones; jet #it{p}_{T}; track #it{p}_{T}, #eta",39,5.,200.,200,0.,20.,200,-1.,1.);
1487
1488   fh1MCMultiplicityPrimary      = new TH1F("fh1MCMultiplicityPrimary", "MC Primary Particles;NPrimary;Count", 201, -0.5, 200.5);
1489   fh1MCMultiplicityTracks       = new TH1F("h1MCMultiplicityTracks", "MC Tracks;Ntracks;Count", 201, -0.5, 200.5);
1490   // fh1MCmotherK0s             = new TH1F("fh1MCmotherK0s","K0s mother pdg codes",10,0.,10.);
1491   fh1MCmotherLa                 = new TH1F("fh1MCmotherLa","Lambdas mother pdg codes",10,0.,10.);
1492   fh1MCmotherLa->GetXaxis()->SetBinLabel(1,"#Sigma^{-}");
1493   fh1MCmotherLa->GetXaxis()->SetBinLabel(2,"#Sigma^{0}");
1494   fh1MCmotherLa->GetXaxis()->SetBinLabel(3,"#Sigma^{+}");
1495   fh1MCmotherLa->GetXaxis()->SetBinLabel(4,"#Omega^{-}");
1496   fh1MCmotherLa->GetXaxis()->SetBinLabel(5,"#Xi^{0}");
1497   fh1MCmotherLa->GetXaxis()->SetBinLabel(6,"#Xi^{-}");
1498   fh1MCmotherLa->GetXaxis()->SetBinLabel(7,"#Xi^{+}");
1499   fh1MCmotherLa->GetXaxis()->SetBinLabel(8,"primary particle");
1500   fh1MCmotherALa                = new TH1F("fh1MCmotherALa","Antilambdas mother pdg codes",10,0.,10.);
1501   fh1MCmotherALa->GetXaxis()->SetBinLabel(1,"#bar{#Sigma^{-}}");
1502   fh1MCmotherALa->GetXaxis()->SetBinLabel(2,"#bar{#Sigma^{0}}");
1503   fh1MCmotherALa->GetXaxis()->SetBinLabel(3,"#bar{#Sigma^{+}}");
1504   fh1MCmotherALa->GetXaxis()->SetBinLabel(4,"#bar{#Omega^{-}}");
1505   fh1MCmotherALa->GetXaxis()->SetBinLabel(5,"#bar{#Xi^{0}}");
1506   fh1MCmotherALa->GetXaxis()->SetBinLabel(6,"#Xi^{-}");
1507   fh1MCmotherALa->GetXaxis()->SetBinLabel(7,"#Xi^{+}");
1508   fh1MCmotherALa->GetXaxis()->SetBinLabel(8,"primary particle");
1509   fh3FeedDownLa                 = new TH3F("fh3FeedDownLa","#Lambda stemming from feeddown from Xi(0/-)", 39, 5., 200., 200, 1.05, 1.25, 200,0.,20.);
1510   fh3FeedDownALa                = new TH3F("fh3FeedDownALa","#bar#Lambda stemming from feeddown from Xibar(0/+)", 39, 5., 200., 200, 1.05, 1.25, 200, 0., 20.);
1511   fh3FeedDownLaCone             = new TH3F("fh3FeedDownLaCone","#Lambda stemming from feeddown from Xi(0/-) in jet cone", 39, 5., 200., 200, 1.05, 1.25, 200,0.,20.);
1512   fh3FeedDownALaCone            = new TH3F("fh3FeedDownALaCone","#bar#Lambda stemming from feeddown from Xibar(0/+) in jet cone", 39, 5., 200., 200, 1.05, 1.25, 200, 0., 20.);
1513   fh1MCProdRadiusK0s            = new TH1F("fh1MCProdRadiusK0s","MC gen. MC K0s prod radius",600,0.,200.);
1514   fh1MCProdRadiusLambda         = new TH1F("fh1MCProdRadiusLambda","MC gen. MC La prod radius",600,0.,200.);
1515   fh1MCProdRadiusAntiLambda     = new TH1F("fh1MCProdRadiusAntiLambda","MC gen. MC ALa prod radius",600,0.,200.);
1516
1517   // Pt and inv mass distributions
1518
1519   fh1MCPtV0s                    = new TH1F("fh1MCPtV0s", "MC gen. V^{0} in rap range;#it{p}_{T} (GeV/#it{c})",200,0,20.);// 0.1 GeV/c steps
1520   fh1MCPtK0s                    = new TH1F("fh1MCPtK0s", "MC gen. K^{0}_{s} in eta range;#it{p}_{T} (GeV/#it{c})",200,0.,20.);
1521   fh1MCPtLambda                 = new TH1F("fh1MCPtLambda", "MC gen. #Lambda in rap range;#it{p}_{T} (GeV/#it{c})",200,0.,20.);
1522   fh1MCPtAntiLambda             = new TH1F("fh1MCPtAntiLambda", "MC gen. #AntiLambda in rap range;#it{p}_{T} (GeV/#it{c})",200,0.,20.);
1523   fh1MCXiPt                     = new TH1F("fh1MCXiPt", "MC gen. #Xi^{-/o};#it{p}_{T} (GeV/#it{c})",200,0.,20.);
1524   fh1MCXibarPt                  = new TH1F("fh1MCXibarPt", "MC gen. #bar{#Xi}^{+/o};#it{p}_{T} (GeV/#it{c})",200,0.,20.);
1525   fh2MCEtaVsPtK0s               = new TH2F("fh2MCEtaVsPtK0s","MC gen. K^{0}_{s} #eta; #it{p}_{T}",200,0.,20.,200,-1.,1.);
1526   fh2MCEtaVsPtLa                = new TH2F("fh2MCEtaVsPtLa","MC gen. #Lambda #eta; #it{p}_{T}",200,0.,20.,200,-1.,1.);
1527   fh2MCEtaVsPtALa               = new TH2F("fh2MCEtaVsPtALa","MC gen. #bar{#Lambda}  #eta; #it{p}_{T}",200,0.,20.,200,-1.,1.);
1528
1529   // Rapidity
1530   fh1MCRapK0s                   = new TH1F("fh1MCRapK0s", "MC gen. K0s;rap with cut",200,-10,10);
1531   fh1MCRapLambda                = new TH1F("fh1MCRapLambda", "MC gen. #Lambda;rap",200,-10,10);
1532   fh1MCRapAntiLambda            = new TH1F("fh1MCRapAntiLambda", "MC gen. #bar{#Lambda};rap",200,-10,10);
1533   fh1MCEtaAllK0s                = new TH1F("fh1MCEtaAllK0s", "MC gen. K0s;#eta",200,-1.,1.);
1534   fh1MCEtaK0s                   = new TH1F("fh1MCEtaK0s", "MC gen. K0s;#eta with cut",200,-1.,1.);
1535   fh1MCEtaLambda                = new TH1F("fh1MCEtaLambda", "MC gen. #Lambda;#eta",200,-1.,1.);
1536   fh1MCEtaAntiLambda            = new TH1F("fh1MCEtaAntiLambda", "MC gen. #bar{#Lambda};#eta",200,-1.,1.);
1537
1538   fV0QAK0->DefineHistos();
1539   fFFHistosRecCuts->DefineHistos();
1540   fFFHistosRecCutsK0Evt->DefineHistos();
1541   fFFHistosIMK0AllEvt->DefineHistos();
1542   fFFHistosIMK0Jet->DefineHistos();
1543   fFFHistosIMK0Cone->DefineHistos();
1544   fFFHistosPhiCorrIMK0->DefineHistos();
1545   fFFHistosIMLaAllEvt->DefineHistos();
1546   fFFHistosIMLaJet->DefineHistos();
1547   fFFHistosIMLaCone->DefineHistos();
1548   fFFHistosPhiCorrIMLa->DefineHistos();
1549   fFFHistosIMALaAllEvt->DefineHistos();
1550   fFFHistosIMALaJet->DefineHistos();
1551   fFFHistosIMALaCone->DefineHistos();
1552   fFFHistosPhiCorrIMALa->DefineHistos();
1553
1554   const Int_t saveLevel = 5;
1555   if(saveLevel>0){
1556
1557     fCommonHistList->Add(fh1EvtAllCent);
1558     fCommonHistList->Add(fh1Evt);
1559     fCommonHistList->Add(fh1EvtSelection);
1560     fCommonHistList->Add(fh1EvtCent);
1561     fCommonHistList->Add(fh1VertexNContributors);
1562     fCommonHistList->Add(fh1VertexZ);
1563     fCommonHistList->Add(fh1Xsec);
1564     fCommonHistList->Add(fh1Trials);
1565     fCommonHistList->Add(fh1PtHard);
1566     fCommonHistList->Add(fh1PtHardTrials);
1567     fCommonHistList->Add(fh1nRecJetsCuts);
1568     fCommonHistList->Add(fh1EvtMult);
1569     fCommonHistList->Add(fh1K0Mult);
1570     fCommonHistList->Add(fh1dPhiJetK0);
1571     fCommonHistList->Add(fh1LaMult);
1572     fCommonHistList->Add(fh1dPhiJetLa);
1573     fCommonHistList->Add(fh1ALaMult);
1574     fCommonHistList->Add(fh1dPhiJetALa);
1575     fCommonHistList->Add(fh1JetEta);
1576     fCommonHistList->Add(fh1JetPhi);
1577     fCommonHistList->Add(fh2JetEtaPhi);
1578     fCommonHistList->Add(fh1V0JetPt);
1579     fCommonHistList->Add(fh2FFJetTrackEta);
1580     fCommonHistList->Add(fh1trackPosNCls);
1581     fCommonHistList->Add(fh1trackNegNCls);
1582     fCommonHistList->Add(fh1trackPosEta);
1583     fCommonHistList->Add(fh1trackNegEta);
1584     fCommonHistList->Add(fh1V0Eta);
1585     fCommonHistList->Add(fh1V0totMom);
1586     fCommonHistList->Add(fh1CosPointAngle);
1587     fCommonHistList->Add(fh1Chi2Pos);
1588     fCommonHistList->Add(fh1Chi2Neg);
1589     fCommonHistList->Add(fh1DecayLengthV0);
1590     fCommonHistList->Add(fh2ProperLifetimeK0sVsPtBeforeCut);
1591     fCommonHistList->Add(fh2ProperLifetimeK0sVsPtAfterCut);
1592     fCommonHistList->Add(fh1ProperLifetimeV0BeforeCut);
1593     fCommonHistList->Add(fh1ProperLifetimeV0AfterCut);
1594     fCommonHistList->Add(fh1V0Radius);
1595     fCommonHistList->Add(fh1DcaV0Daughters);
1596     fCommonHistList->Add(fh1DcaPosToPrimVertex);
1597     fCommonHistList->Add(fh1DcaNegToPrimVertex);
1598     fCommonHistList->Add(fh2ArmenterosBeforeCuts);
1599     fCommonHistList->Add(fh2ArmenterosAfterCuts);
1600     fCommonHistList->Add(fh2BBLaPos);
1601     fCommonHistList->Add(fh2BBLaNeg);
1602     fCommonHistList->Add(fh1CrossedRowsOverFindableNeg);
1603     fCommonHistList->Add(fh1CrossedRowsOverFindablePos);
1604     fCommonHistList->Add(fh1PosDaughterCharge);
1605     fCommonHistList->Add(fh1NegDaughterCharge);
1606     fCommonHistList->Add(fh1PtMCK0s);
1607     fCommonHistList->Add(fh1PtMCLa);
1608     fCommonHistList->Add(fh1PtMCALa);
1609     fCommonHistList->Add(fh1EtaK0s);
1610     fCommonHistList->Add(fh1EtaLa);
1611     fCommonHistList->Add(fh1EtaALa);
1612     fCommonHistList->Add(fh3InvMassEtaTrackPtK0s);
1613     fCommonHistList->Add(fh3InvMassEtaTrackPtLa);
1614     fCommonHistList->Add(fh3InvMassEtaTrackPtALa);
1615     fCommonHistList->Add(fh1TrackMultCone);
1616     fCommonHistList->Add(fh2TrackMultCone);
1617     fCommonHistList->Add(fh2MCgenK0Cone);
1618     fCommonHistList->Add(fh2MCgenLaCone);
1619     fCommonHistList->Add(fh2MCgenALaCone);
1620     fCommonHistList->Add(fh2MCEtagenK0Cone);
1621     fCommonHistList->Add(fh2MCEtagenLaCone);
1622     fCommonHistList->Add(fh2MCEtagenALaCone);
1623     fCommonHistList->Add(fh1FFIMK0ConeSmear);
1624     fCommonHistList->Add(fh1FFIMLaConeSmear);
1625     fCommonHistList->Add(fh1FFIMALaConeSmear);
1626     fCommonHistList->Add(fh3MCrecK0Cone);
1627     fCommonHistList->Add(fh3MCrecLaCone);
1628     fCommonHistList->Add(fh3MCrecALaCone);
1629     fCommonHistList->Add(fh3MCrecK0ConeSmear);
1630     fCommonHistList->Add(fh3MCrecLaConeSmear);
1631     fCommonHistList->Add(fh3MCrecALaConeSmear);
1632     fCommonHistList->Add(fh3SecContinCone);
1633     fCommonHistList->Add(fh3StrContinCone);
1634     fCommonHistList->Add(fh3IMK0PerpCone);
1635     fCommonHistList->Add(fh3IMLaPerpCone);
1636     fCommonHistList->Add(fh3IMALaPerpCone);
1637     fCommonHistList->Add(fh3IMK0MedianCone);
1638     fCommonHistList->Add(fh3IMLaMedianCone);
1639     fCommonHistList->Add(fh3IMALaMedianCone);
1640     fCommonHistList->Add(fh1MCMultiplicityPrimary);
1641     fCommonHistList->Add(fh1MCMultiplicityTracks);
1642     fCommonHistList->Add(fh1MCmotherLa);
1643     fCommonHistList->Add(fh1MCmotherALa);
1644     fCommonHistList->Add(fh3FeedDownLa);
1645     fCommonHistList->Add(fh3FeedDownALa);
1646     fCommonHistList->Add(fh3FeedDownLaCone);
1647     fCommonHistList->Add(fh3FeedDownALaCone);
1648     fCommonHistList->Add(fh1MCProdRadiusK0s);
1649     fCommonHistList->Add(fh1MCProdRadiusLambda);
1650     fCommonHistList->Add(fh1MCProdRadiusAntiLambda);
1651     fCommonHistList->Add(fh1MCPtV0s);
1652     fCommonHistList->Add(fh1MCPtK0s);
1653     fCommonHistList->Add(fh1MCPtLambda);
1654     fCommonHistList->Add(fh1MCPtAntiLambda);
1655     fCommonHistList->Add(fh1MCXiPt);
1656     fCommonHistList->Add(fh1MCXibarPt);
1657     fCommonHistList->Add(fh2MCEtaVsPtK0s);
1658     fCommonHistList->Add(fh2MCEtaVsPtLa);
1659     fCommonHistList->Add(fh2MCEtaVsPtALa);
1660     fCommonHistList->Add(fh1MCRapK0s);
1661     fCommonHistList->Add(fh1MCRapLambda);
1662     fCommonHistList->Add(fh1MCRapAntiLambda);
1663     fCommonHistList->Add(fh1MCEtaAllK0s);
1664     fCommonHistList->Add(fh1MCEtaK0s);
1665     fCommonHistList->Add(fh1MCEtaLambda);
1666     fCommonHistList->Add(fh1MCEtaAntiLambda);
1667
1668     fV0QAK0->AddToOutput(fCommonHistList);
1669     fFFHistosRecCuts->AddToOutput(fCommonHistList);
1670     fFFHistosRecCutsK0Evt->AddToOutput(fCommonHistList);
1671     fFFHistosIMK0AllEvt->AddToOutput(fCommonHistList);
1672     fFFHistosIMK0Jet->AddToOutput(fCommonHistList);
1673     fFFHistosIMK0Cone->AddToOutput(fCommonHistList);
1674     fFFHistosPhiCorrIMK0->AddToOutput(fCommonHistList);
1675     fFFHistosIMLaAllEvt->AddToOutput(fCommonHistList);
1676     fFFHistosIMLaJet->AddToOutput(fCommonHistList);
1677     fFFHistosIMLaCone->AddToOutput(fCommonHistList);
1678     fFFHistosPhiCorrIMLa->AddToOutput(fCommonHistList);
1679     fFFHistosIMALaAllEvt->AddToOutput(fCommonHistList);
1680     fFFHistosIMALaJet->AddToOutput(fCommonHistList);
1681     fFFHistosIMALaCone->AddToOutput(fCommonHistList);
1682     fFFHistosPhiCorrIMALa->AddToOutput(fCommonHistList);
1683
1684
1685   }
1686
1687    // =========== Switch on Sumw2 for all histos ===========
1688    for (Int_t i=0; i<fCommonHistList->GetEntries(); ++i){
1689
1690    TH1 *h1 = dynamic_cast<TH1*>(fCommonHistList->At(i));
1691
1692    if (h1) h1->Sumw2();//The error per bin will be computed as sqrt(sum of squares of weight) for each bin
1693     else{
1694       THnSparse *hnSparse = dynamic_cast<THnSparse*>(fCommonHistList->At(i));
1695       if(hnSparse) hnSparse->Sumw2();
1696     }
1697
1698   }
1699   TH1::AddDirectory(oldStatus);
1700
1701 }
1702
1703 //_______________________________________________
1704 void AliAnalysisTaskJetChem::UserExec(Option_t *)
1705 {
1706   // Main loop
1707   // Called for each event
1708
1709   if(fDebug > 1) Printf("AliAnalysisTaskJetChem::UserExec()");
1710
1711   if(fDebug > 1) Printf("Analysis event #%5d", (Int_t) fEntry);
1712
1713    // Trigger selection
1714   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
1715     ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1716
1717   //std::cout<<"inputHandler->IsEventSelected(): "<<inputHandler->IsEventSelected()<<std::endl;
1718   //std::cout<<"fEvtSelectionMask: "<<fEvtSelectionMask<<std::endl;
1719
1720   if(!(inputHandler->IsEventSelected() & fEvtSelectionMask)){
1721     //std::cout<<"########event rejected!!############"<<std::endl;
1722     fh1EvtSelection->Fill(1.);
1723     if (fDebug > 1 ) Printf(" Trigger Selection: event REJECTED ... ");
1724     PostData(1, fCommonHistList);
1725     return;
1726   }
1727
1728   fESD = dynamic_cast<AliESDEvent*>(InputEvent());//casting of pointers for inherited class, only for ESDs
1729   if(!fESD){
1730     if(fDebug>3) Printf("%s:%d ESDEvent not found in the input", (char*)__FILE__,__LINE__);
1731   }
1732
1733   fMCEvent = MCEvent();
1734   if(!fMCEvent){
1735     if(fDebug>3) Printf("%s:%d MCEvent not found in the input", (char*)__FILE__,__LINE__);
1736   }
1737
1738   // get AOD event from input/output
1739   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
1740   if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
1741     fAOD  =  ((AliAODInputHandler*)handler)->GetEvent();
1742     if(fUseAODInputJets) fAODJets = fAOD;
1743     if (fDebug > 1)  Printf("%s:%d AOD event from input", (char*)__FILE__,__LINE__);
1744   }
1745   else {
1746     handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
1747     if( handler && handler->InheritsFrom("AliAODHandler") ) {
1748       fAOD = ((AliAODHandler*)handler)->GetAOD();
1749       fAODJets = fAOD;
1750       if (fDebug > 1)  Printf("%s:%d AOD event from output", (char*)__FILE__,__LINE__);
1751     }
1752   }
1753
1754   if(!fAODJets && !fUseAODInputJets){ // case we have AOD in input & output and want jets from output
1755     TObject* outHandler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
1756     if( outHandler && outHandler->InheritsFrom("AliAODHandler") ){
1757       fAODJets = ((AliAODHandler*)outHandler)->GetAOD();
1758       if (fDebug > 1)  Printf("%s:%d jets from output AOD", (char*)__FILE__,__LINE__);
1759     }
1760   }
1761
1762   if(fNonStdFile.Length()!=0){
1763     // case we have an AOD extension - fetch the jets from the extended output
1764
1765     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1766     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
1767     if(!fAODExtension){
1768       if(fDebug>1)Printf("AODExtension not found for %s",fNonStdFile.Data());
1769     }
1770   }
1771
1772   if(!fAOD){
1773     Printf("%s:%d AODEvent not found", (char*)__FILE__,__LINE__);
1774     return;
1775   }
1776   if(!fAODJets){
1777     Printf("%s:%d AODEvent with jet branch not found", (char*)__FILE__,__LINE__);
1778     return;
1779   }
1780
1781   //primary vertex position:
1782   AliAODVertex *myPrimaryVertex = NULL;
1783   myPrimaryVertex = (AliAODVertex*)fAOD->GetPrimaryVertex();
1784   if (!myPrimaryVertex) return;
1785   fh1Evt->Fill(1.);//fill in every event that was accessed with InputHandler
1786
1787   // event selection  *****************************************
1788
1789   // *** vertex cut ***
1790   AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
1791   Int_t nTracksPrim = primVtx->GetNContributors();
1792   fh1VertexNContributors->Fill(nTracksPrim);
1793
1794   if (fDebug > 1) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim);
1795   //if(!nTracksPrim){
1796   if(nTracksPrim <= 2){
1797     if (fDebug > 1) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__);
1798     fh1EvtSelection->Fill(3.);
1799     PostData(1, fCommonHistList);
1800     return;
1801   }
1802
1803   fh1VertexZ->Fill(primVtx->GetZ());
1804
1805   if(TMath::Abs(primVtx->GetZ())>fMaxVertexZ){
1806     if (fDebug > 1) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
1807     fh1EvtSelection->Fill(4.);
1808     PostData(1, fCommonHistList);
1809     return;
1810   }
1811
1812   // accepts only events that have same "primary" and SPD vertex, special issue of LHC11h PbPb data
1813
1814   //fAOD: pointer to global primary vertex
1815
1816   const AliAODVertex* spdVtx = fAOD->GetPrimaryVertexSPD();
1817
1818   if (TMath::Abs(spdVtx->GetZ() - primVtx->GetZ())>fDeltaVertexZ) { if (fDebug > 1) Printf("deltaZVertex: event REJECTED..."); return;}
1819
1820
1821   //check for vertex radius to be smaller than 1 cm, (that was first applied by Vit Kucera in his analysis)
1822
1823   Double_t vtxX = primVtx->GetX();
1824   Double_t vtxY = primVtx->GetY();
1825
1826   if(TMath::Sqrt(vtxX*vtxX + vtxY*vtxY)>=1){
1827     if (fDebug > 1) Printf("%s:%d primary vertex r = %f: event REJECTED...",(char*)__FILE__,__LINE__,TMath::Sqrt(vtxX*vtxX + vtxY*vtxY));
1828     return;
1829   }
1830
1831
1832   TString primVtxName(primVtx->GetName());
1833
1834   if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){
1835     if (fDebug > 1) Printf("%s:%d primary vertex selection: TPC vertex, event REJECTED...",(char*)__FILE__,__LINE__);
1836     fh1EvtSelection->Fill(5.);
1837     PostData(1, fCommonHistList);
1838     return;
1839   }
1840
1841   Bool_t selectedHelper = AliAnalysisHelperJetTasks::Selected();
1842   if(!selectedHelper){
1843     fh1EvtSelection->Fill(6.);
1844     PostData(1, fCommonHistList);
1845     return;
1846   }
1847
1848   // event selection  *****************************************
1849
1850   Double_t centPercent = -1;
1851   Int_t cl = 0;
1852   if(fEventClass>0){
1853
1854     if(handler && handler->InheritsFrom("AliAODInputHandler")){
1855
1856       centPercent = fAOD->GetHeader()->GetCentrality();
1857       cl = 1;
1858       //std::cout<<"centPercent: "<<centPercent<<std::endl;
1859
1860       fh1EvtAllCent->Fill(centPercent);
1861       /*
1862       if(centPercent>10) cl = 2; //standard PWG-JE binning
1863       if(centPercent>30) cl = 3;
1864       if(centPercent>50) cl = 4;
1865       */
1866
1867
1868       if(centPercent < 0) cl = -1;
1869       if(centPercent >= 0)  cl = 1;
1870       if(centPercent > 10) cl = 2; //standard PWG-JE binning
1871       if(centPercent > 30) cl = 3;
1872       if(centPercent > 50) cl = 4;
1873       if(centPercent > 80) cl = 5; //takes centralities higher than my upper edge of 80%, not to be used
1874
1875     }
1876     else {
1877
1878       cl = AliAnalysisHelperJetTasks::EventClass();
1879
1880       if(fESD) centPercent = fESD->GetCentrality()->GetCentralityPercentile("V0M"); //ESD JetServices Task has the centrality binning 0-10,10-30,30-50,50-80
1881       fh1EvtAllCent->Fill(centPercent);
1882     }
1883
1884     if(cl!=fEventClass){ // event not in selected event class, reject event#########################################
1885
1886       if (fDebug > 1) Printf("%s:%d event not in selected event class: event REJECTED ...",(char*)__FILE__,__LINE__);
1887       fh1EvtSelection->Fill(2.);
1888       PostData(1, fCommonHistList);
1889       return;
1890     }
1891   }//end if fEventClass > 0
1892
1893
1894   if (fDebug > 1) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__);
1895
1896   //test test
1897   //Printf("Analysis event #%5d", (Int_t) fEntry);
1898
1899   fh1EvtSelection->Fill(0.);
1900   fh1EvtCent->Fill(centPercent);
1901
1902   //___ get MC information __________________________________________________________________
1903
1904
1905   Double_t ptHard = 0.; //parton energy bins -> energy of particle
1906   Double_t nTrials = 1; // trials for MC trigger weight for real data
1907
1908   if(fMCEvent){
1909      AliGenEventHeader* genHeader = fMCEvent->GenEventHeader();
1910      AliGenPythiaEventHeader*  pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);//check usage of Pythia (pp) or Hijing (PbPb)
1911      AliGenHijingEventHeader*  hijingGenHeader = 0x0;
1912
1913      if(pythiaGenHeader){
1914        if(fDebug>3) Printf("%s:%d pythiaGenHeader found", (char*)__FILE__,__LINE__);
1915        nTrials = pythiaGenHeader->Trials();
1916        ptHard  = pythiaGenHeader->GetPtHard();
1917
1918        fh1PtHard->Fill(ptHard);
1919        fh1PtHardTrials->Fill(ptHard,nTrials);
1920
1921
1922      } else { // no pythia, hijing?
1923
1924        if(fDebug>3) Printf("%s:%d no pythiaGenHeader found", (char*)__FILE__,__LINE__);
1925
1926        hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
1927        if(!hijingGenHeader){
1928          if(fDebug>3) Printf("%s:%d no pythiaGenHeader or hjingGenHeader found", (char*)__FILE__,__LINE__);
1929        } else {
1930          if(fDebug>3) Printf("%s:%d hijingGenHeader found", (char*)__FILE__,__LINE__);
1931        }
1932      }
1933
1934      fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
1935   }
1936
1937     //____ fetch jets _______________________________________________________________
1938
1939   Int_t nJCuts = GetListOfJets(fJetsRecCuts, kJetsRecAcceptance);//fetch list with jets
1940
1941   Int_t nRecJetsCuts = 0;                                        //number of reconstructed jets after jet cuts
1942   if(nJCuts>=0) nRecJetsCuts = fJetsRecCuts->GetEntries();
1943   if(fDebug>2)Printf("%s:%d Selected Rec jets after cuts: %d %d",(char*)__FILE__,__LINE__,nJCuts,nRecJetsCuts);
1944   if(nRecJetsCuts != nJCuts) Printf("%s:%d Mismatch selected Rec jets after cuts: %d %d",(char*)__FILE__,__LINE__,nJCuts,nRecJetsCuts);
1945   fh1nRecJetsCuts->Fill(nRecJetsCuts);
1946
1947
1948   //____ fetch background clusters ___________________________________________________
1949   if(fBranchRecBckgClusters.Length() != 0){
1950
1951     Int_t nBJ = GetListOfBckgJets(fBckgJetsRec, kJetsRec);
1952     Int_t nRecBckgJets = 0;
1953     if(nBJ>=0) nRecBckgJets = fBckgJetsRec->GetEntries();
1954     if(fDebug>2)Printf("%s:%d Selected Rec background jets: %d %d",(char*)__FILE__,__LINE__,nBJ,nRecBckgJets);
1955     if(nBJ != nRecBckgJets) Printf("%s:%d Mismatch Selected Rec background jets: %d %d",(char*)__FILE__,__LINE__,nBJ,nRecBckgJets);
1956   }
1957
1958
1959   //____ fetch reconstructed particles __________________________________________________________
1960
1961   Int_t nTCuts = GetListOfTracks(fTracksRecCuts, kTrackAODCuts);//all tracks of event
1962   if(fDebug>2)Printf("%s:%d selected reconstructed tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,fTracksRecCuts->GetEntries());
1963   if(fTracksRecCuts->GetEntries() != nTCuts)
1964     Printf("%s:%d Mismatch selected reconstructed tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,fTracksRecCuts->GetEntries());
1965   fh1EvtMult->Fill(fTracksRecCuts->GetEntries());
1966
1967   Int_t nK0s = GetListOfV0s(fListK0s,fK0Type,kK0,myPrimaryVertex,fAOD);//all V0s in event with K0s assumption
1968
1969   if(fDebug>5){std::cout<<"fK0Type: "<<fK0Type<<" kK0: "<<kK0<<" myPrimaryVertex: "<<myPrimaryVertex<<" fAOD:  "<<fAOD<<std::endl;}
1970
1971   //std::cout<< "nK0s: "<<nK0s<<std::endl;
1972
1973   if(fDebug>2)Printf("%s:%d Selected V0s after cuts: %d %d",(char*)__FILE__,__LINE__,nK0s,fListK0s->GetEntries());
1974   if(nK0s != fListK0s->GetEntries()) Printf("%s:%d Mismatch selected K0s: %d %d",(char*)__FILE__,__LINE__,nK0s,fListK0s->GetEntries());
1975   fh1K0Mult->Fill(fListK0s->GetEntries());
1976
1977
1978   Int_t nLa = GetListOfV0s(fListLa,fLaType,kLambda,myPrimaryVertex,fAOD);//all V0s in event with Lambda particle assumption
1979   if(fDebug>2)Printf("%s:%d Selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nLa,fListLa->GetEntries());
1980   if(nLa != fListLa->GetEntries()) Printf("%s:%d Mismatch selected La: %d %d",(char*)__FILE__,__LINE__,nLa,fListLa->GetEntries());
1981   fh1LaMult->Fill(fListLa->GetEntries());
1982
1983   Int_t nALa = GetListOfV0s(fListALa,fALaType,kAntiLambda,myPrimaryVertex,fAOD);//all V0s in event with Antilambda particle assumption
1984   if(fDebug>2)Printf("%s:%d Selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nALa,fListALa->GetEntries());
1985   if(nALa != fListALa->GetEntries()) Printf("%s:%d Mismatch selected ALa: %d %d",(char*)__FILE__,__LINE__,nALa,fListALa->GetEntries());
1986   fh1ALaMult->Fill(fListALa->GetEntries());
1987
1988
1989
1990   //fetch MC gen particles_______________________________________________________
1991
1992   if(fAnalysisMC){ // here
1993
1994     //fill feeddown histo for associated particles
1995
1996     // Access MC generated particles, fill TLists and histograms :
1997
1998     Int_t nMCgenK0s = GetListOfMCParticles(fListMCgenK0s,kK0,fAOD); //fill TList with MC generated primary true K0s (list to fill, particletype, mc aod event)
1999     if(nMCgenK0s != fListMCgenK0s->GetEntries()) Printf("%s:%d Mismatch selected MCgenK0s: %d %d",(char*)__FILE__,__LINE__,nMCgenK0s,fListMCgenK0s->GetEntries());
2000
2001
2002     for(Int_t it=0; it<fListMCgenK0s->GetSize(); ++it){ // loop MC generated K0s, filling histograms
2003
2004       AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenK0s->At(it));
2005       if(!mcp0) continue;
2006
2007       //MC gen K0s
2008
2009       Double_t fRapCurrentPart   = MyRapidity(mcp0->E(),mcp0->Pz());
2010       Double_t fEtaCurrentPart   = mcp0->Eta();
2011       Double_t fPtCurrentPart    = mcp0->Pt();
2012
2013       fh1MCEtaK0s->Fill(fEtaCurrentPart);
2014       fh1MCRapK0s->Fill(fRapCurrentPart);
2015       fh1MCPtK0s->Fill(fPtCurrentPart);
2016
2017       fh2MCEtaVsPtK0s->Fill(fPtCurrentPart,fEtaCurrentPart);                  //eta cut, physical primary selection and decay mode considered
2018
2019     }//end of the loop
2020
2021
2022     Int_t nMCgenLa = GetListOfMCParticles(fListMCgenLa,kLambda,fAOD); //fill TList with MC generated primary true Lambdas (list to fill, particletype, mc aod event)
2023     if(nMCgenLa != fListMCgenLa->GetEntries()) Printf("%s:%d Mismatch selected MCgenLa: %d %d",(char*)__FILE__,__LINE__,nMCgenLa,fListMCgenLa->GetEntries());
2024
2025
2026     for(Int_t it=0; it<fListMCgenLa->GetSize(); ++it){ // loop MC generated La, filling histograms
2027
2028       AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenLa->At(it));
2029       if(!mcp0) continue;
2030
2031       //MC gen Lambdas
2032
2033       Double_t fRapCurrentPart   = MyRapidity(mcp0->E(),mcp0->Pz());
2034       Double_t fEtaCurrentPart   = mcp0->Eta();
2035       Double_t fPtCurrentPart    = mcp0->Pt();
2036
2037       fh1MCEtaLambda->Fill(fEtaCurrentPart);
2038       fh1MCRapLambda->Fill(fRapCurrentPart);
2039       fh1MCPtLambda->Fill(fPtCurrentPart);
2040       fh2MCEtaVsPtLa->Fill(fPtCurrentPart,fEtaCurrentPart);                  //eta cut, physical primary selection and decay mode considered
2041
2042     }//end of the loop
2043
2044
2045     Int_t nMCgenALa = GetListOfMCParticles(fListMCgenALa,kAntiLambda,fAOD); //fill TList with MC generated primary true Antilambdas (list to fill, particletype, mc aod event)
2046     if(nMCgenALa != fListMCgenALa->GetEntries()) Printf("%s:%d Mismatch selected MCgenALa: %d %d",(char*)__FILE__,__LINE__,nMCgenALa,fListMCgenALa->GetEntries());
2047
2048
2049     for(Int_t it=0; it<fListMCgenALa->GetSize(); ++it){ // loop MC generated ALa, filling histograms
2050
2051       AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenALa->At(it));
2052       if(!mcp0) continue;
2053
2054       //MC gen Antilambdas
2055
2056       Double_t fRapCurrentPart   = MyRapidity(mcp0->E(),mcp0->Pz());
2057       Double_t fEtaCurrentPart   = mcp0->Eta();
2058       Double_t fPtCurrentPart    = mcp0->Pt();
2059
2060       fh1MCEtaAntiLambda->Fill(fEtaCurrentPart);
2061       fh1MCRapAntiLambda->Fill(fRapCurrentPart);
2062       fh1MCPtAntiLambda->Fill(fPtCurrentPart);
2063       fh2MCEtaVsPtALa->Fill(fPtCurrentPart,fEtaCurrentPart);                  //eta cut, physical primary selection and decay mode considered
2064
2065     }//end of the loop
2066
2067
2068
2069   //loop over MC feeddown candidates in TList
2070
2071     //....
2072
2073
2074   } //end MCAnalysis part for gen particles
2075
2076
2077   // ___ V0 QA + K0s + La + ALa pt spectra all events _______________________________________________
2078
2079   Double_t lPrimaryVtxPosition[3];
2080   Double_t lV0Position[3];
2081   lPrimaryVtxPosition[0] = primVtx->GetX();
2082   lPrimaryVtxPosition[1] = primVtx->GetY();
2083   lPrimaryVtxPosition[2] = primVtx->GetZ();
2084
2085   //------------------------------------------
2086
2087   for(Int_t it=0; it<fListK0s->GetSize(); ++it){ // loop all K0s
2088
2089     AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListK0s->At(it));
2090     if(!v0) continue;
2091
2092     // VO's main characteristics to check the reconstruction cuts
2093
2094     Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2095     Double_t invMK0s=0;
2096     Double_t trackPt=0;
2097     Double_t fV0Radius      = -999;
2098     Double_t fDcaV0Daughters = v0->DcaV0Daughters();
2099     Double_t fDcaPosToPrimVertex = v0->DcaPosToPrimVertex();//IP of positive charged daughter
2100     Double_t fDcaNegToPrimVertex = v0->DcaNegToPrimVertex();//IP of negative charged daughter
2101     Int_t negDaughterpdg = 0;
2102     Int_t posDaughterpdg = 0;
2103     Int_t motherType = 0;
2104     Int_t v0Label = -1;
2105     Double_t MCPt = 0;
2106     Bool_t fPhysicalPrimary = kFALSE;//don't use IsPhysicalPrimary() anymore for MC analysis, use instead 2D distance from primary to secondary vertex
2107     Int_t MCv0PdgCode = 0;
2108
2109     AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0));
2110     AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1));
2111
2112     Double_t PosEta = trackPos->AliAODTrack::Eta();//daughter track charge is sometimes wrong here, account for that!!!
2113     Double_t NegEta = trackNeg->AliAODTrack::Eta();
2114
2115     //Double_t trackPosNcls = trackPos->GetTPCNcls();//Get number of clusters for positive charged tracks
2116     //Double_t trackNegNcls = trackNeg->GetTPCNcls();//Get number of clusters for negative charged tracks
2117
2118     CalculateInvMass(v0, kK0, invMK0s, trackPt);  //function to calculate invMass with TLorentzVector class
2119
2120     Double_t jetPt = fFFIMJetPtMin; // assign pro forma jet energy
2121     //Double_t fRap = v0->RapK0Short();
2122     Double_t fEta = v0->PseudoRapV0();
2123
2124     Double_t fV0cosPointAngle = v0->CosPointingAngle(lPrimaryVtxPosition);
2125
2126     lV0Position[0]= v0->DecayVertexV0X();
2127     lV0Position[1]= v0->DecayVertexV0Y();
2128     lV0Position[2]= v0->DecayVertexV0Z();
2129
2130     Double_t fV0mom[3];
2131
2132     fV0mom[0]=v0->MomV0X();
2133     fV0mom[1]=v0->MomV0Y();
2134     fV0mom[2]=v0->MomV0Z();
2135     Double_t fV0TotalMomentum = TMath::Sqrt(fV0mom[0]*fV0mom[0]+fV0mom[1]*fV0mom[1]+fV0mom[2]*fV0mom[2]);
2136     Double_t fV0DecayLength = v0->DecayLengthV0(lPrimaryVtxPosition);
2137     fV0Radius  = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
2138
2139     fV0QAK0->FillTrackQA(v0->Eta(), TVector2::Phi_0_2pi(v0->Phi()), v0->Pt());
2140     fFFHistosIMK0AllEvt->FillFF(trackPt, invMK0s, jetPt, incrementJetPt);
2141     //fh1trackPosNCls->Fill(trackPosNcls);
2142     //fh1trackNegNCls->Fill(trackNegNcls);
2143     fh1EtaK0s->Fill(fEta);
2144     if(fAnalysisMC){
2145       TList *listmc = fAOD->GetList();
2146       Bool_t mclabelcheck = MCLabelCheck(v0, kK0, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode);
2147       //if(fPhysicalPrimary == kFALSE)continue;
2148       //std::cout<<"mclabelcheck: "<<mclabelcheck<<std::endl;
2149       //std::cout<<"IsPhysicalPrimary: "<<fPhysicalPrimary<<std::endl;
2150
2151       if(mclabelcheck == kFALSE)continue;
2152       fh3InvMassEtaTrackPtK0s->Fill(fEta,invMK0s,trackPt);//includes also feeddown particles
2153       fh1PtMCK0s->Fill(MCPt);
2154     }
2155
2156
2157     fh1V0Eta->Fill(fEta);
2158     fh1V0totMom->Fill(fV0TotalMomentum);
2159     fh1CosPointAngle->Fill(fV0cosPointAngle);
2160     fh1DecayLengthV0->Fill(fV0DecayLength);
2161     fh1V0Radius->Fill(fV0Radius);
2162     fh1DcaV0Daughters->Fill(fDcaV0Daughters);
2163     fh1DcaPosToPrimVertex->Fill(fDcaPosToPrimVertex);
2164     fh1DcaNegToPrimVertex->Fill(fDcaNegToPrimVertex);
2165     fh1trackPosEta->Fill(PosEta);
2166     fh1trackNegEta->Fill(NegEta);
2167   }
2168
2169
2170   // __La pt spectra all events _______________________________________________
2171
2172
2173   for(Int_t it=0; it<fListLa->GetSize(); ++it){ // loop all La
2174
2175     AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListLa->At(it));
2176     if(!v0) continue;
2177
2178     // VO's main characteristics to check the reconstruction cuts
2179     Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2180     Double_t invMLa =0;
2181     Double_t trackPt=0;
2182     Double_t fV0Radius      = -999;
2183     Double_t fDcaV0Daughters = v0->DcaV0Daughters();
2184     Double_t fDcaPosToPrimVertex = v0->DcaPosToPrimVertex();//IP of positive charged daughter
2185     Double_t fDcaNegToPrimVertex = v0->DcaNegToPrimVertex();//IP of negative charged daughter
2186     Int_t negDaughterpdg = 0;
2187     Int_t posDaughterpdg = 0;
2188     Int_t motherType = 0;
2189     Int_t v0Label = -1;
2190     Double_t MCPt = 0;
2191     Bool_t fPhysicalPrimary = kFALSE;
2192     Int_t MCv0PdgCode = 0;
2193     AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0));
2194     AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1));
2195
2196     //Double_t trackPosNcls = trackPos->GetTPCNcls();//Get number of clusters for positive charged tracks
2197     //Double_t trackNegNcls = trackNeg->GetTPCNcls();//Get number of clusters for negative charged tracks
2198
2199     Double_t PosEta = trackPos->AliAODTrack::Eta();//daughter track charge is sometimes wrong here, account for that!!!
2200     Double_t NegEta = trackNeg->AliAODTrack::Eta();
2201
2202     CalculateInvMass(v0, kLambda, invMLa, trackPt);//function to calculate invMass with TLorentzVector class
2203
2204
2205     Double_t jetPt = fFFIMJetPtMin; // assign pro forma jet energy
2206     // Double_t fRap = v0->Y(3122);
2207     Double_t fEta = v0->PseudoRapV0();
2208
2209     Double_t fV0mom[3];
2210
2211     fV0mom[0]=v0->MomV0X();
2212     fV0mom[1]=v0->MomV0Y();
2213     fV0mom[2]=v0->MomV0Z();
2214     Double_t fV0TotalMomentum = TMath::Sqrt(fV0mom[0]*fV0mom[0]+fV0mom[1]*fV0mom[1]+fV0mom[2]*fV0mom[2]);
2215     Double_t fV0DecayLength = v0->DecayLengthV0(lPrimaryVtxPosition);
2216     Double_t fV0cosPointAngle = v0->CosPointingAngle(lPrimaryVtxPosition);
2217     lV0Position[0]= v0->DecayVertexV0X();
2218     lV0Position[1]= v0->DecayVertexV0Y();
2219     lV0Position[2]= v0->DecayVertexV0Z();
2220
2221     fV0Radius  = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
2222
2223     fFFHistosIMLaAllEvt->FillFF(trackPt, invMLa, jetPt, incrementJetPt);
2224     //fh1trackPosNCls->Fill(trackPosNcls);
2225     //fh1trackNegNCls->Fill(trackNegNcls);
2226     fh1EtaLa->Fill(fEta);
2227     if(fAnalysisMC){
2228       TList* listmc = fAOD->GetList();
2229       Bool_t mclabelcheck = MCLabelCheck(v0, kLambda, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode);
2230       if(mclabelcheck == kFALSE)continue;
2231       //if(fPhysicalPrimary == kFALSE)continue;
2232        fh3InvMassEtaTrackPtLa->Fill(fEta,invMLa,trackPt);
2233       fh1PtMCLa->Fill(MCPt);
2234     }
2235     fh1V0Eta->Fill(fEta);
2236     fh1V0totMom->Fill(fV0TotalMomentum);
2237     fh1CosPointAngle->Fill(fV0cosPointAngle);
2238     fh1DecayLengthV0->Fill(fV0DecayLength);
2239     fh1V0Radius->Fill(fV0Radius);
2240     fh1DcaV0Daughters->Fill(fDcaV0Daughters);
2241     fh1DcaPosToPrimVertex->Fill(fDcaPosToPrimVertex);
2242     fh1DcaNegToPrimVertex->Fill(fDcaNegToPrimVertex);
2243     fh1trackPosEta->Fill(PosEta);
2244     fh1trackNegEta->Fill(NegEta);
2245   }
2246
2247   // __ALa pt spectra all events _______________________________________________
2248
2249   for(Int_t it=0; it<fListALa->GetSize(); ++it){ // loop all ALa
2250
2251     AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListALa->At(it));
2252     if(!v0) continue;
2253
2254
2255     //VO's main characteristics to check the reconstruction cuts
2256     Double_t invMALa =0;
2257     Double_t trackPt=0;
2258     Double_t fV0Radius      = -999;
2259     Double_t fDcaV0Daughters = v0->DcaV0Daughters();
2260     Double_t fDcaPosToPrimVertex = v0->DcaPosToPrimVertex();//IP of positive charged daughter
2261     Double_t fDcaNegToPrimVertex = v0->DcaNegToPrimVertex();//IP of negative charged daughter
2262     Int_t negDaughterpdg = 0;
2263     Int_t posDaughterpdg = 0;
2264     Int_t motherType = 0;
2265     Int_t v0Label = -1;
2266     Double_t MCPt = 0;
2267     Bool_t fPhysicalPrimary = kFALSE;
2268     Int_t MCv0PdgCode = 0;
2269
2270     AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0));
2271     AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1));
2272
2273     Double_t PosEta = trackPos->AliAODTrack::Eta();//daughter track charge is sometimes wrong here, account for that!!!
2274     Double_t NegEta = trackNeg->AliAODTrack::Eta();
2275
2276     //Double_t trackPosNcls = trackPos->GetTPCNcls();//Get number of clusters for positive charged tracks //not used anymore by Strangeness PAG group
2277     //Double_t trackNegNcls = trackNeg->GetTPCNcls();//Get number of clusters for negative charged tracks
2278
2279     CalculateInvMass(v0, kAntiLambda, invMALa, trackPt);  //function to calculate invMass with TLorentzVector class
2280
2281     Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2282     Double_t jetPt = fFFIMJetPtMin; // assign pro forma jet energy
2283     //      Double_t fRap = v0->Y(-3122);
2284     Double_t fEta = v0->PseudoRapV0();
2285
2286     Double_t fV0mom[3];
2287
2288     fV0mom[0]=v0->MomV0X();
2289     fV0mom[1]=v0->MomV0Y();
2290     fV0mom[2]=v0->MomV0Z();
2291     Double_t fV0TotalMomentum = TMath::Sqrt(fV0mom[0]*fV0mom[0]+fV0mom[1]*fV0mom[1]+fV0mom[2]*fV0mom[2]);
2292
2293     Double_t fV0cosPointAngle = v0->CosPointingAngle(lPrimaryVtxPosition);
2294     lV0Position[0]= v0->DecayVertexV0X();
2295     lV0Position[1]= v0->DecayVertexV0Y();
2296     lV0Position[2]= v0->DecayVertexV0Z();
2297     Double_t fV0DecayLength = v0->DecayLengthV0(lPrimaryVtxPosition);
2298     fV0Radius  = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
2299
2300     fFFHistosIMALaAllEvt->FillFF(trackPt, invMALa, jetPt, incrementJetPt);
2301     //fh1trackPosNCls->Fill(trackPosNcls);
2302     //fh1trackNegNCls->Fill(trackNegNcls);
2303     fh1EtaALa->Fill(fEta);
2304     if(fAnalysisMC){
2305       TList* listmc = fAOD->GetList();
2306       Bool_t mclabelcheck = MCLabelCheck(v0, kAntiLambda, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode);
2307       if(mclabelcheck == kFALSE)continue;
2308       //if(fPhysicalPrimary == kFALSE)continue;
2309       fh3InvMassEtaTrackPtALa->Fill(fEta,invMALa,trackPt);
2310       fh1PtMCALa->Fill(MCPt);
2311     }
2312     fh1V0Eta->Fill(fEta);
2313     fh1V0totMom->Fill(fV0TotalMomentum);
2314     fh1CosPointAngle->Fill(fV0cosPointAngle);
2315     fh1DecayLengthV0->Fill(fV0DecayLength);
2316     fh1V0Radius->Fill(fV0Radius);
2317     fh1DcaV0Daughters->Fill(fDcaV0Daughters);
2318     fh1DcaPosToPrimVertex->Fill(fDcaPosToPrimVertex);
2319     fh1DcaNegToPrimVertex->Fill(fDcaNegToPrimVertex);
2320     fh1trackPosEta->Fill(PosEta);
2321     fh1trackNegEta->Fill(NegEta);
2322   }
2323
2324   //_____no jets events______________________________________________________________________________________________________________________________________
2325
2326   if(nRecJetsCuts == 0){
2327
2328     // std::cout<<"################## found event without any jet!!!!!###################"<<std::endl;
2329
2330
2331   }
2332
2333
2334
2335   //____ fill all jet related histos  ________________________________________________________________________________________________________________________
2336   //##########################jet loop########################################################################################################################
2337
2338   //fill jet histos in general
2339   for(Int_t ij=0; ij<nRecJetsCuts; ++ij){                               // ij is an index running over the list of the reconstructed jets after cuts, all jets in event
2340
2341     AliAODJet* jet = (AliAODJet*) (fJetsRecCuts->At(ij));
2342
2343     Double_t jetPt  = jet->Pt();
2344     Double_t jetEta = jet->Eta();
2345     Double_t jetPhi = jet->Phi();
2346
2347     //if(ij==0){ // loop over leading jets for ij = 0, for ij>= 0 look into all jets
2348
2349     if(ij>=0){//all jets in event
2350
2351       TList* jettracklist = new TList();
2352       Double_t sumPt      = 0.;
2353       Bool_t isBadJet     = kFALSE;
2354       Int_t njetTracks    = 0;
2355
2356       if(GetFFRadius()<=0){
2357         GetJetTracksTrackrefs(jettracklist, jet, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJet);// list of jet tracks from trackrefs
2358       } else {
2359         GetJetTracksPointing(fTracksRecCuts, jettracklist, jet, GetFFRadius(), sumPt, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJet);  // fill list of tracks in cone around jet axis with cone Radius (= 0.4 standard)
2360       }
2361
2362       if(GetFFMinNTracks()>0 && jettracklist->GetSize() <= GetFFMinNTracks()) isBadJet = kTRUE; // reject jets with less tracks than fFFMinNTracks
2363       if(isBadJet) continue; // rejects jets in which no track has a track pt higher than 5 GeV/c (see AddTask macro)
2364
2365
2366       Float_t fJetAreaMin = 0.6*TMath::Pi()*GetFFRadius()*GetFFRadius(); // minimum jet area cut
2367
2368       //std::cout<<"GetFFRadius(): "<<GetFFRadius()<<std::endl;
2369       //std::cout<<"jet->EffectiveAreaCharged()"<<jet->EffectiveAreaCharged()<<std::endl;
2370       //std::cout<<"fJetAreaMin: "<<fJetAreaMin<<std::endl;
2371
2372       if (jet->EffectiveAreaCharged() < fJetAreaMin)continue;// cut on jet area
2373
2374
2375       fh1JetEta->Fill(jetEta);
2376       fh1JetPhi->Fill(jetPhi);
2377       fh2JetEtaPhi->Fill(jetEta,jetPhi);
2378
2379       // printf("pT = %f, eta = %f, phi = %f, leadtr pt = %f\n, ",jetPt,jetEta,jetphi,leadtrack);
2380
2381       for(Int_t it=0; it<jettracklist->GetSize(); ++it){//loop over all particles in jet
2382
2383         AliVParticle* trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));//all tracks in jet cone
2384         if(!trackVP)continue;
2385
2386         Float_t trackPt = trackVP->Pt();//transversal momentum of jet particle
2387         Float_t trackEta = trackVP->Eta();
2388
2389         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2390
2391         fFFHistosRecCuts->FillFF(trackPt, jetPt, incrementJetPt);
2392         if(nK0s>0) fFFHistosRecCutsK0Evt->FillFF(trackPt, jetPt, incrementJetPt);
2393         fh2FFJetTrackEta->Fill(trackEta,jetPt);
2394
2395
2396       }
2397
2398       njetTracks = jettracklist->GetSize();
2399
2400       //____________________________________________________________________________________________________________________
2401       //alternative method to estimate secondary constribution in jet cone (second method you can see below in rec. K0s loop & rec. Lambdas loop & rec. Antilambdas loop)
2402
2403       if(fAnalysisMC){
2404
2405         TList *list = fAOD->GetList();
2406         AliAODMCHeader *mcHeadr=(AliAODMCHeader*)list->FindObject(AliAODMCHeader::StdBranchName());
2407         if(!mcHeadr)continue;
2408
2409         Double_t mcXv=0., mcYv=0., mcZv=0.;//MC primary vertex position
2410
2411         mcXv=mcHeadr->GetVtxX(); mcYv=mcHeadr->GetVtxY(); mcZv=mcHeadr->GetVtxZ(); // position of the MC primary vertex
2412
2413         for(Int_t it=0; it<jettracklist->GetSize(); ++it){//loop over all tracks in the jet
2414
2415           AliVParticle* trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));//track in jet cone
2416           if(!trackVP)continue;
2417           AliAODTrack *tr = dynamic_cast<AliAODTrack*> (trackVP);                   //fetch one jet track from the TList
2418           if(!tr)continue;
2419
2420           //get MC label information
2421           TList *mclist = fAOD->GetList();
2422
2423           //fetch the MC stack
2424           TClonesArray* stackMC = (TClonesArray*)mclist->FindObject(AliAODMCParticle::StdBranchName()); //get MCAOD branch in data
2425           if (!stackMC) {Printf("ERROR: stack not available");}
2426
2427           else {
2428
2429             Int_t trackLabel = TMath::Abs(tr->GetLabel());                       //fetch jet track label in MC stack
2430
2431             AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(stackMC->At(trackLabel));  //fetch MC gen. particle for rec. jet track
2432
2433             if(!part)continue;  //skip non-existing objects
2434
2435
2436             //Bool_t IsPhysicalPrimary = part->IsPhysicalPrimary();//not recommended to check, better use distance between primary vertex and secondary vertex
2437
2438             Float_t fDistPrimaryMax = 0.01;
2439             // Get the distance between production point of the MC mother particle and the primary vertex
2440
2441             Double_t dx = mcXv-part->Xv();//mc primary vertex - mc gen. v0 vertex
2442             Double_t dy = mcYv-part->Yv();
2443             Double_t dz = mcZv-part->Zv();
2444
2445             Float_t fDistPrimary = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
2446             Bool_t fPhysicalPrimary = (fDistPrimary < fDistPrimaryMax);
2447
2448             // std::cout<<"fDistPrimary"<<fDistPrimary<<std::endl;
2449             // std::cout<<"fPhysicalPrimary"<<fPhysicalPrimary<<std::endl;
2450
2451             if(!fPhysicalPrimary)continue;//rejects Kstar and other strong decaying particles from Secondary Contamination
2452
2453             Bool_t isFromStrange = kFALSE;// flag to check whether particle has strange mother
2454
2455             Int_t iMother = part->GetMother(); //get mother MC gen. particle label
2456
2457             if(iMother >= 0){
2458               AliAODMCParticle *partM = dynamic_cast<AliAODMCParticle*>(stackMC->At(iMother)); //fetch mother of MC gen. particle
2459               if(!partM) continue;
2460
2461               Int_t codeM =  TMath::Abs(partM->GetPdgCode());                                 //mothers pdg code
2462
2463               Int_t mfl = Int_t (codeM/ TMath::Power(10, Int_t(TMath::Log10(codeM))));        //asks for first number of mothers pdg code (strange particles always start with 3..)
2464
2465               if  (mfl == 3 && codeM != 3) isFromStrange = kTRUE;
2466             }
2467
2468             //cut on primary particles:
2469
2470
2471
2472
2473             if(isFromStrange == kTRUE){
2474
2475               Double_t trackPt = part->Pt();
2476               Double_t trackEta = part->Eta();
2477               fh3StrContinCone->Fill(jetPt, trackPt, trackEta);//MC gen. particle parameters, but rec. jet pt
2478
2479               }//isFromStrange is kTRUE
2480           } //end else
2481         }//end loop over jet tracks
2482
2483       }// end fAnalysisMC
2484
2485
2486       fh1TrackMultCone->Fill(njetTracks);
2487       fh2TrackMultCone->Fill(njetTracks,jetPt);
2488
2489       // ---- K0s ----
2490
2491       // fQAJetHistosRecCutsLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt() );
2492
2493       for(Int_t it=0; it<fListK0s->GetSize(); ++it){ // loop all K0s
2494
2495         AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListK0s->At(it));
2496         if(!v0) continue;//rejection of events with no V0 vertex
2497
2498         Double_t v0Mom[3];
2499         v0->PxPyPz(v0Mom);
2500         TVector3 v0MomVect(v0Mom);
2501
2502         Double_t dPhiJetK0 = (jet->MomentumVector()->Vect()).DeltaPhi(v0MomVect);
2503         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2504
2505         if(incrementJetPt){fh1V0JetPt->Fill(jetPt);}
2506
2507         Double_t invMK0s =0;
2508         Double_t trackPt=0;
2509         CalculateInvMass(v0, kK0, invMK0s, trackPt);  //function to calculate invMass with TLorentzVector class
2510
2511         fFFHistosIMK0Jet->FillFF(trackPt, invMK0s, jetPt, incrementJetPt);
2512         fFFHistosPhiCorrIMK0->FillPhiCorr(trackPt,TVector2::Phi_0_2pi(dPhiJetK0),invMK0s);  //filled for MC with all particles, no association checks (pdg code, proper decay mode etc..) or matching to MC truth
2513
2514         if(dPhiJetK0<fh1dPhiJetK0->GetXaxis()->GetXmin()) dPhiJetK0 += 2*TMath::Pi();
2515         fh1dPhiJetK0->Fill(dPhiJetK0);
2516
2517       }
2518
2519       if(fListK0s->GetSize() == 0){ // no K0: increment jet pt spectrum
2520
2521         Bool_t incrementJetPt = kTRUE;
2522         fFFHistosIMK0Jet->FillFF(-1, -1, jetPt, incrementJetPt);
2523       }
2524
2525       //____fetch reconstructed K0s in cone around jet axis:_______________________________________________________________________________
2526
2527       TList* jetConeK0list = new TList();
2528
2529       Double_t sumPtK0     = 0.;
2530
2531       Bool_t isBadJetK0    = kFALSE; // dummy, do not use
2532
2533
2534       GetTracksInCone(fListK0s, jetConeK0list, jet, GetFFRadius(), sumPtK0, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetK0); //reconstructed K0s in cone around jet axis
2535
2536       if(fDebug>2)Printf("%s:%d nK0s total: %d, in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nK0s,jetConeK0list->GetEntries(),GetFFRadius());
2537
2538
2539       for(Int_t it=0; it<jetConeK0list->GetSize(); ++it){ // loop for K0s in jet cone
2540
2541         AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetConeK0list->At(it));
2542         if(!v0) continue;
2543
2544         Bool_t   incrementJetPt = (it==0) ? kTRUE : kFALSE;
2545         Double_t invMK0s =0;
2546         Double_t trackPt=0;
2547
2548         CalculateInvMass(v0, kK0, invMK0s, trackPt);  //function to calculate invMass with TLorentzVector class
2549
2550         if(fAnalysisMC){
2551           Double_t jetPtSmear = -1;
2552           SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
2553           if(incrementJetPt == kTRUE){fh1FFIMK0ConeSmear->Fill(jetPtSmear);}                          //fill TH1F for normalization purposes
2554         }
2555
2556         fFFHistosIMK0Cone->FillFF(trackPt, invMK0s, jetPt, incrementJetPt);
2557       }
2558
2559
2560       if(jetConeK0list->GetSize() == 0){ // no K0: increment jet pt spectrum
2561
2562         Bool_t incrementJetPt = kTRUE;
2563         fFFHistosIMK0Cone->FillFF(-1, -1, jetPt, incrementJetPt);
2564         if(fAnalysisMC){
2565           Double_t jetPtSmear = -1;
2566           SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
2567           if(incrementJetPt == kTRUE){fh1FFIMK0ConeSmear->Fill(jetPtSmear);}                          //fill TH1F for normalization purposes
2568         }
2569       }
2570
2571       //fetch particles in perpendicular cone to estimate UE event contribution to particle spectrum
2572       //these perpendicular cone particle spectra serve to subtract the particles in jet cones, that are stemming from the Underlying event, on a statistical basis
2573       //for normalization the common jet pT spectrum is used: fh1FFIMK0Cone, fh1FFIMLaCone and fh1FFIMALaCone
2574
2575       //____fetch reconstructed K0s in cone perpendicular to jet axis:_______________________________________________________________________________
2576
2577       TList* jetPerpConeK0list = new TList();
2578
2579       Double_t sumPerpPtK0     = 0.;
2580
2581       GetTracksInPerpCone(fListK0s, jetPerpConeK0list, jet, GetFFRadius(), sumPerpPtK0); //reconstructed K0s in cone around jet axis
2582
2583       if(fDebug>2)Printf("%s:%d nK0s total: %d, in perp jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nK0s,jetPerpConeK0list->GetEntries(),GetFFRadius());
2584
2585       for(Int_t it=0; it<jetPerpConeK0list->GetSize(); ++it){ // loop for K0s in perpendicular cone
2586
2587         AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetPerpConeK0list->At(it));
2588         if(!v0) continue;
2589
2590         Double_t invMPerpK0s =0;
2591         Double_t trackPt=0;
2592
2593         CalculateInvMass(v0, kK0, invMPerpK0s, trackPt);  //function to calculate invMass with TLorentzVector class
2594
2595         fh3IMK0PerpCone->Fill(jetPt, invMPerpK0s, trackPt); //(x,y,z) //pay attention, this histogram contains the V0 content of both (+/- 90 degrees) perp. cones!!
2596       }
2597
2598
2599       if(jetPerpConeK0list->GetSize() == 0){ // no K0s in jet cone
2600
2601         //Bool_t incrementPerpJetPt = kTRUE;
2602         fh3IMK0PerpCone->Fill(jetPt, -1, -1);
2603       }
2604
2605       // ____ rec K0s in median cluster___________________________________________________________________________________________________________
2606
2607       TList* jetMedianConeK0list = new TList();
2608
2609       AliAODJet* medianCluster = GetMedianCluster();
2610
2611       Double_t sumMedianPtK0     = 0.;
2612
2613       Bool_t isBadJetK0Median    = kFALSE; // dummy, do not use
2614
2615       GetTracksInCone(fListK0s, jetMedianConeK0list, medianCluster, GetFFRadius(), sumMedianPtK0, 0., 0., isBadJetK0Median); //reconstructed K0s in median cone around jet axis
2616       //GetTracksInCone(fListK0s, jetConeK0list, jet, GetFFRadius(), sumPtK0, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetK0); //original use of function
2617
2618       //cut parameters from Fragmentation Function task:
2619       //Float_t fFFMinLTrackPt;   // reject jets with leading track with pt smaller than this value, use GetFFMinLTrackPt()
2620       //Float_t fFFMaxTrackPt;    // reject jetscontaining any track with pt larger than this value, use GetFFMaxTrackPt()
2621
2622       for(Int_t it=0; it<jetMedianConeK0list->GetSize(); ++it){ // loop for K0s in median cone
2623
2624         AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetMedianConeK0list->At(it));
2625         if(!v0) continue;
2626
2627         Double_t invMMedianK0s =0;
2628         Double_t trackPt=0;
2629
2630         CalculateInvMass(v0, kK0, invMMedianK0s, trackPt);  //function to calculate invMass with TLorentzVector class
2631
2632         fh3IMK0MedianCone->Fill(jetPt, invMMedianK0s, trackPt); //(x,y,z)
2633       }
2634
2635       if(jetMedianConeK0list->GetSize() == 0){ // no K0s in median cluster cone
2636
2637         fh3IMK0MedianCone->Fill(jetPt, -1, -1);
2638       }
2639
2640       //_________________________________________________________________________________________________________________________________________
2641
2642       //____fetch reconstructed Lambdas in cone perpendicular to jet axis:__________________________________________________________________________
2643
2644       TList* jetPerpConeLalist = new TList();
2645
2646       Double_t sumPerpPtLa     = 0.;
2647
2648       GetTracksInPerpCone(fListLa, jetPerpConeLalist, jet, GetFFRadius(), sumPerpPtLa); //reconstructed Lambdas in cone around jet axis //pay attention, this histogram contains the V0 content of both (+/- 90 degrees) perp. cones!!
2649
2650       if(fDebug>2)Printf("%s:%d nLa total: %d, in perp jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nLa,jetPerpConeLalist->GetEntries(),GetFFRadius());
2651
2652       for(Int_t it=0; it<jetPerpConeLalist->GetSize(); ++it){ // loop for Lambdas in perpendicular cone
2653
2654         AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetPerpConeLalist->At(it));
2655         if(!v0) continue;
2656
2657         Double_t invMPerpLa =0;
2658         Double_t trackPt=0;
2659
2660         CalculateInvMass(v0, kLambda, invMPerpLa, trackPt);  //function to calculate invMass with TLorentzVector class
2661
2662         fh3IMLaPerpCone->Fill(jetPt, invMPerpLa, trackPt);
2663       }
2664
2665
2666       if(jetPerpConeLalist->GetSize() == 0){ // no Lambdas in jet
2667
2668         fh3IMLaPerpCone->Fill(jetPt, -1, -1);
2669
2670       }
2671
2672       //__________________________________________________________________________________________________________________________________________
2673           // ____ rec Lambdas in median cluster___________________________________________________________________________________________________________
2674
2675       TList* jetMedianConeLalist = new TList();
2676
2677       //AliAODJet* medianCluster = GetMedianCluster(); //already loaded at part for K0s ??
2678
2679       Double_t sumMedianPtLa     = 0.;
2680       Bool_t isBadJetLaMedian    = kFALSE; // dummy, do not use
2681
2682       GetTracksInCone(fListLa, jetMedianConeLalist, medianCluster, GetFFRadius(), sumMedianPtLa, 0, 0, isBadJetLaMedian); //reconstructed Lambdas in median cone around jet axis
2683
2684          //cut parameters from Fragmentation Function task:
2685       //Float_t fFFMinLTrackPt;   // reject jets with leading track with pt smaller than this value, use GetFFMinLTrackPt()
2686       //Float_t fFFMaxTrackPt;    // reject jets containing any track with pt larger than this value, use GetFFMaxTrackPt()
2687
2688       for(Int_t it=0; it<jetMedianConeLalist->GetSize(); ++it){ // loop for Lambdas in perpendicular cone
2689
2690         AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetMedianConeLalist->At(it));
2691         if(!v0) continue;
2692
2693         Double_t invMMedianLa =0;
2694         Double_t trackPt=0;
2695
2696         CalculateInvMass(v0, kLambda, invMMedianLa, trackPt);  //function to calculate invMass with TLorentzVector class
2697
2698         fh3IMLaMedianCone->Fill(jetPt, invMMedianLa, trackPt); //(x,y,z)
2699       }
2700
2701       if(jetMedianConeLalist->GetSize() == 0){ // no Lambdas in median cluster cone
2702
2703         fh3IMLaMedianCone->Fill(jetPt, -1, -1);
2704       }
2705
2706       //_________________________________________________________________________________________________________________________________________
2707
2708
2709       //____fetch reconstructed Antilambdas in cone perpendicular to jet axis:___________________________________________________________________
2710
2711       TList* jetPerpConeALalist = new TList();
2712
2713       Double_t sumPerpPtALa     = 0.;
2714
2715       GetTracksInPerpCone(fListALa, jetPerpConeALalist, jet, GetFFRadius(), sumPerpPtALa); //reconstructed Antilambdas in cone around jet axis //pay attention, this histogram contains the V0 content of both (+/- 90 degrees) perp. cones!!
2716
2717       if(fDebug>2)Printf("%s:%d nALa total: %d, in perp jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nALa,jetPerpConeALalist->GetEntries(),GetFFRadius());
2718
2719       for(Int_t it=0; it<jetPerpConeALalist->GetSize(); ++it){ // loop for ALa in perpendicular cone
2720
2721         AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetPerpConeALalist->At(it));
2722         if(!v0) continue;
2723
2724         Double_t invMPerpALa =0;
2725         Double_t trackPt=0;
2726
2727         CalculateInvMass(v0, kAntiLambda, invMPerpALa, trackPt);  //function to calculate invMass with TLorentzVector class
2728
2729         fh3IMALaPerpCone->Fill(jetPt, invMPerpALa, trackPt);
2730
2731       }
2732
2733
2734       if(jetPerpConeALalist->GetSize() == 0){ // no Antilambda
2735
2736         fh3IMALaPerpCone->Fill(jetPt, -1, -1);
2737
2738       }
2739
2740
2741           // ____ rec Antilambdas in median cluster___________________________________________________________________________________________________________
2742
2743       TList* jetMedianConeALalist = new TList();
2744
2745       //AliAODJet* medianCluster = GetMedianCluster(); //already loaded at part for K0s
2746
2747       Double_t sumMedianPtALa     = 0.;
2748
2749       Bool_t isBadJetALaMedian    = kFALSE; // dummy, do not use
2750
2751       GetTracksInCone(fListALa, jetMedianConeALalist, medianCluster, GetFFRadius(), sumMedianPtALa, 0, 0, isBadJetALaMedian); //reconstructed Antilambdas in median cone around jet axis
2752
2753
2754       //cut parameters from Fragmentation Function task:
2755       //Float_t fFFMinLTrackPt;   // reject jets with leading track with pt smaller than this value, use GetFFMinLTrackPt()
2756       //Float_t fFFMaxTrackPt;    // reject jets containing any track with pt larger than this value, use GetFFMaxTrackPt()
2757
2758       for(Int_t it=0; it<jetMedianConeALalist->GetSize(); ++it){ // loop for Antilambdas in median cluster cone
2759
2760         AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetMedianConeALalist->At(it));
2761         if(!v0) continue;
2762
2763         Double_t invMMedianALa =0;
2764         Double_t trackPt=0;
2765
2766         CalculateInvMass(v0, kAntiLambda, invMMedianALa, trackPt);  //function to calculate invMass with TLorentzVector class
2767
2768         fh3IMALaMedianCone->Fill(jetPt, invMMedianALa, trackPt); //(x,y,z)
2769       }
2770
2771       if(jetMedianConeALalist->GetSize() == 0){ // no Antilambdas in median cluster cone
2772
2773         fh3IMALaMedianCone->Fill(jetPt, -1, -1);
2774       }
2775
2776       //_________________________________________________________________________________________________________________________________________
2777
2778
2779
2780       //MC Analysis
2781       //__________________________________________________________________________________________________________________________________________
2782
2783       if(fAnalysisMC){
2784
2785         //fill feeddown candidates from TList
2786         //std::cout<<"fListFeeddownLaCand entries: "<<fListFeeddownLaCand->GetSize()<<std::endl;
2787
2788         Double_t sumPtFDLa     = 0.;
2789         Bool_t isBadJetFDLa    = kFALSE; // dummy, do not use
2790
2791         GetTracksInCone(fListFeeddownLaCand, jetConeFDLalist, jet, GetFFRadius(), sumPtFDLa, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetFDLa);
2792
2793         Double_t sumPtFDALa     = 0.;
2794         Bool_t isBadJetFDALa    = kFALSE; // dummy, do not use
2795
2796         GetTracksInCone(fListFeeddownALaCand, jetConeFDALalist, jet, GetFFRadius(), sumPtFDALa, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetFDALa);
2797
2798       //_________________________________________________________________
2799         for(Int_t it=0; it<fListFeeddownLaCand->GetSize(); ++it){
2800
2801           AliAODv0* mcfd = dynamic_cast<AliAODv0*>(fListFeeddownLaCand->At(it));
2802           if(!mcfd) continue;
2803
2804           Double_t invMLaFDcand = 0;
2805           Double_t trackPt = 0;//pt of ass. particle, not used for the histos
2806
2807           CalculateInvMass(mcfd, kLambda, invMLaFDcand, trackPt);
2808
2809           //Get MC gen. Lambda transverse momentum
2810           TClonesArray *st = 0x0;
2811
2812           if(!fAOD)continue;
2813           TList *lt = fAOD->GetList();
2814           if(!lt)continue;
2815
2816           st = (TClonesArray*)lt->FindObject(AliAODMCParticle::StdBranchName()); //get MCAOD branch in data
2817           if (!st)continue;
2818
2819           AliAODTrack *daughtertrack = (AliAODTrack *) (mcfd->GetSecondaryVtx()->GetDaughter(0));//fetch the first of the two daughter tracks
2820           Int_t AssLabel = TMath::Abs(daughtertrack->GetLabel());
2821
2822           AliAODMCParticle *mcDaughterPart =(AliAODMCParticle*)st->UncheckedAt(AssLabel);
2823
2824           Int_t v0lab = mcDaughterPart->GetMother();
2825
2826           //  Int_t v0lab= TMath::Abs(mcfd->GetLabel());//GetLabel doesn't work for AliAODv0 class!!! Only for AliAODtrack
2827
2828           if((!v0lab) || (v0lab<0) || (v0lab > st->GetEntriesFast()))continue;//validity checks
2829
2830           AliAODMCParticle *mcp=(AliAODMCParticle*)st->UncheckedAt(v0lab);
2831
2832           Double_t genLaPt = mcp->Pt();
2833
2834           //std::cout<<"Incl FD, genLaPt:"<<genLaPt<<std::endl;
2835
2836           fh3FeedDownLa->Fill(5., invMLaFDcand, genLaPt);
2837
2838         }//end loop over feeddown candidates for Lambda particles in jet cone
2839         //fetch MC truth in jet cones, denominator of rec. efficiency in jet cones
2840         //_________________________________________________________________
2841         for(Int_t it=0; it<jetConeFDLalist->GetSize(); ++it){
2842
2843           AliAODv0* mcfd = dynamic_cast<AliAODv0*>(jetConeFDLalist->At(it));
2844           if(!mcfd) continue;
2845
2846           //std::cout<<"Cone, recLaPt:"<<mcfd->Pt()<<std::endl;
2847
2848           Double_t invMLaFDcand = 0;
2849           Double_t trackPt = mcfd->Pt();//pt of ass. particle, not used for the histos
2850
2851           CalculateInvMass(mcfd, kLambda, invMLaFDcand, trackPt);
2852
2853           //Get MC gen. Lambda transverse momentum
2854           TClonesArray *st = 0x0;
2855
2856           TList *lt = fAOD->GetList();
2857           if(!lt)continue;
2858
2859           st = (TClonesArray*)lt->FindObject(AliAODMCParticle::StdBranchName());
2860
2861           AliAODTrack *daughtertrack = (AliAODTrack *) (mcfd->GetSecondaryVtx()->GetDaughter(0));//fetch the first of the two daughter tracks
2862           Int_t AssLabel = TMath::Abs(daughtertrack->GetLabel());
2863
2864           AliAODMCParticle *mcDaughterPart =(AliAODMCParticle*)st->UncheckedAt(AssLabel);
2865
2866           Int_t v0lab = mcDaughterPart->GetMother();
2867
2868           //std::cout<<"v0lab: "<<v0lab<<std::endl;
2869
2870           if((!v0lab) || (v0lab<0) || (v0lab > st->GetEntriesFast()))continue;//validity checks
2871
2872           AliAODMCParticle *mcp=(AliAODMCParticle*)st->UncheckedAt(v0lab);
2873
2874           Double_t genLaPt = mcp->Pt();
2875
2876
2877           //std::cout<<"Cone FD, genLaPt:"<<genLaPt<<std::endl;
2878
2879           fh3FeedDownLaCone->Fill(jetPt, invMLaFDcand, genLaPt);
2880
2881         }//end loop over feeddown candidates for Lambda particles in jet cone
2882
2883         //_________________________________________________________________
2884         for(Int_t it=0; it<fListFeeddownALaCand->GetSize(); ++it){
2885
2886           AliAODv0* mcfd = dynamic_cast<AliAODv0*>(fListFeeddownALaCand->At(it));
2887           if(!mcfd) continue;
2888
2889           Double_t invMALaFDcand = 0;
2890           Double_t trackPt = 0;//pt of ass. particle, not used for the histos
2891
2892           CalculateInvMass(mcfd, kAntiLambda, invMALaFDcand, trackPt);
2893
2894           //Get MC gen. Antilambda transverse momentum
2895           TClonesArray *st = 0x0;
2896
2897           TList *lt = fAOD->GetList();
2898           if(!lt)continue;
2899
2900           st = (TClonesArray*)lt->FindObject(AliAODMCParticle::StdBranchName());
2901
2902           AliAODTrack *daughtertrack = (AliAODTrack *) (mcfd->GetSecondaryVtx()->GetDaughter(0));//fetch the first of the two daughter tracks
2903           Int_t AssLabel = TMath::Abs(daughtertrack->GetLabel());
2904
2905           AliAODMCParticle *mcDaughterPart =(AliAODMCParticle*)st->UncheckedAt(AssLabel);
2906
2907           Int_t v0lab = mcDaughterPart->GetMother();
2908
2909
2910           if((!v0lab) || (v0lab<0) || (v0lab > st->GetEntriesFast()))continue;//validity checks
2911
2912           AliAODMCParticle *mcp=(AliAODMCParticle*)st->UncheckedAt(v0lab);
2913
2914           Double_t genALaPt = mcp->Pt();
2915
2916           fh3FeedDownALa->Fill(5., invMALaFDcand, genALaPt);
2917
2918         }//end loop over feeddown candidates for Antilambda particles
2919
2920         //_________________________________________________________________
2921         //feeddown for Antilambdas from Xi(bar)+ and Xi(bar)0 in jet cone:
2922
2923         for(Int_t it=0; it<jetConeFDALalist->GetSize(); ++it){
2924
2925           AliAODv0* mcfd = dynamic_cast<AliAODv0*>(jetConeFDALalist->At(it));
2926           if(!mcfd) continue;
2927
2928           Double_t invMALaFDcand = 0;
2929           Double_t trackPt = 0;//pt of ass. particle, not used for the histos
2930
2931           CalculateInvMass(mcfd, kAntiLambda, invMALaFDcand, trackPt);
2932
2933           //Get MC gen. Antilambda transverse momentum
2934           TClonesArray *st = 0x0;
2935
2936           TList *lt = fAOD->GetList();
2937           if(!lt)continue;
2938
2939           st = (TClonesArray*)lt->FindObject(AliAODMCParticle::StdBranchName());
2940
2941           AliAODTrack *daughtertrack = (AliAODTrack *) (mcfd->GetSecondaryVtx()->GetDaughter(0));//fetch the first of the two daughter tracks
2942           Int_t AssLabel = TMath::Abs(daughtertrack->GetLabel());
2943
2944           AliAODMCParticle *mcDaughterPart =(AliAODMCParticle*)st->UncheckedAt(AssLabel);
2945
2946           Int_t v0lab = mcDaughterPart->GetMother();
2947
2948           if((!v0lab) || (v0lab<0) || (v0lab > st->GetEntriesFast()))continue;//validity checks
2949
2950           AliAODMCParticle *mcp=(AliAODMCParticle*)st->UncheckedAt(v0lab);
2951
2952           Double_t genALaPt = mcp->Pt();
2953
2954           fh3FeedDownALaCone->Fill(jetPt, invMALaFDcand, genALaPt);
2955
2956         }//end loop over feeddown candidates for Antilambda particles in jet cone
2957
2958
2959
2960         //____fetch MC generated K0s in cone around jet axis__(note: particles can stem from fragmentation but also from underlying event)________
2961
2962         Double_t sumPtMCgenK0s   = 0.;
2963         Bool_t isBadJetMCgenK0s  = kFALSE; // dummy, do not use
2964
2965
2966         fListMCgenK0sCone = new TList();      //MC generated K0s in (only geometrical) jet cone (these are MC gen K0s falling geometrically into jet cone (R = 0.4) around jet axis, that was found by anti-kt jet finder, particles can stem from fragmentation but also from underlying event!!)
2967
2968         //first: sampling MC gen K0s
2969
2970         GetTracksInCone(fListMCgenK0s, fListMCgenK0sCone, jet, GetFFRadius(), sumPtMCgenK0s, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetMCgenK0s); //MC generated K0s in cone around jet axis
2971
2972         if(fDebug>2)Printf("%s:%d nMCgenK0s in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,fListMCgenK0sCone->GetEntries(),GetFFRadius());
2973
2974
2975         for(Int_t it=0; it<fListMCgenK0sCone->GetSize(); ++it){ // loop MC generated K0s in cone around jet axis
2976
2977           AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenK0sCone->At(it));
2978           if(!mcp0) continue;
2979
2980           //Double_t fRapMCgenK0s   = MyRapidity(mcp0->E(),mcp0->Pz());//get rec. particle in cone information
2981           Double_t fEtaMCgenK0s   = mcp0->Eta();
2982           Double_t fPtMCgenK0s    = mcp0->Pt();
2983
2984           fh2MCgenK0Cone->Fill(jetPt,fPtMCgenK0s);
2985           fh2MCEtagenK0Cone->Fill(jetPt,fEtaMCgenK0s);
2986
2987         }
2988
2989         //check whether the reconstructed K0s in jet cone are stemming from MC gen K0s (on MCgenK0s list):__________________________________________________
2990
2991         for(Int_t ic=0; ic<jetConeK0list->GetSize(); ++ic){    //loop over all reconstructed K0s in jet cone
2992
2993         //for(Int_t ic=0; ic<fListK0s->GetSize(); ++ic){     //loop over all reconstructed K0s -> previous definition of reconstruction efficiency, not sure what is the better one to choose
2994
2995           Int_t negDaughterpdg;
2996           Int_t posDaughterpdg;
2997           Int_t motherType;
2998           Int_t v0Label;
2999           Double_t fPtMCrecK0Match;
3000           Double_t invMK0Match;
3001           Double_t MCPt;
3002           Int_t nnum =-1;
3003           Int_t pnum =-1;
3004           Bool_t fPhysicalPrimary = -1;
3005           Int_t MCv0PDGCode =0;
3006           Double_t jetPtSmear = -1;
3007
3008           AliAODv0* v0c = dynamic_cast<AliAODv0*>(jetConeK0list->At(ic));//pointer to reconstructed K0s inside jet cone (cone is placed around reconstructed jet axis)
3009
3010           //AliAODv0* v0c = dynamic_cast<AliAODv0*>(fListK0s->At(ic));//pointer to reconstructed K0s
3011           if(!v0c) continue;
3012
3013           Bool_t daughtercheck = DaughterTrackCheck(v0c, nnum, pnum);//check daughter tracks have proper sign
3014           if(daughtercheck == kFALSE)continue;
3015
3016           const AliAODTrack *trackMCNeg=(AliAODTrack *)(v0c->GetDaughter(nnum));
3017           const AliAODTrack *trackMCPos=(AliAODTrack *)(v0c->GetDaughter(pnum));
3018
3019           TList *listmc = fAOD->GetList();
3020
3021           Bool_t mclabelcheck = MCLabelCheck(v0c, kK0, trackMCNeg, trackMCPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PDGCode);
3022
3023           if(mclabelcheck == kFALSE)continue;
3024           if(fPhysicalPrimary == kFALSE)continue;  //requirements for rec. V0 associated to MC true primary particle
3025
3026           for(Int_t it=0; it<fListMCgenK0s->GetSize(); ++it){                                    // loop over MC generated K0s in event, check whether associated MC particle is part of it
3027
3028           //for(Int_t it=0; it<fListMCgenK0sCone->GetSize(); ++it){//belongs to previous definition of rec. eff. of V0s within jet cone
3029
3030             //Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
3031             //AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenK0sCone->At(it));
3032             AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenK0s->At(it));
3033             if(!mcp0) continue;
3034
3035             Bool_t particleMatching = IsParticleMatching(mcp0, v0Label);
3036
3037             if(particleMatching == kFALSE)continue;                                              //if reconstructed V0 particle doesn't match to the associated MC particle go to next stack entry
3038             CalculateInvMass(v0c, kK0, invMK0Match, fPtMCrecK0Match);
3039
3040             Double_t fPtMCgenK0s    = mcp0->Pt();
3041
3042             fh3MCrecK0Cone->Fill(jetPt,invMK0Match,fPtMCgenK0s);                                 //fill matching rec. K0s in 3D histogram
3043
3044             SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);           //jetPt, cent, jetRadius, ptmintrack, &jetPtSmear
3045
3046             fh3MCrecK0ConeSmear->Fill(jetPtSmear,invMK0Match,fPtMCgenK0s);    //fill matching rec. K0s in 3D histogram, jet pT smeared according to deltaptjet distribution width
3047
3048
3049           } // end MCgenK0s / MCgenK0sCone loop
3050
3051           //___________
3052           //check the K0s daughters contamination of the jet tracks:
3053
3054           TClonesArray *stackMC = 0x0;
3055
3056           for(Int_t it=0; it<jettracklist->GetSize(); ++it){//loop over all particles in the jet
3057
3058             AliVParticle* trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));//all tracks in jet cone
3059             if(!trackVP)continue;
3060             AliAODTrack *tr = dynamic_cast<AliAODTrack*> (trackVP);                   //fetch one jet track from the TList
3061             if(!tr)continue;
3062
3063             //get MC label information
3064             TList *mclist = fAOD->GetList();                                           //fetch the MC stack
3065             if(!mclist)continue;
3066
3067             stackMC = (TClonesArray*)mclist->FindObject(AliAODMCParticle::StdBranchName()); //get MCAOD branch in data
3068             if (!stackMC) {Printf("ERROR: stack not available");}
3069             else {
3070
3071               Int_t particleLabel = TMath::Abs(tr->GetLabel());                       //fetch jet track label in MC stack
3072               if(!tr)continue;
3073               //v0c is pointer to K0s candidate, is fetched already above, here it is just checked again whether daughters are properly ordered by their charge
3074
3075               Bool_t daughterchecks = DaughterTrackCheck(v0c, nnum, pnum);
3076
3077               if(daughterchecks == kFALSE)continue;                                   //make sure that daughters are properly ordered
3078
3079               const AliAODTrack *trackNeg=(AliAODTrack *)(v0c->GetDaughter(nnum));    //fetch v0 daughters of reconstructed K0s
3080               const AliAODTrack *trackPos=(AliAODTrack *)(v0c->GetDaughter(pnum));
3081
3082               if(!trackNeg)continue;
3083               if(!trackPos)continue;
3084
3085               Int_t negAssLabel = TMath::Abs(trackNeg->GetLabel());                   //negative (reconstructed) charged track label in MC stack
3086               Int_t posAssLabel = TMath::Abs(trackPos->GetLabel());                   //positive (reconstructed) charged track label in MC stack
3087
3088
3089               if(particleLabel == posAssLabel){                                       //check whether jet track and each of the rec. K0s daughters have same MC stack label -> are identical
3090                 AliAODMCParticle* mctrackPos = dynamic_cast<AliAODMCParticle*>(stackMC->At(posAssLabel));
3091                 if(!mctrackPos) continue;
3092                 Double_t trackPosPt = mctrackPos->Pt();
3093                 Double_t trackPosEta = mctrackPos->Eta();
3094                 fh3SecContinCone->Fill(jetPt, trackPosPt, trackPosEta);}              //if it's the case, fill jet pt, daughter track pt and track eta in histo
3095
3096               if(particleLabel == negAssLabel){
3097                 AliAODMCParticle* mctrackNeg = dynamic_cast<AliAODMCParticle*>(stackMC->At(negAssLabel));
3098                 if(!mctrackNeg) continue;
3099                 Double_t trackNegPt = mctrackNeg->Pt();
3100                 Double_t trackNegEta = mctrackNeg->Eta();
3101                 fh3SecContinCone->Fill(jetPt, trackNegPt, trackNegEta);}              //if it's the case, fill jet pt, daughter track pt and track eta in histo
3102             }
3103           }
3104
3105
3106           //_______________
3107
3108
3109         } //end rec-K0-in-cone loop
3110
3111         //________________________________________________________________________________________________________________________________________________________
3112
3113
3114
3115         delete fListMCgenK0sCone;
3116
3117
3118       }//end fAnalysisMC
3119
3120       delete jetConeK0list;
3121       delete jetPerpConeK0list;
3122       delete jetPerpConeLalist;
3123       delete jetPerpConeALalist;
3124       delete jetMedianConeK0list;
3125       delete jetMedianConeLalist;
3126       delete jetMedianConeALalist;
3127
3128       //---------------La--------------------------------------------------------------------------------------------------------------------------------------------
3129
3130       // fQAJetHistosRecCutsLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt() );
3131
3132       for(Int_t it=0; it<fListLa->GetSize(); ++it){ // loop all La
3133
3134         AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListLa->At(it));
3135         if(!v0) continue;
3136
3137         Double_t v0Mom[3];
3138         v0->PxPyPz(v0Mom);
3139         TVector3 v0MomVect(v0Mom);
3140
3141         Double_t dPhiJetLa = (jet->MomentumVector()->Vect()).DeltaPhi(v0MomVect);
3142
3143         Double_t invMLa =0;
3144         Double_t trackPt=0;
3145
3146         CalculateInvMass(v0, kLambda, invMLa, trackPt); //function to calculate invMass with TLorentzVector class
3147         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
3148
3149         //if(incrementJetPt){fh1V0JetPt->Fill(jetPt);}
3150
3151         fFFHistosIMLaJet->FillFF(trackPt, invMLa, jetPt, incrementJetPt);
3152         fFFHistosPhiCorrIMLa->FillPhiCorr(trackPt,TVector2::Phi_0_2pi(dPhiJetLa),invMLa);
3153
3154         if(dPhiJetLa<fh1dPhiJetLa->GetXaxis()->GetXmin()) dPhiJetLa += 2*TMath::Pi();
3155         fh1dPhiJetLa->Fill(dPhiJetLa);
3156       }
3157
3158       if(fListLa->GetSize() == 0){ // no La: increment jet pt spectrum
3159
3160         Bool_t incrementJetPt = kTRUE;
3161         fFFHistosIMLaJet->FillFF(-1, -1, jetPt, incrementJetPt);
3162       }
3163
3164
3165       // ____fetch rec. Lambdas in cone around jet axis_______________________________________________________________________________________
3166
3167       TList* jetConeLalist = new TList();
3168       Double_t sumPtLa     = 0.;
3169       Bool_t isBadJetLa    = kFALSE; // dummy, do not use
3170
3171       GetTracksInCone(fListLa, jetConeLalist, jet, GetFFRadius(), sumPtLa, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetLa);//method inherited from FF
3172
3173       if(fDebug>2)Printf("%s:%d nLa total: %d, in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nLa,jetConeLalist->GetEntries(),GetFFRadius());
3174
3175       for(Int_t it=0; it<jetConeLalist->GetSize(); ++it){ // loop La in jet cone
3176
3177         AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetConeLalist->At(it));
3178         if(!v0) continue;
3179         Double_t invMLa =0;
3180         Double_t trackPt=0;
3181
3182         CalculateInvMass(v0, kLambda, invMLa, trackPt); //function to calculate invMass with TLorentzVector class
3183
3184         Bool_t   incrementJetPt = (it==0) ? kTRUE : kFALSE;
3185
3186         if(fAnalysisMC){
3187           Double_t jetPtSmear = -1;
3188           SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
3189           if(incrementJetPt == kTRUE){fh1FFIMLaConeSmear->Fill(jetPtSmear);}                          //fill TH1F for normalization purposes
3190         }
3191
3192         fFFHistosIMLaCone->FillFF(trackPt, invMLa, jetPt, incrementJetPt);
3193       }
3194
3195       if(jetConeLalist->GetSize() == 0){ // no La: increment jet pt spectrum
3196
3197         Bool_t incrementJetPt = kTRUE;
3198         fFFHistosIMLaCone->FillFF(-1, -1, jetPt, incrementJetPt);
3199
3200         if(fAnalysisMC){
3201           Double_t jetPtSmear;
3202           SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
3203           if(incrementJetPt == kTRUE)fh1FFIMLaConeSmear->Fill(jetPtSmear);}
3204
3205       }
3206
3207       if(fAnalysisMC){
3208
3209         //____fetch MC generated Lambdas in cone around jet axis__(particles can stem from fragmentation but also from underlying event)_____________
3210
3211         Double_t sumPtMCgenLa      = 0.;
3212         Bool_t isBadJetMCgenLa  = kFALSE; // dummy, do not use
3213
3214         //sampling MC gen. Lambdas in cone around reconstructed jet axis
3215         fListMCgenLaCone = new TList();
3216
3217         GetTracksInCone(fListMCgenLa, fListMCgenLaCone, jet, GetFFRadius(), sumPtMCgenLa, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetMCgenLa);//fetch MC generated Lambdas in cone of resolution parameter R around jet axis
3218
3219         if(fDebug>2)Printf("%s:%d nMCgenLa in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,fListMCgenLaCone->GetEntries(),GetFFRadius());
3220
3221         for(Int_t it=0; it<fListMCgenLaCone->GetSize(); ++it){ // loop MC generated La in cone around jet axis
3222
3223           AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenLaCone->At(it));
3224           if(!mcp0) continue;
3225
3226           //Double_t fRapMCgenLa   = MyRapidity(mcp0->E(),mcp0->Pz());
3227           Double_t fEtaMCgenLa   = mcp0->Eta();
3228           Double_t fPtMCgenLa    = mcp0->Pt();
3229
3230           fh2MCgenLaCone->Fill(jetPt,fPtMCgenLa);
3231           fh2MCEtagenLaCone->Fill(jetPt,fEtaMCgenLa);
3232         }
3233
3234
3235         //check whether the reconstructed La are stemming from MC gen La on fListMCgenLa List:__________________________________________________
3236
3237         for(Int_t ic=0; ic<jetConeLalist->GetSize(); ++ic){//loop over all reconstructed La within jet cone, new definition
3238
3239           //for(Int_t ic=0; ic<fListLa->GetSize(); ++ic){//old definition
3240
3241           Int_t negDaughterpdg;
3242           Int_t posDaughterpdg;
3243           Int_t motherType;
3244           Int_t v0Label;
3245           Double_t fPtMCrecLaMatch;
3246           Double_t invMLaMatch;
3247           Double_t MCPt;
3248           Int_t nnum;
3249           Int_t pnum;
3250           Bool_t fPhysicalPrimary = -1;
3251           Int_t MCv0PDGCode =0;
3252           Double_t jetPtSmear = -1;
3253
3254           AliAODv0* v0c = dynamic_cast<AliAODv0*>(jetConeLalist->At(ic));//new definition
3255
3256           //AliAODv0* v0c = dynamic_cast<AliAODv0*>(fListLa->At(ic));//old definition
3257           if(!v0c) continue;
3258
3259           Bool_t daughtercheck = DaughterTrackCheck(v0c, nnum, pnum);
3260           if(daughtercheck == kFALSE)continue;
3261
3262           const AliAODTrack *trackMCNeg=(AliAODTrack *)(v0c->GetDaughter(nnum));
3263           const AliAODTrack *trackMCPos=(AliAODTrack *)(v0c->GetDaughter(pnum));
3264
3265           TList *listmc = fAOD->GetList();
3266
3267           Bool_t mclabelcheck = MCLabelCheck(v0c, kLambda, trackMCNeg, trackMCPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PDGCode);
3268
3269           if(mclabelcheck == kFALSE)continue;
3270           if(fPhysicalPrimary == kFALSE)continue;
3271
3272           for(Int_t it=0; it<fListMCgenLa->GetSize(); ++it){//new definition                                  // loop over MC generated K0s in cone around jet axis
3273
3274             // for(Int_t it=0; it<fListMCgenLaCone->GetSize(); ++it){//old definition                                  // loop over MC generated La in cone around jet axis
3275
3276             //Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
3277
3278             AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenLa->At(it));//new definition
3279             //AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenLaCone->At(it));//old definition
3280
3281             if(!mcp0) continue;
3282
3283             Bool_t particleMatching = IsParticleMatching(mcp0, v0Label);
3284
3285
3286             if(particleMatching == kFALSE)continue; //particle doesn't match on any associated MC gen particle in cone around rec jet axis
3287
3288             CalculateInvMass(v0c, kLambda, invMLaMatch, fPtMCrecLaMatch);
3289
3290             Double_t fPtMCgenLa    = mcp0->Pt();
3291
3292             fh3MCrecLaCone->Fill(jetPt,invMLaMatch,fPtMCgenLa);                        //fill matching rec. K0s 3D histogram
3293
3294             SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
3295
3296             fh3MCrecLaConeSmear->Fill(jetPtSmear,invMLaMatch,fPtMCgenLa);              //fill matching rec. Lambdas in 3D histogram, jet pT smeared according to deltaptjet distribution width
3297
3298
3299           } // end MCgenLa loop
3300
3301             //check the Lambda daughters contamination of the jet tracks://///////////////////////////////////////////////////////////////////////////////////////////
3302
3303           TClonesArray *stackMC = 0x0;
3304
3305           for(Int_t it=0; it<jettracklist->GetSize(); ++it){//loop over all particles in the jet
3306
3307             AliVParticle* trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));//all tracks in jet cone
3308             if(!trackVP)continue;
3309             AliAODTrack *tr = dynamic_cast<AliAODTrack*> (trackVP);                   //fetch one jet track from the TList
3310             if(!tr)continue;
3311
3312             //get MC label information
3313             TList *mclist = fAOD->GetList();                                           //fetch the MC stack
3314
3315             stackMC = (TClonesArray*)mclist->FindObject(AliAODMCParticle::StdBranchName()); //get MCAOD branch in data
3316             if (!stackMC) {Printf("ERROR: stack not available");}
3317             else {
3318
3319               Int_t particleLabel = TMath::Abs(tr->GetLabel());                       //fetch jet track label in MC stack
3320
3321               Bool_t daughterchecks = DaughterTrackCheck(v0c, nnum, pnum);
3322
3323               if(daughterchecks == kFALSE)continue;                                   //make sure that daughters are properly ordered
3324
3325               const AliAODTrack *trackNeg=(AliAODTrack *)(v0c->GetDaughter(nnum));    //fetch v0 daughters of reconstructed K0s
3326               const AliAODTrack *trackPos=(AliAODTrack *)(v0c->GetDaughter(pnum));
3327
3328               Int_t negAssLabel = TMath::Abs(trackNeg->GetLabel());                   //negative (reconstructed) charged track label in MC stack
3329               Int_t posAssLabel = TMath::Abs(trackPos->GetLabel());                   //positive (reconstructed) charged track label in MC stack
3330
3331
3332               if(particleLabel == posAssLabel){                                       //check whether jet track and each of the rec. K0s daughters have same MC stack label -> are identical
3333
3334                 AliAODMCParticle* mctrackPos = dynamic_cast<AliAODMCParticle*>(stackMC->At(posAssLabel));
3335                 if(!mctrackPos) continue;
3336
3337                 Double_t trackPosPt = trackPos->Pt();
3338                 Double_t trackPosEta = trackPos->Eta();
3339                 fh3SecContinCone->Fill(jetPt, trackPosPt, trackPosEta);}              //if it's the case, fill jet pt, daughter track pt and track eta in histo
3340
3341               if(particleLabel == negAssLabel){
3342
3343                 AliAODMCParticle* mctrackNeg = dynamic_cast<AliAODMCParticle*>(stackMC->At(negAssLabel));
3344                 if(!mctrackNeg) continue;
3345
3346                 Double_t trackNegPt = trackNeg->Pt();
3347                 Double_t trackNegEta = trackNeg->Eta();
3348                 fh3SecContinCone->Fill(jetPt, trackNegPt, trackNegEta);}              //if it's the case, fill jet pt, daughter track pt and track eta in histo
3349             }
3350           }
3351
3352
3353
3354         } //end rec-La-in-cone loop
3355         //________________________________________________________________________________________________________________________________________________________
3356
3357         delete fListMCgenLaCone;
3358
3359       }//end fAnalysisMC
3360
3361       delete jetConeLalist;
3362
3363
3364
3365       //---------------ALa-----------
3366
3367
3368       // fQAJetHistosRecCutsLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt() );
3369
3370       for(Int_t it=0; it<fListALa->GetSize(); ++it){ // loop all ALa
3371
3372         AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListALa->At(it));
3373         if(!v0) continue;
3374
3375         Double_t v0Mom[3];
3376         v0->PxPyPz(v0Mom);
3377         TVector3 v0MomVect(v0Mom);
3378
3379         Double_t dPhiJetALa = (jet->MomentumVector()->Vect()).DeltaPhi(v0MomVect);
3380
3381         Double_t invMALa =0;
3382         Double_t trackPt=0;
3383
3384         CalculateInvMass(v0, kAntiLambda, invMALa, trackPt); //function to calculate invMass with TLorentzVector class
3385         Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
3386
3387         //if(incrementJetPt){fh1V0JetPt->Fill(jetPt);}
3388
3389         fFFHistosIMALaJet->FillFF(trackPt, invMALa, jetPt, incrementJetPt);
3390         fFFHistosPhiCorrIMALa->FillPhiCorr(trackPt,TVector2::Phi_0_2pi(dPhiJetALa),invMALa);
3391
3392         if(dPhiJetALa<fh1dPhiJetALa->GetXaxis()->GetXmin()) dPhiJetALa += 2*TMath::Pi();
3393         fh1dPhiJetALa->Fill(dPhiJetALa);
3394       }
3395
3396       if(fListALa->GetSize() == 0){ // no ALa: increment jet pt spectrum
3397
3398         Bool_t incrementJetPt = kTRUE;
3399         fFFHistosIMALaJet->FillFF(-1, -1, jetPt, incrementJetPt);
3400       }
3401
3402
3403       // ____fetch rec. Antilambdas in cone around jet axis_______________________________________________________________________________________
3404
3405       TList* jetConeALalist = new TList();
3406       Double_t sumPtALa     = 0.;
3407       Bool_t isBadJetALa    = kFALSE; // dummy, do not use
3408
3409       GetTracksInCone(fListALa, jetConeALalist, jet, GetFFRadius(), sumPtALa, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetALa);//method inherited from FF
3410
3411       if(fDebug>2)Printf("%s:%d nALa total: %d, in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nALa,jetConeALalist->GetEntries(),GetFFRadius());
3412
3413       for(Int_t it=0; it<jetConeALalist->GetSize(); ++it){ // loop ALa in jet cone
3414
3415         AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetConeALalist->At(it));
3416         if(!v0) continue;
3417         Double_t invMALa =0;
3418         Double_t trackPt=0;
3419
3420         CalculateInvMass(v0, kAntiLambda, invMALa, trackPt); //function to calculate invMass with TLorentzVector class
3421
3422         Bool_t   incrementJetPt = (it==0) ? kTRUE : kFALSE;
3423
3424         if(fAnalysisMC){    //jet pt smearing study for Antilambdas
3425           Double_t jetPtSmear = -1;
3426           SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
3427           if(incrementJetPt == kTRUE){fh1FFIMALaConeSmear->Fill(jetPtSmear);}                          //fill TH1F for normalization purposes
3428         }
3429
3430         fFFHistosIMALaCone->FillFF(trackPt, invMALa, jetPt, incrementJetPt);
3431       }
3432
3433       if(jetConeALalist->GetSize() == 0){ // no ALa: increment jet pt spectrum
3434
3435         Bool_t incrementJetPt = kTRUE;
3436         fFFHistosIMALaCone->FillFF(-1, -1, jetPt, incrementJetPt);
3437
3438         if(fAnalysisMC){
3439           Double_t jetPtSmear;
3440           SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
3441           if(incrementJetPt == kTRUE)fh1FFIMALaConeSmear->Fill(jetPtSmear);}
3442
3443       }
3444
3445       if(fAnalysisMC){
3446
3447         //____fetch MC generated Antilambdas in cone around jet axis__(particles can stem from fragmentation but also from underlying event)_____________
3448
3449         Double_t sumPtMCgenALa      = 0.;
3450         Bool_t isBadJetMCgenALa  = kFALSE; // dummy, do not use
3451
3452         //sampling MC gen Antilambdas in cone around reconstructed jet axis
3453         fListMCgenALaCone = new TList();
3454
3455         GetTracksInCone(fListMCgenALa, fListMCgenALaCone, jet, GetFFRadius(), sumPtMCgenALa, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetMCgenALa);//MC generated K0s in cone around jet axis
3456
3457         if(fDebug>2)Printf("%s:%d nMCgenALa in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,fListMCgenALaCone->GetEntries(),GetFFRadius());
3458
3459         for(Int_t it=0; it<fListMCgenALaCone->GetSize(); ++it){ // loop MC generated La in cone around jet axis
3460
3461           AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenALaCone->At(it));
3462           if(!mcp0) continue;
3463
3464           //Double_t fRapMCgenALa   = MyRapidity(mcp0->E(),mcp0->Pz());
3465           Double_t fEtaMCgenALa   = mcp0->Eta();
3466           Double_t fPtMCgenALa    = mcp0->Pt();
3467
3468           fh2MCgenALaCone->Fill(jetPt,fPtMCgenALa);
3469           fh2MCEtagenALaCone->Fill(jetPt,fEtaMCgenALa);
3470         }
3471
3472
3473         //check whether the reconstructed ALa are stemming from MC gen ALa on MCgenALa List:__________________________________________________
3474
3475         for(Int_t ic=0; ic<jetConeALalist->GetSize(); ++ic){//loop over all reconstructed ALa
3476
3477           Int_t negDaughterpdg;
3478           Int_t posDaughterpdg;
3479           Int_t motherType;
3480           Int_t v0Label;
3481           Double_t fPtMCrecALaMatch;
3482           Double_t invMALaMatch;
3483           Double_t MCPt;
3484           Int_t nnum;
3485           Int_t pnum;
3486           Bool_t fPhysicalPrimary = -1;
3487           Int_t MCv0PDGCode =0;
3488           Double_t jetPtSmear = -1;
3489
3490           AliAODv0* v0c = dynamic_cast<AliAODv0*>(jetConeALalist->At(ic));
3491           if(!v0c) continue;
3492
3493           Bool_t daughtercheck = DaughterTrackCheck(v0c, nnum, pnum);
3494           if(daughtercheck == kFALSE)continue;
3495
3496           const AliAODTrack *trackMCNeg=(AliAODTrack *)(v0c->GetDaughter(nnum));
3497           const AliAODTrack *trackMCPos=(AliAODTrack *)(v0c->GetDaughter(pnum));
3498
3499           TList *listmc = fAOD->GetList();
3500           if(!listmc)continue;
3501
3502           Bool_t mclabelcheck = MCLabelCheck(v0c, kAntiLambda, trackMCNeg, trackMCPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PDGCode);
3503
3504           if(mclabelcheck == kFALSE)continue;
3505           if(fPhysicalPrimary == kFALSE)continue;
3506
3507           for(Int_t it=0; it<fListMCgenALa->GetSize(); ++it){                                  // loop over MC generated Antilambdas in cone around jet axis
3508
3509             //Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
3510
3511             AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenALa->At(it));
3512             if(!mcp0) continue;
3513
3514             Bool_t particleMatching = IsParticleMatching(mcp0, v0Label);
3515
3516             if(particleMatching == kFALSE)continue; //particle doesn't match on any associated MC gen particle in cone around rec jet axis
3517
3518             CalculateInvMass(v0c, kAntiLambda, invMALaMatch, fPtMCrecALaMatch);
3519
3520             Double_t fPtMCgenALa  = mcp0->Pt();
3521
3522             fh3MCrecALaCone->Fill(jetPt,invMALaMatch,fPtMCgenALa);                          //fill matching rec. K0s 3D histogram
3523
3524             SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
3525
3526             fh3MCrecALaConeSmear->Fill(jetPtSmear,invMALaMatch,fPtMCgenALa);
3527
3528
3529           } // end MCgenALa loop
3530
3531
3532            //___________
3533           //check the Antilambda daughters contamination of the jet tracks:
3534
3535           TClonesArray *stackMC = 0x0;
3536
3537           for(Int_t it=0; it<jettracklist->GetSize(); ++it){//loop over all particles in the jet
3538
3539             AliVParticle* trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));//all tracks in jet cone
3540             if(!trackVP)continue;
3541             AliAODTrack *tr = dynamic_cast<AliAODTrack*> (trackVP);                   //fetch one jet track from the TList
3542             if(!tr)continue;
3543
3544             //get MC label information
3545             TList *mclist = fAOD->GetList();                                           //fetch the MC stack
3546             if(!mclist)continue;
3547
3548             stackMC = (TClonesArray*)mclist->FindObject(AliAODMCParticle::StdBranchName()); //get MCAOD branch in data
3549             if (!stackMC) {Printf("ERROR: stack not available");}
3550             else {
3551
3552               Int_t particleLabel = TMath::Abs(tr->GetLabel());                       //fetch jet track label in MC stack
3553
3554               Bool_t daughterchecks = DaughterTrackCheck(v0c, nnum, pnum);
3555
3556               if(daughterchecks == kFALSE)continue;                                   //make sure that daughters are properly ordered
3557
3558               const AliAODTrack *trackNeg=(AliAODTrack *)(v0c->GetDaughter(nnum));    //fetch v0 daughters of reconstructed K0s
3559               const AliAODTrack *trackPos=(AliAODTrack *)(v0c->GetDaughter(pnum));
3560               if(!trackPos)continue;
3561               if(!trackNeg)continue;
3562
3563               Int_t negAssLabel = TMath::Abs(trackNeg->GetLabel());                   //negative (reconstructed) charged track label in MC stack
3564               Int_t posAssLabel = TMath::Abs(trackPos->GetLabel());                   //positive (reconstructed) charged track label in MC stack
3565
3566               if(!negAssLabel)continue;
3567               if(!posAssLabel)continue;
3568
3569               if(particleLabel == posAssLabel){                                       //check whether jet track and each of the rec. K0s daughters have same MC stack label -> are identical
3570                 AliAODMCParticle* mctrackPos = dynamic_cast<AliAODMCParticle*>(stackMC->At(posAssLabel));
3571                 if(!mctrackPos) continue;
3572
3573                 Double_t trackPosPt = trackPos->Pt();
3574                 Double_t trackPosEta = trackPos->Eta();
3575                 if(!trackPosPt)continue;
3576                 if(!trackPosEta)continue;
3577
3578                 fh3SecContinCone->Fill(jetPt, trackPosPt, trackPosEta);}              //if it's the case, fill jet pt, daughter track pt and track eta in histo
3579
3580               if(particleLabel == negAssLabel){
3581
3582                 AliAODMCParticle* mctrackNeg = dynamic_cast<AliAODMCParticle*>(stackMC->At(negAssLabel));
3583                 if(!mctrackNeg) continue;
3584
3585                 Double_t trackNegPt = trackNeg->Pt();
3586                 Double_t trackNegEta = trackNeg->Eta();
3587
3588                 if(!trackNegPt)continue;
3589                 if(!trackNegEta)continue;
3590
3591                 fh3SecContinCone->Fill(jetPt, trackNegPt, trackNegEta);}              //if it's the case, fill jet pt, daughter track pt and track eta in histo
3592             }
3593           }
3594
3595         } //end rec-ALa-in-cone loop
3596         //________________________________________________________________________________________________________________________________________________________
3597
3598         delete fListMCgenALaCone;
3599
3600       }//end fAnalysisMC
3601
3602       delete jetConeALalist;
3603       delete jettracklist; //had been initialised at jet loop beginning
3604
3605
3606       }//end of if 'leading' or 'all jet' requirement
3607   }//end of jet loop
3608
3609
3610
3611
3612   fTracksRecCuts->Clear();
3613   fJetsRecCuts->Clear();
3614   fBckgJetsRec->Clear();
3615   fListK0s->Clear();
3616   fListLa->Clear();
3617   fListALa->Clear();
3618   fListFeeddownLaCand->Clear();
3619   fListFeeddownALaCand->Clear();
3620   jetConeFDLalist->Clear();
3621   jetConeFDALalist->Clear();
3622   fListMCgenK0s->Clear();
3623   fListMCgenLa->Clear();
3624   fListMCgenALa->Clear();
3625
3626   //Post output data.
3627   PostData(1, fCommonHistList);
3628 }
3629
3630 // ____________________________________________________________________________________________
3631 void AliAnalysisTaskJetChem::SetProperties(TH3F* h,const char* x, const char* y, const char* z)
3632 {
3633   //Set properties of histos (x,y and z title)
3634
3635   h->SetXTitle(x);
3636   h->SetYTitle(y);
3637   h->SetZTitle(z);
3638   h->GetXaxis()->SetTitleColor(1);
3639   h->GetYaxis()->SetTitleColor(1);
3640   h->GetZaxis()->SetTitleColor(1);
3641 }
3642
3643
3644 //________________________________________________________________________________________________________________________________________
3645 Bool_t AliAnalysisTaskJetChem::AcceptBetheBloch(AliAODv0 *v0, AliPIDResponse *PIDResponse, const Int_t particletype) //dont use for MC Analysis
3646 {
3647
3648         Int_t nnum = 1;
3649         Int_t pnum = 0;
3650         //---
3651         const AliAODTrack *ntracktest=(AliAODTrack *)v0->GetDaughter(nnum);
3652         if(ntracktest->Charge() > 0){nnum = 0; pnum = 1;}
3653
3654         const AliAODTrack *trackNeg=(AliAODTrack *)(v0->GetDaughter(nnum));
3655         const AliAODTrack *trackPos=(AliAODTrack *)(v0->GetDaughter(pnum));
3656
3657         //Check if both tracks are available
3658         if (!trackPos || !trackNeg) {
3659           Printf("strange analysis::UserExec:: Error:Could not retrieve one of the daughter tracks\n");
3660           return kFALSE;
3661         }
3662
3663         //remove like sign V0s
3664         if ( trackPos->Charge() == trackNeg->Charge() ){
3665           //if(fDebug>1) Printf("%s:%d found like-sign V0", (char*)__FILE__,__LINE__);
3666           return kFALSE;
3667           }
3668         //--
3669
3670         Double_t nsig_p = 0; //number of sigmas that positive daughter track has got in TPC pid information
3671         Double_t nsig_n = 0;
3672
3673         const AliAODPid *pid_p=trackPos->GetDetPid();  // returns fDetPID, more detailed or detector specific pid information
3674         const AliAODPid *pid_n=trackNeg->GetDetPid();
3675
3676         if(!pid_p)return kFALSE;
3677         if(!pid_n)return kFALSE;
3678
3679         if (pid_p)
3680           {
3681             if(particletype == 1) //PID cut on positive charged Lambda daughters (only those with pt < 1 GeV/c)
3682               {
3683
3684                 nsig_p=PIDResponse->NumberOfSigmasTPC(trackPos,AliPID::kProton);
3685                 Double_t protonPt = trackPos->Pt();
3686                 if ((TMath::Abs(nsig_p) >= fCutBetheBloch) && (fCutBetheBloch >0) && (protonPt < 1)) return kFALSE;
3687
3688               }
3689
3690
3691           }
3692
3693         if (pid_n)
3694           {
3695             if(particletype == 2)
3696               {
3697                 nsig_n=PIDResponse->NumberOfSigmasTPC(trackNeg,AliPID::kProton);
3698                 Double_t antiprotonPt = trackNeg->Pt();
3699                 if ((TMath::Abs(nsig_n) >= fCutBetheBloch) && (fCutBetheBloch >0) && (antiprotonPt < 1)) return kFALSE;
3700               }
3701
3702           }
3703
3704         return kTRUE;
3705 }
3706
3707 //___________________________________________________________________
3708 Bool_t AliAnalysisTaskJetChem::IsK0InvMass(const Double_t mass) const
3709 {
3710   // K0 mass ? Use FF histo limits
3711
3712   if(fFFIMInvMMin <= mass && mass < fFFIMInvMMax) return kTRUE;
3713
3714   return kFALSE;
3715 }
3716 //___________________________________________________________________
3717 Bool_t AliAnalysisTaskJetChem::IsLaInvMass(const Double_t mass) const
3718 {
3719   // La mass ? Use FF histo limits
3720
3721
3722   if(fFFIMLaInvMMin <= mass && mass < fFFIMLaInvMMax) return kTRUE;
3723
3724   return kFALSE;
3725 }
3726
3727 //_____________________________________________________________________________________
3728 Int_t AliAnalysisTaskJetChem::GetListOfV0s(TList *list, const Int_t type, const Int_t particletype, AliAODVertex* primVertex, AliAODEvent* aod)
3729 {
3730   // fill list of V0s selected according to type
3731
3732   if(!list){
3733     if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
3734     return -1;
3735   }
3736
3737
3738   if(fDebug>5){std::cout<<"AliAnalysisTaskJetChem::GetListOfV0s(): type: "<<type<<" particletype: "<<particletype<<"aod: "<<aod<<std::endl;
3739     if(type==kTrackUndef){std::cout<<"AliAnalysisTaskJetChem::GetListOfV0s(): kTrackUndef!! "<<std::endl;}
3740   }
3741
3742
3743   if(type==kTrackUndef) return 0;
3744
3745   if(!primVertex) return 0;
3746
3747   Double_t lPrimaryVtxPosition[3];
3748   Double_t lV0Position[3];
3749   lPrimaryVtxPosition[0] = primVertex->GetX();
3750   lPrimaryVtxPosition[1] = primVertex->GetY();
3751   lPrimaryVtxPosition[2] = primVertex->GetZ();
3752
3753   if(fDebug>5){ std::cout<<"AliAnalysisTaskJetChem::GetListOfV0s(): aod->GetNumberOfV0s: "<<aod->GetNumberOfV0s()<<std::endl; }
3754
3755
3756   for(int i=0; i<aod->GetNumberOfV0s(); i++){ // loop over V0s
3757
3758
3759     AliAODv0* v0 = aod->GetV0(i);
3760
3761     if(!v0)
3762       {
3763         std::cout << std::endl
3764                   << "Warning in AliAnalysisTaskJetChem::GetListOfV0s:" << std::endl
3765                   << "v0 = " << v0 << std::endl;
3766         continue;
3767       }
3768
3769     Bool_t isOnFly = v0->GetOnFlyStatus();
3770
3771     if(!isOnFly &&  (type == kOnFly || type == kOnFlyPID || type == kOnFlydEdx || type == kOnFlyPrim)) continue;
3772     if( isOnFly &&  (type == kOffl  || type == kOfflPID  || type == kOffldEdx  || type == kOfflPrim))  continue;
3773
3774     Int_t motherType = -1;
3775      //Double_t v0CalcMass = 0;   //mass of MC v0
3776     Double_t MCPt = 0;         //pt of MC v0
3777
3778     Double_t pp[3]={0,0,0}; //3-momentum positive charged track
3779     Double_t pm[3]={0,0,0}; //3-momentum negative charged track
3780     Double_t v0mom[3]={0,0,0};
3781
3782     Double_t invM = 0;
3783     Double_t invMK0s=0;
3784     Double_t invMLa=0;
3785     Double_t invMALa=0;
3786     Double_t trackPt=0;
3787     Int_t nnum = -1;
3788     Int_t pnum = -1;
3789
3790
3791     Bool_t daughtercheck = DaughterTrackCheck(v0, nnum, pnum);
3792
3793     if(daughtercheck == kFALSE)continue;
3794
3795     const AliAODTrack *trackNeg=(AliAODTrack *)(v0->GetDaughter(nnum));
3796     const AliAODTrack *trackPos=(AliAODTrack *)(v0->GetDaughter(pnum));
3797
3798
3799      ///////////////////////////////////////////////////////////////////////////////////
3800
3801     //calculate InvMass for every V0 particle assumption (Kaon=1,Lambda=2,Antilambda=3)
3802     switch(particletype){
3803     case kK0:
3804       CalculateInvMass(v0, kK0, invM, trackPt); //function to calculate invMass with TLorentzVector class
3805       invMK0s=invM;
3806       break;
3807     case kLambda:
3808       CalculateInvMass(v0, kLambda, invM, trackPt);
3809       invMLa=invM;
3810       break;
3811     case kAntiLambda:
3812       CalculateInvMass(v0, kAntiLambda, invM, trackPt);
3813       invMALa=invM;
3814       break;
3815     default:
3816       std::cout<<"***NO VALID PARTICLETYPE***"<<std::endl;
3817       return 0;
3818     }
3819
3820
3821     /////////////////////////////////////////////////////////////
3822     //V0 and track Cuts:
3823
3824
3825     if(fDebug>6){if(!(IsK0InvMass(invMK0s)) && !(IsLaInvMass(invMLa)) && !(IsLaInvMass(invMALa))){std::cout<<"AliAnalysisTaskJetChem::GetListOfV0s: invM not in selected mass window "<<std::endl;}}
3826
3827     if(!(IsK0InvMass(invMK0s)) && !(IsLaInvMass(invMLa)) && !(IsLaInvMass(invMALa)))continue;
3828
3829     //  Double_t PosEta = trackPos->AliAODTrack::Eta();//daughter track charge is sometimes wrong here, account for that!!!
3830     // Double_t NegEta = trackNeg->AliAODTrack::Eta();
3831
3832    Double_t PosEta = trackPos->Eta();//daughter track charge is sometimes wrong here, account for that!!!
3833    Double_t NegEta = trackNeg->Eta();
3834
3835     Double_t PosCharge = trackPos->Charge();
3836     Double_t NegCharge = trackNeg->Charge();
3837
3838     if((trackPos->Charge() == 1) && (trackNeg->Charge() == -1)) //Fill daughters charge into histo to check if they are symmetric distributed
3839       { fh1PosDaughterCharge->Fill(PosCharge);
3840         fh1NegDaughterCharge->Fill(NegCharge);
3841       }
3842
3843     //DistOverTotMom_in_2D___________
3844
3845     Float_t fMassK0s = TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
3846     Float_t fMassLambda = TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
3847
3848
3849     AliAODVertex* primVtx = fAOD->GetPrimaryVertex(); // get the primary vertex
3850     Double_t dPrimVtxPos[3]; // primary vertex position {x,y,z}
3851     primVtx->GetXYZ(dPrimVtxPos);
3852
3853     Float_t fPtV0 = TMath::Sqrt(v0->Pt2V0()); // transverse momentum of V0
3854     Double_t dSecVtxPos[3]; // V0 vertex position {x,y,z}
3855     v0->GetSecondaryVtx(dSecVtxPos);
3856     Double_t dDecayPath[3];
3857     for (Int_t iPos = 0; iPos < 3; iPos++)
3858       dDecayPath[iPos] = dSecVtxPos[iPos]-dPrimVtxPos[iPos]; // vector of the V0 path
3859     Float_t fDecLen2D = TMath::Sqrt(dDecayPath[0]*dDecayPath[0]+dDecayPath[1]*dDecayPath[1]); //transverse path length R
3860     Float_t fROverPt = fDecLen2D/fPtV0; // R/pT
3861
3862     Float_t fMROverPtK0s = fMassK0s*fROverPt; // m*R/pT
3863     Float_t fMROverPtLambda = fMassLambda*fROverPt; // m*R/pT
3864
3865     //___________________
3866     Double_t fRap = -999;//init values
3867     Double_t fEta = -999;
3868     Double_t fV0cosPointAngle = -999;
3869     Double_t fV0DecayLength = v0->DecayLengthV0(lPrimaryVtxPosition);
3870
3871     Double_t fV0mom[3];
3872
3873     fV0mom[0]=v0->MomV0X();
3874     fV0mom[1]=v0->MomV0Y();
3875     fV0mom[2]=v0->MomV0Z();
3876
3877     Double_t fV0TotalMomentum = TMath::Sqrt(fV0mom[0]*fV0mom[0]+fV0mom[1]*fV0mom[1]+fV0mom[2]*fV0mom[2]);
3878     //  const Double_t K0sPDGmass = 0.497614;
3879     // const Double_t LambdaPDGmass = 1.115683;
3880
3881     const Double_t K0sPDGmass = TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
3882     const Double_t LambdaPDGmass = TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
3883
3884     Double_t fDistOverTotMomK0s = 0;
3885     Double_t fDistOverTotMomLa = 0;
3886
3887     //calculate proper lifetime of particles in 3D (not recommended anymore)
3888
3889     if(particletype == kK0){
3890
3891       fDistOverTotMomK0s = fV0DecayLength * K0sPDGmass;
3892       fDistOverTotMomK0s /= (fV0TotalMomentum+1e-10);
3893     }
3894
3895     if((particletype == kLambda)||(particletype == kAntiLambda)){
3896
3897       fDistOverTotMomLa = fV0DecayLength * LambdaPDGmass;
3898       fDistOverTotMomLa /= (fV0TotalMomentum+1e-10);
3899     }
3900
3901      //TPC cluster (not used anymore) and TPCRefit cuts
3902
3903     //Double_t trackPosNcls = trackPos->GetTPCNcls();//Get number of clusters for positive charged tracks
3904     //Double_t trackNegNcls = trackNeg->GetTPCNcls();//Get number of clusters for negative charged tracks
3905
3906     if(fRequireTPCRefit==kTRUE){//if kTRUE: accept only if daughter track is refitted in TPC!!
3907       Bool_t isPosTPCRefit = (trackPos->AliAODTrack::IsOn(AliESDtrack::kTPCrefit));
3908       Bool_t isNegTPCRefit = (trackNeg->AliAODTrack::IsOn(AliESDtrack::kTPCrefit));
3909       if (!isPosTPCRefit)continue;
3910       if (!isNegTPCRefit)continue;
3911     }
3912
3913     if(fKinkDaughters==kFALSE){//if kFALSE: no acceptance of kink daughters
3914       AliAODVertex* ProdVtxDaughtersPos = (AliAODVertex*) (trackPos->AliAODTrack::GetProdVertex());
3915       Char_t isAcceptKinkDaughtersPos  = ProdVtxDaughtersPos->GetType();
3916       if(isAcceptKinkDaughtersPos==AliAODVertex::kKink)continue;
3917
3918       AliAODVertex* ProdVtxDaughtersNeg = (AliAODVertex*) (trackNeg->AliAODTrack::GetProdVertex());
3919       Char_t isAcceptKinkDaughtersNeg  = ProdVtxDaughtersNeg->GetType();
3920       if(isAcceptKinkDaughtersNeg==AliAODVertex::kKink)continue;
3921
3922     }
3923
3924     Double_t fV0Radius      = -999;
3925     Double_t fDcaV0Daughters = v0->DcaV0Daughters();
3926     Double_t fDcaPosToPrimVertex = v0->DcaPosToPrimVertex();//IP of positive charged daughter
3927     Double_t fDcaNegToPrimVertex = v0->DcaNegToPrimVertex();//IP of negative charged daughter
3928     Double_t avDecayLengthK0s = 2.6844;
3929     Double_t avDecayLengthLa = 7.89;
3930
3931     //Float_t fCTauK0s = 2.6844; // [cm] c tau of K0S
3932     //Float_t fCTauLambda = 7.89; // [cm] c tau of Lambda and Antilambda
3933
3934     fV0cosPointAngle = v0->CosPointingAngle(lPrimaryVtxPosition);
3935     lV0Position[0]= v0->DecayVertexV0X();
3936     lV0Position[1]= v0->DecayVertexV0Y();
3937     lV0Position[2]= v0->DecayVertexV0Z();
3938
3939     fV0Radius  = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
3940
3941     if(particletype == kK0)         {fRap = v0->RapK0Short();
3942                                      fEta = v0->PseudoRapV0();}
3943     if(particletype == kLambda)     {fRap = v0->RapLambda();
3944                                      fEta = v0->PseudoRapV0();}
3945     if(particletype == kAntiLambda) {fRap = v0->Y(-3122);
3946                                      fEta = v0->PseudoRapV0();}
3947
3948
3949     //cut on 3D DistOverTotMom: (not used anymore)
3950     //if((particletype == kLambda)||(particletype == kAntiLambda)){if(fDistOverTotMomLa >= (fCutV0DecayMax * avDecayLengthLa))  continue;}
3951
3952     //cut on K0s applied below after all other cuts for histo fill purposes..
3953
3954     //cut on 2D DistOverTransMom: (recommended from Iouri)
3955     if((particletype == kLambda)||(particletype == kAntiLambda)){if(fMROverPtLambda > (fCutV0DecayMax * avDecayLengthLa))continue;}//fCutV0DecayMax set to 5 in AddTask macro
3956
3957  //Armenteros Podolanski Plot for K0s:////////////////////////////
3958
3959     Double_t ArmenterosAlpha=-999;
3960     Double_t ArmenterosPt=-999;
3961     Double_t PosPl;
3962     Double_t NegPl;
3963     Double_t PosPt;
3964     Double_t NegPt;
3965
3966     if(particletype == kK0){
3967
3968       pp[0]=v0->MomPosX();
3969       pp[1]=v0->MomPosY();
3970       pp[2]=v0->MomPosZ();
3971
3972       pm[0]=v0->MomNegX();
3973       pm[1]=v0->MomNegY();
3974       pm[2]=v0->MomNegZ();
3975
3976
3977       v0mom[0]=v0->MomV0X();
3978       v0mom[1]=v0->MomV0Y();
3979       v0mom[2]=v0->MomV0Z();
3980
3981       TVector3 v0Pos(pp[0],pp[1],pp[2]);
3982       TVector3 v0Neg(pm[0],pm[1],pm[2]);
3983       TVector3 v0totMom(v0mom[0], v0mom[1], v0mom[2]); //vector for tot v0 momentum
3984
3985       PosPt = v0Pos.Perp(v0totMom);             //longitudinal momentum of positive charged daughter track
3986       PosPl = v0Pos.Dot(v0totMom)/v0totMom.Mag();  //transversal momentum of positive charged daughter track
3987
3988       NegPt = v0Neg.Perp(v0totMom);             //longitudinal momentum of negative charged daughter track
3989       NegPl = v0Neg.Dot(v0totMom)/v0totMom.Mag();  //transversal momentum of nergative charged daughter track
3990
3991       ArmenterosAlpha = 1.-2./(1+(PosPl/NegPl));
3992       ArmenterosPt= v0->PtArmV0();
3993
3994     }
3995
3996     if(particletype == kK0){//only cut on K0s histos
3997       if(IsArmenterosSelected == 1){// Armenteros Cut to reject Lambdas contamination in K0s inv. massspectrum
3998         fh2ArmenterosBeforeCuts->Fill(ArmenterosAlpha,ArmenterosPt);
3999       }
4000     }
4001
4002     //some more cuts on v0s and daughter tracks:
4003
4004
4005     if((TMath::Abs(PosEta)>fCutPostrackEta) || (TMath::Abs(NegEta)>fCutNegtrackEta))continue;   //Daughters pseudorapidity cut
4006     if (fV0cosPointAngle < fCutV0cosPointAngle) continue;                                       //cospointangle cut
4007
4008     //if(TMath::Abs(fRap) > fCutRap)continue;                                                     //V0 Rapidity Cut
4009     if(TMath::Abs(fEta) > fCutEta) continue;                                                  //V0 Eta Cut
4010     if (fDcaV0Daughters > fCutDcaV0Daughters)continue;
4011     if ((fDcaPosToPrimVertex < fCutDcaPosToPrimVertex) || (fDcaNegToPrimVertex < fCutDcaNegToPrimVertex))continue;
4012     if ((fV0Radius < fCutV0RadiusMin) || (fV0Radius > fCutV0RadiusMax))continue;
4013
4014     const AliAODPid *pid_p1=trackPos->GetDetPid();
4015     const AliAODPid *pid_n1=trackNeg->GetDetPid();
4016
4017
4018       if(particletype == kLambda){
4019         //      if(AcceptBetheBloch(v0, fPIDResponse, 1) == kFALSE){std::cout<<"******PID cut rejects Lambda!!!************"<<std::endl;}
4020         if(AcceptBetheBloch(v0, fPIDResponse, 1) == kFALSE)continue;
4021         fh2BBLaPos->Fill(pid_p1->GetTPCmomentum(),pid_p1->GetTPCsignal());//positive lambda daughter
4022         fh2BBLaNeg->Fill(pid_n1->GetTPCmomentum(),pid_n1->GetTPCsignal());//negative lambda daughter
4023
4024         //Double_t phi = v0->Phi();
4025         //Double_t massLa = v0->MassLambda();
4026
4027         //printf("La: i = %d, m = %f, pT = %f, eta = %f, phi = %f\n, ",i,massLa,trackPt,fEta,phi);
4028
4029       }
4030
4031       if(particletype == kAntiLambda){
4032
4033         if(AcceptBetheBloch(v0, fPIDResponse, 2) == kFALSE)continue;
4034         fh2BBLaPos->Fill(pid_p1->GetTPCmomentum(),pid_p1->GetTPCsignal());//positive antilambda daughter
4035         fh2BBLaNeg->Fill(pid_n1->GetTPCmomentum(),pid_n1->GetTPCsignal());//negative antilambda daughter
4036
4037       }
4038
4039
4040     //Armenteros cut on K0s:
4041     if(particletype == kK0){
4042       if(IsArmenterosSelected == 1){// Armenteros Cut to reject Lambdas contamination in K0s inv. massspectrum
4043
4044         if((ArmenterosPt<=(TMath::Abs(fCutArmenteros*ArmenterosAlpha))) && (fCutArmenteros!=-999))continue; //Cuts out Lambda contamination in K0s histos
4045         fh2ArmenterosAfterCuts->Fill(ArmenterosAlpha,ArmenterosPt);
4046       }
4047     }
4048
4049     //not used anymore in 3D, z component of total momentum has bad resolution, cut instead in 2D and use pT
4050     //Proper Lifetime Cut: DecayLength3D * PDGmass / |p_tot| < 3*2.68cm (ctau(betagamma=1))  ;  |p|/mass = beta*gamma
4051     //////////////////////////////////////////////
4052
4053     //cut on 3D DistOverTotMom
4054     /*  if(particletype == kK0){
4055
4056       fh2ProperLifetimeK0sVsPtBeforeCut->Fill(trackPt,fDistOverTotMomK0s); //fill these histos after all other cuts
4057       fh1ProperLifetimeV0BeforeCut->Fill(fDistOverTotMomK0s);
4058       if(fDistOverTotMomK0s >= (fCutV0DecayMax * avDecayLengthK0s))continue;
4059       fh1ProperLifetimeV0AfterCut->Fill(fDistOverTotMomK0s);
4060       fh2ProperLifetimeK0sVsPtAfterCut->Fill(trackPt,fDistOverTotMomK0s);
4061     }
4062     */
4063
4064 //cut on 2D DistOverTransMom
4065     if(particletype == kK0){//the cut on Lambdas you can find above
4066
4067       fh2ProperLifetimeK0sVsPtBeforeCut->Fill(trackPt,fMROverPtK0s); //fill these histos after all other cuts
4068       fh1ProperLifetimeV0BeforeCut->Fill(fMROverPtK0s);
4069       if(fMROverPtK0s > (fCutV0DecayMax * avDecayLengthK0s))continue;
4070       fh1ProperLifetimeV0AfterCut->Fill(fMROverPtK0s);
4071       fh2ProperLifetimeK0sVsPtAfterCut->Fill(trackPt,fMROverPtK0s);
4072
4073       //Double_t phi = v0->Phi();
4074       // Double_t massK0s = v0->MassK0Short();
4075       //printf("K0S: i = %d, m = %f, pT = %f, eta = %f, phi = %f\n",i,invMK0s,trackPt,fEta,phi);
4076
4077       //test std::cout<<" Index accepted K0s candidate in list of V0s in event: "<<i<<" m: "<<invMK0s<<" pT: "<<trackPt<<" eta: "<<fEta<<" phi: "<<v0->Phi()<<std::endl;
4078
4079     }
4080     //MC Associated V0 particles: (reconstructed particles associated with MC truth (MC truth: true primary MC generated particle))
4081
4082
4083     if(fAnalysisMC){// begin MC part
4084
4085       Int_t negDaughterpdg = 0;
4086       Int_t posDaughterpdg = 0;
4087       Int_t v0Label = -1;
4088       Bool_t fPhysicalPrimary = -1;   //v0 physical primary check
4089       Int_t MCv0PdgCode = 0;
4090       Bool_t mclabelcheck = kFALSE;
4091
4092       TList *listmc = aod->GetList(); //AliAODEvent* is inherited from AliVEvent*, listmc is pointer to reconstructed event in MC list, member of AliAODEvent
4093
4094       if(!listmc)continue;
4095
4096       if((particletype == kLambda) || (particletype == kAntiLambda)){// at this point the v0 candidates already survived all V0 cuts, for the MC analysis they still have to survive the association checks in the following block
4097
4098         //feeddown-correction for Lambda/Antilambda particles
4099         //feedddown comes mainly from charged and neutral Xi particles
4100         //feeddown from Sigma decays so quickly that it's not possible to distinguish from primary Lambdas with detector
4101         //feeddown for K0s from phi decays is neglectible
4102         //TH2F* fh2FeedDownMatrix = 0x0; //histo for feeddown already decleared above
4103
4104
4105         //first for all Lambda and Antilambda candidates____________________________________________________________________
4106
4107         if(particletype == kLambda){
4108
4109           mclabelcheck = MCLabelCheck(v0, kLambda, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode);
4110
4111
4112           if((motherType == 3312)||(motherType == 3322)){//mother of v0 is neutral or negative Xi
4113
4114             fListFeeddownLaCand->Add(v0); //fill TList with ass. particles, stemming from feeddown from Xi(bar) decays
4115           }
4116         }
4117
4118         if(particletype == kAntiLambda){
4119           mclabelcheck = MCLabelCheck(v0, kAntiLambda, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode);
4120
4121           if((motherType == -3312)||(motherType == -3322)){
4122             fListFeeddownALaCand->Add(v0); //fill TList with ass. particles, stemming from feeddown from Xi(bar) decays
4123           }
4124         }
4125       }
4126
4127       //_only true primary particles survive the following checks_______________________________________________________________________________________________
4128
4129       if(particletype == kK0){
4130         mclabelcheck = MCLabelCheck(v0, kK0, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode);
4131         if(mclabelcheck == kFALSE)continue;
4132       }
4133       if(particletype == kLambda){
4134         mclabelcheck = MCLabelCheck(v0, kLambda, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode);
4135         if(mclabelcheck == kFALSE)continue;
4136       }
4137       if(particletype == kAntiLambda){
4138         mclabelcheck = MCLabelCheck(v0, kAntiLambda, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode);
4139         if(mclabelcheck == kFALSE)continue;
4140       }
4141
4142       if(fPhysicalPrimary != 1)continue; //V0 candidate (K0s, Lambda or Antilambda) must be physical primary, this means there is no mother particle existing
4143
4144     }
4145
4146
4147     list->Add(v0);
4148
4149   }
4150
4151   Int_t nPart=list->GetSize();
4152
4153   return nPart;
4154 } // end GetListOfV0s()
4155
4156 // -------------------------------------------------------------------------------------------------------
4157
4158 void AliAnalysisTaskJetChem::CalculateInvMass(AliAODv0* v0vtx, const Int_t particletype, Double_t& invM, Double_t& trackPt){
4159
4160    //particletype:
4161    // * kaon = 1
4162    // * lambda = 2
4163    // * antilambda = 3
4164
4165      invM    = 0;
4166      trackPt = 0;
4167
4168      Double_t pp[3]={0,0,0}; //3-momentum positive charged track
4169      Double_t pm[3]={0,0,0}; //3-momentum negative charged track
4170
4171      const Double_t massPi = 0.13957018; //better use PDG code at this point
4172      const Double_t massP  = 0.93827203;
4173
4174      Double_t mass1=0;
4175      Double_t mass2=0;
4176
4177      TLorentzVector vector;  //lorentzvector V0 particle
4178      TLorentzVector fourmom1;//lorentzvector positive daughter
4179      TLorentzVector fourmom2;//lorentzvector negative daughter
4180
4181    //--------------------------------------------------------------
4182
4183    AliAODTrack *trackPos = (AliAODTrack *) (v0vtx->GetSecondaryVtx()->GetDaughter(0));//index 0 defined as positive charged track in AliESDFilter
4184
4185    if( trackPos->Charge() == 1 ){
4186
4187      pp[0]=v0vtx->MomPosX();
4188      pp[1]=v0vtx->MomPosY();
4189      pp[2]=v0vtx->MomPosZ();
4190
4191      pm[0]=v0vtx->MomNegX();
4192      pm[1]=v0vtx->MomNegY();
4193      pm[2]=v0vtx->MomNegZ();
4194    }
4195
4196    if( trackPos->Charge() == -1 ){
4197
4198      pm[0]=v0vtx->MomPosX();
4199      pm[1]=v0vtx->MomPosY();
4200      pm[2]=v0vtx->MomPosZ();
4201
4202      pp[0]=v0vtx->MomNegX();
4203      pp[1]=v0vtx->MomNegY();
4204      pp[2]=v0vtx->MomNegZ();
4205    }
4206
4207    if (particletype == kK0){ // case K0s
4208      mass1 = massPi;//positive particle
4209      mass2 = massPi;//negative particle
4210    } else if (particletype == kLambda){ // case Lambda
4211      mass1 = massP;//positive particle
4212      mass2 = massPi;//negative particle
4213    } else if (particletype == kAntiLambda){ //case AntiLambda
4214      mass1 = massPi;//positive particle
4215      mass2 = massP; //negative particle
4216    }
4217
4218    fourmom1.SetXYZM(pp[0],pp[1],pp[2],mass1);//positive track
4219    fourmom2.SetXYZM(pm[0],pm[1],pm[2],mass2);//negative track
4220    vector=fourmom1 + fourmom2;
4221
4222    invM    = vector.M();
4223    trackPt = vector.Pt();
4224
4225    /*// don't apply AliAODv0 methods to get the inv. mass for the OnFly finder, since the daughter labels are sometimes switched!!!! For Offline V0 finder no problem
4226
4227    if(particletype == kK0){
4228      std::cout << "invMK0s: " << invM <<std::endl;
4229      std::cout << "v0vtx->MassK0Short(): " << v0vtx->MassK0Short() << std::endl;
4230      std::cout << "    " <<std::endl;
4231      //invM = v0vtx->MassK0Short();
4232    }
4233
4234    if(particletype == kLambda){
4235      std::cout << "invMLambda: " << invM <<std::endl;
4236      std::cout << "v0vtx->MassMassLambda(): " << v0vtx->MassLambda() << std::endl;
4237      std::cout << "    " <<std::endl;
4238     //invM = v0vtx->MassLambda();
4239    }
4240
4241    if(particletype == kAntiLambda){
4242      std::cout << "invMAntiLambda: " << invM <<std::endl;
4243      std::cout << "v0vtx->MassAntiLambda(): " << v0vtx->MassAntiLambda() << std::endl;
4244      std::cout << "    " <<std::endl;
4245      //invM = v0vtx->MassAntiLambda();
4246    }
4247    */
4248
4249    return;
4250 }
4251
4252
4253 //_____________________________________________________________________________________
4254 Int_t AliAnalysisTaskJetChem::GetListOfMCParticles(TList *outputlist, const Int_t particletype, AliAODEvent *mcaodevent) //(list to fill here e.g. fListMCgenK0s, particle species to search for)
4255 {
4256
4257   outputlist->Clear();
4258
4259   TClonesArray *stack = 0x0;
4260   Double_t mcXv=0., mcYv=0., mcZv=0.;//MC vertex position
4261   Int_t ntrk =0;
4262
4263   // get MC generated particles
4264
4265   Int_t fPdgcodeCurrentPart = 0; //pdg code current particle
4266   Double_t fRapCurrentPart  = 0; //get rapidity
4267   Double_t fPtCurrentPart   = 0; //get transverse momentum
4268   Double_t fEtaCurrentPart = 0;  //get pseudorapidity
4269
4270   //variable for check: physical primary particle
4271   //Bool_t IsPhysicalPrimary = -1;
4272   //Int_t index = 0; //check number of injected particles
4273   //****************************
4274   // Start loop over MC particles
4275
4276   TList *lst = mcaodevent->GetList();
4277
4278   if(!lst){
4279     if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
4280     return -1;
4281   }
4282
4283   stack = (TClonesArray*)lst->FindObject(AliAODMCParticle::StdBranchName());
4284   if (!stack) {
4285     Printf("ERROR: stack not available");
4286     return -1;
4287   }
4288
4289   AliAODMCHeader *mcHdr=(AliAODMCHeader*)lst->FindObject(AliAODMCHeader::StdBranchName());
4290   if(!mcHdr)return -1;
4291
4292   mcXv=mcHdr->GetVtxX(); mcYv=mcHdr->GetVtxY(); mcZv=mcHdr->GetVtxZ(); // position of the MC primary vertex
4293
4294
4295   ntrk=stack->GetEntriesFast();
4296
4297   //if(TMath::Abs(mcZv)>10)return; //i also cut at the reconstructed particles - here i also want to cut for a second time on z vertex (?) -> could be possible bias because of resolution effects on edges of acceptance, also the case for pseudorapidity...
4298
4299
4300   for (Int_t iMc = 0; iMc < ntrk; iMc++) {                                             //loop over mc generated particles
4301
4302
4303     AliAODMCParticle *p0=(AliAODMCParticle*)stack->UncheckedAt(iMc);
4304     if (!p0) {
4305       //Printf("ERROR: particle with label %d not found in stack (mc loop)", iMc);
4306       continue;
4307     }
4308     fPdgcodeCurrentPart = p0->GetPdgCode();
4309
4310     // Keep only K0s, Lambda and AntiLambda, Xi and Phi:
4311     //if ( (fPdgcodeCurrentPart != 310 ) && (fPdgcodeCurrentPart != 3122 ) && (fPdgcodeCurrentPart != -3122 ) && (fPdgcodeCurrentPart != 3312 ) && (fPdgcodeCurrentPart != -3312) && (fPdgcodeCurrentPart != -333) ) continue;
4312
4313
4314
4315       //Rejection of Pythia injected particles with David Chinellatos method - not the latest method, better Method with TString from MC generator in IsInjected() function below!
4316
4317     /*     if( (p0->GetStatus()==21) ||
4318           ((p0->GetPdgCode() == 443) &&
4319            (p0->GetMother() == -1) &&
4320            (p0->GetDaughter(0) == (iMc))) ){ index++; }
4321
4322       if(p0->GetStatus()==21){std::cout<< "hello !!!!" <<std::endl;}
4323
4324       std::cout<< "MC particle status:  " << p0->GetStatus() <<std::endl;
4325     */
4326
4327
4328       //if(index>=1){std::cout<< "MC particle status:  " << p0->GetStatus() <<std::endl;}//if first injected MC particle was found, the Status is printed out for this and every following MC particle
4329
4330
4331         //injected particles could be from GenBox (single high pt tracks) or jet related tracks, both generated from PYTHIA MC generator
4332
4333       //Check: MC particle mother
4334
4335       //for feed-down checks
4336       /* //MC gen particles
4337         Int_t iMother = p0->GetMother(); //Motherparticle of V0 candidate (e.g. phi particle,..)
4338            if(iMother >= 0){
4339         AliAODMCParticle *partM = (AliAODMCParticle*)stack->UncheckedAt(iMother);
4340         Int_t codeM = -1;
4341         if(partM) codeM = TMath::Abs(partM->GetPdgCode());
4342
4343
4344          3312    Xi-                       -3312    Xibar+
4345          3322    Xi0                       -3322    Xibar0
4346
4347
4348         if((codeM == 3312)||(codeM == 3322))// feeddown for Lambda coming from Xi- and Xi0
4349
4350
4351
4352            }
4353         */
4354            /*   //Check: MC gen. particle decays via 2-pion decay? -> only to be done for the rec. particles !! (-> branching ratio ~ 70 % for K0s -> pi+ pi-)
4355
4356         Int_t daughter0Label = p0->GetDaughter(0);
4357         AliAODMCParticle *mcDaughter0 = (AliAODMCParticle *)stack->UncheckedAt(daughter0Label);
4358       if(daughter0Label >= 0)
4359         {daughter0Type = mcDaughter0->GetPdgCode();}
4360
4361       Int_t daughter1Label = p0->GetDaughter(1);
4362       AliAODMCParticle *mcDaughter1 = (AliAODMCParticle *)stack->UncheckedAt(daughter1Label);
4363
4364       if(daughter1Label >= 1)
4365         {daughter1Type = mcDaughter1->GetPdgCode();}    //requirement that daughters are pions is only done for the reconstructed V0s in GetListofV0s() below
4366     }
4367            */
4368
4369
4370       // Keep only K0s, Lambda and AntiLambda:
4371       if ( (fPdgcodeCurrentPart != 310 ) && (fPdgcodeCurrentPart != 3122 ) && (fPdgcodeCurrentPart != -3122 )) continue;
4372       // Check: Is physical primary
4373
4374       //do not use anymore: //IsPhysicalPrimary = p0->IsPhysicalPrimary();
4375       //if(!IsPhysicalPrimary)continue;
4376
4377       Float_t fDistPrimaryMax = 0.01; // [cm] max distance of production point to the primary vertex (criterion for choice of MC particles considered as primary)
4378
4379       // Get the distance between production point of the MC mother particle and the primary vertex
4380
4381       Double_t dx = mcXv-p0->Xv();//mc primary vertex - mc gen. v0 vertex
4382       Double_t dy = mcYv-p0->Yv();
4383       Double_t dz = mcZv-p0->Zv();
4384
4385       Double_t fDistPrimary = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
4386       Bool_t fPhysicalPrimary = (fDistPrimary < fDistPrimaryMax);
4387
4388       if(!fPhysicalPrimary)continue;
4389
4390       //if(fPhysicalPrimary){std::cout<<"hello**********************"<<std::endl;}
4391
4392       /* std::cout<<"dx: "<<dx<<std::endl;
4393       std::cout<<"dy: "<<dy<<std::endl;
4394       std::cout<<"dz: "<<dz<<std::endl;
4395
4396       std::cout<<"start: "<<std::endl;
4397       std::cout<<"mcXv: "<<mcXv<<std::endl;
4398       std::cout<<"mcYv: "<<mcYv<<std::endl;
4399       std::cout<<"mcZv: "<<mcZv<<std::endl;
4400
4401       std::cout<<"p0->Xv(): "<<p0->Xv()<<std::endl;
4402       std::cout<<"p0->Yv(): "<<p0->Yv()<<std::endl;
4403       std::cout<<"p0->Zv(): "<<p0->Zv()<<std::endl;
4404       std::cout<<" "<<std::endl;
4405       std::cout<<"fDistPrimary"<<fDistPrimary<<std::endl;
4406       std::cout<<"fPhysicalPrimary"<<fPhysicalPrimary<<std::endl;
4407       */
4408
4409       //Is close enough to primary vertex to be considered as primary-like?
4410
4411       fRapCurrentPart   = MyRapidity(p0->E(),p0->Pz());
4412       fEtaCurrentPart   = p0->Eta();
4413       fPtCurrentPart    = p0->Pt();
4414
4415       if (TMath::Abs(fEtaCurrentPart) < fCutEta){
4416         // if (TMath::Abs(fRapCurrentPart) > fCutRap)continue;    //rap cut for crosschecks
4417
4418         if(particletype == kK0){                                      //MC gen. K0s
4419           if (fPdgcodeCurrentPart==310){
4420             outputlist->Add(p0);
4421           }
4422         }
4423
4424     if(particletype == kLambda){                                //MC gen. Lambdas
4425       if (fPdgcodeCurrentPart==3122)  {
4426         outputlist->Add(p0);
4427       }
4428     }
4429
4430     if(particletype == kAntiLambda){
4431       if (fPdgcodeCurrentPart==-3122) {                         //MC gen. Antilambdas
4432         outputlist->Add(p0);
4433       }
4434     }
4435   }
4436
4437   }//end  loop over MC generated particle
4438
4439   Int_t nMCPart=outputlist->GetSize();
4440
4441
4442   return nMCPart;
4443
4444 }
4445
4446 //---------------------------------------------------------------------------
4447 /*
4448 Bool_t AliAnalysisTaskJetChem::FillFeeddownMatrix(TList* fListFeeddownCand, Int_t particletype)
4449 {
4450
4451         // Define Feeddown matrix
4452         Double_t lFeedDownMatrix [100][100];
4453         // FeedDownMatrix [Lambda Bin][Xi Bin];
4454
4455         //Initialize entries of matrix:
4456         for(Int_t ilb = 0; ilb<100; ilb++){
4457           for(Int_t ixb = 0; ixb<100; ixb++){
4458             lFeedDownMatrix[ilb][ixb]=0; //first lambda bins, xi bins
4459           }
4460         }
4461 }
4462 */
4463 //----------------------------------------------------------------------------
4464
4465 Double_t AliAnalysisTaskJetChem::MyRapidity(Double_t rE, Double_t rPz) const
4466 {
4467   // Local calculation for rapidity
4468   return 0.5*TMath::Log((rE+rPz)/(rE-rPz+1.e-13));
4469 }
4470 //----------------------------------------------------------------------------
4471
4472 // ________________________________________________________________________________________________________________________//function to get the MC gen. jet particles
4473
4474 void AliAnalysisTaskJetChem::GetTracksInCone(TList* inputlist, TList* outputlist, const AliAODJet* jet,
4475                                                                 const Double_t radius, Double_t& sumPt, const Double_t minPt, const Double_t maxPt, Bool_t& isBadPt)
4476 {
4477   // fill list of tracks in cone around jet axis
4478
4479   sumPt = 0;
4480   Bool_t isBadMaxPt = kFALSE;
4481   Bool_t isBadMinPt = kTRUE;
4482
4483   Double_t jetMom[3];
4484   if(!jet)return;
4485   jet->PxPyPz(jetMom);
4486   TVector3 jet3mom(jetMom);
4487
4488   for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){
4489
4490     AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack));
4491     if(!track)continue;
4492     Double_t trackMom[3];
4493     track->PxPyPz(trackMom);
4494     TVector3 track3mom(trackMom);
4495
4496     Double_t dR = jet3mom.DeltaR(track3mom);
4497
4498     if(dR<radius){
4499
4500       outputlist->Add(track);
4501
4502       sumPt += track->Pt();
4503
4504       if(maxPt>0 && track->Pt()>maxPt) isBadMaxPt = kTRUE;   // reject jets containing any track with pt larger than this value, use GetFFMaxTrackPt()
4505       if(minPt>0 && track->Pt()>minPt) isBadMinPt = kFALSE;  // reject jets with leading track with pt smaller than this value, use GetFFMinLTrackPt()
4506
4507     }
4508   }
4509
4510   isBadPt = kFALSE;
4511   if(minPt>0 && isBadMinPt) isBadPt = kTRUE; //either the jet is bad because of too small leading track pt.. (probability to be purely combinatorial jet is too high to accept it)
4512   if(maxPt>0 && isBadMaxPt) isBadPt = kTRUE; //..or because of leading track with too high pt (could be fake track)
4513
4514   outputlist->Sort();
4515
4516 }
4517
4518 //____________________________________________________________________________________________________________________
4519
4520
4521 void AliAnalysisTaskJetChem::GetTracksInPerpCone(TList* inputlist, TList* outputlist, const AliAODJet* jet,
4522                                                                 const Double_t radius, Double_t& sumPerpPt)
4523 {
4524   // fill list of tracks in two cones around jet axis rotated in phi +/- 90 degrees
4525
4526   Double_t jetMom[3];             //array for entries in TVector3
4527   Double_t perpjetplusMom[3];     //array for entries in TVector3
4528   Double_t perpjetnegMom[3];
4529
4530   if(!jet)return;
4531
4532   jet->PxPyPz(jetMom);   //get 3D jet momentum
4533   Double_t jetPerpPt = jet->Pt(); //original jet pt, invariant under rotations
4534   Double_t jetPhi = jet->Phi(); //original jet phi
4535
4536   Double_t jetPerpposPhi = jetPhi + ((TMath::Pi())*0.5);//get new perp. jet axis phi clockwise
4537   Double_t jetPerpnegPhi = jetPhi - ((TMath::Pi())*0.5);//get new perp. jet axis phi counterclockwise
4538
4539   TVector3 jet3mom(jetMom); //3-Vector for original rec. jet axis
4540
4541   //Double_t phitest  = jet3mom.Phi();
4542
4543   perpjetplusMom[0]=(TMath::Cos(jetPerpposPhi)*jetPerpPt); //x coordinate (sidewards - when looking in beam direction)
4544   perpjetplusMom[1]=(TMath::Sin(jetPerpposPhi)*jetPerpPt); //y coordinate (upwards - when looking in beam direction)
4545   perpjetplusMom[2]=jetMom[2];                          //z coordinate (along beam axis), invariant under azimuthal rotation
4546
4547   perpjetnegMom[0]=(TMath::Cos(jetPerpnegPhi)*jetPerpPt); //x coordinate (sidewards - when looking in beam direction)
4548   perpjetnegMom[1]=(TMath::Sin(jetPerpnegPhi)*jetPerpPt); //y coordinate (upwards - when looking in beam direction)
4549   perpjetnegMom[2]=jetMom[2];                          //z coordinate (along beam axis), invariant under azimuthal rotation
4550
4551
4552   TVector3 perpjetplus3mom(perpjetplusMom);  //3-Vector for new perp. jet axis, clockwise rotated
4553   TVector3 perpjetneg3mom(perpjetnegMom);    //3-Vector for new perp. jet axis, counterclockwise rotated
4554
4555   //for crosscheck TVector3 rotation method
4556
4557   //Double_t jetMomplusTest[3];
4558   //Double_t jetMomminusTest[3];
4559
4560   //jet3mom.RotateZ(TMath::Pi()*0.5);//rotate original jet axis around +90 degrees in phi
4561
4562   //perpjetminus3momTest  = jet3mom.RotateZ((-1)*TMath::Pi()*0.5);
4563
4564   // jet3mom.RotateZ(TMath::Pi()*0.5);
4565   // jet3mom.RotateZ((-1)*TMath::Pi()*0.5);
4566
4567   //jetMomplusTest[0] = jet3mom.X();  //fetching perp. axis coordinates
4568   //jetMomplusTest[1] = jet3mom.Y();
4569   //jetMomplusTest[2] = jet3mom.Z();
4570
4571   //TVector3 perpjetplus3momTest(jetMomplusTest);   //new TVector3 for +90deg rotated jet axis with rotation method from ROOT
4572   //TVector3 perpjetminus3momTest(jetMomminusTest);   //new TVector3 for -90deg rotated jet axis with rotation method from ROOT
4573
4574
4575   for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){ //collect V0 content in perp cone, rotated clockwise
4576
4577     AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack)); //inputlist is fListK0s, all reconstructed K0s in event
4578     if(!track){std::cout<<"K0s track not found!!!"<<std::endl; continue;}
4579
4580     Double_t trackMom[3];//3-mom of V0 particle
4581     track->PxPyPz(trackMom);
4582     TVector3 track3mom(trackMom);
4583
4584     Double_t dR = perpjetplus3mom.DeltaR(track3mom);
4585
4586     if(dR<radius){
4587
4588       outputlist->Add(track); // output list is jetPerpConeK0list
4589
4590       sumPerpPt += track->Pt();
4591
4592
4593     }
4594   }
4595
4596
4597   for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){//collect V0 content in perp cone, rotated counterclockwise
4598
4599     AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack)); //inputlist is fListK0s, all reconstructed K0s in event
4600     if(!track){std::cout<<"K0s track not found!!!"<<std::endl; continue;}
4601
4602     Double_t trackMom[3];//3-mom of V0 particle
4603     track->PxPyPz(trackMom);
4604     TVector3 track3mom(trackMom);
4605
4606     Double_t dR = perpjetneg3mom.DeltaR(track3mom);
4607
4608     if(dR<radius){
4609
4610       outputlist->Add(track); // output list is jetPerpConeK0list
4611
4612       sumPerpPt += track->Pt();
4613
4614
4615     }
4616   }
4617
4618   // pay attention: this list contains the double amount of V0s, found in both cones
4619   // before using it, devide spectra by 2!!!
4620   sumPerpPt = sumPerpPt*0.5; //correct to do this?
4621
4622
4623   outputlist->Sort();
4624
4625 }
4626
4627
4628 // _______________________________________________________________________________________________________________________________________________________
4629
4630 Bool_t AliAnalysisTaskJetChem::MCLabelCheck(AliAODv0* v0, Int_t particletype,const AliAODTrack* trackNeg, const AliAODTrack* trackPos, TList *listmc, Int_t& negDaughterpdg, Int_t& posDaughterpdg, Int_t& motherType, Int_t& v0Label, Double_t& MCPt, Bool_t& fPhysicalPrimary, Int_t& MCv0PDGCode){
4631
4632   TClonesArray *stackmc = 0x0;
4633   stackmc = (TClonesArray*)listmc->FindObject(AliAODMCParticle::StdBranchName()); //get MCAOD branch in data
4634   if (!stackmc)
4635     {
4636       Printf("ERROR: stack not available");
4637       return kFALSE;
4638     }
4639   else
4640     {
4641       Int_t negAssLabel = TMath::Abs(trackNeg->GetLabel());                       //negative (reconstructed) charged track label in MC stack
4642       Int_t posAssLabel = TMath::Abs(trackPos->GetLabel());                       //positive (reconstructed) charged track label in MC stack
4643
4644       //injected particle checks
4645
4646       Double_t mcXv = 0;
4647       Double_t mcYv = 0;
4648       Double_t mcZv = 0;
4649
4650       AliAODMCHeader *header=(AliAODMCHeader*)listmc->FindObject(AliAODMCHeader::StdBranchName());
4651       if(!header)return kFALSE;
4652
4653       mcXv=header->GetVtxX(); mcYv=header->GetVtxY(); mcZv=header->GetVtxZ();
4654
4655       Int_t trackinjected = IsTrackInjected(v0, header, stackmc); //requires AliAODTrack instead of AliVTrack
4656
4657       if(trackinjected == 0){std::cout<<"HIJING track injected!!: "<<trackinjected<<std::endl;}
4658
4659       //mc label checks
4660
4661       if(negAssLabel>=0 && negAssLabel < stackmc->GetEntriesFast() && posAssLabel>=0 && posAssLabel < stackmc->GetEntriesFast()){//safety check if label has valid value of stack
4662
4663         AliAODMCParticle *mcNegPart =(AliAODMCParticle*)stackmc->UncheckedAt(negAssLabel);//fetch the, with one MC truth track associated (reconstructed), negative charged track
4664         v0Label = mcNegPart->GetMother();// mother of negative charged particle is v0, get v0 label here
4665         negDaughterpdg = mcNegPart->GetPdgCode();
4666         AliAODMCParticle *mcPosPart =(AliAODMCParticle*)stackmc->UncheckedAt(posAssLabel);//fetch the, with one MC truth track associated (reconstructed), positive charged track
4667         Int_t v0PosLabel = mcPosPart->GetMother();                                        //get mother label of positive charged track label
4668         posDaughterpdg = mcPosPart->GetPdgCode();
4669
4670         if(v0Label >= 0 && v0Label < stackmc->GetEntriesFast() && v0Label == v0PosLabel){//first v0 mc label check, then: check if both daughters are stemming from same particle
4671
4672           AliAODMCParticle *mcv0 = (AliAODMCParticle *)stackmc->UncheckedAt(v0Label);  //fetch MC ass. particle to v0 (mother of the both charged daughter tracks)
4673
4674           //do not use anymore:
4675           //fPhysicalPrimary = mcv0->IsPhysicalPrimary();
4676
4677
4678           Float_t fDistPrimaryMax = 0.01; // [cm] max distance of production point to the primary vertex (criterion for choice of MC particles considered as primary)
4679
4680           // Get the distance between production point of the MC mother particle and the primary vertex
4681
4682           Double_t dx = mcXv-mcv0->Xv();//mc primary vertex - mc particle production vertex
4683           Double_t dy = mcYv-mcv0->Yv();
4684           Double_t dz = mcZv-mcv0->Zv();
4685
4686           Float_t fDistPrimary = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
4687           fPhysicalPrimary = kFALSE;//init
4688
4689           fPhysicalPrimary = (fDistPrimary < fDistPrimaryMax);
4690
4691           //if(fPhysicalPrimary == kTRUE){std::cout<<"hello*********!!!!!!!!!!!!! "<<std::endl;}
4692           //if(fPhysicalPrimary == kFALSE)return kFALSE;
4693
4694           MCv0PDGCode = mcv0->GetPdgCode();
4695
4696           //std::cout<<"MCv0PDGCode: "<<MCv0PDGCode<<std::endl;
4697
4698           MCPt = mcv0->Pt();//for MC data, always use MC gen. pt for any pt distributions, also for the spectra, used for normalisation
4699           //for feed-down checks later
4700
4701           Int_t motherLabel = mcv0->GetMother();  //get mother particle label of v0 particle
4702           // std::cout<<"motherLabel: "<<motherLabel<<std::endl;
4703
4704           if(motherLabel >= 0 && v0Label < stackmc->GetEntriesFast())                 //do safety check for mother label
4705             {
4706               AliAODMCParticle *mcMother = (AliAODMCParticle *)stackmc->UncheckedAt(motherLabel);  //get mother particle
4707               motherType = mcMother->GetPdgCode(); //get PDG code of mother
4708
4709               Double_t XiPt = 0.;
4710               Double_t XibarPt = 0.;
4711
4712               if(particletype == kLambda){
4713                 if((motherType == 3312)||(motherType == 3322)){ //if v0 mother is Xi0 or Xi- fill MC gen. pt in FD La histogram
4714                   XiPt = mcMother->Pt();
4715                   fh1MCXiPt->Fill(XiPt);
4716                 }
4717               }
4718               if(particletype == kAntiLambda){
4719                 if((motherType == -3312)||(motherType == -3322)){ //if v0 mother is Xibar0 or Xibar+ fill MC gen. pt in FD ALa histogram
4720                   XibarPt = mcMother->Pt();
4721                   fh1MCXibarPt->Fill(XibarPt);
4722                 }
4723                 }
4724
4725             }
4726
4727           //pdg code checks etc..
4728
4729           if(particletype == kK0){
4730
4731                if(TMath::Abs(posDaughterpdg) != 211){return kFALSE;}//one or both of the daughters are not a pion
4732                if(TMath::Abs(negDaughterpdg) != 211){return kFALSE;}
4733
4734                if(MCv0PDGCode != 310)  {return kFALSE;}
4735              }
4736
4737              if(particletype == kLambda){
4738                if(MCv0PDGCode != 3122)return kFALSE;//if particle is not Antilambda, v0 is rejected
4739                if(posDaughterpdg != 2212)return kFALSE;
4740                if(negDaughterpdg != -211)return kFALSE;    //pdg code check for Lambda daughters
4741
4742                //{if((motherType == 3312)||(motherType == 3322)){continue;}//if Xi0 and Xi- is motherparticle of Lambda, particle is rejected, pay attention, most possible Xi-, Xi0 and Omega- are not distributed physically and are much more abundant than expected by physics       //}
4743              }
4744
4745              if(particletype == kAntiLambda){
4746                if(MCv0PDGCode != -3122)return kFALSE;
4747                if(posDaughterpdg != 211)return kFALSE;
4748                if(negDaughterpdg !=-2212)return kFALSE;    //pdg code check for Antilambda daughters
4749                //if(fPhysicalPrimary == 1){fh1MCmotherALa->Fill(7.);}
4750
4751                //{if((motherType == -3312)||(motherType == -3322)){continue;}//if bar{Xi0} and Xi+ is motherparticle of Antilambda, particle is rejected
4752                //}
4753              }
4754
4755              return kTRUE;                     //check was successful
4756            }//end mc v0 label check
4757       }// end of stack label check
4758     }//end of else
4759
4760
4761
4762   return kFALSE;                               //check wasn't successful
4763 }
4764 //________________________________________________________________________________________________________________________________________________________
4765
4766
4767 Bool_t AliAnalysisTaskJetChem::IsParticleMatching(const AliAODMCParticle* mcp0, const Int_t v0Label){
4768
4769   const Int_t mcp0label = mcp0->GetLabel();
4770
4771   if(v0Label == mcp0label)return kTRUE;
4772
4773   return kFALSE;
4774 }
4775
4776 //_______________________________________________________________________________________________________________________________________________________
4777
4778 Bool_t AliAnalysisTaskJetChem::DaughterTrackCheck(AliAODv0* v0, Int_t& nnum, Int_t& pnum){
4779
4780
4781   if(v0->GetNDaughters() != 2) return kFALSE;//case v0 has more or less than 2 daughters, avoids seg. break at some AOD files                                   //reason?
4782
4783
4784   // safety check of input parameters
4785   if(v0 == NULL)
4786     {
4787       if(fDebug > 1){std::cout << std::endl
4788                 << "Warning in AliAnalysisTaskJetChem::DaughterTrackCheck:" << std::endl
4789                            << "v0 = " << v0 << std::endl;}
4790
4791       return kFALSE;
4792     }
4793   else
4794     {
4795       //Daughters track check: its Luke Hanrattys method to check daughters charge
4796
4797       nnum = 1;
4798       pnum = 0;
4799
4800
4801       AliAODTrack *ntracktest =(AliAODTrack*)(v0->GetDaughter(nnum));
4802
4803       if(ntracktest == NULL)
4804         {
4805           if(fDebug > 1){std::cout << std::endl
4806                     << "Warning in AliAnalysisTaskJetChem::DaughterTrackCheck:" << std::endl
4807                          << "ntracktest = " << ntracktest << std::endl;}
4808
4809           return kFALSE;
4810         }
4811
4812       if(ntracktest->Charge() > 0)
4813         {
4814           nnum = 0;
4815           pnum = 1;
4816         }
4817
4818       const AliAODTrack *trackNeg=(AliAODTrack *)(v0->GetDaughter(nnum));
4819       const AliAODTrack *trackPos=(AliAODTrack *)(v0->GetDaughter(pnum));
4820
4821       //Check if both tracks are available
4822       if (!trackPos || !trackNeg) {
4823         if(fDebug > 1) Printf("strange analysis::UserExec:: Error:Could not retrieve one of the daughter tracks\n");
4824         return kFALSE;
4825       }
4826
4827
4828       //remove like sign V0s
4829       if ( trackPos->Charge() == trackNeg->Charge() ){
4830         //if(fDebug>1) Printf("%s:%d found like-sign V0", (char*)__FILE__,__LINE__);
4831         return kFALSE;
4832       }
4833
4834       return kTRUE;
4835     }
4836 }
4837
4838 //_______________________________________________________________________________________________________________________________________________________
4839
4840 Int_t AliAnalysisTaskJetChem::IsTrackInjected(AliAODv0 *v0, AliAODMCHeader *header, TClonesArray *arrayMC){//info in TString should be available from 2011 data productions on..
4841
4842   if(!v0){std::cout << " !part " << std::endl;return 1;}
4843   if(!header){std::cout << " !header " << std::endl;return 1;}
4844   if(!arrayMC){std::cout << " !arrayMC " << std::endl;return 1;}
4845
4846   Int_t lab=v0->GetLabel();
4847   if(lab<0) {return 1;}
4848   TString bbb = GetGenerator(lab,header);
4849   TString empty="";
4850
4851   // std::cout << " TString bbb: " << bbb << std::endl;
4852
4853   //  std::cout << " FIRST CALL " << bbb << std::endl;
4854
4855   while(bbb.IsWhitespace()){
4856     AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab);
4857     if(!mcpart){return 1;}
4858     Int_t mother = mcpart->GetMother();
4859     lab=mother;
4860     bbb = GetGenerator(mother,header);
4861     std::cout << "Testing " << bbb << " "  << std::endl;
4862   }
4863
4864   std::cout << " FINAL CALL " << bbb << std::endl;
4865
4866   //std::transform(bbb.begin(), bbb.end(), bbb.begin(), ::tolower);   //convert TString bbb into lower case, to avoid that TString could be  written in lower or upper case
4867
4868   if(bbb.Contains("ijing")){std::cout << " particle is injected!! " << std::endl; return 0;}//if TString returns something with "ijing" return this method with 0 -> select out all HIJING particles, all others return with "1"
4869
4870
4871   return 1;
4872 }
4873
4874 //______________________________________________________________________
4875 TString AliAnalysisTaskJetChem::GetGenerator(Int_t label, AliAODMCHeader* header){
4876    Int_t nsumpart=0;
4877    TList *lh=header->GetCocktailHeaders();
4878    Int_t nh=lh->GetEntries();
4879    for(Int_t i=0;i<nh;i++){
4880      AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i);
4881      TString genname=gh->GetName();
4882      Int_t npart=gh->NProduced();
4883      if(label>=nsumpart && label<(nsumpart+npart)) return genname;
4884      nsumpart+=npart;
4885    }
4886    TString empty="";
4887    return empty;
4888  }
4889
4890 //_________________________________________________________________________________________________________________________________________
4891 Double_t AliAnalysisTaskJetChem::SmearJetPt(Double_t jetPt, Int_t /*cent*/, Double_t /*jetRadius*/, Double_t /*ptmintrack*/, Double_t& jetPtSmear){
4892
4893   static TF1 fsmear("f1","[0]*exp(-1*(x-[1])*(x-[1])/(2*[2]*[2]))",-100.,100.);   //smearing according to gaussian function in between  +/- 10 GeV/c
4894
4895   //Int_t cl = 1;
4896
4897   /*  if(cent>10) cl = 2;
4898   if(cent>30) cl = 3;
4899   if(cent>50) cl = 4;
4900   */
4901
4902   fsmear.SetParameters(1,0,11.19);//for 2010 PbPb jets, R=0.4, ptmintrack = 0.15 GeV/c, cent 00-10%, delta-pt width estimated via single track embedding
4903   //fsmear->SetParameters(1,0,3.28);//for 2010 PbPb jets, R=0.4, ptmintrack = 0.15 GeV/c, cent 50-60%, delta-pt width estimated via single track embedding
4904
4905   //fsmear->SetParameters(1,0,4.472208);// for 2010 PbPb jets, R=0.2, ptmintrack = 0.15 GeV/c, cent 00-10%
4906
4907   /* //delta-pt width for anti-kt jet finder:
4908
4909   // jet cone R = 0.4
4910   if((cl == 1)&&(jetRadius == 0.4)&&(ptmintrack == 0.15)){
4911   fsmear->SetParameters(1,0,10.178069);//(max.,mean,sigma) of gaussian, needs to be adjusted for every combination of jet cone size, centrality and min. pt constituents cut
4912   }
4913   if((cl == 2)&&(jetRadius == 0.4)&&(ptmintrack == 0.15)){
4914   fsmear->SetParameters(1,0,8.536195);
4915   }
4916   if((cl == 3)&&(jetRadius == 0.4)&&(ptmintrack == 0.15)){
4917   fsmear->SetParameters(1,0,?);
4918   }
4919   if((cl == 4)&&(jetRadius == 0.4)&&(ptmintrack == 0.15)){
4920   fsmear->SetParameters(1,0,5.229839);
4921   }
4922
4923   // jet cone R = 0.3
4924   if((cl == 1)&&(jetRadius == 0.3)&&(ptmintrack == 0.15)){
4925   fsmear->SetParameters(1,0,7.145967);
4926   }
4927   if((cl == 2)&&(jetRadius == 0.3)&&(ptmintrack == 0.15)){
4928   fsmear->SetParameters(1,0,5.844796);
4929   }
4930   if((cl == 3)&&(jetRadius == 0.3)&&(ptmintrack == 0.15)){
4931   fsmear->SetParameters(1,0,?);
4932   }
4933   if((cl == 4)&&(jetRadius == 0.3)&&(ptmintrack == 0.15)){
4934   fsmear->SetParameters(1,0,3.630751);
4935   }
4936
4937   // jet cone R = 0.2
4938   if((cl == 1)&&(jetRadius == 0.2)&&(ptmintrack == 0.15)){
4939   fsmear->SetParameters(1,0,4.472208);
4940   }
4941   if((cl == 2)&&(jetRadius == 0.2)&&(ptmintrack == 0.15)){
4942   fsmear->SetParameters(1,0,3.543938);
4943   }
4944   if((cl == 3)&&(jetRadius == 0.2)&&(ptmintrack == 0.15)){
4945   fsmear->SetParameters(1,0,?);
4946   }
4947   if((cl == 4)&&(jetRadius == 0.2)&&(ptmintrack == 0.15)){
4948   fsmear->SetParameters(1,0,1.037476);
4949   }
4950
4951   */
4952
4953   Double_t r = fsmear.GetRandom();
4954   jetPtSmear = jetPt + r;
4955
4956   //  std::cout<<"jetPt: "<<jetPt<<std::endl;
4957   //  std::cout<<"jetPtSmear: "<<jetPtSmear<<std::endl;
4958   //  std::cout<<"r: "<<r<<std::endl;
4959
4960
4961   return jetPtSmear;
4962 }
4963
4964
4965 // _______________________________________________________________________________________________________________________
4966 AliAODJet* AliAnalysisTaskJetChem::GetMedianCluster()
4967 {
4968   // fill tracks from bckgCluster branch,
4969   // using cluster with median density (odd number of clusters)
4970   // or picking randomly one of the two closest to median (even number)
4971
4972   Int_t nBckgClusters = fBckgJetsRec->GetEntries(); // not 'recCuts': use all clusters in full eta range
4973
4974   if(nBckgClusters<3) return 0; // need at least 3 clusters (skipping 2 highest)
4975
4976   Double_t* bgrDensity = new Double_t[nBckgClusters];
4977   Int_t*    indices    = new Int_t[nBckgClusters];
4978
4979   for(Int_t ij=0; ij<nBckgClusters; ++ij){
4980
4981     AliAODJet* bgrCluster = (AliAODJet*)(fBckgJetsRec->At(ij));
4982     Double_t clusterPt    = bgrCluster->Pt();
4983     Double_t area         = bgrCluster->EffectiveAreaCharged();
4984
4985     Double_t density = 0;
4986     if(area>0) density = clusterPt/area;
4987
4988     bgrDensity[ij] = density;
4989     indices[ij]    = ij;
4990   }
4991
4992   TMath::Sort(nBckgClusters, bgrDensity, indices);
4993
4994   // get median cluster
4995
4996   AliAODJet* medianCluster = 0;
4997   Double_t   medianDensity = 0;
4998
4999   if(TMath::Odd(nBckgClusters)){
5000
5001     //Int_t medianIndex = indices[(Int_t) (0.5*(nBckgClusters-1))];
5002     Int_t medianIndex = indices[(Int_t) (0.5*(nBckgClusters+1))];
5003
5004     medianCluster = (AliAODJet*)(fBckgJetsRec->At(medianIndex));
5005
5006     Double_t clusterPt = medianCluster->Pt();
5007     Double_t area      = medianCluster->EffectiveAreaCharged();
5008
5009     if(area>0) medianDensity = clusterPt/area;
5010   }
5011   else{
5012
5013     //Int_t medianIndex1 = indices[(Int_t) (0.5*nBckgClusters-1)];
5014     //Int_t medianIndex2 = indices[(Int_t) (0.5*nBckgClusters)];
5015
5016     Int_t medianIndex1 = indices[(Int_t) (0.5*nBckgClusters)];
5017     Int_t medianIndex2 = indices[(Int_t) (0.5*nBckgClusters+1)];
5018
5019     AliAODJet* medianCluster1 = (AliAODJet*)(fBckgJetsRec->At(medianIndex1));
5020     AliAODJet* medianCluster2 = (AliAODJet*)(fBckgJetsRec->At(medianIndex2));
5021
5022     Double_t density1 = 0;
5023     Double_t clusterPt1 = medianCluster1->Pt();
5024     Double_t area1      = medianCluster1->EffectiveAreaCharged();
5025     if(area1>0) density1 = clusterPt1/area1;
5026
5027     Double_t density2 = 0;
5028     Double_t clusterPt2 = medianCluster2->Pt();
5029     Double_t area2      = medianCluster2->EffectiveAreaCharged();
5030     if(area2>0) density2 = clusterPt2/area2;
5031
5032     medianDensity = 0.5*(density1+density2);
5033
5034     medianCluster = ( (gRandom->Rndm()>0.5) ? medianCluster1 : medianCluster2 );  // select one randomly to avoid adding areas
5035   }
5036
5037   delete[] bgrDensity;
5038   delete[] indices;
5039
5040   return medianCluster;
5041 }