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