PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoMassBase.cxx

   1 // $Id$
   2 //
   3 // Base class for rho mass calculation.
   4 // Calculates parameterized rho mass for given centrality independent of input.
   5 //
   6 // Author: M. Verweij. Similar to AliAnalysisTaskRhoBase
   7
   8 #include <TF1.h>
   9 #include <TH1F.h>
  10 #include <TH2F.h>
  11 #include <TClonesArray.h>
  12
  13 #include "AliLog.h"
  14 #include "AliRhoParameter.h"
  15 #include "AliEmcalJet.h"
  16
  17 #include "AliAnalysisTaskRhoMassBase.h"
  18
  19 ClassImp(AliAnalysisTaskRhoMassBase)
  20
  21 //________________________________________________________________________
  22 AliAnalysisTaskRhoMassBase::AliAnalysisTaskRhoMassBase() :
  23   AliAnalysisTaskEmcalJet("AliAnalysisTaskRhoMassBase", kFALSE),
  24   fOutRhoMassName(),
  25   fOutRhoMassScaledName(),
  26   fCompareRhoMassName(),
  27   fCompareRhoMassScaledName(),
  28   fRhoMassFunction(0),
  29   fScaleFunction(0),
  30   fAttachToEvent(kTRUE),
  31   fOutRhoMass(0),
  32   fOutRhoMassScaled(0),
  33   fCompareRhoMass(0),
  34   fCompareRhoMassScaled(0),
  35   fHistJetMassvsCent(0),
  36   fHistJetRhoMassvsCent(0),
  37   fHistRhoMassvsCent(0),
  38   fHistRhoMassScaledvsCent(0),
  39   fHistDeltaRhoMassvsCent(0),
  40   fHistDeltaRhoMassScalevsCent(0),
  41   fHistRhoMassvsNtrack(0),
  42   fHistRhoMassScaledvsNtrack(0),
  43   fHistDeltaRhoMassvsNtrack(0),
  44   fHistDeltaRhoMassScalevsNtrack(0),
  45   fHistRhoMassvsNcluster(0),
  46   fHistRhoMassScaledvsNcluster(0),
  47   fHistGammaVsNtrack(0)
  48 {
  49   // Constructor.
  50
  51   for (Int_t i = 0; i < 4; i++)
  52     fHistJetRhoMassvsEta[i] = 0;
  53 }
  54
  55 //________________________________________________________________________
  56 AliAnalysisTaskRhoMassBase::AliAnalysisTaskRhoMassBase(const char *name, Bool_t histo) :
  57   AliAnalysisTaskEmcalJet(name, histo),
  58   fOutRhoMassName(),
  59   fOutRhoMassScaledName(),
  60   fCompareRhoMassName(),
  61   fCompareRhoMassScaledName(),
  62   fRhoMassFunction(0),
  63   fScaleFunction(0),
  64   fAttachToEvent(kTRUE),
  65   fOutRhoMass(0),
  66   fOutRhoMassScaled(0),
  67   fCompareRhoMass(0),
  68   fCompareRhoMassScaled(0),
  69   fHistJetMassvsCent(0),
  70   fHistJetRhoMassvsCent(0),
  71   fHistRhoMassvsCent(0),
  72   fHistRhoMassScaledvsCent(0),
  73   fHistDeltaRhoMassvsCent(0),
  74   fHistDeltaRhoMassScalevsCent(0),
  75   fHistRhoMassvsNtrack(0),
  76   fHistRhoMassScaledvsNtrack(0),
  77   fHistDeltaRhoMassvsNtrack(0),
  78   fHistDeltaRhoMassScalevsNtrack(0),
  79   fHistRhoMassvsNcluster(0),
  80   fHistRhoMassScaledvsNcluster(0),
  81   fHistGammaVsNtrack(0)
  82 {
  83   // Constructor.
  84
  85   for (Int_t i = 0; i < 4; i++)
  86     fHistJetRhoMassvsEta[i] = 0;
  87
  88   SetMakeGeneralHistograms(histo);
  89 }
  90
  91 //________________________________________________________________________
  92 void AliAnalysisTaskRhoMassBase::UserCreateOutputObjects()
  93 {
  94   // User create output objects, called at the beginning of the analysis.
  95
  96   if (!fCreateHisto)
  97     return;
  98
  99   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 100
 101   fHistRhoMassvsCent = new TH2F("fHistRhoMassvsCent", "fHistRhoMassvsCent", 101, -1,  100, fNbins, fMinBinPt, fMaxBinPt/2.);
 102   fHistRhoMassvsCent->GetXaxis()->SetTitle("Centrality (%)");
 103   fHistRhoMassvsCent->GetYaxis()->SetTitle("#rho_{m} (GeV/c * rad^{-1})");
 104   fOutput->Add(fHistRhoMassvsCent);
 105
 106   if (fParticleCollArray.GetEntriesFast()>0) {
 107     fHistRhoMassvsNtrack = new TH2F("fHistRhoMassvsNtrack", "fHistRhoMassvsNtrack", 150, 0, 6000, fNbins, fMinBinPt, fMaxBinPt/2);
 108     fHistRhoMassvsNtrack->GetXaxis()->SetTitle("No. of tracks");
 109     fHistRhoMassvsNtrack->GetYaxis()->SetTitle("#rho_{m} (GeV/c * rad^{-1})");
 110     fOutput->Add(fHistRhoMassvsNtrack);
 111
 112     //fHistGammaVsNtrack
 113     fHistGammaVsNtrack = new TH2F("fHistGammaVsNtrack", "fHistGammaVsNtrack", 150, 0, 6000, 100,0.,10.);
 114     fHistGammaVsNtrack->GetXaxis()->SetTitle("No. of tracks");
 115     fHistGammaVsNtrack->GetYaxis()->SetTitle("#gamma = #LT E #GT / #LT M #GT");
 116     fOutput->Add(fHistGammaVsNtrack);
 117   }
 118
 119   if (fClusterCollArray.GetEntriesFast()>0) {
 120     fHistRhoMassvsNcluster = new TH2F("fHistRhoMassvsNcluster", "fHistRhoMassvsNcluster", 50, 0, 1500, fNbins, fMinBinPt, fMaxBinPt/2);
 121     fHistRhoMassvsNcluster->GetXaxis()->SetTitle("No. of tracks");
 122     fHistRhoMassvsNcluster->GetYaxis()->SetTitle("#rho_{m} (GeV/c * rad^{-1})");
 123     fOutput->Add(fHistRhoMassvsNcluster);
 124   }
 125
 126   if (fJetCollArray.GetEntriesFast()>0) {
 127     fHistJetMassvsCent = new TH2F("fHistJetMassvsCent", "fHistJetMassvsCent", 101, -1,  100, fNbins, fMinBinPt, fMaxBinPt);
 128     fHistJetMassvsCent->GetXaxis()->SetTitle("Centrality (%)");
 129     fHistJetMassvsCent->GetYaxis()->SetTitle("#it{M}_{jet} (GeV/c)");
 130     fOutput->Add(fHistJetMassvsCent);
 131
 132     fHistJetRhoMassvsCent = new TH2F("fHistJetRhoMassvsCent", "fHistJetRhoMassvsCent", 101, -1, 100, fNbins, fMinBinPt, fMaxBinPt*2);
 133     fHistJetRhoMassvsCent->GetXaxis()->SetTitle("Centrality (%)");
 134     fHistJetRhoMassvsCent->GetYaxis()->SetTitle("Jet #rho_{m} (GeV/c)");
 135     fOutput->Add(fHistJetRhoMassvsCent);
 136
 137
 138     TString name;
 139     for (Int_t i = 0; i < 4; i++) {
 140       name = Form("fHistJetRhoMassvsEta_%d",i);
 141       fHistJetRhoMassvsEta[i] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt*2, 16, -0.8, 0.8);
 142       fHistJetRhoMassvsEta[i]->GetXaxis()->SetTitle("#rho_{m} (GeV/c)");
 143       fHistJetRhoMassvsEta[i]->GetYaxis()->SetTitle("#eta");
 144       fOutput->Add(fHistJetRhoMassvsEta[i]);
 145     }
 146   }
 147
 148   if (!fCompareRhoMassName.IsNull()) {
 149     fHistDeltaRhoMassvsCent = new TH2F("fHistDeltaRhoMassvsCent", "fHistDeltaRhoMassvsCent", 101, -1, 100, fNbins, -fMaxBinPt, fMaxBinPt);
 150     fHistDeltaRhoMassvsCent->GetXaxis()->SetTitle("Centrality (%)");
 151     fHistDeltaRhoMassvsCent->GetYaxis()->SetTitle("#Delta#rho (GeV/c * rad^{-1})");
 152     fOutput->Add(fHistDeltaRhoMassvsCent);
 153
 154     if (fParticleCollArray.GetEntriesFast()>0) {
 155       fHistDeltaRhoMassvsNtrack = new TH2F("fHistDeltaRhoMassvsNtrack", "fHistDeltaRhoMassvsNtrack", 150, 0, 6000, fNbins, -fMaxBinPt, fMaxBinPt);
 156       fHistDeltaRhoMassvsNtrack->GetXaxis()->SetTitle("No. of tracks");
 157       fHistDeltaRhoMassvsNtrack->GetYaxis()->SetTitle("#Delta#rho (GeV/c * rad^{-1})");
 158       fOutput->Add(fHistDeltaRhoMassvsNtrack);
 159     }
 160   }
 161
 162   if (fScaleFunction) {
 163     fHistRhoMassScaledvsCent = new TH2F("fHistRhoMassScaledvsCent", "fHistRhoMassScaledvsCent", 101, -1, 100, fNbins, fMinBinPt , fMaxBinPt*2);
 164     fHistRhoMassScaledvsCent->GetXaxis()->SetTitle("Centrality (%)");
 165     fHistRhoMassScaledvsCent->GetYaxis()->SetTitle("#rho_{m,scaled} (GeV/c * rad^{-1})");
 166     fOutput->Add(fHistRhoMassScaledvsCent);
 167
 168     if (fParticleCollArray.GetEntriesFast()>0) {
 169       fHistRhoMassScaledvsNtrack = new TH2F("fHistRhoMassScaledvsNtrack", "fHistRhoMassScaledvsNtrack", 150, 0, 6000, fNbins, fMinBinPt, fMaxBinPt*2);
 170       fHistRhoMassScaledvsNtrack->GetXaxis()->SetTitle("No. of tracks");
 171       fHistRhoMassScaledvsNtrack->GetYaxis()->SetTitle("#rho_{m,scaled} (GeV/c * rad^{-1})");
 172       fOutput->Add(fHistRhoMassScaledvsNtrack);
 173     }
 174
 175     if (fClusterCollArray.GetEntriesFast()>0) {
 176       fHistRhoMassScaledvsNcluster = new TH2F("fHistRhoMassScaledvsNcluster", "fHistRhoMassScaledvsNcluster", 50, 0, 1500, fNbins, fMinBinPt, fMaxBinPt*2);
 177       fHistRhoMassScaledvsNcluster->GetXaxis()->SetTitle("No. of clusters");
 178       fHistRhoMassScaledvsNcluster->GetYaxis()->SetTitle("#rho_{m,scaled} (GeV/c * rad^{-1})");
 179       fOutput->Add(fHistRhoMassScaledvsNcluster);
 180     }
 181
 182     if (!fCompareRhoMassScaledName.IsNull()) {
 183       fHistDeltaRhoMassScalevsCent = new TH2F("fHistDeltaRhoMassScalevsCent", "fHistDeltaRhoMassScalevsCent", 101, -1, 100, fNbins, -fMaxBinPt, fMaxBinPt);
 184       fHistDeltaRhoMassScalevsCent->GetXaxis()->SetTitle("Centrality (%)");
 185       fHistDeltaRhoMassScalevsCent->GetYaxis()->SetTitle("#Delta#rho_{m,scaled} (GeV/c * rad^{-1})");
 186       fOutput->Add(fHistDeltaRhoMassScalevsCent);
 187
 188       if (fParticleCollArray.GetEntriesFast()>0) {
 189         fHistDeltaRhoMassScalevsNtrack = new TH2F("fHistDeltaRhoMassScalevsNtrack", "fHistDeltaRhoMassScalevsNtrack", 150, 0, 6000, fNbins, -fMaxBinPt, fMaxBinPt);
 190         fHistDeltaRhoMassScalevsNtrack->GetXaxis()->SetTitle("No. of tracks");
 191         fHistDeltaRhoMassScalevsNtrack->GetYaxis()->SetTitle("#Delta#rho_{m,scaled} (GeV/c * rad^{-1})");
 192         fOutput->Add(fHistDeltaRhoMassScalevsNtrack);
 193       }
 194     }
 195   }
 196 }
 197
 198 //________________________________________________________________________
 199 Bool_t AliAnalysisTaskRhoMassBase::Run()
 200 {
 201   // Run the analysis.
 202
 203   Double_t rhom = GetRhoMassFactor(fCent);
 204   fOutRhoMass->SetVal(rhom);
 205
 206   if (fScaleFunction) {
 207     Double_t rhomScaled = rhom * GetScaleFactor(fCent);
 208     fOutRhoMassScaled->SetVal(rhomScaled);
 209   }
 210
 211   return kTRUE;
 212 }
 213
 214 //________________________________________________________________________
 215 Bool_t AliAnalysisTaskRhoMassBase::FillHistograms()
 216 {
 217   // Fill histograms.
 218
 219   Int_t Ntracks   = 0;
 220   Int_t Nclusters = 0;
 221
 222   if (fTracks)
 223     Ntracks = fTracks->GetEntries();
 224   if (fCaloClusters)
 225     Nclusters = fCaloClusters->GetEntries();
 226
 227   if (fJets) {
 228     Int_t    Njets         = fJets->GetEntries();
 229     Int_t    NjetAcc       = 0;
 230
 231     for (Int_t i = 0; i < Njets; ++i) {
 232       AliEmcalJet *jet = static_cast<AliEmcalJet*>(fJets->At(i));
 233       if (!jet) {
 234         AliError(Form("%s: Could not receive jet %d", GetName(), i));
 235         continue;
 236       }
 237
 238       if (!AcceptJet(jet))
 239         continue;
 240
 241       fHistJetMassvsCent->Fill(fCent, jet->M());
 242       fHistJetRhoMassvsCent->Fill(fCent, jet->M() / jet->Area());
 243       fHistJetRhoMassvsEta[fCentBin]->Fill(jet->M() / jet->Area(), jet->Eta());
 244
 245       NjetAcc++;
 246     }
 247   }
 248
 249   fHistRhoMassvsCent->Fill(fCent, fOutRhoMass->GetVal());
 250
 251   if (fTracks)
 252     fHistRhoMassvsNtrack->Fill(Ntracks, fOutRhoMass->GetVal());
 253   if (fCaloClusters)
 254     fHistRhoMassvsNcluster->Fill(Nclusters, fOutRhoMass->GetVal());
 255   if (fCompareRhoMass) {
 256     fHistDeltaRhoMassvsCent->Fill(fCent, fOutRhoMass->GetVal() - fCompareRhoMass->GetVal());
 257     if (fTracks)
 258       fHistDeltaRhoMassvsNtrack->Fill(Ntracks, fOutRhoMass->GetVal() - fCompareRhoMass->GetVal());
 259   }
 260
 261   if (fOutRhoMassScaled) {
 262     fHistRhoMassScaledvsCent->Fill(fCent, fOutRhoMassScaled->GetVal());
 263     if (fTracks)
 264       fHistRhoMassScaledvsNtrack->Fill(Ntracks, fOutRhoMassScaled->GetVal());
 265     if (fCaloClusters)
 266       fHistRhoMassScaledvsNcluster->Fill(Nclusters,  fOutRhoMassScaled->GetVal());
 267     if (fCompareRhoMassScaled) {
 268       fHistDeltaRhoMassScalevsCent->Fill(fCent, fOutRhoMassScaled->GetVal() - fCompareRhoMassScaled->GetVal());
 269       if (fTracks)
 270         fHistDeltaRhoMassScalevsNtrack->Fill(Ntracks, fOutRhoMassScaled->GetVal() - fCompareRhoMassScaled->GetVal());
 271     }
 272   }
 273
 274   return kTRUE;
 275 }
 276
 277
 278 //________________________________________________________________________
 279 void AliAnalysisTaskRhoMassBase::ExecOnce()
 280 {
 281   // Init the analysis.
 282
 283   if (!fOutRhoMass) {
 284     fOutRhoMass = new AliRhoParameter(fOutRhoMassName, 0);
 285
 286     if (fAttachToEvent) {
 287       if (!(InputEvent()->FindListObject(fOutRhoMassName))) {
 288         InputEvent()->AddObject(fOutRhoMass);
 289       } else {
 290         AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fOutRhoMassName.Data()));
 291         return;
 292       }
 293     }
 294   }
 295
 296   if (fScaleFunction && !fOutRhoMassScaled) {
 297     fOutRhoMassScaled = new AliRhoParameter(fOutRhoMassScaledName, 0);
 298
 299     if (fAttachToEvent) {
 300       if (!(InputEvent()->FindListObject(fOutRhoMassScaledName))) {
 301         InputEvent()->AddObject(fOutRhoMassScaled);
 302       } else {
 303         AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fOutRhoMassScaledName.Data()));
 304         return;
 305       }
 306     }
 307   }
 308
 309   if (!fCompareRhoMassName.IsNull() && !fCompareRhoMass) {
 310     fCompareRhoMass = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fCompareRhoMassName));
 311     if (!fCompareRhoMass) {
 312       AliWarning(Form("%s: Could not retrieve rho %s!", GetName(), fCompareRhoMassName.Data()));
 313     }
 314   }
 315
 316   if (!fCompareRhoMassScaledName.IsNull() && !fCompareRhoMassScaled) {
 317     fCompareRhoMassScaled = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fCompareRhoMassScaledName));
 318     if (!fCompareRhoMassScaled) {
 319       AliWarning(Form("%s: Could not retrieve rho %s!", GetName(), fCompareRhoMassScaledName.Data()));
 320     }
 321   }
 322
 323   AliAnalysisTaskEmcalJet::ExecOnce();
 324 }
 325
 326 //________________________________________________________________________
 327 Double_t AliAnalysisTaskRhoMassBase::GetRhoMassFactor(Double_t cent)
 328 {
 329   // Return rho per centrality.
 330
 331   Double_t rho = 0;
 332   if (fRhoMassFunction)
 333     rho = fRhoMassFunction->Eval(cent);
 334   return rho;
 335 }
 336
 337 //________________________________________________________________________
 338 Double_t AliAnalysisTaskRhoMassBase::GetScaleFactor(Double_t cent)
 339 {
 340   // Get scale factor.
 341
 342   Double_t scale = 1;
 343   if (fScaleFunction)
 344     scale = fScaleFunction->Eval(cent);
 345   return scale;
 346 }