PWGGA/EMCALJetTasks/AliJetResponseMaker.cxx

   1 // $Id$
   2 //
   3 // Emcal jet response matrix maker task.
   4 //
   5 // Author: S. Aiola
   6
   7 #include "AliJetResponseMaker.h"
   8
   9 #include <TChain.h>
  10 #include <TClonesArray.h>
  11 #include <TH1F.h>
  12 #include <TH2F.h>
  13 #include <TList.h>
  14 #include <TLorentzVector.h>
  15
  16 #include "AliAnalysisManager.h"
  17 #include "AliCentrality.h"
  18 #include "AliFJWrapper.h"
  19 #include "AliVCluster.h"
  20 #include "AliVTrack.h"
  21 #include "AliEmcalJet.h"
  22 #include "AliMCEvent.h"
  23
  24 ClassImp(AliJetResponseMaker)
  25
  26 //________________________________________________________________________
  27 AliJetResponseMaker::AliJetResponseMaker() :
  28   AliAnalysisTaskEmcal("AliJetResponseMaker"),
  29   fMCTracksName("MCParticles"),
  30   fMCJetsName("MCJets"),
  31   fMaxDistance(0.2),
  32   fMCTracks(0),
  33   fMCJets(0),
  34   fHistMCJetPhiEta(0),
  35   fHistMCJetsPt(0),
  36   fHistMCJetsPtNonBias(0),
  37   fHistMCJetsNEFvsPt(0),
  38   fHistMCJetsZvsPt(0),
  39   fHistJetPhiEta(0),
  40   fHistJetsPt(0),
  41   fHistJetsPtNonBias(0),
  42   fHistJetsNEFvsPt(0),
  43   fHistJetsZvsPt(0),
  44   fHistClosestDistance(0),
  45   fHistClosestDeltaPhi(0),
  46   fHistClosestDeltaEta(0),
  47   fHistClosestDeltaPt(0),
  48   fHistPartvsDetecPt(0)
  49 {
  50   // Default constructor.
  51 }
  52
  53 //________________________________________________________________________
  54 AliJetResponseMaker::AliJetResponseMaker(const char *name) :
  55   AliAnalysisTaskEmcal(name),
  56   fMCTracksName("MCParticles"),
  57   fMCJetsName("MCJets"),
  58   fMaxDistance(0.2),
  59   fMCTracks(0),
  60   fMCJets(0),
  61   fHistMCJetPhiEta(0),
  62   fHistMCJetsPt(0),
  63   fHistMCJetsPtNonBias(0),
  64   fHistMCJetsNEFvsPt(0),
  65   fHistMCJetsZvsPt(0),
  66   fHistJetPhiEta(0),
  67   fHistJetsPt(0),
  68   fHistJetsPtNonBias(0),
  69   fHistJetsNEFvsPt(0),
  70   fHistJetsZvsPt(0),
  71   fHistClosestDistance(0),
  72   fHistClosestDeltaPhi(0),
  73   fHistClosestDeltaEta(0),
  74   fHistClosestDeltaPt(0),
  75   fHistPartvsDetecPt(0)
  76 {
  77   // Standard constructor.
  78
  79 }
  80
  81 //________________________________________________________________________
  82 AliJetResponseMaker::~AliJetResponseMaker()
  83 {
  84   // Destructor
  85 }
  86
  87 //________________________________________________________________________
  88 void AliJetResponseMaker::UserCreateOutputObjects()
  89 {
  90   // Create user objects.
  91
  92   OpenFile(1);
  93   fOutput = new TList();
  94   fOutput->SetOwner();
  95
  96   fHistJetPhiEta = new TH2F("fHistJetPhiEta", "fHistJetPhiEta", 20, -2, 2, 32, 0, 6.4);
  97   fHistJetPhiEta->GetXaxis()->SetTitle("#eta");
  98   fHistJetPhiEta->GetYaxis()->SetTitle("#phi");
  99   fOutput->Add(fHistJetPhiEta);
 100
 101   fHistJetsPt = new TH1F("fHistJetsPt", "fHistJetsPt", fNbins, fMinPt, fMaxPt);
 102   fHistJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 103   fHistJetsPt->GetYaxis()->SetTitle("counts");
 104   fOutput->Add(fHistJetsPt);
 105
 106   fHistJetsPtNonBias = new TH1F("fHistJetsPtNonBias", "fHistJetsPtNonBias", fNbins, fMinPt, fMaxPt);
 107   fHistJetsPtNonBias->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 108   fHistJetsPtNonBias->GetYaxis()->SetTitle("counts");
 109   fOutput->Add(fHistJetsPtNonBias);
 110
 111   fHistJetsZvsPt = new TH2F("fHistJetsZvsPt", "fHistJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
 112   fHistJetsZvsPt->GetXaxis()->SetTitle("Z");
 113   fHistJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 114   fOutput->Add(fHistJetsZvsPt);
 115
 116   if (fAnaType == kEMCAL) {
 117     fHistJetsNEFvsPt = new TH2F("fHistJetsNEFvsPt", "fHistJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
 118     fHistJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
 119     fHistJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 120     fOutput->Add(fHistJetsNEFvsPt);
 121   }
 122
 123   fHistMCJetPhiEta = new TH2F("fHistMCJetPhiEta", "fHistMCJetPhiEta", 20, -2, 2, 32, 0, 6.4);
 124   fHistMCJetPhiEta->GetXaxis()->SetTitle("#eta");
 125   fHistMCJetPhiEta->GetYaxis()->SetTitle("#phi");
 126   fOutput->Add(fHistMCJetPhiEta);
 127
 128   fHistMCJetsPt = new TH1F("fHistMCJetsPt", "fHistMCJetsPt", fNbins, fMinPt, fMaxPt);
 129   fHistMCJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 130   fHistMCJetsPt->GetYaxis()->SetTitle("counts");
 131   fOutput->Add(fHistMCJetsPt);
 132
 133   fHistMCJetsPtNonBias = new TH1F("fHistMCJetsPtNonBias", "fHistMCJetsPtNonBias", fNbins, fMinPt, fMaxPt);
 134   fHistMCJetsPtNonBias->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 135   fHistMCJetsPtNonBias->GetYaxis()->SetTitle("counts");
 136   fOutput->Add(fHistMCJetsPtNonBias);
 137
 138   fHistMCJetsZvsPt = new TH2F("fHistMCJetsZvsPt", "fHistMCJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
 139   fHistMCJetsZvsPt->GetXaxis()->SetTitle("Z");
 140   fHistMCJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 141   fOutput->Add(fHistMCJetsZvsPt);
 142
 143   if (fAnaType == kEMCAL) {
 144     fHistMCJetsNEFvsPt = new TH2F("fHistMCJetsNEFvsPt", "fHistMCJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
 145     fHistMCJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
 146     fHistMCJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 147     fOutput->Add(fHistMCJetsNEFvsPt);
 148   }
 149
 150   fHistClosestDistance = new TH1F("fHistClosestDistance", "fHistClosestDistance", 50, 0, fMaxDistance * 1.2);
 151   fHistClosestDistance->GetXaxis()->SetTitle("d");
 152   fHistClosestDistance->GetYaxis()->SetTitle("counts");
 153   fOutput->Add(fHistClosestDistance);
 154
 155   fHistClosestDeltaPhi = new TH1F("fHistClosestDeltaPhi", "fHistClosestDeltaPhi", 64, -1.6, 4.8);
 156   fHistClosestDeltaPhi->GetXaxis()->SetTitle("#Delta#phi");
 157   fHistClosestDeltaPhi->GetYaxis()->SetTitle("counts");
 158   fOutput->Add(fHistClosestDeltaPhi);
 159
 160   fHistClosestDeltaEta = new TH1F("fHistClosestDeltaEta", "fHistClosestDeltaEta", TMath::CeilNint(fMaxEta - fMinEta) * 20, fMinEta * 2, fMaxEta * 2);
 161   fHistClosestDeltaEta->GetXaxis()->SetTitle("#Delta#eta");
 162   fHistClosestDeltaEta->GetYaxis()->SetTitle("counts");
 163   fOutput->Add(fHistClosestDeltaEta);
 164
 165   fHistClosestDeltaPt = new TH1F("fHistClosestDeltaPt", "fHistClosestDeltaPt", fNbins, -fMaxPt / 2, fMaxPt / 2);
 166   fHistClosestDeltaPt->GetXaxis()->SetTitle("#Delta p_{T}");
 167   fHistClosestDeltaPt->GetYaxis()->SetTitle("counts");
 168   fOutput->Add(fHistClosestDeltaPt);
 169
 170   fHistPartvsDetecPt = new TH2F("fHistPartvsDetecPt", "fHistPartvsDetecPt", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
 171   fHistPartvsDetecPt->GetXaxis()->SetTitle("p_{T}^{det}");
 172   fHistPartvsDetecPt->GetYaxis()->SetTitle("p_{T}^{rec}");
 173   fOutput->Add(fHistPartvsDetecPt);
 174
 175   PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
 176 }
 177
 178 //________________________________________________________________________
 179 void AliJetResponseMaker::FillHistograms()
 180 {
 181   // Fill histograms.
 182
 183   // Find the closest jets
 184   DoJetLoop(fJets, fMCJets);
 185   DoJetLoop(fMCJets, fJets);
 186
 187   Int_t nJets = fJets->GetEntriesFast();
 188
 189   for (Int_t i = 0; i < nJets; i++) {
 190
 191     AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fJets->At(i));
 192
 193     if (!jet) {
 194       AliError(Form("Could not receive jet %d", i));
 195       continue;
 196     }
 197
 198     if (!AcceptJet(jet))
 199       continue;
 200
 201     fHistJetsPtNonBias->Fill(jet->Pt());
 202
 203     if (jet->MaxTrackPt() < fPtBiasJetTrack && (fAnaType == kTPC || jet->MaxClusterPt() < fPtBiasJetClus))
 204         continue;
 205
 206     if (jet->ClosestJet() && jet->ClosestJet()->ClosestJet() == jet && jet->ClosestJetDistance() < fMaxDistance) {    // Matched jet found
 207       jet->SetMatchedToClosest();
 208       jet->ClosestJet()->SetMatchedToClosest();
 209       fHistClosestDistance->Fill(jet->ClosestJetDistance());
 210       Double_t deta = jet->Eta() - jet->MatchedJet()->Eta();
 211       fHistClosestDeltaEta->Fill(deta);
 212       Double_t dphi = jet->Phi() - jet->MatchedJet()->Phi();
 213       fHistClosestDeltaPhi->Fill(dphi);
 214       Double_t dpt = jet->Pt() - jet->MatchedJet()->Pt();
 215       fHistClosestDeltaPt->Fill(dpt);
 216       fHistPartvsDetecPt->Fill(jet->Pt(), jet->MatchedJet()->Pt());
 217     }
 218
 219     fHistJetsPt->Fill(jet->Pt());
 220
 221     fHistJetPhiEta->Fill(jet->Eta(), jet->Phi());
 222
 223     if (fAnaType == kEMCAL)
 224       fHistJetsNEFvsPt->Fill(jet->NEF(), jet->Pt());
 225
 226     for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
 227       AliVParticle *track = jet->TrackAt(it, fTracks);
 228       if (track)
 229         fHistJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt());
 230     }
 231
 232     for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
 233       AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
 234       if (cluster) {
 235         TLorentzVector nP;
 236         cluster->GetMomentum(nP, fVertex);
 237         fHistJetsZvsPt->Fill(nP.Pt() / jet->Pt(), jet->Pt());
 238       }
 239     }
 240   }
 241
 242  Int_t nMCJets = fMCJets->GetEntriesFast();
 243
 244   for (Int_t i = 0; i < nMCJets; i++) {
 245
 246     AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fMCJets->At(i));
 247
 248     if (!jet) {
 249       AliError(Form("Could not receive mc jet %d", i));
 250       continue;
 251     }
 252
 253     if (!AcceptJet(jet))
 254       continue;
 255
 256     fHistMCJetsPtNonBias->Fill(jet->Pt());
 257
 258     if (jet->MaxTrackPt() < fPtBiasJetTrack)
 259         continue;
 260
 261     fHistMCJetsPt->Fill(jet->Pt());
 262
 263     fHistMCJetPhiEta->Fill(jet->Eta(), jet->Phi());
 264
 265     if (fAnaType == kEMCAL)
 266       fHistMCJetsNEFvsPt->Fill(jet->NEF(), jet->Pt());
 267
 268     for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
 269       AliVParticle *track = jet->TrackAt(it, fMCTracks);
 270       if (track)
 271         fHistMCJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt());
 272     }
 273   }
 274 }
 275
 276 //________________________________________________________________________
 277 void AliJetResponseMaker::DoJetLoop(TClonesArray *jets1, TClonesArray *jets2)
 278 {
 279   // Do the jet loop.
 280
 281   Int_t nJets1 = jets1->GetEntriesFast();
 282   Int_t nJets2 = jets2->GetEntriesFast();
 283
 284   for (Int_t i = 0; i < nJets1; i++) {
 285
 286     AliEmcalJet* jet1 = dynamic_cast<AliEmcalJet*>(jets1->At(i));
 287
 288     if (!jet1) {
 289       AliError(Form("Could not receive jet %d", i));
 290       continue;
 291     }
 292
 293     if (!AcceptJet(jet1))
 294       continue;
 295
 296     if (jet1->MaxTrackPt() < fPtBiasJetTrack && (fAnaType == kTPC || jet1->IsMC() || jet1->MaxClusterPt() < fPtBiasJetClus))
 297         continue;
 298
 299     for (Int_t j = 0; j < nJets2; j++) {
 300
 301       AliEmcalJet* jet2 = dynamic_cast<AliEmcalJet*>(jets2->At(j));
 302
 303       if (!jet2) {
 304         AliError(Form("Could not receive jet %d", j));
 305         continue;
 306       }
 307
 308       if (!AcceptJet(jet2))
 309         continue;
 310
 311       if (jet2->MaxTrackPt() < fPtBiasJetTrack && (fAnaType == kTPC || jet2->IsMC() || jet2->MaxClusterPt() < fPtBiasJetClus))
 312         continue;
 313
 314       Double_t deta = jet2->Eta() - jet1->Eta();
 315       Double_t dphi = jet2->Phi() - jet1->Phi();
 316       Double_t d = TMath::Sqrt(deta * deta + dphi * dphi);
 317
 318       if (d < jet1->ClosestJetDistance()) {
 319         jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
 320         jet1->SetClosestJet(jet2, d);
 321       }
 322       else if (d < jet1->SecondClosestJetDistance()) {
 323         jet1->SetSecondClosestJet(jet2, d);
 324       }
 325     }
 326   }
 327 }
 328
 329 //________________________________________________________________________
 330 void AliJetResponseMaker::RetrieveEventObjects()
 331 {
 332   // Retrieve event objects.
 333
 334   AliAnalysisTaskEmcal::RetrieveEventObjects();
 335
 336   if (!fMCJetsName.IsNull()) {
 337     fMCJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCJetsName));
 338     if (!fMCJets) {
 339       AliWarning(Form("Could not retrieve MC jets %s!", fMCJetsName.Data()));
 340     }
 341   }
 342
 343   if (!fMCTracksName.IsNull()) {
 344     fMCTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCTracksName));
 345     if (!fMCJets) {
 346       AliWarning(Form("Could not retrieve MC tracks %s!", fMCTracksName.Data()));
 347     }
 348   }
 349 }
 350
 351 //________________________________________________________________________
 352 void AliJetResponseMaker::Terminate(Option_t *)
 353 {
 354   // Called once at the end of the analysis.
 355 }