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