From 213adb92b959bdb359c0403fad7ab647f4de2b07 Mon Sep 17 00:00:00 2001 From: fbellini Date: Wed, 19 Feb 2014 23:07:19 +0100 Subject: [PATCH] Added DCA of daughters and DCA product computation (Massimo) --- PWGLF/RESONANCES/AliRsnCutMiniPair.cxx | 3 ++ PWGLF/RESONANCES/AliRsnCutMiniPair.h | 1 + PWGLF/RESONANCES/AliRsnCutTrackQuality.cxx | 52 ++++++++++++++++++---- PWGLF/RESONANCES/AliRsnCutTrackQuality.h | 15 ++++--- PWGLF/RESONANCES/AliRsnLoopPair.cxx | 12 +++++ PWGLF/RESONANCES/AliRsnLoopPair.h | 4 +- PWGLF/RESONANCES/AliRsnMiniPair.cxx | 27 ++++++++++- PWGLF/RESONANCES/AliRsnMiniPair.h | 29 ++++++++---- PWGLF/RESONANCES/AliRsnMiniParticle.cxx | 22 ++++++++- PWGLF/RESONANCES/AliRsnMiniParticle.h | 6 ++- PWGLF/RESONANCES/AliRsnMiniValue.cxx | 9 ++++ PWGLF/RESONANCES/AliRsnMiniValue.h | 3 ++ PWGLF/RESONANCES/AliRsnMother.cxx | 37 ++++++++++++++- PWGLF/RESONANCES/AliRsnMother.h | 6 ++- PWGLF/RESONANCES/AliRsnValuePair.cxx | 8 ++++ PWGLF/RESONANCES/AliRsnValuePair.h | 1 + 16 files changed, 203 insertions(+), 32 deletions(-) diff --git a/PWGLF/RESONANCES/AliRsnCutMiniPair.cxx b/PWGLF/RESONANCES/AliRsnCutMiniPair.cxx index bdf06c328c1..dbf1b41ce08 100644 --- a/PWGLF/RESONANCES/AliRsnCutMiniPair.cxx +++ b/PWGLF/RESONANCES/AliRsnCutMiniPair.cxx @@ -46,6 +46,9 @@ Bool_t AliRsnCutMiniPair::IsSelected(TObject *obj) case kMomentumComparison: AliWarning("TODO: implement this"); return kTRUE; + case kDCAproduct: + fCutValueD = pair->DCAProduct(); + return OkRangeD(); default: AliWarning("Undefined enum value"); return kTRUE; diff --git a/PWGLF/RESONANCES/AliRsnCutMiniPair.h b/PWGLF/RESONANCES/AliRsnCutMiniPair.h index 2f27c413d91..ba7fbb2faca 100644 --- a/PWGLF/RESONANCES/AliRsnCutMiniPair.h +++ b/PWGLF/RESONANCES/AliRsnCutMiniPair.h @@ -19,6 +19,7 @@ public: kRapidityRange, kRapidityRangeMC, kMomentumComparison, + kDCAproduct, kTypes }; diff --git a/PWGLF/RESONANCES/AliRsnCutTrackQuality.cxx b/PWGLF/RESONANCES/AliRsnCutTrackQuality.cxx index 2b71792d2b5..379e22b2915 100644 --- a/PWGLF/RESONANCES/AliRsnCutTrackQuality.cxx +++ b/PWGLF/RESONANCES/AliRsnCutTrackQuality.cxx @@ -35,9 +35,12 @@ AliRsnCutTrackQuality::AliRsnCutTrackQuality(const char *name) : fFlagsOn(0x0), fFlagsOff(0x0), fRejectKinkDaughters(kTRUE), - fDCARfixed(kTRUE), + fDCARmaxfixed(kTRUE), + fDCARminfixed(kTRUE), fDCARptFormula(""), + fDCARptFormulaMin(""), fDCARmax(1E20), + fDCARmin(0), fDCAZfixed(kTRUE), fDCAZptFormula(""), fDCAZmax(1E20), @@ -68,9 +71,12 @@ AliRsnCutTrackQuality::AliRsnCutTrackQuality(const AliRsnCutTrackQuality ©) fFlagsOn(copy.fFlagsOn), fFlagsOff(copy.fFlagsOff), fRejectKinkDaughters(copy.fRejectKinkDaughters), - fDCARfixed(copy.fDCARfixed), + fDCARmaxfixed(copy.fDCARmaxfixed), + fDCARminfixed(copy.fDCARminfixed), fDCARptFormula(copy.fDCARptFormula), + fDCARptFormulaMin(copy.fDCARptFormulaMin), fDCARmax(copy.fDCARmax), + fDCARmin(copy.fDCARmin), fDCAZfixed(copy.fDCAZfixed), fDCAZptFormula(copy.fDCAZptFormula), fDCAZmax(copy.fDCAZmax), @@ -110,9 +116,12 @@ AliRsnCutTrackQuality &AliRsnCutTrackQuality::operator=(const AliRsnCutTrackQual fFlagsOn = copy.fFlagsOn; fFlagsOff = copy.fFlagsOff; fRejectKinkDaughters = copy.fRejectKinkDaughters; - fDCARfixed = copy.fDCARfixed; + fDCARmaxfixed = copy.fDCARmaxfixed; + fDCARminfixed = copy.fDCARminfixed; fDCARptFormula = copy.fDCARptFormula; + fDCARptFormulaMin = copy.fDCARptFormulaMin; fDCARmax = copy.fDCARmax; + fDCARmin = copy.fDCARmin; fDCAZfixed = copy.fDCAZfixed; fDCAZptFormula = copy.fDCAZptFormula; fDCAZmax = copy.fDCAZmax; @@ -145,9 +154,12 @@ void AliRsnCutTrackQuality::DisableAll() fFlagsOn = 0x0; fFlagsOff = 0x0; fRejectKinkDaughters = kFALSE; - fDCARfixed = kTRUE; + fDCARmaxfixed = kTRUE; + fDCARminfixed = kTRUE; fDCARptFormula = ""; + fDCARptFormulaMin = ""; fDCARmax = 1E20; + fDCARmin = 0; fDCAZfixed = kTRUE; fDCAZptFormula = ""; fDCAZmax = 1E20; @@ -244,10 +256,15 @@ Bool_t AliRsnCutTrackQuality::CheckESD(AliESDtrack *track) cuts.SetEtaRange(fEta[0], fEta[1]); // transverse DCA cuts - if (fDCARfixed) + if (fDCARmaxfixed) cuts.SetMaxDCAToVertexXY(fDCARmax); else cuts.SetMaxDCAToVertexXYPtDep(fDCARptFormula.Data()); + + if (fDCARminfixed) + cuts.SetMinDCAToVertexXY(fDCARmin); + else + cuts.SetMinDCAToVertexXYPtDep(fDCARptFormulaMin.Data()); // longitudinal DCA cuts if (fDCAZfixed) @@ -390,17 +407,28 @@ Bool_t AliRsnCutTrackQuality::CheckAOD(AliAODTrack *track) return kFALSE; } // if the DCA cut is not fixed, compute current value - if (!fDCARfixed) { + if (!fDCARmaxfixed) { TString str(fDCARptFormula); str.ReplaceAll("pt", "x"); TFormula dcaXY(Form("%s_dcaXY", GetName()), str.Data()); fDCARmax = dcaXY.Eval(track->Pt()); } + if (!fDCARminfixed) { + TString str2(fDCARptFormulaMin); + str2.ReplaceAll("pt", "x"); + TFormula dcaXY_2(Form("%s_dcaXY_2", GetName()), str2.Data()); + fDCARmin = dcaXY_2.Eval(track->Pt()); + } // check the cut if (TMath::Abs(b[0]) > fDCARmax) { AliDebug(AliLog::kDebug + 2, "Too large transverse DCA"); return kFALSE; } + + if (TMath::Abs(b[0]) < fDCARmin) { + AliDebug(AliLog::kDebug + 2, "Too short transverse DCA"); + return kFALSE; + } // step #5: DCA cut (longitudinal) // the DCA has already been computed above @@ -449,10 +477,16 @@ void AliRsnCutTrackQuality::Print(const Option_t *) const AliInfo(Form("ITS requirements : min. cluster = %d (all), %d (SPD), max chi2 = %f", fITSminNClusters, fSPDminNClusters, fITSmaxChi2)); - if (fDCARfixed) { - AliInfo(Form("DCA r cut : fixed to %f cm", fDCARmax)); + if (fDCARmaxfixed) { + AliInfo(Form("Max DCA r cut : fixed to %f cm", fDCARmax)); + } else { + AliInfo(Form("Max DCA r cut formula : %s", fDCARptFormula.Data())); + } + + if (fDCARminfixed) { + AliInfo(Form("Min DCA r cut : fixed to %f cm", fDCARmin)); } else { - AliInfo(Form("DCA r cut formula : %s", fDCARptFormula.Data())); + AliInfo(Form("Min DCA r cut formula : %s", fDCARptFormulaMin.Data())); } if (fDCAZfixed) { diff --git a/PWGLF/RESONANCES/AliRsnCutTrackQuality.h b/PWGLF/RESONANCES/AliRsnCutTrackQuality.h index cb27bd1901d..78a4fc539dc 100644 --- a/PWGLF/RESONANCES/AliRsnCutTrackQuality.h +++ b/PWGLF/RESONANCES/AliRsnCutTrackQuality.h @@ -36,8 +36,10 @@ public: void SetPtRange(Double_t a, Double_t b) {fPt[0] = TMath::Min(a, b); fPt[1] = TMath::Max(a, b);} void SetEtaRange(Double_t a, Double_t b) {fEta[0] = TMath::Min(a, b); fEta[1] = TMath::Max(a, b);} - void SetDCARPtFormula(const char *formula) {fDCARptFormula = formula; fDCARfixed = kFALSE;} - void SetDCARmax(Double_t value) {fDCARmax = value; fDCARptFormula = ""; fDCARfixed = kTRUE;} + void SetDCARPtFormula(const char *formula) {fDCARptFormula = formula; fDCARmaxfixed = kFALSE;} + void SetDCARPtFormulaMin(const char *formula) {fDCARptFormulaMin = formula; fDCARminfixed = kFALSE;} + void SetDCARmax(Double_t value) {fDCARmax = value; fDCARptFormula = ""; fDCARmaxfixed = kTRUE;} + void SetDCARmin(Double_t value) {fDCARmin = value; fDCARptFormulaMin = ""; fDCARminfixed = kTRUE;} void SetDCAZPtFormula(const char *formula) {fDCAZptFormula = formula; fDCAZfixed = kFALSE;} void SetDCAZmax(Double_t value) {fDCAZmax = value; fDCAZptFormula = ""; fDCAZfixed = kTRUE;} @@ -79,9 +81,12 @@ protected: Double_t fEta[2]; // eta range Bool_t fRejectKinkDaughters; // switch to kTRUE if daughters of kinks must be rejected - Bool_t fDCARfixed; // flag to switch between fixed and pt-dependent DCA cut - TString fDCARptFormula; // expression to compute transverse DCA sigma w.r. to pt + Bool_t fDCARmaxfixed; // flag to switch between fixed and pt-dependent DCA cut (maximum) + Bool_t fDCARminfixed; // flag to switch between fixed and pt-dependent DCA cut (minimum) + TString fDCARptFormula; // expression to compute transverse DCA sigma w.r. to pt (maximum) + TString fDCARptFormulaMin; // expression to compute transverse DCA sigma w.r. to pt (minimum) Double_t fDCARmax; // maximum value for transverse DCA + Double_t fDCARmin; // minimum value for transverse DCA Bool_t fDCAZfixed; // flag to switch between fixed and pt-dependent DCA cut TString fDCAZptFormula; // expression to compute longitudinal DCA sigma w.r. to pt @@ -105,6 +110,6 @@ protected: Bool_t fCheckOnlyFilterBit; // check only the filter bit AliESDtrackCuts *fESDtrackCuts; // pointer to AliESDtrackCuts object - ClassDef(AliRsnCutTrackQuality, 3) + ClassDef(AliRsnCutTrackQuality, 4) }; #endif diff --git a/PWGLF/RESONANCES/AliRsnLoopPair.cxx b/PWGLF/RESONANCES/AliRsnLoopPair.cxx index 9fe239eda1e..a5d3d1fcdc8 100644 --- a/PWGLF/RESONANCES/AliRsnLoopPair.cxx +++ b/PWGLF/RESONANCES/AliRsnLoopPair.cxx @@ -36,6 +36,7 @@ AliRsnLoopPair::AliRsnLoopPair(const char *name, AliRsnPairDef *def, Bool_t isMi fUseMCRef(kFALSE), fCheckDecay(kFALSE), fRangeY(1E20), + fRangeDCAproduct(1E20), fPairDef(def), fPairCuts(0x0), fMother() @@ -56,6 +57,7 @@ AliRsnLoopPair::AliRsnLoopPair(const AliRsnLoopPair ©) : fUseMCRef(copy.fUseMCRef), fCheckDecay(copy.fCheckDecay), fRangeY(copy.fRangeY), + fRangeDCAproduct(copy.fRangeDCAproduct), fPairDef(copy.fPairDef), fPairCuts(copy.fPairCuts), fMother(copy.fMother) @@ -83,6 +85,7 @@ AliRsnLoopPair &AliRsnLoopPair::operator=(const AliRsnLoopPair ©) fUseMCRef = copy.fUseMCRef; fCheckDecay = copy.fCheckDecay; fRangeY = copy.fRangeY; + fRangeDCAproduct = copy.fRangeDCAproduct; fPairDef = copy.fPairDef; fPairCuts = copy.fPairCuts; fMother = copy.fMother; @@ -205,6 +208,11 @@ Int_t AliRsnLoopPair::DoLoop AliDebugClass(2, Form("[%s]: Outside rapidity range", GetName())); continue; } + // check daughter's DCA product range + if (TMath::Abs(fMother.DCAproduct()) > fRangeDCAproduct) { + AliDebugClass(2, Form("[%s]: Outside daughter's DCA products range", GetName())); + continue; + } // check mother if (fOnlyTrue) { if (!IsTrueMother()) { @@ -500,6 +508,10 @@ Int_t AliRsnLoopPair::LoopTrueMC(AliRsnEvent *rsn) AliDebugClass(2, Form("[%s]: Outside rapidity range", GetName())); continue; } + if (TMath::Abs(fMother.DCAproduct()) > fRangeDCAproduct) { + AliDebugClass(2, Form("[%s]: Outside rapidity range", GetName())); + continue; + } // fill outputs next.Reset(); while ( (out = (AliRsnListOutput *)next()) ) { diff --git a/PWGLF/RESONANCES/AliRsnLoopPair.h b/PWGLF/RESONANCES/AliRsnLoopPair.h index b8c1843f7a4..74a97316018 100644 --- a/PWGLF/RESONANCES/AliRsnLoopPair.h +++ b/PWGLF/RESONANCES/AliRsnLoopPair.h @@ -37,6 +37,7 @@ public: void SetCheckDecay(Bool_t check = kTRUE) {fCheckDecay = check;} void SetListID(Int_t i, Int_t val) {if (i==0||i==1) fListID[i] = val;} void SetRangeY(Double_t range) {fRangeY = range;} + void SetRangeDCAproduct(Double_t range) {fRangeDCAproduct = range;} // methods Bool_t IsTrueMother(); @@ -57,6 +58,7 @@ protected: Bool_t fCheckDecay; // is the decay channel correct in a true pair? Int_t fListID[2]; // indexes of the two entry lists to be used Double_t fRangeY; // range in rapidity (serves always) + Double_t fRangeDCAproduct; // range in the daughters's DCA product AliRsnPairDef *fPairDef; // pair definition AliRsnCutSet *fPairCuts; // collection of all cuts @@ -65,7 +67,7 @@ protected: private: - ClassDef(AliRsnLoopPair, 4) + ClassDef(AliRsnLoopPair, 5) }; #endif diff --git a/PWGLF/RESONANCES/AliRsnMiniPair.cxx b/PWGLF/RESONANCES/AliRsnMiniPair.cxx index cef3a02e36a..ad6e49ee5af 100644 --- a/PWGLF/RESONANCES/AliRsnMiniPair.cxx +++ b/PWGLF/RESONANCES/AliRsnMiniPair.cxx @@ -11,11 +11,13 @@ void AliRsnMiniPair::Fill // Fill this object with data coming // from arguments // - p1->Set4Vector(fP1[0], m1, kFALSE); p2->Set4Vector(fP2[0], m2, kFALSE); p1->Set4Vector(fP1[1], m1, kTRUE ); p2->Set4Vector(fP2[1], m2, kTRUE ); + + fDCA1 = p1->DCA(); + fDCA2 = p2->DCA(); fMother = -1; if (p1->Mother() == p2->Mother()) { @@ -160,6 +162,29 @@ Double_t AliRsnMiniPair::DaughterPt(Int_t daughterId, Bool_t mc) return fP2[ID(mc)].Pt(); } +//__________________________________________________________________________________________________ +Double_t AliRsnMiniPair::DaughterDCA(Int_t daughterId) +{ + // + //returns dca to Primary Vertex of the daughter + // + + if (daughterId==0) + return fDCA1; + else + return fDCA2; +} + +//__________________________________________________________________________________________________ +Double_t AliRsnMiniPair::DCAProduct() +{ + // + //returns products of the DCA of the 2 daughters + // + + return fDCA1*fDCA2; +} + //__________________________________________________________________________________________________ void AliRsnMiniPair::DaughterPxPyPz(Int_t daughterId, Bool_t mc, Double_t *pxpypz) { diff --git a/PWGLF/RESONANCES/AliRsnMiniPair.h b/PWGLF/RESONANCES/AliRsnMiniPair.h index e64a9494778..6fcfcbccd7f 100644 --- a/PWGLF/RESONANCES/AliRsnMiniPair.h +++ b/PWGLF/RESONANCES/AliRsnMiniPair.h @@ -9,21 +9,24 @@ #include #include +#include "TClonesArray.h" +#include "AliRsnListOutput.h" +class AliRsnListOutput; class AliRsnMiniParticle; class AliRsnMiniPair : public TObject { public: - AliRsnMiniPair() : fMother(-1), fMotherPDG(0) { } - + AliRsnMiniPair() : fDCA1(0), fDCA2(0), fMother(-1), fMotherPDG(0), fNsisters(-1) { } + Int_t &Mother() {return fMother;} Int_t &MotherPDG() {return fMotherPDG;} - void Fill(AliRsnMiniParticle *p1, AliRsnMiniParticle *p2, Double_t m1, Double_t m2, Double_t refMass); - void FillRef(Double_t mass); - void InvertP(Bool_t first); + void Fill(AliRsnMiniParticle *p1, AliRsnMiniParticle *p2, Double_t m1, Double_t m2, Double_t refMass); + void FillRef(Double_t mass); + void InvertP(Bool_t first); - Int_t ID(Bool_t mc) const {if (mc) return 1; else return 0;} + Int_t ID(Bool_t mc) const {if (mc) return 1; else return 0;} TLorentzVector &P1 (Bool_t mc) {return fP1 [ID(mc)];} TLorentzVector &P2 (Bool_t mc) {return fP2 [ID(mc)];} @@ -41,9 +44,13 @@ public: Double_t PtRatio(Bool_t mc) const; Double_t DipAngle(Bool_t mc) const; Double_t CosThetaStar(Bool_t mc); - Double_t DaughterPt(Int_t daughterId, Bool_t mc); + Double_t DaughterPt(Int_t daughterId, Bool_t mc); + Double_t DaughterDCA(Int_t daughterId); + Double_t DCAProduct(); void DaughterPxPyPz(Int_t daughterId, Bool_t mc, Double_t *pxpypz); - + Short_t Nsisters() {return fNsisters;} + void SetNsisters(Short_t value) {fNsisters=value;} + private: TLorentzVector fP1 [2]; // 1st daughter momentum @@ -51,10 +58,14 @@ public: TLorentzVector fSum[2]; // sum of momenta TLorentzVector fRef[2]; // same as 'fSum' but with nominal resonance mass + Double_t fDCA1; // 1st daughter DCA + Double_t fDCA2; // 2nd daughter DCA + Int_t fMother; // label of mothers (when common) Int_t fMotherPDG; // PDG code of mother (when common) + Short_t fNsisters; // total number of mother's daughters in the MC stack - ClassDef(AliRsnMiniPair,1) + ClassDef(AliRsnMiniPair,2) }; #endif diff --git a/PWGLF/RESONANCES/AliRsnMiniParticle.cxx b/PWGLF/RESONANCES/AliRsnMiniParticle.cxx index 045f3cf196b..90342d0ca47 100644 --- a/PWGLF/RESONANCES/AliRsnMiniParticle.cxx +++ b/PWGLF/RESONANCES/AliRsnMiniParticle.cxx @@ -9,7 +9,8 @@ #include #include -#include "AliRsnDaughter.h" +#include "AliAODEvent.h" +#include "AliRsnEvent.h" #include "AliRsnMiniParticle.h" ClassImp(AliRsnMiniParticle) @@ -22,6 +23,7 @@ void AliRsnMiniParticle::CopyDaughter(AliRsnDaughter *daughter) // // reset what could not be initialized + fDCA = 0.0; //typically used for D0 analysis fPDG = 0; fMother = -1; fMotherPDG = 0; @@ -52,9 +54,25 @@ void AliRsnMiniParticle::CopyDaughter(AliRsnDaughter *daughter) fMother = daughter->GetMother(); fMotherPDG = daughter->GetMotherPDG(); } + + AliRsnEvent *event = (AliRsnEvent *) daughter->GetOwnerEvent(); + if (event && event->IsAOD()){ + AliAODTrack *track = (AliAODTrack*) daughter->Ref2AODtrack(); + AliAODEvent *aodEvent = (AliAODEvent*) event->GetRefAOD(); + if (track && aodEvent) { + AliVVertex *vertex = (AliVVertex*) aodEvent->GetPrimaryVertex(); + Double_t b[2], cov[3]; + if (vertex) { + track->PropagateToDCA(vertex, aodEvent->GetMagneticField(), kVeryBig, b, cov); + fDCA = b[0]; + } + } + } else { + AliWarning("DCA not implemented for ESDs"); + } + } - //__________________________________________________________________________________________________ Double_t AliRsnMiniParticle::Mass() { diff --git a/PWGLF/RESONANCES/AliRsnMiniParticle.h b/PWGLF/RESONANCES/AliRsnMiniParticle.h index d5abb51dd86..8cedf26ebb5 100644 --- a/PWGLF/RESONANCES/AliRsnMiniParticle.h +++ b/PWGLF/RESONANCES/AliRsnMiniParticle.h @@ -16,7 +16,7 @@ class AliRsnDaughter; class AliRsnMiniParticle : public TObject { public: - AliRsnMiniParticle() : fIndex(-1), fCharge(0), fPDG(0), fMother(0), fMotherPDG(0), fCutBits(0x0) {Int_t i = 3; while (i--) fPsim[i] = fPrec[i] = 0.0;} + AliRsnMiniParticle() : fIndex(-1), fCharge(0), fPDG(0), fMother(0), fMotherPDG(0), fDCA(0), fCutBits(0x0) {Int_t i = 3; while (i--) fPsim[i] = fPrec[i] = 0.0;} Int_t &Index() {return fIndex;} Char_t &Charge() {return fCharge;} @@ -35,6 +35,7 @@ public: Int_t &Mother() {return fMother;} Short_t &MotherPDG() {return fMotherPDG;} UShort_t &CutBits() {return fCutBits;} + Double_t DCA() {return fDCA;} Bool_t HasCutBit(Int_t i) {UShort_t bit = 1 << i; return ((fCutBits & bit) != 0);} void SetCutBit(Int_t i) {UShort_t bit = 1 << i; fCutBits |= bit;} void ClearCutBit(Int_t i) {UShort_t bit = 1 << i; fCutBits &= (~bit);} @@ -51,9 +52,10 @@ private: Short_t fPDG; // particle PDG code Int_t fMother; // index of mother in its container Short_t fMotherPDG; // PDG code of mother + Double_t fDCA; // DCA of the particle UShort_t fCutBits; // list of bits used to know what cuts were passed by this track - ClassDef(AliRsnMiniParticle,2) + ClassDef(AliRsnMiniParticle,3) }; #endif diff --git a/PWGLF/RESONANCES/AliRsnMiniValue.cxx b/PWGLF/RESONANCES/AliRsnMiniValue.cxx index c3db0685ccd..a6d327545a9 100644 --- a/PWGLF/RESONANCES/AliRsnMiniValue.cxx +++ b/PWGLF/RESONANCES/AliRsnMiniValue.cxx @@ -116,6 +116,9 @@ const char *AliRsnMiniValue::TypeName(EType type) case kSecondDaughterPt: return "SecondDaughterPt"; case kFirstDaughterP: return "FirstDaughterP"; case kSecondDaughterP: return "SecondDaughterP"; + case kDCAproduct: return "DaughterDCAproduct"; + case kFirstDaughterDCA: return "FirstDaughterDCA"; + case kSecondDaughterDCA: return "SecondDaughterDCA"; default: return "Undefined"; } } @@ -200,6 +203,12 @@ Float_t AliRsnMiniValue::Eval(AliRsnMiniPair *pair, AliRsnMiniEvent *event) case kSecondDaughterP: pair->DaughterPxPyPz(1,fUseMCInfo, p3); return TMath::Sqrt(p3[0]*p3[0]+p3[1]*p3[1]+p3[2]*p3[2]); + case kDCAproduct: + return pair->DCAProduct(); + case kFirstDaughterDCA: + return pair->DaughterDCA(0); + case kSecondDaughterDCA: + return pair->DaughterDCA(1); default: AliError("Invalid value type"); return 1E20; diff --git a/PWGLF/RESONANCES/AliRsnMiniValue.h b/PWGLF/RESONANCES/AliRsnMiniValue.h index eea1160c2bc..218d9396456 100644 --- a/PWGLF/RESONANCES/AliRsnMiniValue.h +++ b/PWGLF/RESONANCES/AliRsnMiniValue.h @@ -38,6 +38,9 @@ public: kSecondDaughterPt, //pt of the second daughter of the pair kFirstDaughterP, //p of the first daughter of the pair kSecondDaughterP, //p of the second daughter of the pair + kDCAproduct, // product of the daughter's dca to PV (same in AliRsnValuePair) + kFirstDaughterDCA, //DCA to PV of the first daughter of the pair + kSecondDaughterDCA, //DCA to PV of the second daughter of the pair kTypes // -- general limit ---------------------------------------------------------- }; diff --git a/PWGLF/RESONANCES/AliRsnMother.cxx b/PWGLF/RESONANCES/AliRsnMother.cxx index 8f365f36814..f32b7b23e95 100644 --- a/PWGLF/RESONANCES/AliRsnMother.cxx +++ b/PWGLF/RESONANCES/AliRsnMother.cxx @@ -44,7 +44,8 @@ AliRsnMother::AliRsnMother(const AliRsnMother &obj) : fSum(obj.fSum), fSumMC(obj.fSumMC), fRef(obj.fRef), - fRefMC(obj.fRefMC) + fRefMC(obj.fRefMC), + fDCAproduct(obj.fDCAproduct) { // // Copy constructor. @@ -72,6 +73,7 @@ AliRsnMother &AliRsnMother::operator=(const AliRsnMother &obj) fRefEvent = obj.fRefEvent; fDaughter[0] = obj.fDaughter[0]; fDaughter[1] = obj.fDaughter[1]; + fDCAproduct = obj.fDCAproduct; return (*this); } @@ -225,6 +227,39 @@ Double_t AliRsnMother::CosThetaStar(Bool_t first, Bool_t useMC) return cosThetaStar; } +//__________________________________________________________________________________________________ +Double_t AliRsnMother::DCAproduct() +{ + // + // returns product of DCA of the two daughters + // + AliRsnEvent *event1 = fDaughter[0]->GetOwnerEvent(); + AliRsnEvent *event2 = fDaughter[1]->GetOwnerEvent(); + + if(event1 != event2){ + AliError("Attempting to build pair with tracks coming from different events"); + return 0.0; + } + + if (event1->IsAOD()) { + AliAODEvent *aodEvent = (AliAODEvent*)event1->GetRefAOD(); + if (!aodEvent) return 0.0; + AliAODTrack *track1 = (AliAODTrack*)fDaughter[0]->Ref2AODtrack(); + AliAODTrack *track2 = (AliAODTrack*)fDaughter[1]->Ref2AODtrack(); + AliVVertex *vertex = aodEvent->GetPrimaryVertex(); + if (!vertex || !track1 || track2) return 0.0; + + Double_t b1[2], cov1[3], b2[2], cov2[3]; + track1->PropagateToDCA(vertex, aodEvent->GetMagneticField(), kVeryBig, b1, cov1); + track2->PropagateToDCA(vertex, aodEvent->GetMagneticField(), kVeryBig, b2, cov2); + fDCAproduct = b1[0]*b2[0]; + } else { + AliWarning("DCA product not implemented for ESDs"); + fDCAproduct=0.0; + } + return fDCAproduct; +} + //__________________________________________________________________________________________________ void AliRsnMother::PrintInfo(const Option_t * /*option*/) const { diff --git a/PWGLF/RESONANCES/AliRsnMother.h b/PWGLF/RESONANCES/AliRsnMother.h index 834ef3fcb79..d2ebbe7a192 100644 --- a/PWGLF/RESONANCES/AliRsnMother.h +++ b/PWGLF/RESONANCES/AliRsnMother.h @@ -20,7 +20,7 @@ class AliRsnEvent; class AliRsnMother : public TObject { public: - AliRsnMother() : fRefEvent(0), fSum(), fSumMC(), fRef(), fRefMC() {fDaughter[0] = fDaughter[1] = 0;} + AliRsnMother() : fRefEvent(0), fSum(), fSumMC(), fRef(), fRefMC(), fDCAproduct(0) {fDaughter[0] = fDaughter[1] = 0;} AliRsnMother(const AliRsnMother &obj); AliRsnMother &operator=(const AliRsnMother &obj); virtual ~AliRsnMother(); @@ -37,6 +37,7 @@ public: TLorentzVector &Ref(Bool_t mc) {return (mc ? fRefMC : fRef);} Bool_t GetResolution(Double_t &value); Double_t Rapidity(Bool_t mc) {if (mc) return fRefMC.Rapidity(); else return fRef.Rapidity();} + Double_t DCAproduct(); // checks Bool_t IsLabelEqual() const {return TMath::Abs(fDaughter[0]->GetLabel()) == TMath::Abs(fDaughter[1]->GetLabel());} @@ -64,8 +65,9 @@ private: TLorentzVector fSumMC; // sum computed from the two daughters (sim) TLorentzVector fRef; // same to sum, but with fixed mass hypothesis (rec) TLorentzVector fRefMC; // same to sum, but with fixed mass hypothesis (sim) + Double_t fDCAproduct; // product of the daughter's DCA to Primary Vertex - ClassDef(AliRsnMother, 1) + ClassDef(AliRsnMother, 2) }; #endif diff --git a/PWGLF/RESONANCES/AliRsnValuePair.cxx b/PWGLF/RESONANCES/AliRsnValuePair.cxx index 313492aa8ee..3b418abfbda 100644 --- a/PWGLF/RESONANCES/AliRsnValuePair.cxx +++ b/PWGLF/RESONANCES/AliRsnValuePair.cxx @@ -52,6 +52,8 @@ #include "AliRsnPairDef.h" #include "AliRsnDaughterDef.h" +#include "AliRsnMiniPair.h" + #include "AliRsnValuePair.h" ClassImp(AliRsnValuePair) @@ -112,6 +114,7 @@ const char *AliRsnValuePair::GetTypeName() const case kDipAngle: return "PairDipAngle"; case kCosThetaStar: return "PairCosThetaStar"; case kAngleLeading: return "PairAngleToLeading"; + case kDCAproduct: return "DaughterDCAProduct"; default: return "Undefined"; } } @@ -136,6 +139,8 @@ Bool_t AliRsnValuePair::Eval(TObject *object) TLorentzVector &p1 = fMother->GetDaughter(0)->P(fUseMCInfo); TLorentzVector &p2 = fMother->GetDaughter(1)->P(fUseMCInfo); + + // utility variables Bool_t success; @@ -178,6 +183,9 @@ Bool_t AliRsnValuePair::Eval(TObject *object) case kAngleLeading: fComputedValue = fMother->AngleToLeading(success); return success; + case kDCAproduct: + fComputedValue = fMother->DCAproduct(); + return kTRUE; default: AliError(Form("[%s] Invalid value type for this computation", GetName())); return kFALSE; diff --git a/PWGLF/RESONANCES/AliRsnValuePair.h b/PWGLF/RESONANCES/AliRsnValuePair.h index cf15db59861..60e62dbbc78 100644 --- a/PWGLF/RESONANCES/AliRsnValuePair.h +++ b/PWGLF/RESONANCES/AliRsnValuePair.h @@ -27,6 +27,7 @@ public: kDipAngle, // inverse cosine of the angle between daughter vector momenta kCosThetaStar, // polarization angle kAngleLeading, // angle to leading particle + kDCAproduct, // product of the daughter's DCA kTypes }; -- 2.43.0