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 "AliGenPythiaEventHeader.h"
  23 #include "AliMCEvent.h"
  24
  25 ClassImp(AliJetResponseMaker)
  26
  27 //________________________________________________________________________
  28 AliJetResponseMaker::AliJetResponseMaker() :
  29   AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE),
  30   fMCTracksName("MCParticles"),
  31   fMCJetsName("MCJets"),
  32   fMaxDistance(0.25),
  33   fDoWeighting(kFALSE),
  34   fVertexCut(10),
  35   fPythiaHeader(0),
  36   fEventWeight(0),
  37   fPtHardBin(0),
  38   fNTrials(0),
  39   fMCTracks(0),
  40   fMCJets(0),
  41   fHistNTrials(0),
  42   fHistAcceptedEvents(0),
  43   fHistEvents(0),
  44   fHistMCJetPhiEta(0),
  45   fHistMCJetsPt(0),
  46   fHistMCJetsNEFvsPt(0),
  47   fHistMCJetsZvsPt(0),
  48   fHistJetPhiEta(0),
  49   fHistJetsPt(0),
  50   fHistJetsNEFvsPt(0),
  51   fHistJetsZvsPt(0),
  52   fHistClosestDistance(0),
  53   fHistClosestDeltaPhi(0),
  54   fHistClosestDeltaEta(0),
  55   fHistClosestDeltaPt(0),
  56   fHistNonMatchedMCJetPt(0),
  57   fHistNonMatchedJetPt(0),
  58   fHistPartvsDetecPt(0),
  59   fHistMissedMCJets(0)
  60 {
  61   // Default constructor.
  62 }
  63
  64 //________________________________________________________________________
  65 AliJetResponseMaker::AliJetResponseMaker(const char *name) :
  66   AliAnalysisTaskEmcalJet(name, kTRUE),
  67   fMCTracksName("MCParticles"),
  68   fMCJetsName("MCJets"),
  69   fMaxDistance(0.25),
  70   fDoWeighting(kFALSE),
  71   fVertexCut(10),
  72   fPythiaHeader(0),
  73   fEventWeight(0),
  74   fPtHardBin(0),
  75   fNTrials(0),
  76   fMCTracks(0),
  77   fMCJets(0),
  78   fHistNTrials(0),
  79   fHistAcceptedEvents(0),
  80   fHistEvents(0),
  81   fHistMCJetPhiEta(0),
  82   fHistMCJetsPt(0),
  83   fHistMCJetsNEFvsPt(0),
  84   fHistMCJetsZvsPt(0),
  85   fHistJetPhiEta(0),
  86   fHistJetsPt(0),
  87   fHistJetsNEFvsPt(0),
  88   fHistJetsZvsPt(0),
  89   fHistClosestDistance(0),
  90   fHistClosestDeltaPhi(0),
  91   fHistClosestDeltaEta(0),
  92   fHistClosestDeltaPt(0),
  93   fHistNonMatchedMCJetPt(0),
  94   fHistNonMatchedJetPt(0),
  95   fHistPartvsDetecPt(0),
  96   fHistMissedMCJets(0)
  97 {
  98   // Standard constructor.
  99 }
 100
 101 //________________________________________________________________________
 102 AliJetResponseMaker::~AliJetResponseMaker()
 103 {
 104   // Destructor
 105 }
 106
 107 //________________________________________________________________________
 108 void AliJetResponseMaker::UserCreateOutputObjects()
 109 {
 110   // Create user objects.
 111
 112   OpenFile(1);
 113   fOutput = new TList();
 114   fOutput->SetOwner();
 115
 116   const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
 117   const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
 118
 119   fHistNTrials = new TH1F("fHistNTrials", "fHistNTrials", 11, 0, 11);
 120   fHistNTrials->GetXaxis()->SetTitle("p_{T} hard bin");
 121   fHistNTrials->GetYaxis()->SetTitle("trials");
 122   fOutput->Add(fHistNTrials);
 123
 124   fHistAcceptedEvents = new TH1F("fHistAcceptedEvents", "fHistAcceptedEvents", 11, 0, 11);
 125   fHistAcceptedEvents->GetXaxis()->SetTitle("p_{T} hard bin");
 126   fHistAcceptedEvents->GetYaxis()->SetTitle("accepted events");
 127   fOutput->Add(fHistAcceptedEvents);
 128
 129   fHistEvents = new TH1F("fHistEvents", "fHistEvents", 11, 0, 11);
 130   fHistEvents->GetXaxis()->SetTitle("p_{T} hard bin");
 131   fHistEvents->GetYaxis()->SetTitle("total events");
 132   fOutput->Add(fHistEvents);
 133
 134   for (Int_t i = 1; i < 12; i++) {
 135     fHistNTrials->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
 136     fHistAcceptedEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
 137     fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1]));
 138   }
 139
 140   fHistJetPhiEta = new TH2F("fHistJetPhiEta", "fHistJetPhiEta", 20, -2, 2, 32, 0, 6.4);
 141   fHistJetPhiEta->GetXaxis()->SetTitle("#eta");
 142   fHistJetPhiEta->GetYaxis()->SetTitle("#phi");
 143   fOutput->Add(fHistJetPhiEta);
 144
 145   fHistJetsPt = new TH1F("fHistJetsPt", "fHistJetsPt", fNbins, fMinBinPt, fMaxBinPt);
 146   fHistJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 147   fHistJetsPt->GetYaxis()->SetTitle("counts");
 148   fOutput->Add(fHistJetsPt);
 149
 150   fHistJetsZvsPt = new TH2F("fHistJetsZvsPt", "fHistJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 151   fHistJetsZvsPt->GetXaxis()->SetTitle("Z");
 152   fHistJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 153   fOutput->Add(fHistJetsZvsPt);
 154
 155   if (fAnaType == kEMCAL) {
 156     fHistJetsNEFvsPt = new TH2F("fHistJetsNEFvsPt", "fHistJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 157     fHistJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
 158     fHistJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 159     fOutput->Add(fHistJetsNEFvsPt);
 160   }
 161
 162   fHistMCJetPhiEta = new TH2F("fHistMCJetPhiEta", "fHistMCJetPhiEta", 20, -2, 2, 32, 0, 6.4);
 163   fHistMCJetPhiEta->GetXaxis()->SetTitle("#eta");
 164   fHistMCJetPhiEta->GetYaxis()->SetTitle("#phi");
 165   fOutput->Add(fHistMCJetPhiEta);
 166
 167   fHistMCJetsPt = new TH1F("fHistMCJetsPt", "fHistMCJetsPt", fNbins, fMinBinPt, fMaxBinPt);
 168   fHistMCJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 169   fHistMCJetsPt->GetYaxis()->SetTitle("counts");
 170   fOutput->Add(fHistMCJetsPt);
 171
 172   fHistMCJetsZvsPt = new TH2F("fHistMCJetsZvsPt", "fHistMCJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 173   fHistMCJetsZvsPt->GetXaxis()->SetTitle("Z");
 174   fHistMCJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 175   fOutput->Add(fHistMCJetsZvsPt);
 176
 177   if (fAnaType == kEMCAL) {
 178     fHistMCJetsNEFvsPt = new TH2F("fHistMCJetsNEFvsPt", "fHistMCJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
 179     fHistMCJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
 180     fHistMCJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 181     fOutput->Add(fHistMCJetsNEFvsPt);
 182   }
 183
 184   fHistClosestDistance = new TH1F("fHistClosestDistance", "fHistClosestDistance", 50, 0, fMaxDistance * 1.2);
 185   fHistClosestDistance->GetXaxis()->SetTitle("d");
 186   fHistClosestDistance->GetYaxis()->SetTitle("counts");
 187   fOutput->Add(fHistClosestDistance);
 188
 189   fHistClosestDeltaPhi = new TH1F("fHistClosestDeltaPhi", "fHistClosestDeltaPhi", 64, -1.6, 4.8);
 190   fHistClosestDeltaPhi->GetXaxis()->SetTitle("#Delta#phi");
 191   fHistClosestDeltaPhi->GetYaxis()->SetTitle("counts");
 192   fOutput->Add(fHistClosestDeltaPhi);
 193
 194   fHistClosestDeltaEta = new TH1F("fHistClosestDeltaEta", "fHistClosestDeltaEta", TMath::CeilNint(fMaxEta - fMinEta) * 20, fMinEta * 2, fMaxEta * 2);
 195   fHistClosestDeltaEta->GetXaxis()->SetTitle("#Delta#eta");
 196   fHistClosestDeltaEta->GetYaxis()->SetTitle("counts");
 197   fOutput->Add(fHistClosestDeltaEta);
 198
 199   fHistClosestDeltaPt = new TH1F("fHistClosestDeltaPt", "fHistClosestDeltaPt", fNbins, -fMaxBinPt / 2, fMaxBinPt / 2);
 200   fHistClosestDeltaPt->GetXaxis()->SetTitle("#Delta p_{T}");
 201   fHistClosestDeltaPt->GetYaxis()->SetTitle("counts");
 202   fOutput->Add(fHistClosestDeltaPt);
 203
 204   fHistNonMatchedMCJetPt = new TH1F("fHistNonMatchedMCJetPt", "fHistNonMatchedMCJetPt", fNbins, fMinBinPt, fMaxBinPt);
 205   fHistNonMatchedMCJetPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 206   fHistNonMatchedMCJetPt->GetYaxis()->SetTitle("counts");
 207   fOutput->Add(fHistNonMatchedMCJetPt);
 208
 209   fHistNonMatchedJetPt = new TH1F("fHistNonMatchedJetPt", "fHistNonMatchedJetPt", fNbins, fMinBinPt, fMaxBinPt);
 210   fHistNonMatchedJetPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 211   fHistNonMatchedJetPt->GetYaxis()->SetTitle("counts");
 212   fOutput->Add(fHistNonMatchedJetPt);
 213
 214   fHistPartvsDetecPt = new TH2F("fHistPartvsDetecPt", "fHistPartvsDetecPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
 215   fHistPartvsDetecPt->GetXaxis()->SetTitle("p_{T}^{rec}");
 216   fHistPartvsDetecPt->GetYaxis()->SetTitle("p_{T}^{gen}");
 217   fOutput->Add(fHistPartvsDetecPt);
 218
 219   fHistMissedMCJets = new TH1F("fHistMissedMCJets", "fHistMissedMCJets", fNbins, fMinBinPt, fMaxBinPt);
 220   fHistMissedMCJets->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 221   fHistMissedMCJets->GetYaxis()->SetTitle("counts");
 222   fOutput->Add(fHistMissedMCJets);
 223
 224   PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
 225 }
 226
 227 //________________________________________________________________________
 228 void AliJetResponseMaker::ExecOnce()
 229 {
 230   // Execute once.
 231
 232   if (!fMCJetsName.IsNull() && !fMCJets) {
 233     fMCJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCJetsName));
 234     if (!fMCJets) {
 235       AliError(Form("%s: Could not retrieve mc jets %s!", GetName(), fMCJetsName.Data()));
 236       return;
 237     }
 238     else if (!fMCJets->GetClass()->GetBaseClass("AliEmcalJet")) {
 239       AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fMCJetsName.Data()));
 240       fMCJets = 0;
 241       return;
 242     }
 243   }
 244
 245   if (!fMCTracksName.IsNull() && !fMCTracks) {
 246     fMCTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCTracksName));
 247     if (!fMCTracks) {
 248       AliError(Form("%s: Could not retrieve mc tracks %s!", GetName(), fMCTracksName.Data()));
 249       return;
 250     }
 251     else {
 252       TClass *cl = fMCTracks->GetClass();
 253       if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) {
 254         AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fMCTracksName.Data()));
 255         fMCTracks = 0;
 256         return;
 257       }
 258     }
 259   }
 260
 261   AliAnalysisTaskEmcalJet::ExecOnce();
 262 }
 263
 264 //________________________________________________________________________
 265 Bool_t AliJetResponseMaker::RetrieveEventObjects()
 266 {
 267   // Retrieve event objects.
 268
 269   if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
 270     return kFALSE;
 271
 272   fPythiaHeader = dynamic_cast<AliGenPythiaEventHeader*>(MCEvent()->GenEventHeader());
 273
 274   if (!fPythiaHeader)
 275     return kFALSE;
 276
 277   if (fDoWeighting)
 278     fEventWeight = fPythiaHeader->EventWeight();
 279   else
 280     fEventWeight = 1;
 281
 282   const Int_t ptHardLo[11] = { 0, 5,11,21,36,57, 84,117,152,191,234};
 283   const Int_t ptHardHi[11] = { 5,11,21,36,57,84,117,152,191,234,1000000};
 284
 285   Double_t pthard = fPythiaHeader->GetPtHard();
 286
 287   for (fPtHardBin = 0; fPtHardBin < 11; fPtHardBin++) {
 288     if (pthard >= ptHardLo[fPtHardBin] && pthard < ptHardHi[fPtHardBin])
 289       break;
 290   }
 291   fPtHardBin--;
 292
 293   fNTrials = fPythiaHeader->Trials();
 294
 295   return kTRUE;
 296 }
 297
 298 //________________________________________________________________________
 299 Bool_t AliJetResponseMaker::Run()
 300 {
 301   // Find the closest jets
 302
 303   fHistEvents->SetBinContent(fPtHardBin, fHistEvents->GetBinContent(fPtHardBin) + 1);
 304
 305   if (TMath::Abs(fVertex[2]) > fVertexCut)
 306     return kFALSE;
 307
 308   fHistAcceptedEvents->SetBinContent(fPtHardBin, fHistEvents->GetBinContent(fPtHardBin) + 1);
 309   fHistNTrials->SetBinContent(fPtHardBin, fHistNTrials->GetBinContent(fPtHardBin) + fNTrials);
 310
 311   DoJetLoop(fJets, fMCJets, kFALSE);
 312   DoJetLoop(fMCJets, fJets, kTRUE);
 313
 314   return kTRUE;
 315 }
 316
 317 //________________________________________________________________________
 318 void AliJetResponseMaker::DoJetLoop(TClonesArray *jets1, TClonesArray *jets2, Bool_t mc)
 319 {
 320   // Do the jet loop.
 321
 322   Int_t nJets1 = jets1->GetEntriesFast();
 323   Int_t nJets2 = jets2->GetEntriesFast();
 324
 325   for (Int_t i = 0; i < nJets1; i++) {
 326
 327     AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(jets1->At(i));
 328
 329     if (!jet1) {
 330       AliError(Form("Could not receive jet %d", i));
 331       continue;
 332     }
 333
 334     if (!AcceptJet(jet1, kTRUE, !mc))
 335       continue;
 336
 337     for (Int_t j = 0; j < nJets2; j++) {
 338
 339       AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(jets2->At(j));
 340
 341       if (!jet2) {
 342         AliError(Form("Could not receive jet %d", j));
 343         continue;
 344       }
 345
 346       if (!AcceptJet(jet2, kTRUE, mc))
 347         continue;
 348
 349       Double_t deta = jet2->Eta() - jet1->Eta();
 350       Double_t dphi = jet2->Phi() - jet1->Phi();
 351       Double_t d = TMath::Sqrt(deta * deta + dphi * dphi);
 352
 353       if (d < jet1->ClosestJetDistance()) {
 354         jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
 355         jet1->SetClosestJet(jet2, d);
 356       }
 357       else if (d < jet1->SecondClosestJetDistance()) {
 358         jet1->SetSecondClosestJet(jet2, d);
 359       }
 360     }
 361   }
 362 }
 363
 364
 365 //________________________________________________________________________
 366 Bool_t AliJetResponseMaker::FillHistograms()
 367 {
 368   // Fill histograms.
 369
 370   const Int_t nMCJets = fMCJets->GetEntriesFast();
 371
 372   for (Int_t i = 0; i < nMCJets; i++) {
 373
 374     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fMCJets->At(i));
 375
 376     if (!jet) {
 377       AliError(Form("Could not receive jet %d", i));
 378       continue;
 379     }
 380
 381     if (!AcceptJet(jet, kTRUE, kFALSE))
 382       continue;
 383
 384     if (jet->Pt() > fMaxBinPt)
 385       continue;
 386
 387     if (jet->ClosestJet() && jet->ClosestJet()->ClosestJet() == jet &&
 388         jet->ClosestJetDistance() < fMaxDistance) {    // Matched jet found
 389       jet->SetMatchedToClosest();
 390       jet->ClosestJet()->SetMatchedToClosest();
 391       if (jet->MatchedJet()->Pt() > fMaxBinPt) {
 392         fHistMissedMCJets->Fill(jet->Pt(), fEventWeight);
 393       }
 394       else {
 395         fHistClosestDistance->Fill(jet->ClosestJetDistance(), fEventWeight);
 396
 397         Double_t deta = jet->MatchedJet()->Eta() - jet->Eta();
 398         fHistClosestDeltaEta->Fill(deta, fEventWeight);
 399
 400         Double_t dphi = jet->MatchedJet()->Phi() - jet->Phi();
 401         fHistClosestDeltaPhi->Fill(dphi, fEventWeight);
 402
 403         Double_t dpt = jet->MatchedJet()->Pt() - jet->Pt();
 404         fHistClosestDeltaPt->Fill(dpt, fEventWeight);
 405
 406         fHistPartvsDetecPt->Fill(jet->MatchedJet()->Pt(), jet->Pt(), fEventWeight);
 407       }
 408     }
 409     else {
 410       fHistNonMatchedMCJetPt->Fill(jet->Pt(), fEventWeight);
 411       fHistMissedMCJets->Fill(jet->Pt(), fEventWeight);
 412     }
 413
 414     fHistMCJetsPt->Fill(jet->Pt(), fEventWeight);
 415
 416     fHistMCJetPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight);
 417
 418     if (fAnaType == kEMCAL)
 419       fHistMCJetsNEFvsPt->Fill(jet->NEF(), jet->Pt(), fEventWeight);
 420
 421     for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
 422       AliVParticle *track = jet->TrackAt(it, fMCTracks);
 423       if (track)
 424         fHistMCJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt(), fEventWeight);
 425     }
 426   }
 427
 428   const Int_t nJets = fJets->GetEntriesFast();
 429
 430   for (Int_t i = 0; i < nJets; i++) {
 431
 432     AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(i));
 433
 434     if (!jet) {
 435       AliError(Form("Could not receive mc jet %d", i));
 436       continue;
 437     }
 438
 439     if (!AcceptJet(jet))
 440       continue;
 441
 442     if (!jet->MatchedJet()) {
 443       fHistNonMatchedJetPt->Fill(jet->Pt(), fEventWeight);
 444     }
 445
 446     fHistJetsPt->Fill(jet->Pt(), fEventWeight);
 447
 448     fHistJetPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight);
 449
 450     if (fAnaType == kEMCAL)
 451       fHistJetsNEFvsPt->Fill(jet->NEF(), jet->Pt(), fEventWeight);
 452
 453     for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
 454       AliVParticle *track = jet->TrackAt(it, fTracks);
 455       if (track)
 456         fHistJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt(), fEventWeight);
 457     }
 458
 459     for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
 460       AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
 461       if (cluster) {
 462         TLorentzVector nP;
 463         cluster->GetMomentum(nP, fVertex);
 464         fHistJetsZvsPt->Fill(nP.Pt() / jet->Pt(), jet->Pt(), fEventWeight);
 465       }
 466     }
 467   }
 468
 469   return kTRUE;
 470 }