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