Skip to content

Commit 3672228

Browse files
bghanley1995bghanley1995
committed
Added track xy and yz plots
1 parent 32b9733 commit 3672228

File tree

2 files changed

+139
-58
lines changed

2 files changed

+139
-58
lines changed

‎PWGCF/TwoParticleCorrelations/TableProducer/identifiedBfFilter.cxx‎

Lines changed: 90 additions & 56 deletions
Original file line numberDiff line numberDiff line change
@@ -80,6 +80,10 @@ TH1F* fhVertexZB = nullptr;
8080
TH1F* fhVertexZA = nullptr;
8181
TH1F* fhMultB = nullptr;
8282
TH1F* fhMultA = nullptr;
83+
TH2F* fhYZB = nullptr;
84+
TH2F* fhXYB = nullptr;
85+
TH2F* fhYZA = nullptr;
86+
TH2F* fhXYA = nullptr;
8387
TH1F* fhPB = nullptr;
8488
TH1F* fhPA[kIdBfNoOfSpecies + 1] = {nullptr};
8589
TH1F* fhPtB = nullptr;
@@ -134,6 +138,7 @@ TH1F* fhTrueVertexZB = nullptr;
134138
TH1F* fhTrueVertexZA = nullptr;
135139
TH1F* fhTrueVertexZAA = nullptr;
136140
TH1F* fhTruePB = nullptr;
141+
TH1F* fhTrueCharge = nullptr;
137142
TH1F* fhTruePA[kIdBfNoOfSpecies + 1] = {nullptr};
138143
TH1F* fhTruePtB = nullptr;
139144
TH1F* fhTruePtA[kIdBfNoOfSpecies + 1] = {nullptr};
@@ -144,6 +149,7 @@ TH1F* fhTruePtNegA[kIdBfNoOfSpecies + 1] = {nullptr};
144149
TH2F* fhTrueNPosNegA[kIdBfNoOfSpecies + 1] = {nullptr};
145150
TH1F* fhTrueDeltaNA[kIdBfNoOfSpecies + 1] = {nullptr};
146151

152+
147153
TH1F* fhTrueEtaB = nullptr;
148154
TH1F* fhTrueEtaA = nullptr;
149155

@@ -744,6 +750,10 @@ struct IdentifiedBfFilterTracks {
744750

745751
if ((fDataType == kData) || (fDataType == kDataNoEvtSel) || (fDataType == kMC)) {
746752
/* create the reconstructed data histograms */
753+
fhXYB = new TH2F("fHistXYB", "x and y distribution for reconstructed before", 1000, -10.0, 10.0, 1000, -10.0, 10.0);
754+
fhYZB = new TH2F("fHistYZB", "y and z distribution for reconstructed before", 1000, -10.0, 10.0, 1000, -10.0, 10.0);
755+
fhXYA = new TH2F("fHistXYA", "x and y distribution for reconstructed after", 1000, -10.0, 10.0, 1000, -10.0, 10.0);
756+
fhYZA = new TH2F("fHistYZA", "y and z distribution for reconstructed after", 1000, -10.0, 10.0, 1000, -10.0, 10.0);
747757
fhPB = new TH1F("fHistPB", "p distribution for reconstructed before;p (GeV/c);dN/dp (c/GeV)", 100, 0.0, 15.0);
748758
fhPtB = new TH1F("fHistPtB", "p_{T} distribution for reconstructed before;p_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0);
749759
fhPtPosB = new TH1F("fHistPtPosB", "P_{T} distribution for reconstructed (#plus) before;P_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0);
@@ -825,6 +835,10 @@ struct IdentifiedBfFilterTracks {
825835
TString::Format("dE/dx vs P_{IP} reconstructed Wrong Species; P (GeV/c); dE/dx (a.u.)").Data(),
826836
ptbins, ptlow, ptup, 1000, 0.0, 1000.0);
827837
/* add the hstograms to the output list */
838+
fOutputList->Add(fhXYB);
839+
fOutputList->Add(fhYZB);
840+
fOutputList->Add(fhXYA);
841+
fOutputList->Add(fhYZA);
828842
fOutputList->Add(fhPB);
829843
fOutputList->Add(fhPtB);
830844
fOutputList->Add(fhPtPosB);
@@ -877,7 +891,9 @@ struct IdentifiedBfFilterTracks {
877891

878892
if ((fDataType != kData) && (fDataType != kDataNoEvtSel)) {
879893
/* create the true data histograms */
894+
880895
fhTruePB = new TH1F("fTrueHistPB", "p distribution before (truth);p (GeV/c);dN/dp (c/GeV)", 100, 0.0, 15.0);
896+
fhTrueCharge = new TH1F("fTrueHistCharge", "Charge distribution before (truth);charge;count", 3, -1.0, 1.0);
881897
fhTruePtB = new TH1F("fTrueHistPtB", "p_{T} distribution before (truth);p_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0);
882898
fhTruePtPosB = new TH1F("fTrueHistPtPosB", "P_{T} distribution (#plus) before (truth);P_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0);
883899
fhTruePtNegB = new TH1F("fTrueHistPtNegB", "P_{T} distribution (#minus) before (truth);P_{T} (GeV/c);dN/dP_{T} (c/GeV)", 100, 0.0, 15.0);
@@ -921,6 +937,7 @@ struct IdentifiedBfFilterTracks {
921937
fOutputList->Add(fhTruePtB);
922938
fOutputList->Add(fhTruePtPosB);
923939
fOutputList->Add(fhTruePtNegB);
940+
fOutputList->Add(fhTrueCharge);
924941
fOutputList->Add(fhTrueEtaB);
925942
fOutputList->Add(fhTrueEtaA);
926943
fOutputList->Add(fhTruePhiB);
@@ -947,18 +964,18 @@ struct IdentifiedBfFilterTracks {
947964

948965
template <typename TrackObject>
949966
inline MatchRecoGenSpecies IdentifyTrack(TrackObject const& track);
950-
template <typename TrackObject>
951-
int8_t AcceptTrack(TrackObject const& track);
967+
template <typename TrackObject, typename CollisionObject>
968+
int8_t AcceptTrack(TrackObject const& track, CollisionObject const& collision);
952969
template <typename ParticleObject, typename MCCollisionObject>
953970
int8_t AcceptParticle(ParticleObject& particle, MCCollisionObject const& mccollision);
954971
template <typename CollisionObjects, typename TrackObject>
955972
int8_t selectTrackAmbiguousCheck(CollisionObjects const& collisions, TrackObject const& track);
956973
template <typename ParticleObject>
957974
inline MatchRecoGenSpecies IdentifyParticle(ParticleObject const& particle);
958-
template <typename TrackObject>
959-
void fillTrackHistosBeforeSelection(TrackObject const& track);
960-
template <typename TrackObject>
961-
void fillTrackHistosAfterSelection(TrackObject const& track, MatchRecoGenSpecies sp);
975+
template <typename TrackObject, typename CollisionObject>
976+
void fillTrackHistosBeforeSelection(TrackObject const& track, CollisionObject const& collision);
977+
template <typename TrackObject, typename CollisionObject>
978+
void fillTrackHistosAfterSelection(TrackObject const& track, MatchRecoGenSpecies sp, CollisionObject const& collision);
962979
template <typename ParticleObject, typename MCCollisionObject>
963980
void fillParticleHistosBeforeSelection(ParticleObject const& particle,
964981
MCCollisionObject const& collision,
@@ -1044,36 +1061,42 @@ struct IdentifiedBfFilterTracks {
10441061
}
10451062

10461063
for (auto& particle : particles) {
1047-
1048-
10491064
int8_t pid = -1;
1050-
TParticlePDG* pdgpart = fPDG->GetParticle(particle.pdgCode());
1051-
float charge = 0;
1052-
if (pdgpart != nullptr){
1053-
charge = getCharge(pdgpart->Charge());
1054-
//print charge
1055-
}
1056-
1057-
if (charge != 0) {
1058-
if (particle.has_mcCollision() && (particle.template mcCollision_as<soa::Join<aod::McCollisions, aod::IdentifiedBfCFGenCollisionsInfo>>()).collisionaccepted()) {
1059-
auto mccollision = particle.template mcCollision_as<soa::Join<aod::McCollisions, aod::IdentifiedBfCFGenCollisionsInfo>>();
1060-
/* before particle selection */
1061-
fillParticleHistosBeforeSelection(particle, mccollision, charge);
1062-
1063-
/* track selection */
1064-
pid = AcceptParticle(particle, mccollision);
1065-
if (!(pid < 0)) { // if PID isn't negative
1066-
acceptedparticles++;
1065+
if(particle.isPhysicalPrimary()){
1066+
TParticlePDG* pdgpart = fPDG->GetParticle(particle.pdgCode());
1067+
float charge = 0;
1068+
if (pdgpart != nullptr){
1069+
charge = getCharge(pdgpart->Charge());
1070+
//print charge
1071+
}
1072+
fhTrueCharge->Fill(charge);
1073+
if (charge != 0) {
1074+
if (particle.has_mcCollision() && (particle.template mcCollision_as<soa::Join<aod::McCollisions, aod::IdentifiedBfCFGenCollisionsInfo>>()).collisionaccepted()) {
1075+
auto mccollision = particle.template mcCollision_as<soa::Join<aod::McCollisions, aod::IdentifiedBfCFGenCollisionsInfo>>();
1076+
/* before particle selection */
1077+
fillParticleHistosBeforeSelection(particle, mccollision, charge);
1078+
1079+
/* track selection */
1080+
pid = AcceptParticle(particle, mccollision);
1081+
if (!(pid < 0)) { // if PID isn't negative
1082+
acceptedparticles++;
1083+
}
1084+
}
1085+
} else {
1086+
if ((particle.mcCollisionId() == 0) && traceCollId0) {
1087+
LOGF(IDENTIFIEDBFFILTERLOGTRACKS, "Particle %d with fractional charge or equal to zero", particle.globalIndex());
10671088
}
10681089
}
1090+
10691091
} else {
1070-
if ((particle.mcCollisionId() == 0) && traceCollId0) {
1071-
LOGF(IDENTIFIEDBFFILTERLOGTRACKS, "Particle %d with fractional charge or equal to zero", particle.globalIndex());
1092+
if((particle.mcCollisionId() == 0) && traceCollId0){
1093+
LOGF(IDENTIFIEDBFFILTERLOGTRACKS, "Particle %d not Physical Primary", particle.globalIndex());
10721094
}
10731095
}
10741096
if (!fullDerivedData) {
10751097
gentracksinfo(pid);
10761098
}
1099+
10771100
}
10781101
LOGF(info,
10791102
"Processed %d accepted generated collisions out of a total of %d with %d accepted particles out of a "
@@ -1294,7 +1317,7 @@ inline MatchRecoGenSpecies IdentifiedBfFilterTracks::IdentifyTrack(TrackObject c
12941317
float min_nsigma = 999.0f;
12951318
MatchRecoGenSpecies sp_min_nsigma = kWrongSpecies;
12961319
for (int sp = 0; sp < kIdBfNoOfSpecies; ++sp) {
1297-
if (abs(nsigmas[sp]) < abs(min_nsigma)) { // Check if species nsigma is less than current nsigma
1320+
if (fabs(nsigmas[sp]) < fabs(min_nsigma)) { // Check if species nsigma is less than current nsigma
12981321
min_nsigma = nsigmas[sp]; // If yes, set species nsigma to current nsigma
12991322
sp_min_nsigma = MatchRecoGenSpecies(sp); // set current species sp number to current sp
13001323
}
@@ -1344,10 +1367,10 @@ inline MatchRecoGenSpecies IdentifiedBfFilterTracks::IdentifyTrack(TrackObject c
13441367
/// For the time being we keep the convention
13451368
/// - positive track pid even
13461369
/// - negative track pid odd
1347-
template <typename TrackObject>
1348-
inline int8_t IdentifiedBfFilterTracks::AcceptTrack(TrackObject const& track)
1370+
template <typename TrackObject, typename CollisionObject>
1371+
inline int8_t IdentifiedBfFilterTracks::AcceptTrack(TrackObject const& track, CollisionObject const& collision)
13491372
{
1350-
fillTrackHistosBeforeSelection(track); // <Fill "before selection" histo
1373+
fillTrackHistosBeforeSelection(track, collision); // <Fill "before selection" histo
13511374

13521375
/* TODO: incorporate a mask in the scanned tracks table for the rejecting track reason */
13531376
if constexpr (framework::has_type_v<aod::mctracklabel::McParticleId, typename TrackObject::all_columns>) {
@@ -1356,9 +1379,9 @@ inline int8_t IdentifiedBfFilterTracks::AcceptTrack(TrackObject const& track)
13561379
}
13571380
}
13581381

1359-
if (matchTrackType(track)) {
1382+
if (matchTrackType(track, collision)) {
13601383
if (ptlow < track.pt() && track.pt() < ptup && etalow < track.eta() && track.eta() < etaup) {
1361-
fillTrackHistosAfterSelection(track, kIdBfCharged);
1384+
fillTrackHistosAfterSelection(track, kIdBfCharged, collision);
13621385
MatchRecoGenSpecies sp = kWrongSpecies;
13631386
if (recoIdMethod == 0) {
13641387
sp = kIdBfCharged;
@@ -1380,7 +1403,7 @@ inline int8_t IdentifiedBfFilterTracks::AcceptTrack(TrackObject const& track)
13801403
return -1;
13811404
}
13821405
if (!(sp < 0)) {
1383-
fillTrackHistosAfterSelection(track, sp); //<Fill accepted track histo with PID
1406+
fillTrackHistosAfterSelection(track, sp, collision); //<Fill accepted track histo with PID
13841407
if (track.sign() > 0) { // if positive
13851408
trkMultPos[sp]++; //<< Update Particle Multiplicity
13861409
return speciesChargeValue1[sp];
@@ -1504,14 +1527,15 @@ int8_t IdentifiedBfFilterTracks::selectTrackAmbiguousCheck(CollisionObjects cons
15041527
/* feedback of no ambiguous tracks only if checks required */
15051528
fhAmbiguousTrackType->Fill(tracktype, multiplicityclass);
15061529
}
1507-
std::cout<<"Hello Accept Track"<<std::endl;
1508-
return AcceptTrack(track);
1530+
return AcceptTrack(track, collisions.iteratorAt(track.collisionId()));
15091531
}
15101532
}
15111533

1512-
template <typename TrackObject>
1513-
void IdentifiedBfFilterTracks::fillTrackHistosBeforeSelection(TrackObject const& track)
1534+
template <typename TrackObject, typename CollisionObject>
1535+
void IdentifiedBfFilterTracks::fillTrackHistosBeforeSelection(TrackObject const& track, CollisionObject const& collision)
15141536
{
1537+
fhXYB->Fill(track.x(),track.y());
1538+
fhYZB->Fill(track.y(),track.z());
15151539
fhPB->Fill(track.p());
15161540
fhPtB->Fill(track.pt());
15171541
fhEtaB->Fill(track.eta());
@@ -1523,17 +1547,28 @@ void IdentifiedBfFilterTracks::fillTrackHistosBeforeSelection(TrackObject const&
15231547
} else {
15241548
fhPtNegB->Fill(track.pt());
15251549
}
1550+
1551+
//float dcaxy = CalculateDCA(track, collision,0);
1552+
//float dcaz = CalculateDCA(track, collision,1);
1553+
//fhDCAxyB->Fill(dcaxy);
1554+
//fhDCAzB->Fill(dcaz);
15261555
fhDCAxyB->Fill(track.dcaXY());
15271556
fhDCAzB->Fill(track.dcaZ());
15281557
}
15291558

1530-
template <typename TrackObject>
1531-
void IdentifiedBfFilterTracks::fillTrackHistosAfterSelection(TrackObject const& track, MatchRecoGenSpecies sp)
1559+
template <typename TrackObject, typename CollisionObject>
1560+
void IdentifiedBfFilterTracks::fillTrackHistosAfterSelection(TrackObject const& track, MatchRecoGenSpecies sp, CollisionObject const& collision)
15321561
{
15331562
/* the charged species should have been called first so avoid double counting */
15341563
if (sp == kIdBfCharged) {
15351564
fhEtaA->Fill(track.eta());
15361565
fhPhiA->Fill(track.phi());
1566+
fhXYA->Fill(track.x(),track.y());
1567+
fhYZA->Fill(track.y(),track.z());
1568+
//float dcaxy = CalculateDCA(track, collision, 0);
1569+
//float dcaz = CalculateDCA(track, collision, 1);
1570+
//fhDCAxyA->Fill(dcaxy);
1571+
//fhDCAzA->Fill(dcaz);
15371572
fhDCAxyA->Fill(track.dcaXY());
15381573
fhDCAzA->Fill(track.dcaZ());
15391574

@@ -1543,16 +1578,15 @@ void IdentifiedBfFilterTracks::fillTrackHistosAfterSelection(TrackObject const&
15431578
if (track.dcaZ() < 1.0) {
15441579
fhFineDCAzA->Fill(track.dcaZ());
15451580
}
1581+
}
1582+
fhPA[sp]->Fill(track.p());
1583+
fhPtA[sp]->Fill(track.pt());
1584+
fhdEdxA[sp]->Fill(track.p(), track.tpcSignal());
1585+
fhdEdxIPTPCA[sp]->Fill(track.tpcInnerParam(), track.tpcSignal());
1586+
if (track.sign() > 0) {
1587+
fhPtPosA[sp]->Fill(track.pt());
15461588
} else {
1547-
fhPA[sp]->Fill(track.p());
1548-
fhPtA[sp]->Fill(track.pt());
1549-
fhdEdxA[sp]->Fill(track.p(), track.tpcSignal());
1550-
fhdEdxIPTPCA[sp]->Fill(track.tpcInnerParam(), track.tpcSignal());
1551-
if (track.sign() > 0) {
1552-
fhPtPosA[sp]->Fill(track.pt());
1553-
} else {
1554-
fhPtNegA[sp]->Fill(track.pt());
1555-
}
1589+
fhPtNegA[sp]->Fill(track.pt());
15561590
}
15571591
}
15581592

@@ -1598,15 +1632,15 @@ void IdentifiedBfFilterTracks::fillParticleHistosAfterSelection(ParticleObject c
15981632
fhTrueDCAxyA->Fill(TMath::Sqrt((particle.vx() - collision.posX()) * (particle.vx() - collision.posX()) +
15991633
(particle.vy() - collision.posY()) * (particle.vy() - collision.posY())));
16001634
fhTrueDCAzA->Fill((particle.vz() - collision.posZ()));
1635+
}
1636+
fhTruePA[sp]->Fill(particle.p());
1637+
fhTruePtA[sp]->Fill(particle.pt());
1638+
if (charge > 0) {
1639+
fhTruePtPosA[sp]->Fill(particle.pt());
16011640
} else {
1602-
fhTruePA[sp]->Fill(particle.p());
1603-
fhTruePtA[sp]->Fill(particle.pt());
1604-
if (charge > 0) {
1605-
fhTruePtPosA[sp]->Fill(particle.pt());
1606-
} else {
1607-
fhTruePtNegA[sp]->Fill(particle.pt());
1608-
}
1641+
fhTruePtNegA[sp]->Fill(particle.pt());
16091642
}
1643+
16101644
}
16111645

16121646
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

‎PWGCF/TwoParticleCorrelations/TableProducer/identifiedBfFilter.h‎

Lines changed: 49 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,7 @@
2626
#include "Common/Core/TrackSelection.h"
2727
#include "Common/Core/TrackSelectionDefaults.h"
2828
#include "PWGCF/Core/AnalysisConfigurableCuts.h"
29+
#include "MathUtils/Utils.h"
2930
#include <TDatabasePDG.h>
3031

3132
namespace o2
@@ -702,16 +703,62 @@ inline bool IsEvtSelected(CollisionObject const& collision, float& centormult)
702703
//////////////////////////////////////////////////////////////////////////////////
703704
/// Track selection
704705
//////////////////////////////////////////////////////////////////////////////////
706+
template <typename TrackObject, typename CollisionObject>
707+
float CalculateDCA(TrackObject const& track,CollisionObject const& collision, int vtx)
708+
{
709+
float dca;
710+
// propagate track to DCA to the vertex
711+
float sn, cs;
712+
auto alp = track.alpha();
713+
math_utils::detail::sincos<float>(alp, sn, cs);
714+
auto x = track.x(), y = track.y(), snp = track.snp(), csp = math_utils::detail::sqrt<float>((1.f - snp) * (1.f + snp));
715+
auto xv = collision.posX() * cs + collision.posX() * sn, yv = -collision.posX() * sn + collision.posY() * cs, zv = collision.posZ();
716+
717+
x -= xv;
718+
y -= yv;
719+
// Estimate the impact parameter neglecting the track curvature
720+
721+
auto crv = 0;
722+
auto tgfv = -(crv * x - snp) / (crv * y + csp);
723+
sn = tgfv / math_utils::detail::sqrt<float>(1.f + tgfv * tgfv);
724+
cs = math_utils::detail::sqrt<float>((1. - sn) * (1. + sn));
725+
cs = (math_utils::detail::abs<float>(tgfv) > o2::constants::math::Almost0) ? sn / tgfv : o2::constants::math::Almost1;
726+
727+
x = xv * cs + yv * sn;
728+
//yv = -xv * sn + yv * cs;
729+
xv = x;
730+
731+
math_utils::detail::sincos<float>(alp, sn, cs);
732+
733+
if(vtx == 0){
734+
if (fabs(track.x())>fabs(track.y())){
735+
dca = track.y() - yv;
736+
} else{
737+
dca = track.x() - xv;
738+
}
739+
if(dca>10.0){
740+
std::cout<<"Big dca: "<<dca<<std::endl;
741+
std::cout<<"x: "<<track.x()<<std::endl;
742+
std::cout<<"y: "<<track.y()<<std::endl;
743+
std::cout<<"yv: "<<yv<<std::endl;
744+
}
745+
} else if(vtx == 1){
746+
dca = track.z() - zv;
747+
}
748+
return dca;
749+
}
705750

706-
template <typename TrackObject>
707-
inline bool matchTrackType(TrackObject const& track)
751+
template <typename TrackObject, typename CollisionObject>
752+
inline bool matchTrackType(TrackObject const& track, CollisionObject const& collision)
708753
{
709754
if (useOwnTrackSelection) {
710755
return ownTrackSelection.IsSelected(track);
711756
} else {
712757
for (auto filter : trackFilters) {
713758
if (filter->IsSelected(track)) {
714759
if (dca2Dcut) {
760+
//float dcaxy = CalculateDCA(track, collision, 0);
761+
//float dcaz = CalculateDCA(track, collision, 1);
715762
if (track.dcaXY() * track.dcaXY() / maxDCAxy / maxDCAxy + track.dcaZ() * track.dcaZ() / maxDCAz / maxDCAz > 1) {
716763
return false;
717764
} else {

0 commit comments

Comments
 (0)