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