PWGJE/StrangenessInJets/AliAnalysisTaskJetChem.cxx

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