@@ -792,179 +792,183 @@ struct NonPromptCascadeTask {
792792 }
793793 PROCESS_SWITCH (NonPromptCascadeTask, processCascadesData, " process cascades: Data analysis" , false );
794794
795- // colls : Join<aod::Collisions, ...>
796- // tracks: Join<aod::Tracks, aod::McTrackLabels>
797- // mcCollisions: aod::McCollisions
798- // mcParticles : aod::McParticles
799-
800- void processdNdetaMC (CollisionCandidatesRun3MC const & colls,
801- aod::McCollisions const & mcCollisions,
802- aod::McParticles const & mcParticles,
803- TracksWithLabel const & tracks)
804- {
805- // -------------------------------------------------------------
806- // MC mult for all MC coll
807- // --------------------------------------------------------------
808- std::vector<int > mcMult (mcCollisions.size (), 0 );
809- for (auto const & mcp : mcParticles) {
810- int mcid = mcp.mcCollisionId ();
811- if (mcid < 0 || mcid >= (int )mcMult.size ()) continue ;
812-
813- // apply your primary/eta/charge definition here
814- if (!mcp.isPhysicalPrimary ()) continue ;
815- if (std::abs (mcp.eta ()) > 0 .5f ) continue ;
816- int q = 0 ;
817- if (auto pdg = TDatabasePDG::Instance ()->GetParticle (mcp.pdgCode ())) {
818- q = int (std::round (pdg->Charge () / 3.0 ));
819- }
820- if (q == 0 ) continue ;
821-
822- ++mcMult[mcid];
823- }
795+ // colls : Join<aod::Collisions, ...>
796+ // tracks: Join<aod::Tracks, aod::McTrackLabels>
797+ // mcCollisions: aod::McCollisions
798+ // mcParticles : aod::McParticles
799+
800+ void processdNdetaMC (CollisionCandidatesRun3MC const & colls,
801+ aod::McCollisions const & mcCollisions,
802+ aod::McParticles const & mcParticles,
803+ TracksWithLabel const & tracks)
804+ {
805+ // -------------------------------------------------------------
806+ // MC mult for all MC coll
807+ // --------------------------------------------------------------
808+ std::vector<int > mcMult (mcCollisions.size (), 0 );
809+ for (auto const & mcp : mcParticles) {
810+ int mcid = mcp.mcCollisionId ();
811+ if (mcid < 0 || mcid >= (int )mcMult.size ())
812+ continue ;
824813
825- // ------------------------------------------------------------
826- // Build mapping: (aod::Collisions row id used by tracks.collisionId())
827- // -> dense index in 'colls' (0..colls.size()-1)
828- // We assume col.globalIndex() refers to the original aod::Collisions row.
829- // ------------------------------------------------------------
830- int maxCollRowId = -1 ;
831- for (auto const & trk : tracks) {
832- maxCollRowId = std::max (maxCollRowId, (int )trk.collisionId ());
833- }
834- std::vector<int > collRowIdToDense (maxCollRowId + 1 , -1 );
814+ // apply your primary/eta/charge definition here
815+ if (!mcp.isPhysicalPrimary ())
816+ continue ;
817+ if (std::abs (mcp.eta ()) > 0 .5f )
818+ continue ;
819+ int q = 0 ;
820+ if (auto pdg = TDatabasePDG::Instance ()->GetParticle (mcp.pdgCode ())) {
821+ q = int (std::round (pdg->Charge () / 3.0 ));
822+ }
823+ if (q == 0 )
824+ continue ;
835825
836- int dense = 0 ;
837- for (auto const & col : colls) {
838- const int collRowId = col.globalIndex (); // row id in aod::Collisions
839- if (collRowId >= 0 && collRowId < (int )collRowIdToDense.size ()) {
840- collRowIdToDense[collRowId] = dense;
826+ ++mcMult[mcid];
841827 }
842- ++dense;
843- }
844828
845- // ------------------------------------------------------------
846- // Reco multiplicity per *dense collision index in colls*
847- // ------------------------------------------------------------
848- std::vector<int > recoMultDense (colls.size (), 0 );
849- for (auto const & trk : tracks) {
850- if (std::abs (trk.eta ()) > 0 .5f ) {
851- continue ;
852- }
853- const int collRowId = trk.collisionId ();
854- if (collRowId < 0 || collRowId >= (int )collRowIdToDense.size ()) {
855- continue ;
856- }
857- const int dIdx = collRowIdToDense[collRowId];
858- if (dIdx >= 0 ) {
859- ++recoMultDense[dIdx];
829+ // ------------------------------------------------------------
830+ // Build mapping: (aod::Collisions row id used by tracks.collisionId())
831+ // -> dense index in 'colls' (0..colls.size()-1)
832+ // We assume col.globalIndex() refers to the original aod::Collisions row.
833+ // ------------------------------------------------------------
834+ int maxCollRowId = -1 ;
835+ for (auto const & trk : tracks) {
836+ maxCollRowId = std::max (maxCollRowId, (int )trk.collisionId ());
860837 }
861- }
838+ std::vector< int > collRowIdToDense (maxCollRowId + 1 , - 1 );
862839
863- // ------------------------------------------------------------
864- // MC bookkeeping: index by ROW INDEX (0..size-1), not globalIndex()
865- // ------------------------------------------------------------
866- std::vector<char > isReco (mcParticles.size (), 0 );
867- std::vector<float > isRecoMult (mcParticles.size (), 0 .f );
868- std::vector<char > mcReconstructed (mcCollisions.size (), 0 );
869-
870- // Optional cache of MC multiplicity per MC collision
871- std::vector<float > mcMultCache (mcCollisions.size (), -1 .f );
872-
873- // ------------------------------------------------------------
874- // Single pass over tracks: fill RM for tracks whose collision is in colls
875- // ------------------------------------------------------------
876- for (auto const & trk : tracks) {
877- // Accept reco track
878- if (std::abs (trk.eta ()) > 0 .5f ) {
879- continue ;
840+ int dense = 0 ;
841+ for (auto const & col : colls) {
842+ const int collRowId = col.globalIndex (); // row id in aod::Collisions
843+ if (collRowId >= 0 && collRowId < (int )collRowIdToDense.size ()) {
844+ collRowIdToDense[collRowId] = dense;
845+ }
846+ ++dense;
880847 }
881848
882- // Map track's collision row id -> dense colls index
883- const int collRowId = trk.collisionId ();
884- if (collRowId < 0 || collRowId >= (int )collRowIdToDense.size ()) {
885- continue ;
886- }
887- const int dIdx = collRowIdToDense[collRowId];
888- if (dIdx < 0 ) {
889- continue ; // this track's collision is not in our 'colls' view
849+ // ------------------------------------------------------------
850+ // Reco multiplicity per *dense collision index in colls*
851+ // ------------------------------------------------------------
852+ std::vector<int > recoMultDense (colls.size (), 0 );
853+ for (auto const & trk : tracks) {
854+ if (std::abs (trk.eta ()) > 0 .5f ) {
855+ continue ;
856+ }
857+ const int collRowId = trk.collisionId ();
858+ if (collRowId < 0 || collRowId >= (int )collRowIdToDense.size ()) {
859+ continue ;
860+ }
861+ const int dIdx = collRowIdToDense[collRowId];
862+ if (dIdx >= 0 ) {
863+ ++recoMultDense[dIdx];
864+ }
890865 }
891866
892- // Get the collision row (dense index in colls view)
893- auto col = colls.rawIteratorAt (dIdx);
867+ // ------------------------------------------------------------
868+ // MC bookkeeping: index by ROW INDEX (0..size-1), not globalIndex()
869+ // ------------------------------------------------------------
870+ std::vector<char > isReco (mcParticles.size (), 0 );
871+ std::vector<float > isRecoMult (mcParticles.size (), 0 .f );
872+ std::vector<char > mcReconstructed (mcCollisions.size (), 0 );
894873
895- // MC collision id (row index in aod::McCollisions)
896- const int mcCollId = col.mcCollisionId ();
897- if (mcCollId < 0 || mcCollId >= (int )mcCollisions.size ()) {
898- continue ;
899- }
900- mcReconstructed[mcCollId] = 1 ;
874+ // Optional cache of MC multiplicity per MC collision
875+ std::vector<float > mcMultCache (mcCollisions.size (), -1 .f );
901876
902- // MC particle id (row index in aod::McParticles)
903- const int mcPid = trk.mcParticleId ();
904- if (mcPid < 0 || mcPid >= (int )mcParticles.size ()) {
905- continue ;
906- }
877+ // ------------------------------------------------------------
878+ // Single pass over tracks: fill RM for tracks whose collision is in colls
879+ // ------------------------------------------------------------
880+ for (auto const & trk : tracks) {
881+ // Accept reco track
882+ if (std::abs (trk.eta ()) > 0 .5f ) {
883+ continue ;
884+ }
907885
908- // MC multiplicity for that MC collision (cache)
909- float mult = mcMultCache[mcCollId];
910- if (mult < 0 .f ) {
911- std::vector<float > tmp;
912- mult = mcMult[mcCollId];
913- mcMultCache[mcCollId] = mult;
914- }
886+ // Map track's collision row id -> dense colls index
887+ const int collRowId = trk.collisionId ();
888+ if (collRowId < 0 || collRowId >= (int )collRowIdToDense.size ()) {
889+ continue ;
890+ }
891+ const int dIdx = collRowIdToDense[collRowId];
892+ if (dIdx < 0 ) {
893+ continue ; // this track's collision is not in our 'colls' view
894+ }
915895
916- auto mcPar = mcParticles.rawIteratorAt (mcPid);
896+ // Get the collision row (dense index in colls view)
897+ auto col = colls.rawIteratorAt (dIdx);
917898
918- // Apply the same acceptance as in MC multiplicity definition
919- if (!mcPar.isPhysicalPrimary ()) {
920- continue ;
921- }
922- if (std::abs (mcPar.eta ()) > 0 .5f ) {
923- continue ;
924- }
899+ // MC collision id (row index in aod::McCollisions)
900+ const int mcCollId = col.mcCollisionId ();
901+ if (mcCollId < 0 || mcCollId >= (int )mcCollisions.size ()) {
902+ continue ;
903+ }
904+ mcReconstructed[mcCollId] = 1 ;
925905
926- int q = 0 ;
927- if (auto pdgEntry = TDatabasePDG::Instance ()->GetParticle (mcPar.pdgCode ())) {
928- q = int (std::round (pdgEntry->Charge () / 3.0 ));
929- }
930- if (q == 0 ) {
931- continue ;
932- }
906+ // MC particle id (row index in aod::McParticles)
907+ const int mcPid = trk.mcParticleId ();
908+ if (mcPid < 0 || mcPid >= (int )mcParticles.size ()) {
909+ continue ;
910+ }
933911
934- // Mark reconstructed MC particle (now that it truly passed & matched)
935- isReco[mcPid] = 1 ;
936- isRecoMult[mcPid] = mult;
912+ // MC multiplicity for that MC collision (cache)
913+ float mult = mcMultCache[mcCollId];
914+ if (mult < 0 .f ) {
915+ std::vector<float > tmp;
916+ mult = mcMult[mcCollId];
917+ mcMultCache[mcCollId] = mult;
918+ }
937919
938- const float multReco = col.multNTracksGlobal (); // or recoMultDense[dIdx]
939- const float ptReco = trk.pt ();
940- const float ptMC = mcPar.pt ();
920+ auto mcPar = mcParticles.rawIteratorAt (mcPid);
941921
942- mRegistrydNdeta .fill (HIST (" hdNdetaRM/hdNdetaRM" ), mult, multReco, ptMC, ptReco);
943- }
922+ // Apply the same acceptance as in MC multiplicity definition
923+ if (!mcPar.isPhysicalPrimary ()) {
924+ continue ;
925+ }
926+ if (std::abs (mcPar.eta ()) > 0 .5f ) {
927+ continue ;
928+ }
929+
930+ int q = 0 ;
931+ if (auto pdgEntry = TDatabasePDG::Instance ()->GetParticle (mcPar.pdgCode ())) {
932+ q = int (std::round (pdgEntry->Charge () / 3.0 ));
933+ }
934+ if (q == 0 ) {
935+ continue ;
936+ }
937+
938+ // Mark reconstructed MC particle (now that it truly passed & matched)
939+ isReco[mcPid] = 1 ;
940+ isRecoMult[mcPid] = mult;
944941
945- // ------------------------------------------------------------
946- // MC particles with no reco track (iterate by row index)
947- // ------------------------------------------------------------
948- for (int pid = 0 ; pid < (int )mcParticles.size (); ++pid) {
949- if (!isReco[pid]) {
950- auto mcp = mcParticles.rawIteratorAt (pid);
951- mRegistrydNdeta .fill (HIST (" hdNdetaRM/hdNdetaRMNotInRecoTrk" ), isRecoMult[pid], mcp.pt ());
942+ const float multReco = col.multNTracksGlobal (); // or recoMultDense[dIdx]
943+ const float ptReco = trk.pt ();
944+ const float ptMC = mcPar.pt ();
945+
946+ mRegistrydNdeta .fill (HIST (" hdNdetaRM/hdNdetaRM" ), mult, multReco, ptMC, ptReco);
952947 }
953- }
954948
955- // ------------------------------------------------------------
956- // Unreconstructed MC collisions (iterate by row index)
957- // ------------------------------------------------------------
958- for (int mcid = 0 ; mcid < (int )mcCollisions.size (); ++mcid) {
959- if (!mcReconstructed[mcid]) {
960- std::vector<float > mcptvec;
961- const int mult = mcMult[mcid];
962- for (auto const & pt : mcptvec) {
963- mRegistrydNdeta .fill (HIST (" hdNdetaRM/hdNdetaRMNotInRecoCol" ), mult, pt);
949+ // ------------------------------------------------------------
950+ // MC particles with no reco track (iterate by row index)
951+ // ------------------------------------------------------------
952+ for (int pid = 0 ; pid < (int )mcParticles.size (); ++pid) {
953+ if (!isReco[pid]) {
954+ auto mcp = mcParticles.rawIteratorAt (pid);
955+ mRegistrydNdeta .fill (HIST (" hdNdetaRM/hdNdetaRMNotInRecoTrk" ), isRecoMult[pid], mcp.pt ());
956+ }
957+ }
958+
959+ // ------------------------------------------------------------
960+ // Unreconstructed MC collisions (iterate by row index)
961+ // ------------------------------------------------------------
962+ for (int mcid = 0 ; mcid < (int )mcCollisions.size (); ++mcid) {
963+ if (!mcReconstructed[mcid]) {
964+ std::vector<float > mcptvec;
965+ const int mult = mcMult[mcid];
966+ for (auto const & pt : mcptvec) {
967+ mRegistrydNdeta .fill (HIST (" hdNdetaRM/hdNdetaRMNotInRecoCol" ), mult, pt);
968+ }
964969 }
965970 }
966971 }
967- }
968972
969973 PROCESS_SWITCH (NonPromptCascadeTask, processdNdetaMC, " process mc dN/deta" , false );
970974};
0 commit comments