PWGJE/EMCALJetTasks/AliAnalysisTaskEmcalJetSample.cxx

   1 // $Id$
   2 //
   3 // Jet sample analysis task.
   4 //
   5 // Author: S.Aiola
   6
   7 #include <TClonesArray.h>
   8 #include <TH1F.h>
   9 #include <TH2F.h>
  10 #include <TList.h>
  11 #include <TLorentzVector.h>
  12
  13 #include "AliVCluster.h"
  14 #include "AliAODCaloCluster.h"
  15 #include "AliVTrack.h"
  16 #include "AliEmcalJet.h"
  17 #include "AliRhoParameter.h"
  18 #include "AliLog.h"
  19 #include "AliJetContainer.h"
  20 #include "AliParticleContainer.h"
  21 #include "AliClusterContainer.h"
  22
  23 #include "AliAnalysisTaskEmcalJetSample.h"
  24
  25 ClassImp(AliAnalysisTaskEmcalJetSample)
  26
  27 //________________________________________________________________________
  28 AliAnalysisTaskEmcalJetSample::AliAnalysisTaskEmcalJetSample() :
  29   AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetSample", kTRUE),
  30   fHistTracksPt(0),
  31   fHistClustersPt(0),
  32   fHistLeadingJetPt(0),
  33   fHistJetsPhiEta(0),
  34   fHistJetsPtArea(0),
  35   fHistJetsPtLeadHad(0),
  36   fHistJetsCorrPtArea(0),
  37   fJetsCont(0),
  38   fTracksCont(0),
  39   fCaloClustersCont(0)
  40
  41 {
  42   // Default constructor.
  43
  44   fHistTracksPt       = new TH1*[fNcentBins];
  45   fHistClustersPt     = new TH1*[fNcentBins];
  46   fHistLeadingJetPt   = new TH1*[fNcentBins];
  47   fHistJetsPhiEta     = new TH2*[fNcentBins];
  48   fHistJetsPtArea     = new TH2*[fNcentBins];
  49   fHistJetsPtLeadHad  = new TH2*[fNcentBins];
  50   fHistJetsCorrPtArea = new TH2*[fNcentBins];
  51
  52   for (Int_t i = 0; i < fNcentBins; i++) {
  53     fHistTracksPt[i] = 0;
  54     fHistClustersPt[i] = 0;
  55     fHistLeadingJetPt[i] = 0;
  56     fHistJetsPhiEta[i] = 0;
  57     fHistJetsPtArea[i] = 0;
  58     fHistJetsPtLeadHad[i] = 0;
  59     fHistJetsCorrPtArea[i] = 0;
  60   }
  61
  62   SetMakeGeneralHistograms(kTRUE);
  63 }
  64
  65 //________________________________________________________________________
  66 AliAnalysisTaskEmcalJetSample::AliAnalysisTaskEmcalJetSample(const char *name) :
  67   AliAnalysisTaskEmcalJet(name, kTRUE),
  68   fHistTracksPt(0),
  69   fHistClustersPt(0),
  70   fHistLeadingJetPt(0),
  71   fHistJetsPhiEta(0),
  72   fHistJetsPtArea(0),
  73   fHistJetsPtLeadHad(0),
  74   fHistJetsCorrPtArea(0),
  75   fJetsCont(0),
  76   fTracksCont(0),
  77   fCaloClustersCont(0)
  78 {
  79   // Standard constructor.
  80
  81   fHistTracksPt       = new TH1*[fNcentBins];
  82   fHistClustersPt     = new TH1*[fNcentBins];
  83   fHistLeadingJetPt   = new TH1*[fNcentBins];
  84   fHistJetsPhiEta     = new TH2*[fNcentBins];
  85   fHistJetsPtArea     = new TH2*[fNcentBins];
  86   fHistJetsPtLeadHad  = new TH2*[fNcentBins];
  87   fHistJetsCorrPtArea = new TH2*[fNcentBins];
  88
  89   for (Int_t i = 0; i < fNcentBins; i++) {
  90     fHistTracksPt[i] = 0;
  91     fHistClustersPt[i] = 0;
  92     fHistLeadingJetPt[i] = 0;
  93     fHistJetsPhiEta[i] = 0;
  94     fHistJetsPtArea[i] = 0;
  95     fHistJetsPtLeadHad[i] = 0;
  96     fHistJetsCorrPtArea[i] = 0;
  97   }
  98
  99   SetMakeGeneralHistograms(kTRUE);
 100 }
 101
 102 //________________________________________________________________________
 103 AliAnalysisTaskEmcalJetSample::~AliAnalysisTaskEmcalJetSample()
 104 {
 105   // Destructor.
 106 }
 107
 108 //________________________________________________________________________
 109 void AliAnalysisTaskEmcalJetSample::UserCreateOutputObjects()
 110 {
 111   // Create user output.
 112
 113   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 114
 115   fJetsCont           = GetJetContainer(0);
 116   if(fJetsCont) { //get particles and clusters connected to jets
 117     fTracksCont       = fJetsCont->GetParticleContainer();
 118     fCaloClustersCont = fJetsCont->GetClusterContainer();
 119   } else {        //no jets, just analysis tracks and clusters
 120     fTracksCont       = GetParticleContainer(0);
 121     fCaloClustersCont = GetClusterContainer(0);
 122   }
 123   fTracksCont->SetClassName("AliVTrack");
 124   fCaloClustersCont->SetClassName("AliAODCaloCluster");
 125
 126   TString histname;
 127
 128   for (Int_t i = 0; i < fNcentBins; i++) {
 129     if (fParticleCollArray.GetEntriesFast()>0) {
 130       histname = "fHistTracksPt_";
 131       histname += i;
 132       fHistTracksPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
 133       fHistTracksPt[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
 134       fHistTracksPt[i]->GetYaxis()->SetTitle("counts");
 135       fOutput->Add(fHistTracksPt[i]);
 136     }
 137
 138     if (fClusterCollArray.GetEntriesFast()>0) {
 139       histname = "fHistClustersPt_";
 140       histname += i;
 141       fHistClustersPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
 142       fHistClustersPt[i]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
 143       fHistClustersPt[i]->GetYaxis()->SetTitle("counts");
 144       fOutput->Add(fHistClustersPt[i]);
 145     }
 146
 147     if (fJetCollArray.GetEntriesFast()>0) {
 148       histname = "fHistLeadingJetPt_";
 149       histname += i;
 150       fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 151       fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 152       fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
 153       fOutput->Add(fHistLeadingJetPt[i]);
 154
 155       histname = "fHistJetsPhiEta_";
 156       histname += i;
 157       fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 50, -1, 1, 101, 0, TMath::Pi()*2 + TMath::Pi()/200);
 158       fHistJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
 159       fHistJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
 160       fOutput->Add(fHistJetsPhiEta[i]);
 161
 162       histname = "fHistJetsPtArea_";
 163       histname += i;
 164       fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 30, 0, 3);
 165       fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 166       fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
 167       fOutput->Add(fHistJetsPtArea[i]);
 168
 169       histname = "fHistJetsPtLeadHad_";
 170       histname += i;
 171       fHistJetsPtLeadHad[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
 172       fHistJetsPtLeadHad[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 173       fHistJetsPtLeadHad[i]->GetYaxis()->SetTitle("p_{T,lead} (GeV/c)");
 174       fHistJetsPtLeadHad[i]->GetZaxis()->SetTitle("counts");
 175       fOutput->Add(fHistJetsPtLeadHad[i]);
 176
 177       if (!(GetJetContainer()->GetRhoName().IsNull())) {
 178         histname = "fHistJetsCorrPtArea_";
 179         histname += i;
 180         fHistJetsCorrPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins*2, -fMaxBinPt, fMaxBinPt, 30, 0, 3);
 181         fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
 182         fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
 183         fOutput->Add(fHistJetsCorrPtArea[i]);
 184       }
 185     }
 186   }
 187   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 188 }
 189
 190 //________________________________________________________________________
 191 Bool_t AliAnalysisTaskEmcalJetSample::FillHistograms()
 192 {
 193   // Fill histograms.
 194
 195   if (fTracksCont) {
 196     AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0));
 197     while(track) {
 198       fHistTracksPt[fCentBin]->Fill(track->Pt());
 199       Int_t emc1 = track->GetEMCALcluster();
 200       Printf("EMCAL cluster %d",emc1);
 201       track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
 202     }
 203   }
 204
 205   if (fCaloClustersCont) {
 206     //    AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0);
 207     AliAODCaloCluster *cluster = static_cast<AliAODCaloCluster*>(fCaloClustersCont->GetNextAcceptCluster(0));
 208     while(cluster) {
 209
 210       TLorentzVector nPart;
 211       cluster->GetMomentum(nPart, fVertex);
 212       fHistClustersPt[fCentBin]->Fill(nPart.Pt());
 213
 214       AliVTrack *mt = NULL;
 215       Printf("N matched tracks: %d",cluster->GetNTracksMatched());
 216       if(cluster->GetNTracksMatched()>1)
 217         mt = static_cast<AliVTrack*>(cluster->GetTrackMatched(0));
 218       if(mt) Printf("matched track pt: %f eta: %f phi: %f",mt->Pt(),mt->Eta(),mt->Phi());
 219       cluster = static_cast<AliAODCaloCluster*>(fCaloClustersCont->GetNextAcceptCluster());
 220     }
 221   }
 222
 223   if (fJetsCont) {
 224     AliEmcalJet *jet = fJetsCont->GetNextAcceptJet(0);
 225     while(jet) {
 226
 227       fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
 228       fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
 229
 230       Float_t ptLeading = fJetsCont->GetLeadingHadronPt(jet);
 231       fHistJetsPtLeadHad[fCentBin]->Fill(jet->Pt(), ptLeading);
 232
 233       Float_t corrPt = jet->Pt() - fJetsCont->GetRhoVal() * jet->Area();
 234       fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());
 235
 236       jet = fJetsCont->GetNextAcceptJet();
 237     }
 238
 239     jet = fJetsCont->GetLeadingJet();
 240     if(jet) fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
 241   }
 242
 243   return kTRUE;
 244 }
 245
 246 //________________________________________________________________________
 247 void AliAnalysisTaskEmcalJetSample::ExecOnce() {
 248
 249   AliAnalysisTaskEmcalJet::ExecOnce();
 250
 251   if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
 252   if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
 253   if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
 254
 255 }
 256
 257 //________________________________________________________________________
 258 Bool_t AliAnalysisTaskEmcalJetSample::Run()
 259 {
 260   // Run analysis code here, if needed. It will be executed before FillHistograms().
 261
 262   return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
 263 }
 264
 265 //________________________________________________________________________
 266 void AliAnalysisTaskEmcalJetSample::Terminate(Option_t *)
 267 {
 268   // Called once at the end of the analysis.
 269 }