[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskRhoMass.cxx

// $Id$
//
// Calculation of rho mass from a collection of jets.
// If scale function is given the scaled rho will be exported
// with the name as "fOutRhoMassName".Apppend("_Scaled").
//
// Authors: M. Verweij

#include "AliAnalysisTaskRhoMass.h"

#include <TClonesArray.h>
#include <TMath.h>

#include "AliAnalysisManager.h"
#include "AliEmcalJet.h"
#include "AliLog.h"
#include "AliRhoParameter.h"

ClassImp(AliAnalysisTaskRhoMass)

//________________________________________________________________________
AliAnalysisTaskRhoMass::AliAnalysisTaskRhoMass() : 
  AliAnalysisTaskRhoMassBase("AliAnalysisTaskRhoMass"),
  fNExclLeadJets(0),
  fJetRhoMassType(kMd),
  fPionMassClusters(kFALSE),
  fHistMdAreavsCent(0)
{
  // Constructor.
}

//________________________________________________________________________
AliAnalysisTaskRhoMass::AliAnalysisTaskRhoMass(const char *name, Bool_t histo) :
  AliAnalysisTaskRhoMassBase(name, histo),
  fNExclLeadJets(0),
  fJetRhoMassType(kMd),
  fPionMassClusters(kFALSE),
  fHistMdAreavsCent(0)
{
  // Constructor.
}

//________________________________________________________________________
void AliAnalysisTaskRhoMass::UserCreateOutputObjects()
{
  // User create output objects, called at the beginning of the analysis.

  if (!fCreateHisto)
    return;
  
  AliAnalysisTaskRhoMassBase::UserCreateOutputObjects();

  fHistMdAreavsCent = new TH2F("fHistMdAreavsCent", "fHistMdAreavsCent", 101, -1,  100, fNbins, fMinBinPt, fMaxBinPt/2.);
  fHistMdAreavsCent->GetXaxis()->SetTitle("Centrality (%)");
  fHistMdAreavsCent->GetYaxis()->SetTitle("#rho_{m} (GeV/c * rad^{-1})");
  fOutput->Add(fHistMdAreavsCent);
}


//________________________________________________________________________
Bool_t AliAnalysisTaskRhoMass::Run() 
{
  // Run the analysis.

  fOutRhoMass->SetVal(0);
  if (fOutRhoMassScaled)
    fOutRhoMassScaled->SetVal(0);

  if (!fJets)
    return kFALSE;

  const Int_t Njets   = fJets->GetEntries();

  Int_t maxJetIds[]   = {-1, -1};
  Float_t maxJetPts[] = { 0,  0};

  if (fNExclLeadJets > 0) {
    for (Int_t ij = 0; ij < Njets; ++ij) {
      AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
      if (!jet) {
	AliError(Form("%s: Could not receive jet %d", GetName(), ij));
	continue;
      } 

      if (!AcceptJet(jet))
        continue;

      if (jet->Pt() > maxJetPts[0]) {
	maxJetPts[1] = maxJetPts[0];
	maxJetIds[1] = maxJetIds[0];
	maxJetPts[0] = jet->Pt();
	maxJetIds[0] = ij;
      } else if (jet->Pt() > maxJetPts[1]) {
	maxJetPts[1] = jet->Pt();
	maxJetIds[1] = ij;
      }
    }
    if (fNExclLeadJets < 2) {
      maxJetIds[1] = -1;
      maxJetPts[1] = 0;
    }
  }

  static Double_t rhomvec[999];
  static Double_t Evec[999];
  static Double_t Mvec[999];
  Int_t NjetAcc = 0;

  // push all jets within selected acceptance into stack
  for (Int_t iJets = 0; iJets < Njets; ++iJets) {

    // exlcuding lead jets
    if (iJets == maxJetIds[0] || iJets == maxJetIds[1])
      continue;

    AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
    if (!jet) {
      AliError(Form("%s: Could not receive jet %d", GetName(), iJets));
      continue;
    } 

    if (!AcceptJet(jet))
      continue;

    // Double_t sumM = GetSumMConstituents(jet);
    // Double_t sumPt = GetSumPtConstituents(jet);
    if(jet->Area()>0.) {// && (jet->M()*jet->M() + jet->Pt()*jet->Pt())>0.) {
      //rhomvec[NjetAcc] = (TMath::Sqrt(sumM*sumM + sumPt*sumPt) - sumPt ) / jet->Area();
      // rhomvec[NjetAcc] = (TMath::Sqrt(jet->M()*jet->M() + jet->Pt()*jet->Pt()) - jet->Pt() ) / jet->Area();
      rhomvec[NjetAcc] = GetMd(jet) / jet->Area();
      fHistMdAreavsCent->Fill(fCent,rhomvec[NjetAcc]);
      Evec[NjetAcc] = jet->E();
      Mvec[NjetAcc] = jet->M();
      ++NjetAcc;
    }
  }

  if (NjetAcc > 0) {
    //find median value
    Double_t rhom = TMath::Median(NjetAcc, rhomvec);
    fOutRhoMass->SetVal(rhom);

    Int_t Ntracks = fTracks->GetEntries();
    Double_t meanM = TMath::Mean(NjetAcc, Mvec);
    Double_t meanE = TMath::Mean(NjetAcc, Evec);
    Double_t gamma = 0.;
    if(meanM>0.) gamma = meanE/meanM;
    fHistGammaVsNtrack->Fill(Ntracks,gamma);

    if (fOutRhoMassScaled) {
      Double_t rhomScaled = rhom * GetScaleFactor(fCent);
      fOutRhoMassScaled->SetVal(rhomScaled);
    }
  }

  return kTRUE;
} 

//________________________________________________________________________
Double_t AliAnalysisTaskRhoMass::GetSumMConstituents(AliEmcalJet *jet) {
  
  Double_t sum = 0.;
  
  AliVParticle *vp;
  for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
    vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
    if(!vp) continue;
    sum+=vp->M();
  }
  return sum;
}

//________________________________________________________________________
Double_t AliAnalysisTaskRhoMass::GetSumPtConstituents(AliEmcalJet *jet) {
  
  Double_t sum = 0.;
  
  AliVParticle *vp;
  for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
    vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
    if(!vp) continue;
    sum+=vp->Pt();
  }
  return sum;
}

//________________________________________________________________________
Double_t AliAnalysisTaskRhoMass::GetMd(AliEmcalJet *jet) {
  //get md as defined in http://arxiv.org/pdf/1211.2811.pdf
  Double_t sum = 0.;
  Double_t px = 0.;
  Double_t py = 0.;
  Double_t pz = 0.;
  Double_t E = 0.;

  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);
      Double_t m = 0.;
      if(fPionMassClusters) m = 0.13957;
      if(fJetRhoMassType==kMd) sum += TMath::Sqrt(m*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;
    sum = TMath::Sqrt(m2 + pt*pt) - pt;
  }
  return sum;
}
Commit	Line	Data
4d3b366f	1	// $Id$
	2	//
	3	// Calculation of rho mass from a collection of jets.
	4	// If scale function is given the scaled rho will be exported
	5	// with the name as "fOutRhoMassName".Apppend("_Scaled").
	6	//
	7	// Authors: M. Verweij
	8
	9	#include "AliAnalysisTaskRhoMass.h"
	10
	11	#include <TClonesArray.h>
	12	#include <TMath.h>
	13
	14	#include "AliAnalysisManager.h"
	15	#include "AliEmcalJet.h"
	16	#include "AliLog.h"
	17	#include "AliRhoParameter.h"
	18
	19	ClassImp(AliAnalysisTaskRhoMass)
	20
	21	//________________________________________________________________________
	22	AliAnalysisTaskRhoMass::AliAnalysisTaskRhoMass() :
	23	AliAnalysisTaskRhoMassBase("AliAnalysisTaskRhoMass"),
	24	fNExclLeadJets(0),
	25	fJetRhoMassType(kMd),
1e9d5955	26	fPionMassClusters(kFALSE),
4d3b366f	27	fHistMdAreavsCent(0)
	28	{
	29	// Constructor.
	30	}
	31
	32	//________________________________________________________________________
	33	AliAnalysisTaskRhoMass::AliAnalysisTaskRhoMass(const char *name, Bool_t histo) :
	34	AliAnalysisTaskRhoMassBase(name, histo),
	35	fNExclLeadJets(0),
	36	fJetRhoMassType(kMd),
1e9d5955	37	fPionMassClusters(kFALSE),
4d3b366f	38	fHistMdAreavsCent(0)
	39	{
	40	// Constructor.
	41	}
	42
	43	//________________________________________________________________________
	44	void AliAnalysisTaskRhoMass::UserCreateOutputObjects()
	45	{
	46	// User create output objects, called at the beginning of the analysis.
	47
	48	if (!fCreateHisto)
	49	return;
	50
	51	AliAnalysisTaskRhoMassBase::UserCreateOutputObjects();
	52
	53	fHistMdAreavsCent = new TH2F("fHistMdAreavsCent", "fHistMdAreavsCent", 101, -1, 100, fNbins, fMinBinPt, fMaxBinPt/2.);
	54	fHistMdAreavsCent->GetXaxis()->SetTitle("Centrality (%)");
	55	fHistMdAreavsCent->GetYaxis()->SetTitle("#rho_{m} (GeV/c * rad^{-1})");
	56	fOutput->Add(fHistMdAreavsCent);
	57	}
	58
	59
	60	//________________________________________________________________________
	61	Bool_t AliAnalysisTaskRhoMass::Run()
	62	{
	63	// Run the analysis.
	64
	65	fOutRhoMass->SetVal(0);
	66	if (fOutRhoMassScaled)
	67	fOutRhoMassScaled->SetVal(0);
	68
	69	if (!fJets)
	70	return kFALSE;
	71
	72	const Int_t Njets = fJets->GetEntries();
	73
	74	Int_t maxJetIds[] = {-1, -1};
	75	Float_t maxJetPts[] = { 0, 0};
	76
	77	if (fNExclLeadJets > 0) {
	78	for (Int_t ij = 0; ij < Njets; ++ij) {
	79	AliEmcalJet jet = static_cast<AliEmcalJet>(fJets->At(ij));
	80	if (!jet) {
	81	AliError(Form("%s: Could not receive jet %d", GetName(), ij));
	82	continue;
	83	}
	84
	85	if (!AcceptJet(jet))
	86	continue;
	87
	88	if (jet->Pt() > maxJetPts[0]) {
	89	maxJetPts[1] = maxJetPts[0];
	90	maxJetIds[1] = maxJetIds[0];
	91	maxJetPts[0] = jet->Pt();
	92	maxJetIds[0] = ij;
	93	} else if (jet->Pt() > maxJetPts[1]) {
	94	maxJetPts[1] = jet->Pt();
	95	maxJetIds[1] = ij;
	96	}
	97	}
	98	if (fNExclLeadJets < 2) {
	99	maxJetIds[1] = -1;
	100	maxJetPts[1] = 0;
	101	}
102	}
103
104	static Double_t rhomvec[999];
105	static Double_t Evec[999];
106	static Double_t Mvec[999];
107	Int_t NjetAcc = 0;
108
109	// push all jets within selected acceptance into stack
110	for (Int_t iJets = 0; iJets < Njets; ++iJets) {
111
112	// exlcuding lead jets
113	if (iJets == maxJetIds[0] \|\| iJets == maxJetIds[1])
114	continue;
115
116	AliEmcalJet jet = static_cast<AliEmcalJet>(fJets->At(iJets));
117	if (!jet) {
118	AliError(Form("%s: Could not receive jet %d", GetName(), iJets));
119	continue;
120	}
121
122	if (!AcceptJet(jet))
123	continue;
124
125	// Double_t sumM = GetSumMConstituents(jet);
126	// Double_t sumPt = GetSumPtConstituents(jet);
127	if(jet->Area()>0.) {// && (jet->M()jet->M() + jet->Pt()jet->Pt())>0.) {
128	//rhomvec[NjetAcc] = (TMath::Sqrt(sumMsumM + sumPtsumPt) - sumPt ) / jet->Area();
129	// rhomvec[NjetAcc] = (TMath::Sqrt(jet->M()jet->M() + jet->Pt()jet->Pt()) - jet->Pt() ) / jet->Area();
130	rhomvec[NjetAcc] = GetMd(jet) / jet->Area();
131	fHistMdAreavsCent->Fill(fCent,rhomvec[NjetAcc]);
132	Evec[NjetAcc] = jet->E();
133	Mvec[NjetAcc] = jet->M();
134	++NjetAcc;
135	}
136	}
137
138	if (NjetAcc > 0) {
139	//find median value
140	Double_t rhom = TMath::Median(NjetAcc, rhomvec);
141	fOutRhoMass->SetVal(rhom);
142
143	Int_t Ntracks = fTracks->GetEntries();
144	Double_t meanM = TMath::Mean(NjetAcc, Mvec);
145	Double_t meanE = TMath::Mean(NjetAcc, Evec);
146	Double_t gamma = 0.;
147	if(meanM>0.) gamma = meanE/meanM;
148	fHistGammaVsNtrack->Fill(Ntracks,gamma);
149
150	if (fOutRhoMassScaled) {
151	Double_t rhomScaled = rhom * GetScaleFactor(fCent);
152	fOutRhoMassScaled->SetVal(rhomScaled);
153	}
154	}
155
156	return kTRUE;
157	}
158
159	//________________________________________________________________________
160	Double_t AliAnalysisTaskRhoMass::GetSumMConstituents(AliEmcalJet *jet) {
161
162	Double_t sum = 0.;
163
164	AliVParticle *vp;
165	for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
166	vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
167	if(!vp) continue;
168	sum+=vp->M();
169	}
170	return sum;
171	}
172
173	//________________________________________________________________________
174	Double_t AliAnalysisTaskRhoMass::GetSumPtConstituents(AliEmcalJet *jet) {
175
176	Double_t sum = 0.;
177
178	AliVParticle *vp;
179	for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
180	vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
181	if(!vp) continue;
182	sum+=vp->Pt();
183	}
184	return sum;
185	}
186
187	//________________________________________________________________________
188	Double_t AliAnalysisTaskRhoMass::GetMd(AliEmcalJet *jet) {
189	//get md as defined in http://arxiv.org/pdf/1211.2811.pdf
190	Double_t sum = 0.;
191	Double_t px = 0.;
192	Double_t py = 0.;
193	Double_t pz = 0.;
194	Double_t E = 0.;
85a6a0ff	195
	196	if (fTracks) {
	197	AliVParticle *vp;
	198	for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
	199	vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
	200	if(!vp) continue;
	201	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)
	202	else if(fJetRhoMassType==kMdP) sum += TMath::Sqrt(vp->M()vp->M() + vp->P()vp->P()) - vp->P();
	203	else if(fJetRhoMassType==kMd4) {
	204	px+=vp->Px();
	205	py+=vp->Py();
	206	pz+=vp->Pz();
	207	E+=vp->E();
	208	}
4d3b366f	209	}
4d3b366f	210	}
85a6a0ff	211
	212	if (fCaloClusters) {
	213	AliVCluster *vp;
	214	for(Int_t icc=0; icc<jet->GetNumberOfClusters(); icc++) {
	215	vp = static_cast<AliVCluster*>(jet->ClusterAt(icc, fCaloClusters));
	216	if(!vp) continue;
	217	TLorentzVector nPart;
	218	vp->GetMomentum(nPart, fVertex);
1e9d5955	219	Double_t m = 0.;
	220	if(fPionMassClusters) m = 0.13957;
	221	if(fJetRhoMassType==kMd) sum += TMath::Sqrt(mm + nPart.Pt()nPart.Pt()) - nPart.Pt();
85a6a0ff	222	else if(fJetRhoMassType==kMdP) sum += TMath::Sqrt(nPart.M()nPart.M() + nPart.P()nPart.P()) - nPart.P();
	223	else if(fJetRhoMassType==kMd4) {
	224	px+=nPart.Px();
	225	py+=nPart.Py();
	226	pz+=nPart.Pz();
	227	E+=nPart.E();
	228	}
	229	}
	230	}
	231
4d3b366f	232	if(fJetRhoMassType==kMd4) {
	233	Double_t pt = TMath::Sqrt(pxpx + pypy);
	234	Double_t m2 = EE - ptpt - pz*pz;
	235	sum = TMath::Sqrt(m2 + pt*pt) - pt;
	236	}
	237	return sum;
	238	}