From 85a6a0ff5c77c40326cd40ff9fceab18b270b1d1 Mon Sep 17 00:00:00 2001
From: mverweij <marta.verweij@cern.ch>
Date: Tue, 2 Sep 2014 23:35:07 +0200
Subject: [PATCH] include neutral constituents in rho_m calculation

---
 .../UserTasks/AliAnalysisTaskRhoMass.cxx      | 45 ++++++++++++++-----
 1 file changed, 34 insertions(+), 11 deletions(-)

diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoMass.cxx b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoMass.cxx
index 4ea19fbe28d..161790a163e 100644
--- a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoMass.cxx
+++ b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoMass.cxx
@@ -190,19 +190,42 @@ Double_t AliAnalysisTaskRhoMass::GetMd(AliEmcalJet *jet) {
   Double_t py = 0.;
   Double_t pz = 0.;
   Double_t E = 0.;
-  AliVParticle *vp;
-  for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
-    vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
-    if(!vp) continue;
-    if(fJetRhoMassType==kMd) sum += TMath::Sqrt(vp->M()*vp->M() + vp->Pt()*vp->Pt()) - vp->Pt(); //sqrt(E^2-P^2+pt^2)=sqrt(E^2-pz^2)
-    else if(fJetRhoMassType==kMdP) sum += TMath::Sqrt(vp->M()*vp->M() + vp->P()*vp->P()) - vp->P();
-    else if(fJetRhoMassType==kMd4) {
-      px+=vp->Px();
-      py+=vp->Py();
-      pz+=vp->Pz();
-      E+=vp->E();
+
+  if (fTracks) {
+    AliVParticle *vp;
+    for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
+      vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
+      if(!vp) continue;
+      if(fJetRhoMassType==kMd) sum += TMath::Sqrt(vp->M()*vp->M() + vp->Pt()*vp->Pt()) - vp->Pt(); //sqrt(E^2-P^2+pt^2)=sqrt(E^2-pz^2)
+      else if(fJetRhoMassType==kMdP) sum += TMath::Sqrt(vp->M()*vp->M() + vp->P()*vp->P()) - vp->P();
+      else if(fJetRhoMassType==kMd4) {
+	px+=vp->Px();
+	py+=vp->Py();
+	pz+=vp->Pz();
+	E+=vp->E();
+      }
     }
   }
+
+  if (fCaloClusters) {
+    AliVCluster *vp;
+    for(Int_t icc=0; icc<jet->GetNumberOfClusters(); icc++) {
+      vp = static_cast<AliVCluster*>(jet->ClusterAt(icc, fCaloClusters));
+      if(!vp) continue;
+      TLorentzVector nPart;
+      vp->GetMomentum(nPart, fVertex);
+      
+      if(fJetRhoMassType==kMd) sum += TMath::Sqrt(nPart.M()*nPart.M() + nPart.Pt()*nPart.Pt()) - nPart.Pt();
+      else if(fJetRhoMassType==kMdP) sum += TMath::Sqrt(nPart.M()*nPart.M() + nPart.P()*nPart.P()) - nPart.P();
+      else if(fJetRhoMassType==kMd4) {
+	px+=nPart.Px();
+	py+=nPart.Py();
+	pz+=nPart.Pz();
+	E+=nPart.E();
+      }
+    }
+  }
+
   if(fJetRhoMassType==kMd4) {
     Double_t pt = TMath::Sqrt(px*px + py*py);
     Double_t m2 = E*E - pt*pt - pz*pz;
-- 
2.43.0