PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoMass.cxx

   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),
  26   fHistMdAreavsCent(0)
  27 {
  28   // Constructor.
  29 }
  30
  31 //________________________________________________________________________
  32 AliAnalysisTaskRhoMass::AliAnalysisTaskRhoMass(const char *name, Bool_t histo) :
  33   AliAnalysisTaskRhoMassBase(name, histo),
  34   fNExclLeadJets(0),
  35   fJetRhoMassType(kMd),
  36   fHistMdAreavsCent(0)
  37 {
  38   // Constructor.
  39 }
  40
  41 //________________________________________________________________________
  42 void AliAnalysisTaskRhoMass::UserCreateOutputObjects()
  43 {
  44   // User create output objects, called at the beginning of the analysis.
  45
  46   if (!fCreateHisto)
  47     return;
  48
  49   AliAnalysisTaskRhoMassBase::UserCreateOutputObjects();
  50
  51   fHistMdAreavsCent = new TH2F("fHistMdAreavsCent", "fHistMdAreavsCent", 101, -1,  100, fNbins, fMinBinPt, fMaxBinPt/2.);
  52   fHistMdAreavsCent->GetXaxis()->SetTitle("Centrality (%)");
  53   fHistMdAreavsCent->GetYaxis()->SetTitle("#rho_{m} (GeV/c * rad^{-1})");
  54   fOutput->Add(fHistMdAreavsCent);
  55 }
  56
  57
  58 //________________________________________________________________________
  59 Bool_t AliAnalysisTaskRhoMass::Run()
  60 {
  61   // Run the analysis.
  62
  63   fOutRhoMass->SetVal(0);
  64   if (fOutRhoMassScaled)
  65     fOutRhoMassScaled->SetVal(0);
  66
  67   if (!fJets)
  68     return kFALSE;
  69
  70   const Int_t Njets   = fJets->GetEntries();
  71
  72   Int_t maxJetIds[]   = {-1, -1};
  73   Float_t maxJetPts[] = { 0,  0};
  74
  75   if (fNExclLeadJets > 0) {
  76     for (Int_t ij = 0; ij < Njets; ++ij) {
  77       AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(ij));
  78       if (!jet) {
  79         AliError(Form("%s: Could not receive jet %d", GetName(), ij));
  80         continue;
  81       }
  82
  83       if (!AcceptJet(jet))
  84         continue;
  85
  86       if (jet->Pt() > maxJetPts[0]) {
  87         maxJetPts[1] = maxJetPts[0];
  88         maxJetIds[1] = maxJetIds[0];
  89         maxJetPts[0] = jet->Pt();
  90         maxJetIds[0] = ij;
  91       } else if (jet->Pt() > maxJetPts[1]) {
  92         maxJetPts[1] = jet->Pt();
  93         maxJetIds[1] = ij;
  94       }
  95     }
  96     if (fNExclLeadJets < 2) {
  97       maxJetIds[1] = -1;
  98       maxJetPts[1] = 0;
  99     }
 100   }
 101
 102   static Double_t rhomvec[999];
 103   static Double_t Evec[999];
 104   static Double_t Mvec[999];
 105   Int_t NjetAcc = 0;
 106
 107   // push all jets within selected acceptance into stack
 108   for (Int_t iJets = 0; iJets < Njets; ++iJets) {
 109
 110     // exlcuding lead jets
 111     if (iJets == maxJetIds[0] || iJets == maxJetIds[1])
 112       continue;
 113
 114     AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(iJets));
 115     if (!jet) {
 116       AliError(Form("%s: Could not receive jet %d", GetName(), iJets));
 117       continue;
 118     }
 119
 120     if (!AcceptJet(jet))
 121       continue;
 122
 123     // Double_t sumM = GetSumMConstituents(jet);
 124     // Double_t sumPt = GetSumPtConstituents(jet);
 125     if(jet->Area()>0.) {// && (jet->M()*jet->M() + jet->Pt()*jet->Pt())>0.) {
 126       //rhomvec[NjetAcc] = (TMath::Sqrt(sumM*sumM + sumPt*sumPt) - sumPt ) / jet->Area();
 127       // rhomvec[NjetAcc] = (TMath::Sqrt(jet->M()*jet->M() + jet->Pt()*jet->Pt()) - jet->Pt() ) / jet->Area();
 128       rhomvec[NjetAcc] = GetMd(jet) / jet->Area();
 129       fHistMdAreavsCent->Fill(fCent,rhomvec[NjetAcc]);
 130       Evec[NjetAcc] = jet->E();
 131       Mvec[NjetAcc] = jet->M();
 132       ++NjetAcc;
 133     }
 134   }
 135
 136   if (NjetAcc > 0) {
 137     //find median value
 138     Double_t rhom = TMath::Median(NjetAcc, rhomvec);
 139     fOutRhoMass->SetVal(rhom);
 140
 141     Int_t Ntracks = fTracks->GetEntries();
 142     Double_t meanM = TMath::Mean(NjetAcc, Mvec);
 143     Double_t meanE = TMath::Mean(NjetAcc, Evec);
 144     Double_t gamma = 0.;
 145     if(meanM>0.) gamma = meanE/meanM;
 146     fHistGammaVsNtrack->Fill(Ntracks,gamma);
 147
 148     if (fOutRhoMassScaled) {
 149       Double_t rhomScaled = rhom * GetScaleFactor(fCent);
 150       fOutRhoMassScaled->SetVal(rhomScaled);
 151     }
 152   }
 153
 154   return kTRUE;
 155 }
 156
 157 //________________________________________________________________________
 158 Double_t AliAnalysisTaskRhoMass::GetSumMConstituents(AliEmcalJet *jet) {
 159
 160   Double_t sum = 0.;
 161
 162   AliVParticle *vp;
 163   for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
 164     vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
 165     if(!vp) continue;
 166     sum+=vp->M();
 167   }
 168   return sum;
 169 }
 170
 171 //________________________________________________________________________
 172 Double_t AliAnalysisTaskRhoMass::GetSumPtConstituents(AliEmcalJet *jet) {
 173
 174   Double_t sum = 0.;
 175
 176   AliVParticle *vp;
 177   for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
 178     vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
 179     if(!vp) continue;
 180     sum+=vp->Pt();
 181   }
 182   return sum;
 183 }
 184
 185 //________________________________________________________________________
 186 Double_t AliAnalysisTaskRhoMass::GetMd(AliEmcalJet *jet) {
 187   //get md as defined in http://arxiv.org/pdf/1211.2811.pdf
 188   Double_t sum = 0.;
 189   Double_t px = 0.;
 190   Double_t py = 0.;
 191   Double_t pz = 0.;
 192   Double_t E = 0.;
 193   AliVParticle *vp;
 194   for(Int_t icc=0; icc<jet->GetNumberOfTracks(); icc++) {
 195     vp = static_cast<AliVParticle*>(jet->TrackAt(icc, fTracks));
 196     if(!vp) continue;
 197     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)
 198     else if(fJetRhoMassType==kMdP) sum += TMath::Sqrt(vp->M()*vp->M() + vp->P()*vp->P()) - vp->P();
 199     else if(fJetRhoMassType==kMd4) {
 200       px+=vp->Px();
 201       py+=vp->Py();
 202       pz+=vp->Pz();
 203       E+=vp->E();
 204     }
 205   }
 206   if(fJetRhoMassType==kMd4) {
 207     Double_t pt = TMath::Sqrt(px*px + py*py);
 208     Double_t m2 = E*E - pt*pt - pz*pz;
 209     sum = TMath::Sqrt(m2 + pt*pt) - pt;
 210   }
 211   return sum;
 212 }