@@ -37,13 +37,6 @@ namespace femtoDream
3737// / \tparam partOne Type of particle 1 (Track/V0/Cascade/...)
3838// / \tparam partTwo Type of particle 2 (Track/V0/Cascade/...)
3939
40- enum ProngCharmHadron {
41- Prong0 = 0 ,
42- Prong1 = 1 ,
43- Prong2 = 2 ,
44- Nprongs = 3
45- };
46-
4740template <o2::aod::femtodreamparticle::ParticleType partOne, o2::aod::femtodreamparticle::ParticleType partTwo>
4841class FemtoDreamDetaDphiStar
4942{
@@ -118,8 +111,14 @@ class FemtoDreamDetaDphiStar
118111 }
119112 }
120113 }
121- if constexpr (mPartOneType == o2::aod::femtodreamparticle::ParticleType::kTrack && mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kCharmHadron ) {
122- for (int i = 0 ; i < Nprongs; i++) {
114+ if constexpr (mPartOneType == o2::aod::femtodreamparticle::ParticleType::kTrack && (mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kCharmHadron3Prong || mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kCharmHadron2Prong )) {
115+ int nProng = 0 ;
116+ if (mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kCharmHadron3Prong )
117+ nProng = 3 ;
118+ if (mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kCharmHadron2Prong )
119+ nProng = 2 ;
120+
121+ for (int i = 0 ; i < nProng; i++) {
123122 std::string dirName = static_cast <std::string>(dirNames[2 ]);
124123 histdetadpi[i][0 ] = mHistogramRegistry ->add <TH2>((dirName + static_cast <std::string>(histNames[0 ][i]) + static_cast <std::string>(histNameSEorME[meORse])).c_str (), " ; #Delta #eta; #Delta #phi^{*}" kTH2F , {{100 , -0.15 , 0.15 }, {100 , -0.15 , 0.15 }});
125124 histdetadpi[i][1 ] = mHistogramRegistry ->add <TH2>((dirName + static_cast <std::string>(histNames[1 ][i]) + static_cast <std::string>(histNameSEorME[meORse])).c_str (), " ; #Delta #eta; #Delta #phi^{*}" kTH2F , {{100 , -0.15 , 0.15 }, {100 , -0.15 , 0.15 }});
@@ -500,50 +499,102 @@ class FemtoDreamDetaDphiStar
500499 }
501500
502501 return pass;
503- } else if constexpr (mPartOneType == o2::aod::femtodreamparticle::ParticleType::kTrack && mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kCharmHadron ) {
502+ } else if constexpr (mPartOneType == o2::aod::femtodreamparticle::ParticleType::kTrack && ( mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kCharmHadron3Prong || mPartTwoType == o2::aod::femtodreamparticle::ParticleType:: kCharmHadronDstar )) { // to be fixed for Dstar
504503 // check if provided particles are in agreement with the class instantiation
505- if (part2.candidateSelFlag () < o2::aod::fdhf::lcToPKPi ) {
504+ if (! part2.candidateSelFlag ()) {
506505 LOG (fatal) << " FemtoDreamDetaDphiStar: passed arguments don't agree with FemtoDreamDetaDphiStar instantiation! Please provide Charm Hadron candidates."
507506 return false ;
508507 }
509508
510509 bool pass = false ;
511510
512- for (int i = 0 ; i < Nprongs ; ++i) {
511+ for (int i = 0 ; i < 3 ; ++i) {
513512 double deta, dphiAvg, dphi_AT_PV, dphi_AT_SpecificRadii, daughterEta, daughterPhi;
514513 bool sameCharge = false ;
515- daughterEta = -999 .;
516- daughterPhi = -999 .;
517-
518- switch (i) {
519- case Prong0:
520- daughterEta = part2.prong0Eta ();
521- daughterPhi = part2.prong0Phi ();
522- deta = part1.eta () - daughterEta;
523- dphi_AT_PV = part1.phi () - daughterPhi;
524- dphi_AT_SpecificRadii = PhiAtSpecificRadiiTPC (part1, radiiTPC) - PhiAtSpecificRadiiTPC<true , 0 >(part2, radiiTPC);
525- dphiAvg = AveragePhiStar<true >(part1, part2, 0 , &sameCharge);
526- // histdetadpi[0][0]->Fill(deta, dphiAvg);
527- break ;
528- case Prong1:
529- daughterEta = part2.prong1Eta ();
530- daughterPhi = part2.prong1Phi ();
531- deta = part1.eta () - daughterEta;
532- dphi_AT_PV = part1.phi () - daughterPhi;
533- dphi_AT_SpecificRadii = PhiAtSpecificRadiiTPC (part1, radiiTPC) - PhiAtSpecificRadiiTPC<true , 1 >(part2, radiiTPC);
534- dphiAvg = AveragePhiStar<true >(part1, part2, 1 , &sameCharge);
535- // histdetadpi[1][0]->Fill(deta, dphiAvg);
536- break ;
537- case Prong2:
538- daughterEta = part2.prong2Eta ();
539- daughterPhi = part2.prong2Phi ();
540- deta = part1.eta () - daughterEta;
541- dphi_AT_PV = part1.phi () - daughterPhi;
542- dphi_AT_SpecificRadii = PhiAtSpecificRadiiTPC (part1, radiiTPC) - PhiAtSpecificRadiiTPC<true , 2 >(part2, radiiTPC);
543- dphiAvg = AveragePhiStar<true >(part1, part2, 2 , &sameCharge);
544- // histdetadpi[2][0]->Fill(deta, dphiAvg);
545- break ;
514+ daughterEta = part2.prong2Eta ();
515+ daughterPhi = part2.prong2Phi ();
516+ deta = part1.eta () - daughterEta;
517+ dphi_AT_PV = part1.phi () - daughterPhi;
518+ dphi_AT_SpecificRadii = PhiAtSpecificRadiiTPC (part1, radiiTPC) - PhiAtSpecificRadiiTPC<true , 2 >(part2, radiiTPC);
519+ dphiAvg = AveragePhiStar<true >(part1, part2, 2 , &sameCharge);
520+
521+ if (Q3 == 999 ) {
522+ histdetadpi[i][0 ]->Fill (deta, dphiAvg);
523+ histdetadpi[i][2 ]->Fill (deta, dphi_AT_PV);
524+ if (fillQA) {
525+ histdetadpi_eta[i]->Fill (deta, dphiAvg, part1.eta (), daughterEta);
526+ histdetadpi_phi[i]->Fill (deta, dphiAvg, part1.phi (), daughterPhi);
527+ }
528+ } else if (Q3 < upperQ3LimitForPlotting) {
529+ histdetadpi[i][0 ]->Fill (deta, dphiAvg);
530+ histdetadpi[i][2 ]->Fill (deta, dphi_AT_PV);
531+ if (fillQA) {
532+ histdetadpi_eta[i]->Fill (deta, dphiAvg, part1.eta (), daughterEta);
533+ histdetadpi_phi[i]->Fill (deta, dphiAvg, part1.phi (), daughterPhi);
534+ }
546535 }
536+
537+ if (atWhichRadiiToSelect == 1 ) {
538+ if (std::pow (dphiAvg, 2 ) / std::pow (deltaPhiMax, 2 ) + std::pow (deta, 2 ) / std::pow (deltaEtaMax, 2 ) < 1 .) {
539+ pass = true ;
540+ } else {
541+ if (Q3 == 999 ) {
542+ histdetadpi[i][1 ]->Fill (deta, dphiAvg);
543+ histdetadpi[i][3 ]->Fill (deta, dphi_AT_PV);
544+ } else if (Q3 < upperQ3LimitForPlotting) {
545+ histdetadpi[i][1 ]->Fill (deta, dphiAvg);
546+ histdetadpi[i][3 ]->Fill (deta, dphi_AT_PV);
547+ }
548+ }
549+ } else if (atWhichRadiiToSelect == 0 ) {
550+ if (std::pow (dphi_AT_PV, 2 ) / std::pow (deltaPhiMax, 2 ) + std::pow (deta, 2 ) / std::pow (deltaEtaMax, 2 ) < 1 .) {
551+ pass = true ;
552+ } else {
553+ if (Q3 == 999 ) {
554+ histdetadpi[i][1 ]->Fill (deta, dphiAvg);
555+ histdetadpi[i][3 ]->Fill (deta, dphi_AT_PV);
556+ } else if (Q3 < upperQ3LimitForPlotting) {
557+ histdetadpi[i][1 ]->Fill (deta, dphiAvg);
558+ histdetadpi[i][3 ]->Fill (deta, dphi_AT_PV);
559+ }
560+ }
561+ } else if (atWhichRadiiToSelect == 2 ) {
562+ if (std::pow (dphi_AT_SpecificRadii, 2 ) / std::pow (deltaPhiMax, 2 ) + std::pow (deta, 2 ) / std::pow (deltaEtaMax, 2 ) < 1 .) {
563+ pass = true ;
564+ } else {
565+ if (Q3 == 999 ) {
566+ histdetadpi[i][1 ]->Fill (deta, dphiAvg);
567+ histdetadpi[i][3 ]->Fill (deta, dphi_AT_PV);
568+ } else if (Q3 < upperQ3LimitForPlotting) {
569+ histdetadpi[i][1 ]->Fill (deta, dphiAvg);
570+ histdetadpi[i][3 ]->Fill (deta, dphi_AT_PV);
571+ }
572+ }
573+ }
574+ }
575+
576+ return pass;
577+
578+ } else if constexpr (mPartOneType == o2::aod::femtodreamparticle::ParticleType::kTrack && mPartTwoType == o2::aod::femtodreamparticle::ParticleType::kCharmHadron2Prong ) {
579+ // check if provided particles are in agreement with the class instantiation
580+ if (!part2.candidateSelFlag ()) {
581+ LOG (fatal) << " FemtoDreamDetaDphiStar: passed arguments don't agree with FemtoDreamDetaDphiStar instantiation! Please provide Charm Hadron candidates."
582+ return false ;
583+ }
584+
585+ bool pass = false ;
586+
587+ for (int i = 0 ; i < 2 ; ++i) {
588+ double deta, dphiAvg, dphi_AT_PV, dphi_AT_SpecificRadii, daughterEta, daughterPhi;
589+ bool sameCharge = false ;
590+ daughterEta = part2.prong1Eta ();
591+ daughterPhi = part2.prong1Phi ();
592+ deta = part1.eta () - daughterEta;
593+ dphi_AT_PV = part1.phi () - daughterPhi;
594+ dphi_AT_SpecificRadii = PhiAtSpecificRadiiTPC (part1, radiiTPC) - PhiAtSpecificRadiiTPC<true , 1 >(part2, radiiTPC);
595+ dphiAvg = AveragePhiStar<true >(part1, part2, 1 , &sameCharge);
596+ // histdetadpi[1][0]->Fill(deta, dphiAvg);
597+
547598 if (Q3 == 999 ) {
548599 histdetadpi[i][0 ]->Fill (deta, dphiAvg);
549600 histdetadpi[i][2 ]->Fill (deta, dphi_AT_PV);
@@ -780,27 +831,25 @@ class FemtoDreamDetaDphiStar
780831 float PhiAtSpecificRadiiTPC (const T& part, float radii)
781832 {
782833 int charge = 0 ;
783- float phi0, pt;
834+ float phi0 = 0 .f ;
835+ float pt = 0 .f ;
836+
784837 if constexpr (isHF) {
785- switch (prong) {
786- case Prong0:
787- charge = part.charge (); // charge calculation according to 3-prong decay, Lc^+ --> P^+ + K^- + pi^+
788- phi0 = part.prong0Phi ();
789- pt = part.prong0Pt ();
790- break ;
791- case Prong1:
792- charge = -part.charge ();
793- phi0 = part.prong1Phi ();
794- pt = part.prong1Pt ();
795- break ;
796- case Prong2:
797- charge = part.charge ();
798- phi0 = part.prong2Phi ();
799- pt = part.prong2Pt ();
800- break ;
801- default :
802- // Handle invalid prong value if necessary
803- break ;
838+ static_assert (prong == 0 || prong == 1 || prong == 2 ,
839+ " PhiAtSpecificRadiiTPC: invalid prong index for HF"
840+
841+ if constexpr (prong == 0 ) {
842+ charge = part.charge ();
843+ phi0 = part.prong0Phi ();
844+ pt = part.prong0Pt ();
845+ } else if constexpr (prong == 1 ) {
846+ charge = -part.charge ();
847+ phi0 = part.prong1Phi ();
848+ pt = part.prong1Pt ();
849+ } else if constexpr (prong == 2 ) {
850+ charge = part.charge ();
851+ phi0 = part.prong2Phi ();
852+ pt = part.prong2Pt ();
804853 }
805854 } else {
806855 phi0 = part.phi ();
@@ -838,42 +887,42 @@ class FemtoDreamDetaDphiStar
838887 int PhiAtRadiiTPCForHF (const T& part, std::vector<float >& tmpVec, int prong)
839888 {
840889 int charge = 0 ;
841- float pt = -999 .;
842- float phi0 = -999 .;
843- switch (prong) {
844- case Prong0:
845- pt = part.prong0Pt ();
846- phi0 = part.prong0Phi ();
890+ if constexpr (mPartTwoType == o2::aod::femtodreamparticle::kCharmHadron3Prong ) {
891+ if (prong == 0 || prong == 2 ) {
847892 charge = part.charge ();
848- break ;
849- case Prong1:
850- pt = part.prong1Pt ();
851- phi0 = part.prong1Phi ();
893+ } else if (prong == 1 ) {
852894 charge = -part.charge ();
853- break ;
854- case Prong2:
855- pt = part.prong2Pt ();
856- phi0 = part.prong2Phi ();
895+ } else {
896+ return 0 ;
897+ }
898+ for (size_t i = 0 ; i < 9 ; ++i) {
899+ if (prong == 0 ) {
900+ tmpVec.push_back (PhiAtSpecificRadiiTPC<true , 0 >(part, tmpRadiiTPC[i]));
901+ } else if (prong == 1 ) {
902+ tmpVec.push_back (PhiAtSpecificRadiiTPC<true , 1 >(part, tmpRadiiTPC[i]));
903+ } else { // prong == 2
904+ tmpVec.push_back (PhiAtSpecificRadiiTPC<true , 2 >(part, tmpRadiiTPC[i]));
905+ }
906+ }
907+
908+ } else {
909+ if (prong == 0 ) {
857910 charge = part.charge ();
858- break ;
859- default :
860- // Handle invalid prong value
861- break ;
862- }
863- for (size_t i = 0 ; i < 9 ; i++) {
864- if (runOldVersion) {
865- tmpVec.push_back (phi0 - std::asin (0.3 * charge * 0.1 * magfield * tmpRadiiTPC[i] * 0.01 / (2 . * pt)));
911+ } else if (prong == 1 ) {
912+ charge = -part.charge ();
913+ } else {
914+ return 0 ;
866915 }
867- if (!runOldVersion) {
868- auto arg = 0.3 * charge * magfield * tmpRadiiTPC[i] * 0.01 / (2 . * pt);
869- // for very low pT particles, this value goes outside of range -1 to 1 at at large tpc radius; asin fails
870- if (std::fabs (arg) < 1 ) {
871- tmpVec.push_back (phi0 - std::asin (0.3 * charge * magfield * tmpRadiiTPC[i] * 0.01 / (2 . * pt)));
872- } else {
873- tmpVec.push_back (999 );
916+
917+ for (size_t i = 0 ; i < 9 ; ++i) {
918+ if (prong == 0 ) {
919+ tmpVec.push_back (PhiAtSpecificRadiiTPC<true , 0 >(part, tmpRadiiTPC[i]));
920+ } else { // prong == 1
921+ tmpVec.push_back (PhiAtSpecificRadiiTPC<true , 1 >(part, tmpRadiiTPC[i]));
874922 }
875923 }
876924 }
925+
877926 return charge;
878927 }
879928
0 commit comments