5555#include " DataFormatsParameters/GRPMagField.h"
5656#include " CCDB/BasicCCDBManager.h"
5757#include " Common/DataModel/PIDResponseITS.h"
58+ #include " PWGMM/Mult/DataModel/Index.h" // for Particles2Tracks table
5859
5960using namespace o2 ;
6061using namespace o2 ::framework;
@@ -69,12 +70,14 @@ struct phipbpb {
6970 Service<o2::framework::O2DatabasePDG> pdg;
7071
7172 // CCDB options
72- Configurable<std::string> ccdburl{" ccdb-url" , " http://alice-ccdb.cern.ch" , " url of the ccdb repository" };
73- Configurable<std::string> grpPath{" grpPath" , " GLO/GRP/GRP" , " Path of the grp file" };
74- Configurable<std::string> grpmagPath{" grpmagPath" , " GLO/Config/GRPMagField" , " CCDB path of the GRPMagField object" };
75- Configurable<std::string> lutPath{" lutPath" , " GLO/Param/MatLUT" , " Path of the Lut parametrization" };
76- Configurable<std::string> geoPath{" geoPath" , " GLO/Config/GeometryAligned" , " Path of the geometry file" };
77-
73+ // Configurable<std::string> ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
74+ // Configurable<std::string> grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"};
75+ // Configurable<std::string> grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"};
76+ // Configurable<std::string> lutPath{"lutPath", "GLO/Param/MatLUT", "Path of the Lut parametrization"};
77+ // Configurable<std::string> geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"};
78+
79+ // filling
80+ Configurable<bool > fillSA{" fillSA" , false , " fill spin alignment" };
7881 // events
7982 Configurable<float > cfgCutVertex{" cfgCutVertex" , 10 .0f , " Accepted z-vertex range" };
8083 Configurable<float > cfgCutCentrality{" cfgCutCentrality" , 80 .0f , " Accepted maximum Centrality" };
@@ -98,6 +101,7 @@ struct phipbpb {
98101 Configurable<int > cfgTPCcluster{" cfgTPCcluster" , 70 , " Number of TPC cluster" };
99102 Configurable<bool > isDeepAngle{" isDeepAngle" , false , " Deep Angle cut" };
100103 Configurable<bool > ispTdepPID{" ispTdepPID" , true , " pT dependent PID" };
104+ Configurable<bool > checkAllCharge{" checkAllCharge" , true , " check all charge for MC weight" };
101105 Configurable<double > cfgDeepAngle{" cfgDeepAngle" , 0.04 , " Deep Angle cut value" };
102106 Configurable<double > confRapidity{" confRapidity" , 0.5 , " Rapidity cut" };
103107 ConfigurableAxis configThnAxisInvMass{" configThnAxisInvMass" , {120 , 0.98 , 1.1 }, " #it{M} (GeV/#it{c}^{2})" };
@@ -109,6 +113,7 @@ struct phipbpb {
109113 ConfigurableAxis configThnAxisRapidity{" configThnAxisRapidity" , {8 , 0 , 0.8 }, " Rapidity" };
110114 ConfigurableAxis configThnAxisSA{" configThnAxisSA" , {200 , -1 , 1 }, " SA" };
111115 ConfigurableAxis configThnAxiscosthetaSA{" configThnAxiscosthetaSA" , {200 , 0 , 1 }, " costhetaSA" };
116+ ConfigurableAxis axisPtKaonWeight{" axisPtKaonWeight" , {VARIABLE_WIDTH, 0 .0f , 0 .1f , 0 .2f , 0 .3f , 0 .4f , 0 .5f , 0 .6f , 0 .7f , 0 .8f , 0 .9f , 1 .0f , 1 .1f , 1 .2f , 1 .3f , 1 .4f , 1 .5f , 1 .6f , 1 .7f , 1 .8f , 1 .9f , 2 .0f , 2 .2f , 2 .4f , 2 .6f , 2 .8f , 3 .0f , 3 .2f , 3 .4f , 3 .6f , 3 .8f , 4 .0f , 4 .4f , 4 .8f , 5 .2f , 5 .6f , 6 .0f , 6 .5f , 7 .0f , 7 .5f , 8 .0f , 9 .0f , 10 .0f , 11 .0f , 12 .0f }, " pt axis" };
112117 Configurable<bool > isMC{" isMC" , false , " use MC" };
113118 Configurable<bool > genacceptancecut{" genacceptancecut" , true , " use acceptance cut for generated" };
114119 Configurable<bool > avoidsplitrackMC{" avoidsplitrackMC" , false , " avoid split track in MC" };
@@ -192,12 +197,12 @@ struct phipbpb {
192197 histos.add (" hSparseV2SameEventSinDeltaPhi" , " hSparseV2SameEventSinDeltaPhi" , HistType::kTHnSparseF , {thnAxisInvMass, thnAxisPt, thnAxisV2, thnAxisCentrality});
193198 histos.add (" hSparseV2MixedEventSinDeltaPhi" , " hSparseV2MixedEventSinDeltaPhi" , HistType::kTHnSparseF , {thnAxisInvMass, thnAxisPt, thnAxisV2, thnAxisCentrality});
194199
195- histos. add ( " hSparseV2SameEventSA " , " hSparseV2SameEventSA " , HistType:: kTHnSparseF , {thnAxisInvMass, thnAxisPt, thnAxisSA, thnAxisRapidity, thnAxisCentrality});
196- histos.add (" hSparseV2MixedEventSA " , " hSparseV2MixedEventSA " , HistType::kTHnSparseF , {thnAxisInvMass, thnAxisPt, thnAxisSA, thnAxisRapidity, thnAxisCentrality});
197-
198- histos.add (" hSparseV2SameEventCosThetaStar" , " hSparseV2SameEventCosThetaStar" , HistType::kTHnSparseF , {thnAxisInvMass, thnAxisPt, thnAxisCosThetaStar, thnAxisRapidity, thnAxisCentrality});
199- histos.add (" hSparseV2MixedEventCosThetaStar" , " hSparseV2MixedEventCosThetaStar" , HistType::kTHnSparseF , {thnAxisInvMass, thnAxisPt, thnAxisCosThetaStar, thnAxisRapidity, thnAxisCentrality});
200-
200+ if (fillSA) {
201+ histos.add (" hSparseV2SameEventSA " , " hSparseV2SameEventSA " , HistType::kTHnSparseF , {thnAxisInvMass, thnAxisPt, thnAxisSA, thnAxisRapidity, thnAxisCentrality});
202+ histos. add ( " hSparseV2MixedEventSA " , " hSparseV2MixedEventSA " , HistType:: kTHnSparseF , {thnAxisInvMass, thnAxisPt, thnAxisSA, thnAxisRapidity, thnAxisCentrality});
203+ histos.add (" hSparseV2SameEventCosThetaStar" , " hSparseV2SameEventCosThetaStar" , HistType::kTHnSparseF , {thnAxisInvMass, thnAxisPt, thnAxisCosThetaStar, thnAxisRapidity, thnAxisCentrality});
204+ histos.add (" hSparseV2MixedEventCosThetaStar" , " hSparseV2MixedEventCosThetaStar" , HistType::kTHnSparseF , {thnAxisInvMass, thnAxisPt, thnAxisCosThetaStar, thnAxisRapidity, thnAxisCentrality});
205+ }
201206 // histogram for resolution
202207 histos.add (" ResFT0CTPC" , " ResFT0CTPC" , kTH3F , {centAxis, occupancyAxis, resAxis});
203208 histos.add (" ResFT0CTPCR" , " ResFT0CTPCR" , kTH3F , {centAxis, occupancyAxis, resAxis});
@@ -224,6 +229,13 @@ struct phipbpb {
224229
225230 histos.add (" hSparseV2MCRecSA" , " hSparseV2SameEventSA" , HistType::kTHnSparseF , {thnAxisInvMass, thnAxisPt, thnAxisSA, thnAxisRapidity, thnAxisCentrality});
226231 histos.add (" hSparseV2MCRecCosThetaStar_effy" , " hSparseV2SameEventCosThetaStar_effy" , HistType::kTHnSparseF , {thnAxisInvMass, thnAxisPt, thnAxisCosThetaStar, thnAxisRapidity, thnAxisCentrality});
232+
233+ // weight
234+ histos.add (" hImpactParameter" , " Impact parameter" , kTH1F , {{200 , 0 .0f , 20 .0f }});
235+ histos.add (" hEventPlaneAngle" , " hEventPlaneAngle" , kTH1F , {{100 , 0 .0f , 2 .0f * TMath::Pi ()}});
236+ histos.add (" hNchVsImpactParameter" , " hNchVsImpactParameter" , kTH2F , {{200 , 0 .0f , 20 .0f }, {500 , -0 .5f , 5000 .5f }});
237+ histos.add (" hSparseMCGenWeight" , " hSparseMCGenWeight" , HistType::kTHnSparseF , {thnAxisCentrality, {36 , 0 .0f , 2 .0f * TMath::Pi ()}, axisPtKaonWeight, {8 , -0.8 , 0.8 }});
238+ histos.add (" hSparseMCRecWeight" , " hSparseMCRecWeight" , HistType::kTHnSparseF , {thnAxisCentrality, {36 , 0 .0f , 2 .0f * TMath::Pi ()}, axisPtKaonWeight, {8 , -0.8 , 0.8 }});
227239 }
228240 // Event selection cut additional - Alex
229241 if (additionalEvsel) {
@@ -488,17 +500,19 @@ struct phipbpb {
488500 histos.fill (HIST (" hSparseV2SameEventCosPsi" ), PhiMesonMother.M (), PhiMesonMother.Pt (), TMath::Cos (2.0 * psiFT0C), centrality);
489501 histos.fill (HIST (" hSparseV2SameEventSinPsi" ), PhiMesonMother.M (), PhiMesonMother.Pt (), TMath::Sin (2.0 * psiFT0C), centrality);
490502 }
491- ROOT::Math::Boost boost{PhiMesonMother.BoostToCM ()};
492- fourVecDauCM = boost (KaonMinus);
493- threeVecDauCM = fourVecDauCM.Vect ();
494- threeVecDauCMXY = ROOT::Math::XYZVector (threeVecDauCM.X (), threeVecDauCM.Y (), 0 .);
495- eventplaneVec = ROOT::Math::XYZVector (std::cos (2.0 * psiFT0C), std::sin (2.0 * psiFT0C), 0 );
496- eventplaneVecNorm = ROOT::Math::XYZVector (std::sin (2.0 * psiFT0C), -std::cos (2.0 * psiFT0C), 0 );
497- auto cosPhistarminuspsi = GetPhiInRange (fourVecDauCM.Phi () - psiFT0C);
498- auto SA = TMath::Cos (2.0 * cosPhistarminuspsi);
499- auto cosThetaStar = eventplaneVecNorm.Dot (threeVecDauCM) / std::sqrt (threeVecDauCM.Mag2 ()) / std::sqrt (eventplaneVecNorm.Mag2 ());
500- histos.fill (HIST (" hSparseV2SameEventSA" ), PhiMesonMother.M (), PhiMesonMother.Pt (), SA, TMath::Abs (PhiMesonMother.Rapidity ()), centrality);
501- histos.fill (HIST (" hSparseV2SameEventCosThetaStar" ), PhiMesonMother.M (), PhiMesonMother.Pt (), cosThetaStar, TMath::Abs (PhiMesonMother.Rapidity ()), centrality);
503+ if (fillSA) {
504+ ROOT::Math::Boost boost{PhiMesonMother.BoostToCM ()};
505+ fourVecDauCM = boost (KaonMinus);
506+ threeVecDauCM = fourVecDauCM.Vect ();
507+ threeVecDauCMXY = ROOT::Math::XYZVector (threeVecDauCM.X (), threeVecDauCM.Y (), 0 .);
508+ eventplaneVec = ROOT::Math::XYZVector (std::cos (2.0 * psiFT0C), std::sin (2.0 * psiFT0C), 0 );
509+ eventplaneVecNorm = ROOT::Math::XYZVector (std::sin (2.0 * psiFT0C), -std::cos (2.0 * psiFT0C), 0 );
510+ auto cosPhistarminuspsi = GetPhiInRange (fourVecDauCM.Phi () - psiFT0C);
511+ auto SA = TMath::Cos (2.0 * cosPhistarminuspsi);
512+ auto cosThetaStar = eventplaneVecNorm.Dot (threeVecDauCM) / std::sqrt (threeVecDauCM.Mag2 ()) / std::sqrt (eventplaneVecNorm.Mag2 ());
513+ histos.fill (HIST (" hSparseV2SameEventSA" ), PhiMesonMother.M (), PhiMesonMother.Pt (), SA, TMath::Abs (PhiMesonMother.Rapidity ()), centrality);
514+ histos.fill (HIST (" hSparseV2SameEventCosThetaStar" ), PhiMesonMother.M (), PhiMesonMother.Pt (), cosThetaStar, TMath::Abs (PhiMesonMother.Rapidity ()), centrality);
515+ }
502516 }
503517 Npostrack = Npostrack + 1 ;
504518 }
@@ -590,17 +604,19 @@ struct phipbpb {
590604 histos.fill (HIST (" hSparseV2MixedEventCosDeltaPhiSquare" ), PhiMesonMother.M (), PhiMesonMother.Pt (), v2 * v2, centrality);
591605 histos.fill (HIST (" hSparseV2MixedEventSinDeltaPhi" ), PhiMesonMother.M (), PhiMesonMother.Pt (), v2sin * QFT0C, centrality);
592606 }
593- ROOT::Math::Boost boost{PhiMesonMother.BoostToCM ()};
594- fourVecDauCM = boost (KaonMinus);
595- threeVecDauCM = fourVecDauCM.Vect ();
596- threeVecDauCMXY = ROOT::Math::XYZVector (threeVecDauCM.X (), threeVecDauCM.Y (), 0 .);
597- eventplaneVec = ROOT::Math::XYZVector (std::cos (2.0 * psiFT0C), std::sin (2.0 * psiFT0C), 0 );
598- eventplaneVecNorm = ROOT::Math::XYZVector (std::sin (2.0 * psiFT0C), -std::cos (2.0 * psiFT0C), 0 );
599- auto cosPhistarminuspsi = GetPhiInRange (fourVecDauCM.Phi () - psiFT0C);
600- auto SA = TMath::Cos (2.0 * cosPhistarminuspsi);
601- auto cosThetaStar = eventplaneVecNorm.Dot (threeVecDauCM) / std::sqrt (threeVecDauCM.Mag2 ()) / std::sqrt (eventplaneVecNorm.Mag2 ());
602- histos.fill (HIST (" hSparseV2MixedEventSA" ), PhiMesonMother.M (), PhiMesonMother.Pt (), SA, TMath::Abs (PhiMesonMother.Rapidity ()), centrality);
603- histos.fill (HIST (" hSparseV2MixedEventCosThetaStar" ), PhiMesonMother.M (), PhiMesonMother.Pt (), cosThetaStar, TMath::Abs (PhiMesonMother.Rapidity ()), centrality);
607+ if (fillSA) {
608+ ROOT::Math::Boost boost{PhiMesonMother.BoostToCM ()};
609+ fourVecDauCM = boost (KaonMinus);
610+ threeVecDauCM = fourVecDauCM.Vect ();
611+ threeVecDauCMXY = ROOT::Math::XYZVector (threeVecDauCM.X (), threeVecDauCM.Y (), 0 .);
612+ eventplaneVec = ROOT::Math::XYZVector (std::cos (2.0 * psiFT0C), std::sin (2.0 * psiFT0C), 0 );
613+ eventplaneVecNorm = ROOT::Math::XYZVector (std::sin (2.0 * psiFT0C), -std::cos (2.0 * psiFT0C), 0 );
614+ auto cosPhistarminuspsi = GetPhiInRange (fourVecDauCM.Phi () - psiFT0C);
615+ auto SA = TMath::Cos (2.0 * cosPhistarminuspsi);
616+ auto cosThetaStar = eventplaneVecNorm.Dot (threeVecDauCM) / std::sqrt (threeVecDauCM.Mag2 ()) / std::sqrt (eventplaneVecNorm.Mag2 ());
617+ histos.fill (HIST (" hSparseV2MixedEventSA" ), PhiMesonMother.M (), PhiMesonMother.Pt (), SA, TMath::Abs (PhiMesonMother.Rapidity ()), centrality);
618+ histos.fill (HIST (" hSparseV2MixedEventCosThetaStar" ), PhiMesonMother.M (), PhiMesonMother.Pt (), cosThetaStar, TMath::Abs (PhiMesonMother.Rapidity ()), centrality);
619+ }
604620 }
605621 }
606622 }
@@ -791,6 +807,87 @@ struct phipbpb {
791807
792808 } // process MC
793809 PROCESS_SWITCH (phipbpb, processMC, " Process MC" , false );
810+
811+ using recoTracks = soa::Join<aod::TracksIU, aod::TracksExtra>;
812+ void processMCweight (aod::McCollision const & mcCollision, soa::Join<aod::McParticles, aod::ParticlesToTracks> const & mcParticles, recoTracks const &)
813+ {
814+ float imp = mcCollision.impactParameter ();
815+ float evPhi = mcCollision.eventPlaneAngle ();
816+ float centclass = -999 ;
817+ if (imp >= 0 && imp < 3.49 ) {
818+ centclass = 2.5 ;
819+ }
820+ if (imp >= 3.49 && imp < 4.93 ) {
821+ centclass = 7.5 ;
822+ }
823+ if (imp >= 4.93 && imp < 6.98 ) {
824+ centclass = 15.0 ;
825+ }
826+ if (imp >= 6.98 && imp < 8.55 ) {
827+ centclass = 25.0 ;
828+ }
829+ if (imp >= 8.55 && imp < 9.87 ) {
830+ centclass = 35.0 ;
831+ }
832+ if (imp >= 9.87 && imp < 11 ) {
833+ centclass = 45.0 ;
834+ }
835+ if (imp >= 11 && imp < 12.1 ) {
836+ centclass = 55.0 ;
837+ }
838+ if (imp >= 12.1 && imp < 13.1 ) {
839+ centclass = 65.0 ;
840+ }
841+ if (imp >= 13.1 && imp < 14 ) {
842+ centclass = 75.0 ;
843+ }
844+ if (evPhi < 0 )
845+ evPhi += 2 . * TMath::Pi ();
846+
847+ int nCh = 0 ;
848+
849+ if (centclass > 0 && centclass < 80 ) {
850+ // event within range
851+ histos.fill (HIST (" hImpactParameter" ), imp);
852+ histos.fill (HIST (" hEventPlaneAngle" ), evPhi);
853+ for (auto const & mcParticle : mcParticles) {
854+ // focus on bulk: e, mu, pi, k, p
855+ int pdgCode = TMath::Abs (mcParticle.pdgCode ());
856+ if (checkAllCharge && pdgCode != 11 && pdgCode != 13 && pdgCode != 211 && pdgCode != 321 && pdgCode != 2212 )
857+ continue ;
858+ if (pdgCode != 321 )
859+ continue ;
860+ if (!mcParticle.isPhysicalPrimary ())
861+ continue ;
862+ if (TMath::Abs (mcParticle.eta ()) > 0.8 ) // main acceptance
863+ continue ;
864+ float deltaPhi = mcParticle.phi () - mcCollision.eventPlaneAngle ();
865+ if (deltaPhi < 0 )
866+ deltaPhi += 2 . * TMath::Pi ();
867+ if (deltaPhi > 2 . * TMath::Pi ())
868+ deltaPhi -= 2 . * TMath::Pi ();
869+
870+ histos.fill (HIST (" hSparseMCGenWeight" ), centclass, deltaPhi, mcParticle.pt (), mcParticle.eta ());
871+ nCh++;
872+ bool validGlobal = false ;
873+ if (mcParticle.has_tracks ()) {
874+ auto const & tracks = mcParticle.tracks_as <recoTracks>();
875+ for (auto const & track : tracks) {
876+ if (track.hasTPC () && track.hasITS ()) {
877+ validGlobal = true ;
878+ }
879+ }
880+ }
881+ // if valid global, fill
882+ if (validGlobal) {
883+ histos.fill (HIST (" hSparseMCRecWeight" ), centclass, deltaPhi, mcParticle.pt (), mcParticle.eta ());
884+ }
885+ // if any track present, fill
886+ }
887+ }
888+ histos.fill (HIST (" hNchVsImpactParameter" ), imp, nCh);
889+ }
890+ PROCESS_SWITCH (phipbpb, processMCweight, " Process MC Weight" , false );
794891};
795892WorkflowSpec defineDataProcessing (ConfigContext const & cfgc)
796893{
0 commit comments