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