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)
238 ,fhnrecMCHijingLaIncl(0)
239 ,fhnrecMCHijingLaCone(0)
240 ,fhnrecMCHijingALaIncl(0)
241 ,fhnrecMCHijingALaCone(0)
242 ,fhnrecMCInjectLaIncl(0)
243 ,fhnrecMCInjectLaCone(0)
244 ,fhnrecMCInjectALaIncl(0)
245 ,fhnrecMCInjectALaCone(0)
249 ,fhnMCrecK0ConeSmear(0)
250 ,fhnMCrecLaConeSmear(0)
251 ,fhnMCrecALaConeSmear(0)
252 ,fhnK0sSecContinCone(0)
253 ,fhnLaSecContinCone(0)
254 ,fhnALaSecContinCone(0)
276 ,fh1MCMultiplicityPrimary(0)
277 ,fh1MCMultiplicityTracks(0)
280 ,fhnFeedDownLaCone(0)
281 ,fhnFeedDownALaCone(0)
282 ,fh1MCProdRadiusK0s(0)
283 ,fh1MCProdRadiusLambda(0)
284 ,fh1MCProdRadiusAntiLambda(0)
288 ,fh1MCPtAntiLambda(0)
296 //,fh1MCRapAntiLambda(0)
300 ,fh1MCEtaAntiLambda(0)
303 // default constructor
306 //__________________________________________________________________________________________
307 AliAnalysisTaskJetChem::AliAnalysisTaskJetChem(const char *name)
308 : AliAnalysisTaskFragmentationFunction(name)
322 ,fCutV0cosPointAngle(0)
329 ,fCutDcaV0Daughters(0)
330 ,fCutDcaPosToPrimVertex(0)
331 ,fCutDcaNegToPrimVertex(0)
341 ,fFFHistosRecCutsK0Evt(0)
342 //,fFFHistosIMK0AllEvt(0)
343 //,fFFHistosIMK0Jet(0)
344 //,fFFHistosIMK0Cone(0)
348 //,fFFHistosIMLaAllEvt(0)
349 //,fFFHistosIMLaJet(0)
350 //,fFFHistosIMLaCone(0)
354 ,fListFeeddownLaCand(0)
355 ,fListFeeddownALaCand(0)
361 ,fListMCgenK0sCone(0)
363 ,fListMCgenALaCone(0)
364 ,IsArmenterosSelected(0)
365 //,fFFHistosIMALaAllEvt(0)
366 //,fFFHistosIMALaJet(0)
367 // ,fFFHistosIMALaCone(0)
383 ,fFFIMLaNBinsJetPt(0)
414 // ,fh1trackPosNCls(0)
415 // ,fh1trackNegNCls(0)
425 ,fh2ProperLifetimeK0sVsPtBeforeCut(0)
426 ,fh2ProperLifetimeK0sVsPtAfterCut(0)
428 ,fh1DcaV0Daughters(0)
429 ,fh1DcaPosToPrimVertex(0)
430 ,fh1DcaNegToPrimVertex(0)
431 ,fh2ArmenterosBeforeCuts(0)
432 ,fh2ArmenterosAfterCuts(0)
435 ,fh1PosDaughterCharge(0)
436 ,fh1NegDaughterCharge(0)
443 ,fhnInvMassEtaTrackPtK0s(0)
444 ,fhnInvMassEtaTrackPtLa(0)
445 ,fhnInvMassEtaTrackPtALa(0)
454 ,fh2MCEtagenK0Cone(0)
455 ,fh2MCEtagenLaCone(0)
456 ,fh2MCEtagenALaCone(0)
459 ,fh1IMALaConeSmear(0)
460 ,fhnrecMCHijingLaIncl(0)
461 ,fhnrecMCHijingLaCone(0)
462 ,fhnrecMCHijingALaIncl(0)
463 ,fhnrecMCHijingALaCone(0)
464 ,fhnrecMCInjectLaIncl(0)
465 ,fhnrecMCInjectLaCone(0)
466 ,fhnrecMCInjectALaIncl(0)
467 ,fhnrecMCInjectALaCone(0)
471 ,fhnMCrecK0ConeSmear(0)
472 ,fhnMCrecLaConeSmear(0)
473 ,fhnMCrecALaConeSmear(0)
474 ,fhnK0sSecContinCone(0)
475 ,fhnLaSecContinCone(0)
476 ,fhnALaSecContinCone(0)
498 ,fh1MCMultiplicityPrimary(0)
499 ,fh1MCMultiplicityTracks(0)
502 ,fhnFeedDownLaCone(0)
503 ,fhnFeedDownALaCone(0)
504 ,fh1MCProdRadiusK0s(0)
505 ,fh1MCProdRadiusLambda(0)
506 ,fh1MCProdRadiusAntiLambda(0)
510 ,fh1MCPtAntiLambda(0)
518 //,fh1MCRapAntiLambda(0)
522 ,fh1MCEtaAntiLambda(0)
528 DefineOutput(1,TList::Class());
531 //__________________________________________________________________________________________________________________________
532 AliAnalysisTaskJetChem::AliAnalysisTaskJetChem(const AliAnalysisTaskJetChem ©)
533 : AliAnalysisTaskFragmentationFunction()
535 ,fRandom(copy.fRandom)
536 ,fAnalysisMC(copy.fAnalysisMC)
537 ,fDeltaVertexZ(copy.fDeltaVertexZ)
538 ,fCutjetEta(copy.fCutjetEta)
539 ,fCuttrackNegNcls(copy.fCuttrackNegNcls)
540 ,fCuttrackPosNcls(copy.fCuttrackPosNcls)
541 ,fCutPostrackRap(copy.fCutPostrackRap)
542 ,fCutNegtrackRap(copy.fCutNegtrackRap)
543 ,fCutRap(copy.fCutRap)
544 ,fCutPostrackEta(copy.fCutPostrackEta)
545 ,fCutNegtrackEta(copy.fCutNegtrackEta)
546 ,fCutEta(copy.fCutEta)
547 ,fCutV0cosPointAngle(copy.fCutV0cosPointAngle)
548 ,fKinkDaughters(copy.fKinkDaughters)
549 ,fRequireTPCRefit(copy.fRequireTPCRefit)
550 ,fCutArmenteros(copy.fCutArmenteros)
551 ,fCutV0DecayMin(copy.fCutV0DecayMin)
552 ,fCutV0DecayMax(copy.fCutV0DecayMax)
553 ,fCutV0totMom(copy.fCutV0totMom)
554 ,fCutDcaV0Daughters(copy.fCutDcaV0Daughters)
555 ,fCutDcaPosToPrimVertex(copy.fCutDcaPosToPrimVertex)
556 ,fCutDcaNegToPrimVertex(copy.fCutDcaNegToPrimVertex)
557 ,fCutV0RadiusMin(copy.fCutV0RadiusMin)
558 ,fCutV0RadiusMax(copy.fCutV0RadiusMax)
559 ,fCutBetheBloch(copy.fCutBetheBloch)
560 ,fCutRatio(copy.fCutRatio)
561 ,fK0Type(copy.fK0Type)
562 ,fFilterMaskK0(copy.fFilterMaskK0)
563 ,fListK0s(copy.fListK0s)
564 ,fPIDResponse(copy.fPIDResponse)
565 ,fV0QAK0(copy.fV0QAK0)
566 ,fFFHistosRecCutsK0Evt(copy.fFFHistosRecCutsK0Evt)
567 //,fFFHistosIMK0AllEvt(copy.fFFHistosIMK0AllEvt)
568 //,fFFHistosIMK0Jet(copy.fFFHistosIMK0Jet)
569 //,fFFHistosIMK0Cone(copy.fFFHistosIMK0Cone)
570 ,fLaType(copy.fLaType)
571 ,fFilterMaskLa(copy.fFilterMaskLa)
572 ,fListLa(copy.fListLa)
573 //,fFFHistosIMLaAllEvt(copy.fFFHistosIMLaAllEvt)
574 //,fFFHistosIMLaJet(copy.fFFHistosIMLaJet)
575 //,fFFHistosIMLaCone(copy.fFFHistosIMLaCone)
576 ,fALaType(copy.fALaType)
577 ,fFilterMaskALa(copy.fFilterMaskALa)
578 ,fListALa(copy.fListALa)
579 ,fListFeeddownLaCand(copy.fListFeeddownLaCand)
580 ,fListFeeddownALaCand(copy.fListFeeddownALaCand)
581 ,jetConeFDLalist(copy.jetConeFDLalist)
582 ,jetConeFDALalist(copy.jetConeFDALalist)
583 ,fListMCgenK0s(copy.fListMCgenK0s)
584 ,fListMCgenLa(copy.fListMCgenLa)
585 ,fListMCgenALa(copy.fListMCgenALa)
586 ,fListMCgenK0sCone(copy.fListMCgenK0sCone)
587 ,fListMCgenLaCone(copy.fListMCgenLaCone)
588 ,fListMCgenALaCone(copy.fListMCgenALaCone)
589 ,IsArmenterosSelected(copy.IsArmenterosSelected)
590 //,fFFHistosIMALaAllEvt(copy.fFFHistosIMALaAllEvt)
591 //,fFFHistosIMALaJet(copy.fFFHistosIMALaJet)
592 //,fFFHistosIMALaCone(copy.fFFHistosIMALaCone)
593 ,fFFIMNBinsJetPt(copy.fFFIMNBinsJetPt)
594 ,fFFIMJetPtMin(copy.fFFIMJetPtMin)
595 ,fFFIMJetPtMax(copy.fFFIMJetPtMax)
596 ,fFFIMNBinsInvM(copy.fFFIMNBinsInvM)
597 ,fFFIMInvMMin(copy.fFFIMInvMMin)
598 ,fFFIMInvMMax(copy.fFFIMInvMMax)
599 ,fFFIMNBinsPt(copy.fFFIMNBinsPt)
600 ,fFFIMPtMin(copy.fFFIMPtMin)
601 ,fFFIMPtMax(copy.fFFIMPtMax)
602 ,fFFIMNBinsXi(copy.fFFIMNBinsXi)
603 ,fFFIMXiMin(copy.fFFIMXiMin)
604 ,fFFIMXiMax(copy.fFFIMXiMax)
605 ,fFFIMNBinsZ(copy.fFFIMNBinsZ)
606 ,fFFIMZMin(copy.fFFIMZMin)
607 ,fFFIMZMax(copy.fFFIMZMax)
608 ,fFFIMLaNBinsJetPt(copy.fFFIMLaNBinsJetPt)
609 ,fFFIMLaJetPtMin(copy.fFFIMLaJetPtMin)
610 ,fFFIMLaJetPtMax(copy.fFFIMLaJetPtMax)
611 ,fFFIMLaNBinsInvM(copy.fFFIMLaNBinsInvM)
612 ,fFFIMLaInvMMin(copy.fFFIMLaInvMMin)
613 ,fFFIMLaInvMMax(copy.fFFIMLaInvMMax)
614 ,fFFIMLaNBinsPt(copy.fFFIMLaNBinsPt)
615 ,fFFIMLaPtMin(copy.fFFIMLaPtMin)
616 ,fFFIMLaPtMax(copy.fFFIMLaPtMax)
617 ,fFFIMLaNBinsXi(copy.fFFIMLaNBinsXi)
618 ,fFFIMLaXiMin(copy.fFFIMLaXiMin)
619 ,fFFIMLaXiMax(copy.fFFIMLaXiMax)
620 ,fFFIMLaNBinsZ(copy.fFFIMLaNBinsZ)
621 ,fFFIMLaZMin(copy.fFFIMLaZMin)
622 ,fFFIMLaZMax(copy.fFFIMLaZMax)
623 ,fh1EvtAllCent(copy.fh1EvtAllCent)
625 ,fh1K0Mult(copy.fh1K0Mult)
626 ,fh1dPhiJetK0(copy.fh1dPhiJetK0)
627 ,fh1LaMult(copy.fh1LaMult)
628 ,fh1dPhiJetLa(copy.fh1dPhiJetLa)
629 ,fh1ALaMult(copy.fh1ALaMult)
630 ,fh1dPhiJetALa(copy.fh1dPhiJetALa)
631 ,fh1JetEta(copy.fh1JetEta)
632 ,fh1JetPhi(copy.fh1JetPhi)
633 ,fh2JetEtaPhi(copy.fh2JetEtaPhi)
634 //,fh1V0JetPt(copy.fh1V0JetPt)
635 ,fh1IMK0Cone(copy.fh1IMK0Cone)
636 ,fh1IMLaCone(copy.fh1IMLaCone)
637 ,fh1IMALaCone(copy.fh1IMALaCone)
638 ,fh2FFJetTrackEta(copy.fh2FFJetTrackEta)
639 //,fh1trackPosNCls(copy.fh1trackPosNCls)
640 //,fh1trackNegNCls(copy.fh1trackNegNCls)
641 ,fh1trackPosRap(copy.fh1trackPosRap)
642 ,fh1trackNegRap(copy.fh1trackNegRap)
643 //,fh1V0Rap(copy.fh1V0Rap)
644 ,fh1trackPosEta(copy.fh1trackPosEta)
645 ,fh1trackNegEta(copy.fh1trackNegEta)
646 ,fh1V0Eta(copy.fh1V0Eta)
647 //,fh1V0totMom(copy.fh1V0totMom)
648 ,fh1CosPointAngle(copy.fh1CosPointAngle)
649 ,fh1DecayLengthV0(copy.fh1DecayLengthV0)
650 ,fh2ProperLifetimeK0sVsPtBeforeCut(copy.fh2ProperLifetimeK0sVsPtBeforeCut)
651 ,fh2ProperLifetimeK0sVsPtAfterCut(copy.fh2ProperLifetimeK0sVsPtAfterCut)
652 ,fh1V0Radius(copy.fh1V0Radius)
653 ,fh1DcaV0Daughters(copy.fh1DcaV0Daughters)
654 ,fh1DcaPosToPrimVertex(copy.fh1DcaPosToPrimVertex)
655 ,fh1DcaNegToPrimVertex(copy.fh1DcaNegToPrimVertex)
656 ,fh2ArmenterosBeforeCuts(copy.fh2ArmenterosBeforeCuts)
657 ,fh2ArmenterosAfterCuts(copy.fh2ArmenterosAfterCuts)
658 ,fh2BBLaPos(copy.fh2BBLaPos)
659 ,fh2BBLaNeg(copy.fh2BBLaPos)
660 ,fh1PosDaughterCharge(copy.fh1PosDaughterCharge)
661 ,fh1NegDaughterCharge(copy.fh1NegDaughterCharge)
662 ,fh1PtMCK0s(copy.fh1PtMCK0s)
663 ,fh1PtMCLa(copy.fh1PtMCLa)
664 ,fh1PtMCALa(copy.fh1PtMCALa)
665 ,fh1EtaK0s(copy.fh1EtaK0s)
666 ,fh1EtaLa(copy.fh1EtaLa)
667 ,fh1EtaALa(copy.fh1EtaALa)
668 ,fhnInvMassEtaTrackPtK0s(copy.fhnInvMassEtaTrackPtK0s)
669 ,fhnInvMassEtaTrackPtLa(copy.fhnInvMassEtaTrackPtLa)
670 ,fhnInvMassEtaTrackPtALa(copy.fhnInvMassEtaTrackPtALa)
671 ,fh1TrackMultCone(copy.fh1TrackMultCone)
672 ,fh2TrackMultCone(copy.fh2TrackMultCone)
673 ,fhnNJK0(copy.fhnNJK0)
674 ,fhnNJLa(copy.fhnNJLa)
675 ,fhnNJALa(copy.fhnNJALa)
676 ,fh2MCgenK0Cone(copy.fh2MCgenK0Cone)
677 ,fh2MCgenLaCone(copy.fh2MCgenLaCone)
678 ,fh2MCgenALaCone(copy.fh2MCgenALaCone)
679 ,fh2MCEtagenK0Cone(copy.fh2MCEtagenK0Cone)
680 ,fh2MCEtagenLaCone(copy.fh2MCEtagenLaCone)
681 ,fh2MCEtagenALaCone(copy.fh2MCEtagenALaCone)
682 ,fh1IMK0ConeSmear(copy.fh1IMK0ConeSmear)
683 ,fh1IMLaConeSmear(copy.fh1IMLaConeSmear)
684 ,fh1IMALaConeSmear(copy.fh1IMALaConeSmear)
685 ,fhnrecMCHijingLaIncl(copy.fhnrecMCHijingLaIncl)
686 ,fhnrecMCHijingLaCone(copy.fhnrecMCHijingLaCone)
687 ,fhnrecMCHijingALaIncl(copy.fhnrecMCHijingALaIncl)
688 ,fhnrecMCHijingALaCone(copy.fhnrecMCHijingALaCone)
689 ,fhnrecMCInjectLaIncl(copy.fhnrecMCInjectLaIncl)
690 ,fhnrecMCInjectLaCone(copy.fhnrecMCInjectLaCone)
691 ,fhnrecMCInjectALaIncl(copy.fhnrecMCInjectALaIncl)
692 ,fhnrecMCInjectALaCone(copy.fhnrecMCInjectALaCone)
693 ,fhnMCrecK0Cone(copy.fhnMCrecK0Cone)
694 ,fhnMCrecLaCone(copy.fhnMCrecLaCone)
695 ,fhnMCrecALaCone(copy.fhnMCrecALaCone)
696 ,fhnMCrecK0ConeSmear(copy.fhnMCrecK0ConeSmear)
697 ,fhnMCrecLaConeSmear(copy.fhnMCrecLaConeSmear)
698 ,fhnMCrecALaConeSmear(copy.fhnMCrecALaConeSmear)
699 ,fhnK0sSecContinCone(copy.fhnK0sSecContinCone)
700 ,fhnLaSecContinCone(copy.fhnLaSecContinCone)
701 ,fhnALaSecContinCone(copy.fhnALaSecContinCone)
702 ,fhnK0sIncl(copy.fhnK0sIncl)
703 ,fhnK0sCone(copy.fhnK0sCone)
704 ,fhnLaIncl(copy.fhnLaIncl)
705 ,fhnLaCone(copy.fhnLaCone)
706 ,fhnALaIncl(copy.fhnALaIncl)
707 ,fhnALaCone(copy.fhnALaCone)
708 ,fhnK0sPC(copy.fhnK0sPC)
709 ,fhnLaPC(copy.fhnLaPC)
710 ,fhnALaPC(copy.fhnALaPC)
711 ,fhnK0sMCC(copy.fhnK0sMCC)
712 ,fhnLaMCC(copy.fhnLaMCC)
713 ,fhnALaMCC(copy.fhnALaMCC)
714 ,fhnK0sRC(copy.fhnK0sRC)
715 ,fhnLaRC(copy.fhnLaRC)
716 ,fhnALaRC(copy.fhnALaRC)
717 ,fhnK0sOC(copy.fhnK0sOC)
718 ,fhnLaOC(copy.fhnLaOC)
719 ,fhnALaOC(copy.fhnALaOC)
720 ,fh1AreaExcluded(copy.fh1AreaExcluded)
721 ,fh1MedianEta(copy.fh1MedianEta)
722 ,fh1JetPtMedian(copy.fh1JetPtMedian)
723 ,fh1MCMultiplicityPrimary(copy.fh1MCMultiplicityPrimary)
724 ,fh1MCMultiplicityTracks(copy.fh1MCMultiplicityTracks)
725 ,fhnFeedDownLa(copy.fhnFeedDownLa)
726 ,fhnFeedDownALa(copy.fhnFeedDownALa)
727 ,fhnFeedDownLaCone(copy.fhnFeedDownLaCone)
728 ,fhnFeedDownALaCone(copy.fhnFeedDownALaCone)
729 ,fh1MCProdRadiusK0s(copy.fh1MCProdRadiusK0s)
730 ,fh1MCProdRadiusLambda(copy.fh1MCProdRadiusLambda)
731 ,fh1MCProdRadiusAntiLambda(copy.fh1MCProdRadiusAntiLambda)
732 ,fh1MCPtV0s(copy.fh1MCPtV0s)
733 ,fh1MCPtK0s(copy.fh1MCPtK0s)
734 ,fh1MCPtLambda(copy.fh1MCPtLambda)
735 ,fh1MCPtAntiLambda(copy.fh1MCPtAntiLambda)
736 ,fh1MCXiPt(copy.fh1MCXiPt)
737 ,fh1MCXibarPt(copy.fh1MCXibarPt)
738 ,fh2MCEtaVsPtK0s(copy.fh2MCEtaVsPtK0s)
739 ,fh2MCEtaVsPtLa(copy.fh2MCEtaVsPtLa)
740 ,fh2MCEtaVsPtALa(copy.fh2MCEtaVsPtALa)
741 //,fh1MCRapK0s(copy.fh1MCRapK0s)
742 //,fh1MCRapLambda(copy.fh1MCRapLambda)
743 //,fh1MCRapAntiLambda(copy.fh1MCRapAntiLambda)
744 ,fh1MCEtaAllK0s(copy.fh1MCEtaAllK0s)
745 ,fh1MCEtaK0s(copy.fh1MCEtaK0s)
746 ,fh1MCEtaLambda(copy.fh1MCEtaLambda)
747 ,fh1MCEtaAntiLambda(copy.fh1MCEtaAntiLambda)
754 // _________________________________________________________________________________________________________________________________
755 AliAnalysisTaskJetChem& AliAnalysisTaskJetChem::operator=(const AliAnalysisTaskJetChem& o)
760 AliAnalysisTaskFragmentationFunction::operator=(o);
763 fAnalysisMC = o.fAnalysisMC;
764 fDeltaVertexZ = o.fDeltaVertexZ;
765 fCutjetEta = o.fCutjetEta;
766 fCuttrackNegNcls = o.fCuttrackNegNcls;
767 fCuttrackPosNcls = o.fCuttrackPosNcls;
768 fCutPostrackRap = o.fCutPostrackRap;
769 fCutNegtrackRap = o.fCutNegtrackRap;
771 fCutPostrackEta = o.fCutPostrackEta;
772 fCutNegtrackEta = o.fCutNegtrackEta;
774 fCutV0cosPointAngle = o.fCutV0cosPointAngle;
775 fKinkDaughters = o.fKinkDaughters;
776 fRequireTPCRefit = o.fRequireTPCRefit;
777 fCutArmenteros = o.fCutArmenteros;
778 fCutV0DecayMin = o.fCutV0DecayMin;
779 fCutV0DecayMax = o.fCutV0DecayMax;
780 fCutV0totMom = o.fCutV0totMom;
781 fCutDcaV0Daughters = o.fCutDcaV0Daughters;
782 fCutDcaPosToPrimVertex = o.fCutDcaPosToPrimVertex;
783 fCutDcaNegToPrimVertex = o.fCutDcaNegToPrimVertex;
784 fCutV0RadiusMin = o.fCutV0RadiusMin;
785 fCutV0RadiusMax = o.fCutV0RadiusMax;
786 fCutBetheBloch = o.fCutBetheBloch;
787 fCutRatio = o.fCutRatio;
789 fFilterMaskK0 = o.fFilterMaskK0;
790 fListK0s = o.fListK0s;
791 fPIDResponse = o.fPIDResponse;
793 fFFHistosRecCutsK0Evt = o.fFFHistosRecCutsK0Evt;
794 //fFFHistosIMK0AllEvt = o.fFFHistosIMK0AllEvt;
795 //fFFHistosIMK0Jet = o.fFFHistosIMK0Jet;
796 //fFFHistosIMK0Cone = o.fFFHistosIMK0Cone;
798 fFilterMaskLa = o.fFilterMaskLa;
800 //fFFHistosIMLaAllEvt = o.fFFHistosIMLaAllEvt;
801 //fFFHistosIMLaJet = o.fFFHistosIMLaJet;
802 //fFFHistosIMLaCone = o.fFFHistosIMLaCone;
803 fALaType = o.fALaType;
804 fFilterMaskALa = o.fFilterMaskALa;
805 fListFeeddownLaCand = o.fListFeeddownLaCand;
806 fListFeeddownALaCand = o.fListFeeddownALaCand;
807 jetConeFDLalist = o.jetConeFDLalist;
808 jetConeFDALalist = o.jetConeFDALalist;
809 fListMCgenK0s = o.fListMCgenK0s;
810 fListMCgenLa = o.fListMCgenLa;
811 fListMCgenALa = o.fListMCgenALa;
812 fListMCgenK0sCone = o.fListMCgenK0sCone;
813 fListMCgenLaCone = o.fListMCgenLaCone;
814 fListMCgenALaCone = o.fListMCgenALaCone;
815 IsArmenterosSelected = o.IsArmenterosSelected;
816 // fFFHistosIMALaAllEvt = o.fFFHistosIMALaAllEvt;
817 // fFFHistosIMALaJet = o.fFFHistosIMALaJet;
818 // fFFHistosIMALaCone = o.fFFHistosIMALaCone;
819 fFFIMNBinsJetPt = o.fFFIMNBinsJetPt;
820 fFFIMJetPtMin = o.fFFIMJetPtMin;
821 fFFIMJetPtMax = o.fFFIMJetPtMax;
822 fFFIMNBinsPt = o.fFFIMNBinsPt;
823 fFFIMPtMin = o.fFFIMPtMin;
824 fFFIMPtMax = o.fFFIMPtMax;
825 fFFIMNBinsXi = o.fFFIMNBinsXi;
826 fFFIMXiMin = o.fFFIMXiMin;
827 fFFIMXiMax = o.fFFIMXiMax;
828 fFFIMNBinsZ = o.fFFIMNBinsZ;
829 fFFIMZMin = o.fFFIMZMin;
830 fFFIMZMax = o.fFFIMZMax;
831 fFFIMLaNBinsJetPt = o.fFFIMLaNBinsJetPt;
832 fFFIMLaJetPtMin = o.fFFIMLaJetPtMin;
833 fFFIMLaJetPtMax = o.fFFIMLaJetPtMax;
834 fFFIMLaNBinsPt = o.fFFIMLaNBinsPt;
835 fFFIMLaPtMin = o.fFFIMLaPtMin;
836 fFFIMLaPtMax = o.fFFIMLaPtMax;
837 fFFIMLaNBinsXi = o.fFFIMLaNBinsXi;
838 fFFIMLaXiMin = o.fFFIMLaXiMin;
839 fFFIMLaXiMax = o.fFFIMLaXiMax;
840 fFFIMLaNBinsZ = o.fFFIMLaNBinsZ;
841 fFFIMLaZMin = o.fFFIMLaZMin;
842 fFFIMLaZMax = o.fFFIMLaZMax;
843 fh1EvtAllCent = o.fh1EvtAllCent;
845 fh1K0Mult = o.fh1K0Mult;
846 fh1dPhiJetK0 = o.fh1dPhiJetK0;
847 fh1LaMult = o.fh1LaMult;
848 fh1dPhiJetLa = o.fh1dPhiJetLa;
849 fh1ALaMult = o.fh1ALaMult;
850 fh1dPhiJetALa = o.fh1dPhiJetALa;
851 fh1JetEta = o.fh1JetEta;
852 fh1JetPhi = o.fh1JetPhi;
853 fh2JetEtaPhi = o.fh2JetEtaPhi;
854 //fh1V0JetPt = o.fh1V0JetPt;
855 fh1IMK0Cone = o.fh1IMK0Cone;
856 fh1IMLaCone = o.fh1IMLaCone;
857 fh1IMALaCone = o.fh1IMALaCone;
858 fh2FFJetTrackEta = o.fh2FFJetTrackEta;
859 //fh1trackPosNCls = o.fh1trackPosNCls;
860 //fh1trackNegNCls = o.fh1trackNegNCls;
861 fh1trackPosRap = o.fh1trackPosRap;
862 fh1trackNegRap = o.fh1trackNegRap;
863 //fh1V0Rap = o.fh1V0Rap;
864 fh1trackPosEta = o.fh1trackPosEta;
865 fh1trackNegEta = o.fh1trackNegEta;
866 fh1V0Eta = o.fh1V0Eta;
867 // fh1V0totMom = o.fh1V0totMom;
868 fh1CosPointAngle = o.fh1CosPointAngle;
869 fh1DecayLengthV0 = o.fh1DecayLengthV0;
870 fh2ProperLifetimeK0sVsPtBeforeCut = o.fh2ProperLifetimeK0sVsPtBeforeCut;
871 fh2ProperLifetimeK0sVsPtAfterCut= o.fh2ProperLifetimeK0sVsPtAfterCut;
872 fh1V0Radius = o.fh1V0Radius;
873 fh1DcaV0Daughters = o.fh1DcaV0Daughters;
874 fh1DcaPosToPrimVertex = o.fh1DcaPosToPrimVertex;
875 fh1DcaNegToPrimVertex = o.fh1DcaNegToPrimVertex;
876 fh2ArmenterosBeforeCuts = o.fh2ArmenterosBeforeCuts;
877 fh2ArmenterosAfterCuts = o.fh2ArmenterosAfterCuts;
878 fh2BBLaPos = o.fh2BBLaPos;
879 fh2BBLaNeg = o.fh2BBLaPos;
880 fh1PosDaughterCharge = o.fh1PosDaughterCharge;
881 fh1NegDaughterCharge = o.fh1NegDaughterCharge;
882 fh1PtMCK0s = o.fh1PtMCK0s;
883 fh1PtMCLa = o.fh1PtMCLa;
884 fh1PtMCALa = o.fh1PtMCALa;
885 fh1EtaK0s = o.fh1EtaK0s;
886 fh1EtaLa = o.fh1EtaLa;
887 fh1EtaALa = o.fh1EtaALa;
888 fhnInvMassEtaTrackPtK0s = o.fhnInvMassEtaTrackPtK0s;
889 fhnInvMassEtaTrackPtLa = o.fhnInvMassEtaTrackPtLa;
890 fhnInvMassEtaTrackPtALa = o.fhnInvMassEtaTrackPtALa;
891 fh1TrackMultCone = o.fh1TrackMultCone;
892 fh2TrackMultCone = o.fh2TrackMultCone;
895 fhnNJALa = o.fhnNJALa;
896 fh2MCgenK0Cone = o.fh2MCgenK0Cone;
897 fh2MCgenLaCone = o.fh2MCgenLaCone;
898 fh2MCgenALaCone = o.fh2MCgenALaCone;
899 fh2MCEtagenK0Cone = o.fh2MCEtagenK0Cone;
900 fh2MCEtagenLaCone = o.fh2MCEtagenLaCone;
901 fh2MCEtagenALaCone = o.fh2MCEtagenALaCone;
902 fh1IMK0ConeSmear = o.fh1IMK0ConeSmear;
903 fh1IMLaConeSmear = o.fh1IMLaConeSmear;
904 fh1IMALaConeSmear = o.fh1IMALaConeSmear;
905 fhnrecMCHijingLaIncl = o.fhnrecMCHijingLaIncl;
906 fhnrecMCHijingLaCone = o.fhnrecMCHijingLaCone;
907 fhnrecMCHijingALaIncl = o.fhnrecMCHijingALaIncl;
908 fhnrecMCHijingALaCone = o.fhnrecMCHijingALaCone;
909 fhnrecMCInjectLaIncl = o.fhnrecMCInjectLaIncl;
910 fhnrecMCInjectLaCone = o.fhnrecMCInjectLaCone;
911 fhnrecMCInjectALaIncl = o.fhnrecMCInjectALaIncl;
912 fhnrecMCInjectALaCone = o.fhnrecMCInjectALaCone;
913 fhnMCrecK0Cone = o.fhnMCrecK0Cone;
914 fhnMCrecLaCone = o.fhnMCrecLaCone;
915 fhnMCrecALaCone = o.fhnMCrecALaCone;
916 fhnMCrecK0ConeSmear = o.fhnMCrecK0ConeSmear;
917 fhnMCrecLaConeSmear = o.fhnMCrecLaConeSmear;
918 fhnMCrecALaConeSmear = o.fhnMCrecALaConeSmear;
919 fhnK0sSecContinCone = o.fhnK0sSecContinCone;
920 fhnLaSecContinCone = o.fhnLaSecContinCone;
921 fhnALaSecContinCone = o.fhnALaSecContinCone;
922 fhnK0sIncl = o.fhnK0sIncl;
923 fhnK0sCone = o.fhnK0sCone;
924 fhnLaIncl = o.fhnLaIncl;
925 fhnLaCone = o.fhnLaCone;
926 fhnALaIncl = o.fhnALaIncl;
927 fhnALaCone = o.fhnALaCone;
928 fhnK0sPC = o.fhnK0sPC;
930 fhnALaPC = o.fhnALaPC;
931 fhnK0sRC = o.fhnK0sRC;
933 fhnALaRC = o.fhnALaRC;
934 fhnK0sOC = o.fhnK0sOC;
936 fhnALaOC = o.fhnALaOC;
937 fh1AreaExcluded = o.fh1AreaExcluded;
938 fh1MedianEta = o.fh1MedianEta;
939 fh1JetPtMedian = o.fh1JetPtMedian;
940 fh1MCMultiplicityPrimary = o.fh1MCMultiplicityPrimary;
941 fh1MCMultiplicityTracks = o.fh1MCMultiplicityTracks;
942 fhnFeedDownLa = o.fhnFeedDownLa;
943 fhnFeedDownALa = o.fhnFeedDownALa;
944 fhnFeedDownLaCone = o.fhnFeedDownLaCone;
945 fhnFeedDownALaCone = o.fhnFeedDownALaCone;
946 fh1MCProdRadiusK0s = o.fh1MCProdRadiusK0s;
947 fh1MCProdRadiusLambda = o.fh1MCProdRadiusLambda;
948 fh1MCProdRadiusAntiLambda = o.fh1MCProdRadiusAntiLambda;
949 fh1MCPtV0s = o.fh1MCPtV0s;
950 fh1MCPtK0s = o.fh1MCPtK0s;
951 fh1MCPtLambda = o.fh1MCPtLambda;
952 fh1MCPtAntiLambda = o.fh1MCPtAntiLambda;
953 fh1MCXiPt = o.fh1MCXiPt;
954 fh1MCXibarPt = o.fh1MCXibarPt;
955 fh2MCEtaVsPtK0s = o.fh2MCEtaVsPtK0s;
956 fh2MCEtaVsPtLa = o.fh2MCEtaVsPtLa;
957 fh2MCEtaVsPtALa = o.fh2MCEtaVsPtALa;
958 //fh1MCRapK0s = o.fh1MCRapK0s;
959 //fh1MCRapLambda = o.fh1MCRapLambda;
960 //fh1MCRapAntiLambda = o.fh1MCRapAntiLambda;
961 fh1MCEtaAllK0s = o.fh1MCEtaAllK0s;
962 fh1MCEtaK0s = o.fh1MCEtaK0s;
963 fh1MCEtaLambda = o.fh1MCEtaLambda;
964 fh1MCEtaAntiLambda = o.fh1MCEtaAntiLambda;
970 //_______________________________________________
971 AliAnalysisTaskJetChem::~AliAnalysisTaskJetChem()
976 if(fListK0s) delete fListK0s;
977 if(fListLa) delete fListLa;
978 if(fListALa) delete fListALa;
979 if(fListFeeddownLaCand) delete fListFeeddownLaCand;
980 if(fListFeeddownALaCand) delete fListFeeddownALaCand;
981 if(jetConeFDLalist) delete jetConeFDLalist;
982 if(jetConeFDALalist) delete jetConeFDALalist;
983 if(fListMCgenK0s) delete fListMCgenK0s;
984 if(fListMCgenLa) delete fListMCgenLa;
985 if(fListMCgenALa) delete fListMCgenALa;
986 if(fRandom) delete fRandom;
989 //________________________________________________________________________________________________________________________________
990 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AliFragFuncHistosInvMass(const char* name,
991 Int_t nJetPt, Float_t jetPtMin, Float_t jetPtMax,
992 Int_t nInvMass, Float_t invMassMin, Float_t invMassMax,
993 Int_t nPt, Float_t ptMin, Float_t ptMax,
994 Int_t nXi, Float_t xiMin, Float_t xiMax,
995 Int_t nZ , Float_t zMin , Float_t zMax )
1000 ,fNBinsInvMass(nInvMass)
1001 ,fInvMassMin(invMassMin)
1002 ,fInvMassMax(invMassMax)
1018 // default constructor
1022 //______________________________________________________________________________________________________________
1023 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AliFragFuncHistosInvMass(const AliFragFuncHistosInvMass& copy)
1025 ,fNBinsJetPt(copy.fNBinsJetPt)
1026 ,fJetPtMin(copy.fJetPtMin)
1027 ,fJetPtMax(copy.fJetPtMax)
1028 ,fNBinsInvMass(copy.fNBinsInvMass)
1029 ,fInvMassMin(copy.fInvMassMin)
1030 ,fInvMassMax(copy.fInvMassMax)
1031 ,fNBinsPt(copy.fNBinsPt)
1032 ,fPtMin(copy.fPtMin)
1033 ,fPtMax(copy.fPtMax)
1034 ,fNBinsXi(copy.fNBinsXi)
1035 ,fXiMin(copy.fXiMin)
1036 ,fXiMax(copy.fXiMax)
1037 ,fNBinsZ(copy.fNBinsZ)
1040 ,fh3TrackPt(copy.fh3TrackPt)
1043 ,fh1JetPt(copy.fh1JetPt)
1044 ,fNameFF(copy.fNameFF)
1049 //______________________________________________________________________________________________________________________________________________________________________
1050 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass& AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::operator=(const AliAnalysisTaskJetChem::AliFragFuncHistosInvMass& o)
1055 TObject::operator=(o);
1056 fNBinsJetPt = o.fNBinsJetPt;
1057 fJetPtMin = o.fJetPtMin;
1058 fJetPtMax = o.fJetPtMax;
1059 fNBinsInvMass = o.fNBinsInvMass;
1060 fInvMassMin = o.fInvMassMin;
1061 fInvMassMax = o.fInvMassMax;
1062 fNBinsPt = o.fNBinsPt;
1065 fNBinsXi = o.fNBinsXi;
1068 fNBinsZ = o.fNBinsZ;
1071 fh3TrackPt = o.fh3TrackPt;
1074 fh1JetPt = o.fh1JetPt;
1075 fNameFF = o.fNameFF;
1081 //___________________________________________________________________________
1082 AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::~AliFragFuncHistosInvMass()
1086 if(fh1JetPt) delete fh1JetPt;
1087 if(fh3TrackPt) delete fh3TrackPt;
1088 if(fh3Xi) delete fh3Xi;
1089 if(fh3Z) delete fh3Z;
1092 //_________________________________________________________________
1093 void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::DefineHistos()
1097 fh1JetPt = new TH1F(Form("fh1FFJetPtIM%s", fNameFF.Data()),"",fNBinsJetPt,fJetPtMin,fJetPtMax);
1098 fh3TrackPt = new TH3F(Form("fh3FFTrackPtIM%s",fNameFF.Data()),"",fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsPt, fPtMin, fPtMax);
1099 fh3Xi = new TH3F(Form("fh3FFXiIM%s", fNameFF.Data()),"", fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsXi, fXiMin, fXiMax);
1100 fh3Z = new TH3F(Form("fh3FFZIM%s", fNameFF.Data()),"", fNBinsJetPt, fJetPtMin, fJetPtMax, fNBinsInvMass, fInvMassMin, fInvMassMax, fNBinsZ, fZMin, fZMax);
1102 AliAnalysisTaskFragmentationFunction::SetProperties(fh1JetPt, "p_{t} (GeV/c)", "entries");
1103 AliAnalysisTaskJetChem::SetProperties(fh3TrackPt,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","p_{t} (GeV/c)");
1104 AliAnalysisTaskJetChem::SetProperties(fh3Xi,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","#xi");
1105 AliAnalysisTaskJetChem::SetProperties(fh3Z,"jet p_{t} (GeV/c)","inv Mass (GeV/c^2)","z");
1108 //________________________________________________________________________________________________________________________________
1109 void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::FillFF(Float_t trackPt, Float_t invM, Float_t jetPt, Bool_t incrementJetPt)
1111 // fill FF, don't use TH3F anymore use THnSparse instead to save memory
1113 if(incrementJetPt) fh1JetPt->Fill(jetPt);
1114 //fh3TrackPt->Fill(jetPt,invM,trackPt);//Fill(x,y,z)
1117 if(jetPt>0) z = trackPt / jetPt;
1119 if(z>0) xi = TMath::Log(1/z);
1121 //fh3Xi->Fill(jetPt,invM,xi);
1122 //fh3Z->Fill(jetPt,invM,z);
1125 //___________________________________________________________________________________
1126 void AliAnalysisTaskJetChem::AliFragFuncHistosInvMass::AddToOutput(TList* list) const
1128 // add histos to list
1130 list->Add(fh1JetPt);
1131 //list->Add(fh3TrackPt);
1137 //____________________________________________________
1138 void AliAnalysisTaskJetChem::UserCreateOutputObjects()
1140 // create output objects
1142 fRandom = new TRandom3(0);
1143 fRandom->SetSeed(0);
1145 if(fDebug > 1) Printf("AliAnalysisTaskJetChem::UserCreateOutputObjects()");
1147 // create list of tracks and jets
1149 fTracksRecCuts = new TList();
1150 fTracksRecCuts->SetOwner(kFALSE); //objects in TList wont be deleted when TList is deleted
1151 fJetsRecCuts = new TList();
1152 fJetsRecCuts->SetOwner(kFALSE);
1153 fBckgJetsRec = new TList();
1154 fBckgJetsRec->SetOwner(kFALSE);
1155 fListK0s = new TList();
1156 fListK0s->SetOwner(kFALSE);
1157 fListLa = new TList();
1158 fListLa->SetOwner(kFALSE);
1159 fListALa = new TList();
1160 fListALa->SetOwner(kFALSE);
1161 fListFeeddownLaCand = new TList(); //feeddown Lambda candidates
1162 fListFeeddownLaCand->SetOwner(kFALSE);
1163 fListFeeddownALaCand = new TList(); //feeddown Antilambda candidates
1164 fListFeeddownALaCand->SetOwner(kFALSE);
1165 jetConeFDLalist = new TList();
1166 jetConeFDLalist->SetOwner(kFALSE); //feeddown Lambda candidates in jet cone
1167 jetConeFDALalist = new TList();
1168 jetConeFDALalist->SetOwner(kFALSE); //feeddown Antilambda candidates in jet cone
1169 fListMCgenK0s = new TList(); //MC generated K0s
1170 fListMCgenK0s->SetOwner(kFALSE);
1171 fListMCgenLa = new TList(); //MC generated Lambdas
1172 fListMCgenLa->SetOwner(kFALSE);
1173 fListMCgenALa = new TList(); //MC generated Antilambdas
1174 fListMCgenALa->SetOwner(kFALSE);
1177 // Create histograms / output container
1179 fCommonHistList = new TList();
1180 fCommonHistList->SetOwner();
1182 Bool_t oldStatus = TH1::AddDirectoryStatus();
1183 TH1::AddDirectory(kFALSE);//By default (fAddDirectory = kTRUE), histograms are automatically added to the list of objects in memory
1185 // histograms inherited from AliAnalysisTaskFragmentationFunction
1187 fh1EvtSelection = new TH1F("fh1EvtSelection", "Event Selection", 6, -0.5, 5.5);
1188 fh1EvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
1189 fh1EvtSelection->GetXaxis()->SetBinLabel(2,"event trigger selection: rejected");
1190 fh1EvtSelection->GetXaxis()->SetBinLabel(3,"event class: rejected");
1191 fh1EvtSelection->GetXaxis()->SetBinLabel(4,"vertex Ncontr: rejected");
1192 fh1EvtSelection->GetXaxis()->SetBinLabel(5,"vertex z: rejected");
1193 fh1EvtSelection->GetXaxis()->SetBinLabel(6,"vertex type: rejected");
1196 fh1EvtCent = new TH1F("fh1EvtCent","centrality",100,0.,100.);
1197 fh1VertexNContributors = new TH1F("fh1VertexNContributors", "Vertex N contributors", 11,-.5, 10.5);
1198 fh1VertexZ = new TH1F("fh1VertexZ", "Vertex z distribution", 30, -15., 15.);
1199 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
1200 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
1201 fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
1202 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
1203 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
1204 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
1205 fh1nRecJetsCuts = new TH1F("fh1nRecJetsCuts","reconstructed jets per event",100,-0.5,99.5);
1207 // histograms JetChem task
1209 fh1EvtAllCent = new TH1F("fh1EvtAllCent","before centrality selection",100,0.,100.);
1210 fh1Evt = new TH1F("fh1Evt", "All events runned over", 3, 0.,1.);
1211 fh1EvtMult = new TH1F("fh1EvtMult","multiplicity",240,0.,240.);
1212 fh1K0Mult = new TH1F("fh1K0Mult","K0 multiplicity",100,0.,100.);//500. all
1213 fh1dPhiJetK0 = new TH1F("fh1dPhiJetK0","",64,-1,5.4);
1214 fh1LaMult = new TH1F("fh1LaMult","La multiplicity",100,0.,100.);
1215 fh1dPhiJetLa = new TH1F("fh1dPhiJetLa","",64,-1,5.4);
1216 fh1ALaMult = new TH1F("fh1ALaMult","ALa multiplicity",100,0.,100.);
1217 fh1dPhiJetALa = new TH1F("fh1dPhiJetALa","",64,-1,5.4);
1218 fh1JetEta = new TH1F("fh1JetEta","#eta distribution of all jets",40,-2.,2.);
1219 fh1JetPhi = new TH1F("fh1JetPhi","#phi distribution of all jets",63,0.,6.3);
1220 fh2JetEtaPhi = new TH2F("fh2JetEtaPhi","#eta and #phi distribution of all jets",400,-2.,2.,63,0.,6.3);
1223 //fh1V0JetPt = new TH1F("fh1V0JetPt","p_{T} distribution of all jets containing v0s",200,0.,200.);
1224 fh1IMK0Cone = new TH1F("fh1IMK0Cone","p_{T} distribution of all jets containing K0s candidates",19,5.,100.);
1225 fh1IMLaCone = new TH1F("fh1IMLaCone","p_{T} distribution of all jets containing #Lambda candidates",19,5.,100.);
1226 fh1IMALaCone = new TH1F("fh1IMALaCone","p_{T} distribution of all jets containing #bar{#Lambda} candidates",19,5.,100.);
1228 fh2FFJetTrackEta = new TH2F("fh2FFJetTrackEta","charged track eta distr. in jet cone",200,-1.,1.,40,0.,200.);
1229 //fh1trackPosNCls = new TH1F("fh1trackPosNCls","NTPC clusters positive daughters",10,0.,100.);
1230 //fh1trackNegNCls = new TH1F("fh1trackNegNCls","NTPC clusters negative daughters",10,0.,100.);
1231 fh1trackPosEta = new TH1F("fh1trackPosEta","eta positive daughters",100,-2.,2.);
1232 fh1trackNegEta = new TH1F("fh1trackNegEta","eta negative daughters",100,-2.,2.);
1233 fh1V0Eta = new TH1F("fh1V0Eta","V0 eta",60,-1.5,1.5);
1234 //fh1V0totMom = new TH1F("fh1V0totMom","V0 tot mom",100,0.,20.);
1235 fh1CosPointAngle = new TH1F("fh1CosPointAngle", "Cosine of V0's pointing angle",50,0.99,1.0);
1236 fh1DecayLengthV0 = new TH1F("fh1DecayLengthV0", "V0s decay Length;decay length(cm)",1200,0.,120.);
1237 fh2ProperLifetimeK0sVsPtBeforeCut = new TH2F("fh2ProperLifetimeK0sVsPtBeforeCut"," K0s ProperLifetime vs Pt; p_{T} (GeV/#it{c})",150,0.,15.,250,0.,250.);
1238 fh2ProperLifetimeK0sVsPtAfterCut = new TH2F("fh2ProperLifetimeK0sVsPtAfterCut"," K0s ProperLifetime vs Pt; p_{T} (GeV/#it{c})",150,0.,15.,250,0.,250.);
1239 fh1V0Radius = new TH1F("fh1V0Radius", "V0s Radius;Radius(cm)",200,0.,40.);
1240 fh1DcaV0Daughters = new TH1F("fh1DcaV0Daughters", "DCA between daughters;dca(cm)",200,0.,2.);
1241 fh1DcaPosToPrimVertex = new TH1F("fh1DcaPosToPrimVertex", "Positive V0 daughter;dca(cm)",100,0.,10.);
1242 fh1DcaNegToPrimVertex = new TH1F("fh1DcaNegToPrimVertex", "Negative V0 daughter;dca(cm)",100,0.,10.);
1243 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);
1244 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);
1245 fh2BBLaPos = new TH2F("fh2BBLaPos","PID of the positive daughter of La candidates; P (GeV); -dE/dx (keV/cm ?)",100,0,10,200,0,200);
1246 fh2BBLaNeg = new TH2F("fh2BBLaNeg","PID of the negative daughter of La candidates; P (GeV); -dE/dx (keV/cm ?)",100,0,10,200,0,200);
1247 fh1PosDaughterCharge = new TH1F("fh1PosDaughterCharge","charge of V0 positive daughters; V0 daughters",3,-2.,2.);
1248 fh1NegDaughterCharge = new TH1F("fh1NegDaughterCharge","charge of V0 negative daughters; V0 daughters",3,-2.,2.);
1249 fh1PtMCK0s = new TH1F("fh1PtMCK0s","Pt of MC rec K0s; #it{p}_{T} (GeV/#it{c})",200,0.,20.);
1250 fh1PtMCLa = new TH1F("fh1PtMCLa","Pt of MC rec La; #it{p}_{T} (GeV/#it{c})",200,0.,20.);
1251 fh1PtMCALa = new TH1F("fh1PtMCALa","Pt of MC rec ALa; #it{p}_{T} (GeV/#it{c})",200,0.,20.);
1252 fh1EtaK0s = new TH1F("fh1EtaK0s","K^{0}_{s} entries ;#eta",200,-1.,1.);
1253 fh1EtaLa = new TH1F("fh1EtaLa","#Lambda entries ;#eta",200,-1.,1.);
1254 fh1EtaALa = new TH1F("fh1EtaALa","#bar{#Lambda} entries ;#eta",200,-1.,1.);
1256 Int_t binsInvMassEtaTrackPtK0s[3] = {200, 200, 120};//eta,invM,trackPt
1257 Double_t xminInvMassEtaTrackPtK0s[3] = {-1.,0.3,0.};
1258 Double_t xmaxInvMassEtaTrackPtK0s[3] = {1.,0.7,12.};
1260 fhnInvMassEtaTrackPtK0s = new THnSparseF("fhnInvMassEtaTrackPtK0s","#eta; K0s invM (GeV/{#it{c}}^{2}); #it{p}_{T} (GeV/#it{c})",3,binsInvMassEtaTrackPtK0s,xminInvMassEtaTrackPtK0s,xmaxInvMassEtaTrackPtK0s);
1262 Int_t binsInvMassEtaTrackPtLa[3] = {200, 200, 120};//eta,invM,trackPt
1263 Double_t xminInvMassEtaTrackPtLa[3] = {-1.,1.05,0.};
1264 Double_t xmaxInvMassEtaTrackPtLa[3] = {1.,1.25,12.};
1266 fhnInvMassEtaTrackPtLa = new THnSparseF("fhnInvMassEtaTrackPtLa","#eta; #Lambda invM (GeV/{#it{c}}^{2}); #it{p}_{T} (GeV/#it{c})",3,binsInvMassEtaTrackPtLa,xminInvMassEtaTrackPtLa,xmaxInvMassEtaTrackPtLa);
1268 Int_t binsInvMassEtaTrackPtALa[3] = {200, 200, 120};//eta,invM,trackPt
1269 Double_t xminInvMassEtaTrackPtALa[3] = {-1.,1.05,0.};
1270 Double_t xmaxInvMassEtaTrackPtALa[3] = {1.,1.25,12.};
1272 fhnInvMassEtaTrackPtALa = new THnSparseF("fhnInvMassEtaTrackPtALa","#eta; #bar{#Lambda} invM (GeV/{#it{c}}^{2}); #it{p}_{T} (GeV/#it{c})",3,binsInvMassEtaTrackPtALa,xminInvMassEtaTrackPtALa,xmaxInvMassEtaTrackPtALa);
1274 Int_t binsK0sPC[4] = {19, 200, 200, 200};
1275 Double_t xminK0sPC[4] = {5.,0.3, 0., -1.};
1276 Double_t xmaxK0sPC[4] = {100.,0.7, 20., 1.};
1277 fhnK0sPC = new THnSparseF("fhnK0sPC","jet pT; K0s invM; particle pT; particle #eta",4,binsK0sPC,xminK0sPC,xmaxK0sPC);
1279 Int_t binsLaPC[4] = {19, 200, 200, 200};
1280 Double_t xminLaPC[4] = {5.,1.05, 0., -1.};
1281 Double_t xmaxLaPC[4] = {100.,1.25, 20., 1.};
1282 fhnLaPC = new THnSparseF("fhnLaPC","jet pT; #Lambda invM; particle pT; particle #eta",4,binsLaPC,xminLaPC,xmaxLaPC);
1284 Int_t binsALaPC[4] = {19, 200, 200, 200};
1285 Double_t xminALaPC[4] = {5.,1.05, 0., -1.};
1286 Double_t xmaxALaPC[4] = {100.,1.25, 20., 1.};
1287 fhnALaPC = new THnSparseF("fhnALaPC","jet pT; #bar#Lambda invM; particle pT; particle #eta",4,binsALaPC,xminALaPC,xmaxALaPC);
1289 Int_t binsK0sMCC[3] = {200, 200, 200};
1290 Double_t xminK0sMCC[3] = {0.3, 0., -1.};
1291 Double_t xmaxK0sMCC[3] = {0.7, 20., 1.};
1292 fhnK0sMCC = new THnSparseF("fhnK0sMCC","jet pT; K0s invM; particle pT; particle #eta",3,binsK0sMCC,xminK0sMCC,xmaxK0sMCC);
1294 Int_t binsLaMCC[3] = {200, 200, 200};
1295 Double_t xminLaMCC[3] = {1.05, 0., -1.};
1296 Double_t xmaxLaMCC[3] = {1.25, 20., 1.};
1297 fhnLaMCC = new THnSparseF("fhnLaMCC","jet pT; #Lambda invM; particle pT; particle #eta",3,binsLaMCC,xminLaMCC,xmaxLaMCC);
1299 Int_t binsALaMCC[3] = {200, 200, 200};
1300 Double_t xminALaMCC[3] = {1.05, 0., -1.};
1301 Double_t xmaxALaMCC[3] = {1.25, 20., 1.};
1302 fhnALaMCC = new THnSparseF("fhnALaMCC","jet pT; #bara#Lambda invM; particle pT; particle #eta",3,binsALaMCC,xminALaMCC,xmaxALaMCC);
1304 Int_t binsK0sRC[3] = {200, 200, 200};
1305 Double_t xminK0sRC[3] = {0.3, 0., -1.};
1306 Double_t xmaxK0sRC[3] = {0.7, 20., 1.};
1307 fhnK0sRC = new THnSparseF("fhnK0sRC","jet pT; K0s invM; particle pT; particle #eta",3,binsK0sRC,xminK0sRC,xmaxK0sRC);
1309 Int_t binsLaRC[3] = {200, 200, 200};
1310 Double_t xminLaRC[3] = {1.05, 0., -1.};
1311 Double_t xmaxLaRC[3] = {1.25, 20., 1.};
1312 fhnLaRC = new THnSparseF("fhnLaRC","jet pT; #Lambda invM; particle pT; particle #eta",3,binsLaRC,xminLaRC,xmaxLaRC);
1314 Int_t binsALaRC[3] = {200, 200, 200};
1315 Double_t xminALaRC[3] = {1.05, 0., -1.};
1316 Double_t xmaxALaRC[3] = {1.25, 20., 1.};
1317 fhnALaRC = new THnSparseF("fhnALaRC","jet pT; #bara#Lambda invM; particle pT; particle #eta",3,binsALaRC,xminALaRC,xmaxALaRC);
1320 Int_t binsK0sOC[3] = {200, 200, 200};
1321 Double_t xminK0sOC[3] = {0.3, 0., -1.};
1322 Double_t xmaxK0sOC[3] = {0.7, 20., 1.};
1323 fhnK0sOC = new THnSparseF("fhnK0sOC","jet pT; K0s invM; particle pT; particle #eta",3,binsK0sOC,xminK0sOC,xmaxK0sOC);
1325 Int_t binsLaOC[3] = {200, 200, 200};
1326 Double_t xminLaOC[3] = {1.05, 0., -1.};
1327 Double_t xmaxLaOC[3] = {1.25, 20., 1.};
1328 fhnLaOC = new THnSparseF("fhnLaOC","jet pT; #Lambda invM; particle pT; particle #eta",3,binsLaOC,xminLaOC,xmaxLaOC);
1330 Int_t binsALaOC[3] = {200, 200, 200};
1331 Double_t xminALaOC[3] = {1.05, 0., -1.};
1332 Double_t xmaxALaOC[3] = {1.25, 20., 1.};
1334 fhnALaOC = new THnSparseF("fhnALaOC","jet pT; #bara#Lambda invM; particle pT; particle #eta",3,binsALaOC,xminALaOC,xmaxALaOC);
1340 fh1AreaExcluded = new TH1F("fh1AreaExcluded","area excluded for selected jets in event acceptance",100,0.,5.);
1342 fh1MedianEta = new TH1F("fh1MedianEta","Median cluster axis ;#eta",200,-1.,1.);
1343 fh1JetPtMedian = new TH1F("fh1JetPtMedian"," (selected) jet it{p}_{T} distribution for MCC method; #GeV/it{c}",19,5.,100.);
1345 fh1TrackMultCone = new TH1F("fh1TrackMultCone","track multiplicity in jet cone; number of tracks",20,0.,50.);
1347 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.);
1349 Int_t binsNJK0[3] = {200, 200, 200};
1350 Double_t xminNJK0[3] = {0.3, 0., -1.};
1351 Double_t xmaxNJK0[3] = {0.7, 20., 1.};
1352 fhnNJK0 = new THnSparseF("fhnNJK0","K0s candidates in events wo selected jets;",3,binsNJK0,xminNJK0,xmaxNJK0);
1354 Int_t binsNJLa[3] = {200, 200, 200};
1355 Double_t xminNJLa[3] = {1.05, 0., -1.};
1356 Double_t xmaxNJLa[3] = {1.25, 20., 1.};
1357 fhnNJLa = new THnSparseF("fhnNJLa","La candidates in events wo selected jets; ",3,binsNJLa,xminNJLa,xmaxNJLa);
1359 Int_t binsNJALa[3] = {200, 200, 200};
1360 Double_t xminNJALa[3] = {1.05, 0., -1.};
1361 Double_t xmaxNJALa[3] = {1.25, 20., 1.};
1362 fhnNJALa = new THnSparseF("fhnNJALa","ALa candidates in events wo selected jets; ",3,binsNJALa,xminNJALa,xmaxNJALa);
1364 fFFHistosRecCuts = new AliFragFuncHistos("RecCuts", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1365 fFFNBinsPt, fFFPtMin, fFFPtMax,
1366 fFFNBinsXi, fFFXiMin, fFFXiMax,
1367 fFFNBinsZ , fFFZMin , fFFZMax);
1369 fV0QAK0 = new AliFragFuncQATrackHistos("V0QAK0",fQATrackNBinsPt, fQATrackPtMin, fQATrackPtMax,
1370 fQATrackNBinsEta, fQATrackEtaMin, fQATrackEtaMax,
1371 fQATrackNBinsPhi, fQATrackPhiMin, fQATrackPhiMax,
1372 fQATrackHighPtThreshold);
1374 fFFHistosRecCutsK0Evt = new AliFragFuncHistos("RecCutsK0Evt", fFFNBinsJetPt, fFFJetPtMin, fFFJetPtMax,
1375 fFFNBinsPt, fFFPtMin, fFFPtMax,
1376 fFFNBinsXi, fFFXiMin, fFFXiMax,
1377 fFFNBinsZ , fFFZMin , fFFZMax);
1380 fFFHistosIMK0AllEvt = new AliFragFuncHistosInvMass("K0AllEvt", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax,
1381 fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
1382 fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax,
1383 fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,
1384 fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
1386 fFFHistosIMK0Jet = new AliFragFuncHistosInvMass("K0Jet", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax,
1387 fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
1388 fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax,
1389 fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,
1390 fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
1392 fFFHistosIMK0Cone = new AliFragFuncHistosInvMass("K0Cone", fFFIMNBinsJetPt, fFFIMJetPtMin, fFFIMJetPtMax,
1393 fFFIMNBinsInvM,fFFIMInvMMin,fFFIMInvMMax,
1394 fFFIMNBinsPt, fFFIMPtMin, fFFIMPtMax,
1395 fFFIMNBinsXi, fFFIMXiMin, fFFIMXiMax,
1396 fFFIMNBinsZ , fFFIMZMin , fFFIMZMax);
1398 fFFHistosIMLaAllEvt = new AliFragFuncHistosInvMass("LaAllEvt", fFFIMLaNBinsJetPt, fFFIMLaJetPtMin, fFFIMLaJetPtMax,
1399 fFFIMLaNBinsInvM,fFFIMLaInvMMin,fFFIMLaInvMMax,
1400 fFFIMLaNBinsPt, fFFIMLaPtMin, fFFIMLaPtMax,
1401 fFFIMLaNBinsXi, fFFIMLaXiMin, fFFIMLaXiMax,
1402 fFFIMLaNBinsZ , fFFIMLaZMin , fFFIMLaZMax);
1404 fFFHistosIMLaJet = new AliFragFuncHistosInvMass("LaJet", fFFIMLaNBinsJetPt, fFFIMLaJetPtMin, fFFIMLaJetPtMax,
1405 fFFIMLaNBinsInvM,fFFIMLaInvMMin,fFFIMLaInvMMax,
1406 fFFIMLaNBinsPt, fFFIMLaPtMin, fFFIMLaPtMax,
1407 fFFIMLaNBinsXi, fFFIMLaXiMin, fFFIMLaXiMax,
1408 fFFIMLaNBinsZ , fFFIMLaZMin , fFFIMLaZMax);
1411 fFFHistosIMLaCone = new AliFragFuncHistosInvMass("LaCone", fFFIMLaNBinsJetPt, fFFIMLaJetPtMin, fFFIMLaJetPtMax,
1412 fFFIMLaNBinsInvM,fFFIMLaInvMMin,fFFIMLaInvMMax,
1413 fFFIMLaNBinsPt, fFFIMLaPtMin, fFFIMLaPtMax,
1414 fFFIMLaNBinsXi, fFFIMLaXiMin, fFFIMLaXiMax,
1415 fFFIMLaNBinsZ , fFFIMLaZMin , fFFIMLaZMax);
1418 fFFHistosIMALaAllEvt = new AliFragFuncHistosInvMass("ALaAllEvt", fFFIMLaNBinsJetPt, fFFIMLaJetPtMin, fFFIMLaJetPtMax,
1419 fFFIMLaNBinsInvM,fFFIMLaInvMMin,fFFIMLaInvMMax,
1420 fFFIMLaNBinsPt, fFFIMLaPtMin, fFFIMLaPtMax,
1421 fFFIMLaNBinsXi, fFFIMLaXiMin, fFFIMLaXiMax,
1422 fFFIMLaNBinsZ , fFFIMLaZMin , fFFIMLaZMax);
1424 fFFHistosIMALaJet = new AliFragFuncHistosInvMass("ALaJet", fFFIMLaNBinsJetPt, fFFIMLaJetPtMin, fFFIMLaJetPtMax,
1425 fFFIMLaNBinsInvM,fFFIMLaInvMMin,fFFIMLaInvMMax,
1426 fFFIMLaNBinsPt, fFFIMLaPtMin, fFFIMLaPtMax,
1427 fFFIMLaNBinsXi, fFFIMLaXiMin, fFFIMLaXiMax,
1428 fFFIMLaNBinsZ , fFFIMLaZMin , fFFIMLaZMax);
1430 fFFHistosIMALaCone = new AliFragFuncHistosInvMass("ALaCone", fFFIMLaNBinsJetPt, fFFIMLaJetPtMin, fFFIMLaJetPtMax,
1431 fFFIMLaNBinsInvM,fFFIMLaInvMMin,fFFIMLaInvMMax,
1432 fFFIMLaNBinsPt, fFFIMLaPtMin, fFFIMLaPtMax,
1433 fFFIMLaNBinsXi, fFFIMLaXiMin, fFFIMLaXiMax,
1434 fFFIMLaNBinsZ , fFFIMLaZMin , fFFIMLaZMax);
1441 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.);
1442 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.);
1443 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.);
1445 fh2MCgenK0Cone->GetYaxis()->SetTitle("MC gen K^{0}}^{s} #it{p}_{T}");
1446 fh2MCgenLaCone->GetYaxis()->SetTitle("MC gen #Lambda #it{p}_{T}");
1447 fh2MCgenALaCone->GetYaxis()->SetTitle("MC gen #Antilambda #it{p}_{T}");
1449 fh2MCEtagenK0Cone = new TH2F("fh2MCEtagenK0Cone","MC gen {K^{0}}^{s} #it{p}_{T} #eta distribution in jet cone;#eta",19,5.,100.,200,-1.,1.);
1450 fh2MCEtagenLaCone = new TH2F("fh2MCEtagenLaCone","MC gen #Lambda #it{p}_{T} #eta distribution in jet cone;#eta",19,5.,100.,200,-1.,1.);
1451 fh2MCEtagenALaCone = new TH2F("fh2MCEtagenALaCone","MC gen #Antilambda #it{p}_{T} #eta distribution in jet cone;#eta",19,5.,100.,200,-1.,1.);
1452 fh1IMK0ConeSmear = new TH1F("fh1IMK0ConeSmear","Smeared jet pt study for K0s-in-cone-jets; smeared jet #it{p}_{T}", 19,5.,100.);
1453 fh1IMLaConeSmear = new TH1F("fh1IMLaConeSmear","Smeared jet pt study for La-in-cone-jets; smeared jet #it{p}_{T}", 19,5.,100.);
1454 fh1IMALaConeSmear = new TH1F("fh1IMALaConeSmear","Smeared jet pt study for ALa-in-cone-jets; smeared jet #it{p}_{T}", 19,5.,100.);
1456 //8 new histograms: Cone, Incl, Lambda, Antilambda, Hijing, Injected:
1458 Int_t binsrecMCHijingLaIncl[3] = {200, 200, 200};
1459 Double_t xminrecMCHijingLaIncl[3] = {1.05, 0., -1.};
1460 Double_t xmaxrecMCHijingLaIncl[3] = {1.25, 20., 1.};
1461 fhnrecMCHijingLaIncl = new THnSparseF("fhnrecMCHijingLaIncl","La inv. mass; particle pT; particle #eta",3,binsrecMCHijingLaIncl,xminrecMCHijingLaIncl,xmaxrecMCHijingLaIncl);
1463 Int_t binsrecMCHijingLaCone[4] = {19, 200, 200, 200};
1464 Double_t xminrecMCHijingLaCone[4] = {5., 1.05, 0., -1.};
1465 Double_t xmaxrecMCHijingLaCone[4] = {100., 1.25, 20., 1.};
1466 fhnrecMCHijingLaCone = new THnSparseF("fhnrecMCHijingLaCone","La inv. mass; particle pT; particle #eta",4,binsrecMCHijingLaCone,xminrecMCHijingLaCone,xmaxrecMCHijingLaCone);
1468 Int_t binsrecMCHijingALaIncl[3] = {200, 200, 200};
1469 Double_t xminrecMCHijingALaIncl[3] = {1.05, 0., -1.};
1470 Double_t xmaxrecMCHijingALaIncl[3] = {1.25, 20., 1.};
1471 fhnrecMCHijingALaIncl = new THnSparseF("fhnrecMCHijingALaIncl","ALa inv. mass; particle pT; particle #eta",3,binsrecMCHijingALaIncl,xminrecMCHijingALaIncl,xmaxrecMCHijingALaIncl);
1473 Int_t binsrecMCHijingALaCone[4] = {19, 200, 200, 200};
1474 Double_t xminrecMCHijingALaCone[4] = {5., 1.05, 0., -1.};
1475 Double_t xmaxrecMCHijingALaCone[4] = {100., 1.25, 20., 1.};
1476 fhnrecMCHijingALaCone = new THnSparseF("fhnrecMCHijingALaCone","ALa inv. mass; particle pT; particle #eta",4,binsrecMCHijingALaCone,xminrecMCHijingALaCone,xmaxrecMCHijingALaCone);
1478 Int_t binsrecMCInjectLaIncl[3] = {200, 200, 200};
1479 Double_t xminrecMCInjectLaIncl[3] = {1.05, 0., -1.};
1480 Double_t xmaxrecMCInjectLaIncl[3] = {1.25, 20., 1.};
1481 fhnrecMCInjectLaIncl = new THnSparseF("fhnrecMCInjectLaIncl","La inv. mass; particle pT; particle #eta",3,binsrecMCInjectLaIncl,xminrecMCInjectLaIncl,xmaxrecMCInjectLaIncl);
1483 Int_t binsrecMCInjectLaCone[4] = {19, 200, 200, 200};
1484 Double_t xminrecMCInjectLaCone[4] = {5., 1.05, 0., -1.};
1485 Double_t xmaxrecMCInjectLaCone[4] = {100., 1.25, 20., 1.};
1486 fhnrecMCInjectLaCone = new THnSparseF("fhnrecMCInjectLaCone","La jet pT;inv. mass; particle pT; particle #eta",4,binsrecMCInjectLaCone,xminrecMCInjectLaCone,xmaxrecMCInjectLaCone);
1488 Int_t binsrecMCInjectALaIncl[3] = {200, 200, 200};
1489 Double_t xminrecMCInjectALaIncl[3] = {1.05, 0., -1.};
1490 Double_t xmaxrecMCInjectALaIncl[3] = {1.25, 20., 1.};
1491 fhnrecMCInjectALaIncl = new THnSparseF("fhnrecMCInjectALaIncl","ALa inv. mass; particle pT; particle #eta",3,binsrecMCInjectALaIncl,xminrecMCInjectALaIncl,xmaxrecMCInjectALaIncl);
1493 Int_t binsrecMCInjectALaCone[4] = {19, 200, 200, 200};
1494 Double_t xminrecMCInjectALaCone[4] = {5., 1.05, 0., -1.};
1495 Double_t xmaxrecMCInjectALaCone[4] = {100., 1.25, 20., 1.};
1496 fhnrecMCInjectALaCone = new THnSparseF("fhnrecMCInjectALaCone","ALa inv. mass; particle pT; particle #eta",4,binsrecMCInjectALaCone,xminrecMCInjectALaCone,xmaxrecMCInjectALaCone);
1499 Int_t binsMCrecK0Cone[4] = {19, 200, 200, 200};
1500 Double_t xminMCrecK0Cone[4] = {5.,0.3, 0., -1.};
1501 Double_t xmaxMCrecK0Cone[4] = {100.,0.7, 20., 1.};
1502 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);
1504 Int_t binsMCrecLaCone[4] = {19, 200, 200, 200};
1505 Double_t xminMCrecLaCone[4] = {5.,0.3, 0., -1.};
1506 Double_t xmaxMCrecLaCone[4] = {100.,0.7, 20., 1.};
1507 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);
1509 Int_t binsMCrecALaCone[4] = {19, 200, 200, 200};
1510 Double_t xminMCrecALaCone[4] = {5.,0.3, 0., -1.};
1511 Double_t xmaxMCrecALaCone[4] = {100.,0.7, 20., 1.};
1512 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);
1514 Int_t binsMCrecK0ConeSmear[4] = {19, 200, 200, 200};
1515 Double_t xminMCrecK0ConeSmear[4] = {5.,0.3, 0., -1.};
1516 Double_t xmaxMCrecK0ConeSmear[4] = {100.,0.7, 20., 1.};
1517 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);
1519 Int_t binsMCrecLaConeSmear[4] = {19, 200, 200, 200};
1520 Double_t xminMCrecLaConeSmear[4] = {5.,0.3, 0., -1.};
1521 Double_t xmaxMCrecLaConeSmear[4] = {100.,0.7, 20., 1.};
1522 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);
1524 Int_t binsMCrecALaConeSmear[4] = {19, 200, 200, 200};
1525 Double_t xminMCrecALaConeSmear[4] = {5.,0.3, 0., -1.};
1526 Double_t xmaxMCrecALaConeSmear[4] = {100.,0.7, 20., 1.};
1527 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);
1529 Int_t binsK0sSecContinCone[3] = {19, 200, 200};
1530 Double_t xminK0sSecContinCone[3] = {5.,0., -1.};
1531 Double_t xmaxK0sSecContinCone[3] = {100.,20., 1.};
1532 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);
1534 Int_t binsLaSecContinCone[3] = {19, 200, 200};
1535 Double_t xminLaSecContinCone[3] = {5.,0., -1.};
1536 Double_t xmaxLaSecContinCone[3] = {100.,20., 1.};
1537 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);
1539 Int_t binsALaSecContinCone[3] = {19, 200, 200};
1540 Double_t xminALaSecContinCone[3] = {5.,0., -1.};
1541 Double_t xmaxALaSecContinCone[3] = {100.,20., 1.};
1542 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);
1544 Int_t binsK0sIncl[3] = {200, 200, 200};
1545 Double_t xminK0sIncl[3] = {0.3, 0., -1.};
1546 Double_t xmaxK0sIncl[3] = {0.7, 20., 1.};
1547 fhnK0sIncl = new THnSparseF("fhnK0sIncl","K0s inv. mass; particle pT; particle #eta",3,binsK0sIncl,xminK0sIncl,xmaxK0sIncl);
1549 Int_t binsK0sCone[4] = {19, 200, 200, 200};
1550 Double_t xminK0sCone[4] = {5.,0.3, 0., -1.};
1551 Double_t xmaxK0sCone[4] = {100.,0.7, 20., 1.};
1552 fhnK0sCone = new THnSparseF("fhnK0sCone","jet pT; K0s inv. mass; particle pT; particle #eta",4,binsK0sCone,xminK0sCone,xmaxK0sCone);
1554 Int_t binsLaIncl[3] = {200, 200, 200};
1555 Double_t xminLaIncl[3] = {1.05, 0., -1.};
1556 Double_t xmaxLaIncl[3] = {1.25, 20., 1.};
1557 fhnLaIncl = new THnSparseF("fhnLaIncl","La inv. mass; particle pT; particle #eta",3,binsLaIncl,xminLaIncl,xmaxLaIncl);
1559 Int_t binsLaCone[4] = {19, 200, 200, 200};
1560 Double_t xminLaCone[4] = {5.,1.05, 0., -1.};
1561 Double_t xmaxLaCone[4] = {100.,1.25, 20., 1.};
1562 fhnLaCone = new THnSparseF("fhnLaCone","jet pT; La inv. mass; particle pT; particle #eta",4,binsLaCone,xminLaCone,xmaxLaCone);
1564 Int_t binsALaIncl[3] = {200, 200, 200};
1565 Double_t xminALaIncl[3] = {1.05, 0., -1.};
1566 Double_t xmaxALaIncl[3] = {1.25, 20., 1.};
1567 fhnALaIncl = new THnSparseF("fhnALaIncl","ALa inv. mass; particle pT; particle #eta",3,binsALaIncl,xminALaIncl,xmaxALaIncl);
1569 Int_t binsALaCone[4] = {19, 200, 200, 200};
1570 Double_t xminALaCone[4] = {5.,1.05, 0., -1.};
1571 Double_t xmaxALaCone[4] = {100.,1.25, 20., 1.};
1572 fhnALaCone = new THnSparseF("fhnALaCone","jet pT; ALa inv. mass; particle pT; particle #eta",4,binsALaCone,xminALaCone,xmaxALaCone);
1574 fh1MCMultiplicityPrimary = new TH1F("fh1MCMultiplicityPrimary", "MC Primary Particles;NPrimary;Count", 201, -0.5, 200.5);
1575 fh1MCMultiplicityTracks = new TH1F("h1MCMultiplicityTracks", "MC Tracks;Ntracks;Count", 201, -0.5, 200.5);
1578 Int_t binsFeedDownLa[3] = {19, 200, 200};
1579 Double_t xminFeedDownLa[3] = {5.,1.05, 0.};
1580 Double_t xmaxFeedDownLa[3] = {100.,1.25, 20.};
1581 fhnFeedDownLa = new THnSparseF("fhnFeedDownLa","#Lambda stemming from feeddown from Xi(0/-)",3,binsFeedDownLa,xminFeedDownLa,xmaxFeedDownLa);
1583 Int_t binsFeedDownALa[3] = {19, 200, 200};
1584 Double_t xminFeedDownALa[3] = {5.,1.05, 0.};
1585 Double_t xmaxFeedDownALa[3] = {100.,1.25, 20.};
1586 fhnFeedDownALa = new THnSparseF("fhnFeedDownALa","#bar#Lambda stemming from feeddown from Xibar(0/+)",3,binsFeedDownALa,xminFeedDownALa,xmaxFeedDownALa);
1588 Int_t binsFeedDownLaCone[3] = {19, 200, 200};
1589 Double_t xminFeedDownLaCone[3] = {5.,1.05, 0.};
1590 Double_t xmaxFeedDownLaCone[3] = {100.,1.25, 20.};
1591 fhnFeedDownLaCone = new THnSparseF("fhnFeedDownLaCone","#Lambda stemming from feeddown from Xi(0/-) in jet cone",3,binsFeedDownLaCone,xminFeedDownLaCone,xmaxFeedDownLaCone);
1593 Int_t binsFeedDownALaCone[3] = {19, 200, 200};
1594 Double_t xminFeedDownALaCone[3] = {5.,1.05, 0.};
1595 Double_t xmaxFeedDownALaCone[3] = {100.,1.25, 20.};
1596 fhnFeedDownALaCone = new THnSparseF("fhnFeedDownALaCone","#bar#Lambda stemming from feeddown from Xibar(0/+) in jet cone",3,binsFeedDownALaCone,xminFeedDownALaCone,xmaxFeedDownALaCone);
1599 fh1MCProdRadiusK0s = new TH1F("fh1MCProdRadiusK0s","MC gen. MC K0s prod radius",200,0.,200.);
1600 fh1MCProdRadiusLambda = new TH1F("fh1MCProdRadiusLambda","MC gen. MC La prod radius",200,0.,200.);
1601 fh1MCProdRadiusAntiLambda = new TH1F("fh1MCProdRadiusAntiLambda","MC gen. MC ALa prod radius",200,0.,200.);
1603 // Pt and inv mass distributions
1605 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
1606 fh1MCPtK0s = new TH1F("fh1MCPtK0s", "MC gen. K^{0}_{s} in eta range;#it{p}_{T} (GeV/#it{c})",200,0.,20.);
1607 fh1MCPtLambda = new TH1F("fh1MCPtLambda", "MC gen. #Lambda in rap range;#it{p}_{T} (GeV/#it{c})",200,0.,20.);
1608 fh1MCPtAntiLambda = new TH1F("fh1MCPtAntiLambda", "MC gen. #AntiLambda in rap range;#it{p}_{T} (GeV/#it{c})",200,0.,20.);
1609 fh1MCXiPt = new TH1F("fh1MCXiPt", "MC gen. #Xi^{-/o};#it{p}_{T} (GeV/#it{c})",200,0.,20.);
1610 fh1MCXibarPt = new TH1F("fh1MCXibarPt", "MC gen. #bar{#Xi}^{+/o};#it{p}_{T} (GeV/#it{c})",200,0.,20.);
1611 fh2MCEtaVsPtK0s = new TH2F("fh2MCEtaVsPtK0s","MC gen. K^{0}_{s} #eta; #it{p}_{T}",200,0.,20.,200,-1.,1.);
1612 fh2MCEtaVsPtLa = new TH2F("fh2MCEtaVsPtLa","MC gen. #Lambda #eta; #it{p}_{T}",200,0.,20.,200,-1.,1.);
1613 fh2MCEtaVsPtALa = new TH2F("fh2MCEtaVsPtALa","MC gen. #bar{#Lambda} #eta; #it{p}_{T}",200,0.,20.,200,-1.,1.);
1616 //fh1MCRapK0s = new TH1F("fh1MCRapK0s", "MC gen. K0s;rap with cut",200,-10,10);
1617 //fh1MCRapLambda = new TH1F("fh1MCRapLambda", "MC gen. #Lambda;rap",200,-10,10);
1618 //fh1MCRapAntiLambda = new TH1F("fh1MCRapAntiLambda", "MC gen. #bar{#Lambda};rap",200,-10,10);
1619 fh1MCEtaAllK0s = new TH1F("fh1MCEtaAllK0s", "MC gen. K0s;#eta",200,-1.,1.);
1620 fh1MCEtaK0s = new TH1F("fh1MCEtaK0s", "MC gen. K0s;#eta with cut",200,-1.,1.);
1621 fh1MCEtaLambda = new TH1F("fh1MCEtaLambda", "MC gen. #Lambda;#eta",200,-1.,1.);
1622 fh1MCEtaAntiLambda = new TH1F("fh1MCEtaAntiLambda", "MC gen. #bar{#Lambda};#eta",200,-1.,1.);
1624 fV0QAK0->DefineHistos();
1625 fFFHistosRecCuts->DefineHistos();
1626 fFFHistosRecCutsK0Evt->DefineHistos();
1627 /* fFFHistosIMK0AllEvt->DefineHistos();
1628 fFFHistosIMK0Jet->DefineHistos();
1629 fFFHistosIMK0Cone->DefineHistos();
1630 fFFHistosIMLaAllEvt->DefineHistos();
1631 fFFHistosIMLaJet->DefineHistos();
1632 fFFHistosIMLaCone->DefineHistos();
1633 fFFHistosIMALaAllEvt->DefineHistos();
1634 fFFHistosIMALaJet->DefineHistos();
1635 fFFHistosIMALaCone->DefineHistos();
1638 const Int_t saveLevel = 5;
1641 fCommonHistList->Add(fh1EvtAllCent);
1642 fCommonHistList->Add(fh1Evt);
1643 fCommonHistList->Add(fh1EvtSelection);
1644 fCommonHistList->Add(fh1EvtCent);
1645 fCommonHistList->Add(fh1VertexNContributors);
1646 fCommonHistList->Add(fh1VertexZ);
1647 fCommonHistList->Add(fh1Xsec);
1648 fCommonHistList->Add(fh1Trials);
1649 fCommonHistList->Add(fh1PtHard);
1650 fCommonHistList->Add(fh1PtHardTrials);
1651 fCommonHistList->Add(fh1nRecJetsCuts);
1652 fCommonHistList->Add(fh1EvtMult);
1653 fCommonHistList->Add(fh1K0Mult);
1654 fCommonHistList->Add(fh1dPhiJetK0);
1655 fCommonHistList->Add(fh1LaMult);
1656 fCommonHistList->Add(fh1dPhiJetLa);
1657 fCommonHistList->Add(fh1ALaMult);
1658 fCommonHistList->Add(fh1dPhiJetALa);
1659 fCommonHistList->Add(fh1JetEta);
1660 fCommonHistList->Add(fh1JetPhi);
1661 fCommonHistList->Add(fh2JetEtaPhi);
1662 //fCommonHistList->Add(fh1V0JetPt);
1663 fCommonHistList->Add(fh1IMK0Cone);
1664 fCommonHistList->Add(fh1IMLaCone);
1665 fCommonHistList->Add(fh1IMALaCone);
1666 fCommonHistList->Add(fh2FFJetTrackEta);
1667 // fCommonHistList->Add(fh1trackPosNCls);
1668 //fCommonHistList->Add(fh1trackNegNCls);
1669 fCommonHistList->Add(fh1trackPosEta);
1670 fCommonHistList->Add(fh1trackNegEta);
1671 fCommonHistList->Add(fh1V0Eta);
1672 // fCommonHistList->Add(fh1V0totMom);
1673 fCommonHistList->Add(fh1CosPointAngle);
1674 fCommonHistList->Add(fh1DecayLengthV0);
1675 fCommonHistList->Add(fh2ProperLifetimeK0sVsPtBeforeCut);
1676 fCommonHistList->Add(fh2ProperLifetimeK0sVsPtAfterCut);
1677 fCommonHistList->Add(fh1V0Radius);
1678 fCommonHistList->Add(fh1DcaV0Daughters);
1679 fCommonHistList->Add(fh1DcaPosToPrimVertex);
1680 fCommonHistList->Add(fh1DcaNegToPrimVertex);
1681 fCommonHistList->Add(fh2ArmenterosBeforeCuts);
1682 fCommonHistList->Add(fh2ArmenterosAfterCuts);
1683 fCommonHistList->Add(fh2BBLaPos);
1684 fCommonHistList->Add(fh2BBLaNeg);
1685 fCommonHistList->Add(fh1PosDaughterCharge);
1686 fCommonHistList->Add(fh1NegDaughterCharge);
1687 fCommonHistList->Add(fh1PtMCK0s);
1688 fCommonHistList->Add(fh1PtMCLa);
1689 fCommonHistList->Add(fh1PtMCALa);
1690 fCommonHistList->Add(fh1EtaK0s);
1691 fCommonHistList->Add(fh1EtaLa);
1692 fCommonHistList->Add(fh1EtaALa);
1693 fCommonHistList->Add(fhnInvMassEtaTrackPtK0s);
1694 fCommonHistList->Add(fhnInvMassEtaTrackPtLa);
1695 fCommonHistList->Add(fhnInvMassEtaTrackPtALa);
1696 fCommonHistList->Add(fh1TrackMultCone);
1697 fCommonHistList->Add(fh2TrackMultCone);
1698 fCommonHistList->Add(fhnNJK0);
1699 fCommonHistList->Add(fhnNJLa);
1700 fCommonHistList->Add(fhnNJALa);
1701 fCommonHistList->Add(fh2MCgenK0Cone);
1702 fCommonHistList->Add(fh2MCgenLaCone);
1703 fCommonHistList->Add(fh2MCgenALaCone);
1704 fCommonHistList->Add(fh2MCEtagenK0Cone);
1705 fCommonHistList->Add(fh2MCEtagenLaCone);
1706 fCommonHistList->Add(fh2MCEtagenALaCone);
1707 fCommonHistList->Add(fh1IMK0ConeSmear);
1708 fCommonHistList->Add(fh1IMLaConeSmear);
1709 fCommonHistList->Add(fh1IMALaConeSmear);
1710 fCommonHistList->Add(fhnrecMCHijingLaIncl);
1711 fCommonHistList->Add(fhnrecMCHijingLaCone);
1712 fCommonHistList->Add(fhnrecMCHijingALaIncl);
1713 fCommonHistList->Add(fhnrecMCHijingALaCone);
1714 fCommonHistList->Add(fhnrecMCInjectLaIncl);
1715 fCommonHistList->Add(fhnrecMCInjectLaCone);
1716 fCommonHistList->Add(fhnrecMCInjectALaIncl);
1717 fCommonHistList->Add(fhnrecMCInjectALaCone);
1718 fCommonHistList->Add(fhnrecMCHijingLaIncl);
1719 fCommonHistList->Add(fhnrecMCHijingLaCone);
1720 fCommonHistList->Add(fhnrecMCHijingALaIncl);
1721 fCommonHistList->Add(fhnrecMCHijingALaCone);
1722 fCommonHistList->Add(fhnrecMCInjectLaIncl);
1723 fCommonHistList->Add(fhnrecMCInjectLaCone);
1724 fCommonHistList->Add(fhnrecMCInjectALaIncl);
1725 fCommonHistList->Add(fhnrecMCInjectALaCone);
1726 fCommonHistList->Add(fhnMCrecK0Cone);
1727 fCommonHistList->Add(fhnMCrecLaCone);
1728 fCommonHistList->Add(fhnMCrecALaCone);
1729 fCommonHistList->Add(fhnMCrecK0ConeSmear);
1730 fCommonHistList->Add(fhnMCrecLaConeSmear);
1731 fCommonHistList->Add(fhnMCrecALaConeSmear);
1732 fCommonHistList->Add(fhnK0sSecContinCone);
1733 fCommonHistList->Add(fhnLaSecContinCone);
1734 fCommonHistList->Add(fhnALaSecContinCone);
1735 fCommonHistList->Add(fhnK0sIncl);
1736 fCommonHistList->Add(fhnK0sCone);
1737 fCommonHistList->Add(fhnLaIncl);
1738 fCommonHistList->Add(fhnLaCone);
1739 fCommonHistList->Add(fhnALaIncl);
1740 fCommonHistList->Add(fhnALaCone);
1741 fCommonHistList->Add(fhnK0sPC);
1742 fCommonHistList->Add(fhnLaPC);
1743 fCommonHistList->Add(fhnALaPC);
1744 fCommonHistList->Add(fhnK0sMCC);
1745 fCommonHistList->Add(fhnLaMCC);
1746 fCommonHistList->Add(fhnALaMCC);
1747 fCommonHistList->Add(fhnK0sRC);
1748 fCommonHistList->Add(fhnLaRC);
1749 fCommonHistList->Add(fhnALaRC);
1750 fCommonHistList->Add(fhnK0sOC);
1751 fCommonHistList->Add(fhnLaOC);
1752 fCommonHistList->Add(fhnALaOC);
1753 fCommonHistList->Add(fh1AreaExcluded);
1754 fCommonHistList->Add(fh1MedianEta);
1755 fCommonHistList->Add(fh1JetPtMedian);
1756 fCommonHistList->Add(fh1MCMultiplicityPrimary);
1757 fCommonHistList->Add(fh1MCMultiplicityTracks);
1758 fCommonHistList->Add(fhnFeedDownLa);
1759 fCommonHistList->Add(fhnFeedDownALa);
1760 fCommonHistList->Add(fhnFeedDownLaCone);
1761 fCommonHistList->Add(fhnFeedDownALaCone);
1762 fCommonHistList->Add(fh1MCProdRadiusK0s);
1763 fCommonHistList->Add(fh1MCProdRadiusLambda);
1764 fCommonHistList->Add(fh1MCProdRadiusAntiLambda);
1765 fCommonHistList->Add(fh1MCPtV0s);
1766 fCommonHistList->Add(fh1MCPtK0s);
1767 fCommonHistList->Add(fh1MCPtLambda);
1768 fCommonHistList->Add(fh1MCPtAntiLambda);
1769 fCommonHistList->Add(fh1MCXiPt);
1770 fCommonHistList->Add(fh1MCXibarPt);
1771 fCommonHistList->Add(fh2MCEtaVsPtK0s);
1772 fCommonHistList->Add(fh2MCEtaVsPtLa);
1773 fCommonHistList->Add(fh2MCEtaVsPtALa);
1774 //fCommonHistList->Add(fh1MCRapK0s);
1775 //fCommonHistList->Add(fh1MCRapLambda);
1776 //fCommonHistList->Add(fh1MCRapAntiLambda);
1777 fCommonHistList->Add(fh1MCEtaAllK0s);
1778 fCommonHistList->Add(fh1MCEtaK0s);
1779 fCommonHistList->Add(fh1MCEtaLambda);
1780 fCommonHistList->Add(fh1MCEtaAntiLambda);
1784 fV0QAK0->AddToOutput(fCommonHistList);
1785 fFFHistosRecCuts->AddToOutput(fCommonHistList);
1786 fFFHistosRecCutsK0Evt->AddToOutput(fCommonHistList);
1787 // fFFHistosIMK0AllEvt->AddToOutput(fCommonHistList);
1788 // fFFHistosIMK0Jet->AddToOutput(fCommonHistList);
1789 // fFFHistosIMK0Cone->AddToOutput(fCommonHistList);
1790 // fFFHistosIMLaAllEvt->AddToOutput(fCommonHistList);
1791 // fFFHistosIMLaJet->AddToOutput(fCommonHistList);
1792 // fFFHistosIMLaCone->AddToOutput(fCommonHistList);
1793 // fFFHistosIMALaAllEvt->AddToOutput(fCommonHistList);
1794 // fFFHistosIMALaJet->AddToOutput(fCommonHistList);
1795 // fFFHistosIMALaCone->AddToOutput(fCommonHistList);
1800 // =========== Switch on Sumw2 for all histos ===========
1801 for (Int_t i=0; i<fCommonHistList->GetEntries(); ++i){
1803 TH1 *h1 = dynamic_cast<TH1*>(fCommonHistList->At(i));
1805 if (h1) h1->Sumw2();//The error per bin will be computed as sqrt(sum of squares of weight) for each bin
1807 THnSparse *hnSparse = dynamic_cast<THnSparse*>(fCommonHistList->At(i));
1808 if(hnSparse) hnSparse->Sumw2();
1812 TH1::AddDirectory(oldStatus);
1813 PostData(1, fCommonHistList);
1816 //_______________________________________________
1817 void AliAnalysisTaskJetChem::UserExec(Option_t *)
1820 // Called for each event
1822 if(fDebug > 1) Printf("AliAnalysisTaskJetChem::UserExec()");
1824 if(fDebug > 1) Printf("Analysis event #%5d", (Int_t) fEntry);
1826 // Trigger selection
1827 AliInputEventHandler* inputHandler = (AliInputEventHandler*)
1828 ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1831 //for AliPIDResponse:
1832 //AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
1833 //AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1834 fPIDResponse = inputHandler->GetPIDResponse();
1836 if (!fPIDResponse){if(fDebug > 1) Printf("AliAnalysisTaskJetChem::UserExec(): fPIDResponse does not exist!"); return;}
1838 //std::cout<<"inputHandler->IsEventSelected(): "<<inputHandler->IsEventSelected()<<std::endl;
1839 //std::cout<<"fEvtSelectionMask: "<<fEvtSelectionMask<<std::endl;
1841 if(!(inputHandler->IsEventSelected() & fEvtSelectionMask)){
1842 //std::cout<<"########event rejected!!############"<<std::endl;
1843 fh1EvtSelection->Fill(1.);
1844 if (fDebug > 1 ) Printf(" Trigger Selection: event REJECTED ... ");
1845 PostData(1, fCommonHistList);
1849 fESD = dynamic_cast<AliESDEvent*>(InputEvent());//casting of pointers for inherited class, only for ESDs
1851 if(fDebug>3) Printf("%s:%d ESDEvent not found in the input", (char*)__FILE__,__LINE__);
1854 fMCEvent = MCEvent();
1856 if(fDebug>3) Printf("%s:%d MCEvent not found in the input", (char*)__FILE__,__LINE__);
1859 // get AOD event from input/output
1860 TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
1861 if( handler && handler->InheritsFrom("AliAODInputHandler") ) {
1862 fAOD = ((AliAODInputHandler*)handler)->GetEvent();
1863 if(fUseAODInputJets) fAODJets = fAOD;
1864 if (fDebug > 1) Printf("%s:%d AOD event from input", (char*)__FILE__,__LINE__);
1867 handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
1868 if( handler && handler->InheritsFrom("AliAODHandler") ) {
1869 fAOD = ((AliAODHandler*)handler)->GetAOD();
1871 if (fDebug > 1) Printf("%s:%d AOD event from output", (char*)__FILE__,__LINE__);
1875 if(!fAODJets && !fUseAODInputJets){ // case we have AOD in input & output and want jets from output
1876 TObject* outHandler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
1877 if( outHandler && outHandler->InheritsFrom("AliAODHandler") ){
1878 fAODJets = ((AliAODHandler*)outHandler)->GetAOD();
1879 if (fDebug > 1) Printf("%s:%d jets from output AOD", (char*)__FILE__,__LINE__);
1883 if(fNonStdFile.Length()!=0){
1884 // case we have an AOD extension - fetch the jets from the extended output
1886 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1887 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
1889 if(fDebug>1)Printf("AODExtension not found for %s",fNonStdFile.Data());
1894 Printf("%s:%d AODEvent not found", (char*)__FILE__,__LINE__);
1898 Printf("%s:%d AODEvent with jet branch not found", (char*)__FILE__,__LINE__);
1902 //primary vertex position:
1903 AliAODVertex *myPrimaryVertex = NULL;
1904 myPrimaryVertex = (AliAODVertex*)fAOD->GetPrimaryVertex();
1905 if (!myPrimaryVertex) return;
1906 fh1Evt->Fill(1.);//fill in every event that was accessed with InputHandler
1908 // event selection *****************************************
1910 // *** vertex cut ***
1911 AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
1912 Int_t nTracksPrim = primVtx->GetNContributors();
1913 fh1VertexNContributors->Fill(nTracksPrim);
1915 if (fDebug > 1) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim);
1917 if(nTracksPrim <= 2){
1918 if (fDebug > 1) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__);
1919 fh1EvtSelection->Fill(3.);
1920 PostData(1, fCommonHistList);
1924 fh1VertexZ->Fill(primVtx->GetZ());
1926 if(TMath::Abs(primVtx->GetZ())>fMaxVertexZ){
1927 if (fDebug > 1) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
1928 fh1EvtSelection->Fill(4.);
1929 PostData(1, fCommonHistList);
1933 // accepts only events that have same "primary" and SPD vertex, special issue of LHC11h PbPb data
1935 //fAOD: pointer to global primary vertex
1937 const AliAODVertex* spdVtx = fAOD->GetPrimaryVertexSPD();
1939 if (TMath::Abs(spdVtx->GetZ() - primVtx->GetZ())>fDeltaVertexZ) { if (fDebug > 1) Printf("deltaZVertex: event REJECTED..."); return;}
1942 //check for vertex radius to be smaller than 1 cm, (that was first applied by Vit Kucera in his analysis)
1944 Double_t vtxX = primVtx->GetX();
1945 Double_t vtxY = primVtx->GetY();
1947 if(TMath::Sqrt(vtxX*vtxX + vtxY*vtxY)>=1){
1948 if (fDebug > 1) Printf("%s:%d primary vertex r = %f: event REJECTED...",(char*)__FILE__,__LINE__,TMath::Sqrt(vtxX*vtxX + vtxY*vtxY));
1953 TString primVtxName(primVtx->GetName());
1955 if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){
1956 if (fDebug > 1) Printf("%s:%d primary vertex selection: TPC vertex, event REJECTED...",(char*)__FILE__,__LINE__);
1957 fh1EvtSelection->Fill(5.);
1958 PostData(1, fCommonHistList);
1962 Bool_t selectedHelper = AliAnalysisHelperJetTasks::Selected();
1963 if(!selectedHelper){
1964 fh1EvtSelection->Fill(6.);
1965 PostData(1, fCommonHistList);
1969 // event selection *****************************************
1971 Double_t centPercent = -1;
1975 if(handler && handler->InheritsFrom("AliAODInputHandler")){
1977 centPercent = fAOD->GetHeader()->GetCentrality();
1979 //std::cout<<"centPercent: "<<centPercent<<std::endl;
1981 fh1EvtAllCent->Fill(centPercent);
1983 if(centPercent>10) cl = 2; //standard PWG-JE binning
1984 if(centPercent>30) cl = 3;
1985 if(centPercent>50) cl = 4;
1989 if(centPercent < 0) cl = -1;
1990 if(centPercent >= 0) cl = 1;
1991 if(centPercent > 10) cl = 2; //standard PWG-JE binning
1992 if(centPercent > 30) cl = 3;
1993 if(centPercent > 50) cl = 4;
1994 if(centPercent > 80) cl = 5; //takes centralities higher than my upper edge of 80%, not to be used
1999 cl = AliAnalysisHelperJetTasks::EventClass();
2001 if(fESD) centPercent = fESD->GetCentrality()->GetCentralityPercentile("V0M"); //ESD JetServices Task has the centrality binning 0-10,10-30,30-50,50-80
2002 fh1EvtAllCent->Fill(centPercent);
2005 if(cl!=fEventClass){ // event not in selected event class, reject event#########################################
2007 if (fDebug > 1) Printf("%s:%d event not in selected event class: event REJECTED ...",(char*)__FILE__,__LINE__);
2008 fh1EvtSelection->Fill(2.);
2009 PostData(1, fCommonHistList);
2012 }//end if fEventClass > 0
2015 if (fDebug > 1) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__);
2018 //Printf("Analysis event #%5d", (Int_t) fEntry);
2020 fh1EvtSelection->Fill(0.);
2021 fh1EvtCent->Fill(centPercent);
2023 //___ get MC information __________________________________________________________________
2026 Double_t ptHard = 0.; //parton energy bins -> energy of particle
2027 Double_t nTrials = 1; // trials for MC trigger weight for real data
2030 AliGenEventHeader* genHeader = fMCEvent->GenEventHeader();
2031 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);//check usage of Pythia (pp) or Hijing (PbPb)
2032 AliGenHijingEventHeader* hijingGenHeader = 0x0;
2034 if(pythiaGenHeader){
2035 if(fDebug>3) Printf("%s:%d pythiaGenHeader found", (char*)__FILE__,__LINE__);
2036 nTrials = pythiaGenHeader->Trials();
2037 ptHard = pythiaGenHeader->GetPtHard();
2039 fh1PtHard->Fill(ptHard);
2040 fh1PtHardTrials->Fill(ptHard,nTrials);
2043 } else { // no pythia, hijing?
2045 if(fDebug>3) Printf("%s:%d no pythiaGenHeader found", (char*)__FILE__,__LINE__);
2047 hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
2048 if(!hijingGenHeader){
2049 if(fDebug>3) Printf("%s:%d no pythiaGenHeader or hjingGenHeader found", (char*)__FILE__,__LINE__);
2051 if(fDebug>3) Printf("%s:%d hijingGenHeader found", (char*)__FILE__,__LINE__);
2055 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
2058 //____ fetch jets _______________________________________________________________
2060 Int_t nJCuts = GetListOfJets(fJetsRecCuts, kJetsRecAcceptance);//fetch list with jets that survived all jet cuts: fJetsRecCuts
2062 Int_t nRecJetsCuts = 0; //number of reconstructed jets after jet cuts
2063 if(nJCuts>=0) nRecJetsCuts = fJetsRecCuts->GetEntries();
2064 if(fDebug>2)Printf("%s:%d Selected Rec jets after cuts: %d %d",(char*)__FILE__,__LINE__,nJCuts,nRecJetsCuts);
2065 if(nRecJetsCuts != nJCuts) Printf("%s:%d Mismatch selected Rec jets after cuts: %d %d",(char*)__FILE__,__LINE__,nJCuts,nRecJetsCuts);
2066 fh1nRecJetsCuts->Fill(nRecJetsCuts);
2069 //____ fetch background clusters ___________________________________________________
2070 if(fBranchRecBckgClusters.Length() != 0){
2072 Int_t nBJ = GetListOfBckgJets(fBckgJetsRec, kJetsRec);
2073 Int_t nRecBckgJets = 0;
2074 if(nBJ>=0) nRecBckgJets = fBckgJetsRec->GetEntries();
2075 if(fDebug>2)Printf("%s:%d Selected Rec background jets: %d %d",(char*)__FILE__,__LINE__,nBJ,nRecBckgJets);
2076 if(nBJ != nRecBckgJets) Printf("%s:%d Mismatch Selected Rec background jets: %d %d",(char*)__FILE__,__LINE__,nBJ,nRecBckgJets);
2080 //____ fetch reconstructed particles __________________________________________________________
2082 Int_t nTCuts = GetListOfTracks(fTracksRecCuts, kTrackAODCuts);//all tracks of event
2083 if(fDebug>2)Printf("%s:%d selected reconstructed tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,fTracksRecCuts->GetEntries());
2084 if(fTracksRecCuts->GetEntries() != nTCuts)
2085 Printf("%s:%d Mismatch selected reconstructed tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nTCuts,fTracksRecCuts->GetEntries());
2086 fh1EvtMult->Fill(fTracksRecCuts->GetEntries());
2088 Int_t nK0s = GetListOfV0s(fListK0s,fK0Type,kK0,myPrimaryVertex,fAOD);//all V0s in event with K0s assumption
2090 if(fDebug>5){std::cout<<"fK0Type: "<<fK0Type<<" kK0: "<<kK0<<" myPrimaryVertex: "<<myPrimaryVertex<<" fAOD: "<<fAOD<<std::endl;}
2092 //std::cout<< "nK0s: "<<nK0s<<std::endl;
2094 if(fDebug>2)Printf("%s:%d Selected V0s after cuts: %d %d",(char*)__FILE__,__LINE__,nK0s,fListK0s->GetEntries());
2095 if(nK0s != fListK0s->GetEntries()) Printf("%s:%d Mismatch selected K0s: %d %d",(char*)__FILE__,__LINE__,nK0s,fListK0s->GetEntries());
2096 fh1K0Mult->Fill(fListK0s->GetEntries());
2099 Int_t nLa = GetListOfV0s(fListLa,fLaType,kLambda,myPrimaryVertex,fAOD);//all V0s in event with Lambda particle assumption
2100 if(fDebug>2)Printf("%s:%d Selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nLa,fListLa->GetEntries());
2101 if(nLa != fListLa->GetEntries()) Printf("%s:%d Mismatch selected La: %d %d",(char*)__FILE__,__LINE__,nLa,fListLa->GetEntries());
2102 fh1LaMult->Fill(fListLa->GetEntries());
2104 Int_t nALa = GetListOfV0s(fListALa,fALaType,kAntiLambda,myPrimaryVertex,fAOD);//all V0s in event with Antilambda particle assumption
2105 if(fDebug>2)Printf("%s:%d Selected Rec tracks after cuts: %d %d",(char*)__FILE__,__LINE__,nALa,fListALa->GetEntries());
2106 if(nALa != fListALa->GetEntries()) Printf("%s:%d Mismatch selected ALa: %d %d",(char*)__FILE__,__LINE__,nALa,fListALa->GetEntries());
2107 fh1ALaMult->Fill(fListALa->GetEntries());
2111 //fetch MC gen particles_______________________________________________________
2113 if(fAnalysisMC){ // here
2115 //fill feeddown histo for associated particles
2117 // Access MC generated particles, fill TLists and histograms :
2119 Int_t nMCgenK0s = GetListOfMCParticles(fListMCgenK0s,kK0,fAOD); //fill TList with MC generated primary true K0s (list to fill, particletype, mc aod event)
2120 if(nMCgenK0s != fListMCgenK0s->GetEntries()) Printf("%s:%d Mismatch selected MCgenK0s: %d %d",(char*)__FILE__,__LINE__,nMCgenK0s,fListMCgenK0s->GetEntries());
2123 for(Int_t it=0; it<fListMCgenK0s->GetSize(); ++it){ // loop MC generated K0s, filling histograms
2125 AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenK0s->At(it));
2130 //Double_t fRapCurrentPart = MyRapidity(mcp0->E(),mcp0->Pz());
2131 Double_t fEtaCurrentPart = mcp0->Eta();
2132 Double_t fPtCurrentPart = mcp0->Pt();
2134 fh1MCEtaK0s->Fill(fEtaCurrentPart);
2135 //fh1MCRapK0s->Fill(fRapCurrentPart);
2136 fh1MCPtK0s->Fill(fPtCurrentPart);
2138 fh2MCEtaVsPtK0s->Fill(fPtCurrentPart,fEtaCurrentPart); //eta cut, physical primary selection and decay mode considered
2143 Int_t nMCgenLa = GetListOfMCParticles(fListMCgenLa,kLambda,fAOD); //fill TList with MC generated primary true Lambdas (list to fill, particletype, mc aod event)
2144 if(nMCgenLa != fListMCgenLa->GetEntries()) Printf("%s:%d Mismatch selected MCgenLa: %d %d",(char*)__FILE__,__LINE__,nMCgenLa,fListMCgenLa->GetEntries());
2147 for(Int_t it=0; it<fListMCgenLa->GetSize(); ++it){ // loop MC generated La, filling histograms
2149 AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenLa->At(it));
2154 //Double_t fRapCurrentPart = MyRapidity(mcp0->E(),mcp0->Pz());
2155 Double_t fEtaCurrentPart = mcp0->Eta();
2156 Double_t fPtCurrentPart = mcp0->Pt();
2158 fh1MCEtaLambda->Fill(fEtaCurrentPart);
2159 //fh1MCRapLambda->Fill(fRapCurrentPart);
2160 fh1MCPtLambda->Fill(fPtCurrentPart);
2161 fh2MCEtaVsPtLa->Fill(fPtCurrentPart,fEtaCurrentPart); //eta cut, physical primary selection and decay mode considered
2166 Int_t nMCgenALa = GetListOfMCParticles(fListMCgenALa,kAntiLambda,fAOD); //fill TList with MC generated primary true Antilambdas (list to fill, particletype, mc aod event)
2167 if(nMCgenALa != fListMCgenALa->GetEntries()) Printf("%s:%d Mismatch selected MCgenALa: %d %d",(char*)__FILE__,__LINE__,nMCgenALa,fListMCgenALa->GetEntries());
2170 for(Int_t it=0; it<fListMCgenALa->GetSize(); ++it){ // loop MC generated ALa, filling histograms
2172 AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenALa->At(it));
2175 //MC gen Antilambdas
2177 //Double_t fRapCurrentPart = MyRapidity(mcp0->E(),mcp0->Pz());
2178 Double_t fEtaCurrentPart = mcp0->Eta();
2179 Double_t fPtCurrentPart = mcp0->Pt();
2181 fh1MCEtaAntiLambda->Fill(fEtaCurrentPart);
2182 //fh1MCRapAntiLambda->Fill(fRapCurrentPart);
2183 fh1MCPtAntiLambda->Fill(fPtCurrentPart);
2184 fh2MCEtaVsPtALa->Fill(fPtCurrentPart,fEtaCurrentPart); //eta cut, physical primary selection and decay mode considered
2190 //loop over MC feeddown candidates in TList
2195 } //end MCAnalysis part for gen particles
2198 // ___ V0 QA + K0s + La + ALa pt spectra all events _______________________________________________
2200 Double_t lPrimaryVtxPosition[3];
2201 Double_t lV0Position[3];
2202 lPrimaryVtxPosition[0] = primVtx->GetX();
2203 lPrimaryVtxPosition[1] = primVtx->GetY();
2204 lPrimaryVtxPosition[2] = primVtx->GetZ();
2205 Double_t dRadiusExcludeCone = 2*GetFFRadius(); //2 times jet radius
2206 //------------------------------------------
2207 for(Int_t it=0; it<fListK0s->GetSize(); ++it){
2209 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListK0s->At(it));
2214 // VO's main characteristics to check the reconstruction cuts
2216 //Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2219 Double_t fV0Radius = -999;
2220 Double_t fDcaV0Daughters = v0->DcaV0Daughters();
2221 Double_t fDcaPosToPrimVertex = v0->DcaPosToPrimVertex();//IP of positive charged daughter
2222 Double_t fDcaNegToPrimVertex = v0->DcaNegToPrimVertex();//IP of negative charged daughter
2223 Int_t negDaughterpdg = 0;
2224 Int_t posDaughterpdg = 0;
2225 Int_t motherType = 0;
2228 Bool_t fPhysicalPrimary = kFALSE;//don't use IsPhysicalPrimary() anymore for MC analysis, use instead 2D distance from primary to secondary vertex
2229 Int_t MCv0PdgCode = 0;
2231 AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0));
2232 AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1));
2234 Double_t PosEta = trackPos->AliAODTrack::Eta();//daughter track charge is sometimes wrong here, account for that!!!
2235 Double_t NegEta = trackNeg->AliAODTrack::Eta();
2237 //Double_t trackPosNcls = trackPos->GetTPCNcls();//Get number of clusters for positive charged tracks
2238 //Double_t trackNegNcls = trackNeg->GetTPCNcls();//Get number of clusters for negative charged tracks
2240 CalculateInvMass(v0, kK0, invMK0s, trackPt); //function to calculate invMass with TLorentzVector class
2242 //Double_t jetPt = fFFIMJetPtMin; // assign pro forma jet energy
2243 //Double_t fRap = v0->RapK0Short();
2244 Double_t fEta = v0->PseudoRapV0();
2245 Bool_t bIsInCone = kFALSE;//init boolean, is not in any cone (OC)
2247 for(Int_t ij=0; ij<nRecJetsCuts; ++ij){ // loop over all jets in event
2249 AliAODJet* jet = (AliAODJet*) (fJetsRecCuts->At(ij));
2251 if(IsParticleInCone(jet, v0, dRadiusExcludeCone) == kTRUE) {bIsInCone = kTRUE;}
2255 if(bIsInCone==kFALSE){//K0s is not part of any selected jet in event
2256 Double_t vK0sOC[3] = {invMK0s,trackPt,fEta};
2257 fhnK0sOC->Fill(vK0sOC);
2260 Double_t fV0cosPointAngle = v0->CosPointingAngle(lPrimaryVtxPosition);
2262 lV0Position[0]= v0->DecayVertexV0X();
2263 lV0Position[1]= v0->DecayVertexV0Y();
2264 lV0Position[2]= v0->DecayVertexV0Z();
2266 Double_t fV0DecayLength = v0->DecayLengthV0(lPrimaryVtxPosition);
2267 fV0Radius = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
2270 //fetch V0s outside of jet cones (outside of 2R):
2280 fV0QAK0->FillTrackQA(v0->Eta(), TVector2::Phi_0_2pi(v0->Phi()), v0->Pt());
2281 //fFFHistosIMK0AllEvt->FillFF(trackPt, invMK0s, jetPt, incrementJetPt);
2282 //fh1trackPosNCls->Fill(trackPosNcls);
2283 //fh1trackNegNCls->Fill(trackNegNcls);
2284 fh1EtaK0s->Fill(fEta);
2286 Double_t vK0sIncl[3] = {invMK0s,trackPt,fEta}; //fill all K0s in event into THnSparse of 3 dimensions
2287 fhnK0sIncl->Fill(vK0sIncl);
2291 TString generatorName;
2293 TList *listmc = fAOD->GetList();
2294 Bool_t mclabelcheck = MCLabelCheck(v0, kK0, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode, generatorName, isinjected);
2295 //if(fPhysicalPrimary == kFALSE)continue;
2296 //std::cout<<"mclabelcheck: "<<mclabelcheck<<std::endl;
2297 //std::cout<<"IsPhysicalPrimary: "<<fPhysicalPrimary<<std::endl;
2299 if(mclabelcheck == kFALSE)continue;
2301 Double_t vInvMassEtaTrackPtK0s[3] = {fEta,invMK0s,trackPt};
2302 fhnInvMassEtaTrackPtK0s->Fill(vInvMassEtaTrackPtK0s);//includes also feeddown particles, mainly phi particles whose decay products are considered here as primary
2305 fh1PtMCK0s->Fill(MCPt);
2309 fh1V0Eta->Fill(fEta);
2310 //fh1V0totMom->Fill(fV0TotalMomentum);
2311 fh1CosPointAngle->Fill(fV0cosPointAngle);
2312 fh1DecayLengthV0->Fill(fV0DecayLength);
2313 fh1V0Radius->Fill(fV0Radius);
2314 fh1DcaV0Daughters->Fill(fDcaV0Daughters);
2315 fh1DcaPosToPrimVertex->Fill(fDcaPosToPrimVertex);
2316 fh1DcaNegToPrimVertex->Fill(fDcaNegToPrimVertex);
2317 fh1trackPosEta->Fill(PosEta);
2318 fh1trackNegEta->Fill(NegEta);
2322 // __La pt spectra all events _______________________________________________
2325 for(Int_t it=0; it<fListLa->GetSize(); ++it){ // loop all La
2327 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListLa->At(it));
2330 // VO's main characteristics to check the reconstruction cuts
2331 // Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2334 Double_t fV0Radius = -999;
2335 Double_t fDcaV0Daughters = v0->DcaV0Daughters();
2336 Double_t fDcaPosToPrimVertex = v0->DcaPosToPrimVertex();//IP of positive charged daughter
2337 Double_t fDcaNegToPrimVertex = v0->DcaNegToPrimVertex();//IP of negative charged daughter
2338 Int_t negDaughterpdg = 0;
2339 Int_t posDaughterpdg = 0;
2340 Int_t motherType = 0;
2343 Bool_t fPhysicalPrimary = kFALSE;
2344 Int_t MCv0PdgCode = 0;
2345 AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0));
2346 AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1));
2348 //Double_t trackPosNcls = trackPos->GetTPCNcls();//Get number of clusters for positive charged tracks
2349 //Double_t trackNegNcls = trackNeg->GetTPCNcls();//Get number of clusters for negative charged tracks
2351 Double_t PosEta = trackPos->AliAODTrack::Eta();//daughter track charge is sometimes wrong here, account for that!!!
2352 Double_t NegEta = trackNeg->AliAODTrack::Eta();
2354 Double_t fEta = v0->PseudoRapV0();
2355 Bool_t bIsInCone = kFALSE;//init boolean, is not in any cone (OC)
2357 CalculateInvMass(v0, kLambda, invMLa, trackPt);//function to calculate invMass with TLorentzVector class
2360 for(Int_t ij=0; ij<nRecJetsCuts; ++ij){ // loop over all jets in event
2362 AliAODJet* jet = (AliAODJet*) (fJetsRecCuts->At(ij));
2364 if(IsParticleInCone(jet, v0, dRadiusExcludeCone) == kTRUE) {bIsInCone = kTRUE;
2368 if(bIsInCone == kFALSE){//success!
2369 Double_t vLaOC[3] = {invMLa, trackPt,fEta};
2370 fhnLaOC->Fill(vLaOC);
2373 // Double_t jetPt = fFFIMJetPtMin; // assign pro forma jet energy
2374 // Double_t fRap = v0->Y(3122);
2376 Double_t fV0DecayLength = v0->DecayLengthV0(lPrimaryVtxPosition);
2377 Double_t fV0cosPointAngle = v0->CosPointingAngle(lPrimaryVtxPosition);
2378 lV0Position[0]= v0->DecayVertexV0X();
2379 lV0Position[1]= v0->DecayVertexV0Y();
2380 lV0Position[2]= v0->DecayVertexV0Z();
2382 fV0Radius = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
2384 //fFFHistosIMLaAllEvt->FillFF(trackPt, invMLa, jetPt, incrementJetPt);
2385 //fh1trackPosNCls->Fill(trackPosNcls);
2386 //fh1trackNegNCls->Fill(trackNegNcls);
2387 fh1EtaLa->Fill(fEta);
2389 Double_t vLaIncl[3] = {invMLa,trackPt,fEta};
2390 fhnLaIncl->Fill(vLaIncl);
2394 TString generatorName;
2396 TList* listmc = fAOD->GetList();
2397 Bool_t mclabelcheck = MCLabelCheck(v0, kLambda, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode, generatorName, isinjected);
2398 if(mclabelcheck == kFALSE)continue;
2399 //if(fPhysicalPrimary == kFALSE)continue;
2401 if(generatorName == "Hijing"){
2402 Double_t vrecMCHijingLaIncl[3] = {invMLa,trackPt,fEta};
2403 fhnrecMCHijingLaIncl->Fill(vrecMCHijingLaIncl);
2406 if(isinjected == kTRUE){
2407 Double_t vrecMCInjectLaIncl[3] = {invMLa,trackPt,fEta};
2408 fhnrecMCInjectLaIncl->Fill(vrecMCInjectLaIncl);
2411 Double_t vInvMassEtaTrackPtLa[3] = {fEta,invMLa,trackPt};
2412 fhnInvMassEtaTrackPtLa->Fill(vInvMassEtaTrackPtLa);//includes also feed-down particles
2413 fh1PtMCLa->Fill(MCPt);
2416 fh1PtMCLa->Fill(MCPt);
2418 fh1V0Eta->Fill(fEta);
2419 //fh1V0totMom->Fill(fV0TotalMomentum);
2420 fh1CosPointAngle->Fill(fV0cosPointAngle);
2421 fh1DecayLengthV0->Fill(fV0DecayLength);
2422 fh1V0Radius->Fill(fV0Radius);
2423 fh1DcaV0Daughters->Fill(fDcaV0Daughters);
2424 fh1DcaPosToPrimVertex->Fill(fDcaPosToPrimVertex);
2425 fh1DcaNegToPrimVertex->Fill(fDcaNegToPrimVertex);
2426 fh1trackPosEta->Fill(PosEta);
2427 fh1trackNegEta->Fill(NegEta);
2430 // __ALa pt spectra all events _______________________________________________
2432 for(Int_t it=0; it<fListALa->GetSize(); ++it){ // loop all ALa
2434 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListALa->At(it));
2438 //VO's main characteristics to check the reconstruction cuts
2439 Double_t invMALa =0;
2441 Double_t fV0Radius = -999;
2442 Double_t fDcaV0Daughters = v0->DcaV0Daughters();
2443 Double_t fDcaPosToPrimVertex = v0->DcaPosToPrimVertex();//IP of positive charged daughter
2444 Double_t fDcaNegToPrimVertex = v0->DcaNegToPrimVertex();//IP of negative charged daughter
2445 Int_t negDaughterpdg = 0;
2446 Int_t posDaughterpdg = 0;
2447 Int_t motherType = 0;
2450 Bool_t fPhysicalPrimary = kFALSE;
2451 Int_t MCv0PdgCode = 0;
2453 AliAODTrack *trackPos = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(0));
2454 AliAODTrack *trackNeg = (AliAODTrack *) (v0->GetSecondaryVtx()->GetDaughter(1));
2456 Double_t PosEta = trackPos->AliAODTrack::Eta();//daughter track charge is sometimes wrong here, account for that!!!
2457 Double_t NegEta = trackNeg->AliAODTrack::Eta();
2459 Double_t fEta = v0->PseudoRapV0();
2460 Bool_t bIsInCone = kFALSE;//init boolean for OC
2463 CalculateInvMass(v0, kAntiLambda, invMALa, trackPt); //function to calculate invMass with TLorentzVector class
2465 for(Int_t ij=0; ij<nRecJetsCuts; ++ij){ // loop over all jets in event
2467 AliAODJet* jet = (AliAODJet*) (fJetsRecCuts->At(ij));
2469 if(IsParticleInCone(jet, v0, dRadiusExcludeCone) == kTRUE){
2474 if(bIsInCone == kFALSE){//success!
2475 Double_t vALaOC[3] = {invMALa, trackPt,fEta};
2476 fhnALaOC->Fill(vALaOC);
2479 //Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2480 //Double_t jetPt = fFFIMJetPtMin; // assign pro forma jet energy
2481 // Double_t fRap = v0->Y(-3122);
2484 Double_t fV0cosPointAngle = v0->CosPointingAngle(lPrimaryVtxPosition);
2485 lV0Position[0]= v0->DecayVertexV0X();
2486 lV0Position[1]= v0->DecayVertexV0Y();
2487 lV0Position[2]= v0->DecayVertexV0Z();
2488 Double_t fV0DecayLength = v0->DecayLengthV0(lPrimaryVtxPosition);
2489 fV0Radius = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
2491 //fFFHistosIMALaAllEvt->FillFF(trackPt, invMALa, jetPt, incrementJetPt);
2492 //fh1trackPosNCls->Fill(trackPosNcls);
2493 //fh1trackNegNCls->Fill(trackNegNcls);
2494 fh1EtaALa->Fill(fEta);
2496 Double_t vALaIncl[3] = {invMALa,trackPt,fEta};
2497 fhnALaIncl->Fill(vALaIncl);
2500 TString generatorName;
2502 TList* listmc = fAOD->GetList();
2503 Bool_t mclabelcheck = MCLabelCheck(v0, kAntiLambda, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode, generatorName, isinjected);
2504 if(mclabelcheck == kFALSE)continue;
2505 //if(fPhysicalPrimary == kFALSE)continue;
2507 if(generatorName == "Hijing"){
2508 Double_t vrecMCHijingALaIncl[3] = {invMALa,trackPt,fEta};
2509 fhnrecMCHijingALaIncl->Fill(vrecMCHijingALaIncl);
2512 if(isinjected == kTRUE){
2513 Double_t vrecMCInjectALaIncl[3] = {invMALa,trackPt,fEta};
2514 fhnrecMCInjectALaIncl->Fill(vrecMCInjectALaIncl);
2518 Double_t vInvMassEtaTrackPtALa[3] = {fEta,invMALa,trackPt};
2519 fhnInvMassEtaTrackPtALa->Fill(vInvMassEtaTrackPtALa);
2520 fh1PtMCALa->Fill(MCPt);
2523 fh1V0Eta->Fill(fEta);
2524 //fh1V0totMom->Fill(fV0TotalMomentum);
2525 fh1CosPointAngle->Fill(fV0cosPointAngle);
2526 fh1DecayLengthV0->Fill(fV0DecayLength);
2527 fh1V0Radius->Fill(fV0Radius);
2528 fh1DcaV0Daughters->Fill(fDcaV0Daughters);
2529 fh1DcaPosToPrimVertex->Fill(fDcaPosToPrimVertex);
2530 fh1DcaNegToPrimVertex->Fill(fDcaNegToPrimVertex);
2531 fh1trackPosEta->Fill(PosEta);
2532 fh1trackNegEta->Fill(NegEta);
2535 //_____no jets events______________________________________________________________________________________________________________________________________
2537 if(nRecJetsCuts == 0){//no jet events
2539 if(fDebug>6) { std::cout<<"################## nRecJetsCuts == 0 ###################"<<std::endl;
2540 std::cout<<"fListK0s->GetSize() in NJ event: "<<fListK0s->GetSize()<<std::endl;
2543 for(Int_t it=0; it<fListK0s->GetSize(); ++it){ // loop all K0s
2545 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListK0s->At(it));
2548 Double_t invMK0s =0;
2550 CalculateInvMass(v0, kK0, invMK0s, trackPt);
2551 Double_t fEta = v0->Eta();
2553 Double_t vNJK0[3] = {invMK0s,trackPt,fEta}; //fill all K0s in events wo selected jets
2554 fhnNJK0->Fill(vNJK0);
2558 for(Int_t it=0; it<fListLa->GetSize(); ++it){ // loop all La
2560 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListLa->At(it));
2565 CalculateInvMass(v0, kLambda, invMLa, trackPt);
2566 Double_t fEta = v0->Eta();
2568 Double_t vNJLa[3] = {invMLa,trackPt,fEta}; //fill all K0s in events wo selected jets
2569 fhnNJLa->Fill(vNJLa);
2573 for(Int_t it=0; it<fListALa->GetSize(); ++it){ // loop all ALa
2575 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListALa->At(it));
2578 Double_t invMALa =0;
2580 CalculateInvMass(v0, kAntiLambda, invMALa, trackPt);
2582 Double_t fEta = v0->Eta();
2584 Double_t vNJALa[3] = {invMALa,trackPt,fEta}; //fill all K0s in events wo selected jets
2585 fhnNJALa->Fill(vNJALa);
2592 //____ fill all jet related histos ________________________________________________________________________________________________________________________
2593 //##########################jet loop########################################################################################################################
2595 //fill jet histos in general
2596 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
2598 AliAODJet* jet = (AliAODJet*) (fJetsRecCuts->At(ij));
2600 Double_t jetPt = jet->Pt();
2601 Double_t jetEta = jet->Eta();
2602 Double_t jetPhi = jet->Phi();
2604 //if(ij==0){ // loop over leading jets for ij = 0, for ij>= 0 look into all jets
2606 if(ij>=0){//all jets in event
2608 TList* jettracklist = new TList();
2609 Double_t sumPt = 0.;
2610 Bool_t isBadJet = kFALSE;
2611 Int_t njetTracks = 0;
2613 if(GetFFRadius()<=0){
2614 GetJetTracksTrackrefs(jettracklist, jet, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJet);// list of jet tracks from trackrefs
2616 GetJetTracksPointing(fTracksRecCuts, jettracklist, jet, GetFFRadius(), sumPt, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJet); // fill list of tracks in cone around jet axis with cone Radius (= 0.4 standard)
2619 if(GetFFMinNTracks()>0 && jettracklist->GetSize() <= GetFFMinNTracks()) isBadJet = kTRUE; // reject jets with less tracks than fFFMinNTracks
2620 if(isBadJet) continue; // rejects jets in which no track has a track pt higher than 5 GeV/c (see AddTask macro)
2622 //Float_t fJetAreaMin = 0.6*TMath::Pi()*GetFFRadius()*GetFFRadius(); // minimum jet area cut
2624 //std::cout<<"GetFFRadius(): "<<GetFFRadius()<<std::endl;
2625 //std::cout<<"jet->EffectiveAreaCharged()"<<jet->EffectiveAreaCharged()<<std::endl;
2626 //std::cout<<"fJetAreaMin: "<<fJetAreaMin<<std::endl;
2628 //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
2630 Double_t dAreaExcluded = TMath::Pi()*dRadiusExcludeCone*dRadiusExcludeCone; // area of the cone
2631 dAreaExcluded -= AreaCircSegment(dRadiusExcludeCone,fCutjetEta-jet->Eta()); // positive eta overhang
2632 dAreaExcluded -= AreaCircSegment(dRadiusExcludeCone,fCutjetEta+jet->Eta()); // negative eta overhang
2633 fh1AreaExcluded->Fill(dAreaExcluded);//histo contains all areas that are jet related and have to be excluded concerning OC UE pt spectrum normalisation by area
2635 fh1JetEta->Fill(jetEta);
2636 fh1JetPhi->Fill(jetPhi);
2637 fh2JetEtaPhi->Fill(jetEta,jetPhi);
2639 // printf("pT = %f, eta = %f, phi = %f, leadtr pt = %f\n, ",jetPt,jetEta,jetphi,leadtrack);
2641 for(Int_t it=0; it<jettracklist->GetSize(); ++it){//loop over all particles in jet
2643 AliVParticle* trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));//all tracks in jet cone
2644 if(!trackVP)continue;
2646 Float_t trackPt = trackVP->Pt();//transversal momentum of jet particle
2647 Float_t trackEta = trackVP->Eta();
2649 Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2651 fFFHistosRecCuts->FillFF(trackPt, jetPt, incrementJetPt);//histo with tracks/jets after cut selection, for all events
2652 if(nK0s>0) fFFHistosRecCutsK0Evt->FillFF(trackPt, jetPt, incrementJetPt);//only for K0s events
2653 fh2FFJetTrackEta->Fill(trackEta,jetPt);
2658 njetTracks = jettracklist->GetSize();
2660 //____________________________________________________________________________________________________________________
2661 //strangeness constribution to jet cone
2665 TList *list = fAOD->GetList();
2666 AliAODMCHeader *mcHeadr=(AliAODMCHeader*)list->FindObject(AliAODMCHeader::StdBranchName());
2667 if(!mcHeadr)continue;
2669 Double_t mcXv=0., mcYv=0., mcZv=0.;//MC primary vertex position
2671 mcXv=mcHeadr->GetVtxX(); mcYv=mcHeadr->GetVtxY(); mcZv=mcHeadr->GetVtxZ(); // position of the MC primary vertex
2673 for(Int_t it=0; it<jettracklist->GetSize(); ++it){//loop over all tracks in the jet
2675 AliVParticle* trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));//track in jet cone
2676 if(!trackVP)continue;
2677 AliAODTrack *tr = dynamic_cast<AliAODTrack*> (trackVP); //fetch one jet track from the TList
2680 //get MC label information
2681 TList *mclist = fAOD->GetList();
2683 //fetch the MC stack
2684 TClonesArray* stackMC = (TClonesArray*)mclist->FindObject(AliAODMCParticle::StdBranchName()); //get MCAOD branch in data
2685 if (!stackMC) {Printf("ERROR: stack not available");}
2689 Int_t trackLabel = TMath::Abs(tr->GetLabel()); //fetch jet track label in MC stack
2691 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(stackMC->At(trackLabel)); //fetch MC gen. particle for rec. jet track
2693 if(!part)continue; //skip non-existing objects
2696 //Bool_t IsPhysicalPrimary = part->IsPhysicalPrimary();//not recommended to check, better use distance between primary vertex and secondary vertex
2698 Float_t fDistPrimaryMax = 0.01;
2699 // Get the distance between production point of the MC mother particle and the primary vertex
2701 Double_t dx = mcXv-part->Xv();//mc primary vertex - mc gen. v0 vertex
2702 Double_t dy = mcYv-part->Yv();
2703 Double_t dz = mcZv-part->Zv();
2705 Float_t fDistPrimary = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
2706 Bool_t fPhysicalPrimary = (fDistPrimary < fDistPrimaryMax);
2708 // std::cout<<"fDistPrimary"<<fDistPrimary<<std::endl;
2709 // std::cout<<"fPhysicalPrimary"<<fPhysicalPrimary<<std::endl;
2711 if(!fPhysicalPrimary)continue;//rejects Kstar and other strong decaying particles from Secondary Contamination
2713 Bool_t isFromStrange = kFALSE;// flag to check whether particle has strange mother
2715 Int_t iMother = part->GetMother(); //get mother MC gen. particle label
2718 AliAODMCParticle *partM = dynamic_cast<AliAODMCParticle*>(stackMC->At(iMother)); //fetch mother of MC gen. particle
2719 if(!partM) continue;
2721 Int_t codeM = TMath::Abs(partM->GetPdgCode()); //mothers pdg code
2723 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..)
2725 if (mfl == 3 && codeM != 3) isFromStrange = kTRUE;
2728 if(isFromStrange == kTRUE){
2730 Double_t trackPt = part->Pt();
2731 Double_t trackEta = part->Eta();
2732 //fh3StrContinCone->Fill(jetPt, trackPt, trackEta);//MC gen. particle parameters, but rec. jet pt
2734 }//isFromStrange is kTRUE
2736 }//end loop over jet tracks
2741 fh1TrackMultCone->Fill(njetTracks);
2742 fh2TrackMultCone->Fill(njetTracks,jetPt);
2746 // fQAJetHistosRecCutsLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt() );
2748 for(Int_t it=0; it<fListK0s->GetSize(); ++it){ // loop all K0s
2750 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListK0s->At(it));
2751 if(!v0) continue;//rejection of events with no V0 vertex
2755 TVector3 v0MomVect(v0Mom);
2757 Double_t dPhiJetK0 = (jet->MomentumVector()->Vect()).DeltaPhi(v0MomVect);
2758 // Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2760 // if(incrementJetPt){fh1V0JetPt->Fill(jetPt);}
2762 Double_t invMK0s =0;
2764 CalculateInvMass(v0, kK0, invMK0s, trackPt); //function to calculate invMass with TLorentzVector class
2766 // fFFHistosIMK0Jet->FillFF(trackPt, invMK0s, jetPt, incrementJetPt);
2769 if(dPhiJetK0<fh1dPhiJetK0->GetXaxis()->GetXmin()) dPhiJetK0 += 2*TMath::Pi();
2770 fh1dPhiJetK0->Fill(dPhiJetK0);
2774 // if(fListK0s->GetSize() == 0){ // no K0: increment jet pt spectrum
2776 // Bool_t incrementJetPt = kTRUE;
2777 // fFFHistosIMK0Jet->FillFF(-1, -1, jetPt, incrementJetPt);
2780 //____fetch reconstructed K0s in cone around jet axis:_______________________________________________________________________________
2782 TList* jetConeK0list = new TList();
2784 Double_t sumPtK0 = 0.;
2786 Bool_t isBadJetK0 = kFALSE; // dummy, do not use
2789 GetTracksInCone(fListK0s, jetConeK0list, jet, GetFFRadius(), sumPtK0, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetK0); //reconstructed K0s in cone around jet axis
2791 if(fDebug>2)Printf("%s:%d nK0s total: %d, in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nK0s,jetConeK0list->GetEntries(),GetFFRadius());
2794 for(Int_t it=0; it<jetConeK0list->GetSize(); ++it){ // loop for K0s in jet cone
2796 AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetConeK0list->At(it));
2799 Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
2800 Double_t invMK0s =0;
2805 CalculateInvMass(v0, kK0, invMK0s, trackPt); //function to calculate invMass with TLorentzVector class
2809 Double_t jetPtSmear = -1;
2810 SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
2811 if(incrementJetPt == kTRUE){fh1IMK0ConeSmear->Fill(jetPtSmear);} //fill TH1F for normalization purposes
2814 if(incrementJetPt==kTRUE){
2815 fh1IMK0Cone->Fill(jetPt);}//normalisation by number of selected jets
2817 //fFFHistosIMK0Cone->FillFF(trackPt, invMK0s, jetPt, incrementJetPt);
2819 Double_t vK0sCone[4] = {jetPt, invMK0s,trackPt,fEta};
2820 fhnK0sCone->Fill(vK0sCone);
2824 if(jetConeK0list->GetSize() == 0){ // no K0: increment jet pt spectrum
2827 Bool_t incrementJetPt = kTRUE;//jets without K0s will be only filled in TH1F only once, so no increment needed
2828 //fFFHistosIMK0Cone->FillFF(-1, -1, jetPt, incrementJetPt);
2829 Double_t vK0sCone[4] = {jetPt, -1, -1, -1};
2830 fhnK0sCone->Fill(vK0sCone);
2832 if(incrementJetPt==kTRUE){
2833 fh1IMK0Cone->Fill(jetPt);}//normalisation by number of selected jets
2836 Double_t jetPtSmear = -1;
2837 SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
2838 if(incrementJetPt == kTRUE){fh1IMK0ConeSmear->Fill(jetPtSmear);} //fill TH1F for normalization purposes
2842 //Random cones________________________________________________________________________
2844 if(ij==0){//fetch random cone V0s only once per event
2846 //______fetch random cones___________________________________________________________
2849 AliAODJet* jetRC = 0;
2850 jetRC = GetRandomCone(fJetsRecCuts, fCutjetEta, 2*GetFFRadius());//fetch one random cone for each event
2851 TList* fListK0sRC = new TList();//list for K0s in random cone (RC), one RC per event
2852 TList* fListLaRC = new TList();
2853 TList* fListALaRC = new TList();
2855 Double_t sumPtK0sRC = 0;
2856 Double_t sumPtLaRC = 0;
2857 Double_t sumPtALaRC = 0;
2858 Bool_t isBadJetK0sRC = kFALSE;
2859 Bool_t isBadJetLaRC = kFALSE;
2860 Bool_t isBadJetALaRC = kFALSE;
2865 GetTracksInCone(fListK0s, fListK0sRC, jetRC, GetFFRadius(), sumPtK0sRC, 0, 0, isBadJetK0sRC);
2867 for(Int_t it=0; it<fListK0sRC->GetSize(); ++it){ // loop for K0s in random cone
2869 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListK0sRC->At(it));
2872 Double_t invMK0s =0;
2877 CalculateInvMass(v0, kK0, invMK0s, trackPt); //function to calculate invMass with TLorentzVector class
2879 Double_t vK0sRC[3] = {invMK0s,trackPt,fEta};
2880 fhnK0sRC->Fill(vK0sRC);
2885 GetTracksInCone(fListLa, fListLaRC, jetRC, GetFFRadius(), sumPtLaRC, 0, 0, isBadJetLaRC);
2887 for(Int_t it=0; it<fListLaRC->GetSize(); ++it){ // loop for Lambdas in random cone
2889 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListLaRC->At(it));
2897 CalculateInvMass(v0, kLambda, invMLa, trackPt); //function to calculate invMass with TLorentzVector class
2899 Double_t vLaRC[3] = {invMLa,trackPt,fEta};
2900 fhnLaRC->Fill(vLaRC);
2905 GetTracksInCone(fListALa, fListALaRC, jetRC, GetFFRadius(), sumPtALaRC, 0, 0, isBadJetALaRC);
2907 for(Int_t it=0; it<fListALaRC->GetSize(); ++it){ // loop for Lambdas in random cone
2909 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListALaRC->At(it));
2912 Double_t invMALa =0;
2917 CalculateInvMass(v0, kAntiLambda, invMALa, trackPt); //function to calculate invMass with TLorentzVector class
2919 Double_t vALaRC[3] = {invMALa,trackPt,fEta};
2920 fhnALaRC->Fill(vALaRC);
2930 //fetch particles in perpendicular cone to estimate UE event contribution to particle spectrum
2931 //these perpendicular cone particle spectra serve to subtract the particles in jet cones, that are stemming from the Underlying event, on a statistical basis
2932 //for normalization the common jet pT spectrum is used: fh1IMK0Cone, fh1IMLaCone and fh1IMALaCone
2934 //____fetch reconstructed K0s in cone perpendicular to jet axis:_______________________________________________________________________________
2936 TList* jetPerpConeK0list = new TList();
2938 Double_t sumPerpPtK0 = 0.;
2940 GetTracksInPerpCone(fListK0s, jetPerpConeK0list, jet, GetFFRadius(), sumPerpPtK0); //reconstructed K0s in cone around jet axis
2942 if(fDebug>2)Printf("%s:%d nK0s total: %d, in perp jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nK0s,jetPerpConeK0list->GetEntries(),GetFFRadius());
2944 for(Int_t it=0; it<jetPerpConeK0list->GetSize(); ++it){ // loop for K0s in perpendicular cone
2946 AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetPerpConeK0list->At(it));
2949 Double_t invMPerpK0s =0;
2954 CalculateInvMass(v0, kK0, invMPerpK0s, trackPt); //function to calculate invMass with TLorentzVector class
2955 Double_t vK0sPC[4] = {jetPt, invMPerpK0s,trackPt,fEta};
2957 fhnK0sPC->Fill(vK0sPC); //(x,y,z) //pay attention, this histogram contains the V0 content of both (+/- 90 degrees) perp. cones!!
2962 if(jetPerpConeK0list->GetSize() == 0){ // no K0s in jet cone
2964 Double_t vK0sPC[4] = {jetPt, -1, -1 , -999};//default values for case: no K0s is found in PC
2965 fhnK0sPC->Fill(vK0sPC);
2969 if(ij==0){//median cluster only once for event
2971 AliAODJet* medianCluster = GetMedianCluster();
2974 // ____ rec K0s in median cluster___________________________________________________________________________________________________________
2976 TList* jetMedianConeK0list = new TList();
2977 TList* jetMedianConeLalist = new TList();
2978 TList* jetMedianConeALalist = new TList();
2981 Double_t medianEta = medianCluster->Eta();
2983 if(TMath::Abs(medianEta)<=fCutjetEta){
2985 fh1MedianEta->Fill(medianEta);
2986 fh1JetPtMedian->Fill(jetPt); //for normalisation by total number of median cluster jets
2988 Double_t sumMedianPtK0 = 0.;
2990 Bool_t isBadJetK0Median = kFALSE; // dummy, do not use
2992 GetTracksInCone(fListK0s, jetMedianConeK0list, medianCluster, GetFFRadius(), sumMedianPtK0, 0., 0., isBadJetK0Median); //reconstructed K0s in median cone around jet axis
2993 //GetTracksInCone(fListK0s, jetConeK0list, jet, GetFFRadius(), sumPtK0, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetK0); //original use of function
2995 //cut parameters from Fragmentation Function task:
2996 //Float_t fFFMinLTrackPt; // reject jets with leading track with pt smaller than this value, use GetFFMinLTrackPt()
2997 //Float_t fFFMaxTrackPt; // reject jetscontaining any track with pt larger than this value, use GetFFMaxTrackPt()
2999 for(Int_t it=0; it<jetMedianConeK0list->GetSize(); ++it){ // loop for K0s in median cone
3001 AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetMedianConeK0list->At(it));
3004 Double_t invMMedianK0s =0;
3009 CalculateInvMass(v0, kK0, invMMedianK0s, trackPt); //function to calculate invMass with TLorentzVector class
3010 Double_t vK0sMCC[3] = {invMMedianK0s,trackPt,fEta};
3011 fhnK0sMCC->Fill(vK0sMCC);
3015 if(jetMedianConeK0list->GetSize() == 0){ // no K0s in median cluster cone
3017 Double_t vK0sMCC[3] = {-1, -1, -999};
3018 fhnK0sMCC->Fill(vK0sMCC);
3022 //__________________________________________________________________________________________________________________________________________
3023 // ____ rec Lambdas in median cluster___________________________________________________________________________________________________________
3025 Double_t sumMedianPtLa = 0.;
3026 Bool_t isBadJetLaMedian = kFALSE; // dummy, do not use
3028 GetTracksInCone(fListLa, jetMedianConeLalist, medianCluster, GetFFRadius(), sumMedianPtLa, 0, 0, isBadJetLaMedian); //reconstructed Lambdas in median cone around jet axis
3030 //cut parameters from Fragmentation Function task:
3031 //Float_t fFFMinLTrackPt; // reject jets with leading track with pt smaller than this value, use GetFFMinLTrackPt()
3032 //Float_t fFFMaxTrackPt; // reject jets containing any track with pt larger than this value, use GetFFMaxTrackPt()
3034 for(Int_t it=0; it<jetMedianConeLalist->GetSize(); ++it){ // loop for Lambdas in perpendicular cone
3036 AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetMedianConeLalist->At(it));
3039 Double_t invMMedianLa =0;
3044 CalculateInvMass(v0, kLambda, invMMedianLa, trackPt); //function to calculate invMass with TLorentzVector class
3046 Double_t vLaMCC[4] = {jetPt, invMMedianLa,trackPt,fEta};
3047 fhnLaMCC->Fill(vLaMCC);
3050 if(jetMedianConeLalist->GetSize() == 0){ // no Lambdas in median cluster cone
3052 Double_t vLaMCC[4] = {jetPt, -1, -1, -999};
3053 fhnLaMCC->Fill(vLaMCC);
3058 // ____ rec Antilambdas in median cluster___________________________________________________________________________________________________________
3061 Double_t sumMedianPtALa = 0.;
3063 Bool_t isBadJetALaMedian = kFALSE; // dummy, do not use
3065 GetTracksInCone(fListALa, jetMedianConeALalist, medianCluster, GetFFRadius(), sumMedianPtALa, 0, 0, isBadJetALaMedian); //reconstructed Antilambdas in median cone around jet axis
3068 //cut parameters from Fragmentation Function task:
3069 //Float_t fFFMinLTrackPt; // reject jets with leading track with pt smaller than this value, use GetFFMinLTrackPt()
3070 //Float_t fFFMaxTrackPt; // reject jets containing any track with pt larger than this value, use GetFFMaxTrackPt()
3072 for(Int_t it=0; it<jetMedianConeALalist->GetSize(); ++it){ // loop for Antilambdas in median cluster cone
3074 AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetMedianConeALalist->At(it));
3077 Double_t invMMedianALa =0;
3083 CalculateInvMass(v0, kAntiLambda, invMMedianALa, trackPt); //function to calculate invMass with TLorentzVector class
3084 Double_t vALaMCC[4] = {jetPt, invMMedianALa,trackPt,fEta};
3085 fhnALaMCC->Fill(vALaMCC);
3089 if(jetMedianConeALalist->GetSize() == 0){ // no Antilambdas in median cluster cone
3091 Double_t vALaMCC[4] = {jetPt, -1, -1, -999};
3092 fhnALaMCC->Fill(vALaMCC);
3095 }//median cluster eta cut
3097 delete jetMedianConeK0list;
3098 delete jetMedianConeLalist;
3099 delete jetMedianConeALalist;
3101 }//if mediancluster is existing
3103 //_________________________________________________________________________________________________________________________________________
3105 //____fetch reconstructed Lambdas in cone perpendicular to jet axis:__________________________________________________________________________
3107 TList* jetPerpConeLalist = new TList();
3109 Double_t sumPerpPtLa = 0.;
3111 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!!
3113 if(fDebug>2)Printf("%s:%d nLa total: %d, in perp jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nLa,jetPerpConeLalist->GetEntries(),GetFFRadius());
3115 for(Int_t it=0; it<jetPerpConeLalist->GetSize(); ++it){ // loop for Lambdas in perpendicular cone
3117 AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetPerpConeLalist->At(it));
3120 Double_t invMPerpLa =0;
3125 CalculateInvMass(v0, kLambda, invMPerpLa, trackPt); //function to calculate invMass with TLorentzVector class
3126 Double_t vLaPC[4] = {jetPt, invMPerpLa,trackPt,fEta};
3127 fhnLaPC->Fill(vLaPC); //(x,y,z) //pay attention, this histogram contains the V0 content of both (+/- 90 degrees) perp. cones!!
3132 if(jetPerpConeLalist->GetSize() == 0){ // no Lambdas in jet
3134 Double_t vLaPC[4] = {jetPt, -1, -1 , -999};//default values for case: no K0s is found in PC
3135 fhnLaPC->Fill(vLaPC);
3141 //____fetch reconstructed Antilambdas in cone perpendicular to jet axis:___________________________________________________________________
3143 TList* jetPerpConeALalist = new TList();
3145 Double_t sumPerpPtALa = 0.;
3147 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!!
3149 if(fDebug>2)Printf("%s:%d nALa total: %d, in perp jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nALa,jetPerpConeALalist->GetEntries(),GetFFRadius());
3151 for(Int_t it=0; it<jetPerpConeALalist->GetSize(); ++it){ // loop for ALa in perpendicular cone
3153 AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetPerpConeALalist->At(it));
3156 Double_t invMPerpALa =0;
3161 CalculateInvMass(v0, kAntiLambda, invMPerpALa, trackPt); //function to calculate invMass with TLorentzVector class
3162 Double_t vALaPC[4] = {jetPt, invMPerpALa,trackPt,fEta};
3163 fhnALaPC->Fill(vALaPC);
3168 if(jetPerpConeALalist->GetSize() == 0){ // no Antilambda
3170 Double_t vALaPC[4] = {jetPt, -1, -1, -999};
3171 fhnALaPC->Fill(vALaPC);
3191 //###########################################################################################################
3193 //__________________________________________________________________________________________________________________________________________
3197 //fill feeddown candidates from TList
3198 //std::cout<<"fListFeeddownLaCand entries: "<<fListFeeddownLaCand->GetSize()<<std::endl;
3200 Double_t sumPtFDLa = 0.;
3201 Bool_t isBadJetFDLa = kFALSE; // dummy, do not use
3203 GetTracksInCone(fListFeeddownLaCand, jetConeFDLalist, jet, GetFFRadius(), sumPtFDLa, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetFDLa);
3205 Double_t sumPtFDALa = 0.;
3206 Bool_t isBadJetFDALa = kFALSE; // dummy, do not use
3208 GetTracksInCone(fListFeeddownALaCand, jetConeFDALalist, jet, GetFFRadius(), sumPtFDALa, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetFDALa);
3210 //_________________________________________________________________
3211 for(Int_t it=0; it<fListFeeddownLaCand->GetSize(); ++it){
3213 AliAODv0* mcfd = dynamic_cast<AliAODv0*>(fListFeeddownLaCand->At(it));
3216 Double_t invMLaFDcand = 0;
3217 Double_t trackPt = 0;//pt of ass. particle, not used for the histos
3219 CalculateInvMass(mcfd, kLambda, invMLaFDcand, trackPt);
3221 //Get MC gen. Lambda transverse momentum
3222 TClonesArray *st = 0x0;
3225 TList *lt = fAOD->GetList();
3228 st = (TClonesArray*)lt->FindObject(AliAODMCParticle::StdBranchName()); //get MCAOD branch in data
3231 AliAODTrack *daughtertrack = (AliAODTrack *) (mcfd->GetSecondaryVtx()->GetDaughter(0));//fetch the first of the two daughter tracks
3232 Int_t AssLabel = TMath::Abs(daughtertrack->GetLabel());
3234 AliAODMCParticle *mcDaughterPart =(AliAODMCParticle*)st->UncheckedAt(AssLabel);
3236 Int_t v0lab = mcDaughterPart->GetMother();
3238 // Int_t v0lab= TMath::Abs(mcfd->GetLabel());//GetLabel doesn't work for AliAODv0 class!!! Only for AliAODtrack
3240 if((!v0lab) || (v0lab<0) || (v0lab > st->GetEntriesFast()))continue;//validity checks
3242 AliAODMCParticle *mcp=(AliAODMCParticle*)st->UncheckedAt(v0lab);
3244 Double_t genLaPt = mcp->Pt();
3246 //std::cout<<"Incl FD, genLaPt:"<<genLaPt<<std::endl;
3248 Double_t vFeedDownLa[3] = {5., invMLaFDcand, genLaPt};
3249 fhnFeedDownLa->Fill(vFeedDownLa);
3252 }//end loop over feeddown candidates for Lambda particles in jet cone
3253 //fetch MC truth in jet cones, denominator of rec. efficiency in jet cones
3254 //_________________________________________________________________
3255 for(Int_t it=0; it<jetConeFDLalist->GetSize(); ++it){
3257 AliAODv0* mcfd = dynamic_cast<AliAODv0*>(jetConeFDLalist->At(it));
3260 //std::cout<<"Cone, recLaPt:"<<mcfd->Pt()<<std::endl;
3262 Double_t invMLaFDcand = 0;
3263 Double_t trackPt = mcfd->Pt();//pt of ass. particle, not used for the histos
3265 CalculateInvMass(mcfd, kLambda, invMLaFDcand, trackPt);
3267 //Get MC gen. Lambda transverse momentum
3268 TClonesArray *st = 0x0;
3270 TList *lt = fAOD->GetList();
3273 st = (TClonesArray*)lt->FindObject(AliAODMCParticle::StdBranchName());
3275 AliAODTrack *daughtertrack = (AliAODTrack *) (mcfd->GetSecondaryVtx()->GetDaughter(0));//fetch the first of the two daughter tracks
3276 Int_t AssLabel = TMath::Abs(daughtertrack->GetLabel());
3278 AliAODMCParticle *mcDaughterPart =(AliAODMCParticle*)st->UncheckedAt(AssLabel);
3280 Int_t v0lab = mcDaughterPart->GetMother();
3282 //std::cout<<"v0lab: "<<v0lab<<std::endl;
3284 if((!v0lab) || (v0lab<0) || (v0lab > st->GetEntriesFast()))continue;//validity checks
3286 AliAODMCParticle *mcp=(AliAODMCParticle*)st->UncheckedAt(v0lab);
3288 Double_t genLaPt = mcp->Pt();
3291 //std::cout<<"Cone FD, genLaPt:"<<genLaPt<<std::endl;
3293 Double_t vFeedDownLaCone[3] = {jetPt, invMLaFDcand, genLaPt};
3294 fhnFeedDownLaCone->Fill(vFeedDownLaCone);
3297 }//end loop over feeddown candidates for Lambda particles in jet cone
3299 //_________________________________________________________________
3300 for(Int_t it=0; it<fListFeeddownALaCand->GetSize(); ++it){
3302 AliAODv0* mcfd = dynamic_cast<AliAODv0*>(fListFeeddownALaCand->At(it));
3305 Double_t invMALaFDcand = 0;
3306 Double_t trackPt = 0;//pt of ass. particle, not used for the histos
3308 CalculateInvMass(mcfd, kAntiLambda, invMALaFDcand, trackPt);
3310 //Get MC gen. Antilambda transverse momentum
3311 TClonesArray *st = 0x0;
3313 TList *lt = fAOD->GetList();
3316 st = (TClonesArray*)lt->FindObject(AliAODMCParticle::StdBranchName());
3318 AliAODTrack *daughtertrack = (AliAODTrack *) (mcfd->GetSecondaryVtx()->GetDaughter(0));//fetch the first of the two daughter tracks
3319 Int_t AssLabel = TMath::Abs(daughtertrack->GetLabel());
3321 AliAODMCParticle *mcDaughterPart =(AliAODMCParticle*)st->UncheckedAt(AssLabel);
3323 Int_t v0lab = mcDaughterPart->GetMother();
3326 if((!v0lab) || (v0lab<0) || (v0lab > st->GetEntriesFast()))continue;//validity checks
3328 AliAODMCParticle *mcp=(AliAODMCParticle*)st->UncheckedAt(v0lab);
3330 Double_t genALaPt = mcp->Pt();
3332 Double_t vFeedDownALa[3] = {5., invMALaFDcand, genALaPt};
3333 fhnFeedDownALa->Fill(vFeedDownALa);
3336 }//end loop over feeddown candidates for Antilambda particles
3338 //_________________________________________________________________
3339 //feeddown for Antilambdas from Xi(bar)+ and Xi(bar)0 in jet cone:
3341 for(Int_t it=0; it<jetConeFDALalist->GetSize(); ++it){
3343 AliAODv0* mcfd = dynamic_cast<AliAODv0*>(jetConeFDALalist->At(it));
3346 Double_t invMALaFDcand = 0;
3347 Double_t trackPt = 0;//pt of ass. particle, not used for the histos
3349 CalculateInvMass(mcfd, kAntiLambda, invMALaFDcand, trackPt);
3351 //Get MC gen. Antilambda transverse momentum
3352 TClonesArray *st = 0x0;
3354 TList *lt = fAOD->GetList();
3357 st = (TClonesArray*)lt->FindObject(AliAODMCParticle::StdBranchName());
3359 AliAODTrack *daughtertrack = (AliAODTrack *) (mcfd->GetSecondaryVtx()->GetDaughter(0));//fetch the first of the two daughter tracks
3360 Int_t AssLabel = TMath::Abs(daughtertrack->GetLabel());
3362 AliAODMCParticle *mcDaughterPart =(AliAODMCParticle*)st->UncheckedAt(AssLabel);
3364 Int_t v0lab = mcDaughterPart->GetMother();
3366 if((!v0lab) || (v0lab<0) || (v0lab > st->GetEntriesFast()))continue;//validity checks
3368 AliAODMCParticle *mcp=(AliAODMCParticle*)st->UncheckedAt(v0lab);
3370 Double_t genALaPt = mcp->Pt();
3372 Double_t vFeedDownALaCone[3] = {jetPt, invMALaFDcand, genALaPt};
3373 fhnFeedDownALaCone->Fill(vFeedDownALaCone);
3376 }//end loop over feeddown candidates for Antilambda particles in jet cone
3380 //____fetch MC generated K0s in cone around jet axis__(note: particles can stem from fragmentation but also from underlying event)________
3382 Double_t sumPtMCgenK0s = 0.;
3383 Bool_t isBadJetMCgenK0s = kFALSE; // dummy, do not use
3386 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!!)
3388 //first: sampling MC gen K0s
3390 GetTracksInCone(fListMCgenK0s, fListMCgenK0sCone, jet, GetFFRadius(), sumPtMCgenK0s, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetMCgenK0s); //MC generated K0s in cone around jet axis
3392 if(fDebug>2)Printf("%s:%d nMCgenK0s in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,fListMCgenK0sCone->GetEntries(),GetFFRadius());
3395 for(Int_t it=0; it<fListMCgenK0sCone->GetSize(); ++it){ // loop MC generated K0s in cone around jet axis
3397 AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenK0sCone->At(it));
3400 //Double_t fRapMCgenK0s = MyRapidity(mcp0->E(),mcp0->Pz());//get rec. particle in cone information
3401 Double_t fEtaMCgenK0s = mcp0->Eta();
3402 Double_t fPtMCgenK0s = mcp0->Pt();
3404 fh2MCgenK0Cone->Fill(jetPt,fPtMCgenK0s);
3405 fh2MCEtagenK0Cone->Fill(jetPt,fEtaMCgenK0s);
3409 //check whether the reconstructed K0s in jet cone are stemming from MC gen K0s (on MCgenK0s list):__________________________________________________
3411 for(Int_t ic=0; ic<jetConeK0list->GetSize(); ++ic){ //loop over all reconstructed K0s in jet cone
3413 //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
3415 Int_t negDaughterpdg;
3416 Int_t posDaughterpdg;
3419 Double_t fPtMCrecK0Match;
3420 Double_t invMK0Match;
3424 Bool_t fPhysicalPrimary = -1;
3425 Int_t MCv0PDGCode =0;
3426 Double_t jetPtSmear = -1;
3428 AliAODv0* v0c = dynamic_cast<AliAODv0*>(jetConeK0list->At(ic));//pointer to reconstructed K0s inside jet cone (cone is placed around reconstructed jet axis)
3430 //AliAODv0* v0c = dynamic_cast<AliAODv0*>(fListK0s->At(ic));//pointer to reconstructed K0s
3433 Bool_t daughtercheck = DaughterTrackCheck(v0c, nnum, pnum);//check daughter tracks have proper sign
3434 if(daughtercheck == kFALSE)continue;
3436 const AliAODTrack *trackMCNeg=(AliAODTrack *)(v0c->GetDaughter(nnum));
3437 const AliAODTrack *trackMCPos=(AliAODTrack *)(v0c->GetDaughter(pnum));
3439 TString generatorName;
3440 TList *listmc = fAOD->GetList();
3442 Bool_t mclabelcheck = MCLabelCheck(v0c, kK0, trackMCNeg, trackMCPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PDGCode, generatorName, isinjected);
3444 if(mclabelcheck == kFALSE)continue;
3445 if(fPhysicalPrimary == kFALSE)continue; //requirements for rec. V0 associated to MC true primary particle
3447 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
3449 //for(Int_t it=0; it<fListMCgenK0sCone->GetSize(); ++it){//belongs to previous definition of rec. eff. of V0s within jet cone
3451 //Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
3452 //AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenK0sCone->At(it));
3453 AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenK0s->At(it));
3456 Bool_t particleMatching = IsParticleMatching(mcp0, v0Label);
3458 if(particleMatching == kFALSE)continue; //if reconstructed V0 particle doesn't match to the associated MC particle go to next stack entry
3459 CalculateInvMass(v0c, kK0, invMK0Match, fPtMCrecK0Match);
3460 Double_t fEta = v0c->Eta();
3461 Double_t fPtMCgenK0s = mcp0->Pt();//pt has to be always MC truth value!
3463 Double_t vMCrecK0Cone[4] = {jetPt, invMK0Match,fPtMCgenK0s,fEta};
3464 fhnMCrecK0Cone->Fill(vMCrecK0Cone); //fill matching rec. K0s in 3D histogram
3466 SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear); //jetPt, cent, jetRadius, ptmintrack, &jetPtSmear
3468 Double_t vMCrecK0ConeSmear[4] = {jetPtSmear, invMK0Match,fPtMCgenK0s,fEta};
3469 fhnMCrecK0ConeSmear->Fill(vMCrecK0ConeSmear);
3471 //fill matching rec. K0s in 3D histogram, jet pT smeared according to deltaptjet distribution width
3474 } // end MCgenK0s / MCgenK0sCone loop
3477 //check the K0s daughters contamination of the jet tracks:
3479 TClonesArray *stackMC = 0x0;
3481 for(Int_t it=0; it<jettracklist->GetSize(); ++it){//loop over all particles in the jet
3483 AliVParticle* trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));//all tracks in jet cone
3484 if(!trackVP)continue;
3485 AliAODTrack *tr = dynamic_cast<AliAODTrack*> (trackVP); //fetch one jet track from the TList
3488 //get MC label information
3489 TList *mclist = fAOD->GetList(); //fetch the MC stack
3490 if(!mclist)continue;
3492 stackMC = (TClonesArray*)mclist->FindObject(AliAODMCParticle::StdBranchName()); //get MCAOD branch in data
3493 if (!stackMC) {Printf("ERROR: stack not available");}
3496 Int_t particleLabel = TMath::Abs(tr->GetLabel()); //fetch jet track label in MC stack
3498 //v0c is pointer to K0s candidate, is fetched already above, here it is just checked again whether daughters are properly ordered by their charge
3500 Bool_t daughterchecks = DaughterTrackCheck(v0c, nnum, pnum);
3502 if(daughterchecks == kFALSE)continue; //make sure that daughters are properly ordered
3504 const AliAODTrack *trackNeg=(AliAODTrack *)(v0c->GetDaughter(nnum)); //fetch v0 daughters of reconstructed K0s
3505 const AliAODTrack *trackPos=(AliAODTrack *)(v0c->GetDaughter(pnum));
3507 if(!trackNeg)continue;
3508 if(!trackPos)continue;
3510 Int_t negAssLabel = TMath::Abs(trackNeg->GetLabel()); //negative (reconstructed) charged track label in MC stack
3511 Int_t posAssLabel = TMath::Abs(trackPos->GetLabel()); //positive (reconstructed) charged track label in MC stack
3514 if(particleLabel == posAssLabel){ //check whether jet track and each of the rec. K0s daughters have same MC stack label -> are identical
3515 AliAODMCParticle* mctrackPos = dynamic_cast<AliAODMCParticle*>(stackMC->At(posAssLabel));
3516 if(!mctrackPos) continue;
3517 Double_t trackPosPt = mctrackPos->Pt();
3518 Double_t trackPosEta = mctrackPos->Eta();
3520 Double_t vK0sSecContinCone[3] = {jetPt, trackPosPt, trackPosEta};
3521 fhnK0sSecContinCone->Fill(vK0sSecContinCone);} //if it's the case, fill jet pt, daughter track pt and track eta in histo
3523 if(particleLabel == negAssLabel){
3524 AliAODMCParticle* mctrackNeg = dynamic_cast<AliAODMCParticle*>(stackMC->At(negAssLabel));
3525 if(!mctrackNeg) continue;
3526 Double_t trackNegPt = mctrackNeg->Pt();
3527 Double_t trackNegEta = mctrackNeg->Eta();
3529 Double_t vK0sSecContinCone[3] = {jetPt, trackNegPt, trackNegEta};
3530 fhnK0sSecContinCone->Fill(vK0sSecContinCone);} //if it's the case, fill jet pt, daughter track pt and track eta in histo
3538 } //end rec-K0-in-cone loop
3540 //________________________________________________________________________________________________________________________________________________________
3542 delete fListMCgenK0sCone;
3547 delete jetConeK0list;
3548 delete jetPerpConeK0list;
3549 delete jetPerpConeLalist;
3550 delete jetPerpConeALalist;
3553 //---------------La--------------------------------------------------------------------------------------------------------------------------------------------
3555 // fQAJetHistosRecCutsLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt() );
3557 for(Int_t it=0; it<fListLa->GetSize(); ++it){ // loop all La
3559 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListLa->At(it));
3564 TVector3 v0MomVect(v0Mom);
3566 Double_t dPhiJetLa = (jet->MomentumVector()->Vect()).DeltaPhi(v0MomVect);
3571 CalculateInvMass(v0, kLambda, invMLa, trackPt); //function to calculate invMass with TLorentzVector class
3572 // Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
3574 //if(incrementJetPt){fh1V0JetPt->Fill(jetPt);}
3576 //fFFHistosIMLaJet->FillFF(trackPt, invMLa, jetPt, incrementJetPt);
3578 if(dPhiJetLa<fh1dPhiJetLa->GetXaxis()->GetXmin()) dPhiJetLa += 2*TMath::Pi();
3579 fh1dPhiJetLa->Fill(dPhiJetLa);
3582 /* if(fListLa->GetSize() == 0){ // no La: increment jet pt spectrum
3584 Bool_t incrementJetPt = kTRUE;
3585 fFFHistosIMLaJet->FillFF(-1, -1, jetPt, incrementJetPt);
3589 // ____fetch rec. Lambdas in cone around jet axis_______________________________________________________________________________________
3591 TList* jetConeLalist = new TList();
3592 Double_t sumPtLa = 0.;
3593 Bool_t isBadJetLa = kFALSE; // dummy, do not use
3595 GetTracksInCone(fListLa, jetConeLalist, jet, GetFFRadius(), sumPtLa, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetLa);//method inherited from FF
3597 if(fDebug>2)Printf("%s:%d nLa total: %d, in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nLa,jetConeLalist->GetEntries(),GetFFRadius());
3599 for(Int_t it=0; it<jetConeLalist->GetSize(); ++it){ // loop La in jet cone
3601 AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetConeLalist->At(it));
3607 Bool_t daughtercheck = DaughterTrackCheck(v0, nnum, pnum);
3608 if(daughtercheck == kFALSE)continue;
3615 CalculateInvMass(v0, kLambda, invMLa, trackPt); //function to calculate invMass with TLorentzVector class
3617 Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;//needed for all histos, which serve for normalisation
3621 Int_t negDaughterpdg;
3622 Int_t posDaughterpdg;
3625 Double_t jetPtSmear = -1;
3627 Bool_t fPhysicalPrimary = -1;
3628 Int_t MCv0PDGCode =0;
3629 TString generatorName;
3631 SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
3632 if(incrementJetPt == kTRUE){fh1IMLaConeSmear->Fill(jetPtSmear);
3634 const AliAODTrack *trackNeg=(AliAODTrack *)(v0->GetDaughter(nnum));
3635 const AliAODTrack *trackPos=(AliAODTrack *)(v0->GetDaughter(pnum));
3637 TList *listmc = fAOD->GetList();
3639 Bool_t mclabelcheck = MCLabelCheck(v0, kLambda, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PDGCode, generatorName, isinjected);
3640 if(mclabelcheck == kFALSE)continue;
3642 //std::cout<<"generatorName: "<<generatorName<<std::endl;
3644 if(generatorName == "Hijing"){
3645 Double_t vrecMCHijingLaCone[4] = {jetPt, invMLa,trackPt,fEta};
3646 fhnrecMCHijingLaCone->Fill(vrecMCHijingLaCone);
3649 if(isinjected == kTRUE){
3650 Double_t vrecMCInjectLaCone[4] = {jetPt, invMLa,trackPt,fEta};
3651 fhnrecMCInjectLaCone->Fill(vrecMCInjectLaCone);
3654 }//fill TH1F for normalization purposes
3657 if(incrementJetPt==kTRUE){
3658 fh1IMLaCone->Fill(jetPt);}//normalisation by number of selected jets
3660 //fFFHistosIMLaCone->FillFF(trackPt, invMLa, jetPt, incrementJetPt);
3664 if(jetConeLalist->GetSize() == 0){ // no La: increment jet pt spectrum
3666 Bool_t incrementJetPt = kTRUE;
3667 // fFFHistosIMLaCone->FillFF(-1, -1, jetPt, incrementJetPt);
3668 Double_t vLaCone[4] = {jetPt, -1, -1, -1};
3669 fhnLaCone->Fill(vLaCone);
3671 if(incrementJetPt==kTRUE){
3672 fh1IMLaCone->Fill(jetPt);}//normalisation by number of selected jets
3675 Double_t jetPtSmear;
3676 SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
3677 if(incrementJetPt == kTRUE){
3678 fh1IMLaConeSmear->Fill(jetPtSmear);
3687 //____fetch MC generated Lambdas in cone around jet axis__(particles can stem from fragmentation but also from underlying event)_____________
3689 Double_t sumPtMCgenLa = 0.;
3690 Bool_t isBadJetMCgenLa = kFALSE; // dummy, do not use
3692 //sampling MC gen. Lambdas in cone around reconstructed jet axis
3693 fListMCgenLaCone = new TList();
3695 GetTracksInCone(fListMCgenLa, fListMCgenLaCone, jet, GetFFRadius(), sumPtMCgenLa, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetMCgenLa);//fetch MC generated Lambdas in cone of resolution parameter R around jet axis
3697 if(fDebug>2)Printf("%s:%d nMCgenLa in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,fListMCgenLaCone->GetEntries(),GetFFRadius());
3699 for(Int_t it=0; it<fListMCgenLaCone->GetSize(); ++it){ // loop MC generated La in cone around jet axis
3701 AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenLaCone->At(it));
3704 //Double_t fRapMCgenLa = MyRapidity(mcp0->E(),mcp0->Pz());
3705 Double_t fEtaMCgenLa = mcp0->Eta();
3706 Double_t fPtMCgenLa = mcp0->Pt();
3708 fh2MCgenLaCone->Fill(jetPt,fPtMCgenLa);
3709 fh2MCEtagenLaCone->Fill(jetPt,fEtaMCgenLa);
3713 //check whether the reconstructed La are stemming from MC gen La on fListMCgenLa List:__________________________________________________
3715 for(Int_t ic=0; ic<jetConeLalist->GetSize(); ++ic){//loop over all reconstructed La within jet cone, new definition
3717 //for(Int_t ic=0; ic<fListLa->GetSize(); ++ic){//old definition
3719 Int_t negDaughterpdg;
3720 Int_t posDaughterpdg;
3723 Double_t fPtMCrecLaMatch;
3724 Double_t invMLaMatch;
3728 Bool_t fPhysicalPrimary = -1;
3729 Int_t MCv0PDGCode =0;
3730 Double_t jetPtSmear = -1;
3731 TString generatorName;
3733 AliAODv0* v0c = dynamic_cast<AliAODv0*>(jetConeLalist->At(ic));//new definition
3738 Bool_t daughtercheck = DaughterTrackCheck(v0c, nnum, pnum);
3739 if(daughtercheck == kFALSE)continue;
3741 const AliAODTrack *trackMCNeg=(AliAODTrack *)(v0c->GetDaughter(nnum));
3742 const AliAODTrack *trackMCPos=(AliAODTrack *)(v0c->GetDaughter(pnum));
3744 TList *listmc = fAOD->GetList();
3746 Bool_t mclabelcheck = MCLabelCheck(v0c, kLambda, trackMCNeg, trackMCPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PDGCode, generatorName, isinjected);
3748 if(mclabelcheck == kFALSE)continue;
3749 if(fPhysicalPrimary == kFALSE)continue;
3751 for(Int_t it=0; it<fListMCgenLa->GetSize(); ++it){//new definition // loop over MC generated K0s in cone around jet axis
3753 // for(Int_t it=0; it<fListMCgenLaCone->GetSize(); ++it){//old definition // loop over MC generated La in cone around jet axis
3755 //Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
3757 AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenLa->At(it));//new definition
3758 //AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenLaCone->At(it));//old definition
3762 Bool_t particleMatching = IsParticleMatching(mcp0, v0Label);
3765 if(particleMatching == kFALSE)continue; //particle doesn't match on any associated MC gen particle in cone around rec jet axis
3767 CalculateInvMass(v0c, kLambda, invMLaMatch, fPtMCrecLaMatch);
3769 Double_t fPtMCgenLa = mcp0->Pt();
3770 Double_t fEta = v0c->Eta();//rec. MC particle
3771 Double_t vMCrecLaCone[4] = {jetPt, invMLaMatch,fPtMCgenLa,fEta};
3772 fhnMCrecLaCone->Fill(vMCrecLaCone);
3774 SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
3776 Double_t vMCrecLaConeSmear[4] = {jetPtSmear, invMLaMatch,fPtMCgenLa,fEta};
3777 fhnMCrecLaConeSmear->Fill(vMCrecLaConeSmear); //fill matching rec. Lambdas in 3D histogram, jet pT smeared according to deltaptjet distribution width
3780 } // end MCgenLa loop
3782 //check the Lambda daughters contamination of the jet tracks://///////////////////////////////////////////////////////////////////////////////////////////
3784 TClonesArray *stackMC = 0x0;
3786 for(Int_t it=0; it<jettracklist->GetSize(); ++it){//loop over all particles in the jet
3788 AliVParticle* trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));//all tracks in jet cone
3789 if(!trackVP)continue;
3790 AliAODTrack *tr = dynamic_cast<AliAODTrack*> (trackVP); //fetch one jet track from the TList
3793 //get MC label information
3794 TList *mclist = fAOD->GetList(); //fetch the MC stack
3796 stackMC = (TClonesArray*)mclist->FindObject(AliAODMCParticle::StdBranchName()); //get MCAOD branch in data
3797 if (!stackMC) {Printf("ERROR: stack not available");}
3800 Int_t particleLabel = TMath::Abs(tr->GetLabel()); //fetch jet track label in MC stack
3802 Bool_t daughterchecks = DaughterTrackCheck(v0c, nnum, pnum);
3804 if(daughterchecks == kFALSE)continue; //make sure that daughters are properly ordered
3806 const AliAODTrack *trackNeg=(AliAODTrack *)(v0c->GetDaughter(nnum)); //fetch v0 daughters of reconstructed K0s
3807 const AliAODTrack *trackPos=(AliAODTrack *)(v0c->GetDaughter(pnum));
3809 Int_t negAssLabel = TMath::Abs(trackNeg->GetLabel()); //negative (reconstructed) charged track label in MC stack
3810 Int_t posAssLabel = TMath::Abs(trackPos->GetLabel()); //positive (reconstructed) charged track label in MC stack
3813 if(particleLabel == posAssLabel){ //check whether jet track and each of the rec. K0s daughters have same MC stack label -> are identical
3815 AliAODMCParticle* mctrackPos = dynamic_cast<AliAODMCParticle*>(stackMC->At(posAssLabel));
3816 if(!mctrackPos) continue;
3818 Double_t trackPosPt = trackPos->Pt();
3819 Double_t trackPosEta = trackPos->Eta();
3820 Double_t vLaSecContinCone[3] = {jetPt, trackPosPt, trackPosEta};
3821 fhnLaSecContinCone->Fill(vLaSecContinCone);
3823 } //if it's the case, fill jet pt, daughter track pt and track eta in histo
3826 if(particleLabel == negAssLabel){
3828 AliAODMCParticle* mctrackNeg = dynamic_cast<AliAODMCParticle*>(stackMC->At(negAssLabel));
3829 if(!mctrackNeg) continue;
3831 Double_t trackNegPt = trackNeg->Pt();
3832 Double_t trackNegEta = trackNeg->Eta();
3834 Double_t vLaSecContinCone[3] = {jetPt, trackNegPt, trackNegEta};
3835 fhnLaSecContinCone->Fill(vLaSecContinCone);
3838 } //if it's the case, fill jet pt, daughter track pt and track eta in histo
3843 } //end rec-La-in-cone loop
3844 //________________________________________________________________________________________________________________________________________________________
3846 delete fListMCgenLaCone;
3850 delete jetConeLalist;
3854 //---------------ALa-----------
3857 // fQAJetHistosRecCutsLeading->FillJetQA( jet->Eta(), TVector2::Phi_0_2pi(jet->Phi()), jet->Pt() );
3859 for(Int_t it=0; it<fListALa->GetSize(); ++it){ // loop all ALa
3861 AliAODv0* v0 = dynamic_cast<AliAODv0*>(fListALa->At(it));
3866 TVector3 v0MomVect(v0Mom);
3868 Double_t dPhiJetALa = (jet->MomentumVector()->Vect()).DeltaPhi(v0MomVect);
3870 Double_t invMALa =0;
3873 CalculateInvMass(v0, kAntiLambda, invMALa, trackPt); //function to calculate invMass with TLorentzVector class
3874 //Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
3876 //if(incrementJetPt){fh1V0JetPt->Fill(jetPt);}
3878 //fFFHistosIMALaJet->FillFF(trackPt, invMALa, jetPt, incrementJetPt);
3880 if(dPhiJetALa<fh1dPhiJetALa->GetXaxis()->GetXmin()) dPhiJetALa += 2*TMath::Pi();
3881 fh1dPhiJetALa->Fill(dPhiJetALa);
3884 // if(fListALa->GetSize() == 0){ // no ALa: increment jet pt spectrum
3886 // Bool_t incrementJetPt = kTRUE;
3887 //fFFHistosIMALaJet->FillFF(-1, -1, jetPt, incrementJetPt);
3891 // ____fetch rec. Antilambdas in cone around jet axis_______________________________________________________________________________________
3893 TList* jetConeALalist = new TList();
3894 Double_t sumPtALa = 0.;
3895 Bool_t isBadJetALa = kFALSE; // dummy, do not use
3897 GetTracksInCone(fListALa, jetConeALalist, jet, GetFFRadius(), sumPtALa, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetALa);//method inherited from FF
3899 if(fDebug>2)Printf("%s:%d nALa total: %d, in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,nALa,jetConeALalist->GetEntries(),GetFFRadius());
3901 for(Int_t it=0; it<jetConeALalist->GetSize(); ++it){ // loop ALa in jet cone
3903 AliAODv0* v0 = dynamic_cast<AliAODv0*>(jetConeALalist->At(it));
3910 Bool_t daughtercheck = DaughterTrackCheck(v0, nnum, pnum);
3911 if(daughtercheck == kFALSE)continue;
3914 Double_t invMALa =0;
3920 CalculateInvMass(v0, kAntiLambda, invMALa, trackPt); //function to calculate invMass with TLorentzVector class
3922 Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
3924 if(fAnalysisMC){ //jet pt smearing study for Antilambdas
3926 Int_t negDaughterpdg;
3927 Int_t posDaughterpdg;
3930 Double_t jetPtSmear = -1;
3932 Bool_t fPhysicalPrimary = -1;
3933 Int_t MCv0PDGCode =0;
3934 TString generatorName;
3936 SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
3937 const AliAODTrack *trackNeg=(AliAODTrack *)(v0->GetDaughter(nnum));
3938 const AliAODTrack *trackPos=(AliAODTrack *)(v0->GetDaughter(pnum));
3940 TList *listmc = fAOD->GetList();
3942 Bool_t mclabelcheck = MCLabelCheck(v0, kAntiLambda, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PDGCode, generatorName, isinjected);
3943 if(mclabelcheck == kFALSE)continue;
3945 //std::cout<<"generatorName: "<<generatorName<<std::endl;
3947 if(generatorName == "Hijing"){
3948 Double_t vrecMCHijingALaCone[4] = {jetPt, invMALa,trackPt,fEta};
3949 fhnrecMCHijingALaCone->Fill(vrecMCHijingALaCone);
3952 if(isinjected == kTRUE){
3953 Double_t vrecMCInjectALaCone[4] = {jetPt, invMALa,trackPt,fEta};
3954 fhnrecMCInjectALaCone->Fill(vrecMCInjectALaCone);
3957 if(incrementJetPt == kTRUE){fh1IMALaConeSmear->Fill(jetPtSmear);} //fill TH1F for normalization purposes
3960 if(incrementJetPt==kTRUE){
3961 fh1IMALaCone->Fill(jetPt);}//normalisation by number of selected jets
3963 //fFFHistosIMALaCone->FillFF(trackPt, invMALa, jetPt, incrementJetPt);
3964 Double_t vALaCone[4] = {jetPt, invMALa,trackPt,fEta};
3965 fhnALaCone->Fill(vALaCone);
3968 if(jetConeALalist->GetSize() == 0){ // no ALa: increment jet pt spectrum
3970 Bool_t incrementJetPt = kTRUE;
3972 if(incrementJetPt==kTRUE){
3973 fh1IMALaCone->Fill(jetPt);}//normalisation by number of selected jets
3975 //fFFHistosIMALaCone->FillFF(-1, -1, jetPt, incrementJetPt);
3976 Double_t vALaCone[4] = {jetPt, -1, -1, -1};
3977 fhnALaCone->Fill(vALaCone);
3980 Double_t jetPtSmear;
3981 SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
3982 if(incrementJetPt == kTRUE)fh1IMALaConeSmear->Fill(jetPtSmear);}
3988 //____fetch MC generated Antilambdas in cone around jet axis__(particles can stem from fragmentation but also from underlying event)_____________
3990 Double_t sumPtMCgenALa = 0.;
3991 Bool_t isBadJetMCgenALa = kFALSE; // dummy, do not use
3993 //sampling MC gen Antilambdas in cone around reconstructed jet axis
3994 fListMCgenALaCone = new TList();
3996 GetTracksInCone(fListMCgenALa, fListMCgenALaCone, jet, GetFFRadius(), sumPtMCgenALa, GetFFMinLTrackPt(), GetFFMaxTrackPt(), isBadJetMCgenALa);//MC generated K0s in cone around jet axis
3998 if(fDebug>2)Printf("%s:%d nMCgenALa in jet cone: %d,FFRadius %f ",(char*)__FILE__,__LINE__,fListMCgenALaCone->GetEntries(),GetFFRadius());
4000 for(Int_t it=0; it<fListMCgenALaCone->GetSize(); ++it){ // loop MC generated La in cone around jet axis
4002 AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenALaCone->At(it));
4005 //Double_t fRapMCgenALa = MyRapidity(mcp0->E(),mcp0->Pz());
4006 Double_t fEtaMCgenALa = mcp0->Eta();
4007 Double_t fPtMCgenALa = mcp0->Pt();
4009 fh2MCgenALaCone->Fill(jetPt,fPtMCgenALa);
4010 fh2MCEtagenALaCone->Fill(jetPt,fEtaMCgenALa);
4014 //check whether the reconstructed ALa are stemming from MC gen ALa on MCgenALa List:__________________________________________________
4016 for(Int_t ic=0; ic<jetConeALalist->GetSize(); ++ic){//loop over all reconstructed ALa
4018 Int_t negDaughterpdg;
4019 Int_t posDaughterpdg;
4022 Double_t fPtMCrecALaMatch;
4023 Double_t invMALaMatch;
4027 Bool_t fPhysicalPrimary = -1;
4028 Int_t MCv0PDGCode =0;
4029 Double_t jetPtSmear = -1;
4030 TString generatorName;
4032 AliAODv0* v0c = dynamic_cast<AliAODv0*>(jetConeALalist->At(ic));
4035 Bool_t daughtercheck = DaughterTrackCheck(v0c, nnum, pnum);
4036 if(daughtercheck == kFALSE)continue;
4038 const AliAODTrack *trackMCNeg=(AliAODTrack *)(v0c->GetDaughter(nnum));
4039 const AliAODTrack *trackMCPos=(AliAODTrack *)(v0c->GetDaughter(pnum));
4041 TList *listmc = fAOD->GetList();
4042 if(!listmc)continue;
4044 Bool_t mclabelcheck = MCLabelCheck(v0c, kAntiLambda, trackMCNeg, trackMCPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PDGCode, generatorName, isinjected);
4046 if(mclabelcheck == kFALSE)continue;
4047 if(fPhysicalPrimary == kFALSE)continue;
4049 for(Int_t it=0; it<fListMCgenALa->GetSize(); ++it){ // loop over MC generated Antilambdas in cone around jet axis
4051 //Bool_t incrementJetPt = (it==0) ? kTRUE : kFALSE;
4053 AliAODMCParticle* mcp0 = dynamic_cast<AliAODMCParticle*>(fListMCgenALa->At(it));
4056 Bool_t particleMatching = IsParticleMatching(mcp0, v0Label);
4058 if(particleMatching == kFALSE)continue; //particle doesn't match on any associated MC gen particle in cone around rec jet axis
4060 CalculateInvMass(v0c, kAntiLambda, invMALaMatch, fPtMCrecALaMatch);
4062 Double_t fPtMCgenALa = mcp0->Pt();
4063 Double_t fEta = v0c->Eta();
4064 Double_t vMCrecALaCone[4] = {jetPt, invMALaMatch,fPtMCgenALa,fEta};
4065 fhnMCrecALaCone->Fill(vMCrecALaCone); //fill matching rec. Antilambda in 3D histogram
4067 SmearJetPt(jetPt,cl,GetFFRadius(),GetFFMinLTrackPt(),jetPtSmear);
4069 Double_t vMCrecALaConeSmear[4] = {jetPtSmear, invMALaMatch,fPtMCgenALa,fEta};
4070 fhnMCrecALaConeSmear->Fill(vMCrecALaConeSmear); //fill matching rec. Antilambda in 3D histogram
4072 } // end MCgenALa loop
4076 //check the Antilambda daughters contamination of the jet tracks:
4078 TClonesArray *stackMC = 0x0;
4080 for(Int_t it=0; it<jettracklist->GetSize(); ++it){//loop over all particles in the jet
4082 AliVParticle* trackVP = dynamic_cast<AliVParticle*>(jettracklist->At(it));//all tracks in jet cone
4083 if(!trackVP)continue;
4084 AliAODTrack *tr = dynamic_cast<AliAODTrack*> (trackVP); //fetch one jet track from the TList
4087 //get MC label information
4088 TList *mclist = fAOD->GetList(); //fetch the MC stack
4089 if(!mclist)continue;
4091 stackMC = (TClonesArray*)mclist->FindObject(AliAODMCParticle::StdBranchName()); //get MCAOD branch in data
4092 if (!stackMC) {Printf("ERROR: stack not available");}
4095 Int_t particleLabel = TMath::Abs(tr->GetLabel()); //fetch jet track label in MC stack
4097 Bool_t daughterchecks = DaughterTrackCheck(v0c, nnum, pnum);
4099 if(daughterchecks == kFALSE)continue; //make sure that daughters are properly ordered
4101 const AliAODTrack *trackNeg=(AliAODTrack *)(v0c->GetDaughter(nnum)); //fetch v0 daughters of reconstructed K0s
4102 const AliAODTrack *trackPos=(AliAODTrack *)(v0c->GetDaughter(pnum));
4103 if(!trackPos)continue;
4104 if(!trackNeg)continue;
4106 Int_t negAssLabel = TMath::Abs(trackNeg->GetLabel()); //negative (reconstructed) charged track label in MC stack
4107 Int_t posAssLabel = TMath::Abs(trackPos->GetLabel()); //positive (reconstructed) charged track label in MC stack
4109 if(!negAssLabel)continue;
4110 if(!posAssLabel)continue;
4112 if(particleLabel == posAssLabel){ //check whether jet track and each of the rec. K0s daughters have same MC stack label -> are identical
4113 AliAODMCParticle* mctrackPos = dynamic_cast<AliAODMCParticle*>(stackMC->At(posAssLabel));
4114 if(!mctrackPos) continue;
4116 Double_t trackPosPt = trackPos->Pt();
4117 Double_t trackPosEta = trackPos->Eta();
4118 if(!trackPosPt)continue;
4119 if(!trackPosEta)continue;
4121 Double_t vLaSecContinCone[3] = {jetPt, trackPosPt, trackPosEta};
4122 fhnLaSecContinCone->Fill(vLaSecContinCone);
4126 //fh3SecContinCone->Fill(jetPt, trackPosPt, trackPosEta);
4127 } //if it's the case, fill jet pt, daughter track pt and track eta in histo
4129 if(particleLabel == negAssLabel){
4131 AliAODMCParticle* mctrackNeg = dynamic_cast<AliAODMCParticle*>(stackMC->At(negAssLabel));
4132 if(!mctrackNeg) continue;
4134 Double_t trackNegPt = trackNeg->Pt();
4135 Double_t trackNegEta = trackNeg->Eta();
4137 if(!trackNegPt)continue;
4138 if(!trackNegEta)continue;
4140 Double_t vLaSecContinCone[3] = {jetPt, trackNegPt, trackNegEta};
4141 fhnLaSecContinCone->Fill(vLaSecContinCone);
4143 //fh3SecContinCone->Fill(jetPt, trackNegPt, trackNegEta);
4144 } //if it's the case, fill jet pt, daughter track pt and track eta in histo
4148 } //end rec-ALa-in-cone loop
4149 //________________________________________________________________________________________________________________________________________________________
4151 delete fListMCgenALaCone;
4155 delete jetConeALalist;
4156 delete jettracklist; //had been initialised at jet loop beginning
4158 }//end of if 'leading' or 'all jet' requirement
4164 fTracksRecCuts->Clear();
4165 fJetsRecCuts->Clear();
4166 fBckgJetsRec->Clear();
4170 fListFeeddownLaCand->Clear();
4171 fListFeeddownALaCand->Clear();
4172 jetConeFDLalist->Clear();
4173 jetConeFDALalist->Clear();
4174 fListMCgenK0s->Clear();
4175 fListMCgenLa->Clear();
4176 fListMCgenALa->Clear();
4180 PostData(1, fCommonHistList);
4185 // ____________________________________________________________________________________________
4186 void AliAnalysisTaskJetChem::SetProperties(TH3F* h,const char* x, const char* y, const char* z)
4188 //Set properties of histos (x,y and z title)
4193 h->GetXaxis()->SetTitleColor(1);
4194 h->GetYaxis()->SetTitleColor(1);
4195 h->GetZaxis()->SetTitleColor(1);
4199 //________________________________________________________________________________________________________________________________________
4200 Bool_t AliAnalysisTaskJetChem::AcceptBetheBloch(AliAODv0 *v0, AliPIDResponse *PIDResponse, const Int_t particletype) //dont use for MC Analysis
4206 const AliAODTrack *ntracktest=(AliAODTrack *)v0->GetDaughter(nnum);
4207 if(ntracktest->Charge() > 0){nnum = 0; pnum = 1;}
4209 const AliAODTrack *trackNeg=(AliAODTrack *)(v0->GetDaughter(nnum));
4210 const AliAODTrack *trackPos=(AliAODTrack *)(v0->GetDaughter(pnum));
4212 //Check if both tracks are available
4213 if (!trackPos || !trackNeg) {
4214 Printf("strange analysis::UserExec:: Error:Could not retrieve one of the daughter tracks\n");
4218 //remove like sign V0s
4219 if ( trackPos->Charge() == trackNeg->Charge() ){
4220 //if(fDebug>1) Printf("%s:%d found like-sign V0", (char*)__FILE__,__LINE__);
4225 Double_t nsig_p = 0; //number of sigmas that positive daughter track has got in TPC pid information
4226 Double_t nsig_n = 0;
4228 const AliAODPid *pid_p=trackPos->GetDetPid(); // returns fDetPID, more detailed or detector specific pid information
4229 const AliAODPid *pid_n=trackNeg->GetDetPid();
4231 if(!pid_p)return kFALSE;
4232 if(!pid_n)return kFALSE;
4236 if(particletype == 1) //PID cut on positive charged Lambda daughters (only those with pt < 1 GeV/c)
4239 nsig_p=PIDResponse->NumberOfSigmasTPC(trackPos,AliPID::kProton);
4240 Double_t protonPt = trackPos->Pt();
4241 if ((TMath::Abs(nsig_p) >= fCutBetheBloch) && (fCutBetheBloch >0) && (protonPt < 1)) return kFALSE;
4250 if(particletype == 2)
4252 nsig_n=PIDResponse->NumberOfSigmasTPC(trackNeg,AliPID::kProton);
4253 Double_t antiprotonPt = trackNeg->Pt();
4254 if ((TMath::Abs(nsig_n) >= fCutBetheBloch) && (fCutBetheBloch >0) && (antiprotonPt < 1)) return kFALSE;
4262 //___________________________________________________________________
4263 Bool_t AliAnalysisTaskJetChem::IsK0InvMass(const Double_t mass) const
4265 // K0 mass ? Use FF histo limits
4267 if(fFFIMInvMMin <= mass && mass < fFFIMInvMMax) return kTRUE;
4271 //___________________________________________________________________
4272 Bool_t AliAnalysisTaskJetChem::IsLaInvMass(const Double_t mass) const
4274 // La mass ? Use FF histo limits
4277 if(fFFIMLaInvMMin <= mass && mass < fFFIMLaInvMMax) return kTRUE;
4282 //_____________________________________________________________________________________
4283 Int_t AliAnalysisTaskJetChem::GetListOfV0s(TList *list, const Int_t type, const Int_t particletype, AliAODVertex* primVertex, AliAODEvent* aod)
4285 // fill list of V0s selected according to type
4288 if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
4293 if(fDebug>5){std::cout<<"AliAnalysisTaskJetChem::GetListOfV0s(): type: "<<type<<" particletype: "<<particletype<<"aod: "<<aod<<std::endl;
4294 if(type==kTrackUndef){std::cout<<"AliAnalysisTaskJetChem::GetListOfV0s(): kTrackUndef!! "<<std::endl;}
4298 if(type==kTrackUndef) return 0;
4300 if(!primVertex) return 0;
4302 Double_t lPrimaryVtxPosition[3];
4303 Double_t lV0Position[3];
4304 lPrimaryVtxPosition[0] = primVertex->GetX();
4305 lPrimaryVtxPosition[1] = primVertex->GetY();
4306 lPrimaryVtxPosition[2] = primVertex->GetZ();
4308 if(fDebug>5){ std::cout<<"AliAnalysisTaskJetChem::GetListOfV0s(): aod->GetNumberOfV0s: "<<aod->GetNumberOfV0s()<<std::endl; }
4311 for(int i=0; i<aod->GetNumberOfV0s(); i++){ // loop over V0s
4314 AliAODv0* v0 = aod->GetV0(i);
4318 std::cout << std::endl
4319 << "Warning in AliAnalysisTaskJetChem::GetListOfV0s:" << std::endl
4320 << "v0 = " << v0 << std::endl;
4324 Bool_t isOnFly = v0->GetOnFlyStatus();
4326 if(!isOnFly && (type == kOnFly || type == kOnFlyPID || type == kOnFlydEdx || type == kOnFlyPrim)) continue;
4327 if( isOnFly && (type == kOffl || type == kOfflPID || type == kOffldEdx || type == kOfflPrim)) continue;
4329 Int_t motherType = -1;
4330 //Double_t v0CalcMass = 0; //mass of MC v0
4331 Double_t MCPt = 0; //pt of MC v0
4333 Double_t pp[3]={0,0,0}; //3-momentum positive charged track
4334 Double_t pm[3]={0,0,0}; //3-momentum negative charged track
4335 Double_t v0mom[3]={0,0,0};
4346 Bool_t daughtercheck = DaughterTrackCheck(v0, nnum, pnum);
4348 if(daughtercheck == kFALSE)continue;
4350 const AliAODTrack *trackNeg=(AliAODTrack *)(v0->GetDaughter(nnum));
4351 const AliAODTrack *trackPos=(AliAODTrack *)(v0->GetDaughter(pnum));
4354 ///////////////////////////////////////////////////////////////////////////////////
4356 //calculate InvMass for every V0 particle assumption (Kaon=1,Lambda=2,Antilambda=3)
4357 switch(particletype){
4359 CalculateInvMass(v0, kK0, invM, trackPt); //function to calculate invMass with TLorentzVector class
4363 CalculateInvMass(v0, kLambda, invM, trackPt);
4367 CalculateInvMass(v0, kAntiLambda, invM, trackPt);
4371 std::cout<<"***NO VALID PARTICLETYPE***"<<std::endl;
4376 /////////////////////////////////////////////////////////////
4377 //V0 and track Cuts:
4380 if(fDebug>7){if(!(IsK0InvMass(invMK0s)) && !(IsLaInvMass(invMLa)) && !(IsLaInvMass(invMALa))){std::cout<<"AliAnalysisTaskJetChem::GetListOfV0s: invM not in selected mass window "<<std::endl;}}
4382 if(!(IsK0InvMass(invMK0s)) && !(IsLaInvMass(invMLa)) && !(IsLaInvMass(invMALa)))continue;
4384 // Double_t PosEta = trackPos->AliAODTrack::Eta();//daughter track charge is sometimes wrong here, account for that!!!
4385 // Double_t NegEta = trackNeg->AliAODTrack::Eta();
4387 Double_t PosEta = trackPos->Eta();//daughter track charge is sometimes wrong here, account for that!!!
4388 Double_t NegEta = trackNeg->Eta();
4390 Double_t PosCharge = trackPos->Charge();
4391 Double_t NegCharge = trackNeg->Charge();
4393 if((trackPos->Charge() == 1) && (trackNeg->Charge() == -1)) //Fill daughters charge into histo to check if they are symmetric distributed
4394 { fh1PosDaughterCharge->Fill(PosCharge);
4395 fh1NegDaughterCharge->Fill(NegCharge);
4398 //DistOverTotMom_in_2D___________
4400 Float_t fMassK0s = TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
4401 Float_t fMassLambda = TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
4404 AliAODVertex* primVtx = fAOD->GetPrimaryVertex(); // get the primary vertex
4405 Double_t dPrimVtxPos[3]; // primary vertex position {x,y,z}
4406 primVtx->GetXYZ(dPrimVtxPos);
4408 Float_t fPtV0 = TMath::Sqrt(v0->Pt2V0()); // transverse momentum of V0
4409 Double_t dSecVtxPos[3]; // V0 vertex position {x,y,z}
4410 v0->GetSecondaryVtx(dSecVtxPos);
4411 Double_t dDecayPath[3];
4412 for (Int_t iPos = 0; iPos < 3; iPos++)
4413 dDecayPath[iPos] = dSecVtxPos[iPos]-dPrimVtxPos[iPos]; // vector of the V0 path
4414 Float_t fDecLen2D = TMath::Sqrt(dDecayPath[0]*dDecayPath[0]+dDecayPath[1]*dDecayPath[1]); //transverse path length R
4415 Float_t fROverPt = fDecLen2D/fPtV0; // R/pT
4417 Float_t fMROverPtK0s = fMassK0s*fROverPt; // m*R/pT
4418 Float_t fMROverPtLambda = fMassLambda*fROverPt; // m*R/pT
4420 //___________________
4421 //Double_t fRap = -999;//init values
4422 Double_t fEta = -999;
4423 Double_t fV0cosPointAngle = -999;
4424 Double_t fV0DecayLength = v0->DecayLengthV0(lPrimaryVtxPosition);
4428 fV0mom[0]=v0->MomV0X();
4429 fV0mom[1]=v0->MomV0Y();
4430 fV0mom[2]=v0->MomV0Z();
4432 Double_t fV0TotalMomentum = TMath::Sqrt(fV0mom[0]*fV0mom[0]+fV0mom[1]*fV0mom[1]+fV0mom[2]*fV0mom[2]);
4433 // const Double_t K0sPDGmass = 0.497614;
4434 // const Double_t LambdaPDGmass = 1.115683;
4436 const Double_t K0sPDGmass = TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
4437 const Double_t LambdaPDGmass = TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
4439 Double_t fDistOverTotMomK0s = 0;
4440 Double_t fDistOverTotMomLa = 0;
4442 //calculate proper lifetime of particles in 3D (not recommended anymore)
4444 if(particletype == kK0){
4446 fDistOverTotMomK0s = fV0DecayLength * K0sPDGmass;
4447 fDistOverTotMomK0s /= (fV0TotalMomentum+1e-10);
4450 if((particletype == kLambda)||(particletype == kAntiLambda)){
4452 fDistOverTotMomLa = fV0DecayLength * LambdaPDGmass;
4453 fDistOverTotMomLa /= (fV0TotalMomentum+1e-10);
4456 //TPC cluster (not used anymore) and TPCRefit cuts
4458 //Double_t trackPosNcls = trackPos->GetTPCNcls();//Get number of clusters for positive charged tracks
4459 //Double_t trackNegNcls = trackNeg->GetTPCNcls();//Get number of clusters for negative charged tracks
4461 if(fRequireTPCRefit==kTRUE){//if kTRUE: accept only if daughter track is refitted in TPC!!
4462 Bool_t isPosTPCRefit = (trackPos->AliAODTrack::IsOn(AliESDtrack::kTPCrefit));
4463 Bool_t isNegTPCRefit = (trackNeg->AliAODTrack::IsOn(AliESDtrack::kTPCrefit));
4464 if (!isPosTPCRefit)continue;
4465 if (!isNegTPCRefit)continue;
4468 if(fKinkDaughters==kFALSE){//if kFALSE: no acceptance of kink daughters
4469 AliAODVertex* ProdVtxDaughtersPos = (AliAODVertex*) (trackPos->AliAODTrack::GetProdVertex());
4470 Char_t isAcceptKinkDaughtersPos = ProdVtxDaughtersPos->GetType();
4471 if(isAcceptKinkDaughtersPos==AliAODVertex::kKink)continue;
4473 AliAODVertex* ProdVtxDaughtersNeg = (AliAODVertex*) (trackNeg->AliAODTrack::GetProdVertex());
4474 Char_t isAcceptKinkDaughtersNeg = ProdVtxDaughtersNeg->GetType();
4475 if(isAcceptKinkDaughtersNeg==AliAODVertex::kKink)continue;
4479 Double_t fV0Radius = -999;
4480 Double_t fDcaV0Daughters = v0->DcaV0Daughters();
4481 Double_t fDcaPosToPrimVertex = v0->DcaPosToPrimVertex();//IP of positive charged daughter
4482 Double_t fDcaNegToPrimVertex = v0->DcaNegToPrimVertex();//IP of negative charged daughter
4483 Double_t avDecayLengthK0s = 2.6844;
4484 Double_t avDecayLengthLa = 7.89;
4486 //Float_t fCTauK0s = 2.6844; // [cm] c tau of K0S
4487 //Float_t fCTauLambda = 7.89; // [cm] c tau of Lambda and Antilambda
4489 fV0cosPointAngle = v0->CosPointingAngle(lPrimaryVtxPosition);
4490 lV0Position[0]= v0->DecayVertexV0X();
4491 lV0Position[1]= v0->DecayVertexV0Y();
4492 lV0Position[2]= v0->DecayVertexV0Z();
4494 fV0Radius = TMath::Sqrt(lV0Position[0]*lV0Position[0]+lV0Position[1]*lV0Position[1]);
4496 if(particletype == kK0) {//fRap = v0->RapK0Short();
4497 fEta = v0->PseudoRapV0();}
4498 if(particletype == kLambda) {//fRap = v0->RapLambda();
4499 fEta = v0->PseudoRapV0();}
4500 if(particletype == kAntiLambda) {//fRap = v0->Y(-3122);
4501 fEta = v0->PseudoRapV0();}
4504 //cut on 3D DistOverTotMom: (not used anymore)
4505 //if((particletype == kLambda)||(particletype == kAntiLambda)){if(fDistOverTotMomLa >= (fCutV0DecayMax * avDecayLengthLa)) continue;}
4507 //cut on K0s applied below after all other cuts for histo fill purposes..
4509 //cut on 2D DistOverTransMom: (recommended from Iouri)
4510 if((particletype == kLambda)||(particletype == kAntiLambda)){if(fMROverPtLambda > (fCutV0DecayMax * avDecayLengthLa))continue;}//fCutV0DecayMax set to 5 in AddTask macro
4512 //Armenteros Podolanski Plot for K0s:////////////////////////////
4514 Double_t ArmenterosAlpha=-999;
4515 Double_t ArmenterosPt=-999;
4521 if(particletype == kK0){
4523 pp[0]=v0->MomPosX();
4524 pp[1]=v0->MomPosY();
4525 pp[2]=v0->MomPosZ();
4527 pm[0]=v0->MomNegX();
4528 pm[1]=v0->MomNegY();
4529 pm[2]=v0->MomNegZ();
4532 v0mom[0]=v0->MomV0X();
4533 v0mom[1]=v0->MomV0Y();
4534 v0mom[2]=v0->MomV0Z();
4536 TVector3 v0Pos(pp[0],pp[1],pp[2]);
4537 TVector3 v0Neg(pm[0],pm[1],pm[2]);
4538 TVector3 v0totMom(v0mom[0], v0mom[1], v0mom[2]); //vector for tot v0 momentum
4540 //PosPt = v0Pos.Perp(v0totMom); //longitudinal momentum of positive charged daughter track
4541 PosPl = v0Pos.Dot(v0totMom)/v0totMom.Mag(); //transversal momentum of positive charged daughter track
4543 //NegPt = v0Neg.Perp(v0totMom); //longitudinal momentum of negative charged daughter track
4544 NegPl = v0Neg.Dot(v0totMom)/v0totMom.Mag(); //transversal momentum of nergative charged daughter track
4546 ArmenterosAlpha = 1.-2./(1+(PosPl/NegPl));
4547 ArmenterosPt= v0->PtArmV0();
4551 if(particletype == kK0){//only cut on K0s histos
4552 if(IsArmenterosSelected == 1){// Armenteros Cut to reject Lambdas contamination in K0s inv. massspectrum
4553 fh2ArmenterosBeforeCuts->Fill(ArmenterosAlpha,ArmenterosPt);
4557 //some more cuts on v0s and daughter tracks:
4560 if((TMath::Abs(PosEta)>fCutPostrackEta) || (TMath::Abs(NegEta)>fCutNegtrackEta))continue; //Daughters pseudorapidity cut
4561 if (fV0cosPointAngle < fCutV0cosPointAngle) continue; //cospointangle cut
4563 //if(TMath::Abs(fRap) > fCutRap)continue; //V0 Rapidity Cut
4564 if(TMath::Abs(fEta) > fCutEta) continue; //V0 Eta Cut
4565 if (fDcaV0Daughters > fCutDcaV0Daughters)continue;
4566 if ((fDcaPosToPrimVertex < fCutDcaPosToPrimVertex) || (fDcaNegToPrimVertex < fCutDcaNegToPrimVertex))continue;
4567 if ((fV0Radius < fCutV0RadiusMin) || (fV0Radius > fCutV0RadiusMax))continue;
4569 const AliAODPid *pid_p1=trackPos->GetDetPid();
4570 const AliAODPid *pid_n1=trackNeg->GetDetPid();
4573 if(particletype == kLambda){
4574 // if(AcceptBetheBloch(v0, fPIDResponse, 1) == kFALSE){std::cout<<"******PID cut rejects Lambda!!!************"<<std::endl;}
4575 if(AcceptBetheBloch(v0, fPIDResponse, 1) == kFALSE)continue;
4576 fh2BBLaPos->Fill(pid_p1->GetTPCmomentum(),pid_p1->GetTPCsignal());//positive lambda daughter
4577 fh2BBLaNeg->Fill(pid_n1->GetTPCmomentum(),pid_n1->GetTPCsignal());//negative lambda daughter
4579 //Double_t phi = v0->Phi();
4580 //Double_t massLa = v0->MassLambda();
4582 //printf("La: i = %d, m = %f, pT = %f, eta = %f, phi = %f\n, ",i,massLa,trackPt,fEta,phi);
4586 if(particletype == kAntiLambda){
4588 if(AcceptBetheBloch(v0, fPIDResponse, 2) == kFALSE)continue;
4589 fh2BBLaPos->Fill(pid_p1->GetTPCmomentum(),pid_p1->GetTPCsignal());//positive antilambda daughter
4590 fh2BBLaNeg->Fill(pid_n1->GetTPCmomentum(),pid_n1->GetTPCsignal());//negative antilambda daughter
4595 //Armenteros cut on K0s:
4596 if(particletype == kK0){
4597 if(IsArmenterosSelected == 1){// Armenteros Cut to reject Lambdas contamination in K0s inv. massspectrum
4599 if((ArmenterosPt<=(TMath::Abs(fCutArmenteros*ArmenterosAlpha))) && (fCutArmenteros!=-999))continue; //Cuts out Lambda contamination in K0s histos
4600 fh2ArmenterosAfterCuts->Fill(ArmenterosAlpha,ArmenterosPt);
4604 //not used anymore in 3D, z component of total momentum has bad resolution, cut instead in 2D and use pT
4605 //Proper Lifetime Cut: DecayLength3D * PDGmass / |p_tot| < 3*2.68cm (ctau(betagamma=1)) ; |p|/mass = beta*gamma
4606 //////////////////////////////////////////////
4609 //cut on 2D DistOverTransMom
4610 if(particletype == kK0){//the cut on Lambdas you can find above
4612 fh2ProperLifetimeK0sVsPtBeforeCut->Fill(trackPt,fMROverPtK0s); //fill these histos after all other cuts
4613 if(fMROverPtK0s > (fCutV0DecayMax * avDecayLengthK0s))continue;
4614 fh2ProperLifetimeK0sVsPtAfterCut->Fill(trackPt,fMROverPtK0s);
4616 //Double_t phi = v0->Phi();
4617 // Double_t massK0s = v0->MassK0Short();
4618 //printf("K0S: i = %d, m = %f, pT = %f, eta = %f, phi = %f\n",i,invMK0s,trackPt,fEta,phi);
4620 //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;
4623 //MC Associated V0 particles: (reconstructed particles associated with MC truth (MC truth: true primary MC generated particle))
4626 if(fAnalysisMC){// begin MC part
4628 Int_t negDaughterpdg = 0;
4629 Int_t posDaughterpdg = 0;
4631 Bool_t fPhysicalPrimary = -1; //v0 physical primary check
4632 Int_t MCv0PdgCode = 0;
4633 Bool_t mclabelcheck = kFALSE;
4635 TList *listmc = aod->GetList(); //AliAODEvent* is inherited from AliVEvent*, listmc is pointer to reconstructed event in MC list, member of AliAODEvent
4637 if(!listmc)continue;
4639 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
4641 //feeddown-correction for Lambda/Antilambda particles
4642 //feedddown comes mainly from charged and neutral Xi particles
4643 //feeddown from Sigma decays so quickly that it's not possible to distinguish from primary Lambdas with detector
4644 //feeddown for K0s from phi decays is neglectible
4645 //TH2F* fh2FeedDownMatrix = 0x0; //histo for feeddown already decleared above
4648 //first for all Lambda and Antilambda candidates____________________________________________________________________
4649 TString generatorName;
4651 if(particletype == kLambda){
4653 mclabelcheck = MCLabelCheck(v0, kLambda, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode, generatorName, isinjected);
4656 if((motherType == 3312)||(motherType == 3322)){//mother of v0 is neutral or negative Xi
4658 fListFeeddownLaCand->Add(v0); //fill TList with ass. particles, stemming from feeddown from Xi(bar) decays
4662 if(particletype == kAntiLambda){
4664 mclabelcheck = MCLabelCheck(v0, kAntiLambda, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode, generatorName, isinjected);
4666 if((motherType == -3312)||(motherType == -3322)){
4667 fListFeeddownALaCand->Add(v0); //fill TList with ass. particles, stemming from feeddown from Xi(bar) decays
4673 //_only true primary particles survive the following checks_______________________________________________________________________________________________
4674 TString generatorName;
4676 if(particletype == kK0){
4678 mclabelcheck = MCLabelCheck(v0, kK0, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode, generatorName, isinjected);
4679 if(mclabelcheck == kFALSE)continue;
4681 if(particletype == kLambda){
4683 mclabelcheck = MCLabelCheck(v0, kLambda, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode, generatorName, isinjected);
4684 if(mclabelcheck == kFALSE)continue;
4686 if(particletype == kAntiLambda){
4688 mclabelcheck = MCLabelCheck(v0, kAntiLambda, trackNeg, trackPos, listmc, negDaughterpdg, posDaughterpdg, motherType, v0Label, MCPt, fPhysicalPrimary, MCv0PdgCode, generatorName, isinjected);
4689 if(mclabelcheck == kFALSE)continue;
4692 if(fPhysicalPrimary != 1)continue; //V0 candidate (K0s, Lambda or Antilambda) must be physical primary, this means there is no mother particle existing
4701 Int_t nPart=list->GetSize();
4704 } // end GetListOfV0s()
4706 // -------------------------------------------------------------------------------------------------------
4708 void AliAnalysisTaskJetChem::CalculateInvMass(AliAODv0* v0vtx, const Int_t particletype, Double_t& invM, Double_t& trackPt){
4718 Double_t pp[3]={0,0,0}; //3-momentum positive charged track
4719 Double_t pm[3]={0,0,0}; //3-momentum negative charged track
4721 const Double_t massPi = 0.13957018; //better use PDG code at this point
4722 const Double_t massP = 0.93827203;
4727 TLorentzVector vector; //lorentzvector V0 particle
4728 TLorentzVector fourmom1;//lorentzvector positive daughter
4729 TLorentzVector fourmom2;//lorentzvector negative daughter
4731 //--------------------------------------------------------------
4733 AliAODTrack *trackPos = (AliAODTrack *) (v0vtx->GetSecondaryVtx()->GetDaughter(0));//index 0 defined as positive charged track in AliESDFilter
4735 if( trackPos->Charge() == 1 ){
4737 pp[0]=v0vtx->MomPosX();
4738 pp[1]=v0vtx->MomPosY();
4739 pp[2]=v0vtx->MomPosZ();
4741 pm[0]=v0vtx->MomNegX();
4742 pm[1]=v0vtx->MomNegY();
4743 pm[2]=v0vtx->MomNegZ();
4746 if( trackPos->Charge() == -1 ){
4748 pm[0]=v0vtx->MomPosX();
4749 pm[1]=v0vtx->MomPosY();
4750 pm[2]=v0vtx->MomPosZ();
4752 pp[0]=v0vtx->MomNegX();
4753 pp[1]=v0vtx->MomNegY();
4754 pp[2]=v0vtx->MomNegZ();
4757 if (particletype == kK0){ // case K0s
4758 mass1 = massPi;//positive particle
4759 mass2 = massPi;//negative particle
4760 } else if (particletype == kLambda){ // case Lambda
4761 mass1 = massP;//positive particle
4762 mass2 = massPi;//negative particle
4763 } else if (particletype == kAntiLambda){ //case AntiLambda
4764 mass1 = massPi;//positive particle
4765 mass2 = massP; //negative particle
4768 fourmom1.SetXYZM(pp[0],pp[1],pp[2],mass1);//positive track
4769 fourmom2.SetXYZM(pm[0],pm[1],pm[2],mass2);//negative track
4770 vector=fourmom1 + fourmom2;
4773 trackPt = vector.Pt();
4775 /*// 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
4777 if(particletype == kK0){
4778 std::cout << "invMK0s: " << invM <<std::endl;
4779 std::cout << "v0vtx->MassK0Short(): " << v0vtx->MassK0Short() << std::endl;
4780 std::cout << " " <<std::endl;
4781 //invM = v0vtx->MassK0Short();
4784 if(particletype == kLambda){
4785 std::cout << "invMLambda: " << invM <<std::endl;
4786 std::cout << "v0vtx->MassMassLambda(): " << v0vtx->MassLambda() << std::endl;
4787 std::cout << " " <<std::endl;
4788 //invM = v0vtx->MassLambda();
4791 if(particletype == kAntiLambda){
4792 std::cout << "invMAntiLambda: " << invM <<std::endl;
4793 std::cout << "v0vtx->MassAntiLambda(): " << v0vtx->MassAntiLambda() << std::endl;
4794 std::cout << " " <<std::endl;
4795 //invM = v0vtx->MassAntiLambda();
4803 //_____________________________________________________________________________________
4804 Int_t AliAnalysisTaskJetChem::GetListOfMCParticles(TList *outputlist, const Int_t particletype, AliAODEvent *mcaodevent) //(list to fill here e.g. fListMCgenK0s, particle species to search for)
4807 outputlist->Clear();
4809 TClonesArray *stack = 0x0;
4810 Double_t mcXv=0., mcYv=0., mcZv=0.;//MC vertex position
4813 // get MC generated particles
4815 Int_t fPdgcodeCurrentPart = 0; //pdg code current particle
4816 //Double_t fRapCurrentPart = 0; //get rapidity
4817 //Double_t fPtCurrentPart = 0; //get transverse momentum
4818 Double_t fEtaCurrentPart = 0; //get pseudorapidity
4820 //variable for check: physical primary particle
4821 //Bool_t IsPhysicalPrimary = -1;
4822 //Int_t index = 0; //check number of injected particles
4823 //****************************
4824 // Start loop over MC particles
4826 TList *lst = mcaodevent->GetList();
4829 if(fDebug>1) Printf("%s:%d no input list", (char*)__FILE__,__LINE__);
4833 stack = (TClonesArray*)lst->FindObject(AliAODMCParticle::StdBranchName());
4835 Printf("ERROR: stack not available");
4839 AliAODMCHeader *mcHdr=(AliAODMCHeader*)lst->FindObject(AliAODMCHeader::StdBranchName());
4840 if(!mcHdr)return -1;
4842 mcXv=mcHdr->GetVtxX(); mcYv=mcHdr->GetVtxY(); mcZv=mcHdr->GetVtxZ(); // position of the MC primary vertex
4845 ntrk=stack->GetEntriesFast();
4847 //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...
4850 for (Int_t iMc = 0; iMc < ntrk; iMc++) { //loop over mc generated particles
4853 AliAODMCParticle *p0=(AliAODMCParticle*)stack->UncheckedAt(iMc);
4855 //Printf("ERROR: particle with label %d not found in stack (mc loop)", iMc);
4858 fPdgcodeCurrentPart = p0->GetPdgCode();
4860 // Keep only K0s, Lambda and AntiLambda, Xi and Phi:
4861 //if ( (fPdgcodeCurrentPart != 310 ) && (fPdgcodeCurrentPart != 3122 ) && (fPdgcodeCurrentPart != -3122 ) && (fPdgcodeCurrentPart != 3312 ) && (fPdgcodeCurrentPart != -3312) && (fPdgcodeCurrentPart != -333) ) continue;
4865 //Rejection of Pythia injected particles with David Chinellatos method - not the latest method, better Method with TString from MC generator in IsInjected() function below!
4867 /* if( (p0->GetStatus()==21) ||
4868 ((p0->GetPdgCode() == 443) &&
4869 (p0->GetMother() == -1) &&
4870 (p0->GetDaughter(0) == (iMc))) ){ index++; }
4872 if(p0->GetStatus()==21){std::cout<< "hello !!!!" <<std::endl;}
4874 std::cout<< "MC particle status: " << p0->GetStatus() <<std::endl;
4878 //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
4881 //injected particles could be from GenBox (single high pt tracks) or jet related tracks, both generated from PYTHIA MC generator
4883 //Check: MC particle mother
4885 //for feed-down checks
4886 /* //MC gen particles
4887 Int_t iMother = p0->GetMother(); //Motherparticle of V0 candidate (e.g. phi particle,..)
4889 AliAODMCParticle *partM = (AliAODMCParticle*)stack->UncheckedAt(iMother);
4891 if(partM) codeM = TMath::Abs(partM->GetPdgCode());
4894 3312 Xi- -3312 Xibar+
4895 3322 Xi0 -3322 Xibar0
4898 if((codeM == 3312)||(codeM == 3322))// feeddown for Lambda coming from Xi- and Xi0
4904 /* //Check: MC gen. particle decays via 2-pion decay? -> only to be done for the rec. particles !! (-> branching ratio ~ 70 % for K0s -> pi+ pi-)
4906 Int_t daughter0Label = p0->GetDaughter(0);
4907 AliAODMCParticle *mcDaughter0 = (AliAODMCParticle *)stack->UncheckedAt(daughter0Label);
4908 if(daughter0Label >= 0)
4909 {daughter0Type = mcDaughter0->GetPdgCode();}
4911 Int_t daughter1Label = p0->GetDaughter(1);
4912 AliAODMCParticle *mcDaughter1 = (AliAODMCParticle *)stack->UncheckedAt(daughter1Label);
4914 if(daughter1Label >= 1)
4915 {daughter1Type = mcDaughter1->GetPdgCode();} //requirement that daughters are pions is only done for the reconstructed V0s in GetListofV0s() below
4920 // Keep only K0s, Lambda and AntiLambda:
4921 if ( (fPdgcodeCurrentPart != 310 ) && (fPdgcodeCurrentPart != 3122 ) && (fPdgcodeCurrentPart != -3122 )) continue;
4922 // Check: Is physical primary
4924 //do not use anymore: //IsPhysicalPrimary = p0->IsPhysicalPrimary();
4925 //if(!IsPhysicalPrimary)continue;
4927 Float_t fDistPrimaryMax = 0.01; // [cm] max distance of production point to the primary vertex (criterion for choice of MC particles considered as primary)
4929 // Get the distance between production point of the MC mother particle and the primary vertex
4931 Double_t dx = mcXv-p0->Xv();//mc primary vertex - mc gen. v0 vertex
4932 Double_t dy = mcYv-p0->Yv();
4933 Double_t dz = mcZv-p0->Zv();
4935 Double_t fDistPrimary = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
4936 Bool_t fPhysicalPrimary = (fDistPrimary < fDistPrimaryMax);
4938 if(!fPhysicalPrimary)continue;
4940 //if(fPhysicalPrimary){std::cout<<"hello**********************"<<std::endl;}
4942 /* std::cout<<"dx: "<<dx<<std::endl;
4943 std::cout<<"dy: "<<dy<<std::endl;
4944 std::cout<<"dz: "<<dz<<std::endl;
4946 std::cout<<"start: "<<std::endl;
4947 std::cout<<"mcXv: "<<mcXv<<std::endl;
4948 std::cout<<"mcYv: "<<mcYv<<std::endl;
4949 std::cout<<"mcZv: "<<mcZv<<std::endl;
4951 std::cout<<"p0->Xv(): "<<p0->Xv()<<std::endl;
4952 std::cout<<"p0->Yv(): "<<p0->Yv()<<std::endl;
4953 std::cout<<"p0->Zv(): "<<p0->Zv()<<std::endl;
4954 std::cout<<" "<<std::endl;
4955 std::cout<<"fDistPrimary"<<fDistPrimary<<std::endl;
4956 std::cout<<"fPhysicalPrimary"<<fPhysicalPrimary<<std::endl;
4959 //Is close enough to primary vertex to be considered as primary-like?
4961 //fRapCurrentPart = MyRapidity(p0->E(),p0->Pz());
4962 fEtaCurrentPart = p0->Eta();
4963 //fPtCurrentPart = p0->Pt();
4965 if (TMath::Abs(fEtaCurrentPart) < fCutEta){
4966 // if (TMath::Abs(fRapCurrentPart) > fCutRap)continue; //rap cut for crosschecks
4968 if(particletype == kK0){ //MC gen. K0s
4969 if (fPdgcodeCurrentPart==310){
4970 outputlist->Add(p0);
4974 if(particletype == kLambda){ //MC gen. Lambdas
4975 if (fPdgcodeCurrentPart==3122) {
4976 outputlist->Add(p0);
4980 if(particletype == kAntiLambda){
4981 if (fPdgcodeCurrentPart==-3122) { //MC gen. Antilambdas
4982 outputlist->Add(p0);
4987 }//end loop over MC generated particle
4989 Int_t nMCPart=outputlist->GetSize();
4996 //---------------------------------------------------------------------------
4998 Bool_t AliAnalysisTaskJetChem::FillFeeddownMatrix(TList* fListFeeddownCand, Int_t particletype)
5001 // Define Feeddown matrix
5002 Double_t lFeedDownMatrix [100][100];
5003 // FeedDownMatrix [Lambda Bin][Xi Bin];
5005 //Initialize entries of matrix:
5006 for(Int_t ilb = 0; ilb<100; ilb++){
5007 for(Int_t ixb = 0; ixb<100; ixb++){
5008 lFeedDownMatrix[ilb][ixb]=0; //first lambda bins, xi bins
5013 //----------------------------------------------------------------------------
5015 Double_t AliAnalysisTaskJetChem::MyRapidity(Double_t rE, Double_t rPz) const
5017 // Local calculation for rapidity
5018 return 0.5*TMath::Log((rE+rPz)/(rE-rPz+1.e-13));
5020 //----------------------------------------------------------------------------
5022 // ________________________________________________________________________________________________________________________//function to get the MC gen. jet particles
5024 void AliAnalysisTaskJetChem::GetTracksInCone(TList* inputlist, TList* outputlist, const AliAODJet* jet,
5025 const Double_t radius, Double_t& sumPt, const Double_t minPt, const Double_t maxPt, Bool_t& isBadPt)
5027 // fill list of tracks in cone around jet axis
5030 Bool_t isBadMaxPt = kFALSE;
5031 Bool_t isBadMinPt = kTRUE;
5035 jet->PxPyPz(jetMom);
5036 TVector3 jet3mom(jetMom);
5038 //if(jetets < jetetscutr)continue;
5040 for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){
5042 AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack));
5044 Double_t trackMom[3];
5045 track->PxPyPz(trackMom);
5046 TVector3 track3mom(trackMom);
5048 Double_t dR = jet3mom.DeltaR(track3mom);
5052 outputlist->Add(track);
5054 sumPt += track->Pt();
5056 if(maxPt>0 && track->Pt()>maxPt) isBadMaxPt = kTRUE; // reject jets containing any track with pt larger than this value, use GetFFMaxTrackPt()
5057 if(minPt>0 && track->Pt()>minPt) isBadMinPt = kFALSE; // reject jets with leading track with pt smaller than this value, use GetFFMinLTrackPt()
5063 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)
5064 if(maxPt>0 && isBadMaxPt) isBadPt = kTRUE; //..or because of leading track with too high pt (could be fake track)
5070 //____________________________________________________________________________________________________________________
5073 void AliAnalysisTaskJetChem::GetTracksInPerpCone(TList* inputlist, TList* outputlist, const AliAODJet* jet,
5074 const Double_t radius, Double_t& sumPerpPt)
5076 // fill list of tracks in two cones around jet axis rotated in phi +/- 90 degrees
5078 Double_t jetMom[3]; //array for entries in TVector3
5079 Double_t perpjetplusMom[3]; //array for entries in TVector3
5080 Double_t perpjetnegMom[3];
5084 jet->PxPyPz(jetMom); //get 3D jet momentum
5085 Double_t jetPerpPt = jet->Pt(); //original jet pt, invariant under rotations
5086 Double_t jetPhi = jet->Phi(); //original jet phi
5088 Double_t jetPerpposPhi = jetPhi + ((TMath::Pi())*0.5);//get new perp. jet axis phi clockwise
5089 Double_t jetPerpnegPhi = jetPhi - ((TMath::Pi())*0.5);//get new perp. jet axis phi counterclockwise
5091 TVector3 jet3mom(jetMom); //3-Vector for original rec. jet axis
5093 //Double_t phitest = jet3mom.Phi();
5095 perpjetplusMom[0]=(TMath::Cos(jetPerpposPhi)*jetPerpPt); //x coordinate (sidewards - when looking in beam direction)
5096 perpjetplusMom[1]=(TMath::Sin(jetPerpposPhi)*jetPerpPt); //y coordinate (upwards - when looking in beam direction)
5097 perpjetplusMom[2]=jetMom[2]; //z coordinate (along beam axis), invariant under azimuthal rotation
5099 perpjetnegMom[0]=(TMath::Cos(jetPerpnegPhi)*jetPerpPt); //x coordinate (sidewards - when looking in beam direction)
5100 perpjetnegMom[1]=(TMath::Sin(jetPerpnegPhi)*jetPerpPt); //y coordinate (upwards - when looking in beam direction)
5101 perpjetnegMom[2]=jetMom[2]; //z coordinate (along beam axis), invariant under azimuthal rotation
5104 TVector3 perpjetplus3mom(perpjetplusMom); //3-Vector for new perp. jet axis, clockwise rotated
5105 TVector3 perpjetneg3mom(perpjetnegMom); //3-Vector for new perp. jet axis, counterclockwise rotated
5107 //for crosscheck TVector3 rotation method
5109 //Double_t jetMomplusTest[3];
5110 //Double_t jetMomminusTest[3];
5112 //jet3mom.RotateZ(TMath::Pi()*0.5);//rotate original jet axis around +90 degrees in phi
5114 //perpjetminus3momTest = jet3mom.RotateZ((-1)*TMath::Pi()*0.5);
5116 // jet3mom.RotateZ(TMath::Pi()*0.5);
5117 // jet3mom.RotateZ((-1)*TMath::Pi()*0.5);
5119 //jetMomplusTest[0] = jet3mom.X(); //fetching perp. axis coordinates
5120 //jetMomplusTest[1] = jet3mom.Y();
5121 //jetMomplusTest[2] = jet3mom.Z();
5123 //TVector3 perpjetplus3momTest(jetMomplusTest); //new TVector3 for +90deg rotated jet axis with rotation method from ROOT
5124 //TVector3 perpjetminus3momTest(jetMomminusTest); //new TVector3 for -90deg rotated jet axis with rotation method from ROOT
5127 for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){ //collect V0 content in perp cone, rotated clockwise
5129 AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack)); //inputlist is fListK0s, all reconstructed K0s in event
5130 if(!track){std::cout<<"K0s track not found!!!"<<std::endl; continue;}
5132 Double_t trackMom[3];//3-mom of V0 particle
5133 track->PxPyPz(trackMom);
5134 TVector3 track3mom(trackMom);
5136 Double_t dR = perpjetplus3mom.DeltaR(track3mom);
5140 outputlist->Add(track); // output list is jetPerpConeK0list
5142 sumPerpPt += track->Pt();
5149 for (Int_t itrack=0; itrack<inputlist->GetSize(); itrack++){//collect V0 content in perp cone, rotated counterclockwise
5151 AliVParticle* track = dynamic_cast<AliVParticle*>(inputlist->At(itrack)); //inputlist is fListK0s, all reconstructed K0s in event
5152 if(!track){std::cout<<"K0s track not found!!!"<<std::endl; continue;}
5154 Double_t trackMom[3];//3-mom of V0 particle
5155 track->PxPyPz(trackMom);
5156 TVector3 track3mom(trackMom);
5158 Double_t dR = perpjetneg3mom.DeltaR(track3mom);
5162 outputlist->Add(track); // output list is jetPerpConeK0list
5164 sumPerpPt += track->Pt();
5170 // pay attention: this list contains the double amount of V0s, found in both cones
5171 // before using it, devide spectra by 2!!!
5172 sumPerpPt = sumPerpPt*0.5; //correct to do this?
5180 // _______________________________________________________________________________________________________________________________________________________
5182 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, TString& generatorName, Bool_t& isinjected){
5184 if(!v0)return kFALSE;
5186 TClonesArray *stackmc = 0x0;
5187 stackmc = (TClonesArray*)listmc->FindObject(AliAODMCParticle::StdBranchName()); //get MCAOD branch in data
5190 Printf("ERROR: stack not available");
5195 Int_t negAssLabel = TMath::Abs(trackNeg->GetLabel()); //negative (reconstructed) charged track label in MC stack
5196 Int_t posAssLabel = TMath::Abs(trackPos->GetLabel()); //positive (reconstructed) charged track label in MC stack
5198 //injected particle checks
5204 AliAODMCHeader *header=(AliAODMCHeader*)listmc->FindObject(AliAODMCHeader::StdBranchName());
5205 if(!header)return kFALSE;
5207 mcXv=header->GetVtxX(); mcYv=header->GetVtxY(); mcZv=header->GetVtxZ();
5213 if(negAssLabel>=0 && negAssLabel < stackmc->GetEntriesFast() && posAssLabel>=0 && posAssLabel < stackmc->GetEntriesFast()){//safety check if label has valid value of stack
5215 AliAODMCParticle *mcNegPart =(AliAODMCParticle*)stackmc->UncheckedAt(negAssLabel);//fetch the, with one MC truth track associated (reconstructed), negative charged track
5216 v0Label = mcNegPart->GetMother();// mother of negative charged particle is v0, get v0 label here
5217 negDaughterpdg = mcNegPart->GetPdgCode();
5218 AliAODMCParticle *mcPosPart =(AliAODMCParticle*)stackmc->UncheckedAt(posAssLabel);//fetch the, with one MC truth track associated (reconstructed), positive charged track
5219 Int_t v0PosLabel = mcPosPart->GetMother(); //get mother label of positive charged track label
5220 posDaughterpdg = mcPosPart->GetPdgCode();
5222 if(v0Label >= 0 && v0Label < stackmc->GetEntriesFast() && v0Label == v0PosLabel){//first v0 mc label check, then: check if both daughters are stemming from same particle
5224 AliAODMCParticle *mcv0 = (AliAODMCParticle *)stackmc->UncheckedAt(v0Label); //fetch MC ass. particle to v0 (mother of the both charged daughter tracks)
5226 //do not use anymore:
5227 //fPhysicalPrimary = mcv0->IsPhysicalPrimary();
5229 Float_t fDistPrimaryMax = 0.01; // [cm] max distance of production point to the primary vertex (criterion for choice of MC particles considered as primary)
5231 // Get the distance between production point of the MC mother particle and the primary vertex
5233 Double_t dx = mcXv-mcv0->Xv();//mc primary vertex - mc particle production vertex
5234 Double_t dy = mcYv-mcv0->Yv();
5235 Double_t dz = mcZv-mcv0->Zv();
5237 Float_t fDistPrimary = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
5238 fPhysicalPrimary = kFALSE;//init
5240 fPhysicalPrimary = (fDistPrimary < fDistPrimaryMax);
5241 MCv0PDGCode = mcv0->GetPdgCode();
5243 //if(fPhysicalPrimary == kTRUE){//look only at physical primary particles
5245 isinjected = IsTrackInjected(v0Label, header, stackmc, generatorName); //requires AliAODv0 instead of AliVTrack
5247 //trackinjected is kFALSE if it is either Hijing or has no generator name
5248 // std::cout<<" "<<std::endl;
5249 // std::cout<<"#### next particle: ####"<<std::endl;
5250 //std::cout<<"Is track injected: "<< trackinjected <<std::endl;
5251 // std::cout<<"pdg code: "<<MCv0PDGCode<<std::endl;
5252 // std::cout<<"v0Label: "<<v0Label<<std::endl;
5254 MCPt = mcv0->Pt();//for MC data, always use MC gen. pt for any pt distributions, also for the spectra, used for normalisation
5255 //for feed-down checks later
5257 Int_t motherLabel = mcv0->GetMother(); //get mother particle label of v0 particle
5258 // std::cout<<"motherLabel: "<<motherLabel<<std::endl;
5260 if(motherLabel >= 0 && v0Label < stackmc->GetEntriesFast()) //do safety check for mother label
5262 AliAODMCParticle *mcMother = (AliAODMCParticle *)stackmc->UncheckedAt(motherLabel); //get mother particle
5263 motherType = mcMother->GetPdgCode(); //get PDG code of mother
5266 Double_t XibarPt = 0.;
5268 if(particletype == kLambda){
5269 if((motherType == 3312)||(motherType == 3322)){ //if v0 mother is Xi0 or Xi- fill MC gen. pt in FD La histogram
5270 XiPt = mcMother->Pt();
5271 fh1MCXiPt->Fill(XiPt);
5274 if(particletype == kAntiLambda){
5275 if((motherType == -3312)||(motherType == -3322)){ //if v0 mother is Xibar0 or Xibar+ fill MC gen. pt in FD ALa histogram
5276 XibarPt = mcMother->Pt();
5277 fh1MCXibarPt->Fill(XibarPt);
5283 //pdg code checks etc..
5285 if(particletype == kK0){
5287 if(TMath::Abs(posDaughterpdg) != 211){return kFALSE;}//one or both of the daughters are not a pion
5288 if(TMath::Abs(negDaughterpdg) != 211){return kFALSE;}
5290 if(MCv0PDGCode != 310) {return kFALSE;}
5293 if(particletype == kLambda){
5294 if(MCv0PDGCode != 3122)return kFALSE;//if particle is not Antilambda, v0 is rejected
5295 if(posDaughterpdg != 2212)return kFALSE;
5296 if(negDaughterpdg != -211)return kFALSE; //pdg code check for Lambda daughters
5298 //{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 //}
5301 if(particletype == kAntiLambda){
5302 if(MCv0PDGCode != -3122)return kFALSE;
5303 if(posDaughterpdg != 211)return kFALSE;
5304 if(negDaughterpdg !=-2212)return kFALSE; //pdg code check for Antilambda daughters
5307 //{if((motherType == -3312)||(motherType == -3322)){continue;}//if bar{Xi0} and Xi+ is motherparticle of Antilambda, particle is rejected
5311 return kTRUE; //check was successful
5312 }//end mc v0 label check
5313 }// end of stack label check
5318 return kFALSE; //check wasn't successful
5320 //________________________________________________________________________________________________________________________________________________________
5323 Bool_t AliAnalysisTaskJetChem::IsParticleMatching(const AliAODMCParticle* mcp0, const Int_t v0Label){
5325 const Int_t mcp0label = mcp0->GetLabel();
5327 if(v0Label == mcp0label)return kTRUE;
5332 //_______________________________________________________________________________________________________________________________________________________
5334 Bool_t AliAnalysisTaskJetChem::DaughterTrackCheck(AliAODv0* v0, Int_t& nnum, Int_t& pnum){
5337 if(v0->GetNDaughters() != 2) return kFALSE;//case v0 has more or less than 2 daughters, avoids seg. break at some AOD files //reason?
5340 // safety check of input parameters
5343 if(fDebug > 1){std::cout << std::endl
5344 << "Warning in AliAnalysisTaskJetChem::DaughterTrackCheck:" << std::endl
5345 << "v0 = " << v0 << std::endl;}
5351 //Daughters track check: its Luke Hanrattys method to check daughters charge
5357 AliAODTrack *ntracktest =(AliAODTrack*)(v0->GetDaughter(nnum));
5359 if(ntracktest == NULL)
5361 if(fDebug > 1){std::cout << std::endl
5362 << "Warning in AliAnalysisTaskJetChem::DaughterTrackCheck:" << std::endl
5363 << "ntracktest = " << ntracktest << std::endl;}
5368 if(ntracktest->Charge() > 0)
5374 const AliAODTrack *trackNeg=(AliAODTrack *)(v0->GetDaughter(nnum));
5375 const AliAODTrack *trackPos=(AliAODTrack *)(v0->GetDaughter(pnum));
5377 //Check if both tracks are available
5378 if (!trackPos || !trackNeg) {
5379 if(fDebug > 1) Printf("strange analysis::UserExec:: Error:Could not retrieve one of the daughter tracks\n");
5384 //remove like sign V0s
5385 if ( trackPos->Charge() == trackNeg->Charge() ){
5386 //if(fDebug>1) Printf("%s:%d found like-sign V0", (char*)__FILE__,__LINE__);
5394 //______________________________________________________________________
5395 TString AliAnalysisTaskJetChem::GetGenerator(Int_t label, AliAODMCHeader* header){
5396 Int_t nsumpart=0;//number of particles
5397 TList *lh=header->GetCocktailHeaders();//TList with all generator headers
5398 Int_t nh=lh->GetEntries();//number of entries in TList with all headers
5400 for(Int_t i=0;i<nh;i++){
5401 AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(i);
5402 TString genname=gh->GetName();//name of particle generator
5403 Int_t npart=gh->NProduced();//number of stable or undecayed particles in MC stack block (?)
5404 if(label>=nsumpart && label<(nsumpart+npart)) return genname;
5411 //____________________________________________________________________________________________________________-
5413 Int_t AliAnalysisTaskJetChem::SplitCocktail(AliAODMCParticle *mcv0, Int_t v0Label, AliAODMCHeader *header, TClonesArray *arrayMC){//info in TString should be available from 2011 data productions on..
5415 if(!mcv0){std::cout << " !mcv0 " << std::endl;return 1;}
5416 if(!header){std::cout << " !header " << std::endl;return 1;}
5417 if(!arrayMC){std::cout << " !arrayMC " << std::endl;return 1;}
5419 //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
5420 //the complete amount of MC truth particles produced by several sources is named 'cocktail'
5422 //std::cout<<"v0 label: "<<v0Label<<std::endl;
5424 if(v0Label < 0) {std::cout<<"v0Label is negative!"<<std::endl;return 1;} //if label is negative, this particle is in the first part of the MC stack, it can be either Hijing or injected, but it is a primary particle
5426 TString generatorName = GetGenerator(v0Label,header);//this function returns a string with the generator name, used to produce this certain particle
5427 TString empty="";//no header was found
5429 std::cout << " TString generatorName: " << generatorName << std::endl;
5431 std::cout << " FIRST CALL " << generatorName << std::endl;
5433 while(generatorName.IsWhitespace()){
5434 AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(v0Label);
5435 if(!mcpart){return 1;}
5436 Int_t mother = mcpart->GetMother();
5438 generatorName = GetGenerator(mother,header);//see function directly below..
5439 std::cout << "Testing " << generatorName << " " << std::endl;
5442 std::cout << " FINAL CALL " << generatorName << std::endl;
5444 //std::transform(generatorName.begin(), generatorName.end(), generatorName.begin(), ::tolower); //convert TString bbb into lower case, to avoid that TString could be written in lower or upper case
5446 if(generatorName.Contains("ijing")){std::cout << " particle is Hijing!! " << 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"
5451 //_____________________________________________________________________
5452 void AliAnalysisTaskJetChem::GetTrackPrimaryGenerator(Int_t lab, AliAODMCHeader *header,TClonesArray *arrayMC,TString &nameGen){
5454 // method to check if a v0 comes from a given generator
5456 nameGen=GetGenerator(lab,header);
5458 // Int_t countControl=0;
5460 while(nameGen.IsWhitespace()){
5461 AliAODMCParticle *mcpart= (AliAODMCParticle*)arrayMC->At(lab);//get MC generated particle for AliAODv0 particle
5463 printf("AliAnalysisTaskJetChem::IsTrackInjected - BREAK: No valid AliAODMCParticle at label %i\n",lab);
5466 Int_t mother = mcpart->GetMother();
5469 printf("AliAnalysisTaskJetChem::IsTrackInjected - BREAK: Reached primary particle without valid mother\n");
5473 nameGen=GetGenerator(mother,header);
5476 // if(countControl>=10){ // 10 = arbitrary number; protection from infinite loops
5477 // printf("AliVertexingHFUtils::IsTrackInjected - BREAK: Protection from infinite loop active\n");
5487 //---------------------------------------------------------------------------------------------------------------------
5488 Bool_t AliAnalysisTaskJetChem::IsTrackInjected(Int_t lab, AliAODMCHeader *header,TClonesArray *arrayMC, TString& nameGen){
5489 // method to check if a v0 particle comes from the signal event or from the underlying Hijing event
5492 GetTrackPrimaryGenerator(lab, header, arrayMC, nameGen);
5494 if(nameGen.IsWhitespace() || nameGen.Contains("ijing")) return kFALSE;//particle has either no info about generator or is Hijing particle, so it is not injected
5496 //std::cout<<"generator name: "<<nameGen<<std::endl;
5501 //_________________________________________________________________________________________________________________________________________
5502 Double_t AliAnalysisTaskJetChem::SmearJetPt(Double_t jetPt, Int_t /*cent*/, Double_t /*jetRadius*/, Double_t /*ptmintrack*/, Double_t& jetPtSmear){
5504 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
5508 /* if(cent>10) cl = 2;
5513 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
5514 //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
5516 //fsmear->SetParameters(1,0,4.472208);// for 2010 PbPb jets, R=0.2, ptmintrack = 0.15 GeV/c, cent 00-10%
5518 /* //delta-pt width for anti-kt jet finder:
5521 if((cl == 1)&&(jetRadius == 0.4)&&(ptmintrack == 0.15)){
5522 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
5524 if((cl == 2)&&(jetRadius == 0.4)&&(ptmintrack == 0.15)){
5525 fsmear->SetParameters(1,0,8.536195);
5527 if((cl == 3)&&(jetRadius == 0.4)&&(ptmintrack == 0.15)){
5528 fsmear->SetParameters(1,0,?);
5530 if((cl == 4)&&(jetRadius == 0.4)&&(ptmintrack == 0.15)){
5531 fsmear->SetParameters(1,0,5.229839);
5535 if((cl == 1)&&(jetRadius == 0.3)&&(ptmintrack == 0.15)){
5536 fsmear->SetParameters(1,0,7.145967);
5538 if((cl == 2)&&(jetRadius == 0.3)&&(ptmintrack == 0.15)){
5539 fsmear->SetParameters(1,0,5.844796);
5541 if((cl == 3)&&(jetRadius == 0.3)&&(ptmintrack == 0.15)){
5542 fsmear->SetParameters(1,0,?);
5544 if((cl == 4)&&(jetRadius == 0.3)&&(ptmintrack == 0.15)){
5545 fsmear->SetParameters(1,0,3.630751);
5549 if((cl == 1)&&(jetRadius == 0.2)&&(ptmintrack == 0.15)){
5550 fsmear->SetParameters(1,0,4.472208);
5552 if((cl == 2)&&(jetRadius == 0.2)&&(ptmintrack == 0.15)){
5553 fsmear->SetParameters(1,0,3.543938);
5555 if((cl == 3)&&(jetRadius == 0.2)&&(ptmintrack == 0.15)){
5556 fsmear->SetParameters(1,0,?);
5558 if((cl == 4)&&(jetRadius == 0.2)&&(ptmintrack == 0.15)){
5559 fsmear->SetParameters(1,0,1.037476);
5564 Double_t r = fsmear.GetRandom();
5565 jetPtSmear = jetPt + r;
5567 // std::cout<<"jetPt: "<<jetPt<<std::endl;
5568 // std::cout<<"jetPtSmear: "<<jetPtSmear<<std::endl;
5569 // std::cout<<"r: "<<r<<std::endl;
5576 //______________________________________________________________________________________________________________________
5577 //____________________________________________________________________________________________________________________
5579 Bool_t AliAnalysisTaskJetChem::IsParticleInCone(const AliVParticle* part1, const AliVParticle* part2, Double_t dRMax) const
5581 // decides whether a particle is inside a jet cone
5582 if (!part1 || !part2)
5585 TVector3 vecMom2(part2->Px(),part2->Py(),part2->Pz());
5586 TVector3 vecMom1(part1->Px(),part1->Py(),part1->Pz());
5587 Double_t dR = vecMom2.DeltaR(vecMom1); // = sqrt(dEta*dEta+dPhi*dPhi)
5588 if(dR<dRMax) // momentum vectors of part1 and part2 are closer than dRMax
5592 //__________________________________________________________________________________________________________________
5595 Bool_t AliAnalysisTaskJetChem::IsRCJCOverlap(TList* recjetlist, const AliVParticle* part, Double_t dDistance) const{
5597 if(!recjetlist) return kFALSE;
5598 if(!part) return kFALSE;
5599 if(!dDistance) return kFALSE;
5600 Int_t nRecJetsCuts = fJetsRecCuts->GetEntries();
5602 for(Int_t i=0; i<nRecJetsCuts; ++i){ //loop over all reconstructed jets in events
5603 AliAODJet* jet = (AliAODJet*) (recjetlist->At(i));
5604 if(!jet){if(fDebug>2)std::cout<<"AliAnalysisTaskJetChem::IsRCJCOverlap jet pointer invalid!"<<std::endl;continue;}
5605 if(IsParticleInCone(jet, part, dDistance) == kTRUE)return kTRUE;//RC and JC are overlapping
5607 }//end loop testing RC-JC overlap
5608 return kFALSE;//RC and JC are not overlapping -> good!
5611 //_______________________________________________________________________________________________________________________
5612 AliAODJet* AliAnalysisTaskJetChem::GetRandomCone(TList* jetlist, Double_t dEtaConeMax, Double_t dDistance) const
5614 TLorentzVector vecRdCone;
5615 AliAODJet* jetRC = 0;//random cone candidate
5616 Double_t dEta, dPhi; //random eta and phi value for RC
5617 Bool_t IsRCoutJC = kFALSE;//check whether RC is not overlapping with any selected jet cone in event
5618 Int_t iRCTrials = 10;//search at maximum 10 times for random cone that doesn't overlap with jet cone
5620 for(Int_t i=0; i<iRCTrials; iRCTrials++){
5622 dEta = dEtaConeMax*(2*fRandom->Rndm()-1.); //random eta value in range: [-dEtaConeMax,+dEtaConeMax]
5623 dPhi = TMath::TwoPi()*fRandom->Rndm(); //random phi value in range: [0,2*Pi]
5624 vecRdCone.SetPtEtaPhiM(1.,dEta,dPhi,0.);
5625 jetRC = new AliAODJet(vecRdCone);//new RC candidate
5627 if (!IsRCJCOverlap(jetlist,jetRC,dDistance))
5629 IsRCoutJC = kTRUE; //std::cout<<"RC and JC are not overlapping!!!"<<std::endl;
5633 delete jetRC; //RC is overlapping with JC, delete this RC candidate
5636 if(!IsRCoutJC) {jetRC = 0;}//in case no random cone was selected
5642 // _______________________________________________________________________________________________________________________
5643 AliAODJet* AliAnalysisTaskJetChem::GetMedianCluster()
5645 // fill tracks from bckgCluster branch,
5646 // using cluster with median density (odd number of clusters)
5647 // or picking randomly one of the two closest to median (even number)
5649 Int_t nBckgClusters = fBckgJetsRec->GetEntries(); // not 'recCuts': use all clusters in full eta range
5651 if(nBckgClusters<3) return 0; // need at least 3 clusters (skipping 2 highest)
5653 Double_t* bgrDensity = new Double_t[nBckgClusters];
5654 Int_t* indices = new Int_t[nBckgClusters];
5656 for(Int_t ij=0; ij<nBckgClusters; ++ij){
5658 AliAODJet* bgrCluster = (AliAODJet*)(fBckgJetsRec->At(ij));
5659 Double_t clusterPt = bgrCluster->Pt();
5660 Double_t area = bgrCluster->EffectiveAreaCharged();
5662 Double_t density = 0;
5663 if(area>0) density = clusterPt/area;
5665 bgrDensity[ij] = density;
5669 TMath::Sort(nBckgClusters, bgrDensity, indices);
5671 // get median cluster
5673 AliAODJet* medianCluster = 0;
5675 if(TMath::Odd(nBckgClusters)){
5677 Int_t medianIndex = indices[(Int_t) (0.5*(nBckgClusters+1))];
5679 medianCluster = (AliAODJet*)(fBckgJetsRec->At(medianIndex));
5681 Double_t clusterPt = medianCluster->Pt();
5682 Double_t area = medianCluster->EffectiveAreaCharged();
5686 Int_t medianIndex1 = indices[(Int_t) (0.5*nBckgClusters)];
5687 Int_t medianIndex2 = indices[(Int_t) (0.5*nBckgClusters+1)];
5689 AliAODJet* medianCluster1 = (AliAODJet*)(fBckgJetsRec->At(medianIndex1));
5690 AliAODJet* medianCluster2 = (AliAODJet*)(fBckgJetsRec->At(medianIndex2));
5692 Double_t density1 = 0;
5693 Double_t clusterPt1 = medianCluster1->Pt();
5694 Double_t area1 = medianCluster1->EffectiveAreaCharged();
5695 if(area1>0) density1 = clusterPt1/area1;
5697 Double_t density2 = 0;
5698 Double_t clusterPt2 = medianCluster2->Pt();
5699 Double_t area2 = medianCluster2->EffectiveAreaCharged();
5700 if(area2>0) density2 = clusterPt2/area2;
5702 medianCluster = ( (gRandom->Rndm()>0.5) ? medianCluster1 : medianCluster2 ); // select one randomly to avoid adding areas
5705 delete[] bgrDensity;
5708 return medianCluster;
5711 //____________________________________________________________________________________________
5713 Double_t AliAnalysisTaskJetChem::AreaCircSegment(Double_t dRadius, Double_t dDistance) const
5715 // calculate area of a circular segment defined by the circle radius and the (oriented) distance between the secant line and the circle centre
5716 Double_t dEpsilon = 1e-2;
5717 Double_t dR = dRadius;
5718 Double_t dD = dDistance;
5719 if (TMath::Abs(dR)<dEpsilon)
5721 if(fDebug>0) printf("AliAnalysisTaskJetChem::AreaCircSegment: Error: Too small radius: %f < %f\n",dR,dEpsilon);
5727 return TMath::Pi()*dR*dR;
5728 return dR*dR*TMath::ACos(dD/dR)-dD*TMath::Sqrt(dR*dR-dD*dD);