1 /*************************************************************************
4 * Task for Jet Chemistry Analysis in PWG4 Jet Task Force Train *
5 * Analysis of K0s, Lambda and Antilambda with and without Jetevents *
7 *************************************************************************/
9 /**************************************************************************
10 * Copyright(c) 1998-2012, ALICE Experiment at CERN, All rights reserved. *
12 * Author: The ALICE Off-line Project. *
13 * Contributors are mentioned in the code where appropriate. *
15 * Permission to use, copy, modify and distribute this software and its *
16 * documentation strictly for non-commercial purposes is hereby grante *
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 //Task for K0s, Lambda and Antilambda analysis in jets
25 //Author: Alice Zimmermann (zimmermann@physi.uni-heidelberg.de)
43 #include "THnSparse.h"
46 #include "AliAnalysisHelperJetTasks.h"
47 #include "TDatabasePDG.h"
49 #include "AliAnalysisManager.h"
50 #include "AliAODHandler.h"
51 #include "AliAODInputHandler.h"
52 #include "AliESDEvent.h"
53 #include "AliGenPythiaEventHeader.h"
54 #include "AliGenHijingEventHeader.h"
55 #include "AliGenEventHeader.h"
56 #include "TLorentzVector.h"
57 #include "AliAODEvent.h"
58 #include "AliAODJet.h"
60 #include "AliAODTrack.h"
61 #include "AliCentrality.h"
62 #include "AliAnalysisTaskSE.h"
63 #include "AliESDtrack.h"
64 #include "AliESDtrackCuts.h"
65 #include "AliESDEvent.h"
66 #include "AliESDInputHandler.h"
68 #include "AliPIDResponse.h"
69 #include "AliAODPid.h"
70 #include "AliExternalTrackParam.h"
71 #include "AliAnalysisTaskJetChem.h"
72 #include "AliPhysicsSelection.h"
73 #include "AliBackgroundSelection.h"
74 #include "AliInputEventHandler.h"
75 #include "AliAODMCHeader.h"
76 #include "AliAODPid.h"
77 #include "AliVEvent.h"
78 #include "AliAODMCParticle.h"
82 ClassImp(AliAnalysisTaskJetChem)
84 //____________________________________________________________________________
85 AliAnalysisTaskJetChem::AliAnalysisTaskJetChem()
86 : AliAnalysisTaskFragmentationFunction()
100 ,fCutV0cosPointAngle(0)
107 ,fCutDcaV0Daughters(0)
108 ,fCutDcaPosToPrimVertex(0)
109 ,fCutDcaNegToPrimVertex(0)
119 ,fFFHistosRecCutsK0Evt(0)
120 //,fFFHistosIMK0AllEvt(0)
121 //,fFFHistosIMK0Jet(0)
122 //,fFFHistosIMK0Cone(0)
126 // ,fFFHistosIMLaAllEvt(0)
127 // ,fFFHistosIMLaJet(0)
128 //,fFFHistosIMLaCone(0)
132 ,fListFeeddownLaCand(0)
133 ,fListFeeddownALaCand(0)
139 ,fListMCgenK0sCone(0)
141 ,fListMCgenALaCone(0)
142 ,IsArmenterosSelected(0)
143 // ,fFFHistosIMALaAllEvt(0)
144 // ,fFFHistosIMALaJet(0)
145 // ,fFFHistosIMALaCone(0)
161 ,fFFIMLaNBinsJetPt(0)
192 //,fh1trackPosNCls(0)
193 //,fh1trackNegNCls(0)
203 ,fh2ProperLifetimeK0sVsPtBeforeCut(0)
204 ,fh2ProperLifetimeK0sVsPtAfterCut(0)
206 ,fh1DcaV0Daughters(0)
207 ,fh1DcaPosToPrimVertex(0)
208 ,fh1DcaNegToPrimVertex(0)
209 ,fh2ArmenterosBeforeCuts(0)
210 ,fh2ArmenterosAfterCuts(0)
213 ,fh1PosDaughterCharge(0)
214 ,fh1NegDaughterCharge(0)
221 ,fhnInvMassEtaTrackPtK0s(0)
222 ,fhnInvMassEtaTrackPtLa(0)
223 ,fhnInvMassEtaTrackPtALa(0)
232 ,fh2MCEtagenK0Cone(0)
233 ,fh2MCEtagenLaCone(0)
234 ,fh2MCEtagenALaCone(0)
237 ,fh1IMALaConeSmear(0)
241 ,fhnMCrecK0ConeSmear(0)
242 ,fhnMCrecLaConeSmear(0)
243 ,fhnMCrecALaConeSmear(0)
244 ,fhnK0sSecContinCone(0)
245 ,fhnLaSecContinCone(0)
246 ,fhnALaSecContinCone(0)
268 ,fh1MCMultiplicityPrimary(0)
269 ,fh1MCMultiplicityTracks(0)
272 ,fhnFeedDownLaCone(0)
273 ,fhnFeedDownALaCone(0)
274 ,fh1MCProdRadiusK0s(0)
275 ,fh1MCProdRadiusLambda(0)
276 ,fh1MCProdRadiusAntiLambda(0)
280 ,fh1MCPtAntiLambda(0)
288 ,fh1MCRapAntiLambda(0)
292 ,fh1MCEtaAntiLambda(0)
295 // default constructor
298 //__________________________________________________________________________________________
299 AliAnalysisTaskJetChem::AliAnalysisTaskJetChem(const char *name)
300 : AliAnalysisTaskFragmentationFunction(name)
314 ,fCutV0cosPointAngle(0)
321 ,fCutDcaV0Daughters(0)
322 ,fCutDcaPosToPrimVertex(0)
323 ,fCutDcaNegToPrimVertex(0)
333 ,fFFHistosRecCutsK0Evt(0)
334 //,fFFHistosIMK0AllEvt(0)
335 //,fFFHistosIMK0Jet(0)
336 //,fFFHistosIMK0Cone(0)
340 //,fFFHistosIMLaAllEvt(0)
341 //,fFFHistosIMLaJet(0)
342 //,fFFHistosIMLaCone(0)
346 ,fListFeeddownLaCand(0)
347 ,fListFeeddownALaCand(0)
353 ,fListMCgenK0sCone(0)
355 ,fListMCgenALaCone(0)
356 ,IsArmenterosSelected(0)
357 //,fFFHistosIMALaAllEvt(0)
358 //,fFFHistosIMALaJet(0)
359 // ,fFFHistosIMALaCone(0)
375 ,fFFIMLaNBinsJetPt(0)
406 // ,fh1trackPosNCls(0)
407 // ,fh1trackNegNCls(0)
417 ,fh2ProperLifetimeK0sVsPtBeforeCut(0)
418 ,fh2ProperLifetimeK0sVsPtAfterCut(0)
420 ,fh1DcaV0Daughters(0)
421 ,fh1DcaPosToPrimVertex(0)
422 ,fh1DcaNegToPrimVertex(0)
423 ,fh2ArmenterosBeforeCuts(0)
424 ,fh2ArmenterosAfterCuts(0)
427 ,fh1PosDaughterCharge(0)
428 ,fh1NegDaughterCharge(0)
435 ,fhnInvMassEtaTrackPtK0s(0)
436 ,fhnInvMassEtaTrackPtLa(0)
437 ,fhnInvMassEtaTrackPtALa(0)
446 ,fh2MCEtagenK0Cone(0)
447 ,fh2MCEtagenLaCone(0)
448 ,fh2MCEtagenALaCone(0)
451 ,fh1IMALaConeSmear(0)
455 ,fhnMCrecK0ConeSmear(0)
456 ,fhnMCrecLaConeSmear(0)
457 ,fhnMCrecALaConeSmear(0)
458 ,fhnK0sSecContinCone(0)
459 ,fhnLaSecContinCone(0)
460 ,fhnALaSecContinCone(0)
482 ,fh1MCMultiplicityPrimary(0)
483 ,fh1MCMultiplicityTracks(0)
486 ,fhnFeedDownLaCone(0)
487 ,fhnFeedDownALaCone(0)
488 ,fh1MCProdRadiusK0s(0)
489 ,fh1MCProdRadiusLambda(0)
490 ,fh1MCProdRadiusAntiLambda(0)
494 ,fh1MCPtAntiLambda(0)
502 ,fh1MCRapAntiLambda(0)
506 ,fh1MCEtaAntiLambda(0)
512 DefineOutput(1,TList::Class());
515 //__________________________________________________________________________________________________________________________
516 AliAnalysisTaskJetChem::AliAnalysisTaskJetChem(const AliAnalysisTaskJetChem ©)
517 : AliAnalysisTaskFragmentationFunction()
519 ,fRandom(copy.fRandom)
520 ,fAnalysisMC(copy.fAnalysisMC)
521 ,fDeltaVertexZ(copy.fDeltaVertexZ)
522 ,fCutjetEta(copy.fCutjetEta)
523 ,fCuttrackNegNcls(copy.fCuttrackNegNcls)
524 ,fCuttrackPosNcls(copy.fCuttrackPosNcls)
525 ,fCutPostrackRap(copy.fCutPostrackRap)
526 ,fCutNegtrackRap(copy.fCutNegtrackRap)
527 ,fCutRap(copy.fCutRap)
528 ,fCutPostrackEta(copy.fCutPostrackEta)
529 ,fCutNegtrackEta(copy.fCutNegtrackEta)
530 ,fCutEta(copy.fCutEta)
531 ,fCutV0cosPointAngle(copy.fCutV0cosPointAngle)
532 ,fKinkDaughters(copy.fKinkDaughters)
533 ,fRequireTPCRefit(copy.fRequireTPCRefit)
534 ,fCutArmenteros(copy.fCutArmenteros)
535 ,fCutV0DecayMin(copy.fCutV0DecayMin)
536 ,fCutV0DecayMax(copy.fCutV0DecayMax)
537 ,fCutV0totMom(copy.fCutV0totMom)
538 ,fCutDcaV0Daughters(copy.fCutDcaV0Daughters)
539 ,fCutDcaPosToPrimVertex(copy.fCutDcaPosToPrimVertex)
540 ,fCutDcaNegToPrimVertex(copy.fCutDcaNegToPrimVertex)
541 ,fCutV0RadiusMin(copy.fCutV0RadiusMin)
542 ,fCutV0RadiusMax(copy.fCutV0RadiusMax)
543 ,fCutBetheBloch(copy.fCutBetheBloch)
544 ,fCutRatio(copy.fCutRatio)
545 ,fK0Type(copy.fK0Type)
546 ,fFilterMaskK0(copy.fFilterMaskK0)
547 ,fListK0s(copy.fListK0s)
548 ,fPIDResponse(copy.fPIDResponse)
549 ,fV0QAK0(copy.fV0QAK0)
550 ,fFFHistosRecCutsK0Evt(copy.fFFHistosRecCutsK0Evt)
551 //,fFFHistosIMK0AllEvt(copy.fFFHistosIMK0AllEvt)
552 //,fFFHistosIMK0Jet(copy.fFFHistosIMK0Jet)
553 //,fFFHistosIMK0Cone(copy.fFFHistosIMK0Cone)
554 ,fLaType(copy.fLaType)
555 ,fFilterMaskLa(copy.fFilterMaskLa)
556 ,fListLa(copy.fListLa)
557 //,fFFHistosIMLaAllEvt(copy.fFFHistosIMLaAllEvt)
558 //,fFFHistosIMLaJet(copy.fFFHistosIMLaJet)
559 //,fFFHistosIMLaCone(copy.fFFHistosIMLaCone)
560 ,fALaType(copy.fALaType)
561 ,fFilterMaskALa(copy.fFilterMaskALa)
562 ,fListALa(copy.fListALa)
563 ,fListFeeddownLaCand(copy.fListFeeddownLaCand)
564 ,fListFeeddownALaCand(copy.fListFeeddownALaCand)
565 ,jetConeFDLalist(copy.jetConeFDLalist)
566 ,jetConeFDALalist(copy.jetConeFDALalist)
567 ,fListMCgenK0s(copy.fListMCgenK0s)
568 ,fListMCgenLa(copy.fListMCgenLa)
569 ,fListMCgenALa(copy.fListMCgenALa)
570 ,fListMCgenK0sCone(copy.fListMCgenK0sCone)
571 ,fListMCgenLaCone(copy.fListMCgenLaCone)
572 ,fListMCgenALaCone(copy.fListMCgenALaCone)
573 ,IsArmenterosSelected(copy.IsArmenterosSelected)
574 //,fFFHistosIMALaAllEvt(copy.fFFHistosIMALaAllEvt)
575 //,fFFHistosIMALaJet(copy.fFFHistosIMALaJet)
576 //,fFFHistosIMALaCone(copy.fFFHistosIMALaCone)
577 ,fFFIMNBinsJetPt(copy.fFFIMNBinsJetPt)
578 ,fFFIMJetPtMin(copy.fFFIMJetPtMin)
579 ,fFFIMJetPtMax(copy.fFFIMJetPtMax)
580 ,fFFIMNBinsInvM(copy.fFFIMNBinsInvM)
581 ,fFFIMInvMMin(copy.fFFIMInvMMin)
582 ,fFFIMInvMMax(copy.fFFIMInvMMax)
583 ,fFFIMNBinsPt(copy.fFFIMNBinsPt)
584 ,fFFIMPtMin(copy.fFFIMPtMin)
585 ,fFFIMPtMax(copy.fFFIMPtMax)
586 ,fFFIMNBinsXi(copy.fFFIMNBinsXi)
587 ,fFFIMXiMin(copy.fFFIMXiMin)
588 ,fFFIMXiMax(copy.fFFIMXiMax)
589 ,fFFIMNBinsZ(copy.fFFIMNBinsZ)
590 ,fFFIMZMin(copy.fFFIMZMin)
591 ,fFFIMZMax(copy.fFFIMZMax)
592 ,fFFIMLaNBinsJetPt(copy.fFFIMLaNBinsJetPt)
593 ,fFFIMLaJetPtMin(copy.fFFIMLaJetPtMin)
594 ,fFFIMLaJetPtMax(copy.fFFIMLaJetPtMax)
595 ,fFFIMLaNBinsInvM(copy.fFFIMLaNBinsInvM)
596 ,fFFIMLaInvMMin(copy.fFFIMLaInvMMin)
597 ,fFFIMLaInvMMax(copy.fFFIMLaInvMMax)
598 ,fFFIMLaNBinsPt(copy.fFFIMLaNBinsPt)
599 ,fFFIMLaPtMin(copy.fFFIMLaPtMin)
600 ,fFFIMLaPtMax(copy.fFFIMLaPtMax)
601 ,fFFIMLaNBinsXi(copy.fFFIMLaNBinsXi)
602 ,fFFIMLaXiMin(copy.fFFIMLaXiMin)
603 ,fFFIMLaXiMax(copy.fFFIMLaXiMax)
604 ,fFFIMLaNBinsZ(copy.fFFIMLaNBinsZ)
605 ,fFFIMLaZMin(copy.fFFIMLaZMin)
606 ,fFFIMLaZMax(copy.fFFIMLaZMax)
607 ,fh1EvtAllCent(copy.fh1EvtAllCent)
609 ,fh1K0Mult(copy.fh1K0Mult)
610 ,fh1dPhiJetK0(copy.fh1dPhiJetK0)
611 ,fh1LaMult(copy.fh1LaMult)
612 ,fh1dPhiJetLa(copy.fh1dPhiJetLa)
613 ,fh1ALaMult(copy.fh1ALaMult)
614 ,fh1dPhiJetALa(copy.fh1dPhiJetALa)
615 ,fh1JetEta(copy.fh1JetEta)
616 ,fh1JetPhi(copy.fh1JetPhi)
617 ,fh2JetEtaPhi(copy.fh2JetEtaPhi)
618 //,fh1V0JetPt(copy.fh1V0JetPt)
619 ,fh1IMK0Cone(copy.fh1IMK0Cone)
620 ,fh1IMLaCone(copy.fh1IMLaCone)
621 ,fh1IMALaCone(copy.fh1IMALaCone)
622 ,fh2FFJetTrackEta(copy.fh2FFJetTrackEta)
623 //,fh1trackPosNCls(copy.fh1trackPosNCls)
624 //,fh1trackNegNCls(copy.fh1trackNegNCls)
625 ,fh1trackPosRap(copy.fh1trackPosRap)
626 ,fh1trackNegRap(copy.fh1trackNegRap)
627 //,fh1V0Rap(copy.fh1V0Rap)
628 ,fh1trackPosEta(copy.fh1trackPosEta)
629 ,fh1trackNegEta(copy.fh1trackNegEta)
630 ,fh1V0Eta(copy.fh1V0Eta)
631 //,fh1V0totMom(copy.fh1V0totMom)
632 ,fh1CosPointAngle(copy.fh1CosPointAngle)
633 ,fh1DecayLengthV0(copy.fh1DecayLengthV0)
634 ,fh2ProperLifetimeK0sVsPtBeforeCut(copy.fh2ProperLifetimeK0sVsPtBeforeCut)
635 ,fh2ProperLifetimeK0sVsPtAfterCut(copy.fh2ProperLifetimeK0sVsPtAfterCut)
636 ,fh1V0Radius(copy.fh1V0Radius)
637 ,fh1DcaV0Daughters(copy.fh1DcaV0Daughters)
638 ,fh1DcaPosToPrimVertex(copy.fh1DcaPosToPrimVertex)
639 ,fh1DcaNegToPrimVertex(copy.fh1DcaNegToPrimVertex)
640 ,fh2ArmenterosBeforeCuts(copy.fh2ArmenterosBeforeCuts)
641 ,fh2ArmenterosAfterCuts(copy.fh2ArmenterosAfterCuts)
642 ,fh2BBLaPos(copy.fh2BBLaPos)
643 ,fh2BBLaNeg(copy.fh2BBLaPos)
644 ,fh1PosDaughterCharge(copy.fh1PosDaughterCharge)
645 ,fh1NegDaughterCharge(copy.fh1NegDaughterCharge)
646 ,fh1PtMCK0s(copy.fh1PtMCK0s)
647 ,fh1PtMCLa(copy.fh1PtMCLa)
648 ,fh1PtMCALa(copy.fh1PtMCALa)
649 ,fh1EtaK0s(copy.fh1EtaK0s)
650 ,fh1EtaLa(copy.fh1EtaLa)
651 ,fh1EtaALa(copy.fh1EtaALa)
652 ,fhnInvMassEtaTrackPtK0s(copy.fhnInvMassEtaTrackPtK0s)
653 ,fhnInvMassEtaTrackPtLa(copy.fhnInvMassEtaTrackPtLa)
654 ,fhnInvMassEtaTrackPtALa(copy.fhnInvMassEtaTrackPtALa)
655 ,fh1TrackMultCone(copy.fh1TrackMultCone)
656 ,fh2TrackMultCone(copy.fh2TrackMultCone)
657 ,fh2NJK0(copy.fh2NJK0)
658 ,fh2NJLa(copy.fh2NJLa)
659 ,fh2NJALa(copy.fh2NJALa)
660 ,fh2MCgenK0Cone(copy.fh2MCgenK0Cone)
661 ,fh2MCgenLaCone(copy.fh2MCgenLaCone)
662 ,fh2MCgenALaCone(copy.fh2MCgenALaCone)
663 ,fh2MCEtagenK0Cone(copy.fh2MCEtagenK0Cone)
664 ,fh2MCEtagenLaCone(copy.fh2MCEtagenLaCone)
665 ,fh2MCEtagenALaCone(copy.fh2MCEtagenALaCone)
666 ,fh1IMK0ConeSmear(copy.fh1IMK0ConeSmear)
667 ,fh1IMLaConeSmear(copy.fh1IMLaConeSmear)
668 ,fh1IMALaConeSmear(copy.fh1IMALaConeSmear)
669 ,fhnMCrecK0Cone(copy.fhnMCrecK0Cone)
670 ,fhnMCrecLaCone(copy.fhnMCrecLaCone)
671 ,fhnMCrecALaCone(copy.fhnMCrecALaCone)
672 ,fhnMCrecK0ConeSmear(copy.fhnMCrecK0ConeSmear)
673 ,fhnMCrecLaConeSmear(copy.fhnMCrecLaConeSmear)
674 ,fhnMCrecALaConeSmear(copy.fhnMCrecALaConeSmear)
675 ,fhnK0sSecContinCone(copy.fhnK0sSecContinCone)
676 ,fhnLaSecContinCone(copy.fhnLaSecContinCone)
677 ,fhnALaSecContinCone(copy.fhnALaSecContinCone)
678 ,fhnK0sIncl(copy.fhnK0sIncl)
679 ,fhnK0sCone(copy.fhnK0sCone)
680 ,fhnLaIncl(copy.fhnLaIncl)
681 ,fhnLaCone(copy.fhnLaCone)
682 ,fhnALaIncl(copy.fhnALaIncl)
683 ,fhnALaCone(copy.fhnALaCone)
684 ,fhnK0sPC(copy.fhnK0sPC)
685 ,fhnLaPC(copy.fhnLaPC)
686 ,fhnALaPC(copy.fhnALaPC)
687 ,fhnK0sMCC(copy.fhnK0sMCC)
688 ,fhnLaMCC(copy.fhnLaMCC)
689 ,fhnALaMCC(copy.fhnALaMCC)
690 ,fhnK0sRC(copy.fhnK0sRC)
691 ,fhnLaRC(copy.fhnLaRC)
692 ,fhnALaRC(copy.fhnALaRC)
693 ,fhnK0sOC(copy.fhnK0sOC)
694 ,fhnLaOC(copy.fhnLaOC)
695 ,fhnALaOC(copy.fhnALaOC)
696 ,fh1AreaExcluded(copy.fh1AreaExcluded)
697 ,fh1MedianEta(copy.fh1MedianEta)
698 ,fh1JetPtMedian(copy.fh1JetPtMedian)
699 ,fh1MCMultiplicityPrimary(copy.fh1MCMultiplicityPrimary)
700 ,fh1MCMultiplicityTracks(copy.fh1MCMultiplicityTracks)
701 ,fhnFeedDownLa(copy.fhnFeedDownLa)
702 ,fhnFeedDownALa(copy.fhnFeedDownALa)
703 ,fhnFeedDownLaCone(copy.fhnFeedDownLaCone)
704 ,fhnFeedDownALaCone(copy.fhnFeedDownALaCone)
705 ,fh1MCProdRadiusK0s(copy.fh1MCProdRadiusK0s)
706 ,fh1MCProdRadiusLambda(copy.fh1MCProdRadiusLambda)
707 ,fh1MCProdRadiusAntiLambda(copy.fh1MCProdRadiusAntiLambda)
708 ,fh1MCPtV0s(copy.fh1MCPtV0s)
709 ,fh1MCPtK0s(copy.fh1MCPtK0s)
710 ,fh1MCPtLambda(copy.fh1MCPtLambda)
711 ,fh1MCPtAntiLambda(copy.fh1MCPtAntiLambda)
712 ,fh1MCXiPt(copy.fh1MCXiPt)
713 ,fh1MCXibarPt(copy.fh1MCXibarPt)
714 ,fh2MCEtaVsPtK0s(copy.fh2MCEtaVsPtK0s)
715 ,fh2MCEtaVsPtLa(copy.fh2MCEtaVsPtLa)
716 ,fh2MCEtaVsPtALa(copy.fh2MCEtaVsPtALa)
717 ,fh1MCRapK0s(copy.fh1MCRapK0s)
718 ,fh1MCRapLambda(copy.fh1MCRapLambda)
719 ,fh1MCRapAntiLambda(copy.fh1MCRapAntiLambda)
720 ,fh1MCEtaAllK0s(copy.fh1MCEtaAllK0s)
721 ,fh1MCEtaK0s(copy.fh1MCEtaK0s)
722 ,fh1MCEtaLambda(copy.fh1MCEtaLambda)
723 ,fh1MCEtaAntiLambda(copy.fh1MCEtaAntiLambda)
730 // _________________________________________________________________________________________________________________________________
731 AliAnalysisTaskJetChem& AliAnalysisTaskJetChem::operator=(const AliAnalysisTaskJetChem& o)
736 AliAnalysisTaskFragmentationFunction::operator=(o);
739 fAnalysisMC = o.fAnalysisMC;
740 fDeltaVertexZ = o.fDeltaVertexZ;
741 fCutjetEta = o.fCutjetEta;
742 fCuttrackNegNcls = o.fCuttrackNegNcls;
743 fCuttrackPosNcls = o.fCuttrackPosNcls;
744 fCutPostrackRap = o.fCutPostrackRap;
745 fCutNegtrackRap = o.fCutNegtrackRap;
747 fCutPostrackEta = o.fCutPostrackEta;
748 fCutNegtrackEta = o.fCutNegtrackEta;
750 fCutV0cosPointAngle = o.fCutV0cosPointAngle;
751 fKinkDaughters = o.fKinkDaughters;
752 fRequireTPCRefit = o.fRequireTPCRefit;
753 fCutArmenteros = o.fCutArmenteros;
754 fCutV0DecayMin = o.fCutV0DecayMin;
755 fCutV0DecayMax = o.fCutV0DecayMax;
756 fCutV0totMom = o.fCutV0totMom;
757 fCutDcaV0Daughters = o.fCutDcaV0Daughters;
758 fCutDcaPosToPrimVertex = o.fCutDcaPosToPrimVertex;
759 fCutDcaNegToPrimVertex = o.fCutDcaNegToPrimVertex;
760 fCutV0RadiusMin = o.fCutV0RadiusMin;
761 fCutV0RadiusMax = o.fCutV0RadiusMax;
762 fCutBetheBloch = o.fCutBetheBloch;
763 fCutRatio = o.fCutRatio;
765 fFilterMaskK0 = o.fFilterMaskK0;
766 fListK0s = o.fListK0s;
767 fPIDResponse = o.fPIDResponse;
769 fFFHistosRecCutsK0Evt = o.fFFHistosRecCutsK0Evt;
770 //fFFHistosIMK0AllEvt = o.fFFHistosIMK0AllEvt;
771 //fFFHistosIMK0Jet = o.fFFHistosIMK0Jet;
772 //fFFHistosIMK0Cone = o.fFFHistosIMK0Cone;
774 fFilterMaskLa = o.fFilterMaskLa;
776 //fFFHistosIMLaAllEvt = o.fFFHistosIMLaAllEvt;
777 //fFFHistosIMLaJet = o.fFFHistosIMLaJet;
778 //fFFHistosIMLaCone = o.fFFHistosIMLaCone;
779 fALaType = o.fALaType;
780 fFilterMaskALa = o.fFilterMaskALa;
781 fListFeeddownLaCand = o.fListFeeddownLaCand;
782 fListFeeddownALaCand = o.fListFeeddownALaCand;
783 jetConeFDLalist = o.jetConeFDLalist;
784 jetConeFDALalist = o.jetConeFDALalist;
785 fListMCgenK0s = o.fListMCgenK0s;
786 fListMCgenLa = o.fListMCgenLa;
787 fListMCgenALa = o.fListMCgenALa;
788 fListMCgenK0sCone = o.fListMCgenK0sCone;
789 fListMCgenLaCone = o.fListMCgenLaCone;
790 fListMCgenALaCone = o.fListMCgenALaCone;
791 IsArmenterosSelected = o.IsArmenterosSelected;
792 // fFFHistosIMALaAllEvt = o.fFFHistosIMALaAllEvt;
793 // fFFHistosIMALaJet = o.fFFHistosIMALaJet;
794 // fFFHistosIMALaCone = o.fFFHistosIMALaCone;
795 fFFIMNBinsJetPt = o.fFFIMNBinsJetPt;
796 fFFIMJetPtMin = o.fFFIMJetPtMin;
797 fFFIMJetPtMax = o.fFFIMJetPtMax;
798 fFFIMNBinsPt = o.fFFIMNBinsPt;
799 fFFIMPtMin = o.fFFIMPtMin;
800 fFFIMPtMax = o.fFFIMPtMax;
801 fFFIMNBinsXi = o.fFFIMNBinsXi;
802 fFFIMXiMin = o.fFFIMXiMin;
803 fFFIMXiMax = o.fFFIMXiMax;
804 fFFIMNBinsZ = o.fFFIMNBinsZ;
805 fFFIMZMin = o.fFFIMZMin;
806 fFFIMZMax = o.fFFIMZMax;
807 fFFIMLaNBinsJetPt = o.fFFIMLaNBinsJetPt;
808 fFFIMLaJetPtMin = o.fFFIMLaJetPtMin;
809 fFFIMLaJetPtMax = o.fFFIMLaJetPtMax;
810 fFFIMLaNBinsPt = o.fFFIMLaNBinsPt;
811 fFFIMLaPtMin = o.fFFIMLaPtMin;
812 fFFIMLaPtMax = o.fFFIMLaPtMax;
813 fFFIMLaNBinsXi = o.fFFIMLaNBinsXi;
814 fFFIMLaXiMin = o.fFFIMLaXiMin;
815 fFFIMLaXiMax = o.fFFIMLaXiMax;
816 fFFIMLaNBinsZ = o.fFFIMLaNBinsZ;
817 fFFIMLaZMin = o.fFFIMLaZMin;
818 fFFIMLaZMax = o.fFFIMLaZMax;
819 fh1EvtAllCent = o.fh1EvtAllCent;
821 fh1K0Mult = o.fh1K0Mult;
822 fh1dPhiJetK0 = o.fh1dPhiJetK0;
823 fh1LaMult = o.fh1LaMult;
824 fh1dPhiJetLa = o.fh1dPhiJetLa;
825 fh1ALaMult = o.fh1ALaMult;
826 fh1dPhiJetALa = o.fh1dPhiJetALa;
827 fh1JetEta = o.fh1JetEta;
828 fh1JetPhi = o.fh1JetPhi;
829 fh2JetEtaPhi = o.fh2JetEtaPhi;
830 //fh1V0JetPt = o.fh1V0JetPt;
831 fh1IMK0Cone = o.fh1IMK0Cone;
832 fh1IMLaCone = o.fh1IMLaCone;
833 fh1IMALaCone = o.fh1IMALaCone;
834 fh2FFJetTrackEta = o.fh2FFJetTrackEta;
835 //fh1trackPosNCls = o.fh1trackPosNCls;
836 //fh1trackNegNCls = o.fh1trackNegNCls;
837 fh1trackPosRap = o.fh1trackPosRap;
838 fh1trackNegRap = o.fh1trackNegRap;
839 //fh1V0Rap = o.fh1V0Rap;
840 fh1trackPosEta = o.fh1trackPosEta;
841 fh1trackNegEta = o.fh1trackNegEta;
842 fh1V0Eta = o.fh1V0Eta;
843 // fh1V0totMom = o.fh1V0totMom;
844 fh1CosPointAngle = o.fh1CosPointAngle;
845 fh1DecayLengthV0 = o.fh1DecayLengthV0;
846 fh2ProperLifetimeK0sVsPtBeforeCut = o.fh2ProperLifetimeK0sVsPtBeforeCut;
847 fh2ProperLifetimeK0sVsPtAfterCut= o.fh2ProperLifetimeK0sVsPtAfterCut;
848 fh1V0Radius = o.fh1V0Radius;
849 fh1DcaV0Daughters = o.fh1DcaV0Daughters;
850 fh1DcaPosToPrimVertex = o.fh1DcaPosToPrimVertex;
851 fh1DcaNegToPrimVertex = o.fh1DcaNegToPrimVertex;
852 fh2ArmenterosBeforeCuts = o.fh2ArmenterosBeforeCuts;
853 fh2ArmenterosAfterCuts = o.fh2ArmenterosAfterCuts;
854 fh2BBLaPos = o.fh2BBLaPos;
855 fh2BBLaNeg = o.fh2BBLaPos;
856 fh1PosDaughterCharge = o.fh1PosDaughterCharge;
857 fh1NegDaughterCharge = o.fh1NegDaughterCharge;
858 fh1PtMCK0s = o.fh1PtMCK0s;
859 fh1PtMCLa = o.fh1PtMCLa;
860 fh1PtMCALa = o.fh1PtMCALa;
861 fh1EtaK0s = o.fh1EtaK0s;
862 fh1EtaLa = o.fh1EtaLa;
863 fh1EtaALa = o.fh1EtaALa;
864 fhnInvMassEtaTrackPtK0s = o.fhnInvMassEtaTrackPtK0s;
865 fhnInvMassEtaTrackPtLa = o.fhnInvMassEtaTrackPtLa;
866 fhnInvMassEtaTrackPtALa = o.fhnInvMassEtaTrackPtALa;
867 fh1TrackMultCone = o.fh1TrackMultCone;
868 fh2TrackMultCone = o.fh2TrackMultCone;
871 fh2NJALa = o.fh2NJALa;
872 fh2MCgenK0Cone = o.fh2MCgenK0Cone;
873 fh2MCgenLaCone = o.fh2MCgenLaCone;
874 fh2MCgenALaCone = o.fh2MCgenALaCone;
875 fh2MCEtagenK0Cone = o.fh2MCEtagenK0Cone;
876 fh2MCEtagenLaCone = o.fh2MCEtagenLaCone;
877 fh2MCEtagenALaCone = o.fh2MCEtagenALaCone;
878 fh1IMK0ConeSmear = o.fh1IMK0ConeSmear;
879 fh1IMLaConeSmear = o.fh1IMLaConeSmear;
880 fh1IMALaConeSmear = o.fh1IMALaConeSmear;
881 fhnMCrecK0Cone = o.fhnMCrecK0Cone;
882 fhnMCrecLaCone = o.fhnMCrecLaCone;
883 fhnMCrecALaCone = o.fhnMCrecALaCone;
884 fhnMCrecK0ConeSmear = o.fhnMCrecK0ConeSmear;
885 fhnMCrecLaConeSmear = o.fhnMCrecLaConeSmear;
886 fhnMCrecALaConeSmear = o.fhnMCrecALaConeSmear;
887 fhnK0sSecContinCone = o.fhnK0sSecContinCone;
888 fhnLaSecContinCone = o.fhnLaSecContinCone;
889 fhnALaSecContinCone = o.fhnALaSecContinCone;
890 fhnK0sIncl = o.fhnK0sIncl;
891 fhnK0sCone = o.fhnK0sCone;
892 fhnLaIncl = o.fhnLaIncl;
893 fhnLaCone = o.fhnLaCone;
894 fhnALaIncl = o.fhnALaIncl;
895 fhnALaCone = o.fhnALaCone;
896 fhnK0sPC = o.fhnK0sPC;
898 fhnALaPC = o.fhnALaPC;
899 fhnK0sRC = o.fhnK0sRC;
901 fhnALaRC = o.fhnALaRC;
902 fhnK0sOC = o.fhnK0sOC;
904 fhnALaOC = o.fhnALaOC;
905 fh1AreaExcluded = o.fh1AreaExcluded;
906 fh1MedianEta = o.fh1MedianEta;
907 fh1JetPtMedian = o.fh1JetPtMedian;
908 fh1MCMultiplicityPrimary = o.fh1MCMultiplicityPrimary;
909 fh1MCMultiplicityTracks = o.fh1MCMultiplicityTracks;
910 fhnFeedDownLa = o.fhnFeedDownLa;
911 fhnFeedDownALa = o.fhnFeedDownALa;
912 fhnFeedDownLaCone = o.fhnFeedDownLaCone;
913 fhnFeedDownALaCone = o.fhnFeedDownALaCone;
914 fh1MCProdRadiusK0s = o.fh1MCProdRadiusK0s;
915 fh1MCProdRadiusLambda = o.fh1MCProdRadiusLambda;
916 fh1MCProdRadiusAntiLambda = o.fh1MCProdRadiusAntiLambda;
917 fh1MCPtV0s = o.fh1MCPtV0s;
918 fh1MCPtK0s = o.fh1MCPtK0s;
919 fh1MCPtLambda = o.fh1MCPtLambda;
920 fh1MCPtAntiLambda = o.fh1MCPtAntiLambda;
921 fh1MCXiPt = o.fh1MCXiPt;
922 fh1MCXibarPt = o.fh1MCXibarPt;
923 fh2MCEtaVsPtK0s = o.fh2MCEtaVsPtK0s;
924 fh2MCEtaVsPtLa = o.fh2MCEtaVsPtLa;
925 fh2MCEtaVsPtALa = o.fh2MCEtaVsPtALa;
926 fh1MCRapK0s = o.fh1MCRapK0s;
927 fh1MCRapLambda = o.fh1MCRapLambda;
928 fh1MCRapAntiLambda = o.fh1MCRapAntiLambda;
929 fh1MCEtaAllK0s = o.fh1MCEtaAllK0s;
930 fh1MCEtaK0s = o.fh1MCEtaK0s;
931 fh1MCEtaLambda = o.fh1MCEtaLambda;
932 fh1MCEtaAntiLambda = o.fh1MCEtaAntiLambda;
938 //_______________________________________________
939 AliAnalysisTaskJetChem::~AliAnalysisTaskJetChem()
944 if(fListK0s) delete fListK0s;
945 if(fListLa) delete fListLa;
946 if(fListALa) delete fListALa;
947 if(fListFeeddownLaCand) delete fListFeeddownLaCand;
948 if(fListFeeddownALaCand) delete fListFeeddownALaCand;
949 if(jetConeFDLalist) delete jetConeFDLalist;
950 if(jetConeFDALalist) delete jetConeFDALalist;
951 if(fListMCgenK0s) delete fListMCgenK0s;
952 if(fListMCgenLa) delete fListMCgenLa;
953 if(fListMCgenALa) delete fListMCgenALa;
954 if(fRandom) delete fRandom;
957 //________________________________________________________________________________________________________________________________
958 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AliFragFuncHistosInvMass(const char* name,
959 Int_t nJetPt, Float_t jetPtMin, Float_t jetPtMax,
960 Int_t nInvMass, Float_t invMassMin, Float_t invMassMax,
961 Int_t nPt, Float_t ptMin, Float_t ptMax,
962 Int_t nXi, Float_t xiMin, Float_t xiMax,
963 Int_t nZ , Float_t zMin , Float_t zMax )
968 ,fNBinsInvMass(nInvMass)
969 ,fInvMassMin(invMassMin)
970 ,fInvMassMax(invMassMax)
986 // default constructor
990 //______________________________________________________________________________________________________________
991 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AliFragFuncHistosInvMass(const AliFragFuncHistosInvMass& copy)
993 ,fNBinsJetPt(copy.fNBinsJetPt)
994 ,fJetPtMin(copy.fJetPtMin)
995 ,fJetPtMax(copy.fJetPtMax)
996 ,fNBinsInvMass(copy.fNBinsInvMass)
997 ,fInvMassMin(copy.fInvMassMin)
998 ,fInvMassMax(copy.fInvMassMax)
999 ,fNBinsPt(copy.fNBinsPt)
1000 ,fPtMin(copy.fPtMin)
1001 ,fPtMax(copy.fPtMax)
1002 ,fNBinsXi(copy.fNBinsXi)
1003 ,fXiMin(copy.fXiMin)
1004 ,fXiMax(copy.fXiMax)
1005 ,fNBinsZ(copy.fNBinsZ)
1008 ,fh3TrackPt(copy.fh3TrackPt)
1011 ,fh1JetPt(copy.fh1JetPt)
1012 ,fNameFF(copy.fNameFF)
1017 //______________________________________________________________________________________________________________________________________________________________________
1018 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass& AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::operator=(const AliAnalysisTaskJetChem::AliFragFuncHistosInvMass& o)
1023 TObject::operator=(o);
1024 fNBinsJetPt = o.fNBinsJetPt;
1025 fJetPtMin = o.fJetPtMin;
1026 fJetPtMax = o.fJetPtMax;
1027 fNBinsInvMass = o.fNBinsInvMass;
1028 fInvMassMin = o.fInvMassMin;
1029 fInvMassMax = o.fInvMassMax;
1030 fNBinsPt = o.fNBinsPt;
1033 fNBinsXi = o.fNBinsXi;
1036 fNBinsZ = o.fNBinsZ;
1039 fh3TrackPt = o.fh3TrackPt;
1042 fh1JetPt = o.fh1JetPt;
1043 fNameFF = o.fNameFF;
1049 //___________________________________________________________________________
1050 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::~AliFragFuncHistosInvMass()
1054 if(fh1JetPt) delete fh1JetPt;
1055 if(fh3TrackPt) delete fh3TrackPt;
1056 if(fh3Xi) delete fh3Xi;
1057 if(fh3Z) delete fh3Z;
1060 //_________________________________________________________________
1061 void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::DefineHistos()
1065 fh1JetPt = new TH1F(Form("fh1FFJetPtIM%s", fNameFF.Data()),"",fNBinsJetPt,fJetPtMin,fJetPtMax);
1066 fh3TrackPt = new TH3F(Form("fh3FFTrackPtIM%s",fNameFF.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsPt, fPtMin, fPtMax);
1067 fh3Xi = new TH3F(Form("fh3FFXiIM%s", fNameFF.Data()),"", fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsXi, fXiMin, fXiMax);
1068 fh3Z = new TH3F(Form("fh3FFZIM%s", fNameFF.Data()),"", fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsZ, fZMin, fZMax);
1070 AliAnalysisTaskFragmentationFunction::SetProperties(fh1JetPt, "p_{t} (GeV/c)", "entries");
1071 AliAnalysisTaskJetChem::SetProperties(fh3TrackPt,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","p_{t} (GeV/c)");
1072 AliAnalysisTaskJetChem::SetProperties(fh3Xi,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","#xi");
1073 AliAnalysisTaskJetChem::SetProperties(fh3Z,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","z");
1076 //________________________________________________________________________________________________________________________________
1077 void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::FillFF(Float_t trackPt, Float_t invM, Float_t jetPt, Bool_t incrementJetPt)
1079 // fill FF, don't use TH3F anymore use THnSparse instead to save memory
1081 if(incrementJetPt) fh1JetPt->Fill(jetPt);
1082 //fh3TrackPt->Fill(jetPt,invM,trackPt);//Fill(x,y,z)
1085 if(jetPt>0) z = trackPt / jetPt;
1087 if(z>0) xi = TMath::Log(1/z);
1089 //fh3Xi->Fill(jetPt,invM,xi);
1090 //fh3Z->Fill(jetPt,invM,z);
1093 //___________________________________________________________________________________
1094 void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AddToOutput(TList* list) const
1096 // add histos to list
1098 list->Add(fh1JetPt);
1099 //list->Add(fh3TrackPt);
1105 //____________________________________________________
1106 void AliAnalysisTaskJetChem::UserCreateOutputObjects()
1108 // create output objects
1110 fRandom = new TRandom3(0);
1111 fRandom->SetSeed(0);
1113 if(fDebug > 1) Printf("AliAnalysisTaskJetChem::UserCreateOutputObjects()");
1115 // create list of tracks and jets
1117 fTracksRecCuts = new TList();
1118 fTracksRecCuts->SetOwner(kFALSE); //objects in TList wont be deleted when TList is deleted
1119 fJetsRecCuts = new TList();
1120 fJetsRecCuts->SetOwner(kFALSE);
1121 fBckgJetsRec = new TList();
1122 fBckgJetsRec->SetOwner(kFALSE);
1123 fListK0s = new TList();
1124 fListK0s->SetOwner(kFALSE);
1125 fListLa = new TList();
1126 fListLa->SetOwner(kFALSE);
1127 fListALa = new TList();
1128 fListALa->SetOwner(kFALSE);
1129 fListFeeddownLaCand = new TList(); //feeddown Lambda candidates
1130 fListFeeddownLaCand->SetOwner(kFALSE);
1131 fListFeeddownALaCand = new TList(); //feeddown Antilambda candidates
1132 fListFeeddownALaCand->SetOwner(kFALSE);
1133 jetConeFDLalist = new TList();
1134 jetConeFDLalist->SetOwner(kFALSE); //feeddown Lambda candidates in jet cone
1135 jetConeFDALalist = new TList();
1136 jetConeFDALalist->SetOwner(kFALSE); //feeddown Antilambda candidates in jet cone
1137 fListMCgenK0s = new TList(); //MC generated K0s
1138 fListMCgenK0s->SetOwner(kFALSE);
1139 fListMCgenLa = new TList(); //MC generated Lambdas
1140 fListMCgenLa->SetOwner(kFALSE);
1141 fListMCgenALa = new TList(); //MC generated Antilambdas
1142 fListMCgenALa->SetOwner(kFALSE);
1145 // Create histograms / output container
1147 fCommonHistList = new TList();
1148 fCommonHistList->SetOwner();
1150 Bool_t oldStatus = TH1::AddDirectoryStatus();
1151 TH1::AddDirectory(kFALSE);//By default (fAddDirectory = kTRUE), histograms are automatically added to the list of objects in memory
1153 // histograms inherited from AliAnalysisTaskFragmentationFunction
1155 fh1EvtSelection = new TH1F("fh1EvtSelection", "Event Selection", 6, -0.5, 5.5);
1156 fh1EvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
1157 fh1EvtSelection->GetXaxis()->SetBinLabel(2,"event trigger selection: rejected");
1158 fh1EvtSelection->GetXaxis()->SetBinLabel(3,"event class: rejected");
1159 fh1EvtSelection->GetXaxis()->SetBinLabel(4,"vertex Ncontr: rejected");
1160 fh1EvtSelection->GetXaxis()->SetBinLabel(5,"vertex z: rejected");
1161 fh1EvtSelection->GetXaxis()->SetBinLabel(6,"vertex type: rejected");
1164 fh1EvtCent = new TH1F("fh1EvtCent","centrality",100,0.,100.);
1165 fh1VertexNContributors = new TH1F("fh1VertexNContributors", "Vertex N contributors", 11,-.5, 10.5);
1166 fh1VertexZ = new TH1F("fh1VertexZ", "Vertex z distribution", 30, -15., 15.);
1167 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
1168 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
1169 fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
1170 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
1171 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
1172 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
1173 fh1nRecJetsCuts = new TH1F("fh1nRecJetsCuts","reconstructed jets per event",100,-0.5,99.5);
1175 // histograms JetChem task
1177 fh1EvtAllCent = new TH1F("fh1EvtAllCent","before centrality selection",100,0.,100.);
1178 fh1Evt = new TH1F("fh1Evt", "All events runned over", 3, 0.,1.);
1179 fh1EvtMult = new TH1F("fh1EvtMult","multiplicity",240,0.,240.);
1180 fh1K0Mult = new TH1F("fh1K0Mult","K0 multiplicity",100,0.,100.);//500. all
1181 fh1dPhiJetK0 = new TH1F("fh1dPhiJetK0","",64,-1,5.4);
1182 fh1LaMult = new TH1F("fh1LaMult","La multiplicity",100,0.,100.);
1183 fh1dPhiJetLa = new TH1F("fh1dPhiJetLa","",64,-1,5.4);
1184 fh1ALaMult = new TH1F("fh1ALaMult","ALa multiplicity",100,0.,100.);
1185 fh1dPhiJetALa = new TH1F("fh1dPhiJetALa","",64,-1,5.4);
1186 fh1JetEta = new TH1F("fh1JetEta","#eta distribution of all jets",40,-2.,2.);
1187 fh1JetPhi = new TH1F("fh1JetPhi","#phi distribution of all jets",63,0.,6.3);
1188 fh2JetEtaPhi = new TH2F("fh2JetEtaPhi","#eta and #phi distribution of all jets",400,-2.,2.,63,0.,6.3);
1191 //fh1V0JetPt = new TH1F("fh1V0JetPt","p_{T} distribution of all jets containing v0s",200,0.,200.);
1192 fh1IMK0Cone = new TH1F("fh1IMK0Cone","p_{T} distribution of all jets containing K0s candidates",19,5.,100.);
1193 fh1IMLaCone = new TH1F("fh1IMLaCone","p_{T} distribution of all jets containing #Lambda candidates",19,5.,100.);
1194 fh1IMALaCone = new TH1F("fh1IMALaCone","p_{T} distribution of all jets containing #bar{#Lambda} candidates",19,5.,100.);
1196 fh2FFJetTrackEta = new TH2F("fh2FFJetTrackEta","charged track eta distr. in jet cone",200,-1.,1.,40,0.,200.);
1197 //fh1trackPosNCls = new TH1F("fh1trackPosNCls","NTPC clusters positive daughters",10,0.,100.);
1198 //fh1trackNegNCls = new TH1F("fh1trackNegNCls","NTPC clusters negative daughters",10,0.,100.);
1199 fh1trackPosEta = new TH1F("fh1trackPosEta","eta positive daughters",100,-2.,2.);
1200 fh1trackNegEta = new TH1F("fh1trackNegEta","eta negative daughters",100,-2.,2.);
1201 fh1V0Eta = new TH1F("fh1V0Eta","V0 eta",60,-1.5,1.5);
1202 //fh1V0totMom = new TH1F("fh1V0totMom","V0 tot mom",100,0.,20.);
1203 fh1CosPointAngle = new TH1F("fh1CosPointAngle", "Cosine of V0's pointing angle",50,0.99,1.0);
1204 fh1DecayLengthV0 = new TH1F("fh1DecayLengthV0", "V0s decay Length;decay length(cm)",1200,0.,120.);
1205 fh2ProperLifetimeK0sVsPtBeforeCut = new TH2F("fh2ProperLifetimeK0sVsPtBeforeCut"," K0s ProperLifetime vs Pt; p_{T} (GeV/#it{c})",150,0.,15.,250,0.,250.);
1206 fh2ProperLifetimeK0sVsPtAfterCut = new TH2F("fh2ProperLifetimeK0sVsPtAfterCut"," K0s ProperLifetime vs Pt; p_{T} (GeV/#it{c})",150,0.,15.,250,0.,250.);
1207 fh1V0Radius = new TH1F("fh1V0Radius", "V0s Radius;Radius(cm)",200,0.,40.);
1208 fh1DcaV0Daughters = new TH1F("fh1DcaV0Daughters", "DCA between daughters;dca(cm)",200,0.,2.);
1209 fh1DcaPosToPrimVertex = new TH1F("fh1DcaPosToPrimVertex", "Positive V0 daughter;dca(cm)",100,0.,10.);
1210 fh1DcaNegToPrimVertex = new TH1F("fh1DcaNegToPrimVertex", "Negative V0 daughter;dca(cm)",100,0.,10.);
1211 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);
1212 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);
1213 fh2BBLaPos = new TH2F("fh2BBLaPos","PID of the positive daughter of La candidates; P (GeV); -dE/dx (keV/cm ?)",100,0,10,200,0,200);
1214 fh2BBLaNeg = new TH2F("fh2BBLaNeg","PID of the negative daughter of La candidates; P (GeV); -dE/dx (keV/cm ?)",100,0,10,200,0,200);
1215 fh1PosDaughterCharge = new TH1F("fh1PosDaughterCharge","charge of V0 positive daughters; V0 daughters",3,-2.,2.);
1216 fh1NegDaughterCharge = new TH1F("fh1NegDaughterCharge","charge of V0 negative daughters; V0 daughters",3,-2.,2.);
1217 fh1PtMCK0s = new TH1F("fh1PtMCK0s","Pt of MC rec K0s; #it{p}_{T} (GeV/#it{c})",200,0.,20.);
1218 fh1PtMCLa = new TH1F("fh1PtMCLa","Pt of MC rec La; #it{p}_{T} (GeV/#it{c})",200,0.,20.);
1219 fh1PtMCALa = new TH1F("fh1PtMCALa","Pt of MC rec ALa; #it{p}_{T} (GeV/#it{c})",200,0.,20.);
1220 fh1EtaK0s = new TH1F("fh1EtaK0s","K^{0}_{s} entries ;#eta",200,-1.,1.);
1221 fh1EtaLa = new TH1F("fh1EtaLa","#Lambda entries ;#eta",200,-1.,1.);
1222 fh1EtaALa = new TH1F("fh1EtaALa","#bar{#Lambda} entries ;#eta",200,-1.,1.);
1224 Int_t binsInvMassEtaTrackPtK0s[3] = {200, 200, 120};//eta,invM,trackPt
1225 Double_t xminInvMassEtaTrackPtK0s[3] = {-1.,0.3,0.};
1226 Double_t xmaxInvMassEtaTrackPtK0s[3] = {1.,0.7,12.};
1228 fhnInvMassEtaTrackPtK0s = new THnSparseF("fhnInvMassEtaTrackPtK0s","#eta; K0s invM (GeV/{#it{c}}^{2}); #it{p}_{T} (GeV/#it{c})",3,binsInvMassEtaTrackPtK0s,xminInvMassEtaTrackPtK0s,xmaxInvMassEtaTrackPtK0s);
1230 Int_t binsInvMassEtaTrackPtLa[3] = {200, 200, 120};//eta,invM,trackPt
1231 Double_t xminInvMassEtaTrackPtLa[3] = {-1.,1.05,0.};
1232 Double_t xmaxInvMassEtaTrackPtLa[3] = {1.,1.25,12.};
1234 fhnInvMassEtaTrackPtLa = new THnSparseF("fhnInvMassEtaTrackPtLa","#eta; #Lambda invM (GeV/{#it{c}}^{2}); #it{p}_{T} (GeV/#it{c})",3,binsInvMassEtaTrackPtLa,xminInvMassEtaTrackPtLa,xmaxInvMassEtaTrackPtLa);
1236 Int_t binsInvMassEtaTrackPtALa[3] = {200, 200, 120};//eta,invM,trackPt
1237 Double_t xminInvMassEtaTrackPtALa[3] = {-1.,1.05,0.};
1238 Double_t xmaxInvMassEtaTrackPtALa[3] = {1.,1.25,12.};
1240 fhnInvMassEtaTrackPtALa = new THnSparseF("fhnInvMassEtaTrackPtALa","#eta; #bar{#Lambda} invM (GeV/{#it{c}}^{2}); #it{p}_{T} (GeV/#it{c})",3,binsInvMassEtaTrackPtALa,xminInvMassEtaTrackPtALa,xmaxInvMassEtaTrackPtALa);
1242 Int_t binsK0sPC[4] = {19, 200, 200, 200};
1243 Double_t xminK0sPC[4] = {5.,0.3, 0., -1.};
1244 Double_t xmaxK0sPC[4] = {100.,0.7, 20., 1.};
1245 fhnK0sPC = new THnSparseF("fhnK0sPC","jet pT; K0s invM; particle pT; particle #eta",4,binsK0sPC,xminK0sPC,xmaxK0sPC);
1247 Int_t binsLaPC[4] = {19, 200, 200, 200};
1248 Double_t xminLaPC[4] = {5.,1.05, 0., -1.};
1249 Double_t xmaxLaPC[4] = {100.,1.25, 20., 1.};
1250 fhnLaPC = new THnSparseF("fhnLaPC","jet pT; #Lambda invM; particle pT; particle #eta",4,binsLaPC,xminLaPC,xmaxLaPC);
1252 Int_t binsALaPC[4] = {19, 200, 200, 200};
1253 Double_t xminALaPC[4] = {5.,1.05, 0., -1.};
1254 Double_t xmaxALaPC[4] = {100.,1.25, 20., 1.};
1255 fhnALaPC = new THnSparseF("fhnALaPC","jet pT; #bar#Lambda invM; particle pT; particle #eta",4,binsALaPC,xminALaPC,xmaxALaPC);
1257 Int_t binsK0sMCC[3] = {200, 200, 200};
1258 Double_t xminK0sMCC[3] = {0.3, 0., -1.};
1259 Double_t xmaxK0sMCC[3] = {0.7, 20., 1.};
1260 fhnK0sMCC = new THnSparseF("fhnK0sMCC","jet pT; K0s invM; particle pT; particle #eta",3,binsK0sMCC,xminK0sMCC,xmaxK0sMCC);
1262 Int_t binsLaMCC[3] = {200, 200, 200};
1263 Double_t xminLaMCC[3] = {1.05, 0., -1.};
1264 Double_t xmaxLaMCC[3] = {1.25, 20., 1.};
1265 fhnLaMCC = new THnSparseF("fhnLaMCC","jet pT; #Lambda invM; particle pT; particle #eta",3,binsLaMCC,xminLaMCC,xmaxLaMCC);
1267 Int_t binsALaMCC[3] = {200, 200, 200};
1268 Double_t xminALaMCC[3] = {1.05, 0., -1.};
1269 Double_t xmaxALaMCC[3] = {1.25, 20., 1.};
1270 fhnALaMCC = new THnSparseF("fhnALaMCC","jet pT; #bara#Lambda invM; particle pT; particle #eta",3,binsALaMCC,xminALaMCC,xmaxALaMCC);
1272 Int_t binsK0sRC[3] = {200, 200, 200};
1273 Double_t xminK0sRC[3] = {0.3, 0., -1.};
1274 Double_t xmaxK0sRC[3] = {0.7, 20., 1.};
1275 fhnK0sRC = new THnSparseF("fhnK0sRC","jet pT; K0s invM; particle pT; particle #eta",3,binsK0sRC,xminK0sRC,xmaxK0sRC);
1277 Int_t binsLaRC[3] = {200, 200, 200};
1278 Double_t xminLaRC[3] = {1.05, 0., -1.};
1279 Double_t xmaxLaRC[3] = {1.25, 20., 1.};
1280 fhnLaRC = new THnSparseF("fhnLaRC","jet pT; #Lambda invM; particle pT; particle #eta",3,binsLaRC,xminLaRC,xmaxLaRC);
1282 Int_t binsALaRC[3] = {200, 200, 200};
1283 Double_t xminALaRC[3] = {1.05, 0., -1.};
1284 Double_t xmaxALaRC[3] = {1.25, 20., 1.};
1285 fhnALaRC = new THnSparseF("fhnALaRC","jet pT; #bara#Lambda invM; particle pT; particle #eta",3,binsALaRC,xminALaRC,xmaxALaRC);
1288 Int_t binsK0sOC[3] = {200, 200, 200};
1289 Double_t xminK0sOC[3] = {0.3, 0., -1.};
1290 Double_t xmaxK0sOC[3] = {0.7, 20., 1.};
1291 fhnK0sOC = new THnSparseF("fhnK0sOC","jet pT; K0s invM; particle pT; particle #eta",3,binsK0sOC,xminK0sOC,xmaxK0sOC);
1293 Int_t binsLaOC[3] = {200, 200, 200};
1294 Double_t xminLaOC[3] = {1.05, 0., -1.};
1295 Double_t xmaxLaOC[3] = {1.25, 20., 1.};
1296 fhnLaOC = new THnSparseF("fhnLaOC","jet pT; #Lambda invM; particle pT; particle #eta",3,binsLaOC,xminLaOC,xmaxLaOC);
1298 Int_t binsALaOC[3] = {200, 200, 200};
1299 Double_t xminALaOC[3] = {1.05, 0., -1.};
1300 Double_t xmaxALaOC[3] = {1.25, 20., 1.};
1302 fhnALaOC = new THnSparseF("fhnALaOC","jet pT; #bara#Lambda invM; particle pT; particle #eta",3,binsALaOC,xminALaOC,xmaxALaOC);
1308 fh1AreaExcluded = new TH1F("fh1AreaExcluded","area excluded for selected jets in event acceptance",100,0.,5.);
1310 fh1MedianEta = new TH1F("fh1MedianEta","Median cluster axis ;#eta",200,-1.,1.);
1311 fh1JetPtMedian = new TH1F("fh1JetPtMedian"," (selected) jet it{p}_{T} distribution for MCC method; #GeV/it{c}",19,5.,100.);
1313 fh1TrackMultCone = new TH1F("fh1TrackMultCone","track multiplicity in jet cone; number of tracks",20,0.,50.);
1315 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.);
1317 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.);
1319 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.);
1321 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.);
1323 fFFHistosRecCuts = new AliFragFuncHistos("RecCuts", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1324 fFFNBinsPt, fFFPtMin, fFFPtMax,
1325 fFFNBinsXi, fFFXiMin, fFFXiMax,
1326 fFFNBinsZ , fFFZMin , fFFZMax);
1328 fV0QAK0 = new AliFragFuncQATrackHistos("V0QAK0",fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
1329 fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1330 fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
1331 fQATrackHighPtThreshold);
1333 fFFHistosRecCutsK0Evt = new AliFragFuncHistos("RecCutsK0Evt", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1334 fFFNBinsPt, fFFPtMin, fFFPtMax,
1335 fFFNBinsXi, fFFXiMin, fFFXiMax,
1336 fFFNBinsZ , fFFZMin , fFFZMax);
1339 fFFHistosIMK0AllEvt = new AliFragFuncHistosInvMass("K0AllEvt", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax,
1340 fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
1341 fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax,
1342 fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,
1343 fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
1345 fFFHistosIMK0Jet = new AliFragFuncHistosInvMass("K0Jet", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax,
1346 fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
1347 fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax,
1348 fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,
1349 fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
1351 fFFHistosIMK0Cone = new AliFragFuncHistosInvMass("K0Cone", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax,
1352 fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
1353 fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax,
1354 fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,
1355 fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
1357 fFFHistosIMLaAllEvt = new AliFragFuncHistosInvMass("LaAllEvt", fFFIMLaNBinsJetPt, fFFIMLaJetPtMin, fFFIMLaJetPtMax,
1358 fFFIMLaNBinsInvM,fFFIMLaInvMMin,fFFIMLaInvMMax,
1359 fFFIMLaNBinsPt, fFFIMLaPtMin, fFFIMLaPtMax,
1360 fFFIMLaNBinsXi, fFFIMLaXiMin, fFFIMLaXiMax,
1361 fFFIMLaNBinsZ , fFFIMLaZMin , fFFIMLaZMax);
1363 fFFHistosIMLaJet = new AliFragFuncHistosInvMass("LaJet", fFFIMLaNBinsJetPt, fFFIMLaJetPtMin, fFFIMLaJetPtMax,
1364 fFFIMLaNBinsInvM,fFFIMLaInvMMin,fFFIMLaInvMMax,
1365 fFFIMLaNBinsPt, fFFIMLaPtMin, fFFIMLaPtMax,
1366 fFFIMLaNBinsXi, fFFIMLaXiMin, fFFIMLaXiMax,
1367 fFFIMLaNBinsZ , fFFIMLaZMin , fFFIMLaZMax);
1370 fFFHistosIMLaCone = new AliFragFuncHistosInvMass("LaCone", fFFIMLaNBinsJetPt, fFFIMLaJetPtMin, fFFIMLaJetPtMax,
1371 fFFIMLaNBinsInvM,fFFIMLaInvMMin,fFFIMLaInvMMax,
1372 fFFIMLaNBinsPt, fFFIMLaPtMin, fFFIMLaPtMax,
1373 fFFIMLaNBinsXi, fFFIMLaXiMin, fFFIMLaXiMax,
1374 fFFIMLaNBinsZ , fFFIMLaZMin , fFFIMLaZMax);
1377 fFFHistosIMALaAllEvt = new AliFragFuncHistosInvMass("ALaAllEvt", fFFIMLaNBinsJetPt, fFFIMLaJetPtMin, fFFIMLaJetPtMax,
1378 fFFIMLaNBinsInvM,fFFIMLaInvMMin,fFFIMLaInvMMax,
1379 fFFIMLaNBinsPt, fFFIMLaPtMin, fFFIMLaPtMax,
1380 fFFIMLaNBinsXi, fFFIMLaXiMin, fFFIMLaXiMax,
1381 fFFIMLaNBinsZ , fFFIMLaZMin , fFFIMLaZMax);
1383 fFFHistosIMALaJet = new AliFragFuncHistosInvMass("ALaJet", fFFIMLaNBinsJetPt, fFFIMLaJetPtMin, fFFIMLaJetPtMax,
1384 fFFIMLaNBinsInvM,fFFIMLaInvMMin,fFFIMLaInvMMax,
1385 fFFIMLaNBinsPt, fFFIMLaPtMin, fFFIMLaPtMax,
1386 fFFIMLaNBinsXi, fFFIMLaXiMin, fFFIMLaXiMax,
1387 fFFIMLaNBinsZ , fFFIMLaZMin , fFFIMLaZMax);
1389 fFFHistosIMALaCone = new AliFragFuncHistosInvMass("ALaCone", fFFIMLaNBinsJetPt, fFFIMLaJetPtMin, fFFIMLaJetPtMax,
1390 fFFIMLaNBinsInvM,fFFIMLaInvMMin,fFFIMLaInvMMax,
1391 fFFIMLaNBinsPt, fFFIMLaPtMin, fFFIMLaPtMax,
1392 fFFIMLaNBinsXi, fFFIMLaXiMin, fFFIMLaXiMax,
1393 fFFIMLaNBinsZ , fFFIMLaZMin , fFFIMLaZMax);
1400 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.);
1401 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.);
1402 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.);
1404 fh2MCgenK0Cone->GetYaxis()->SetTitle("MC gen K^{0}}^{s} #it{p}_{T}");
1405 fh2MCgenLaCone->GetYaxis()->SetTitle("MC gen #Lambda #it{p}_{T}");
1406 fh2MCgenALaCone->GetYaxis()->SetTitle("MC gen #Antilambda #it{p}_{T}");
1408 fh2MCEtagenK0Cone = new TH2F("fh2MCEtagenK0Cone","MC gen {K^{0}}^{s} #it{p}_{T} #eta distribution in jet cone;#eta",19,5.,100.,200,-1.,1.);
1409 fh2MCEtagenLaCone = new TH2F("fh2MCEtagenLaCone","MC gen #Lambda #it{p}_{T} #eta distribution in jet cone;#eta",19,5.,100.,200,-1.,1.);
1410 fh2MCEtagenALaCone = new TH2F("fh2MCEtagenALaCone","MC gen #Antilambda #it{p}_{T} #eta distribution in jet cone;#eta",19,5.,100.,200,-1.,1.);
1411 fh1IMK0ConeSmear = new TH1F("fh1IMK0ConeSmear","Smeared jet pt study for K0s-in-cone-jets; smeared jet #it{p}_{T}", 19,5.,100.);
1412 fh1IMLaConeSmear = new TH1F("fh1IMLaConeSmear","Smeared jet pt study for La-in-cone-jets; smeared jet #it{p}_{T}", 19,5.,100.);
1413 fh1IMALaConeSmear = new TH1F("fh1IMALaConeSmear","Smeared jet pt study for ALa-in-cone-jets; smeared jet #it{p}_{T}", 19,5.,100.);
1415 Int_t binsMCrecK0Cone[4] = {19, 200, 200, 200};
1416 Double_t xminMCrecK0Cone[4] = {5.,0.3, 0., -1.};
1417 Double_t xmaxMCrecK0Cone[4] = {100.,0.7, 20., 1.};
1418 fhnMCrecK0Cone = new THnSparseF("fhnMCrecK0Cone", "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}",4,binsMCrecK0Cone,xminMCrecK0Cone,xmaxMCrecK0Cone);
1420 Int_t binsMCrecLaCone[4] = {19, 200, 200, 200};
1421 Double_t xminMCrecLaCone[4] = {5.,0.3, 0., -1.};
1422 Double_t xmaxMCrecLaCone[4] = {100.,0.7, 20., 1.};
1423 fhnMCrecLaCone = new THnSparseF("fhnMCrecLaCone", "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}",4,binsMCrecLaCone,xminMCrecLaCone,xmaxMCrecLaCone);
1425 Int_t binsMCrecALaCone[4] = {19, 200, 200, 200};
1426 Double_t xminMCrecALaCone[4] = {5.,0.3, 0., -1.};
1427 Double_t xmaxMCrecALaCone[4] = {100.,0.7, 20., 1.};
1428 fhnMCrecALaCone = new THnSparseF("fhnMCrecALaCone", "MC rec {#bar{#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}",4,binsMCrecALaCone,xminMCrecALaCone,xmaxMCrecALaCone);
1430 Int_t binsMCrecK0ConeSmear[4] = {19, 200, 200, 200};
1431 Double_t xminMCrecK0ConeSmear[4] = {5.,0.3, 0., -1.};
1432 Double_t xmaxMCrecK0ConeSmear[4] = {100.,0.7, 20., 1.};
1433 fhnMCrecK0ConeSmear = new THnSparseF("fhnMCrecK0ConeSmear", "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}",4,binsMCrecK0ConeSmear,xminMCrecK0ConeSmear,xmaxMCrecK0ConeSmear);
1435 Int_t binsMCrecLaConeSmear[4] = {19, 200, 200, 200};
1436 Double_t xminMCrecLaConeSmear[4] = {5.,0.3, 0., -1.};
1437 Double_t xmaxMCrecLaConeSmear[4] = {100.,0.7, 20., 1.};
1438 fhnMCrecLaConeSmear = new THnSparseF("fhnMCrecLaConeSmear", "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}",4,binsMCrecLaConeSmear,xminMCrecLaConeSmear,xmaxMCrecLaConeSmear);
1440 Int_t binsMCrecALaConeSmear[4] = {19, 200, 200, 200};
1441 Double_t xminMCrecALaConeSmear[4] = {5.,0.3, 0., -1.};
1442 Double_t xmaxMCrecALaConeSmear[4] = {100.,0.7, 20., 1.};
1443 fhnMCrecALaConeSmear = new THnSparseF("fhnMCrecALaConeSmear", "MC rec {#bar{#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}",4,binsMCrecALaConeSmear,xminMCrecALaConeSmear,xmaxMCrecALaConeSmear);
1445 Int_t binsK0sSecContinCone[3] = {19, 200, 200};
1446 Double_t xminK0sSecContinCone[3] = {5.,0., -1.};
1447 Double_t xmaxK0sSecContinCone[3] = {100.,20., 1.};
1448 fhnK0sSecContinCone = new THnSparseF("fhnK0sSecContinCone", "Secondary contamination {K^{0}}^{s} #it{p}_{T} in cone around jet axis; jet #it{p}_{T}; inv mass (GeV/#it{c}^{2};#it{p}_{T}",3,binsK0sSecContinCone,xminK0sSecContinCone,xmaxK0sSecContinCone);
1450 Int_t binsLaSecContinCone[3] = {19, 200, 200};
1451 Double_t xminLaSecContinCone[3] = {5.,0., -1.};
1452 Double_t xmaxLaSecContinCone[3] = {100.,20., 1.};
1453 fhnLaSecContinCone = new THnSparseF("fhnLaSecContinCone", "Secondary contamination {#Lambda #it{p}_{T} in cone around jet axis; jet #it{p}_{T}; inv mass (GeV/#it{c}^{2};#it{p}_{T}",3,binsLaSecContinCone,xminLaSecContinCone,xmaxLaSecContinCone);
1455 Int_t binsALaSecContinCone[3] = {19, 200, 200};
1456 Double_t xminALaSecContinCone[3] = {5.,0., -1.};
1457 Double_t xmaxALaSecContinCone[3] = {100.,20., 1.};
1458 fhnALaSecContinCone = new THnSparseF("fhnALaSecContinCone", "Secondary contamination {#bar{#Lambda} #it{p}_{T} in cone around jet axis; jet #it{p}_{T}; inv mass (GeV/#it{c}^{2};#it{p}_{T}",3,binsALaSecContinCone,xminALaSecContinCone,xmaxALaSecContinCone);
1460 Int_t binsK0sIncl[3] = {200, 200, 200};
1461 Double_t xminK0sIncl[3] = {0.3, 0., -1.};
1462 Double_t xmaxK0sIncl[3] = {0.7, 20., 1.};
1463 fhnK0sIncl = new THnSparseF("fhnK0sIncl","K0s inv. mass; particle pT; particle #eta",3,binsK0sIncl,xminK0sIncl,xmaxK0sIncl);
1465 Int_t binsK0sCone[4] = {19, 200, 200, 200};
1466 Double_t xminK0sCone[4] = {5.,0.3, 0., -1.};
1467 Double_t xmaxK0sCone[4] = {100.,0.7, 20., 1.};
1468 fhnK0sCone = new THnSparseF("fhnK0sCone","jet pT; K0s inv. mass; particle pT; particle #eta",4,binsK0sCone,xminK0sCone,xmaxK0sCone);
1470 Int_t binsLaIncl[3] = {200, 200, 200};
1471 Double_t xminLaIncl[3] = {1.05, 0., -1.};
1472 Double_t xmaxLaIncl[3] = {1.25, 20., 1.};
1473 fhnLaIncl = new THnSparseF("fhnLaIncl","La inv. mass; particle pT; particle #eta",3,binsLaIncl,xminLaIncl,xmaxLaIncl);
1475 Int_t binsLaCone[4] = {19, 200, 200, 200};
1476 Double_t xminLaCone[4] = {5.,1.05, 0., -1.};
1477 Double_t xmaxLaCone[4] = {100.,1.25, 20., 1.};
1478 fhnLaCone = new THnSparseF("fhnLaCone","jet pT; La inv. mass; particle pT; particle #eta",4,binsLaCone,xminLaCone,xmaxLaCone);
1480 Int_t binsALaIncl[3] = {200, 200, 200};
1481 Double_t xminALaIncl[3] = {1.05, 0., -1.};
1482 Double_t xmaxALaIncl[3] = {1.25, 20., 1.};
1483 fhnALaIncl = new THnSparseF("fhnALaIncl","ALa inv. mass; particle pT; particle #eta",3,binsALaIncl,xminALaIncl,xmaxALaIncl);
1485 Int_t binsALaCone[4] = {19, 200, 200, 200};
1486 Double_t xminALaCone[4] = {5.,1.05, 0., -1.};
1487 Double_t xmaxALaCone[4] = {100.,1.25, 20., 1.};
1488 fhnALaCone = new THnSparseF("fhnALaCone","jet pT; ALa inv. mass; particle pT; particle #eta",4,binsALaCone,xminALaCone,xmaxALaCone);
1490 fh1MCMultiplicityPrimary = new TH1F("fh1MCMultiplicityPrimary", "MC Primary Particles;NPrimary;Count", 201, -0.5, 200.5);
1491 fh1MCMultiplicityTracks = new TH1F("h1MCMultiplicityTracks", "MC Tracks;Ntracks;Count", 201, -0.5, 200.5);
1494 Int_t binsFeedDownLa[3] = {19, 200, 200};
1495 Double_t xminFeedDownLa[3] = {5.,1.05, 0.};
1496 Double_t xmaxFeedDownLa[3] = {100.,1.25, 20.};
1497 fhnFeedDownLa = new THnSparseF("fhnFeedDownLa","#Lambda stemming from feeddown from Xi(0/-)",3,binsFeedDownLa,xminFeedDownLa,xmaxFeedDownLa);
1499 Int_t binsFeedDownALa[3] = {19, 200, 200};
1500 Double_t xminFeedDownALa[3] = {5.,1.05, 0.};
1501 Double_t xmaxFeedDownALa[3] = {100.,1.25, 20.};
1502 fhnFeedDownALa = new THnSparseF("fhnFeedDownALa","#bar#Lambda stemming from feeddown from Xibar(0/+)",3,binsFeedDownALa,xminFeedDownALa,xmaxFeedDownALa);
1504 Int_t binsFeedDownLaCone[3] = {19, 200, 200};
1505 Double_t xminFeedDownLaCone[3] = {5.,1.05, 0.};
1506 Double_t xmaxFeedDownLaCone[3] = {100.,1.25, 20.};
1507 fhnFeedDownLaCone = new THnSparseF("fhnFeedDownLaCone","#Lambda stemming from feeddown from Xi(0/-) in jet cone",3,binsFeedDownLaCone,xminFeedDownLaCone,xmaxFeedDownLaCone);
1509 Int_t binsFeedDownALaCone[3] = {19, 200, 200};
1510 Double_t xminFeedDownALaCone[3] = {5.,1.05, 0.};
1511 Double_t xmaxFeedDownALaCone[3] = {100.,1.25, 20.};
1512 fhnFeedDownALaCone = new THnSparseF("fhnFeedDownALaCone","#bar#Lambda stemming from feeddown from Xibar(0/+) in jet cone",3,binsFeedDownALaCone,xminFeedDownALaCone,xmaxFeedDownALaCone);
1515 fh1MCProdRadiusK0s = new TH1F("fh1MCProdRadiusK0s","MC gen. MC K0s prod radius",200,0.,200.);
1516 fh1MCProdRadiusLambda = new TH1F("fh1MCProdRadiusLambda","MC gen. MC La prod radius",200,0.,200.);
1517 fh1MCProdRadiusAntiLambda = new TH1F("fh1MCProdRadiusAntiLambda","MC gen. MC ALa prod radius",200,0.,200.);
1519 // Pt and inv mass distributions
1521 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
1522 fh1MCPtK0s = new TH1F("fh1MCPtK0s", "MC gen. K^{0}_{s} in eta range;#it{p}_{T} (GeV/#it{c})",200,0.,20.);
1523 fh1MCPtLambda = new TH1F("fh1MCPtLambda", "MC gen. #Lambda in rap range;#it{p}_{T} (GeV/#it{c})",200,0.,20.);
1524 fh1MCPtAntiLambda = new TH1F("fh1MCPtAntiLambda", "MC gen. #AntiLambda in rap range;#it{p}_{T} (GeV/#it{c})",200,0.,20.);
1525 fh1MCXiPt = new TH1F("fh1MCXiPt", "MC gen. #Xi^{-/o};#it{p}_{T} (GeV/#it{c})",200,0.,20.);
1526 fh1MCXibarPt = new TH1F("fh1MCXibarPt", "MC gen. #bar{#Xi}^{+/o};#it{p}_{T} (GeV/#it{c})",200,0.,20.);
1527 fh2MCEtaVsPtK0s = new TH2F("fh2MCEtaVsPtK0s","MC gen. K^{0}_{s} #eta; #it{p}_{T}",200,0.,20.,200,-1.,1.);
1528 fh2MCEtaVsPtLa = new TH2F("fh2MCEtaVsPtLa","MC gen. #Lambda #eta; #it{p}_{T}",200,0.,20.,200,-1.,1.);
1529 fh2MCEtaVsPtALa = new TH2F("fh2MCEtaVsPtALa","MC gen. #bar{#Lambda} #eta; #it{p}_{T}",200,0.,20.,200,-1.,1.);
1532 fh1MCRapK0s = new TH1F("fh1MCRapK0s", "MC gen. K0s;rap with cut",200,-10,10);
1533 fh1MCRapLambda = new TH1F("fh1MCRapLambda", "MC gen. #Lambda;rap",200,-10,10);
1534 fh1MCRapAntiLambda = new TH1F("fh1MCRapAntiLambda", "MC gen. #bar{#Lambda};rap",200,-10,10);
1535 fh1MCEtaAllK0s = new TH1F("fh1MCEtaAllK0s", "MC gen. K0s;#eta",200,-1.,1.);
1536 fh1MCEtaK0s = new TH1F("fh1MCEtaK0s", "MC gen. K0s;#eta with cut",200,-1.,1.);
1537 fh1MCEtaLambda = new TH1F("fh1MCEtaLambda", "MC gen. #Lambda;#eta",200,-1.,1.);
1538 fh1MCEtaAntiLambda = new TH1F("fh1MCEtaAntiLambda", "MC gen. #bar{#Lambda};#eta",200,-1.,1.);
1540 fV0QAK0->DefineHistos();
1541 fFFHistosRecCuts->DefineHistos();
1542 fFFHistosRecCutsK0Evt->DefineHistos();
1543 /* fFFHistosIMK0AllEvt->DefineHistos();
1544 fFFHistosIMK0Jet->DefineHistos();
1545 fFFHistosIMK0Cone->DefineHistos();
1546 fFFHistosIMLaAllEvt->DefineHistos();
1547 fFFHistosIMLaJet->DefineHistos();
1548 fFFHistosIMLaCone->DefineHistos();
1549 fFFHistosIMALaAllEvt->DefineHistos();
1550 fFFHistosIMALaJet->DefineHistos();
1551 fFFHistosIMALaCone->DefineHistos();
1554 const Int_t saveLevel = 5;
1557 fCommonHistList->Add(fh1EvtAllCent);
1558 fCommonHistList->Add(fh1Evt);
1559 fCommonHistList->Add(fh1EvtSelection);
1560 fCommonHistList->Add(fh1EvtCent);
1561 fCommonHistList->Add(fh1VertexNContributors);
1562 fCommonHistList->Add(fh1VertexZ);
1563 fCommonHistList->Add(fh1Xsec);
1564 fCommonHistList->Add(fh1Trials);
1565 fCommonHistList->Add(fh1PtHard);
1566 fCommonHistList->Add(fh1PtHardTrials);
1567 fCommonHistList->Add(fh1nRecJetsCuts);
1568 fCommonHistList->Add(fh1EvtMult);
1569 fCommonHistList->Add(fh1K0Mult);
1570 fCommonHistList->Add(fh1dPhiJetK0);
1571 fCommonHistList->Add(fh1LaMult);
1572 fCommonHistList->Add(fh1dPhiJetLa);
1573 fCommonHistList->Add(fh1ALaMult);
1574 fCommonHistList->Add(fh1dPhiJetALa);
1575 fCommonHistList->Add(fh1JetEta);
1576 fCommonHistList->Add(fh1JetPhi);
1577 fCommonHistList->Add(fh2JetEtaPhi);
1578 //fCommonHistList->Add(fh1V0JetPt);
1579 fCommonHistList->Add(fh1IMK0Cone);
1580 fCommonHistList->Add(fh1IMLaCone);
1581 fCommonHistList->Add(fh1IMALaCone);
1582 fCommonHistList->Add(fh2FFJetTrackEta);
1583 // fCommonHistList->Add(fh1trackPosNCls);
1584 //fCommonHistList->Add(fh1trackNegNCls);
1585 fCommonHistList->Add(fh1trackPosEta);
1586 fCommonHistList->Add(fh1trackNegEta);
1587 fCommonHistList->Add(fh1V0Eta);
1588 // fCommonHistList->Add(fh1V0totMom);
1589 fCommonHistList->Add(fh1CosPointAngle);
1590 fCommonHistList->Add(fh1DecayLengthV0);
1591 fCommonHistList->Add(fh2ProperLifetimeK0sVsPtBeforeCut);
1592 fCommonHistList->Add(fh2ProperLifetimeK0sVsPtAfterCut);
1593 fCommonHistList->Add(fh1V0Radius);
1594 fCommonHistList->Add(fh1DcaV0Daughters);
1595 fCommonHistList->Add(fh1DcaPosToPrimVertex);
1596 fCommonHistList->Add(fh1DcaNegToPrimVertex);
1597 fCommonHistList->Add(fh2ArmenterosBeforeCuts);
1598 fCommonHistList->Add(fh2ArmenterosAfterCuts);
1599 fCommonHistList->Add(fh2BBLaPos);
1600 fCommonHistList->Add(fh2BBLaNeg);
1601 fCommonHistList->Add(fh1PosDaughterCharge);
1602 fCommonHistList->Add(fh1NegDaughterCharge);
1603 fCommonHistList->Add(fh1PtMCK0s);
1604 fCommonHistList->Add(fh1PtMCLa);
1605 fCommonHistList->Add(fh1PtMCALa);
1606 fCommonHistList->Add(fh1EtaK0s);
1607 fCommonHistList->Add(fh1EtaLa);
1608 fCommonHistList->Add(fh1EtaALa);
1609 fCommonHistList->Add(fhnInvMassEtaTrackPtK0s);
1610 fCommonHistList->Add(fhnInvMassEtaTrackPtLa);
1611 fCommonHistList->Add(fhnInvMassEtaTrackPtALa);
1612 fCommonHistList->Add(fh1TrackMultCone);
1613 fCommonHistList->Add(fh2TrackMultCone);
1614 fCommonHistList->Add(fh2NJK0);
1615 fCommonHistList->Add(fh2NJLa);
1616 fCommonHistList->Add(fh2NJALa);
1617 fCommonHistList->Add(fh2MCgenK0Cone);
1618 fCommonHistList->Add(fh2MCgenLaCone);
1619 fCommonHistList->Add(fh2MCgenALaCone);
1620 fCommonHistList->Add(fh2MCEtagenK0Cone);
1621 fCommonHistList->Add(fh2MCEtagenLaCone);
1622 fCommonHistList->Add(fh2MCEtagenALaCone);
1623 fCommonHistList->Add(fh1IMK0ConeSmear);
1624 fCommonHistList->Add(fh1IMLaConeSmear);
1625 fCommonHistList->Add(fh1IMALaConeSmear);
1626 fCommonHistList->Add(fhnMCrecK0Cone);
1627 fCommonHistList->Add(fhnMCrecLaCone);
1628 fCommonHistList->Add(fhnMCrecALaCone);
1629 fCommonHistList->Add(fhnMCrecK0ConeSmear);
1630 fCommonHistList->Add(fhnMCrecLaConeSmear);
1631 fCommonHistList->Add(fhnMCrecALaConeSmear);
1632 fCommonHistList->Add(fhnK0sSecContinCone);
1633 fCommonHistList->Add(fhnLaSecContinCone);
1634 fCommonHistList->Add(fhnALaSecContinCone);
1635 fCommonHistList->Add(fhnK0sIncl);
1636 fCommonHistList->Add(fhnK0sCone);
1637 fCommonHistList->Add(fhnLaIncl);
1638 fCommonHistList->Add(fhnLaCone);
1639 fCommonHistList->Add(fhnALaIncl);
1640 fCommonHistList->Add(fhnALaCone);
1641 fCommonHistList->Add(fhnK0sPC);
1642 fCommonHistList->Add(fhnLaPC);
1643 fCommonHistList->Add(fhnALaPC);
1644 fCommonHistList->Add(fhnK0sMCC);
1645 fCommonHistList->Add(fhnLaMCC);
1646 fCommonHistList->Add(fhnALaMCC);
1647 fCommonHistList->Add(fhnK0sRC);
1648 fCommonHistList->Add(fhnLaRC);
1649 fCommonHistList->Add(fhnALaRC);
1650 fCommonHistList->Add(fhnK0sOC);
1651 fCommonHistList->Add(fhnLaOC);
1652 fCommonHistList->Add(fhnALaOC);
1653 fCommonHistList->Add(fh1AreaExcluded);
1654 fCommonHistList->Add(fh1MedianEta);
1655 fCommonHistList->Add(fh1JetPtMedian);
1656 fCommonHistList->Add(fh1MCMultiplicityPrimary);
1657 fCommonHistList->Add(fh1MCMultiplicityTracks);
1658 fCommonHistList->Add(fhnFeedDownLa);
1659 fCommonHistList->Add(fhnFeedDownALa);
1660 fCommonHistList->Add(fhnFeedDownLaCone);
1661 fCommonHistList->Add(fhnFeedDownALaCone);
1662 fCommonHistList->Add(fh1MCProdRadiusK0s);
1663 fCommonHistList->Add(fh1MCProdRadiusLambda);
1664 fCommonHistList->Add(fh1MCProdRadiusAntiLambda);
1665 fCommonHistList->Add(fh1MCPtV0s);
1666 fCommonHistList->Add(fh1MCPtK0s);
1667 fCommonHistList->Add(fh1MCPtLambda);
1668 fCommonHistList->Add(fh1MCPtAntiLambda);
1669 fCommonHistList->Add(fh1MCXiPt);
1670 fCommonHistList->Add(fh1MCXibarPt);
1671 fCommonHistList->Add(fh2MCEtaVsPtK0s);
1672 fCommonHistList->Add(fh2MCEtaVsPtLa);
1673 fCommonHistList->Add(fh2MCEtaVsPtALa);
1674 fCommonHistList->Add(fh1MCRapK0s);
1675 fCommonHistList->Add(fh1MCRapLambda);
1676 fCommonHistList->Add(fh1MCRapAntiLambda);
1677 fCommonHistList->Add(fh1MCEtaAllK0s);
1678 fCommonHistList->Add(fh1MCEtaK0s);
1679 fCommonHistList->Add(fh1MCEtaLambda);
1680 fCommonHistList->Add(fh1MCEtaAntiLambda);
1684 fV0QAK0->AddToOutput(fCommonHistList);
1685 fFFHistosRecCuts->AddToOutput(fCommonHistList);
1686 fFFHistosRecCutsK0Evt->AddToOutput(fCommonHistList);
1687 // fFFHistosIMK0AllEvt->AddToOutput(fCommonHistList);
1688 // fFFHistosIMK0Jet->AddToOutput(fCommonHistList);
1689 // fFFHistosIMK0Cone->AddToOutput(fCommonHistList);
1690 // fFFHistosIMLaAllEvt->AddToOutput(fCommonHistList);
1691 // fFFHistosIMLaJet->AddToOutput(fCommonHistList);
1692 // fFFHistosIMLaCone->AddToOutput(fCommonHistList);
1693 // fFFHistosIMALaAllEvt->AddToOutput(fCommonHistList);
1694 // fFFHistosIMALaJet->AddToOutput(fCommonHistList);
1695 // fFFHistosIMALaCone->AddToOutput(fCommonHistList);
1700 // =========== Switch on Sumw2 for all histos ===========
1701 for (Int_t i=0; i<fCommonHistList->GetEntries(); ++i){
1703 TH1 *h1 = dynamic_cast<TH1*>(fCommonHistList->At(i));
1705 if (h1) h1->Sumw2();//The error per bin will be computed as sqrt(sum of squares of weight) for each bin
1707 THnSparse *hnSparse = dynamic_cast<THnSparse*>(fCommonHistList->At(i));
1708 if(hnSparse) hnSparse->Sumw2();
1712 TH1::AddDirectory(oldStatus);
1713 PostData(1, fCommonHistList);
1716 //_______________________________________________
1717 void AliAnalysisTaskJetChem::UserExec(Option_t *)
1720 // Called for each event
1722 if(fDebug > 1) Printf("AliAnalysisTaskJetChem::UserExec()");
1724 if(fDebug > 1) Printf("Analysis event #%5d", (Int_t) fEntry);
1726 // Trigger selection
1727 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
1728 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1731 //for AliPIDResponse:
1732 //AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
1733 //AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1734 fPIDResponse = inputHandler->GetPIDResponse();
1736 if (!fPIDResponse){if(fDebug > 1) Printf("AliAnalysisTaskJetChem::UserExec(): fPIDResponse does not exist!"); return;}
1738 //std::cout<<"inputHandler->IsEventSelected(): "<<inputHandler->IsEventSelected()<<std::endl;
1739 //std::cout<<"fEvtSelectionMask: "<<fEvtSelectionMask<<std::endl;
1741 if(!(inputHandler->IsEventSelected() & fEvtSelectionMask)){
1742 //std::cout<<"########event rejected!!############"<<std::endl;
1743 fh1EvtSelection->Fill(1.);
1744 if (fDebug > 1 ) Printf(" Trigger Selection: event REJECTED ... ");
1745 PostData(1, fCommonHistList);
1749 fESD = dynamic_cast<AliESDEvent*>(InputEvent());//casting of pointers for inherited class, only for ESDs
1751 if(fDebug>3) Printf("%s:%d ESDEvent not found in the input", (char*)__FILE__,__LINE__);
1754 fMCEvent = MCEvent();
1756 if(fDebug>3) Printf("%s:%d MCEvent not found in the input", (char*)__FILE__,__LINE__);
1759 // get AOD event from input/output
1760 TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
1761 if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
1762 fAOD = ((AliAODInputHandler*)handler)->GetEvent();
1763 if(fUseAODInputJets) fAODJets = fAOD;
1764 if (fDebug > 1) Printf("%s:%d AOD event from input", (char*)__FILE__,__LINE__);
1767 handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
1768 if( handler && handler->InheritsFrom("AliAODHandler") ) {
1769 fAOD = ((AliAODHandler*)handler)->GetAOD();
1771 if (fDebug > 1) Printf("%s:%d AOD event from output", (char*)__FILE__,__LINE__);
1775 if(!fAODJets && !fUseAODInputJets){ // case we have AOD in input & output and want jets from output
1776 TObject* outHandler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
1777 if( outHandler && outHandler->InheritsFrom("AliAODHandler") ){
1778 fAODJets = ((AliAODHandler*)outHandler)->GetAOD();
1779 if (fDebug > 1) Printf("%s:%d jets from output AOD", (char*)__FILE__,__LINE__);
1783 if(fNonStdFile.Length()!=0){
1784 // case we have an AOD extension - fetch the jets from the extended output
1786 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1787 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
1789 if(fDebug>1)Printf("AODExtension not found for %s",fNonStdFile.Data());
1794 Printf("%s:%d AODEvent not found", (char*)__FILE__,__LINE__);
1798 Printf("%s:%d AODEvent with jet branch not found", (char*)__FILE__,__LINE__);
1802 //primary vertex position:
1803 AliAODVertex *myPrimaryVertex = NULL;
1804 myPrimaryVertex = (AliAODVertex*)fAOD->GetPrimaryVertex();
1805 if (!myPrimaryVertex) return;
1806 fh1Evt->Fill(1.);//fill in every event that was accessed with InputHandler
1808 // event selection *****************************************
1810 // *** vertex cut ***
1811 AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
1812 Int_t nTracksPrim = primVtx->GetNContributors();
1813 fh1VertexNContributors->Fill(nTracksPrim);
1815 if (fDebug > 1) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim);
1817 if(nTracksPrim <= 2){
1818 if (fDebug > 1) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__);
1819 fh1EvtSelection->Fill(3.);
1820 PostData(1, fCommonHistList);
1824 fh1VertexZ->Fill(primVtx->GetZ());
1826 if(TMath::Abs(primVtx->GetZ())>fMaxVertexZ){
1827 if (fDebug > 1) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
1828 fh1EvtSelection->Fill(4.);
1829 PostData(1, fCommonHistList);
1833 // accepts only events that have same "primary" and SPD vertex, special issue of LHC11h PbPb data
1835 //fAOD: pointer to global primary vertex
1837 const AliAODVertex* spdVtx = fAOD->GetPrimaryVertexSPD();
1839 if (TMath::Abs(spdVtx->GetZ() - primVtx->GetZ())>fDeltaVertexZ) { if (fDebug > 1) Printf("deltaZVertex: event REJECTED..."); return;}
1842 //check for vertex radius to be smaller than 1 cm, (that was first applied by Vit Kucera in his analysis)
1844 Double_t vtxX = primVtx->GetX();
1845 Double_t vtxY = primVtx->GetY();
1847 if(TMath::Sqrt(vtxX*vtxX + vtxY*vtxY)>=1){
1848 if (fDebug > 1) Printf("%s:%d primary vertex r = %f: event REJECTED...",(char*)__FILE__,__LINE__,TMath::Sqrt(vtxX*vtxX + vtxY*vtxY));
1853 TString primVtxName(primVtx->GetName());
1855 if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){
1856 if (fDebug > 1) Printf("%s:%d primary vertex selection: TPC vertex, event REJECTED...",(char*)__FILE__,__LINE__);
1857 fh1EvtSelection->Fill(5.);
1858 PostData(1, fCommonHistList);
1862 Bool_t selectedHelper = AliAnalysisHelperJetTasks::Selected();
1863 if(!selectedHelper){
1864 fh1EvtSelection->Fill(6.);
1865 PostData(1, fCommonHistList);
1869 // event selection *****************************************
1871 Double_t centPercent = -1;
1875 if(handler && handler->InheritsFrom("AliAODInputHandler")){
1877 centPercent = fAOD->GetHeader()->GetCentrality();
1879 //std::cout<<"centPercent: "<<centPercent<<std::endl;
1881 fh1EvtAllCent->Fill(centPercent);
1883 if(centPercent>10) cl = 2; //standard PWG-JE binning
1884 if(centPercent>30) cl = 3;
1885 if(centPercent>50) cl = 4;
1889 if(centPercent < 0) cl = -1;
1890 if(centPercent >= 0) cl = 1;
1891 if(centPercent > 10) cl = 2; //standard PWG-JE binning
1892 if(centPercent > 30) cl = 3;
1893 if(centPercent > 50) cl = 4;
1894 if(centPercent > 80) cl = 5; //takes centralities higher than my upper edge of 80%, not to be used
1899 cl = AliAnalysisHelperJetTasks::EventClass();
1901 if(fESD) centPercent = fESD->GetCentrality()->GetCentralityPercentile("V0M"); //ESD JetServices Task has the centrality binning 0-10,10-30,30-50,50-80
1902 fh1EvtAllCent->Fill(centPercent);
1905 if(cl!=fEventClass){ // event not in selected event class, reject event#########################################
1907 if (fDebug > 1) Printf("%s:%d event not in selected event class: event REJECTED ...",(char*)__FILE__,__LINE__);
1908 fh1EvtSelection->Fill(2.);
1909 PostData(1, fCommonHistList);
1912 }//end if fEventClass > 0
1915 if (fDebug > 1) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__);
1918 //Printf("Analysis event #%5d", (Int_t) fEntry);
1920 fh1EvtSelection->Fill(0.);
1921 fh1EvtCent->Fill(centPercent);
1923 //___ get MC information __________________________________________________________________
1926 Double_t ptHard = 0.; //parton energy bins -> energy of particle
1927 Double_t nTrials = 1; // trials for MC trigger weight for real data
1930 AliGenEventHeader* genHeader = fMCEvent->GenEventHeader();
1931 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);//check usage of Pythia (pp) or Hijing (PbPb)
1932 AliGenHijingEventHeader* hijingGenHeader = 0x0;
1934 if(pythiaGenHeader){
1935 if(fDebug>3) Printf("%s:%d pythiaGenHeader found", (char*)__FILE__,__LINE__);
1936 nTrials = pythiaGenHeader->Trials();
1937 ptHard = pythiaGenHeader->GetPtHard();
1939 fh1PtHard->Fill(ptHard);
1940 fh1PtHardTrials->Fill(ptHard,nTrials);
1943 } else { // no pythia, hijing?
1945 if(fDebug>3) Printf("%s:%d no pythiaGenHeader found", (char*)__FILE__,__LINE__);
1947 hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
1948 if(!hijingGenHeader){
1949 if(fDebug>3) Printf("%s:%d no pythiaGenHeader or hjingGenHeader found", (char*)__FILE__,__LINE__);
1951 if(fDebug>3) Printf("%s:%d hijingGenHeader found", (char*)__FILE__,__LINE__);
1955 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
1958 //____ fetch jets _______________________________________________________________
1960 Int_t nJCuts = GetListOfJets(fJetsRecCuts, kJetsRecAcceptance);//fetch list with jets that survived all jet cuts: fJetsRecCuts
1962 Int_t nRecJetsCuts = 0; //number of reconstructed jets after jet cuts
1963 if(nJCuts>=0) nRecJetsCuts = fJetsRecCuts->GetEntries();
1964 if(fDebug>2)Printf("%s:%d Selected Rec jets after cuts: %d %d",(char*)__FILE__,__LINE__,nJCuts,nRecJetsCuts);
1965 if(nRecJetsCuts != nJCuts) Printf("%s:%d Mismatch selected Rec jets after cuts: %d %d",(char*)__FILE__,__LINE__,nJCuts,nRecJetsCuts);
1966 fh1nRecJetsCuts->Fill(nRecJetsCuts);
1969 //____ fetch background clusters ___________________________________________________
1970 if(fBranchRecBckgClusters.Length() != 0){
1972 Int_t nBJ = GetListOfBckgJets(fBckgJetsRec, kJetsRec);
1973 Int_t nRecBckgJets = 0;
1974 if(nBJ>=0) nRecBckgJets = fBckgJetsRec->GetEntries();
1975 if(fDebug>2)Printf("%s:%d Selected Rec background jets: %d %d",(char*)__FILE__,__LINE__,nBJ,nRecBckgJets);
1976 if(nBJ != nRecBckgJets) Printf("%s:%d Mismatch Selected Rec background jets: %d %d",(char*)__FILE__,__LINE__,nBJ,nRecBckgJets);
1980 //____ fetch reconstructed particles __________________________________________________________
1982 Int_t nTCuts = GetListOfTracks(fTracksRecCuts, kTrackAODCuts);//all tracks of event
1983 if(fDebug>2)Printf("%s:%d selected reconstructed tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,fTracksRecCuts->GetEntries());
1984 if(fTracksRecCuts->GetEntries() != nTCuts)
1985 Printf("%s:%d Mismatch selected reconstructed tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,fTracksRecCuts->GetEntries());
1986 fh1EvtMult->Fill(fTracksRecCuts->GetEntries());
1988 Int_t nK0s = GetListOfV0s(fListK0s,fK0Type,kK0,myPrimaryVertex,fAOD);//all V0s in event with K0s assumption
1990 if(fDebug>5){std::cout<<"fK0Type: "<<fK0Type<<" kK0: "<<kK0<<" myPrimaryVertex: "<<myPrimaryVertex<<" fAOD: "<<fAOD<<std::endl;}
1992 //std::cout<< "nK0s: "<<nK0s<<std::endl;
1994 if(fDebug>2)Printf("%s:%d Selected V0s after cuts: %d %d",(char*)__FILE__,__LINE__,nK0s,fListK0s->GetEntries());
1995 if(nK0s != fListK0s->GetEntries()) Printf("%s:%d Mismatch selected K0s: %d %d",(char*)__FILE__,__LINE__,nK0s,fListK0s->GetEntries());
1996 fh1K0Mult->Fill(fListK0s->GetEntries());
1999 Int_t nLa = GetListOfV0s(fListLa,fLaType,kLambda,myPrimaryVertex,fAOD);//all V0s in event with Lambda particle assumption
2000 if(fDebug>2)Printf("%s:%d Selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nLa,fListLa->GetEntries());
2001 if(nLa != fListLa->GetEntries()) Printf("%s:%d Mismatch selected La: %d %d",(char*)__FILE__,__LINE__,nLa,fListLa->GetEntries());
2002 fh1LaMult->Fill(fListLa->GetEntries());
2004 Int_t nALa = GetListOfV0s(fListALa,fALaType,kAntiLambda,myPrimaryVertex,fAOD);//all V0s in event with Antilambda particle assumption
2005 if(fDebug>2)Printf("%s:%d Selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nALa,fListALa->GetEntries());
2006 if(nALa != fListALa->GetEntries()) Printf("%s:%d Mismatch selected ALa: %d %d",(char*)__FILE__,__LINE__,nALa,fListALa->GetEntries());
2007 fh1ALaMult->Fill(fListALa->GetEntries());
2011 //fetch MC gen particles_______________________________________________________
2013 if(fAnalysisMC){ // here
2015 //fill feeddown histo for associated particles
2017 // Access MC generated particles, fill TLists and histograms :
2019 Int_t nMCgenK0s = GetListOfMCParticles(fListMCgenK0s,kK0,fAOD); //fill TList with MC generated primary true K0s (list to fill, particletype, mc aod event)
2020 if(nMCgenK0s != fListMCgenK0s->GetEntries()) Printf("%s:%d Mismatch selected MCgenK0s: %d %d",(char*)__FILE__,__LINE__,nMCgenK0s,fListMCgenK0s->GetEntries());
2023 for(Int_t it=0; it<fListMCgenK0s->GetSize(); ++it){ // loop MC generated K0s, filling histograms
2025 AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenK0s->At(it));
2030 Double_t fRapCurrentPart = MyRapidity(mcp0->E(),mcp0->Pz());
2031 Double_t fEtaCurrentPart = mcp0->Eta();
2032 Double_t fPtCurrentPart = mcp0->Pt();
2034 fh1MCEtaK0s->Fill(fEtaCurrentPart);
2035 fh1MCRapK0s->Fill(fRapCurrentPart);
2036 fh1MCPtK0s->Fill(fPtCurrentPart);
2038 fh2MCEtaVsPtK0s->Fill(fPtCurrentPart,fEtaCurrentPart); //eta cut, physical primary selection and decay mode considered
2043 Int_t nMCgenLa = GetListOfMCParticles(fListMCgenLa,kLambda,fAOD); //fill TList with MC generated primary true Lambdas (list to fill, particletype, mc aod event)
2044 if(nMCgenLa != fListMCgenLa->GetEntries()) Printf("%s:%d Mismatch selected MCgenLa: %d %d",(char*)__FILE__,__LINE__,nMCgenLa,fListMCgenLa->GetEntries());
2047 for(Int_t it=0; it<fListMCgenLa->GetSize(); ++it){ // loop MC generated La, filling histograms
2049 AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenLa->At(it));
2054 Double_t fRapCurrentPart = MyRapidity(mcp0->E(),mcp0->Pz());
2055 Double_t fEtaCurrentPart = mcp0->Eta();
2056 Double_t fPtCurrentPart = mcp0->Pt();
2058 fh1MCEtaLambda->Fill(fEtaCurrentPart);
2059 fh1MCRapLambda->Fill(fRapCurrentPart);
2060 fh1MCPtLambda->Fill(fPtCurrentPart);
2061 fh2MCEtaVsPtLa->Fill(fPtCurrentPart,fEtaCurrentPart); //eta cut, physical primary selection and decay mode considered
2066 Int_t nMCgenALa = GetListOfMCParticles(fListMCgenALa,kAntiLambda,fAOD); //fill TList with MC generated primary true Antilambdas (list to fill, particletype, mc aod event)
2067 if(nMCgenALa != fListMCgenALa->GetEntries()) Printf("%s:%d Mismatch selected MCgenALa: %d %d",(char*)__FILE__,__LINE__,nMCgenALa,fListMCgenALa->GetEntries());
2070 for(Int_t it=0; it<fListMCgenALa->GetSize(); ++it){ // loop MC generated ALa, filling histograms
2072 AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenALa->At(it));
2075 //MC gen Antilambdas
2077 Double_t fRapCurrentPart = MyRapidity(mcp0->E(),mcp0->Pz());
2078 Double_t fEtaCurrentPart = mcp0->Eta();
2079 Double_t fPtCurrentPart = mcp0->Pt();
2081 fh1MCEtaAntiLambda->Fill(fEtaCurrentPart);
2082 fh1MCRapAntiLambda->Fill(fRapCurrentPart);
2083 fh1MCPtAntiLambda->Fill(fPtCurrentPart);
2084 fh2MCEtaVsPtALa->Fill(fPtCurrentPart,fEtaCurrentPart); //eta cut, physical primary selection and decay mode considered
2090 //loop over MC feeddown candidates in TList
2095 } //end MCAnalysis part for gen particles
2098 // ___ V0 QA + K0s + La + ALa pt spectra all events _______________________________________________
2100 Double_t lPrimaryVtxPosition[3];
2101 Double_t lV0Position[3];
2102 lPrimaryVtxPosition[0] = primVtx->GetX();
2103 lPrimaryVtxPosition[1] = primVtx->GetY();
2104 lPrimaryVtxPosition[2] = primVtx->GetZ();
2105 Double_t dRadiusExcludeCone = 2*GetFFRadius(); //2 times jet radius
2106 //------------------------------------------
2107 for(Int_t it=0; it<fListK0s->GetSize(); ++it){
2109 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListK0s->At(it));
2114 // VO's main characteristics to check the reconstruction cuts
2116 //Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2119 Double_t fV0Radius = -999;
2120 Double_t fDcaV0Daughters = v0->DcaV0Daughters();
2121 Double_t fDcaPosToPrimVertex = v0->DcaPosToPrimVertex();//IP of positive charged daughter
2122 Double_t fDcaNegToPrimVertex = v0->DcaNegToPrimVertex();//IP of negative charged daughter
2123 Int_t negDaughterpdg = 0;
2124 Int_t posDaughterpdg = 0;
2125 Int_t motherType = 0;
2128 Bool_t fPhysicalPrimary = kFALSE;//don't use IsPhysicalPrimary() anymore for MC analysis, use instead 2D distance from primary to secondary vertex
2129 Int_t MCv0PdgCode = 0;
2131 AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0));
2132 AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1));
2134 Double_t PosEta = trackPos->AliAODTrack::Eta();//daughter track charge is sometimes wrong here, account for that!!!
2135 Double_t NegEta = trackNeg->AliAODTrack::Eta();
2137 //Double_t trackPosNcls = trackPos->GetTPCNcls();//Get number of clusters for positive charged tracks
2138 //Double_t trackNegNcls = trackNeg->GetTPCNcls();//Get number of clusters for negative charged tracks
2140 CalculateInvMass(v0, kK0, invMK0s, trackPt); //function to calculate invMass with TLorentzVector class
2142 //Double_t jetPt = fFFIMJetPtMin; // assign pro forma jet energy
2143 //Double_t fRap = v0->RapK0Short();
2144 Double_t fEta = v0->PseudoRapV0();
2145 Bool_t bIsInCone = kFALSE;//init boolean, is not in any cone (OC)
2147 for(Int_t ij=0; ij<nRecJetsCuts; ++ij){ // loop over all jets in event
2149 AliAODJet* jet = (AliAODJet*) (fJetsRecCuts->At(ij));
2151 if(IsParticleInCone(jet, v0, dRadiusExcludeCone) == kTRUE) {bIsInCone = kTRUE;}
2155 if(bIsInCone==kFALSE){//K0s is not part of any selected jet in event
2156 Double_t vK0sOC[3] = {invMK0s,trackPt,fEta};
2157 fhnK0sOC->Fill(vK0sOC);
2160 Double_t fV0cosPointAngle = v0->CosPointingAngle(lPrimaryVtxPosition);
2162 lV0Position[0]= v0->DecayVertexV0X();
2163 lV0Position[1]= v0->DecayVertexV0Y();
2164 lV0Position[2]= v0->DecayVertexV0Z();
2168 fV0mom[0]=v0->MomV0X();
2169 fV0mom[1]=v0->MomV0Y();
2170 fV0mom[2]=v0->MomV0Z();
2171 // Double_t fV0TotalMomentum = TMath::Sqrt(fV0mom[0]*fV0mom[0]+fV0mom[1]*fV0mom[1]+fV0mom[2]*fV0mom[2]);
2172 Double_t fV0DecayLength = v0->DecayLengthV0(lPrimaryVtxPosition);
2173 fV0Radius = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
2176 //fetch V0s outside of jet cones (outside of 2R):
2186 fV0QAK0->FillTrackQA(v0->Eta(), TVector2::Phi_0_2pi(v0->Phi()), v0->Pt());
2187 //fFFHistosIMK0AllEvt->FillFF(trackPt, invMK0s, jetPt, incrementJetPt);
2188 //fh1trackPosNCls->Fill(trackPosNcls);
2189 //fh1trackNegNCls->Fill(trackNegNcls);
2190 fh1EtaK0s->Fill(fEta);
2192 Double_t vK0sIncl[3] = {invMK0s,trackPt,fEta}; //fill all K0s in event into THnSparse of 3 dimensions
2193 fhnK0sIncl->Fill(vK0sIncl);
2196 TList *listmc = fAOD->GetList();
2197 Bool_t mclabelcheck = MCLabelCheck(v0, kK0, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode);
2198 //if(fPhysicalPrimary == kFALSE)continue;
2199 //std::cout<<"mclabelcheck: "<<mclabelcheck<<std::endl;
2200 //std::cout<<"IsPhysicalPrimary: "<<fPhysicalPrimary<<std::endl;
2202 if(mclabelcheck == kFALSE)continue;
2204 Double_t vInvMassEtaTrackPtK0s[3] = {fEta,invMK0s,trackPt};
2205 fhnInvMassEtaTrackPtK0s->Fill(vInvMassEtaTrackPtK0s);//includes also feeddown particles, mainly phi particles whose decay products are considered here as primary
2208 fh1PtMCK0s->Fill(MCPt);
2212 fh1V0Eta->Fill(fEta);
2213 //fh1V0totMom->Fill(fV0TotalMomentum);
2214 fh1CosPointAngle->Fill(fV0cosPointAngle);
2215 fh1DecayLengthV0->Fill(fV0DecayLength);
2216 fh1V0Radius->Fill(fV0Radius);
2217 fh1DcaV0Daughters->Fill(fDcaV0Daughters);
2218 fh1DcaPosToPrimVertex->Fill(fDcaPosToPrimVertex);
2219 fh1DcaNegToPrimVertex->Fill(fDcaNegToPrimVertex);
2220 fh1trackPosEta->Fill(PosEta);
2221 fh1trackNegEta->Fill(NegEta);
2225 // __La pt spectra all events _______________________________________________
2228 for(Int_t it=0; it<fListLa->GetSize(); ++it){ // loop all La
2230 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListLa->At(it));
2233 // VO's main characteristics to check the reconstruction cuts
2234 // Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2237 Double_t fV0Radius = -999;
2238 Double_t fDcaV0Daughters = v0->DcaV0Daughters();
2239 Double_t fDcaPosToPrimVertex = v0->DcaPosToPrimVertex();//IP of positive charged daughter
2240 Double_t fDcaNegToPrimVertex = v0->DcaNegToPrimVertex();//IP of negative charged daughter
2241 Int_t negDaughterpdg = 0;
2242 Int_t posDaughterpdg = 0;
2243 Int_t motherType = 0;
2246 Bool_t fPhysicalPrimary = kFALSE;
2247 Int_t MCv0PdgCode = 0;
2248 AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0));
2249 AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1));
2251 //Double_t trackPosNcls = trackPos->GetTPCNcls();//Get number of clusters for positive charged tracks
2252 //Double_t trackNegNcls = trackNeg->GetTPCNcls();//Get number of clusters for negative charged tracks
2254 Double_t PosEta = trackPos->AliAODTrack::Eta();//daughter track charge is sometimes wrong here, account for that!!!
2255 Double_t NegEta = trackNeg->AliAODTrack::Eta();
2257 Double_t fEta = v0->PseudoRapV0();
2258 Bool_t bIsInCone = kFALSE;//init boolean, is not in any cone (OC)
2260 CalculateInvMass(v0, kLambda, invMLa, trackPt);//function to calculate invMass with TLorentzVector class
2263 for(Int_t ij=0; ij<nRecJetsCuts; ++ij){ // loop over all jets in event
2265 AliAODJet* jet = (AliAODJet*) (fJetsRecCuts->At(ij));
2267 if(IsParticleInCone(jet, v0, dRadiusExcludeCone) == kTRUE) {bIsInCone = kTRUE;
2271 if(bIsInCone == kFALSE){//success!
2272 Double_t vLaOC[3] = {invMLa, trackPt,fEta};
2273 fhnLaOC->Fill(vLaOC);
2276 // Double_t jetPt = fFFIMJetPtMin; // assign pro forma jet energy
2277 // Double_t fRap = v0->Y(3122);
2282 fV0mom[0]=v0->MomV0X();
2283 fV0mom[1]=v0->MomV0Y();
2284 fV0mom[2]=v0->MomV0Z();
2285 //Double_t fV0TotalMomentum = TMath::Sqrt(fV0mom[0]*fV0mom[0]+fV0mom[1]*fV0mom[1]+fV0mom[2]*fV0mom[2]);
2286 Double_t fV0DecayLength = v0->DecayLengthV0(lPrimaryVtxPosition);
2287 Double_t fV0cosPointAngle = v0->CosPointingAngle(lPrimaryVtxPosition);
2288 lV0Position[0]= v0->DecayVertexV0X();
2289 lV0Position[1]= v0->DecayVertexV0Y();
2290 lV0Position[2]= v0->DecayVertexV0Z();
2292 fV0Radius = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
2294 //fFFHistosIMLaAllEvt->FillFF(trackPt, invMLa, jetPt, incrementJetPt);
2295 //fh1trackPosNCls->Fill(trackPosNcls);
2296 //fh1trackNegNCls->Fill(trackNegNcls);
2297 fh1EtaLa->Fill(fEta);
2299 Double_t vLaIncl[3] = {invMLa,trackPt,fEta};
2300 fhnLaIncl->Fill(vLaIncl);
2303 TList* listmc = fAOD->GetList();
2304 Bool_t mclabelcheck = MCLabelCheck(v0, kLambda, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode);
2305 if(mclabelcheck == kFALSE)continue;
2306 //if(fPhysicalPrimary == kFALSE)continue;
2308 Double_t vInvMassEtaTrackPtLa[3] = {fEta,invMLa,trackPt};
2309 fhnInvMassEtaTrackPtLa->Fill(vInvMassEtaTrackPtLa);//includes also feed-down particles
2310 fh1PtMCLa->Fill(MCPt);
2313 fh1PtMCLa->Fill(MCPt);
2315 fh1V0Eta->Fill(fEta);
2316 //fh1V0totMom->Fill(fV0TotalMomentum);
2317 fh1CosPointAngle->Fill(fV0cosPointAngle);
2318 fh1DecayLengthV0->Fill(fV0DecayLength);
2319 fh1V0Radius->Fill(fV0Radius);
2320 fh1DcaV0Daughters->Fill(fDcaV0Daughters);
2321 fh1DcaPosToPrimVertex->Fill(fDcaPosToPrimVertex);
2322 fh1DcaNegToPrimVertex->Fill(fDcaNegToPrimVertex);
2323 fh1trackPosEta->Fill(PosEta);
2324 fh1trackNegEta->Fill(NegEta);
2327 // __ALa pt spectra all events _______________________________________________
2329 for(Int_t it=0; it<fListALa->GetSize(); ++it){ // loop all ALa
2331 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListALa->At(it));
2335 //VO's main characteristics to check the reconstruction cuts
2336 Double_t invMALa =0;
2338 Double_t fV0Radius = -999;
2339 Double_t fDcaV0Daughters = v0->DcaV0Daughters();
2340 Double_t fDcaPosToPrimVertex = v0->DcaPosToPrimVertex();//IP of positive charged daughter
2341 Double_t fDcaNegToPrimVertex = v0->DcaNegToPrimVertex();//IP of negative charged daughter
2342 Int_t negDaughterpdg = 0;
2343 Int_t posDaughterpdg = 0;
2344 Int_t motherType = 0;
2347 Bool_t fPhysicalPrimary = kFALSE;
2348 Int_t MCv0PdgCode = 0;
2350 AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0));
2351 AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1));
2353 Double_t PosEta = trackPos->AliAODTrack::Eta();//daughter track charge is sometimes wrong here, account for that!!!
2354 Double_t NegEta = trackNeg->AliAODTrack::Eta();
2356 Double_t fEta = v0->PseudoRapV0();
2357 Bool_t bIsInCone = kFALSE;//init boolean for OC
2360 CalculateInvMass(v0, kAntiLambda, invMALa, trackPt); //function to calculate invMass with TLorentzVector class
2362 for(Int_t ij=0; ij<nRecJetsCuts; ++ij){ // loop over all jets in event
2364 AliAODJet* jet = (AliAODJet*) (fJetsRecCuts->At(ij));
2366 if(IsParticleInCone(jet, v0, dRadiusExcludeCone) == kTRUE){
2371 if(bIsInCone == kFALSE){//success!
2372 Double_t vALaOC[3] = {invMALa, trackPt,fEta};
2373 fhnALaOC->Fill(vALaOC);
2376 //Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2377 //Double_t jetPt = fFFIMJetPtMin; // assign pro forma jet energy
2378 // Double_t fRap = v0->Y(-3122);
2383 fV0mom[0]=v0->MomV0X();
2384 fV0mom[1]=v0->MomV0Y();
2385 fV0mom[2]=v0->MomV0Z();
2386 //Double_t fV0TotalMomentum = TMath::Sqrt(fV0mom[0]*fV0mom[0]+fV0mom[1]*fV0mom[1]+fV0mom[2]*fV0mom[2]);
2388 Double_t fV0cosPointAngle = v0->CosPointingAngle(lPrimaryVtxPosition);
2389 lV0Position[0]= v0->DecayVertexV0X();
2390 lV0Position[1]= v0->DecayVertexV0Y();
2391 lV0Position[2]= v0->DecayVertexV0Z();
2392 Double_t fV0DecayLength = v0->DecayLengthV0(lPrimaryVtxPosition);
2393 fV0Radius = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
2395 //fFFHistosIMALaAllEvt->FillFF(trackPt, invMALa, jetPt, incrementJetPt);
2396 //fh1trackPosNCls->Fill(trackPosNcls);
2397 //fh1trackNegNCls->Fill(trackNegNcls);
2398 fh1EtaALa->Fill(fEta);
2400 Double_t vALaIncl[3] = {invMALa,trackPt,fEta};
2401 fhnALaIncl->Fill(vALaIncl);
2404 TList* listmc = fAOD->GetList();
2405 Bool_t mclabelcheck = MCLabelCheck(v0, kAntiLambda, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode);
2406 if(mclabelcheck == kFALSE)continue;
2407 //if(fPhysicalPrimary == kFALSE)continue;
2409 Double_t vInvMassEtaTrackPtALa[3] = {fEta,invMALa,trackPt};
2410 fhnInvMassEtaTrackPtALa->Fill(vInvMassEtaTrackPtALa);
2411 fh1PtMCALa->Fill(MCPt);
2414 fh1V0Eta->Fill(fEta);
2415 //fh1V0totMom->Fill(fV0TotalMomentum);
2416 fh1CosPointAngle->Fill(fV0cosPointAngle);
2417 fh1DecayLengthV0->Fill(fV0DecayLength);
2418 fh1V0Radius->Fill(fV0Radius);
2419 fh1DcaV0Daughters->Fill(fDcaV0Daughters);
2420 fh1DcaPosToPrimVertex->Fill(fDcaPosToPrimVertex);
2421 fh1DcaNegToPrimVertex->Fill(fDcaNegToPrimVertex);
2422 fh1trackPosEta->Fill(PosEta);
2423 fh1trackNegEta->Fill(NegEta);
2426 //_____no jets events______________________________________________________________________________________________________________________________________
2428 if(nRecJetsCuts == 0){//no jet events
2430 if(fDebug>6) { std::cout<<"################## nRecJetsCuts == 0 ###################"<<std::endl;
2431 std::cout<<"fListK0s->GetSize() in NJ event: "<<fListK0s->GetSize()<<std::endl;
2434 for(Int_t it=0; it<fListK0s->GetSize(); ++it){ // loop all K0s
2436 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListK0s->At(it));
2439 Double_t invMK0s =0;
2441 CalculateInvMass(v0, kK0, invMK0s, trackPt);
2442 fh2NJK0->Fill(invMK0s, trackPt);
2446 for(Int_t it=0; it<fListLa->GetSize(); ++it){ // loop all La
2448 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListLa->At(it));
2453 CalculateInvMass(v0, kLambda, invMLa, trackPt);
2454 fh2NJLa->Fill(invMLa, trackPt);
2458 for(Int_t it=0; it<fListALa->GetSize(); ++it){ // loop all ALa
2460 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListALa->At(it));
2463 Double_t invMALa =0;
2465 CalculateInvMass(v0, kAntiLambda, invMALa, trackPt);
2466 fh2NJALa->Fill(invMALa, trackPt);
2472 //____ fill all jet related histos ________________________________________________________________________________________________________________________
2473 //##########################jet loop########################################################################################################################
2475 //fill jet histos in general
2476 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
2478 AliAODJet* jet = (AliAODJet*) (fJetsRecCuts->At(ij));
2480 Double_t jetPt = jet->Pt();
2481 Double_t jetEta = jet->Eta();
2482 Double_t jetPhi = jet->Phi();
2484 //if(ij==0){ // loop over leading jets for ij = 0, for ij>= 0 look into all jets
2486 if(ij>=0){//all jets in event
2488 TList* jettracklist = new TList();
2489 Double_t sumPt = 0.;
2490 Bool_t isBadJet = kFALSE;
2491 Int_t njetTracks = 0;
2493 if(GetFFRadius()<=0){
2494 GetJetTracksTrackrefs(jettracklist, jet, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJet);// list of jet tracks from trackrefs
2496 GetJetTracksPointing(fTracksRecCuts, jettracklist, jet, GetFFRadius(), sumPt, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJet); // fill list of tracks in cone around jet axis with cone Radius (= 0.4 standard)
2499 if(GetFFMinNTracks()>0 && jettracklist->GetSize() <= GetFFMinNTracks()) isBadJet = kTRUE; // reject jets with less tracks than fFFMinNTracks
2500 if(isBadJet) continue; // rejects jets in which no track has a track pt higher than 5 GeV/c (see AddTask macro)
2502 //Float_t fJetAreaMin = 0.6*TMath::Pi()*GetFFRadius()*GetFFRadius(); // minimum jet area cut
2504 //std::cout<<"GetFFRadius(): "<<GetFFRadius()<<std::endl;
2505 //std::cout<<"jet->EffectiveAreaCharged()"<<jet->EffectiveAreaCharged()<<std::endl;
2506 //std::cout<<"fJetAreaMin: "<<fJetAreaMin<<std::endl;
2508 // if(fDebug > 2) {if (jet->EffectiveAreaCharged() < fJetAreaMin) {std::cout<<" fCutjetArea cut removed a jet!!!!! Should not have to be done again!!"<<std::endl;}}// cut on jet area, already done by jet selection in FF task
2510 Double_t dAreaExcluded = TMath::Pi()*dRadiusExcludeCone*dRadiusExcludeCone; // area of the cone
2511 dAreaExcluded -= AreaCircSegment(dRadiusExcludeCone,fCutjetEta-jet->Eta()); // positive eta overhang
2512 dAreaExcluded -= AreaCircSegment(dRadiusExcludeCone,fCutjetEta+jet->Eta()); // negative eta overhang
2513 fh1AreaExcluded->Fill(dAreaExcluded);//histo contains all areas that are jet related and have to be excluded concerning OC UE pt spectrum normalisation by area
2515 fh1JetEta->Fill(jetEta);
2516 fh1JetPhi->Fill(jetPhi);
2517 fh2JetEtaPhi->Fill(jetEta,jetPhi);
2519 // printf("pT = %f, eta = %f, phi = %f, leadtr pt = %f\n, ",jetPt,jetEta,jetphi,leadtrack);
2521 for(Int_t it=0; it<jettracklist->GetSize(); ++it){//loop over all particles in jet
2523 AliVParticle* trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));//all tracks in jet cone
2524 if(!trackVP)continue;
2526 Float_t trackPt = trackVP->Pt();//transversal momentum of jet particle
2527 Float_t trackEta = trackVP->Eta();
2529 Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2531 fFFHistosRecCuts->FillFF(trackPt, jetPt, incrementJetPt);//histo with tracks/jets after cut selection, for all events
2532 if(nK0s>0) fFFHistosRecCutsK0Evt->FillFF(trackPt, jetPt, incrementJetPt);//only for K0s events
2533 fh2FFJetTrackEta->Fill(trackEta,jetPt);
2538 njetTracks = jettracklist->GetSize();
2540 //____________________________________________________________________________________________________________________
2541 //strangeness constribution to jet cone
2545 TList *list = fAOD->GetList();
2546 AliAODMCHeader *mcHeadr=(AliAODMCHeader*)list->FindObject(AliAODMCHeader::StdBranchName());
2547 if(!mcHeadr)continue;
2549 Double_t mcXv=0., mcYv=0., mcZv=0.;//MC primary vertex position
2551 mcXv=mcHeadr->GetVtxX(); mcYv=mcHeadr->GetVtxY(); mcZv=mcHeadr->GetVtxZ(); // position of the MC primary vertex
2553 for(Int_t it=0; it<jettracklist->GetSize(); ++it){//loop over all tracks in the jet
2555 AliVParticle* trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));//track in jet cone
2556 if(!trackVP)continue;
2557 AliAODTrack *tr = dynamic_cast<AliAODTrack*> (trackVP); //fetch one jet track from the TList
2560 //get MC label information
2561 TList *mclist = fAOD->GetList();
2563 //fetch the MC stack
2564 TClonesArray* stackMC = (TClonesArray*)mclist->FindObject(AliAODMCParticle::StdBranchName()); //get MCAOD branch in data
2565 if (!stackMC) {Printf("ERROR: stack not available");}
2569 Int_t trackLabel = TMath::Abs(tr->GetLabel()); //fetch jet track label in MC stack
2571 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(stackMC->At(trackLabel)); //fetch MC gen. particle for rec. jet track
2573 if(!part)continue; //skip non-existing objects
2576 //Bool_t IsPhysicalPrimary = part->IsPhysicalPrimary();//not recommended to check, better use distance between primary vertex and secondary vertex
2578 Float_t fDistPrimaryMax = 0.01;
2579 // Get the distance between production point of the MC mother particle and the primary vertex
2581 Double_t dx = mcXv-part->Xv();//mc primary vertex - mc gen. v0 vertex
2582 Double_t dy = mcYv-part->Yv();
2583 Double_t dz = mcZv-part->Zv();
2585 Float_t fDistPrimary = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
2586 Bool_t fPhysicalPrimary = (fDistPrimary < fDistPrimaryMax);
2588 // std::cout<<"fDistPrimary"<<fDistPrimary<<std::endl;
2589 // std::cout<<"fPhysicalPrimary"<<fPhysicalPrimary<<std::endl;
2591 if(!fPhysicalPrimary)continue;//rejects Kstar and other strong decaying particles from Secondary Contamination
2593 Bool_t isFromStrange = kFALSE;// flag to check whether particle has strange mother
2595 Int_t iMother = part->GetMother(); //get mother MC gen. particle label
2598 AliAODMCParticle *partM = dynamic_cast<AliAODMCParticle*>(stackMC->At(iMother)); //fetch mother of MC gen. particle
2599 if(!partM) continue;
2601 Int_t codeM = TMath::Abs(partM->GetPdgCode()); //mothers pdg code
2603 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..)
2605 if (mfl == 3 && codeM != 3) isFromStrange = kTRUE;
2608 if(isFromStrange == kTRUE){
2610 Double_t trackPt = part->Pt();
2611 Double_t trackEta = part->Eta();
2612 //fh3StrContinCone->Fill(jetPt, trackPt, trackEta);//MC gen. particle parameters, but rec. jet pt
2614 }//isFromStrange is kTRUE
2616 }//end loop over jet tracks
2621 fh1TrackMultCone->Fill(njetTracks);
2622 fh2TrackMultCone->Fill(njetTracks,jetPt);
2626 // fQAJetHistosRecCutsLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt() );
2628 for(Int_t it=0; it<fListK0s->GetSize(); ++it){ // loop all K0s
2630 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListK0s->At(it));
2631 if(!v0) continue;//rejection of events with no V0 vertex
2635 TVector3 v0MomVect(v0Mom);
2637 Double_t dPhiJetK0 = (jet->MomentumVector()->Vect()).DeltaPhi(v0MomVect);
2638 // Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2640 // if(incrementJetPt){fh1V0JetPt->Fill(jetPt);}
2642 Double_t invMK0s =0;
2644 CalculateInvMass(v0, kK0, invMK0s, trackPt); //function to calculate invMass with TLorentzVector class
2646 // fFFHistosIMK0Jet->FillFF(trackPt, invMK0s, jetPt, incrementJetPt);
2649 if(dPhiJetK0<fh1dPhiJetK0->GetXaxis()->GetXmin()) dPhiJetK0 += 2*TMath::Pi();
2650 fh1dPhiJetK0->Fill(dPhiJetK0);
2654 // if(fListK0s->GetSize() == 0){ // no K0: increment jet pt spectrum
2656 // Bool_t incrementJetPt = kTRUE;
2657 // fFFHistosIMK0Jet->FillFF(-1, -1, jetPt, incrementJetPt);
2660 //____fetch reconstructed K0s in cone around jet axis:_______________________________________________________________________________
2662 TList* jetConeK0list = new TList();
2664 Double_t sumPtK0 = 0.;
2666 Bool_t isBadJetK0 = kFALSE; // dummy, do not use
2669 GetTracksInCone(fListK0s, jetConeK0list, jet, GetFFRadius(), sumPtK0, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetK0); //reconstructed K0s in cone around jet axis
2671 if(fDebug>2)Printf("%s:%d nK0s total: %d, in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nK0s,jetConeK0list->GetEntries(),GetFFRadius());
2674 for(Int_t it=0; it<jetConeK0list->GetSize(); ++it){ // loop for K0s in jet cone
2676 AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetConeK0list->At(it));
2679 Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2680 Double_t invMK0s =0;
2685 CalculateInvMass(v0, kK0, invMK0s, trackPt); //function to calculate invMass with TLorentzVector class
2689 Double_t jetPtSmear = -1;
2690 SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
2691 if(incrementJetPt == kTRUE){fh1IMK0ConeSmear->Fill(jetPtSmear);} //fill TH1F for normalization purposes
2694 if(incrementJetPt==kTRUE){
2695 fh1IMK0Cone->Fill(jetPt);}//normalisation by number of selected jets
2697 //fFFHistosIMK0Cone->FillFF(trackPt, invMK0s, jetPt, incrementJetPt);
2699 Double_t vK0sCone[4] = {jetPt, invMK0s,trackPt,fEta};
2700 fhnK0sCone->Fill(vK0sCone);
2704 if(jetConeK0list->GetSize() == 0){ // no K0: increment jet pt spectrum
2707 Bool_t incrementJetPt = kTRUE;//jets without K0s will be only filled in TH1F only once, so no increment needed
2708 //fFFHistosIMK0Cone->FillFF(-1, -1, jetPt, incrementJetPt);
2709 Double_t vK0sCone[4] = {jetPt, -1, -1, -1};
2710 fhnK0sCone->Fill(vK0sCone);
2712 if(incrementJetPt==kTRUE){
2713 fh1IMK0Cone->Fill(jetPt);}//normalisation by number of selected jets
2716 Double_t jetPtSmear = -1;
2717 SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
2718 if(incrementJetPt == kTRUE){fh1IMK0ConeSmear->Fill(jetPtSmear);} //fill TH1F for normalization purposes
2722 //Random cones________________________________________________________________________
2724 if(ij==0){//fetch random cone V0s only once per event
2726 //______fetch random cones___________________________________________________________
2729 AliAODJet* jetRC = 0;
2730 jetRC = GetRandomCone(fJetsRecCuts, fCutjetEta, 2*GetFFRadius());//fetch one random cone for each event
2731 TList* fListK0sRC = new TList();//list for K0s in random cone (RC), one RC per event
2732 TList* fListLaRC = new TList();
2733 TList* fListALaRC = new TList();
2735 Double_t sumPtK0sRC = 0;
2736 Double_t sumPtLaRC = 0;
2737 Double_t sumPtALaRC = 0;
2738 Bool_t isBadJetK0sRC = kFALSE;
2739 Bool_t isBadJetLaRC = kFALSE;
2740 Bool_t isBadJetALaRC = kFALSE;
2745 GetTracksInCone(fListK0s, fListK0sRC, jetRC, GetFFRadius(), sumPtK0sRC, 0, 0, isBadJetK0sRC);
2747 for(Int_t it=0; it<fListK0sRC->GetSize(); ++it){ // loop for K0s in random cone
2749 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListK0sRC->At(it));
2752 Double_t invMK0s =0;
2757 CalculateInvMass(v0, kK0, invMK0s, trackPt); //function to calculate invMass with TLorentzVector class
2759 Double_t vK0sRC[3] = {invMK0s,trackPt,fEta};
2760 fhnK0sRC->Fill(vK0sRC);
2765 GetTracksInCone(fListLa, fListLaRC, jetRC, GetFFRadius(), sumPtLaRC, 0, 0, isBadJetLaRC);
2767 for(Int_t it=0; it<fListLaRC->GetSize(); ++it){ // loop for Lambdas in random cone
2769 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListLaRC->At(it));
2777 CalculateInvMass(v0, kLambda, invMLa, trackPt); //function to calculate invMass with TLorentzVector class
2779 Double_t vLaRC[3] = {invMLa,trackPt,fEta};
2780 fhnLaRC->Fill(vLaRC);
2785 GetTracksInCone(fListALa, fListALaRC, jetRC, GetFFRadius(), sumPtALaRC, 0, 0, isBadJetALaRC);
2787 for(Int_t it=0; it<fListALaRC->GetSize(); ++it){ // loop for Lambdas in random cone
2789 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListALaRC->At(it));
2792 Double_t invMALa =0;
2797 CalculateInvMass(v0, kAntiLambda, invMALa, trackPt); //function to calculate invMass with TLorentzVector class
2799 Double_t vALaRC[3] = {invMALa,trackPt,fEta};
2800 fhnALaRC->Fill(vALaRC);
2810 //fetch particles in perpendicular cone to estimate UE event contribution to particle spectrum
2811 //these perpendicular cone particle spectra serve to subtract the particles in jet cones, that are stemming from the Underlying event, on a statistical basis
2812 //for normalization the common jet pT spectrum is used: fh1IMK0Cone, fh1IMLaCone and fh1IMALaCone
2814 //____fetch reconstructed K0s in cone perpendicular to jet axis:_______________________________________________________________________________
2816 TList* jetPerpConeK0list = new TList();
2818 Double_t sumPerpPtK0 = 0.;
2820 GetTracksInPerpCone(fListK0s, jetPerpConeK0list, jet, GetFFRadius(), sumPerpPtK0); //reconstructed K0s in cone around jet axis
2822 if(fDebug>2)Printf("%s:%d nK0s total: %d, in perp jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nK0s,jetPerpConeK0list->GetEntries(),GetFFRadius());
2824 for(Int_t it=0; it<jetPerpConeK0list->GetSize(); ++it){ // loop for K0s in perpendicular cone
2826 AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetPerpConeK0list->At(it));
2829 Double_t invMPerpK0s =0;
2834 CalculateInvMass(v0, kK0, invMPerpK0s, trackPt); //function to calculate invMass with TLorentzVector class
2835 Double_t vK0sPC[4] = {jetPt, invMPerpK0s,trackPt,fEta};
2837 fhnK0sPC->Fill(vK0sPC); //(x,y,z) //pay attention, this histogram contains the V0 content of both (+/- 90 degrees) perp. cones!!
2842 if(jetPerpConeK0list->GetSize() == 0){ // no K0s in jet cone
2844 Double_t vK0sPC[4] = {jetPt, -1, -1 , -999};//default values for case: no K0s is found in PC
2845 fhnK0sPC->Fill(vK0sPC);
2849 if(ij==0){//median cluster only once for event
2851 AliAODJet* medianCluster = GetMedianCluster();
2854 // ____ rec K0s in median cluster___________________________________________________________________________________________________________
2856 TList* jetMedianConeK0list = new TList();
2857 TList* jetMedianConeLalist = new TList();
2858 TList* jetMedianConeALalist = new TList();
2861 Double_t medianEta = medianCluster->Eta();
2863 if(TMath::Abs(medianEta)<=fCutjetEta){
2865 fh1MedianEta->Fill(medianEta);
2866 fh1JetPtMedian->Fill(jetPt); //for normalisation by total number of median cluster jets
2868 Double_t sumMedianPtK0 = 0.;
2870 Bool_t isBadJetK0Median = kFALSE; // dummy, do not use
2872 GetTracksInCone(fListK0s, jetMedianConeK0list, medianCluster, GetFFRadius(), sumMedianPtK0, 0., 0., isBadJetK0Median); //reconstructed K0s in median cone around jet axis
2873 //GetTracksInCone(fListK0s, jetConeK0list, jet, GetFFRadius(), sumPtK0, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetK0); //original use of function
2875 //cut parameters from Fragmentation Function task:
2876 //Float_t fFFMinLTrackPt; // reject jets with leading track with pt smaller than this value, use GetFFMinLTrackPt()
2877 //Float_t fFFMaxTrackPt; // reject jetscontaining any track with pt larger than this value, use GetFFMaxTrackPt()
2879 for(Int_t it=0; it<jetMedianConeK0list->GetSize(); ++it){ // loop for K0s in median cone
2881 AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetMedianConeK0list->At(it));
2884 Double_t invMMedianK0s =0;
2889 CalculateInvMass(v0, kK0, invMMedianK0s, trackPt); //function to calculate invMass with TLorentzVector class
2890 Double_t vK0sMCC[3] = {invMMedianK0s,trackPt,fEta};
2891 fhnK0sMCC->Fill(vK0sMCC);
2895 if(jetMedianConeK0list->GetSize() == 0){ // no K0s in median cluster cone
2897 Double_t vK0sMCC[3] = {-1, -1, -999};
2898 fhnK0sMCC->Fill(vK0sMCC);
2902 //__________________________________________________________________________________________________________________________________________
2903 // ____ rec Lambdas in median cluster___________________________________________________________________________________________________________
2905 Double_t sumMedianPtLa = 0.;
2906 Bool_t isBadJetLaMedian = kFALSE; // dummy, do not use
2908 GetTracksInCone(fListLa, jetMedianConeLalist, medianCluster, GetFFRadius(), sumMedianPtLa, 0, 0, isBadJetLaMedian); //reconstructed Lambdas in median cone around jet axis
2910 //cut parameters from Fragmentation Function task:
2911 //Float_t fFFMinLTrackPt; // reject jets with leading track with pt smaller than this value, use GetFFMinLTrackPt()
2912 //Float_t fFFMaxTrackPt; // reject jets containing any track with pt larger than this value, use GetFFMaxTrackPt()
2914 for(Int_t it=0; it<jetMedianConeLalist->GetSize(); ++it){ // loop for Lambdas in perpendicular cone
2916 AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetMedianConeLalist->At(it));
2919 Double_t invMMedianLa =0;
2924 CalculateInvMass(v0, kLambda, invMMedianLa, trackPt); //function to calculate invMass with TLorentzVector class
2926 Double_t vLaMCC[4] = {jetPt, invMMedianLa,trackPt,fEta};
2927 fhnLaMCC->Fill(vLaMCC);
2930 if(jetMedianConeLalist->GetSize() == 0){ // no Lambdas in median cluster cone
2932 Double_t vLaMCC[4] = {jetPt, -1, -1, -999};
2933 fhnLaMCC->Fill(vLaMCC);
2938 // ____ rec Antilambdas in median cluster___________________________________________________________________________________________________________
2941 Double_t sumMedianPtALa = 0.;
2943 Bool_t isBadJetALaMedian = kFALSE; // dummy, do not use
2945 GetTracksInCone(fListALa, jetMedianConeALalist, medianCluster, GetFFRadius(), sumMedianPtALa, 0, 0, isBadJetALaMedian); //reconstructed Antilambdas in median cone around jet axis
2948 //cut parameters from Fragmentation Function task:
2949 //Float_t fFFMinLTrackPt; // reject jets with leading track with pt smaller than this value, use GetFFMinLTrackPt()
2950 //Float_t fFFMaxTrackPt; // reject jets containing any track with pt larger than this value, use GetFFMaxTrackPt()
2952 for(Int_t it=0; it<jetMedianConeALalist->GetSize(); ++it){ // loop for Antilambdas in median cluster cone
2954 AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetMedianConeALalist->At(it));
2957 Double_t invMMedianALa =0;
2963 CalculateInvMass(v0, kAntiLambda, invMMedianALa, trackPt); //function to calculate invMass with TLorentzVector class
2964 Double_t vALaMCC[4] = {jetPt, invMMedianALa,trackPt,fEta};
2965 fhnALaMCC->Fill(vALaMCC);
2969 if(jetMedianConeALalist->GetSize() == 0){ // no Antilambdas in median cluster cone
2971 Double_t vALaMCC[4] = {jetPt, -1, -1, -999};
2972 fhnALaMCC->Fill(vALaMCC);
2975 }//median cluster eta cut
2977 delete jetMedianConeK0list;
2978 delete jetMedianConeLalist;
2979 delete jetMedianConeALalist;
2981 }//if mediancluster is existing
2983 //_________________________________________________________________________________________________________________________________________
2985 //____fetch reconstructed Lambdas in cone perpendicular to jet axis:__________________________________________________________________________
2987 TList* jetPerpConeLalist = new TList();
2989 Double_t sumPerpPtLa = 0.;
2991 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!!
2993 if(fDebug>2)Printf("%s:%d nLa total: %d, in perp jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nLa,jetPerpConeLalist->GetEntries(),GetFFRadius());
2995 for(Int_t it=0; it<jetPerpConeLalist->GetSize(); ++it){ // loop for Lambdas in perpendicular cone
2997 AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetPerpConeLalist->At(it));
3000 Double_t invMPerpLa =0;
3005 CalculateInvMass(v0, kLambda, invMPerpLa, trackPt); //function to calculate invMass with TLorentzVector class
3006 Double_t vLaPC[4] = {jetPt, invMPerpLa,trackPt,fEta};
3007 fhnLaPC->Fill(vLaPC); //(x,y,z) //pay attention, this histogram contains the V0 content of both (+/- 90 degrees) perp. cones!!
3012 if(jetPerpConeLalist->GetSize() == 0){ // no Lambdas in jet
3014 Double_t vLaPC[4] = {jetPt, -1, -1 , -999};//default values for case: no K0s is found in PC
3015 fhnLaPC->Fill(vLaPC);
3021 //____fetch reconstructed Antilambdas in cone perpendicular to jet axis:___________________________________________________________________
3023 TList* jetPerpConeALalist = new TList();
3025 Double_t sumPerpPtALa = 0.;
3027 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!!
3029 if(fDebug>2)Printf("%s:%d nALa total: %d, in perp jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nALa,jetPerpConeALalist->GetEntries(),GetFFRadius());
3031 for(Int_t it=0; it<jetPerpConeALalist->GetSize(); ++it){ // loop for ALa in perpendicular cone
3033 AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetPerpConeALalist->At(it));
3036 Double_t invMPerpALa =0;
3041 CalculateInvMass(v0, kAntiLambda, invMPerpALa, trackPt); //function to calculate invMass with TLorentzVector class
3042 Double_t vALaPC[4] = {jetPt, invMPerpALa,trackPt,fEta};
3043 fhnALaPC->Fill(vALaPC);
3048 if(jetPerpConeALalist->GetSize() == 0){ // no Antilambda
3050 Double_t vALaPC[4] = {jetPt, -1, -1, -999};
3051 fhnALaPC->Fill(vALaPC);
3071 //###########################################################################################################
3073 //__________________________________________________________________________________________________________________________________________
3077 //fill feeddown candidates from TList
3078 //std::cout<<"fListFeeddownLaCand entries: "<<fListFeeddownLaCand->GetSize()<<std::endl;
3080 Double_t sumPtFDLa = 0.;
3081 Bool_t isBadJetFDLa = kFALSE; // dummy, do not use
3083 GetTracksInCone(fListFeeddownLaCand, jetConeFDLalist, jet, GetFFRadius(), sumPtFDLa, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetFDLa);
3085 Double_t sumPtFDALa = 0.;
3086 Bool_t isBadJetFDALa = kFALSE; // dummy, do not use
3088 GetTracksInCone(fListFeeddownALaCand, jetConeFDALalist, jet, GetFFRadius(), sumPtFDALa, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetFDALa);
3090 //_________________________________________________________________
3091 for(Int_t it=0; it<fListFeeddownLaCand->GetSize(); ++it){
3093 AliAODv0* mcfd = dynamic_cast<AliAODv0*>(fListFeeddownLaCand->At(it));
3096 Double_t invMLaFDcand = 0;
3097 Double_t trackPt = 0;//pt of ass. particle, not used for the histos
3099 CalculateInvMass(mcfd, kLambda, invMLaFDcand, trackPt);
3101 //Get MC gen. Lambda transverse momentum
3102 TClonesArray *st = 0x0;
3105 TList *lt = fAOD->GetList();
3108 st = (TClonesArray*)lt->FindObject(AliAODMCParticle::StdBranchName()); //get MCAOD branch in data
3111 AliAODTrack *daughtertrack = (AliAODTrack *) (mcfd->GetSecondaryVtx()->GetDaughter(0));//fetch the first of the two daughter tracks
3112 Int_t AssLabel = TMath::Abs(daughtertrack->GetLabel());
3114 AliAODMCParticle *mcDaughterPart =(AliAODMCParticle*)st->UncheckedAt(AssLabel);
3116 Int_t v0lab = mcDaughterPart->GetMother();
3118 // Int_t v0lab= TMath::Abs(mcfd->GetLabel());//GetLabel doesn't work for AliAODv0 class!!! Only for AliAODtrack
3120 if((!v0lab) || (v0lab<0) || (v0lab > st->GetEntriesFast()))continue;//validity checks
3122 AliAODMCParticle *mcp=(AliAODMCParticle*)st->UncheckedAt(v0lab);
3124 Double_t genLaPt = mcp->Pt();
3126 //std::cout<<"Incl FD, genLaPt:"<<genLaPt<<std::endl;
3128 Double_t vFeedDownLa[3] = {5., invMLaFDcand, genLaPt};
3129 fhnFeedDownLa->Fill(vFeedDownLa);
3132 }//end loop over feeddown candidates for Lambda particles in jet cone
3133 //fetch MC truth in jet cones, denominator of rec. efficiency in jet cones
3134 //_________________________________________________________________
3135 for(Int_t it=0; it<jetConeFDLalist->GetSize(); ++it){
3137 AliAODv0* mcfd = dynamic_cast<AliAODv0*>(jetConeFDLalist->At(it));
3140 //std::cout<<"Cone, recLaPt:"<<mcfd->Pt()<<std::endl;
3142 Double_t invMLaFDcand = 0;
3143 Double_t trackPt = mcfd->Pt();//pt of ass. particle, not used for the histos
3145 CalculateInvMass(mcfd, kLambda, invMLaFDcand, trackPt);
3147 //Get MC gen. Lambda transverse momentum
3148 TClonesArray *st = 0x0;
3150 TList *lt = fAOD->GetList();
3153 st = (TClonesArray*)lt->FindObject(AliAODMCParticle::StdBranchName());
3155 AliAODTrack *daughtertrack = (AliAODTrack *) (mcfd->GetSecondaryVtx()->GetDaughter(0));//fetch the first of the two daughter tracks
3156 Int_t AssLabel = TMath::Abs(daughtertrack->GetLabel());
3158 AliAODMCParticle *mcDaughterPart =(AliAODMCParticle*)st->UncheckedAt(AssLabel);
3160 Int_t v0lab = mcDaughterPart->GetMother();
3162 //std::cout<<"v0lab: "<<v0lab<<std::endl;
3164 if((!v0lab) || (v0lab<0) || (v0lab > st->GetEntriesFast()))continue;//validity checks
3166 AliAODMCParticle *mcp=(AliAODMCParticle*)st->UncheckedAt(v0lab);
3168 Double_t genLaPt = mcp->Pt();
3171 //std::cout<<"Cone FD, genLaPt:"<<genLaPt<<std::endl;
3173 Double_t vFeedDownLaCone[3] = {jetPt, invMLaFDcand, genLaPt};
3174 fhnFeedDownLaCone->Fill(vFeedDownLaCone);
3177 }//end loop over feeddown candidates for Lambda particles in jet cone
3179 //_________________________________________________________________
3180 for(Int_t it=0; it<fListFeeddownALaCand->GetSize(); ++it){
3182 AliAODv0* mcfd = dynamic_cast<AliAODv0*>(fListFeeddownALaCand->At(it));
3185 Double_t invMALaFDcand = 0;
3186 Double_t trackPt = 0;//pt of ass. particle, not used for the histos
3188 CalculateInvMass(mcfd, kAntiLambda, invMALaFDcand, trackPt);
3190 //Get MC gen. Antilambda transverse momentum
3191 TClonesArray *st = 0x0;
3193 TList *lt = fAOD->GetList();
3196 st = (TClonesArray*)lt->FindObject(AliAODMCParticle::StdBranchName());
3198 AliAODTrack *daughtertrack = (AliAODTrack *) (mcfd->GetSecondaryVtx()->GetDaughter(0));//fetch the first of the two daughter tracks
3199 Int_t AssLabel = TMath::Abs(daughtertrack->GetLabel());
3201 AliAODMCParticle *mcDaughterPart =(AliAODMCParticle*)st->UncheckedAt(AssLabel);
3203 Int_t v0lab = mcDaughterPart->GetMother();
3206 if((!v0lab) || (v0lab<0) || (v0lab > st->GetEntriesFast()))continue;//validity checks
3208 AliAODMCParticle *mcp=(AliAODMCParticle*)st->UncheckedAt(v0lab);
3210 Double_t genALaPt = mcp->Pt();
3212 Double_t vFeedDownALa[3] = {5., invMALaFDcand, genALaPt};
3213 fhnFeedDownALa->Fill(vFeedDownALa);
3216 }//end loop over feeddown candidates for Antilambda particles
3218 //_________________________________________________________________
3219 //feeddown for Antilambdas from Xi(bar)+ and Xi(bar)0 in jet cone:
3221 for(Int_t it=0; it<jetConeFDALalist->GetSize(); ++it){
3223 AliAODv0* mcfd = dynamic_cast<AliAODv0*>(jetConeFDALalist->At(it));
3226 Double_t invMALaFDcand = 0;
3227 Double_t trackPt = 0;//pt of ass. particle, not used for the histos
3229 CalculateInvMass(mcfd, kAntiLambda, invMALaFDcand, trackPt);
3231 //Get MC gen. Antilambda transverse momentum
3232 TClonesArray *st = 0x0;
3234 TList *lt = fAOD->GetList();
3237 st = (TClonesArray*)lt->FindObject(AliAODMCParticle::StdBranchName());
3239 AliAODTrack *daughtertrack = (AliAODTrack *) (mcfd->GetSecondaryVtx()->GetDaughter(0));//fetch the first of the two daughter tracks
3240 Int_t AssLabel = TMath::Abs(daughtertrack->GetLabel());
3242 AliAODMCParticle *mcDaughterPart =(AliAODMCParticle*)st->UncheckedAt(AssLabel);
3244 Int_t v0lab = mcDaughterPart->GetMother();
3246 if((!v0lab) || (v0lab<0) || (v0lab > st->GetEntriesFast()))continue;//validity checks
3248 AliAODMCParticle *mcp=(AliAODMCParticle*)st->UncheckedAt(v0lab);
3250 Double_t genALaPt = mcp->Pt();
3252 Double_t vFeedDownALaCone[3] = {jetPt, invMALaFDcand, genALaPt};
3253 fhnFeedDownALaCone->Fill(vFeedDownALaCone);
3256 }//end loop over feeddown candidates for Antilambda particles in jet cone
3260 //____fetch MC generated K0s in cone around jet axis__(note: particles can stem from fragmentation but also from underlying event)________
3262 Double_t sumPtMCgenK0s = 0.;
3263 Bool_t isBadJetMCgenK0s = kFALSE; // dummy, do not use
3266 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!!)
3268 //first: sampling MC gen K0s
3270 GetTracksInCone(fListMCgenK0s, fListMCgenK0sCone, jet, GetFFRadius(), sumPtMCgenK0s, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetMCgenK0s); //MC generated K0s in cone around jet axis
3272 if(fDebug>2)Printf("%s:%d nMCgenK0s in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,fListMCgenK0sCone->GetEntries(),GetFFRadius());
3275 for(Int_t it=0; it<fListMCgenK0sCone->GetSize(); ++it){ // loop MC generated K0s in cone around jet axis
3277 AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenK0sCone->At(it));
3280 //Double_t fRapMCgenK0s = MyRapidity(mcp0->E(),mcp0->Pz());//get rec. particle in cone information
3281 Double_t fEtaMCgenK0s = mcp0->Eta();
3282 Double_t fPtMCgenK0s = mcp0->Pt();
3284 fh2MCgenK0Cone->Fill(jetPt,fPtMCgenK0s);
3285 fh2MCEtagenK0Cone->Fill(jetPt,fEtaMCgenK0s);
3289 //check whether the reconstructed K0s in jet cone are stemming from MC gen K0s (on MCgenK0s list):__________________________________________________
3291 for(Int_t ic=0; ic<jetConeK0list->GetSize(); ++ic){ //loop over all reconstructed K0s in jet cone
3293 //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
3295 Int_t negDaughterpdg;
3296 Int_t posDaughterpdg;
3299 Double_t fPtMCrecK0Match;
3300 Double_t invMK0Match;
3304 Bool_t fPhysicalPrimary = -1;
3305 Int_t MCv0PDGCode =0;
3306 Double_t jetPtSmear = -1;
3308 AliAODv0* v0c = dynamic_cast<AliAODv0*>(jetConeK0list->At(ic));//pointer to reconstructed K0s inside jet cone (cone is placed around reconstructed jet axis)
3310 //AliAODv0* v0c = dynamic_cast<AliAODv0*>(fListK0s->At(ic));//pointer to reconstructed K0s
3313 Bool_t daughtercheck = DaughterTrackCheck(v0c, nnum, pnum);//check daughter tracks have proper sign
3314 if(daughtercheck == kFALSE)continue;
3316 const AliAODTrack *trackMCNeg=(AliAODTrack *)(v0c->GetDaughter(nnum));
3317 const AliAODTrack *trackMCPos=(AliAODTrack *)(v0c->GetDaughter(pnum));
3319 TList *listmc = fAOD->GetList();
3321 Bool_t mclabelcheck = MCLabelCheck(v0c, kK0, trackMCNeg, trackMCPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PDGCode);
3323 if(mclabelcheck == kFALSE)continue;
3324 if(fPhysicalPrimary == kFALSE)continue; //requirements for rec. V0 associated to MC true primary particle
3326 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
3328 //for(Int_t it=0; it<fListMCgenK0sCone->GetSize(); ++it){//belongs to previous definition of rec. eff. of V0s within jet cone
3330 //Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
3331 //AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenK0sCone->At(it));
3332 AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenK0s->At(it));
3335 Bool_t particleMatching = IsParticleMatching(mcp0, v0Label);
3337 if(particleMatching == kFALSE)continue; //if reconstructed V0 particle doesn't match to the associated MC particle go to next stack entry
3338 CalculateInvMass(v0c, kK0, invMK0Match, fPtMCrecK0Match);
3339 Double_t fEta = v0c->Eta();
3340 Double_t fPtMCgenK0s = mcp0->Pt();//pt has to be always MC truth value!
3342 Double_t vMCrecK0Cone[4] = {jetPt, invMK0Match,fPtMCgenK0s,fEta};
3343 fhnMCrecK0Cone->Fill(vMCrecK0Cone); //fill matching rec. K0s in 3D histogram
3345 SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear); //jetPt, cent, jetRadius, ptmintrack, &jetPtSmear
3347 Double_t vMCrecK0ConeSmear[4] = {jetPtSmear, invMK0Match,fPtMCgenK0s,fEta};
3348 fhnMCrecK0ConeSmear->Fill(vMCrecK0ConeSmear);
3350 //fill matching rec. K0s in 3D histogram, jet pT smeared according to deltaptjet distribution width
3353 } // end MCgenK0s / MCgenK0sCone loop
3356 //check the K0s daughters contamination of the jet tracks:
3358 TClonesArray *stackMC = 0x0;
3360 for(Int_t it=0; it<jettracklist->GetSize(); ++it){//loop over all particles in the jet
3362 AliVParticle* trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));//all tracks in jet cone
3363 if(!trackVP)continue;
3364 AliAODTrack *tr = dynamic_cast<AliAODTrack*> (trackVP); //fetch one jet track from the TList
3367 //get MC label information
3368 TList *mclist = fAOD->GetList(); //fetch the MC stack
3369 if(!mclist)continue;
3371 stackMC = (TClonesArray*)mclist->FindObject(AliAODMCParticle::StdBranchName()); //get MCAOD branch in data
3372 if (!stackMC) {Printf("ERROR: stack not available");}
3375 Int_t particleLabel = TMath::Abs(tr->GetLabel()); //fetch jet track label in MC stack
3377 //v0c is pointer to K0s candidate, is fetched already above, here it is just checked again whether daughters are properly ordered by their charge
3379 Bool_t daughterchecks = DaughterTrackCheck(v0c, nnum, pnum);
3381 if(daughterchecks == kFALSE)continue; //make sure that daughters are properly ordered
3383 const AliAODTrack *trackNeg=(AliAODTrack *)(v0c->GetDaughter(nnum)); //fetch v0 daughters of reconstructed K0s
3384 const AliAODTrack *trackPos=(AliAODTrack *)(v0c->GetDaughter(pnum));
3386 if(!trackNeg)continue;
3387 if(!trackPos)continue;
3389 Int_t negAssLabel = TMath::Abs(trackNeg->GetLabel()); //negative (reconstructed) charged track label in MC stack
3390 Int_t posAssLabel = TMath::Abs(trackPos->GetLabel()); //positive (reconstructed) charged track label in MC stack
3393 if(particleLabel == posAssLabel){ //check whether jet track and each of the rec. K0s daughters have same MC stack label -> are identical
3394 AliAODMCParticle* mctrackPos = dynamic_cast<AliAODMCParticle*>(stackMC->At(posAssLabel));
3395 if(!mctrackPos) continue;
3396 Double_t trackPosPt = mctrackPos->Pt();
3397 Double_t trackPosEta = mctrackPos->Eta();
3399 Double_t vK0sSecContinCone[3] = {jetPt, trackPosPt, trackPosEta};
3400 fhnK0sSecContinCone->Fill(vK0sSecContinCone);} //if it's the case, fill jet pt, daughter track pt and track eta in histo
3402 if(particleLabel == negAssLabel){
3403 AliAODMCParticle* mctrackNeg = dynamic_cast<AliAODMCParticle*>(stackMC->At(negAssLabel));
3404 if(!mctrackNeg) continue;
3405 Double_t trackNegPt = mctrackNeg->Pt();
3406 Double_t trackNegEta = mctrackNeg->Eta();
3408 Double_t vK0sSecContinCone[3] = {jetPt, trackNegPt, trackNegEta};
3409 fhnK0sSecContinCone->Fill(vK0sSecContinCone);} //if it's the case, fill jet pt, daughter track pt and track eta in histo
3417 } //end rec-K0-in-cone loop
3419 //________________________________________________________________________________________________________________________________________________________
3421 delete fListMCgenK0sCone;
3426 delete jetConeK0list;
3427 delete jetPerpConeK0list;
3428 delete jetPerpConeLalist;
3429 delete jetPerpConeALalist;
3432 //---------------La--------------------------------------------------------------------------------------------------------------------------------------------
3434 // fQAJetHistosRecCutsLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt() );
3436 for(Int_t it=0; it<fListLa->GetSize(); ++it){ // loop all La
3438 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListLa->At(it));
3443 TVector3 v0MomVect(v0Mom);
3445 Double_t dPhiJetLa = (jet->MomentumVector()->Vect()).DeltaPhi(v0MomVect);
3450 CalculateInvMass(v0, kLambda, invMLa, trackPt); //function to calculate invMass with TLorentzVector class
3451 // Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
3453 //if(incrementJetPt){fh1V0JetPt->Fill(jetPt);}
3455 //fFFHistosIMLaJet->FillFF(trackPt, invMLa, jetPt, incrementJetPt);
3457 if(dPhiJetLa<fh1dPhiJetLa->GetXaxis()->GetXmin()) dPhiJetLa += 2*TMath::Pi();
3458 fh1dPhiJetLa->Fill(dPhiJetLa);
3461 /* if(fListLa->GetSize() == 0){ // no La: increment jet pt spectrum
3463 Bool_t incrementJetPt = kTRUE;
3464 fFFHistosIMLaJet->FillFF(-1, -1, jetPt, incrementJetPt);
3468 // ____fetch rec. Lambdas in cone around jet axis_______________________________________________________________________________________
3470 TList* jetConeLalist = new TList();
3471 Double_t sumPtLa = 0.;
3472 Bool_t isBadJetLa = kFALSE; // dummy, do not use
3474 GetTracksInCone(fListLa, jetConeLalist, jet, GetFFRadius(), sumPtLa, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetLa);//method inherited from FF
3476 if(fDebug>2)Printf("%s:%d nLa total: %d, in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nLa,jetConeLalist->GetEntries(),GetFFRadius());
3478 for(Int_t it=0; it<jetConeLalist->GetSize(); ++it){ // loop La in jet cone
3480 AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetConeLalist->At(it));
3488 CalculateInvMass(v0, kLambda, invMLa, trackPt); //function to calculate invMass with TLorentzVector class
3490 Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;//needed for all histos, which serve for normalisation
3493 Double_t jetPtSmear = -1;
3494 SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
3495 if(incrementJetPt == kTRUE){fh1IMLaConeSmear->Fill(jetPtSmear);} //fill TH1F for normalization purposes
3498 if(incrementJetPt==kTRUE){
3499 fh1IMLaCone->Fill(jetPt);}//normalisation by number of selected jets
3501 //fFFHistosIMLaCone->FillFF(trackPt, invMLa, jetPt, incrementJetPt);
3502 Double_t vLaCone[4] = {jetPt, invMLa,trackPt,fEta};
3503 fhnLaCone->Fill(vLaCone);
3507 if(jetConeLalist->GetSize() == 0){ // no La: increment jet pt spectrum
3509 Bool_t incrementJetPt = kTRUE;
3510 // fFFHistosIMLaCone->FillFF(-1, -1, jetPt, incrementJetPt);
3511 Double_t vLaCone[4] = {jetPt, -1, -1, -1};
3512 fhnLaCone->Fill(vLaCone);
3514 if(incrementJetPt==kTRUE){
3515 fh1IMLaCone->Fill(jetPt);}//normalisation by number of selected jets
3518 Double_t jetPtSmear;
3519 SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
3520 if(incrementJetPt == kTRUE)fh1IMLaConeSmear->Fill(jetPtSmear);}
3526 //____fetch MC generated Lambdas in cone around jet axis__(particles can stem from fragmentation but also from underlying event)_____________
3528 Double_t sumPtMCgenLa = 0.;
3529 Bool_t isBadJetMCgenLa = kFALSE; // dummy, do not use
3531 //sampling MC gen. Lambdas in cone around reconstructed jet axis
3532 fListMCgenLaCone = new TList();
3534 GetTracksInCone(fListMCgenLa, fListMCgenLaCone, jet, GetFFRadius(), sumPtMCgenLa, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetMCgenLa);//fetch MC generated Lambdas in cone of resolution parameter R around jet axis
3536 if(fDebug>2)Printf("%s:%d nMCgenLa in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,fListMCgenLaCone->GetEntries(),GetFFRadius());
3538 for(Int_t it=0; it<fListMCgenLaCone->GetSize(); ++it){ // loop MC generated La in cone around jet axis
3540 AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenLaCone->At(it));
3543 //Double_t fRapMCgenLa = MyRapidity(mcp0->E(),mcp0->Pz());
3544 Double_t fEtaMCgenLa = mcp0->Eta();
3545 Double_t fPtMCgenLa = mcp0->Pt();
3547 fh2MCgenLaCone->Fill(jetPt,fPtMCgenLa);
3548 fh2MCEtagenLaCone->Fill(jetPt,fEtaMCgenLa);
3552 //check whether the reconstructed La are stemming from MC gen La on fListMCgenLa List:__________________________________________________
3554 for(Int_t ic=0; ic<jetConeLalist->GetSize(); ++ic){//loop over all reconstructed La within jet cone, new definition
3556 //for(Int_t ic=0; ic<fListLa->GetSize(); ++ic){//old definition
3558 Int_t negDaughterpdg;
3559 Int_t posDaughterpdg;
3562 Double_t fPtMCrecLaMatch;
3563 Double_t invMLaMatch;
3567 Bool_t fPhysicalPrimary = -1;
3568 Int_t MCv0PDGCode =0;
3569 Double_t jetPtSmear = -1;
3571 AliAODv0* v0c = dynamic_cast<AliAODv0*>(jetConeLalist->At(ic));//new definition
3573 //AliAODv0* v0c = dynamic_cast<AliAODv0*>(fListLa->At(ic));//old definition
3576 Bool_t daughtercheck = DaughterTrackCheck(v0c, nnum, pnum);
3577 if(daughtercheck == kFALSE)continue;
3579 const AliAODTrack *trackMCNeg=(AliAODTrack *)(v0c->GetDaughter(nnum));
3580 const AliAODTrack *trackMCPos=(AliAODTrack *)(v0c->GetDaughter(pnum));
3582 TList *listmc = fAOD->GetList();
3584 Bool_t mclabelcheck = MCLabelCheck(v0c, kLambda, trackMCNeg, trackMCPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PDGCode);
3586 if(mclabelcheck == kFALSE)continue;
3587 if(fPhysicalPrimary == kFALSE)continue;
3589 for(Int_t it=0; it<fListMCgenLa->GetSize(); ++it){//new definition // loop over MC generated K0s in cone around jet axis
3591 // for(Int_t it=0; it<fListMCgenLaCone->GetSize(); ++it){//old definition // loop over MC generated La in cone around jet axis
3593 //Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
3595 AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenLa->At(it));//new definition
3596 //AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenLaCone->At(it));//old definition
3600 Bool_t particleMatching = IsParticleMatching(mcp0, v0Label);
3603 if(particleMatching == kFALSE)continue; //particle doesn't match on any associated MC gen particle in cone around rec jet axis
3605 CalculateInvMass(v0c, kLambda, invMLaMatch, fPtMCrecLaMatch);
3607 Double_t fPtMCgenLa = mcp0->Pt();
3608 Double_t fEta = v0c->Eta();//rec. MC particle
3609 Double_t vMCrecLaCone[4] = {jetPt, invMLaMatch,fPtMCgenLa,fEta};
3610 fhnMCrecLaCone->Fill(vMCrecLaCone);
3612 SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
3614 Double_t vMCrecLaConeSmear[4] = {jetPtSmear, invMLaMatch,fPtMCgenLa,fEta};
3615 fhnMCrecLaConeSmear->Fill(vMCrecLaConeSmear); //fill matching rec. Lambdas in 3D histogram, jet pT smeared according to deltaptjet distribution width
3618 } // end MCgenLa loop
3620 //check the Lambda daughters contamination of the jet tracks://///////////////////////////////////////////////////////////////////////////////////////////
3622 TClonesArray *stackMC = 0x0;
3624 for(Int_t it=0; it<jettracklist->GetSize(); ++it){//loop over all particles in the jet
3626 AliVParticle* trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));//all tracks in jet cone
3627 if(!trackVP)continue;
3628 AliAODTrack *tr = dynamic_cast<AliAODTrack*> (trackVP); //fetch one jet track from the TList
3631 //get MC label information
3632 TList *mclist = fAOD->GetList(); //fetch the MC stack
3634 stackMC = (TClonesArray*)mclist->FindObject(AliAODMCParticle::StdBranchName()); //get MCAOD branch in data
3635 if (!stackMC) {Printf("ERROR: stack not available");}
3638 Int_t particleLabel = TMath::Abs(tr->GetLabel()); //fetch jet track label in MC stack
3640 Bool_t daughterchecks = DaughterTrackCheck(v0c, nnum, pnum);
3642 if(daughterchecks == kFALSE)continue; //make sure that daughters are properly ordered
3644 const AliAODTrack *trackNeg=(AliAODTrack *)(v0c->GetDaughter(nnum)); //fetch v0 daughters of reconstructed K0s
3645 const AliAODTrack *trackPos=(AliAODTrack *)(v0c->GetDaughter(pnum));
3647 Int_t negAssLabel = TMath::Abs(trackNeg->GetLabel()); //negative (reconstructed) charged track label in MC stack
3648 Int_t posAssLabel = TMath::Abs(trackPos->GetLabel()); //positive (reconstructed) charged track label in MC stack
3651 if(particleLabel == posAssLabel){ //check whether jet track and each of the rec. K0s daughters have same MC stack label -> are identical
3653 AliAODMCParticle* mctrackPos = dynamic_cast<AliAODMCParticle*>(stackMC->At(posAssLabel));
3654 if(!mctrackPos) continue;
3656 Double_t trackPosPt = trackPos->Pt();
3657 Double_t trackPosEta = trackPos->Eta();
3658 Double_t vLaSecContinCone[3] = {jetPt, trackPosPt, trackPosEta};
3659 fhnLaSecContinCone->Fill(vLaSecContinCone);
3661 } //if it's the case, fill jet pt, daughter track pt and track eta in histo
3664 if(particleLabel == negAssLabel){
3666 AliAODMCParticle* mctrackNeg = dynamic_cast<AliAODMCParticle*>(stackMC->At(negAssLabel));
3667 if(!mctrackNeg) continue;
3669 Double_t trackNegPt = trackNeg->Pt();
3670 Double_t trackNegEta = trackNeg->Eta();
3672 Double_t vLaSecContinCone[3] = {jetPt, trackNegPt, trackNegEta};
3673 fhnLaSecContinCone->Fill(vLaSecContinCone);
3676 } //if it's the case, fill jet pt, daughter track pt and track eta in histo
3681 } //end rec-La-in-cone loop
3682 //________________________________________________________________________________________________________________________________________________________
3684 delete fListMCgenLaCone;
3688 delete jetConeLalist;
3692 //---------------ALa-----------
3695 // fQAJetHistosRecCutsLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt() );
3697 for(Int_t it=0; it<fListALa->GetSize(); ++it){ // loop all ALa
3699 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListALa->At(it));
3704 TVector3 v0MomVect(v0Mom);
3706 Double_t dPhiJetALa = (jet->MomentumVector()->Vect()).DeltaPhi(v0MomVect);
3708 Double_t invMALa =0;
3711 CalculateInvMass(v0, kAntiLambda, invMALa, trackPt); //function to calculate invMass with TLorentzVector class
3712 //Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
3714 //if(incrementJetPt){fh1V0JetPt->Fill(jetPt);}
3716 //fFFHistosIMALaJet->FillFF(trackPt, invMALa, jetPt, incrementJetPt);
3718 if(dPhiJetALa<fh1dPhiJetALa->GetXaxis()->GetXmin()) dPhiJetALa += 2*TMath::Pi();
3719 fh1dPhiJetALa->Fill(dPhiJetALa);
3722 // if(fListALa->GetSize() == 0){ // no ALa: increment jet pt spectrum
3724 // Bool_t incrementJetPt = kTRUE;
3725 //fFFHistosIMALaJet->FillFF(-1, -1, jetPt, incrementJetPt);
3729 // ____fetch rec. Antilambdas in cone around jet axis_______________________________________________________________________________________
3731 TList* jetConeALalist = new TList();
3732 Double_t sumPtALa = 0.;
3733 Bool_t isBadJetALa = kFALSE; // dummy, do not use
3735 GetTracksInCone(fListALa, jetConeALalist, jet, GetFFRadius(), sumPtALa, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetALa);//method inherited from FF
3737 if(fDebug>2)Printf("%s:%d nALa total: %d, in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nALa,jetConeALalist->GetEntries(),GetFFRadius());
3739 for(Int_t it=0; it<jetConeALalist->GetSize(); ++it){ // loop ALa in jet cone
3741 AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetConeALalist->At(it));
3743 Double_t invMALa =0;
3749 CalculateInvMass(v0, kAntiLambda, invMALa, trackPt); //function to calculate invMass with TLorentzVector class
3751 Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
3753 if(fAnalysisMC){ //jet pt smearing study for Antilambdas
3754 Double_t jetPtSmear = -1;
3755 SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
3756 if(incrementJetPt == kTRUE){fh1IMALaConeSmear->Fill(jetPtSmear);} //fill TH1F for normalization purposes
3759 if(incrementJetPt==kTRUE){
3760 fh1IMALaCone->Fill(jetPt);}//normalisation by number of selected jets
3762 //fFFHistosIMALaCone->FillFF(trackPt, invMALa, jetPt, incrementJetPt);
3763 Double_t vALaCone[4] = {jetPt, invMALa,trackPt,fEta};
3764 fhnALaCone->Fill(vALaCone);
3767 if(jetConeALalist->GetSize() == 0){ // no ALa: increment jet pt spectrum
3769 Bool_t incrementJetPt = kTRUE;
3771 if(incrementJetPt==kTRUE){
3772 fh1IMALaCone->Fill(jetPt);}//normalisation by number of selected jets
3774 //fFFHistosIMALaCone->FillFF(-1, -1, jetPt, incrementJetPt);
3775 Double_t vALaCone[4] = {jetPt, -1, -1, -1};
3776 fhnALaCone->Fill(vALaCone);
3779 Double_t jetPtSmear;
3780 SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
3781 if(incrementJetPt == kTRUE)fh1IMALaConeSmear->Fill(jetPtSmear);}
3787 //____fetch MC generated Antilambdas in cone around jet axis__(particles can stem from fragmentation but also from underlying event)_____________
3789 Double_t sumPtMCgenALa = 0.;
3790 Bool_t isBadJetMCgenALa = kFALSE; // dummy, do not use
3792 //sampling MC gen Antilambdas in cone around reconstructed jet axis
3793 fListMCgenALaCone = new TList();
3795 GetTracksInCone(fListMCgenALa, fListMCgenALaCone, jet, GetFFRadius(), sumPtMCgenALa, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetMCgenALa);//MC generated K0s in cone around jet axis
3797 if(fDebug>2)Printf("%s:%d nMCgenALa in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,fListMCgenALaCone->GetEntries(),GetFFRadius());
3799 for(Int_t it=0; it<fListMCgenALaCone->GetSize(); ++it){ // loop MC generated La in cone around jet axis
3801 AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenALaCone->At(it));
3804 //Double_t fRapMCgenALa = MyRapidity(mcp0->E(),mcp0->Pz());
3805 Double_t fEtaMCgenALa = mcp0->Eta();
3806 Double_t fPtMCgenALa = mcp0->Pt();
3808 fh2MCgenALaCone->Fill(jetPt,fPtMCgenALa);
3809 fh2MCEtagenALaCone->Fill(jetPt,fEtaMCgenALa);
3813 //check whether the reconstructed ALa are stemming from MC gen ALa on MCgenALa List:__________________________________________________
3815 for(Int_t ic=0; ic<jetConeALalist->GetSize(); ++ic){//loop over all reconstructed ALa
3817 Int_t negDaughterpdg;
3818 Int_t posDaughterpdg;
3821 Double_t fPtMCrecALaMatch;
3822 Double_t invMALaMatch;
3826 Bool_t fPhysicalPrimary = -1;
3827 Int_t MCv0PDGCode =0;
3828 Double_t jetPtSmear = -1;
3830 AliAODv0* v0c = dynamic_cast<AliAODv0*>(jetConeALalist->At(ic));
3833 Bool_t daughtercheck = DaughterTrackCheck(v0c, nnum, pnum);
3834 if(daughtercheck == kFALSE)continue;
3836 const AliAODTrack *trackMCNeg=(AliAODTrack *)(v0c->GetDaughter(nnum));
3837 const AliAODTrack *trackMCPos=(AliAODTrack *)(v0c->GetDaughter(pnum));
3839 TList *listmc = fAOD->GetList();
3840 if(!listmc)continue;
3842 Bool_t mclabelcheck = MCLabelCheck(v0c, kAntiLambda, trackMCNeg, trackMCPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PDGCode);
3844 if(mclabelcheck == kFALSE)continue;
3845 if(fPhysicalPrimary == kFALSE)continue;
3847 for(Int_t it=0; it<fListMCgenALa->GetSize(); ++it){ // loop over MC generated Antilambdas in cone around jet axis
3849 //Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
3851 AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenALa->At(it));
3854 Bool_t particleMatching = IsParticleMatching(mcp0, v0Label);
3856 if(particleMatching == kFALSE)continue; //particle doesn't match on any associated MC gen particle in cone around rec jet axis
3858 CalculateInvMass(v0c, kAntiLambda, invMALaMatch, fPtMCrecALaMatch);
3860 Double_t fPtMCgenALa = mcp0->Pt();
3861 Double_t fEta = v0c->Eta();
3862 Double_t vMCrecALaCone[4] = {jetPt, invMALaMatch,fPtMCgenALa,fEta};
3863 fhnMCrecALaCone->Fill(vMCrecALaCone); //fill matching rec. Antilambda in 3D histogram
3865 SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
3867 Double_t vMCrecALaConeSmear[4] = {jetPtSmear, invMALaMatch,fPtMCgenALa,fEta};
3868 fhnMCrecALaConeSmear->Fill(vMCrecALaConeSmear); //fill matching rec. Antilambda in 3D histogram
3870 } // end MCgenALa loop
3874 //check the Antilambda daughters contamination of the jet tracks:
3876 TClonesArray *stackMC = 0x0;
3878 for(Int_t it=0; it<jettracklist->GetSize(); ++it){//loop over all particles in the jet
3880 AliVParticle* trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));//all tracks in jet cone
3881 if(!trackVP)continue;
3882 AliAODTrack *tr = dynamic_cast<AliAODTrack*> (trackVP); //fetch one jet track from the TList
3885 //get MC label information
3886 TList *mclist = fAOD->GetList(); //fetch the MC stack
3887 if(!mclist)continue;
3889 stackMC = (TClonesArray*)mclist->FindObject(AliAODMCParticle::StdBranchName()); //get MCAOD branch in data
3890 if (!stackMC) {Printf("ERROR: stack not available");}
3893 Int_t particleLabel = TMath::Abs(tr->GetLabel()); //fetch jet track label in MC stack
3895 Bool_t daughterchecks = DaughterTrackCheck(v0c, nnum, pnum);
3897 if(daughterchecks == kFALSE)continue; //make sure that daughters are properly ordered
3899 const AliAODTrack *trackNeg=(AliAODTrack *)(v0c->GetDaughter(nnum)); //fetch v0 daughters of reconstructed K0s
3900 const AliAODTrack *trackPos=(AliAODTrack *)(v0c->GetDaughter(pnum));
3901 if(!trackPos)continue;
3902 if(!trackNeg)continue;
3904 Int_t negAssLabel = TMath::Abs(trackNeg->GetLabel()); //negative (reconstructed) charged track label in MC stack
3905 Int_t posAssLabel = TMath::Abs(trackPos->GetLabel()); //positive (reconstructed) charged track label in MC stack
3907 if(!negAssLabel)continue;
3908 if(!posAssLabel)continue;
3910 if(particleLabel == posAssLabel){ //check whether jet track and each of the rec. K0s daughters have same MC stack label -> are identical
3911 AliAODMCParticle* mctrackPos = dynamic_cast<AliAODMCParticle*>(stackMC->At(posAssLabel));
3912 if(!mctrackPos) continue;
3914 Double_t trackPosPt = trackPos->Pt();
3915 Double_t trackPosEta = trackPos->Eta();
3916 if(!trackPosPt)continue;
3917 if(!trackPosEta)continue;
3919 Double_t vLaSecContinCone[3] = {jetPt, trackPosPt, trackPosEta};
3920 fhnLaSecContinCone->Fill(vLaSecContinCone);
3924 //fh3SecContinCone->Fill(jetPt, trackPosPt, trackPosEta);
3925 } //if it's the case, fill jet pt, daughter track pt and track eta in histo
3927 if(particleLabel == negAssLabel){
3929 AliAODMCParticle* mctrackNeg = dynamic_cast<AliAODMCParticle*>(stackMC->At(negAssLabel));
3930 if(!mctrackNeg) continue;
3932 Double_t trackNegPt = trackNeg->Pt();
3933 Double_t trackNegEta = trackNeg->Eta();
3935 if(!trackNegPt)continue;
3936 if(!trackNegEta)continue;
3938 Double_t vLaSecContinCone[3] = {jetPt, trackNegPt, trackNegEta};
3939 fhnLaSecContinCone->Fill(vLaSecContinCone);
3941 //fh3SecContinCone->Fill(jetPt, trackNegPt, trackNegEta);
3942 } //if it's the case, fill jet pt, daughter track pt and track eta in histo
3946 } //end rec-ALa-in-cone loop
3947 //________________________________________________________________________________________________________________________________________________________
3949 delete fListMCgenALaCone;
3953 delete jetConeALalist;
3954 delete jettracklist; //had been initialised at jet loop beginning
3956 }//end of if 'leading' or 'all jet' requirement
3962 fTracksRecCuts->Clear();
3963 fJetsRecCuts->Clear();
3964 fBckgJetsRec->Clear();
3968 fListFeeddownLaCand->Clear();
3969 fListFeeddownALaCand->Clear();
3970 jetConeFDLalist->Clear();
3971 jetConeFDALalist->Clear();
3972 fListMCgenK0s->Clear();
3973 fListMCgenLa->Clear();
3974 fListMCgenALa->Clear();
3978 PostData(1, fCommonHistList);
3983 // ____________________________________________________________________________________________
3984 void AliAnalysisTaskJetChem::SetProperties(TH3F* h,const char* x, const char* y, const char* z)
3986 //Set properties of histos (x,y and z title)
3991 h->GetXaxis()->SetTitleColor(1);
3992 h->GetYaxis()->SetTitleColor(1);
3993 h->GetZaxis()->SetTitleColor(1);
3997 //________________________________________________________________________________________________________________________________________
3998 Bool_t AliAnalysisTaskJetChem::AcceptBetheBloch(AliAODv0 *v0, AliPIDResponse *PIDResponse, const Int_t particletype) //dont use for MC Analysis
4004 const AliAODTrack *ntracktest=(AliAODTrack *)v0->GetDaughter(nnum);
4005 if(ntracktest->Charge() > 0){nnum = 0; pnum = 1;}
4007 const AliAODTrack *trackNeg=(AliAODTrack *)(v0->GetDaughter(nnum));
4008 const AliAODTrack *trackPos=(AliAODTrack *)(v0->GetDaughter(pnum));
4010 //Check if both tracks are available
4011 if (!trackPos || !trackNeg) {
4012 Printf("strange analysis::UserExec:: Error:Could not retrieve one of the daughter tracks\n");
4016 //remove like sign V0s
4017 if ( trackPos->Charge() == trackNeg->Charge() ){
4018 //if(fDebug>1) Printf("%s:%d found like-sign V0", (char*)__FILE__,__LINE__);
4023 Double_t nsig_p = 0; //number of sigmas that positive daughter track has got in TPC pid information
4024 Double_t nsig_n = 0;
4026 const AliAODPid *pid_p=trackPos->GetDetPid(); // returns fDetPID, more detailed or detector specific pid information
4027 const AliAODPid *pid_n=trackNeg->GetDetPid();
4029 if(!pid_p)return kFALSE;
4030 if(!pid_n)return kFALSE;
4034 if(particletype == 1) //PID cut on positive charged Lambda daughters (only those with pt < 1 GeV/c)
4037 nsig_p=PIDResponse->NumberOfSigmasTPC(trackPos,AliPID::kProton);
4038 Double_t protonPt = trackPos->Pt();
4039 if ((TMath::Abs(nsig_p) >= fCutBetheBloch) && (fCutBetheBloch >0) && (protonPt < 1)) return kFALSE;
4048 if(particletype == 2)
4050 nsig_n=PIDResponse->NumberOfSigmasTPC(trackNeg,AliPID::kProton);
4051 Double_t antiprotonPt = trackNeg->Pt();
4052 if ((TMath::Abs(nsig_n) >= fCutBetheBloch) && (fCutBetheBloch >0) && (antiprotonPt < 1)) return kFALSE;
4060 //___________________________________________________________________
4061 Bool_t AliAnalysisTaskJetChem::IsK0InvMass(const Double_t mass) const
4063 // K0 mass ? Use FF histo limits
4065 if(fFFIMInvMMin <= mass && mass < fFFIMInvMMax) return kTRUE;
4069 //___________________________________________________________________
4070 Bool_t AliAnalysisTaskJetChem::IsLaInvMass(const Double_t mass) const
4072 // La mass ? Use FF histo limits
4075 if(fFFIMLaInvMMin <= mass && mass < fFFIMLaInvMMax) return kTRUE;
4080 //_____________________________________________________________________________________
4081 Int_t AliAnalysisTaskJetChem::GetListOfV0s(TList *list, const Int_t type, const Int_t particletype, AliAODVertex* primVertex, AliAODEvent* aod)
4083 // fill list of V0s selected according to type
4086 if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
4091 if(fDebug>5){std::cout<<"AliAnalysisTaskJetChem::GetListOfV0s(): type: "<<type<<" particletype: "<<particletype<<"aod: "<<aod<<std::endl;
4092 if(type==kTrackUndef){std::cout<<"AliAnalysisTaskJetChem::GetListOfV0s(): kTrackUndef!! "<<std::endl;}
4096 if(type==kTrackUndef) return 0;
4098 if(!primVertex) return 0;
4100 Double_t lPrimaryVtxPosition[3];
4101 Double_t lV0Position[3];
4102 lPrimaryVtxPosition[0] = primVertex->GetX();
4103 lPrimaryVtxPosition[1] = primVertex->GetY();
4104 lPrimaryVtxPosition[2] = primVertex->GetZ();
4106 if(fDebug>5){ std::cout<<"AliAnalysisTaskJetChem::GetListOfV0s(): aod->GetNumberOfV0s: "<<aod->GetNumberOfV0s()<<std::endl; }
4109 for(int i=0; i<aod->GetNumberOfV0s(); i++){ // loop over V0s
4112 AliAODv0* v0 = aod->GetV0(i);
4116 std::cout << std::endl
4117 << "Warning in AliAnalysisTaskJetChem::GetListOfV0s:" << std::endl
4118 << "v0 = " << v0 << std::endl;
4122 Bool_t isOnFly = v0->GetOnFlyStatus();
4124 if(!isOnFly && (type == kOnFly || type == kOnFlyPID || type == kOnFlydEdx || type == kOnFlyPrim)) continue;
4125 if( isOnFly && (type == kOffl || type == kOfflPID || type == kOffldEdx || type == kOfflPrim)) continue;
4127 Int_t motherType = -1;
4128 //Double_t v0CalcMass = 0; //mass of MC v0
4129 Double_t MCPt = 0; //pt of MC v0
4131 Double_t pp[3]={0,0,0}; //3-momentum positive charged track
4132 Double_t pm[3]={0,0,0}; //3-momentum negative charged track
4133 Double_t v0mom[3]={0,0,0};
4144 Bool_t daughtercheck = DaughterTrackCheck(v0, nnum, pnum);
4146 if(daughtercheck == kFALSE)continue;
4148 const AliAODTrack *trackNeg=(AliAODTrack *)(v0->GetDaughter(nnum));
4149 const AliAODTrack *trackPos=(AliAODTrack *)(v0->GetDaughter(pnum));
4152 ///////////////////////////////////////////////////////////////////////////////////
4154 //calculate InvMass for every V0 particle assumption (Kaon=1,Lambda=2,Antilambda=3)
4155 switch(particletype){
4157 CalculateInvMass(v0, kK0, invM, trackPt); //function to calculate invMass with TLorentzVector class
4161 CalculateInvMass(v0, kLambda, invM, trackPt);
4165 CalculateInvMass(v0, kAntiLambda, invM, trackPt);
4169 std::cout<<"***NO VALID PARTICLETYPE***"<<std::endl;
4174 /////////////////////////////////////////////////////////////
4175 //V0 and track Cuts:
4178 if(fDebug>7){if(!(IsK0InvMass(invMK0s)) && !(IsLaInvMass(invMLa)) && !(IsLaInvMass(invMALa))){std::cout<<"AliAnalysisTaskJetChem::GetListOfV0s: invM not in selected mass window "<<std::endl;}}
4180 if(!(IsK0InvMass(invMK0s)) && !(IsLaInvMass(invMLa)) && !(IsLaInvMass(invMALa)))continue;
4182 // Double_t PosEta = trackPos->AliAODTrack::Eta();//daughter track charge is sometimes wrong here, account for that!!!
4183 // Double_t NegEta = trackNeg->AliAODTrack::Eta();
4185 Double_t PosEta = trackPos->Eta();//daughter track charge is sometimes wrong here, account for that!!!
4186 Double_t NegEta = trackNeg->Eta();
4188 Double_t PosCharge = trackPos->Charge();
4189 Double_t NegCharge = trackNeg->Charge();
4191 if((trackPos->Charge() == 1) && (trackNeg->Charge() == -1)) //Fill daughters charge into histo to check if they are symmetric distributed
4192 { fh1PosDaughterCharge->Fill(PosCharge);
4193 fh1NegDaughterCharge->Fill(NegCharge);
4196 //DistOverTotMom_in_2D___________
4198 Float_t fMassK0s = TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
4199 Float_t fMassLambda = TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
4202 AliAODVertex* primVtx = fAOD->GetPrimaryVertex(); // get the primary vertex
4203 Double_t dPrimVtxPos[3]; // primary vertex position {x,y,z}
4204 primVtx->GetXYZ(dPrimVtxPos);
4206 Float_t fPtV0 = TMath::Sqrt(v0->Pt2V0()); // transverse momentum of V0
4207 Double_t dSecVtxPos[3]; // V0 vertex position {x,y,z}
4208 v0->GetSecondaryVtx(dSecVtxPos);
4209 Double_t dDecayPath[3];
4210 for (Int_t iPos = 0; iPos < 3; iPos++)
4211 dDecayPath[iPos] = dSecVtxPos[iPos]-dPrimVtxPos[iPos]; // vector of the V0 path
4212 Float_t fDecLen2D = TMath::Sqrt(dDecayPath[0]*dDecayPath[0]+dDecayPath[1]*dDecayPath[1]); //transverse path length R
4213 Float_t fROverPt = fDecLen2D/fPtV0; // R/pT
4215 Float_t fMROverPtK0s = fMassK0s*fROverPt; // m*R/pT
4216 Float_t fMROverPtLambda = fMassLambda*fROverPt; // m*R/pT
4218 //___________________
4219 Double_t fRap = -999;//init values
4220 Double_t fEta = -999;
4221 Double_t fV0cosPointAngle = -999;
4222 Double_t fV0DecayLength = v0->DecayLengthV0(lPrimaryVtxPosition);
4226 fV0mom[0]=v0->MomV0X();
4227 fV0mom[1]=v0->MomV0Y();
4228 fV0mom[2]=v0->MomV0Z();
4230 Double_t fV0TotalMomentum = TMath::Sqrt(fV0mom[0]*fV0mom[0]+fV0mom[1]*fV0mom[1]+fV0mom[2]*fV0mom[2]);
4231 // const Double_t K0sPDGmass = 0.497614;
4232 // const Double_t LambdaPDGmass = 1.115683;
4234 const Double_t K0sPDGmass = TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
4235 const Double_t LambdaPDGmass = TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
4237 Double_t fDistOverTotMomK0s = 0;
4238 Double_t fDistOverTotMomLa = 0;
4240 //calculate proper lifetime of particles in 3D (not recommended anymore)
4242 if(particletype == kK0){
4244 fDistOverTotMomK0s = fV0DecayLength * K0sPDGmass;
4245 fDistOverTotMomK0s /= (fV0TotalMomentum+1e-10);
4248 if((particletype == kLambda)||(particletype == kAntiLambda)){
4250 fDistOverTotMomLa = fV0DecayLength * LambdaPDGmass;
4251 fDistOverTotMomLa /= (fV0TotalMomentum+1e-10);
4254 //TPC cluster (not used anymore) and TPCRefit cuts
4256 //Double_t trackPosNcls = trackPos->GetTPCNcls();//Get number of clusters for positive charged tracks
4257 //Double_t trackNegNcls = trackNeg->GetTPCNcls();//Get number of clusters for negative charged tracks
4259 if(fRequireTPCRefit==kTRUE){//if kTRUE: accept only if daughter track is refitted in TPC!!
4260 Bool_t isPosTPCRefit = (trackPos->AliAODTrack::IsOn(AliESDtrack::kTPCrefit));
4261 Bool_t isNegTPCRefit = (trackNeg->AliAODTrack::IsOn(AliESDtrack::kTPCrefit));
4262 if (!isPosTPCRefit)continue;
4263 if (!isNegTPCRefit)continue;
4266 if(fKinkDaughters==kFALSE){//if kFALSE: no acceptance of kink daughters
4267 AliAODVertex* ProdVtxDaughtersPos = (AliAODVertex*) (trackPos->AliAODTrack::GetProdVertex());
4268 Char_t isAcceptKinkDaughtersPos = ProdVtxDaughtersPos->GetType();
4269 if(isAcceptKinkDaughtersPos==AliAODVertex::kKink)continue;
4271 AliAODVertex* ProdVtxDaughtersNeg = (AliAODVertex*) (trackNeg->AliAODTrack::GetProdVertex());
4272 Char_t isAcceptKinkDaughtersNeg = ProdVtxDaughtersNeg->GetType();
4273 if(isAcceptKinkDaughtersNeg==AliAODVertex::kKink)continue;
4277 Double_t fV0Radius = -999;
4278 Double_t fDcaV0Daughters = v0->DcaV0Daughters();
4279 Double_t fDcaPosToPrimVertex = v0->DcaPosToPrimVertex();//IP of positive charged daughter
4280 Double_t fDcaNegToPrimVertex = v0->DcaNegToPrimVertex();//IP of negative charged daughter
4281 Double_t avDecayLengthK0s = 2.6844;
4282 Double_t avDecayLengthLa = 7.89;
4284 //Float_t fCTauK0s = 2.6844; // [cm] c tau of K0S
4285 //Float_t fCTauLambda = 7.89; // [cm] c tau of Lambda and Antilambda
4287 fV0cosPointAngle = v0->CosPointingAngle(lPrimaryVtxPosition);
4288 lV0Position[0]= v0->DecayVertexV0X();
4289 lV0Position[1]= v0->DecayVertexV0Y();
4290 lV0Position[2]= v0->DecayVertexV0Z();
4292 fV0Radius = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
4294 if(particletype == kK0) {fRap = v0->RapK0Short();
4295 fEta = v0->PseudoRapV0();}
4296 if(particletype == kLambda) {fRap = v0->RapLambda();
4297 fEta = v0->PseudoRapV0();}
4298 if(particletype == kAntiLambda) {fRap = v0->Y(-3122);
4299 fEta = v0->PseudoRapV0();}
4302 //cut on 3D DistOverTotMom: (not used anymore)
4303 //if((particletype == kLambda)||(particletype == kAntiLambda)){if(fDistOverTotMomLa >= (fCutV0DecayMax * avDecayLengthLa)) continue;}
4305 //cut on K0s applied below after all other cuts for histo fill purposes..
4307 //cut on 2D DistOverTransMom: (recommended from Iouri)
4308 if((particletype == kLambda)||(particletype == kAntiLambda)){if(fMROverPtLambda > (fCutV0DecayMax * avDecayLengthLa))continue;}//fCutV0DecayMax set to 5 in AddTask macro
4310 //Armenteros Podolanski Plot for K0s:////////////////////////////
4312 Double_t ArmenterosAlpha=-999;
4313 Double_t ArmenterosPt=-999;
4319 if(particletype == kK0){
4321 pp[0]=v0->MomPosX();
4322 pp[1]=v0->MomPosY();
4323 pp[2]=v0->MomPosZ();
4325 pm[0]=v0->MomNegX();
4326 pm[1]=v0->MomNegY();
4327 pm[2]=v0->MomNegZ();
4330 v0mom[0]=v0->MomV0X();
4331 v0mom[1]=v0->MomV0Y();
4332 v0mom[2]=v0->MomV0Z();
4334 TVector3 v0Pos(pp[0],pp[1],pp[2]);
4335 TVector3 v0Neg(pm[0],pm[1],pm[2]);
4336 TVector3 v0totMom(v0mom[0], v0mom[1], v0mom[2]); //vector for tot v0 momentum
4338 PosPt = v0Pos.Perp(v0totMom); //longitudinal momentum of positive charged daughter track
4339 PosPl = v0Pos.Dot(v0totMom)/v0totMom.Mag(); //transversal momentum of positive charged daughter track
4341 NegPt = v0Neg.Perp(v0totMom); //longitudinal momentum of negative charged daughter track
4342 NegPl = v0Neg.Dot(v0totMom)/v0totMom.Mag(); //transversal momentum of nergative charged daughter track
4344 ArmenterosAlpha = 1.-2./(1+(PosPl/NegPl));
4345 ArmenterosPt= v0->PtArmV0();
4349 if(particletype == kK0){//only cut on K0s histos
4350 if(IsArmenterosSelected == 1){// Armenteros Cut to reject Lambdas contamination in K0s inv. massspectrum
4351 fh2ArmenterosBeforeCuts->Fill(ArmenterosAlpha,ArmenterosPt);
4355 //some more cuts on v0s and daughter tracks:
4358 if((TMath::Abs(PosEta)>fCutPostrackEta) || (TMath::Abs(NegEta)>fCutNegtrackEta))continue; //Daughters pseudorapidity cut
4359 if (fV0cosPointAngle < fCutV0cosPointAngle) continue; //cospointangle cut
4361 //if(TMath::Abs(fRap) > fCutRap)continue; //V0 Rapidity Cut
4362 if(TMath::Abs(fEta) > fCutEta) continue; //V0 Eta Cut
4363 if (fDcaV0Daughters > fCutDcaV0Daughters)continue;
4364 if ((fDcaPosToPrimVertex < fCutDcaPosToPrimVertex) || (fDcaNegToPrimVertex < fCutDcaNegToPrimVertex))continue;
4365 if ((fV0Radius < fCutV0RadiusMin) || (fV0Radius > fCutV0RadiusMax))continue;
4367 const AliAODPid *pid_p1=trackPos->GetDetPid();
4368 const AliAODPid *pid_n1=trackNeg->GetDetPid();
4371 if(particletype == kLambda){
4372 // if(AcceptBetheBloch(v0, fPIDResponse, 1) == kFALSE){std::cout<<"******PID cut rejects Lambda!!!************"<<std::endl;}
4373 if(AcceptBetheBloch(v0, fPIDResponse, 1) == kFALSE)continue;
4374 fh2BBLaPos->Fill(pid_p1->GetTPCmomentum(),pid_p1->GetTPCsignal());//positive lambda daughter
4375 fh2BBLaNeg->Fill(pid_n1->GetTPCmomentum(),pid_n1->GetTPCsignal());//negative lambda daughter
4377 //Double_t phi = v0->Phi();
4378 //Double_t massLa = v0->MassLambda();
4380 //printf("La: i = %d, m = %f, pT = %f, eta = %f, phi = %f\n, ",i,massLa,trackPt,fEta,phi);
4384 if(particletype == kAntiLambda){
4386 if(AcceptBetheBloch(v0, fPIDResponse, 2) == kFALSE)continue;
4387 fh2BBLaPos->Fill(pid_p1->GetTPCmomentum(),pid_p1->GetTPCsignal());//positive antilambda daughter
4388 fh2BBLaNeg->Fill(pid_n1->GetTPCmomentum(),pid_n1->GetTPCsignal());//negative antilambda daughter
4393 //Armenteros cut on K0s:
4394 if(particletype == kK0){
4395 if(IsArmenterosSelected == 1){// Armenteros Cut to reject Lambdas contamination in K0s inv. massspectrum
4397 if((ArmenterosPt<=(TMath::Abs(fCutArmenteros*ArmenterosAlpha))) && (fCutArmenteros!=-999))continue; //Cuts out Lambda contamination in K0s histos
4398 fh2ArmenterosAfterCuts->Fill(ArmenterosAlpha,ArmenterosPt);
4402 //not used anymore in 3D, z component of total momentum has bad resolution, cut instead in 2D and use pT
4403 //Proper Lifetime Cut: DecayLength3D * PDGmass / |p_tot| < 3*2.68cm (ctau(betagamma=1)) ; |p|/mass = beta*gamma
4404 //////////////////////////////////////////////
4407 //cut on 2D DistOverTransMom
4408 if(particletype == kK0){//the cut on Lambdas you can find above
4410 fh2ProperLifetimeK0sVsPtBeforeCut->Fill(trackPt,fMROverPtK0s); //fill these histos after all other cuts
4411 if(fMROverPtK0s > (fCutV0DecayMax * avDecayLengthK0s))continue;
4412 fh2ProperLifetimeK0sVsPtAfterCut->Fill(trackPt,fMROverPtK0s);
4414 //Double_t phi = v0->Phi();
4415 // Double_t massK0s = v0->MassK0Short();
4416 //printf("K0S: i = %d, m = %f, pT = %f, eta = %f, phi = %f\n",i,invMK0s,trackPt,fEta,phi);
4418 //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;
4421 //MC Associated V0 particles: (reconstructed particles associated with MC truth (MC truth: true primary MC generated particle))
4424 if(fAnalysisMC){// begin MC part
4426 Int_t negDaughterpdg = 0;
4427 Int_t posDaughterpdg = 0;
4429 Bool_t fPhysicalPrimary = -1; //v0 physical primary check
4430 Int_t MCv0PdgCode = 0;
4431 Bool_t mclabelcheck = kFALSE;
4433 TList *listmc = aod->GetList(); //AliAODEvent* is inherited from AliVEvent*, listmc is pointer to reconstructed event in MC list, member of AliAODEvent
4435 if(!listmc)continue;
4437 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
4439 //feeddown-correction for Lambda/Antilambda particles
4440 //feedddown comes mainly from charged and neutral Xi particles
4441 //feeddown from Sigma decays so quickly that it's not possible to distinguish from primary Lambdas with detector
4442 //feeddown for K0s from phi decays is neglectible
4443 //TH2F* fh2FeedDownMatrix = 0x0; //histo for feeddown already decleared above
4446 //first for all Lambda and Antilambda candidates____________________________________________________________________
4448 if(particletype == kLambda){
4450 mclabelcheck = MCLabelCheck(v0, kLambda, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode);
4453 if((motherType == 3312)||(motherType == 3322)){//mother of v0 is neutral or negative Xi
4455 fListFeeddownLaCand->Add(v0); //fill TList with ass. particles, stemming from feeddown from Xi(bar) decays
4459 if(particletype == kAntiLambda){
4460 mclabelcheck = MCLabelCheck(v0, kAntiLambda, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode);
4462 if((motherType == -3312)||(motherType == -3322)){
4463 fListFeeddownALaCand->Add(v0); //fill TList with ass. particles, stemming from feeddown from Xi(bar) decays
4468 //_only true primary particles survive the following checks_______________________________________________________________________________________________
4470 if(particletype == kK0){
4471 mclabelcheck = MCLabelCheck(v0, kK0, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode);
4472 if(mclabelcheck == kFALSE)continue;
4474 if(particletype == kLambda){
4475 mclabelcheck = MCLabelCheck(v0, kLambda, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode);
4476 if(mclabelcheck == kFALSE)continue;
4478 if(particletype == kAntiLambda){
4479 mclabelcheck = MCLabelCheck(v0, kAntiLambda, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode);
4480 if(mclabelcheck == kFALSE)continue;
4483 if(fPhysicalPrimary != 1)continue; //V0 candidate (K0s, Lambda or Antilambda) must be physical primary, this means there is no mother particle existing
4492 Int_t nPart=list->GetSize();
4495 } // end GetListOfV0s()
4497 // -------------------------------------------------------------------------------------------------------
4499 void AliAnalysisTaskJetChem::CalculateInvMass(AliAODv0* v0vtx, const Int_t particletype, Double_t& invM, Double_t& trackPt){
4509 Double_t pp[3]={0,0,0}; //3-momentum positive charged track
4510 Double_t pm[3]={0,0,0}; //3-momentum negative charged track
4512 const Double_t massPi = 0.13957018; //better use PDG code at this point
4513 const Double_t massP = 0.93827203;
4518 TLorentzVector vector; //lorentzvector V0 particle
4519 TLorentzVector fourmom1;//lorentzvector positive daughter
4520 TLorentzVector fourmom2;//lorentzvector negative daughter
4522 //--------------------------------------------------------------
4524 AliAODTrack *trackPos = (AliAODTrack *) (v0vtx->GetSecondaryVtx()->GetDaughter(0));//index 0 defined as positive charged track in AliESDFilter
4526 if( trackPos->Charge() == 1 ){
4528 pp[0]=v0vtx->MomPosX();
4529 pp[1]=v0vtx->MomPosY();
4530 pp[2]=v0vtx->MomPosZ();
4532 pm[0]=v0vtx->MomNegX();
4533 pm[1]=v0vtx->MomNegY();
4534 pm[2]=v0vtx->MomNegZ();
4537 if( trackPos->Charge() == -1 ){
4539 pm[0]=v0vtx->MomPosX();
4540 pm[1]=v0vtx->MomPosY();
4541 pm[2]=v0vtx->MomPosZ();
4543 pp[0]=v0vtx->MomNegX();
4544 pp[1]=v0vtx->MomNegY();
4545 pp[2]=v0vtx->MomNegZ();
4548 if (particletype == kK0){ // case K0s
4549 mass1 = massPi;//positive particle
4550 mass2 = massPi;//negative particle
4551 } else if (particletype == kLambda){ // case Lambda
4552 mass1 = massP;//positive particle
4553 mass2 = massPi;//negative particle
4554 } else if (particletype == kAntiLambda){ //case AntiLambda
4555 mass1 = massPi;//positive particle
4556 mass2 = massP; //negative particle
4559 fourmom1.SetXYZM(pp[0],pp[1],pp[2],mass1);//positive track
4560 fourmom2.SetXYZM(pm[0],pm[1],pm[2],mass2);//negative track
4561 vector=fourmom1 + fourmom2;
4564 trackPt = vector.Pt();
4566 /*// 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
4568 if(particletype == kK0){
4569 std::cout << "invMK0s: " << invM <<std::endl;
4570 std::cout << "v0vtx->MassK0Short(): " << v0vtx->MassK0Short() << std::endl;
4571 std::cout << " " <<std::endl;
4572 //invM = v0vtx->MassK0Short();
4575 if(particletype == kLambda){
4576 std::cout << "invMLambda: " << invM <<std::endl;
4577 std::cout << "v0vtx->MassMassLambda(): " << v0vtx->MassLambda() << std::endl;
4578 std::cout << " " <<std::endl;
4579 //invM = v0vtx->MassLambda();
4582 if(particletype == kAntiLambda){
4583 std::cout << "invMAntiLambda: " << invM <<std::endl;
4584 std::cout << "v0vtx->MassAntiLambda(): " << v0vtx->MassAntiLambda() << std::endl;
4585 std::cout << " " <<std::endl;
4586 //invM = v0vtx->MassAntiLambda();
4594 //_____________________________________________________________________________________
4595 Int_t AliAnalysisTaskJetChem::GetListOfMCParticles(TList *outputlist, const Int_t particletype, AliAODEvent *mcaodevent) //(list to fill here e.g. fListMCgenK0s, particle species to search for)
4598 outputlist->Clear();
4600 TClonesArray *stack = 0x0;
4601 Double_t mcXv=0., mcYv=0., mcZv=0.;//MC vertex position
4604 // get MC generated particles
4606 Int_t fPdgcodeCurrentPart = 0; //pdg code current particle
4607 Double_t fRapCurrentPart = 0; //get rapidity
4608 Double_t fPtCurrentPart = 0; //get transverse momentum
4609 Double_t fEtaCurrentPart = 0; //get pseudorapidity
4611 //variable for check: physical primary particle
4612 //Bool_t IsPhysicalPrimary = -1;
4613 //Int_t index = 0; //check number of injected particles
4614 //****************************
4615 // Start loop over MC particles
4617 TList *lst = mcaodevent->GetList();
4620 if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
4624 stack = (TClonesArray*)lst->FindObject(AliAODMCParticle::StdBranchName());
4626 Printf("ERROR: stack not available");
4630 AliAODMCHeader *mcHdr=(AliAODMCHeader*)lst->FindObject(AliAODMCHeader::StdBranchName());
4631 if(!mcHdr)return -1;
4633 mcXv=mcHdr->GetVtxX(); mcYv=mcHdr->GetVtxY(); mcZv=mcHdr->GetVtxZ(); // position of the MC primary vertex
4636 ntrk=stack->GetEntriesFast();
4638 //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...
4641 for (Int_t iMc = 0; iMc < ntrk; iMc++) { //loop over mc generated particles
4644 AliAODMCParticle *p0=(AliAODMCParticle*)stack->UncheckedAt(iMc);
4646 //Printf("ERROR: particle with label %d not found in stack (mc loop)", iMc);
4649 fPdgcodeCurrentPart = p0->GetPdgCode();
4651 // Keep only K0s, Lambda and AntiLambda, Xi and Phi:
4652 //if ( (fPdgcodeCurrentPart != 310 ) && (fPdgcodeCurrentPart != 3122 ) && (fPdgcodeCurrentPart != -3122 ) && (fPdgcodeCurrentPart != 3312 ) && (fPdgcodeCurrentPart != -3312) && (fPdgcodeCurrentPart != -333) ) continue;
4656 //Rejection of Pythia injected particles with David Chinellatos method - not the latest method, better Method with TString from MC generator in IsInjected() function below!
4658 /* if( (p0->GetStatus()==21) ||
4659 ((p0->GetPdgCode() == 443) &&
4660 (p0->GetMother() == -1) &&
4661 (p0->GetDaughter(0) == (iMc))) ){ index++; }
4663 if(p0->GetStatus()==21){std::cout<< "hello !!!!" <<std::endl;}
4665 std::cout<< "MC particle status: " << p0->GetStatus() <<std::endl;
4669 //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
4672 //injected particles could be from GenBox (single high pt tracks) or jet related tracks, both generated from PYTHIA MC generator
4674 //Check: MC particle mother
4676 //for feed-down checks
4677 /* //MC gen particles
4678 Int_t iMother = p0->GetMother(); //Motherparticle of V0 candidate (e.g. phi particle,..)
4680 AliAODMCParticle *partM = (AliAODMCParticle*)stack->UncheckedAt(iMother);
4682 if(partM) codeM = TMath::Abs(partM->GetPdgCode());
4685 3312 Xi- -3312 Xibar+
4686 3322 Xi0 -3322 Xibar0
4689 if((codeM == 3312)||(codeM == 3322))// feeddown for Lambda coming from Xi- and Xi0
4695 /* //Check: MC gen. particle decays via 2-pion decay? -> only to be done for the rec. particles !! (-> branching ratio ~ 70 % for K0s -> pi+ pi-)
4697 Int_t daughter0Label = p0->GetDaughter(0);
4698 AliAODMCParticle *mcDaughter0 = (AliAODMCParticle *)stack->UncheckedAt(daughter0Label);
4699 if(daughter0Label >= 0)
4700 {daughter0Type = mcDaughter0->GetPdgCode();}
4702 Int_t daughter1Label = p0->GetDaughter(1);
4703 AliAODMCParticle *mcDaughter1 = (AliAODMCParticle *)stack->UncheckedAt(daughter1Label);
4705 if(daughter1Label >= 1)
4706 {daughter1Type = mcDaughter1->GetPdgCode();} //requirement that daughters are pions is only done for the reconstructed V0s in GetListofV0s() below
4711 // Keep only K0s, Lambda and AntiLambda:
4712 if ( (fPdgcodeCurrentPart != 310 ) && (fPdgcodeCurrentPart != 3122 ) && (fPdgcodeCurrentPart != -3122 )) continue;
4713 // Check: Is physical primary
4715 //do not use anymore: //IsPhysicalPrimary = p0->IsPhysicalPrimary();
4716 //if(!IsPhysicalPrimary)continue;
4718 Float_t fDistPrimaryMax = 0.01; // [cm] max distance of production point to the primary vertex (criterion for choice of MC particles considered as primary)
4720 // Get the distance between production point of the MC mother particle and the primary vertex
4722 Double_t dx = mcXv-p0->Xv();//mc primary vertex - mc gen. v0 vertex
4723 Double_t dy = mcYv-p0->Yv();
4724 Double_t dz = mcZv-p0->Zv();
4726 Double_t fDistPrimary = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
4727 Bool_t fPhysicalPrimary = (fDistPrimary < fDistPrimaryMax);
4729 if(!fPhysicalPrimary)continue;
4731 //if(fPhysicalPrimary){std::cout<<"hello**********************"<<std::endl;}
4733 /* std::cout<<"dx: "<<dx<<std::endl;
4734 std::cout<<"dy: "<<dy<<std::endl;
4735 std::cout<<"dz: "<<dz<<std::endl;
4737 std::cout<<"start: "<<std::endl;
4738 std::cout<<"mcXv: "<<mcXv<<std::endl;
4739 std::cout<<"mcYv: "<<mcYv<<std::endl;
4740 std::cout<<"mcZv: "<<mcZv<<std::endl;
4742 std::cout<<"p0->Xv(): "<<p0->Xv()<<std::endl;
4743 std::cout<<"p0->Yv(): "<<p0->Yv()<<std::endl;
4744 std::cout<<"p0->Zv(): "<<p0->Zv()<<std::endl;
4745 std::cout<<" "<<std::endl;
4746 std::cout<<"fDistPrimary"<<fDistPrimary<<std::endl;
4747 std::cout<<"fPhysicalPrimary"<<fPhysicalPrimary<<std::endl;
4750 //Is close enough to primary vertex to be considered as primary-like?
4752 fRapCurrentPart = MyRapidity(p0->E(),p0->Pz());
4753 fEtaCurrentPart = p0->Eta();
4754 fPtCurrentPart = p0->Pt();
4756 if (TMath::Abs(fEtaCurrentPart) < fCutEta){
4757 // if (TMath::Abs(fRapCurrentPart) > fCutRap)continue; //rap cut for crosschecks
4759 if(particletype == kK0){ //MC gen. K0s
4760 if (fPdgcodeCurrentPart==310){
4761 outputlist->Add(p0);
4765 if(particletype == kLambda){ //MC gen. Lambdas
4766 if (fPdgcodeCurrentPart==3122) {
4767 outputlist->Add(p0);
4771 if(particletype == kAntiLambda){
4772 if (fPdgcodeCurrentPart==-3122) { //MC gen. Antilambdas
4773 outputlist->Add(p0);
4778 }//end loop over MC generated particle
4780 Int_t nMCPart=outputlist->GetSize();
4787 //---------------------------------------------------------------------------
4789 Bool_t AliAnalysisTaskJetChem::FillFeeddownMatrix(TList* fListFeeddownCand, Int_t particletype)
4792 // Define Feeddown matrix
4793 Double_t lFeedDownMatrix [100][100];
4794 // FeedDownMatrix [Lambda Bin][Xi Bin];
4796 //Initialize entries of matrix:
4797 for(Int_t ilb = 0; ilb<100; ilb++){
4798 for(Int_t ixb = 0; ixb<100; ixb++){
4799 lFeedDownMatrix[ilb][ixb]=0; //first lambda bins, xi bins
4804 //----------------------------------------------------------------------------
4806 Double_t AliAnalysisTaskJetChem::MyRapidity(Double_t rE, Double_t rPz) const
4808 // Local calculation for rapidity
4809 return 0.5*TMath::Log((rE+rPz)/(rE-rPz+1.e-13));
4811 //----------------------------------------------------------------------------
4813 // ________________________________________________________________________________________________________________________//function to get the MC gen. jet particles
4815 void AliAnalysisTaskJetChem::GetTracksInCone(TList* inputlist, TList* outputlist, const AliAODJet* jet,
4816 const Double_t radius, Double_t& sumPt, const Double_t minPt, const Double_t maxPt, Bool_t& isBadPt)
4818 // fill list of tracks in cone around jet axis
4821 Bool_t isBadMaxPt = kFALSE;
4822 Bool_t isBadMinPt = kTRUE;
4826 jet->PxPyPz(jetMom);
4827 TVector3 jet3mom(jetMom);
4829 //if(jetets < jetetscutr)continue;
4831 for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){
4833 AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack));
4835 Double_t trackMom[3];
4836 track->PxPyPz(trackMom);
4837 TVector3 track3mom(trackMom);
4839 Double_t dR = jet3mom.DeltaR(track3mom);
4843 outputlist->Add(track);
4845 sumPt += track->Pt();
4847 if(maxPt>0 && track->Pt()>maxPt) isBadMaxPt = kTRUE; // reject jets containing any track with pt larger than this value, use GetFFMaxTrackPt()
4848 if(minPt>0 && track->Pt()>minPt) isBadMinPt = kFALSE; // reject jets with leading track with pt smaller than this value, use GetFFMinLTrackPt()
4854 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)
4855 if(maxPt>0 && isBadMaxPt) isBadPt = kTRUE; //..or because of leading track with too high pt (could be fake track)
4861 //____________________________________________________________________________________________________________________
4864 void AliAnalysisTaskJetChem::GetTracksInPerpCone(TList* inputlist, TList* outputlist, const AliAODJet* jet,
4865 const Double_t radius, Double_t& sumPerpPt)
4867 // fill list of tracks in two cones around jet axis rotated in phi +/- 90 degrees
4869 Double_t jetMom[3]; //array for entries in TVector3
4870 Double_t perpjetplusMom[3]; //array for entries in TVector3
4871 Double_t perpjetnegMom[3];
4875 jet->PxPyPz(jetMom); //get 3D jet momentum
4876 Double_t jetPerpPt = jet->Pt(); //original jet pt, invariant under rotations
4877 Double_t jetPhi = jet->Phi(); //original jet phi
4879 Double_t jetPerpposPhi = jetPhi + ((TMath::Pi())*0.5);//get new perp. jet axis phi clockwise
4880 Double_t jetPerpnegPhi = jetPhi - ((TMath::Pi())*0.5);//get new perp. jet axis phi counterclockwise
4882 TVector3 jet3mom(jetMom); //3-Vector for original rec. jet axis
4884 //Double_t phitest = jet3mom.Phi();
4886 perpjetplusMom[0]=(TMath::Cos(jetPerpposPhi)*jetPerpPt); //x coordinate (sidewards - when looking in beam direction)
4887 perpjetplusMom[1]=(TMath::Sin(jetPerpposPhi)*jetPerpPt); //y coordinate (upwards - when looking in beam direction)
4888 perpjetplusMom[2]=jetMom[2]; //z coordinate (along beam axis), invariant under azimuthal rotation
4890 perpjetnegMom[0]=(TMath::Cos(jetPerpnegPhi)*jetPerpPt); //x coordinate (sidewards - when looking in beam direction)
4891 perpjetnegMom[1]=(TMath::Sin(jetPerpnegPhi)*jetPerpPt); //y coordinate (upwards - when looking in beam direction)
4892 perpjetnegMom[2]=jetMom[2]; //z coordinate (along beam axis), invariant under azimuthal rotation
4895 TVector3 perpjetplus3mom(perpjetplusMom); //3-Vector for new perp. jet axis, clockwise rotated
4896 TVector3 perpjetneg3mom(perpjetnegMom); //3-Vector for new perp. jet axis, counterclockwise rotated
4898 //for crosscheck TVector3 rotation method
4900 //Double_t jetMomplusTest[3];
4901 //Double_t jetMomminusTest[3];
4903 //jet3mom.RotateZ(TMath::Pi()*0.5);//rotate original jet axis around +90 degrees in phi
4905 //perpjetminus3momTest = jet3mom.RotateZ((-1)*TMath::Pi()*0.5);
4907 // jet3mom.RotateZ(TMath::Pi()*0.5);
4908 // jet3mom.RotateZ((-1)*TMath::Pi()*0.5);
4910 //jetMomplusTest[0] = jet3mom.X(); //fetching perp. axis coordinates
4911 //jetMomplusTest[1] = jet3mom.Y();
4912 //jetMomplusTest[2] = jet3mom.Z();
4914 //TVector3 perpjetplus3momTest(jetMomplusTest); //new TVector3 for +90deg rotated jet axis with rotation method from ROOT
4915 //TVector3 perpjetminus3momTest(jetMomminusTest); //new TVector3 for -90deg rotated jet axis with rotation method from ROOT
4918 for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){ //collect V0 content in perp cone, rotated clockwise
4920 AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack)); //inputlist is fListK0s, all reconstructed K0s in event
4921 if(!track){std::cout<<"K0s track not found!!!"<<std::endl; continue;}
4923 Double_t trackMom[3];//3-mom of V0 particle
4924 track->PxPyPz(trackMom);
4925 TVector3 track3mom(trackMom);
4927 Double_t dR = perpjetplus3mom.DeltaR(track3mom);
4931 outputlist->Add(track); // output list is jetPerpConeK0list
4933 sumPerpPt += track->Pt();
4940 for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){//collect V0 content in perp cone, rotated counterclockwise
4942 AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack)); //inputlist is fListK0s, all reconstructed K0s in event
4943 if(!track){std::cout<<"K0s track not found!!!"<<std::endl; continue;}
4945 Double_t trackMom[3];//3-mom of V0 particle
4946 track->PxPyPz(trackMom);
4947 TVector3 track3mom(trackMom);
4949 Double_t dR = perpjetneg3mom.DeltaR(track3mom);
4953 outputlist->Add(track); // output list is jetPerpConeK0list
4955 sumPerpPt += track->Pt();
4961 // pay attention: this list contains the double amount of V0s, found in both cones
4962 // before using it, devide spectra by 2!!!
4963 sumPerpPt = sumPerpPt*0.5; //correct to do this?
4971 // _______________________________________________________________________________________________________________________________________________________
4973 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){
4975 TClonesArray *stackmc = 0x0;
4976 stackmc = (TClonesArray*)listmc->FindObject(AliAODMCParticle::StdBranchName()); //get MCAOD branch in data
4979 Printf("ERROR: stack not available");
4984 Int_t negAssLabel = TMath::Abs(trackNeg->GetLabel()); //negative (reconstructed) charged track label in MC stack
4985 Int_t posAssLabel = TMath::Abs(trackPos->GetLabel()); //positive (reconstructed) charged track label in MC stack
4987 //injected particle checks
4993 AliAODMCHeader *header=(AliAODMCHeader*)listmc->FindObject(AliAODMCHeader::StdBranchName());
4994 if(!header)return kFALSE;
4996 mcXv=header->GetVtxX(); mcYv=header->GetVtxY(); mcZv=header->GetVtxZ();
4998 Int_t trackinjected = IsTrackInjected(v0, header, stackmc); //requires AliAODTrack instead of AliVTrack
5000 // std::cout<<"is track injected: "<<trackinjected<<std::endl;
5002 if(trackinjected == 0){std::cout<<"HIJING track injected!!: "<<trackinjected<<std::endl;}
5006 if(negAssLabel>=0 && negAssLabel < stackmc->GetEntriesFast() && posAssLabel>=0 && posAssLabel < stackmc->GetEntriesFast()){//safety check if label has valid value of stack
5008 AliAODMCParticle *mcNegPart =(AliAODMCParticle*)stackmc->UncheckedAt(negAssLabel);//fetch the, with one MC truth track associated (reconstructed), negative charged track
5009 v0Label = mcNegPart->GetMother();// mother of negative charged particle is v0, get v0 label here
5010 negDaughterpdg = mcNegPart->GetPdgCode();
5011 AliAODMCParticle *mcPosPart =(AliAODMCParticle*)stackmc->UncheckedAt(posAssLabel);//fetch the, with one MC truth track associated (reconstructed), positive charged track
5012 Int_t v0PosLabel = mcPosPart->GetMother(); //get mother label of positive charged track label
5013 posDaughterpdg = mcPosPart->GetPdgCode();
5015 if(v0Label >= 0 && v0Label < stackmc->GetEntriesFast() && v0Label == v0PosLabel){//first v0 mc label check, then: check if both daughters are stemming from same particle
5017 AliAODMCParticle *mcv0 = (AliAODMCParticle *)stackmc->UncheckedAt(v0Label); //fetch MC ass. particle to v0 (mother of the both charged daughter tracks)
5019 //do not use anymore:
5020 //fPhysicalPrimary = mcv0->IsPhysicalPrimary();
5023 Float_t fDistPrimaryMax = 0.01; // [cm] max distance of production point to the primary vertex (criterion for choice of MC particles considered as primary)
5025 // Get the distance between production point of the MC mother particle and the primary vertex
5027 Double_t dx = mcXv-mcv0->Xv();//mc primary vertex - mc particle production vertex
5028 Double_t dy = mcYv-mcv0->Yv();
5029 Double_t dz = mcZv-mcv0->Zv();
5031 Float_t fDistPrimary = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
5032 fPhysicalPrimary = kFALSE;//init
5034 fPhysicalPrimary = (fDistPrimary < fDistPrimaryMax);
5036 //if(fPhysicalPrimary == kTRUE){std::cout<<"hello*********!!!!!!!!!!!!! "<<std::endl;}
5037 //if(fPhysicalPrimary == kFALSE)return kFALSE;
5039 MCv0PDGCode = mcv0->GetPdgCode();
5041 //std::cout<<"MCv0PDGCode: "<<MCv0PDGCode<<std::endl;
5043 MCPt = mcv0->Pt();//for MC data, always use MC gen. pt for any pt distributions, also for the spectra, used for normalisation
5044 //for feed-down checks later
5046 Int_t motherLabel = mcv0->GetMother(); //get mother particle label of v0 particle
5047 // std::cout<<"motherLabel: "<<motherLabel<<std::endl;
5049 if(motherLabel >= 0 && v0Label < stackmc->GetEntriesFast()) //do safety check for mother label
5051 AliAODMCParticle *mcMother = (AliAODMCParticle *)stackmc->UncheckedAt(motherLabel); //get mother particle
5052 motherType = mcMother->GetPdgCode(); //get PDG code of mother
5055 Double_t XibarPt = 0.;
5057 if(particletype == kLambda){
5058 if((motherType == 3312)||(motherType == 3322)){ //if v0 mother is Xi0 or Xi- fill MC gen. pt in FD La histogram
5059 XiPt = mcMother->Pt();
5060 fh1MCXiPt->Fill(XiPt);
5063 if(particletype == kAntiLambda){
5064 if((motherType == -3312)||(motherType == -3322)){ //if v0 mother is Xibar0 or Xibar+ fill MC gen. pt in FD ALa histogram
5065 XibarPt = mcMother->Pt();
5066 fh1MCXibarPt->Fill(XibarPt);
5072 //pdg code checks etc..
5074 if(particletype == kK0){
5076 if(TMath::Abs(posDaughterpdg) != 211){return kFALSE;}//one or both of the daughters are not a pion
5077 if(TMath::Abs(negDaughterpdg) != 211){return kFALSE;}
5079 if(MCv0PDGCode != 310) {return kFALSE;}
5082 if(particletype == kLambda){
5083 if(MCv0PDGCode != 3122)return kFALSE;//if particle is not Antilambda, v0 is rejected
5084 if(posDaughterpdg != 2212)return kFALSE;
5085 if(negDaughterpdg != -211)return kFALSE; //pdg code check for Lambda daughters
5087 //{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 //}
5090 if(particletype == kAntiLambda){
5091 if(MCv0PDGCode != -3122)return kFALSE;
5092 if(posDaughterpdg != 211)return kFALSE;
5093 if(negDaughterpdg !=-2212)return kFALSE; //pdg code check for Antilambda daughters
5096 //{if((motherType == -3312)||(motherType == -3322)){continue;}//if bar{Xi0} and Xi+ is motherparticle of Antilambda, particle is rejected
5100 return kTRUE; //check was successful
5101 }//end mc v0 label check
5102 }// end of stack label check
5107 return kFALSE; //check wasn't successful
5109 //________________________________________________________________________________________________________________________________________________________
5112 Bool_t AliAnalysisTaskJetChem::IsParticleMatching(const AliAODMCParticle* mcp0, const Int_t v0Label){
5114 const Int_t mcp0label = mcp0->GetLabel();
5116 if(v0Label == mcp0label)return kTRUE;
5121 //_______________________________________________________________________________________________________________________________________________________
5123 Bool_t AliAnalysisTaskJetChem::DaughterTrackCheck(AliAODv0* v0, Int_t& nnum, Int_t& pnum){
5126 if(v0->GetNDaughters() != 2) return kFALSE;//case v0 has more or less than 2 daughters, avoids seg. break at some AOD files //reason?
5129 // safety check of input parameters
5132 if(fDebug > 1){std::cout << std::endl
5133 << "Warning in AliAnalysisTaskJetChem::DaughterTrackCheck:" << std::endl
5134 << "v0 = " << v0 << std::endl;}
5140 //Daughters track check: its Luke Hanrattys method to check daughters charge
5146 AliAODTrack *ntracktest =(AliAODTrack*)(v0->GetDaughter(nnum));
5148 if(ntracktest == NULL)
5150 if(fDebug > 1){std::cout << std::endl
5151 << "Warning in AliAnalysisTaskJetChem::DaughterTrackCheck:" << std::endl
5152 << "ntracktest = " << ntracktest << std::endl;}
5157 if(ntracktest->Charge() > 0)
5163 const AliAODTrack *trackNeg=(AliAODTrack *)(v0->GetDaughter(nnum));
5164 const AliAODTrack *trackPos=(AliAODTrack *)(v0->GetDaughter(pnum));
5166 //Check if both tracks are available
5167 if (!trackPos || !trackNeg) {
5168 if(fDebug > 1) Printf("strange analysis::UserExec:: Error:Could not retrieve one of the daughter tracks\n");
5173 //remove like sign V0s
5174 if ( trackPos->Charge() == trackNeg->Charge() ){
5175 //if(fDebug>1) Printf("%s:%d found like-sign V0", (char*)__FILE__,__LINE__);
5183 //______________________________________________________________________
5184 TString AliAnalysisTaskJetChem::GetGenerator(Int_t label, AliAODMCHeader* header){
5185 Int_t nsumpart=0;//number of particles
5186 TList *lh=header->GetCocktailHeaders();//TList with all generator headers
5187 Int_t nh=lh->GetEntries();//number of entries in TList with all headers
5189 for(Int_t i=0;i<nh;i++){
5190 AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i);
5191 TString genname=gh->GetName();//name of particle generator
5192 Int_t npart=gh->NProduced();//number of stable or undecayed particles in MC stack block (?)
5193 if(label>=nsumpart && label<(nsumpart+npart)) return genname;
5200 //____________________________________________________________________________________________________________-
5202 Int_t AliAnalysisTaskJetChem::IsTrackInjected(AliAODv0 *v0, AliAODMCHeader *header, TClonesArray *arrayMC){//info in TString should be available from 2011 data productions on..
5204 if(!v0){std::cout << " !part " << std::endl;return 1;}
5205 if(!header){std::cout << " !header " << std::endl;return 1;}
5206 if(!arrayMC){std::cout << " !arrayMC " << std::endl;return 1;}
5208 //comment: all MC truth particles are sorted in the MC stack, first comes a block with all hijing produced particles, then several blocks for particletype specific injected particles, after this comes a (not-ordered) block with all decay-products
5209 //the complete amount of MC truth particles produced by several sources is named 'cocktail'
5211 Int_t lab=v0->GetLabel();//get particle label in MC stack
5212 if(lab<0) {return 1;} //if label is negative, this particle is HIJING produced and not injected
5213 TString generatorName = GetGenerator(lab,header);//this function returns a string with the generator name, used to produce this certain particle
5214 TString empty="";//no header was found
5216 std::cout << " TString generatorName: " << generatorName << std::endl;
5218 //std::cout << " FIRST CALL " << bbb << std::endl;
5220 while(generatorName.IsWhitespace()){
5221 AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab);
5222 if(!mcpart){return 1;}
5223 Int_t mother = mcpart->GetMother();
5225 generatorName = GetGenerator(mother,header);//see function direclty below..
5226 std::cout << "Testing " << generatorName << " " << std::endl;
5229 std::cout << " FINAL CALL " << generatorName << std::endl;
5231 //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
5233 if(generatorName.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"
5240 //_________________________________________________________________________________________________________________________________________
5241 Double_t AliAnalysisTaskJetChem::SmearJetPt(Double_t jetPt, Int_t /*cent*/, Double_t /*jetRadius*/, Double_t /*ptmintrack*/, Double_t& jetPtSmear){
5243 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
5247 /* if(cent>10) cl = 2;
5252 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
5253 //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
5255 //fsmear->SetParameters(1,0,4.472208);// for 2010 PbPb jets, R=0.2, ptmintrack = 0.15 GeV/c, cent 00-10%
5257 /* //delta-pt width for anti-kt jet finder:
5260 if((cl == 1)&&(jetRadius == 0.4)&&(ptmintrack == 0.15)){
5261 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
5263 if((cl == 2)&&(jetRadius == 0.4)&&(ptmintrack == 0.15)){
5264 fsmear->SetParameters(1,0,8.536195);
5266 if((cl == 3)&&(jetRadius == 0.4)&&(ptmintrack == 0.15)){
5267 fsmear->SetParameters(1,0,?);
5269 if((cl == 4)&&(jetRadius == 0.4)&&(ptmintrack == 0.15)){
5270 fsmear->SetParameters(1,0,5.229839);
5274 if((cl == 1)&&(jetRadius == 0.3)&&(ptmintrack == 0.15)){
5275 fsmear->SetParameters(1,0,7.145967);
5277 if((cl == 2)&&(jetRadius == 0.3)&&(ptmintrack == 0.15)){
5278 fsmear->SetParameters(1,0,5.844796);
5280 if((cl == 3)&&(jetRadius == 0.3)&&(ptmintrack == 0.15)){
5281 fsmear->SetParameters(1,0,?);
5283 if((cl == 4)&&(jetRadius == 0.3)&&(ptmintrack == 0.15)){
5284 fsmear->SetParameters(1,0,3.630751);
5288 if((cl == 1)&&(jetRadius == 0.2)&&(ptmintrack == 0.15)){
5289 fsmear->SetParameters(1,0,4.472208);
5291 if((cl == 2)&&(jetRadius == 0.2)&&(ptmintrack == 0.15)){
5292 fsmear->SetParameters(1,0,3.543938);
5294 if((cl == 3)&&(jetRadius == 0.2)&&(ptmintrack == 0.15)){
5295 fsmear->SetParameters(1,0,?);
5297 if((cl == 4)&&(jetRadius == 0.2)&&(ptmintrack == 0.15)){
5298 fsmear->SetParameters(1,0,1.037476);
5303 Double_t r = fsmear.GetRandom();
5304 jetPtSmear = jetPt + r;
5306 // std::cout<<"jetPt: "<<jetPt<<std::endl;
5307 // std::cout<<"jetPtSmear: "<<jetPtSmear<<std::endl;
5308 // std::cout<<"r: "<<r<<std::endl;
5315 //______________________________________________________________________________________________________________________
5316 //____________________________________________________________________________________________________________________
5318 Bool_t AliAnalysisTaskJetChem::IsParticleInCone(const AliVParticle* part1, const AliVParticle* part2, Double_t dRMax) const
5320 // decides whether a particle is inside a jet cone
5321 if (!part1 || !part2)
5324 TVector3 vecMom2(part2->Px(),part2->Py(),part2->Pz());
5325 TVector3 vecMom1(part1->Px(),part1->Py(),part1->Pz());
5326 Double_t dR = vecMom2.DeltaR(vecMom1); // = sqrt(dEta*dEta+dPhi*dPhi)
5327 if(dR<dRMax) // momentum vectors of part1 and part2 are closer than dRMax
5331 //__________________________________________________________________________________________________________________
5334 Bool_t AliAnalysisTaskJetChem::IsRCJCOverlap(TList* recjetlist, const AliVParticle* part, Double_t dDistance) const{
5336 if(!recjetlist) return kFALSE;
5337 if(!part) return kFALSE;
5338 if(!dDistance) return kFALSE;
5339 Int_t nRecJetsCuts = fJetsRecCuts->GetEntries();
5341 for(Int_t i=0; i<nRecJetsCuts; ++i){ //loop over all reconstructed jets in events
5342 AliAODJet* jet = (AliAODJet*) (recjetlist->At(i));
5343 if(!jet){if(fDebug>2)std::cout<<"AliAnalysisTaskJetChem::IsRCJCOverlap jet pointer invalid!"<<std::endl;continue;}
5344 if(IsParticleInCone(jet, part, dDistance) == kTRUE)return kTRUE;//RC and JC are overlapping
5346 }//end loop testing RC-JC overlap
5347 return kFALSE;//RC and JC are not overlapping -> good!
5350 //_______________________________________________________________________________________________________________________
5351 AliAODJet* AliAnalysisTaskJetChem::GetRandomCone(TList* jetlist, Double_t dEtaConeMax, Double_t dDistance) const
5353 TLorentzVector vecRdCone;
5354 AliAODJet* jetRC = 0;//random cone candidate
5355 Double_t dEta, dPhi; //random eta and phi value for RC
5356 Bool_t IsRCoutJC;//check whether RC is not overlapping with any selected jet cone in event
5357 Int_t iRCTrials = 10;//search at maximum 10 times for random cone that doesn't overlap with jet cone
5359 for(Int_t i=0; i<iRCTrials; iRCTrials++){
5361 dEta = dEtaConeMax*(2*fRandom->Rndm()-1.); //random eta value in range: [-dEtaConeMax,+dEtaConeMax]
5362 dPhi = TMath::TwoPi()*fRandom->Rndm(); //random phi value in range: [0,2*Pi]
5363 vecRdCone.SetPtEtaPhiM(1.,dEta,dPhi,0.);
5364 jetRC = new AliAODJet(vecRdCone);//new RC candidate
5366 if (IsRCJCOverlap(jetlist,jetRC,dDistance))
5368 IsRCoutJC = kTRUE; std::cout<<"RC and JC are not overlapping!!!"<<std::endl;
5372 delete jetRC; //RC is overlapping with JC, delete this RC candidate
5375 if(!IsRCoutJC) {jetRC = 0;}//in case no random cone was selected
5381 // _______________________________________________________________________________________________________________________
5382 AliAODJet* AliAnalysisTaskJetChem::GetMedianCluster()
5384 // fill tracks from bckgCluster branch,
5385 // using cluster with median density (odd number of clusters)
5386 // or picking randomly one of the two closest to median (even number)
5388 Int_t nBckgClusters = fBckgJetsRec->GetEntries(); // not 'recCuts': use all clusters in full eta range
5390 if(nBckgClusters<3) return 0; // need at least 3 clusters (skipping 2 highest)
5392 Double_t* bgrDensity = new Double_t[nBckgClusters];
5393 Int_t* indices = new Int_t[nBckgClusters];
5395 for(Int_t ij=0; ij<nBckgClusters; ++ij){
5397 AliAODJet* bgrCluster = (AliAODJet*)(fBckgJetsRec->At(ij));
5398 Double_t clusterPt = bgrCluster->Pt();
5399 Double_t area = bgrCluster->EffectiveAreaCharged();
5401 Double_t density = 0;
5402 if(area>0) density = clusterPt/area;
5404 bgrDensity[ij] = density;
5408 TMath::Sort(nBckgClusters, bgrDensity, indices);
5410 // get median cluster
5412 AliAODJet* medianCluster = 0;
5413 Double_t medianDensity = 0;
5415 if(TMath::Odd(nBckgClusters)){
5417 //Int_t medianIndex = indices[(Int_t) (0.5*(nBckgClusters-1))];
5418 Int_t medianIndex = indices[(Int_t) (0.5*(nBckgClusters+1))];
5420 medianCluster = (AliAODJet*)(fBckgJetsRec->At(medianIndex));
5422 Double_t clusterPt = medianCluster->Pt();
5423 Double_t area = medianCluster->EffectiveAreaCharged();
5425 if(area>0) medianDensity = clusterPt/area;
5429 //Int_t medianIndex1 = indices[(Int_t) (0.5*nBckgClusters-1)];
5430 //Int_t medianIndex2 = indices[(Int_t) (0.5*nBckgClusters)];
5432 Int_t medianIndex1 = indices[(Int_t) (0.5*nBckgClusters)];
5433 Int_t medianIndex2 = indices[(Int_t) (0.5*nBckgClusters+1)];
5435 AliAODJet* medianCluster1 = (AliAODJet*)(fBckgJetsRec->At(medianIndex1));
5436 AliAODJet* medianCluster2 = (AliAODJet*)(fBckgJetsRec->At(medianIndex2));
5438 Double_t density1 = 0;
5439 Double_t clusterPt1 = medianCluster1->Pt();
5440 Double_t area1 = medianCluster1->EffectiveAreaCharged();
5441 if(area1>0) density1 = clusterPt1/area1;
5443 Double_t density2 = 0;
5444 Double_t clusterPt2 = medianCluster2->Pt();
5445 Double_t area2 = medianCluster2->EffectiveAreaCharged();
5446 if(area2>0) density2 = clusterPt2/area2;
5448 medianDensity = 0.5*(density1+density2);
5450 medianCluster = ( (gRandom->Rndm()>0.5) ? medianCluster1 : medianCluster2 ); // select one randomly to avoid adding areas
5453 delete[] bgrDensity;
5456 return medianCluster;
5459 //____________________________________________________________________________________________
5461 Double_t AliAnalysisTaskJetChem::AreaCircSegment(Double_t dRadius, Double_t dDistance) const
5463 // calculate area of a circular segment defined by the circle radius and the (oriented) distance between the secant line and the circle centre
5464 Double_t dEpsilon = 1e-2;
5465 Double_t dR = dRadius;
5466 Double_t dD = dDistance;
5467 if (TMath::Abs(dR)<dEpsilon)
5469 if(fDebug>0) printf("AliAnalysisTaskJetChem::AreaCircSegment: Error: Too small radius: %f < %f\n",dR,dEpsilon);
5475 return TMath::Pi()*dR*dR;
5476 return dR*dR*TMath::ACos(dD/dR)-dD*TMath::Sqrt(dR*dR-dD*dD);