]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG2/RESONANCES/AliRsnMother.cxx
Protection against negative argument of log, add histograms to study shape/size of...
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnMother.cxx
index ba09a00dcf457745edd330f97fdc21c6b0131974..3700074e6c1bf64eafc15954baaa0a63d7905562 100644 (file)
@@ -1,83 +1,80 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+////////////////////////////////////////////////////////////////////////////////
 //
-// Class AliRsnMother
-//
-// Implementation of a pair of tracks, for several purposes
-// - computing the total 4-momentum & inv. mass for output histos filling
-// - evaluating cut checks on the pair of particles
-// - evaluating any kind of kinematic value over their sum
+//  This class implements a candidate resonance. It has two pointers to its
+//  two candidate daughters, whose 4-momenta are combined to obtain the mother
+//  invariant mass and other kinematical quantities.
+//  This class contains also some methods used to compute kinematical relations
+//  between the candidate resonance and other particles.
 //
-// authors: Martin Vala (martin.vala@cern.ch)
-//          Alberto Pulvirenti (alberto.pulvirenti@ct.infn.it)
+//  authors: A. Pulvirenti (alberto.pulvirenti@ct.infn.it)
+//           M. Vala (martin.vala@cern.ch)
 //
+////////////////////////////////////////////////////////////////////////////////
+
 #include <Riostream.h>
 #include <TVector3.h>
-#include "AliRsnDaughter.h"
-#include "AliRsnPairDef.h"
-#include "AliRsnMother.h"
 
-ClassImp(AliRsnMother)
+#include "AliAODMCParticle.h"
+#include "AliMCParticle.h"
+#include "AliRsnEvent.h"
 
-//_____________________________________________________________________________
-AliRsnMother::AliRsnMother() : 
-  fUseMC(kFALSE),
-  fDefaultMass(0.0),
-  fSum(),
-  fSumMC(),
-  fRef(),
-  fRefMC()
-{
-//
-// Constructor.
-// Initializes all variables to meaningless values.
-//
+#include "AliRsnMother.h"
 
-  Int_t i;
-  for (i = 0; i < 2; i++) fDaughter[i] = 0x0;
-}
+ClassImp(AliRsnMother)
 
-//_____________________________________________________________________________
-AliRsnMother::AliRsnMother(const AliRsnMother &obj) : 
-  TObject(obj), 
-  fUseMC(obj.fUseMC),
-  fDefaultMass(obj.fDefaultMass),
-  fSum(obj.fSum),
-  fSumMC(obj.fSumMC),
-  fRef(obj.fRef),
-  fRefMC(obj.fRefMC)
+//__________________________________________________________________________________________________
+AliRsnMother::AliRsnMother(const AliRsnMother &obj) :
+   TObject(obj),
+   fRefEvent(obj.fRefEvent),
+   fSum(obj.fSum),
+   fSumMC(obj.fSumMC),
+   fRef(obj.fRef),
+   fRefMC(obj.fRefMC)
 {
 //
 // Copy constructor.
-// Initializes all variables to copy values.
 // Does not duplicate pointers.
 //
 
-  Int_t i;
-  for (i = 0; i < 2; i++) fDaughter[i] = obj.fDaughter[i];
+   Int_t i;
+   for (i = 0; i < 2; i++) fDaughter[i] = obj.fDaughter[i];
 }
 
-//_____________________________________________________________________________
+//__________________________________________________________________________________________________
 AliRsnMother& AliRsnMother::operator=(const AliRsnMother &obj)
 {
 //
 // Assignment operator.
-// Initializes all variables to copy values.
 // Does not duplicate pointers.
 //
 
-  Int_t i;
-  
-  fDefaultMass = obj.fDefaultMass;
-  fSum = obj.fSum;
-  fRef = obj.fRef;
-  fSumMC = obj.fSumMC;
-  fRefMC = obj.fRefMC;
-  
-  for (i = 0; i < 2; i++) fDaughter[i] = obj.fDaughter[i];
-
-  return (*this);
+   fSum = obj.fSum;
+   fRef = obj.fRef;
+   fSumMC = obj.fSumMC;
+   fRefMC = obj.fRefMC;
+   fRefEvent = obj.fRefEvent;
+   fDaughter[0] = obj.fDaughter[0];
+   fDaughter[1] = obj.fDaughter[1];
+   
+   return (*this);
 }
 
-//_____________________________________________________________________________
+//__________________________________________________________________________________________________
 AliRsnMother::~AliRsnMother()
 {
 //
@@ -86,108 +83,147 @@ AliRsnMother::~AliRsnMother()
 //
 }
 
-//_____________________________________________________________________________
-Int_t AliRsnMother::CommonMother() const
+//_______________________________________________________________________________________________________________________
+void AliRsnMother::Reset()
 {
 //
-// Checks if the two tracks in the pair have the same mother.
-// This can be known if MC info is present.
-// If the mother label is the same, rhe PDG code of the mother is returned,
-// otherwise the method returns 0.
+// Resets the mother, zeroing all data members.
 //
 
-  // if MC info is not available, the pairs is not true by default
-  if (!fDaughter[0]->GetRefMC() || !fDaughter[1]->GetRefMC()) 
-  {
-    AliInfo("Cannot know if the pairs is true or not because MC Info is not present");
-    return 0;
-  }
+   Int_t i;
+   for (i = 0; i < 2; i++) fDaughter[i] = 0x0;
+   fRefEvent = 0x0;
+   fSum.SetXYZM(0.0, 0.0, 0.0, 0.0);
+   fRef.SetXYZM(0.0, 0.0, 0.0, 0.0);
+}
 
-  // check that labels are the same
-  if (fDaughter[0]->GetParticle()->GetFirstMother() != fDaughter[1]->GetParticle()->GetFirstMother())
-    return 0;
+//__________________________________________________________________________________________________
+Int_t AliRsnMother::CommonMother() const
+{
+//
+// If MC info is available, checks if the two tracks in the pair have the same mother.
+// If the mother label is the same, the function returns the PDG code of mother,
+// otherwise it returns 0.
+// The two arguments passed by reference contain the GEANT labels of the mother
+// of the two particles to which the two daughters point. This is for being able 
+// to check if they are really coming from a resonance (indexes >= 0) or not.
+//
 
-  // if we reach this point, the two tracks have the same mother
-  // let's check now the PDG code of this common mother
-  return TMath::Abs(fDaughter[0]->GetMotherPDG());
+   Int_t m1  = fDaughter[0]->GetMother();
+   Int_t m2  = fDaughter[1]->GetMother();
+   Int_t out = 0;
+   
+   // a true mother makes sense only if both mothers
+   // are not-negative and equal
+   if (m1 >= 0 && m2 >= 0 && m1 == m2) {
+      out = TMath::Abs(fDaughter[0]->GetMotherPDG());
+      AliDebugClass(1, Form("Mothers are: %d %d --> EQUAL (PDG = %d)", m1, m2, out));
+   } 
+   
+   return out;
 }
 
-//_____________________________________________________________________________
-void AliRsnMother::SetDaughters
-(AliRsnDaughter *d0, Double_t mass0, AliRsnDaughter *d1, Double_t mass1)
+//__________________________________________________________________________________________________
+Double_t AliRsnMother::AngleToLeading(Bool_t &success)
 {
 //
-// Sets the pair defined in this usind tso passed daughters and two masses
-// which will be assigned to them, in order to recompute their 4-momenta
-// and sum them into the datamembers of this object.
+// Compute the angle betwee this and the leading particls
+// of the reference event (if this was set properly).
+// In case one of the two daughters is the leading, return
+// a meaningless value, in order to skip this pair.
+// if second argument is kTRUE, use MC values.
 //
 
-  if (d0) fDaughter[0] = d0;
-  if (d1) fDaughter[1] = d1;
-  
-  if (!d0 || !d1) return;
-  
-  fDaughter[0]->SetMass(mass0);
-  fDaughter[1]->SetMass(mass1);
-  
-  fSum   = fDaughter[0]->P(kFALSE) + fDaughter[1]->P(kFALSE);
-  fSumMC = fDaughter[0]->P(kTRUE)  + fDaughter[1]->P(kTRUE);
-  
-  fRef  .SetXYZM(fSum  .X(), fSum  .Y(), fSum  .Z(), fDefaultMass);
-  fRefMC.SetXYZM(fSumMC.X(), fSumMC.Y(), fSumMC.Z(), fDefaultMass);
+   if (!fRefEvent) {
+      success = kFALSE;
+      return -99.0;
+   }
+   
+   Int_t id1 = fDaughter[0]->GetID();
+   Int_t id2 = fDaughter[1]->GetID();
+   Int_t idL = fRefEvent->GetLeadingIndex();
+   
+   if (id1 == idL || id2 == idL) {
+      success = kFALSE;
+      return -99.0;
+   }
+   
+   AliRsnDaughter leading = fRefEvent->GetDaughter(idL, kFALSE);
+   AliVParticle  *ref     = leading.GetRef();
+   Double_t       angle   = ref->Phi() - fSum.Phi();
+   
+   //return angle w.r.t. leading particle in the range -pi/2, 3/2pi
+   while (angle >= 1.5 * TMath::Pi()) angle -= 2 * TMath::Pi();
+   while (angle < -0.5 * TMath::Pi()) angle += 2 * TMath::Pi();
+   success = kTRUE;
+   
+   return angle;
 }
 
-//_____________________________________________________________________________
-void AliRsnMother::ResetPair()
+//__________________________________________________________________________________________________
+void AliRsnMother::ComputeSum(Double_t mass0, Double_t mass1, Double_t motherMass)
 {
 //
-// Resets the mother, nullifying all data members
+// Sets the masses for the 4-momenta of the daughters and then
+// sums them, taking into account that the space part is set to
+// each of them when the reference object is set (see AliRsnDaughter::SetRef)
 //
 
-  Int_t i;
-  for (i = 0; i < 2; i++) fDaughter[i] = 0x0;
-  
-  fSum  .SetXYZM(0.0, 0.0, 0.0, 0.0);
-  fRef  .SetXYZM(0.0, 0.0, 0.0, 0.0);
-  fSumMC.SetXYZM(0.0, 0.0, 0.0, 0.0);
-  fRefMC.SetXYZM(0.0, 0.0, 0.0, 0.0);
+   fDaughter[0]->FillP(mass0);
+   fDaughter[1]->FillP(mass1);
+
+   // sum
+   fSum   = fDaughter[0]->Prec() + fDaughter[1]->Prec();
+   fSumMC = fDaughter[0]->Psim() + fDaughter[1]->Psim();
+   
+   // reference
+   fRef.SetXYZM(fSum.X(), fSum.Y(), fSum.Z(), motherMass);
+   fRefMC.SetXYZM(fSumMC.X(), fSumMC.Y(), fSumMC.Z(), motherMass);
 }
 
-//_____________________________________________________________________________
-Double_t AliRsnMother::ThetaStar(Bool_t first, Bool_t useMC)
+//__________________________________________________________________________________________________
+Double_t AliRsnMother::CosThetaStar(Bool_t first, Bool_t useMC)
 {
 //
-// Returns the theta* as the angle of the first daughter
-// w.r. to the mother momentum, in its rest frame
+// Computes the cosine of theta*, which is the angle of one of the daughters
+// with respect to the total momentum of the resonance, in its rest frame.
+// The arguments are needed to choose which of the daughters one want to use
+// and if reconstructed or MC momentum must be used.
+// [Contribution from Z. Feckova]
 //
 
-  TLorentzVector &mother   = (useMC ? fSumMC : fSum);
-  TLorentzVector &daughter = (first ? fDaughter[0]->P() : fDaughter[1]->P());
-  
-  Double_t theta = mother.Vect().Theta();
-  Double_t phi   = mother.Vect().Phi();
-  Double_t beta  = mother.Beta();
-  
-  Double_t betax = beta * TMath::Sin(theta) * TMath::Cos(phi);
-  Double_t betay = beta * TMath::Sin(theta) * TMath::Sin(phi);
-  Double_t betaz = beta * TMath::Cos(theta);
-  Double_t gamx  = 1.0 / TMath::Sqrt(1.0 - betax*betax);
-  Double_t gamy  = 1.0 / TMath::Sqrt(1.0 - betay*betay);
-  Double_t gamz  = 1.0 / TMath::Sqrt(1.0 - betaz*betaz);
-  
-  Double_t pstarx = gamx * (daughter.Vect().X() - betax * daughter.E());
-  Double_t pstary = gamy * (daughter.Vect().Y() - betay * daughter.E());
-  Double_t pstarz = gamz * (daughter.Vect().Z() - betaz * daughter.E());
-  
-  TVector3 dstar(pstarx, pstary, pstarz);
-  TVector3 vnorm(mother.Vect().Y() / mother.Perp(), mother.Vect().X() / mother.Perp(), 0.0);
-  
-  Double_t cosThetaStar = dstar.Dot(vnorm) / dstar.Mag();
-  
-  return cosThetaStar;
+   TLorentzVector &mother    = Sum(useMC);
+   TLorentzVector &daughter0 = (first ? fDaughter[0]->P(useMC) : fDaughter[1]->P(useMC));
+   TLorentzVector &daughter1 = (first ? fDaughter[1]->P(useMC) : fDaughter[0]->P(useMC));
+   TVector3 momentumM(mother.Vect());
+   TVector3 normal(mother.Y() / momentumM.Mag(), -mother.X() / momentumM.Mag(), 0.0);
+
+   // Computes first the invariant mass of the mother
+   Double_t mass0      = fDaughter[0]->P(useMC).M();
+   Double_t mass1      = fDaughter[1]->P(useMC).M();
+   Double_t p0         = daughter0.Vect().Mag();
+   Double_t p1         = daughter1.Vect().Mag();
+   Double_t E0         = TMath::Sqrt(mass0 * mass0 + p0 * p0);
+   Double_t E1         = TMath::Sqrt(mass1 * mass1 + p1 * p1);
+   Double_t MotherMass = TMath::Sqrt((E0 + E1) * (E0 + E1) - (p0 * p0 + 2.0 * daughter0.Vect().Dot(daughter1.Vect()) + p1 * p1));
+   MotherMass = fSum.M();
+
+   // Computes components of beta
+   Double_t betaX = -mother.X() / mother.E();
+   Double_t betaY = -mother.Y() / mother.E();
+   Double_t betaZ = -mother.Z() / mother.E();
+
+   // Computes Lorentz transformation of the momentum of the first daughter
+   // into the rest frame of the mother and theta*
+   daughter0.Boost(betaX, betaY, betaZ);
+   TVector3 momentumD = daughter0.Vect();
+
+   Double_t cosThetaStar = normal.Dot(momentumD) / momentumD.Mag();
+
+   return cosThetaStar;
 }
 
-//_____________________________________________________________________________
+//__________________________________________________________________________________________________
 void AliRsnMother::PrintInfo(const Option_t * /*option*/) const
 {
 //
@@ -195,16 +231,16 @@ void AliRsnMother::PrintInfo(const Option_t * /*option*/) const
 // The options are passed to the AliRsnDaughter::Print() method
 //
 
-  AliInfo("======== BEGIN PAIR INFO ===========");
-  AliInfo("Track #1");
-  fDaughter[0]->Print();
-  AliInfo("Track #2");
-  fDaughter[1]->Print();
-  AliInfo("========= END PAIR INFO ===========");
+   AliInfo("======== BEGIN PAIR INFO ===========");
+   AliInfo("Track #1");
+   fDaughter[0]->Print();
+   AliInfo("Track #2");
+   fDaughter[1]->Print();
+   AliInfo("========= END PAIR INFO ===========");
 }
 
-//_____________________________________________________________________________
-Bool_t AliRsnMother::CheckPair() const
+//__________________________________________________________________________________________________
+Bool_t AliRsnMother::CheckPair(Bool_t checkMC) const
 {
 //
 // Checks that the pair is well initialized:
@@ -212,61 +248,17 @@ Bool_t AliRsnMother::CheckPair() const
 // - if MC is required, both daughters have a MC reference
 //
 
-  if (!fDaughter[0] || !fDaughter[1]) 
-  {
-    AliError("One of the two tracks is NULL in this pair!");
-    return kFALSE;
-  }
-  
-  if (fUseMC)
-  {
-    if (fDaughter[0]->GetRefMC() == 0x0 || fDaughter[1]->GetRefMC() == 0x0)
-    {
-      AliError("Required MC info but not all MC refs are available");
+   if (!fDaughter[0] || !fDaughter[1]) {
+      AliError("One of the two tracks is NULL in this pair!");
       return kFALSE;
-    }
-  }
-  
-  return kTRUE;
-}
+   }
 
-//_____________________________________________________________________________
-Bool_t AliRsnMother::MatchesDef(AliRsnPairDef *def)
-{
-//
-// Checks if the daughters, in any order, do match a given decay channel,
-// using the specified identification method, which can be the 'true' one
-// or the 'realistic' one only.
-//
+   if (checkMC) {
+      if (fDaughter[0]->GetRefMC() == 0x0 || fDaughter[1]->GetRefMC() == 0x0) {
+         AliError("Required MC info but not all MC refs are available");
+         return kFALSE;
+      }
+   }
 
-  if (!def) return kFALSE;
-  if (!fDaughter[0]->GetRefMC()) return kFALSE;
-  if (!fDaughter[1]->GetRefMC()) return kFALSE;
-
-  Bool_t decayMatch = kFALSE;
-  Int_t  pdg[2], ref[2];
-  pdg[0] = fDaughter[0]->GetRefMC()->Particle()->GetPdgCode();
-  pdg[1] = fDaughter[1]->GetRefMC()->Particle()->GetPdgCode();
-  ref[0] = TMath::Abs(AliPID::ParticleCode(def->GetPID(0)));
-  ref[1] = TMath::Abs(AliPID::ParticleCode(def->GetPID(1)));
-
-  // check #1:
-  // if first member of pairDef has same sign as first member of this,
-  // daughter[0] perfect PID must match first member of pairDef
-  // daughter[1] perfect PID must march second member of pairDef
-  if (fDaughter[0]->IsSign(def->GetCharge(0)) && fDaughter[1]->IsSign(def->GetCharge(1))) 
-  {
-    decayMatch = (pdg[0] == ref[0] && pdg[1] == ref[1]);
-  }
-
-  // check #2:
-  // if first member of pairDef has same sign as second member of this,
-  // daughter[0] perfect PID must match second member of pairDef
-  // daughter[1] perfect PID must march first member of pairDef
-  if (fDaughter[1]->IsSign(def->GetCharge(0)) && fDaughter[0]->IsSign(def->GetCharge(1))) 
-  {
-    decayMatch = (pdg[0] == ref[1] && pdg[1] == ref[0]);
-  }
-
-  return (decayMatch && (CommonMother() == def->GetMotherPDG()));
+   return kTRUE;
 }