PWGJE/EMCALJetTasks/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   fPionMassClusters(kFALSE),
  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),
  37   fPionMassClusters(kFALSE),
  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(sumM*sumM + sumPt*sumPt) - 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.;
 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       }
 209     }
 210   }
 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);
 219       Double_t m = 0.;
 220       if(fPionMassClusters) m = 0.13957;
 221       if(fJetRhoMassType==kMd) sum += TMath::Sqrt(m*m + nPart.Pt()*nPart.Pt()) - nPart.Pt();
 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
 232   if(fJetRhoMassType==kMd4) {
 233     Double_t pt = TMath::Sqrt(px*px + py*py);
 234     Double_t m2 = E*E - pt*pt - pz*pz;
 235     sum = TMath::Sqrt(m2 + pt*pt) - pt;
 236   }
 237   return sum;
 238 }