2424#include < array>
2525#include < cmath>
2626#include < cstdlib>
27+ #include < vector>
2728#include " TF1.h"
2829#include " TRandom3.h"
2930#include " Math/Vector3D.h"
@@ -79,6 +80,8 @@ struct strangeness_tutorial {
7980 Configurable<bool > goodzvertex{" goodzvertex" , false , " removes collisions with large differences between z of PV by tracks and z of PV from FT0 A-C time difference." };
8081 Configurable<bool > itstpctracks{" itstpctracks" , false , " selects collisions with at least one ITS-TPC track," };
8182 Configurable<bool > additionalEvsel{" additionalEvsel" , false , " Additional event selcection" };
83+ Configurable<bool > applyOccupancyCut{" applyOccupancyCut" , false , " Apply occupancy cut" };
84+ Configurable<int > OccupancyCut{" OccupancyCut" , 1000 , " Mimimum Occupancy cut" };
8285
8386 // Configurable parameters for V0 selection
8487 Configurable<float > ConfV0DCADaughMax{" ConfV0DCADaughMax" , 1 .0f , " DCA b/w V0 daughters" };
@@ -140,6 +143,7 @@ struct strangeness_tutorial {
140143 ConfigurableAxis axisdEdx{" axisdEdx" , {20000 , 0 .0f , 200 .0f }, " dE/dx (a.u.)" };
141144 ConfigurableAxis axisPtfordEbydx{" axisPtfordEbydx" , {2000 , 0 , 20 }, " pT (GeV/c)" };
142145 ConfigurableAxis axisMultdist{" axisMultdist" , {3500 , 0 , 70000 }, " Multiplicity distribution" };
146+ ConfigurableAxis occupancy_bins{" occupancy_bins" , {VARIABLE_WIDTH, 0.0 , 100 , 500 , 600 , 1000 , 1100 , 1500 , 1600 , 2000 , 2100 , 2500 , 2600 , 3000 , 3100 , 3500 , 3600 , 4000 , 4100 , 4500 , 4600 , 5000 , 5100 , 9999 }, " Binning of the occupancy axis" };
143147
144148 // Event selection cuts - Alex (Temporary, need to fix!)
145149 TF1* fMultPVCutLow = nullptr ;
@@ -160,6 +164,7 @@ struct strangeness_tutorial {
160164 // AxisSpec multiplicityAxis = {110, 0.0f, 150.0f, "Multiplicity Axis"};
161165 AxisSpec multiplicityAxis = {binsCent, " Multiplicity Axis" };
162166 AxisSpec thnAxisPOL{configThnAxisPOL, " Configurabel theta axis" };
167+ AxisSpec occupancy_axis = {occupancy_bins, " Occupancy [-40,100]" };
163168
164169 // THnSparses
165170 std::array<int , 4 > sparses = {activateTHnSparseCosThStarHelicity, activateTHnSparseCosThStarProduction, activateTHnSparseCosThStarBeam, activateTHnSparseCosThStarRandom};
@@ -199,9 +204,9 @@ struct strangeness_tutorial {
199204 hglue.add (" h1glueInvMassRot" , " h1glueInvMassRot" , kTH1F , {glueballMassAxis});
200205 }
201206
202- hglue.add (" h3glueInvMassDS" , " h3glueInvMassDS" , kTHnSparseF , {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL}, true );
203- hglue.add (" h3glueInvMassME" , " h3glueInvMassME" , kTHnSparseF , {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL}, true );
204- hglue.add (" h3glueInvMassRot" , " h3glueInvMassRot" , kTHnSparseF , {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL}, true );
207+ hglue.add (" h3glueInvMassDS" , " h3glueInvMassDS" , kTHnSparseF , {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL, occupancy_axis }, true );
208+ hglue.add (" h3glueInvMassME" , " h3glueInvMassME" , kTHnSparseF , {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL, occupancy_axis }, true );
209+ hglue.add (" h3glueInvMassRot" , " h3glueInvMassRot" , kTHnSparseF , {multiplicityAxis, ptAxis, glueballMassAxis, thnAxisPOL, occupancy_axis }, true );
205210 hglue.add (" heventscheck" , " heventscheck" , kTH1I , {{10 , 0 , 10 }});
206211 hglue.add (" htrackscheck_v0" , " htrackscheck_v0" , kTH1I , {{15 , 0 , 15 }});
207212 hglue.add (" htrackscheck_v0_daughters" , " htrackscheck_v0_daughters" , kTH1I , {{15 , 0 , 15 }});
@@ -540,6 +545,11 @@ struct strangeness_tutorial {
540545 return ;
541546 }
542547
548+ auto occupancy_no = collision.trackOccupancyInTimeRange ();
549+ if (applyOccupancyCut && occupancy_no < OccupancyCut) {
550+ return ;
551+ }
552+
543553 if (QAevents) {
544554 rEventSelection.fill (HIST (" hVertexZRec" ), collision.posZ ());
545555 rEventSelection.fill (HIST (" hmultiplicity" ), multiplicity);
@@ -644,45 +654,45 @@ struct strangeness_tutorial {
644654 if (activateTHnSparseCosThStarHelicity) {
645655 ROOT::Math::XYZVector helicityVec = fourVecMother.Vect (); // 3 vector of mother in COM frame
646656 auto cosThetaStarHelicity = helicityVec.Dot (threeVecDauCM) / (std::sqrt (threeVecDauCM.Mag2 ()) * std::sqrt (helicityVec.Mag2 ()));
647- hglue.fill (HIST (" h3glueInvMassDS" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarHelicity);
657+ hglue.fill (HIST (" h3glueInvMassDS" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarHelicity, occupancy_no );
648658 for (int i = 0 ; i < c_nof_rotations; i++) {
649659 float theta2 = rn->Uniform (TMath::Pi () - TMath::Pi () / rotational_cut, TMath::Pi () + TMath::Pi () / rotational_cut);
650660 lv4.SetPtEtaPhiM (v1.pt (), v1.eta (), v1.phi () + theta2, massK0s); // for rotated background
651661 lv5 = lv2 + lv4;
652- hglue.fill (HIST (" h3glueInvMassRot" ), multiplicity, lv5.Pt (), lv5.M (), cosThetaStarHelicity);
662+ hglue.fill (HIST (" h3glueInvMassRot" ), multiplicity, lv5.Pt (), lv5.M (), cosThetaStarHelicity, occupancy_no );
653663 }
654664
655665 } else if (activateTHnSparseCosThStarProduction) {
656666 ROOT::Math::XYZVector normalVec = ROOT::Math::XYZVector (lv3.Py (), -lv3.Px (), 0 .f );
657667 auto cosThetaStarProduction = normalVec.Dot (threeVecDauCM) / (std::sqrt (threeVecDauCM.Mag2 ()) * std::sqrt (normalVec.Mag2 ()));
658- hglue.fill (HIST (" h3glueInvMassDS" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarProduction);
668+ hglue.fill (HIST (" h3glueInvMassDS" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarProduction, occupancy_no );
659669 for (int i = 0 ; i < c_nof_rotations; i++) {
660670 float theta2 = rn->Uniform (TMath::Pi () - TMath::Pi () / rotational_cut, TMath::Pi () + TMath::Pi () / rotational_cut);
661671 lv4.SetPtEtaPhiM (v1.pt (), v1.eta (), v1.phi () + theta2, massK0s); // for rotated background
662672 lv5 = lv2 + lv4;
663- hglue.fill (HIST (" h3glueInvMassRot" ), multiplicity, lv5.Pt (), lv5.M (), cosThetaStarProduction);
673+ hglue.fill (HIST (" h3glueInvMassRot" ), multiplicity, lv5.Pt (), lv5.M (), cosThetaStarProduction, occupancy_no );
664674 }
665675 } else if (activateTHnSparseCosThStarBeam) {
666676 ROOT::Math::XYZVector beamVec = ROOT::Math::XYZVector (0 .f , 0 .f , 1 .f );
667677 auto cosThetaStarBeam = beamVec.Dot (threeVecDauCM) / std::sqrt (threeVecDauCM.Mag2 ());
668- hglue.fill (HIST (" h3glueInvMassDS" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarBeam);
678+ hglue.fill (HIST (" h3glueInvMassDS" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarBeam, occupancy_no );
669679 for (int i = 0 ; i < c_nof_rotations; i++) {
670680 float theta2 = rn->Uniform (TMath::Pi () - TMath::Pi () / rotational_cut, TMath::Pi () + TMath::Pi () / rotational_cut);
671681 lv4.SetPtEtaPhiM (v1.pt (), v1.eta (), v1.phi () + theta2, massK0s); // for rotated background
672682 lv5 = lv2 + lv4;
673- hglue.fill (HIST (" h3glueInvMassRot" ), multiplicity, lv5.Pt (), lv5.M (), cosThetaStarBeam);
683+ hglue.fill (HIST (" h3glueInvMassRot" ), multiplicity, lv5.Pt (), lv5.M (), cosThetaStarBeam, occupancy_no );
674684 }
675685 } else if (activateTHnSparseCosThStarRandom) {
676686 auto phiRandom = gRandom ->Uniform (0 .f , constants::math::TwoPI);
677687 auto thetaRandom = gRandom ->Uniform (0 .f , constants::math::PI);
678688 ROOT::Math::XYZVector randomVec = ROOT::Math::XYZVector (std::sin (thetaRandom) * std::cos (phiRandom), std::sin (thetaRandom) * std::sin (phiRandom), std::cos (thetaRandom));
679689 auto cosThetaStarRandom = randomVec.Dot (threeVecDauCM) / std::sqrt (threeVecDauCM.Mag2 ());
680- hglue.fill (HIST (" h3glueInvMassDS" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarRandom);
690+ hglue.fill (HIST (" h3glueInvMassDS" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarRandom, occupancy_no );
681691 for (int i = 0 ; i < c_nof_rotations; i++) {
682692 float theta2 = rn->Uniform (TMath::Pi () - TMath::Pi () / rotational_cut, TMath::Pi () + TMath::Pi () / rotational_cut);
683693 lv4.SetPtEtaPhiM (v1.pt (), v1.eta (), v1.phi () + theta2, massK0s); // for rotated background
684694 lv5 = lv2 + lv4;
685- hglue.fill (HIST (" h3glueInvMassRot" ), multiplicity, lv5.Pt (), lv5.M (), cosThetaStarRandom);
695+ hglue.fill (HIST (" h3glueInvMassRot" ), multiplicity, lv5.Pt (), lv5.M (), cosThetaStarRandom, occupancy_no );
686696 }
687697 }
688698 }
@@ -725,6 +735,11 @@ struct strangeness_tutorial {
725735 if (!eventselection (c1, multiplicity) || !eventselection (c2, multiplicity)) {
726736 continue ;
727737 }
738+ auto occupancy_no = c1.trackOccupancyInTimeRange ();
739+ auto occupancy_no2 = c2.trackOccupancyInTimeRange ();
740+ if (applyOccupancyCut && (occupancy_no < OccupancyCut || occupancy_no2 < OccupancyCut)) {
741+ return ;
742+ }
728743
729744 for (auto & [t1, t2] : o2::soa::combinations (o2::soa::CombinationsFullIndexPolicy (tracks1, tracks2))) {
730745
@@ -782,21 +797,21 @@ struct strangeness_tutorial {
782797 if (activateTHnSparseCosThStarHelicity) {
783798 ROOT::Math::XYZVector helicityVec = fourVecMother.Vect (); // 3 vector of mother in COM frame
784799 auto cosThetaStarHelicity = helicityVec.Dot (threeVecDauCM) / (std::sqrt (threeVecDauCM.Mag2 ()) * std::sqrt (helicityVec.Mag2 ()));
785- hglue.fill (HIST (" h3glueInvMassME" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarHelicity);
800+ hglue.fill (HIST (" h3glueInvMassME" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarHelicity, occupancy_no );
786801 } else if (activateTHnSparseCosThStarProduction) {
787802 ROOT::Math::XYZVector normalVec = ROOT::Math::XYZVector (lv3.Py (), -lv3.Px (), 0 .f );
788803 auto cosThetaStarProduction = normalVec.Dot (threeVecDauCM) / (std::sqrt (threeVecDauCM.Mag2 ()) * std::sqrt (normalVec.Mag2 ()));
789- hglue.fill (HIST (" h3glueInvMassME" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarProduction);
804+ hglue.fill (HIST (" h3glueInvMassME" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarProduction, occupancy_no );
790805 } else if (activateTHnSparseCosThStarBeam) {
791806 ROOT::Math::XYZVector beamVec = ROOT::Math::XYZVector (0 .f , 0 .f , 1 .f );
792807 auto cosThetaStarBeam = beamVec.Dot (threeVecDauCM) / std::sqrt (threeVecDauCM.Mag2 ());
793- hglue.fill (HIST (" h3glueInvMassME" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarBeam);
808+ hglue.fill (HIST (" h3glueInvMassME" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarBeam, occupancy_no );
794809 } else if (activateTHnSparseCosThStarRandom) {
795810 auto phiRandom = gRandom ->Uniform (0 .f , constants::math::TwoPI);
796811 auto thetaRandom = gRandom ->Uniform (0 .f , constants::math::PI);
797812 ROOT::Math::XYZVector randomVec = ROOT::Math::XYZVector (std::sin (thetaRandom) * std::cos (phiRandom), std::sin (thetaRandom) * std::sin (phiRandom), std::cos (thetaRandom));
798813 auto cosThetaStarRandom = randomVec.Dot (threeVecDauCM) / std::sqrt (threeVecDauCM.Mag2 ());
799- hglue.fill (HIST (" h3glueInvMassME" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarRandom);
814+ hglue.fill (HIST (" h3glueInvMassME" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarRandom, occupancy_no );
800815 }
801816 }
802817
@@ -817,6 +832,11 @@ struct strangeness_tutorial {
817832 if (!eventselection (c1, multiplicity) || !eventselection (c2, multiplicity)) {
818833 continue ;
819834 }
835+ auto occupancy_no = c1.trackOccupancyInTimeRange ();
836+ auto occupancy_no2 = c2.trackOccupancyInTimeRange ();
837+ if (applyOccupancyCut && (occupancy_no < OccupancyCut || occupancy_no2 < OccupancyCut)) {
838+ return ;
839+ }
820840
821841 for (auto & [t1, t2] : o2::soa::combinations (o2::soa::CombinationsFullIndexPolicy (tracks1, tracks2))) {
822842 if (t1.size () == 0 || t2.size () == 0 ) {
@@ -873,21 +893,21 @@ struct strangeness_tutorial {
873893 if (activateTHnSparseCosThStarHelicity) {
874894 ROOT::Math::XYZVector helicityVec = fourVecMother.Vect (); // 3 vector of mother in COM frame
875895 auto cosThetaStarHelicity = helicityVec.Dot (threeVecDauCM) / (std::sqrt (threeVecDauCM.Mag2 ()) * std::sqrt (helicityVec.Mag2 ()));
876- hglue.fill (HIST (" h3glueInvMassME" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarHelicity);
896+ hglue.fill (HIST (" h3glueInvMassME" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarHelicity, occupancy_no );
877897 } else if (activateTHnSparseCosThStarProduction) {
878898 ROOT::Math::XYZVector normalVec = ROOT::Math::XYZVector (lv3.Py (), -lv3.Px (), 0 .f );
879899 auto cosThetaStarProduction = normalVec.Dot (threeVecDauCM) / (std::sqrt (threeVecDauCM.Mag2 ()) * std::sqrt (normalVec.Mag2 ()));
880- hglue.fill (HIST (" h3glueInvMassME" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarProduction);
900+ hglue.fill (HIST (" h3glueInvMassME" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarProduction, occupancy_no );
881901 } else if (activateTHnSparseCosThStarBeam) {
882902 ROOT::Math::XYZVector beamVec = ROOT::Math::XYZVector (0 .f , 0 .f , 1 .f );
883903 auto cosThetaStarBeam = beamVec.Dot (threeVecDauCM) / std::sqrt (threeVecDauCM.Mag2 ());
884- hglue.fill (HIST (" h3glueInvMassME" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarBeam);
904+ hglue.fill (HIST (" h3glueInvMassME" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarBeam, occupancy_no );
885905 } else if (activateTHnSparseCosThStarRandom) {
886906 auto phiRandom = gRandom ->Uniform (0 .f , constants::math::TwoPI);
887907 auto thetaRandom = gRandom ->Uniform (0 .f , constants::math::PI);
888908 ROOT::Math::XYZVector randomVec = ROOT::Math::XYZVector (std::sin (thetaRandom) * std::cos (phiRandom), std::sin (thetaRandom) * std::sin (phiRandom), std::cos (thetaRandom));
889909 auto cosThetaStarRandom = randomVec.Dot (threeVecDauCM) / std::sqrt (threeVecDauCM.Mag2 ());
890- hglue.fill (HIST (" h3glueInvMassME" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarRandom);
910+ hglue.fill (HIST (" h3glueInvMassME" ), multiplicity, lv3.Pt (), lv3.M (), cosThetaStarRandom, occupancy_no );
891911 }
892912 }
893913
0 commit comments