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 <TH3F.h>
  11 #include <TList.h>
  12 #include <TLorentzVector.h>
  13
  14 #include "AliVCluster.h"
  15 #include "AliAODCaloCluster.h"
  16 #include "AliESDCaloCluster.h"
  17 #include "AliVTrack.h"
  18 #include "AliEmcalJet.h"
  19 #include "AliRhoParameter.h"
  20 #include "AliLog.h"
  21 #include "AliJetContainer.h"
  22 #include "AliParticleContainer.h"
  23 #include "AliClusterContainer.h"
  24 #include "AliPicoTrack.h"
  25
  26 #include "AliAnalysisTaskEmcalJetSample.h"
  27
  28 ClassImp(AliAnalysisTaskEmcalJetSample)
  29
  30 //________________________________________________________________________
  31 AliAnalysisTaskEmcalJetSample::AliAnalysisTaskEmcalJetSample() :
  32   AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetSample", kTRUE),
  33   fHistTracksPt(0),
  34   fHistClustersPt(0),
  35   fHistLeadingJetPt(0),
  36   fHistJetsPhiEta(0),
  37   fHistJetsPtArea(0),
  38   fHistJetsPtLeadHad(0),
  39   fHistJetsCorrPtArea(0),
  40   fHistPtDEtaDPhiTrackClus(0),
  41   fHistPtDEtaDPhiClusTrack(0),
  42   fJetsCont(0),
  43   fTracksCont(0),
  44   fCaloClustersCont(0)
  45
  46 {
  47   // Default constructor.
  48
  49   fHistTracksPt       = new TH1*[fNcentBins];
  50   fHistClustersPt     = new TH1*[fNcentBins];
  51   fHistLeadingJetPt   = new TH1*[fNcentBins];
  52   fHistJetsPhiEta     = new TH2*[fNcentBins];
  53   fHistJetsPtArea     = new TH2*[fNcentBins];
  54   fHistJetsPtLeadHad  = new TH2*[fNcentBins];
  55   fHistJetsCorrPtArea = new TH2*[fNcentBins];
  56
  57   for (Int_t i = 0; i < fNcentBins; i++) {
  58     fHistTracksPt[i] = 0;
  59     fHistClustersPt[i] = 0;
  60     fHistLeadingJetPt[i] = 0;
  61     fHistJetsPhiEta[i] = 0;
  62     fHistJetsPtArea[i] = 0;
  63     fHistJetsPtLeadHad[i] = 0;
  64     fHistJetsCorrPtArea[i] = 0;
  65   }
  66
  67   SetMakeGeneralHistograms(kTRUE);
  68 }
  69
  70 //________________________________________________________________________
  71 AliAnalysisTaskEmcalJetSample::AliAnalysisTaskEmcalJetSample(const char *name) :
  72   AliAnalysisTaskEmcalJet(name, kTRUE),
  73   fHistTracksPt(0),
  74   fHistClustersPt(0),
  75   fHistLeadingJetPt(0),
  76   fHistJetsPhiEta(0),
  77   fHistJetsPtArea(0),
  78   fHistJetsPtLeadHad(0),
  79   fHistJetsCorrPtArea(0),
  80   fHistPtDEtaDPhiTrackClus(0),
  81   fHistPtDEtaDPhiClusTrack(0),
  82   fJetsCont(0),
  83   fTracksCont(0),
  84   fCaloClustersCont(0)
  85 {
  86   // Standard constructor.
  87
  88   fHistTracksPt       = new TH1*[fNcentBins];
  89   fHistClustersPt     = new TH1*[fNcentBins];
  90   fHistLeadingJetPt   = new TH1*[fNcentBins];
  91   fHistJetsPhiEta     = new TH2*[fNcentBins];
  92   fHistJetsPtArea     = new TH2*[fNcentBins];
  93   fHistJetsPtLeadHad  = new TH2*[fNcentBins];
  94   fHistJetsCorrPtArea = new TH2*[fNcentBins];
  95
  96   for (Int_t i = 0; i < fNcentBins; i++) {
  97     fHistTracksPt[i] = 0;
  98     fHistClustersPt[i] = 0;
  99     fHistLeadingJetPt[i] = 0;
 100     fHistJetsPhiEta[i] = 0;
 101     fHistJetsPtArea[i] = 0;
 102     fHistJetsPtLeadHad[i] = 0;
 103     fHistJetsCorrPtArea[i] = 0;
 104   }
 105
 106   SetMakeGeneralHistograms(kTRUE);
 107 }
 108
 109 //________________________________________________________________________
 110 AliAnalysisTaskEmcalJetSample::~AliAnalysisTaskEmcalJetSample()
 111 {
 112   // Destructor.
 113 }
 114
 115 //________________________________________________________________________
 116 void AliAnalysisTaskEmcalJetSample::UserCreateOutputObjects()
 117 {
 118   // Create user output.
 119
 120   AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
 121
 122   fJetsCont           = GetJetContainer(0);
 123   if(fJetsCont) { //get particles and clusters connected to jets
 124     fTracksCont       = fJetsCont->GetParticleContainer();
 125     fCaloClustersCont = fJetsCont->GetClusterContainer();
 126   } else {        //no jets, just analysis tracks and clusters
 127     fTracksCont       = GetParticleContainer(0);
 128     fCaloClustersCont = GetClusterContainer(0);
 129   }
 130   fTracksCont->SetClassName("AliVTrack");
 131   fCaloClustersCont->SetClassName("AliAODCaloCluster");
 132
 133   TString histname;
 134
 135   for (Int_t i = 0; i < fNcentBins; i++) {
 136     if (fParticleCollArray.GetEntriesFast()>0) {
 137       histname = "fHistTracksPt_";
 138       histname += i;
 139       fHistTracksPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
 140       fHistTracksPt[i]->GetXaxis()->SetTitle("p_{T,track} (GeV/c)");
 141       fHistTracksPt[i]->GetYaxis()->SetTitle("counts");
 142       fOutput->Add(fHistTracksPt[i]);
 143     }
 144
 145     if (fClusterCollArray.GetEntriesFast()>0) {
 146       histname = "fHistClustersPt_";
 147       histname += i;
 148       fHistClustersPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins / 2, fMinBinPt, fMaxBinPt / 2);
 149       fHistClustersPt[i]->GetXaxis()->SetTitle("p_{T,clus} (GeV/c)");
 150       fHistClustersPt[i]->GetYaxis()->SetTitle("counts");
 151       fOutput->Add(fHistClustersPt[i]);
 152     }
 153
 154     if (fJetCollArray.GetEntriesFast()>0) {
 155       histname = "fHistLeadingJetPt_";
 156       histname += i;
 157       fHistLeadingJetPt[i] = new TH1F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt);
 158       fHistLeadingJetPt[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 159       fHistLeadingJetPt[i]->GetYaxis()->SetTitle("counts");
 160       fOutput->Add(fHistLeadingJetPt[i]);
 161
 162       histname = "fHistJetsPhiEta_";
 163       histname += i;
 164       fHistJetsPhiEta[i] = new TH2F(histname.Data(), histname.Data(), 50, -1, 1, 101, 0, TMath::Pi()*2 + TMath::Pi()/200);
 165       fHistJetsPhiEta[i]->GetXaxis()->SetTitle("#eta");
 166       fHistJetsPhiEta[i]->GetYaxis()->SetTitle("#phi");
 167       fOutput->Add(fHistJetsPhiEta[i]);
 168
 169       histname = "fHistJetsPtArea_";
 170       histname += i;
 171       fHistJetsPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, 30, 0, 3);
 172       fHistJetsPtArea[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 173       fHistJetsPtArea[i]->GetYaxis()->SetTitle("area");
 174       fOutput->Add(fHistJetsPtArea[i]);
 175
 176       histname = "fHistJetsPtLeadHad_";
 177       histname += i;
 178       fHistJetsPtLeadHad[i] = new TH2F(histname.Data(), histname.Data(), fNbins, fMinBinPt, fMaxBinPt, fNbins / 2, fMinBinPt, fMaxBinPt / 2);
 179       fHistJetsPtLeadHad[i]->GetXaxis()->SetTitle("p_{T}^{raw} (GeV/c)");
 180       fHistJetsPtLeadHad[i]->GetYaxis()->SetTitle("p_{T,lead} (GeV/c)");
 181       fHistJetsPtLeadHad[i]->GetZaxis()->SetTitle("counts");
 182       fOutput->Add(fHistJetsPtLeadHad[i]);
 183
 184       if (!(GetJetContainer()->GetRhoName().IsNull())) {
 185         histname = "fHistJetsCorrPtArea_";
 186         histname += i;
 187         fHistJetsCorrPtArea[i] = new TH2F(histname.Data(), histname.Data(), fNbins*2, -fMaxBinPt, fMaxBinPt, 30, 0, 3);
 188         fHistJetsCorrPtArea[i]->GetXaxis()->SetTitle("p_{T}^{corr} [GeV/c]");
 189         fHistJetsCorrPtArea[i]->GetYaxis()->SetTitle("area");
 190         fOutput->Add(fHistJetsCorrPtArea[i]);
 191       }
 192     }
 193   }
 194
 195   histname = "fHistPtDEtaDPhiTrackClus";
 196   fHistPtDEtaDPhiTrackClus = new TH3F(histname.Data(),Form("%s;#it{p}_{T}^{track};#Delta#eta;#Delta#varphi",histname.Data()),100,0.,100.,100,-0.1,0.1,100,-0.1,0.1);
 197   fOutput->Add(fHistPtDEtaDPhiTrackClus);
 198
 199   histname = "fHistPtDEtaDPhiClusTrack";
 200   fHistPtDEtaDPhiClusTrack = new TH3F(histname.Data(),Form("%s;#it{p}_{T}^{clus};#Delta#eta;#Delta#varphi",histname.Data()),100,0.,100.,100,-0.1,0.1,100,-0.1,0.1);
 201   fOutput->Add(fHistPtDEtaDPhiClusTrack);
 202
 203   PostData(1, fOutput); // Post data for ALL output slots > 0 here.
 204 }
 205
 206 //________________________________________________________________________
 207 Bool_t AliAnalysisTaskEmcalJetSample::FillHistograms()
 208 {
 209   // Fill histograms.
 210
 211   if (fTracksCont) {
 212     AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0));
 213     while(track) {
 214       fHistTracksPt[fCentBin]->Fill(track->Pt());
 215       track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
 216     }
 217   }
 218
 219   if (fCaloClustersCont) {
 220     AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0);
 221     while(cluster) {
 222       TLorentzVector nPart;
 223       cluster->GetMomentum(nPart, fVertex);
 224       fHistClustersPt[fCentBin]->Fill(nPart.Pt());
 225
 226       cluster = fCaloClustersCont->GetNextAcceptCluster();
 227     }
 228   }
 229
 230   if (fJetsCont) {
 231     AliEmcalJet *jet = fJetsCont->GetNextAcceptJet(0);
 232     while(jet) {
 233
 234       fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
 235       fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
 236
 237       Float_t ptLeading = fJetsCont->GetLeadingHadronPt(jet);
 238       fHistJetsPtLeadHad[fCentBin]->Fill(jet->Pt(), ptLeading);
 239
 240       if (fHistJetsCorrPtArea[fCentBin]) {
 241         Float_t corrPt = jet->Pt() - fJetsCont->GetRhoVal() * jet->Area();
 242         fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());
 243       }
 244       jet = fJetsCont->GetNextAcceptJet();
 245     }
 246
 247     jet = fJetsCont->GetLeadingJet();
 248     if(jet) fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
 249   }
 250
 251   CheckClusTrackMatching();
 252
 253   return kTRUE;
 254 }
 255
 256 //________________________________________________________________________
 257 void AliAnalysisTaskEmcalJetSample::CheckClusTrackMatching()
 258 {
 259
 260   if(!fTracksCont || !fCaloClustersCont)
 261     return;
 262
 263   Double_t deta = 999;
 264   Double_t dphi = 999;
 265
 266   //Get closest cluster to track
 267   AliVTrack *track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle(0));
 268   while(track) {
 269     //Get matched cluster
 270     Int_t emc1 = track->GetEMCALcluster();
 271     if(fCaloClustersCont && emc1>=0) {
 272       AliVCluster *clusMatch = fCaloClustersCont->GetCluster(emc1);
 273       if(clusMatch) {
 274         AliPicoTrack::GetEtaPhiDiff(track, clusMatch, dphi, deta);
 275         fHistPtDEtaDPhiTrackClus->Fill(track->Pt(),deta,dphi);
 276       }
 277     }
 278     track = static_cast<AliVTrack*>(fTracksCont->GetNextAcceptParticle());
 279   }
 280
 281   //Get closest track to cluster
 282   AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0);
 283   while(cluster) {
 284     TLorentzVector nPart;
 285     cluster->GetMomentum(nPart, fVertex);
 286     fHistClustersPt[fCentBin]->Fill(nPart.Pt());
 287
 288     //Get matched track
 289     AliVTrack *mt = NULL;
 290     AliAODCaloCluster *acl = dynamic_cast<AliAODCaloCluster*>(cluster);
 291     if(acl) {
 292       if(acl->GetNTracksMatched()>1)
 293         mt = static_cast<AliVTrack*>(acl->GetTrackMatched(0));
 294     }
 295     else {
 296       AliESDCaloCluster *ecl = dynamic_cast<AliESDCaloCluster*>(cluster);
 297       Int_t im = ecl->GetTrackMatchedIndex();
 298       if(fTracksCont && im>=0) {
 299         mt = static_cast<AliVTrack*>(fTracksCont->GetParticle(im));
 300       }
 301     }
 302     if(mt) {
 303       AliPicoTrack::GetEtaPhiDiff(mt, cluster, dphi, deta);
 304       fHistPtDEtaDPhiClusTrack->Fill(nPart.Pt(),deta,dphi);
 305
 306       /* //debugging
 307          if(mt->IsEMCAL()) {
 308          Int_t emc1 = mt->GetEMCALcluster();
 309          Printf("current id: %d  emc1: %d",fCaloClustersCont->GetCurrentID(),emc1);
 310          AliVCluster *clm = fCaloClustersCont->GetCluster(emc1);
 311          AliPicoTrack::GetEtaPhiDiff(mt, clm, dphi, deta);
 312          Printf("deta: %f dphi: %f",deta,dphi);
 313          }
 314       */
 315     }
 316     cluster = fCaloClustersCont->GetNextAcceptCluster();
 317   }
 318 }
 319
 320 //________________________________________________________________________
 321 void AliAnalysisTaskEmcalJetSample::ExecOnce() {
 322
 323   AliAnalysisTaskEmcalJet::ExecOnce();
 324
 325   if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
 326   if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
 327   if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
 328
 329 }
 330
 331 //________________________________________________________________________
 332 Bool_t AliAnalysisTaskEmcalJetSample::Run()
 333 {
 334   // Run analysis code here, if needed. It will be executed before FillHistograms().
 335
 336   return kTRUE;  // If return kFALSE FillHistogram() will NOT be executed.
 337 }
 338
 339 //________________________________________________________________________
 340 void AliAnalysisTaskEmcalJetSample::Terminate(Option_t *)
 341 {
 342   // Called once at the end of the analysis.
 343 }