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