ROOT logo
// $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;
}
 AliAnalysisTaskRhoMass.cxx:1
 AliAnalysisTaskRhoMass.cxx:2
 AliAnalysisTaskRhoMass.cxx:3
 AliAnalysisTaskRhoMass.cxx:4
 AliAnalysisTaskRhoMass.cxx:5
 AliAnalysisTaskRhoMass.cxx:6
 AliAnalysisTaskRhoMass.cxx:7
 AliAnalysisTaskRhoMass.cxx:8
 AliAnalysisTaskRhoMass.cxx:9
 AliAnalysisTaskRhoMass.cxx:10
 AliAnalysisTaskRhoMass.cxx:11
 AliAnalysisTaskRhoMass.cxx:12
 AliAnalysisTaskRhoMass.cxx:13
 AliAnalysisTaskRhoMass.cxx:14
 AliAnalysisTaskRhoMass.cxx:15
 AliAnalysisTaskRhoMass.cxx:16
 AliAnalysisTaskRhoMass.cxx:17
 AliAnalysisTaskRhoMass.cxx:18
 AliAnalysisTaskRhoMass.cxx:19
 AliAnalysisTaskRhoMass.cxx:20
 AliAnalysisTaskRhoMass.cxx:21
 AliAnalysisTaskRhoMass.cxx:22
 AliAnalysisTaskRhoMass.cxx:23
 AliAnalysisTaskRhoMass.cxx:24
 AliAnalysisTaskRhoMass.cxx:25
 AliAnalysisTaskRhoMass.cxx:26
 AliAnalysisTaskRhoMass.cxx:27
 AliAnalysisTaskRhoMass.cxx:28
 AliAnalysisTaskRhoMass.cxx:29
 AliAnalysisTaskRhoMass.cxx:30
 AliAnalysisTaskRhoMass.cxx:31
 AliAnalysisTaskRhoMass.cxx:32
 AliAnalysisTaskRhoMass.cxx:33
 AliAnalysisTaskRhoMass.cxx:34
 AliAnalysisTaskRhoMass.cxx:35
 AliAnalysisTaskRhoMass.cxx:36
 AliAnalysisTaskRhoMass.cxx:37
 AliAnalysisTaskRhoMass.cxx:38
 AliAnalysisTaskRhoMass.cxx:39
 AliAnalysisTaskRhoMass.cxx:40
 AliAnalysisTaskRhoMass.cxx:41
 AliAnalysisTaskRhoMass.cxx:42
 AliAnalysisTaskRhoMass.cxx:43
 AliAnalysisTaskRhoMass.cxx:44
 AliAnalysisTaskRhoMass.cxx:45
 AliAnalysisTaskRhoMass.cxx:46
 AliAnalysisTaskRhoMass.cxx:47
 AliAnalysisTaskRhoMass.cxx:48
 AliAnalysisTaskRhoMass.cxx:49
 AliAnalysisTaskRhoMass.cxx:50
 AliAnalysisTaskRhoMass.cxx:51
 AliAnalysisTaskRhoMass.cxx:52
 AliAnalysisTaskRhoMass.cxx:53
 AliAnalysisTaskRhoMass.cxx:54
 AliAnalysisTaskRhoMass.cxx:55
 AliAnalysisTaskRhoMass.cxx:56
 AliAnalysisTaskRhoMass.cxx:57
 AliAnalysisTaskRhoMass.cxx:58
 AliAnalysisTaskRhoMass.cxx:59
 AliAnalysisTaskRhoMass.cxx:60
 AliAnalysisTaskRhoMass.cxx:61
 AliAnalysisTaskRhoMass.cxx:62
 AliAnalysisTaskRhoMass.cxx:63
 AliAnalysisTaskRhoMass.cxx:64
 AliAnalysisTaskRhoMass.cxx:65
 AliAnalysisTaskRhoMass.cxx:66
 AliAnalysisTaskRhoMass.cxx:67
 AliAnalysisTaskRhoMass.cxx:68
 AliAnalysisTaskRhoMass.cxx:69
 AliAnalysisTaskRhoMass.cxx:70
 AliAnalysisTaskRhoMass.cxx:71
 AliAnalysisTaskRhoMass.cxx:72
 AliAnalysisTaskRhoMass.cxx:73
 AliAnalysisTaskRhoMass.cxx:74
 AliAnalysisTaskRhoMass.cxx:75
 AliAnalysisTaskRhoMass.cxx:76
 AliAnalysisTaskRhoMass.cxx:77
 AliAnalysisTaskRhoMass.cxx:78
 AliAnalysisTaskRhoMass.cxx:79
 AliAnalysisTaskRhoMass.cxx:80
 AliAnalysisTaskRhoMass.cxx:81
 AliAnalysisTaskRhoMass.cxx:82
 AliAnalysisTaskRhoMass.cxx:83
 AliAnalysisTaskRhoMass.cxx:84
 AliAnalysisTaskRhoMass.cxx:85
 AliAnalysisTaskRhoMass.cxx:86
 AliAnalysisTaskRhoMass.cxx:87
 AliAnalysisTaskRhoMass.cxx:88
 AliAnalysisTaskRhoMass.cxx:89
 AliAnalysisTaskRhoMass.cxx:90
 AliAnalysisTaskRhoMass.cxx:91
 AliAnalysisTaskRhoMass.cxx:92
 AliAnalysisTaskRhoMass.cxx:93
 AliAnalysisTaskRhoMass.cxx:94
 AliAnalysisTaskRhoMass.cxx:95
 AliAnalysisTaskRhoMass.cxx:96
 AliAnalysisTaskRhoMass.cxx:97
 AliAnalysisTaskRhoMass.cxx:98
 AliAnalysisTaskRhoMass.cxx:99
 AliAnalysisTaskRhoMass.cxx:100
 AliAnalysisTaskRhoMass.cxx:101
 AliAnalysisTaskRhoMass.cxx:102
 AliAnalysisTaskRhoMass.cxx:103
 AliAnalysisTaskRhoMass.cxx:104
 AliAnalysisTaskRhoMass.cxx:105
 AliAnalysisTaskRhoMass.cxx:106
 AliAnalysisTaskRhoMass.cxx:107
 AliAnalysisTaskRhoMass.cxx:108
 AliAnalysisTaskRhoMass.cxx:109
 AliAnalysisTaskRhoMass.cxx:110
 AliAnalysisTaskRhoMass.cxx:111
 AliAnalysisTaskRhoMass.cxx:112
 AliAnalysisTaskRhoMass.cxx:113
 AliAnalysisTaskRhoMass.cxx:114
 AliAnalysisTaskRhoMass.cxx:115
 AliAnalysisTaskRhoMass.cxx:116
 AliAnalysisTaskRhoMass.cxx:117
 AliAnalysisTaskRhoMass.cxx:118
 AliAnalysisTaskRhoMass.cxx:119
 AliAnalysisTaskRhoMass.cxx:120
 AliAnalysisTaskRhoMass.cxx:121
 AliAnalysisTaskRhoMass.cxx:122
 AliAnalysisTaskRhoMass.cxx:123
 AliAnalysisTaskRhoMass.cxx:124
 AliAnalysisTaskRhoMass.cxx:125
 AliAnalysisTaskRhoMass.cxx:126
 AliAnalysisTaskRhoMass.cxx:127
 AliAnalysisTaskRhoMass.cxx:128
 AliAnalysisTaskRhoMass.cxx:129
 AliAnalysisTaskRhoMass.cxx:130
 AliAnalysisTaskRhoMass.cxx:131
 AliAnalysisTaskRhoMass.cxx:132
 AliAnalysisTaskRhoMass.cxx:133
 AliAnalysisTaskRhoMass.cxx:134
 AliAnalysisTaskRhoMass.cxx:135
 AliAnalysisTaskRhoMass.cxx:136
 AliAnalysisTaskRhoMass.cxx:137
 AliAnalysisTaskRhoMass.cxx:138
 AliAnalysisTaskRhoMass.cxx:139
 AliAnalysisTaskRhoMass.cxx:140
 AliAnalysisTaskRhoMass.cxx:141
 AliAnalysisTaskRhoMass.cxx:142
 AliAnalysisTaskRhoMass.cxx:143
 AliAnalysisTaskRhoMass.cxx:144
 AliAnalysisTaskRhoMass.cxx:145
 AliAnalysisTaskRhoMass.cxx:146
 AliAnalysisTaskRhoMass.cxx:147
 AliAnalysisTaskRhoMass.cxx:148
 AliAnalysisTaskRhoMass.cxx:149
 AliAnalysisTaskRhoMass.cxx:150
 AliAnalysisTaskRhoMass.cxx:151
 AliAnalysisTaskRhoMass.cxx:152
 AliAnalysisTaskRhoMass.cxx:153
 AliAnalysisTaskRhoMass.cxx:154
 AliAnalysisTaskRhoMass.cxx:155
 AliAnalysisTaskRhoMass.cxx:156
 AliAnalysisTaskRhoMass.cxx:157
 AliAnalysisTaskRhoMass.cxx:158
 AliAnalysisTaskRhoMass.cxx:159
 AliAnalysisTaskRhoMass.cxx:160
 AliAnalysisTaskRhoMass.cxx:161
 AliAnalysisTaskRhoMass.cxx:162
 AliAnalysisTaskRhoMass.cxx:163
 AliAnalysisTaskRhoMass.cxx:164
 AliAnalysisTaskRhoMass.cxx:165
 AliAnalysisTaskRhoMass.cxx:166
 AliAnalysisTaskRhoMass.cxx:167
 AliAnalysisTaskRhoMass.cxx:168
 AliAnalysisTaskRhoMass.cxx:169
 AliAnalysisTaskRhoMass.cxx:170
 AliAnalysisTaskRhoMass.cxx:171
 AliAnalysisTaskRhoMass.cxx:172
 AliAnalysisTaskRhoMass.cxx:173
 AliAnalysisTaskRhoMass.cxx:174
 AliAnalysisTaskRhoMass.cxx:175
 AliAnalysisTaskRhoMass.cxx:176
 AliAnalysisTaskRhoMass.cxx:177
 AliAnalysisTaskRhoMass.cxx:178
 AliAnalysisTaskRhoMass.cxx:179
 AliAnalysisTaskRhoMass.cxx:180
 AliAnalysisTaskRhoMass.cxx:181
 AliAnalysisTaskRhoMass.cxx:182
 AliAnalysisTaskRhoMass.cxx:183
 AliAnalysisTaskRhoMass.cxx:184
 AliAnalysisTaskRhoMass.cxx:185
 AliAnalysisTaskRhoMass.cxx:186
 AliAnalysisTaskRhoMass.cxx:187
 AliAnalysisTaskRhoMass.cxx:188
 AliAnalysisTaskRhoMass.cxx:189
 AliAnalysisTaskRhoMass.cxx:190
 AliAnalysisTaskRhoMass.cxx:191
 AliAnalysisTaskRhoMass.cxx:192
 AliAnalysisTaskRhoMass.cxx:193
 AliAnalysisTaskRhoMass.cxx:194
 AliAnalysisTaskRhoMass.cxx:195
 AliAnalysisTaskRhoMass.cxx:196
 AliAnalysisTaskRhoMass.cxx:197
 AliAnalysisTaskRhoMass.cxx:198
 AliAnalysisTaskRhoMass.cxx:199
 AliAnalysisTaskRhoMass.cxx:200
 AliAnalysisTaskRhoMass.cxx:201
 AliAnalysisTaskRhoMass.cxx:202
 AliAnalysisTaskRhoMass.cxx:203
 AliAnalysisTaskRhoMass.cxx:204
 AliAnalysisTaskRhoMass.cxx:205
 AliAnalysisTaskRhoMass.cxx:206
 AliAnalysisTaskRhoMass.cxx:207
 AliAnalysisTaskRhoMass.cxx:208
 AliAnalysisTaskRhoMass.cxx:209
 AliAnalysisTaskRhoMass.cxx:210
 AliAnalysisTaskRhoMass.cxx:211
 AliAnalysisTaskRhoMass.cxx:212
 AliAnalysisTaskRhoMass.cxx:213
 AliAnalysisTaskRhoMass.cxx:214
 AliAnalysisTaskRhoMass.cxx:215
 AliAnalysisTaskRhoMass.cxx:216
 AliAnalysisTaskRhoMass.cxx:217
 AliAnalysisTaskRhoMass.cxx:218
 AliAnalysisTaskRhoMass.cxx:219
 AliAnalysisTaskRhoMass.cxx:220
 AliAnalysisTaskRhoMass.cxx:221
 AliAnalysisTaskRhoMass.cxx:222
 AliAnalysisTaskRhoMass.cxx:223
 AliAnalysisTaskRhoMass.cxx:224
 AliAnalysisTaskRhoMass.cxx:225
 AliAnalysisTaskRhoMass.cxx:226
 AliAnalysisTaskRhoMass.cxx:227
 AliAnalysisTaskRhoMass.cxx:228
 AliAnalysisTaskRhoMass.cxx:229
 AliAnalysisTaskRhoMass.cxx:230
 AliAnalysisTaskRhoMass.cxx:231
 AliAnalysisTaskRhoMass.cxx:232
 AliAnalysisTaskRhoMass.cxx:233
 AliAnalysisTaskRhoMass.cxx:234
 AliAnalysisTaskRhoMass.cxx:235
 AliAnalysisTaskRhoMass.cxx:236
 AliAnalysisTaskRhoMass.cxx:237
 AliAnalysisTaskRhoMass.cxx:238