]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGGA/EMCALJetTasks/AliJetResponseMaker.cxx
name of task
[u/mrichter/AliRoot.git] / PWGGA / EMCALJetTasks / AliJetResponseMaker.cxx
index d1abe0648272967ed22c65c09a1f2a45287eef78..6607b8366fb7ee1984dd0d5c08e41e3c6aff9230 100644 (file)
@@ -25,57 +25,58 @@ ClassImp(AliJetResponseMaker)
 
 //________________________________________________________________________
 AliJetResponseMaker::AliJetResponseMaker() : 
-  AliAnalysisTaskEmcal("AliJetResponseMaker"),
+  AliAnalysisTaskEmcalJet("AliJetResponseMaker", kTRUE),
   fMCTracksName("MCParticles"),
   fMCJetsName("MCJets"),
-  fMaxDistance(0.2),
+  fMaxDistance(0.25),
   fMCTracks(0),
   fMCJets(0),
   fHistMCJetPhiEta(0),
   fHistMCJetsPt(0),
-  fHistMCJetsPtNonBias(0),
   fHistMCJetsNEFvsPt(0),
   fHistMCJetsZvsPt(0),
   fHistJetPhiEta(0),
   fHistJetsPt(0),
-  fHistJetsPtNonBias(0),
   fHistJetsNEFvsPt(0),
   fHistJetsZvsPt(0),
   fHistClosestDistance(0),
   fHistClosestDeltaPhi(0),
   fHistClosestDeltaEta(0),
   fHistClosestDeltaPt(0),
-  fHistPartvsDetecPt(0)
+  fHistNonMatchedMCJetPt(0),
+  fHistNonMatchedJetPt(0),
+  fHistPartvsDetecPt(0),
+  fHistMissedMCJets(0)
 {
   // Default constructor.
 }
 
 //________________________________________________________________________
 AliJetResponseMaker::AliJetResponseMaker(const char *name) : 
-  AliAnalysisTaskEmcal(name),
+  AliAnalysisTaskEmcalJet(name, kTRUE),
   fMCTracksName("MCParticles"),
   fMCJetsName("MCJets"),
-  fMaxDistance(0.2),
+  fMaxDistance(0.25),
   fMCTracks(0),
   fMCJets(0),
   fHistMCJetPhiEta(0),
   fHistMCJetsPt(0),
-  fHistMCJetsPtNonBias(0),
   fHistMCJetsNEFvsPt(0),
   fHistMCJetsZvsPt(0),
   fHistJetPhiEta(0),
   fHistJetsPt(0),
-  fHistJetsPtNonBias(0),
   fHistJetsNEFvsPt(0),
   fHistJetsZvsPt(0),
   fHistClosestDistance(0),
   fHistClosestDeltaPhi(0),
   fHistClosestDeltaEta(0),
   fHistClosestDeltaPt(0),
-  fHistPartvsDetecPt(0)
+  fHistNonMatchedMCJetPt(0),
+  fHistNonMatchedJetPt(0),
+  fHistPartvsDetecPt(0),
+  fHistMissedMCJets(0)
 {
   // Standard constructor.
-
 }
 
 //________________________________________________________________________
@@ -98,23 +99,18 @@ void AliJetResponseMaker::UserCreateOutputObjects()
   fHistJetPhiEta->GetYaxis()->SetTitle("#phi");
   fOutput->Add(fHistJetPhiEta);
   
-  fHistJetsPt = new TH1F("fHistJetsPt", "fHistJetsPt", fNbins, fMinPt, fMaxPt);
+  fHistJetsPt = new TH1F("fHistJetsPt", "fHistJetsPt", fNbins, fMinBinPt, fMaxBinPt);
   fHistJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
   fHistJetsPt->GetYaxis()->SetTitle("counts");
   fOutput->Add(fHistJetsPt);
   
-  fHistJetsPtNonBias = new TH1F("fHistJetsPtNonBias", "fHistJetsPtNonBias", fNbins, fMinPt, fMaxPt);
-  fHistJetsPtNonBias->GetXaxis()->SetTitle("p_{T} [GeV/c]");
-  fHistJetsPtNonBias->GetYaxis()->SetTitle("counts");
-  fOutput->Add(fHistJetsPtNonBias);
-  
-  fHistJetsZvsPt = new TH2F("fHistJetsZvsPt", "fHistJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
+  fHistJetsZvsPt = new TH2F("fHistJetsZvsPt", "fHistJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
   fHistJetsZvsPt->GetXaxis()->SetTitle("Z");
   fHistJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
   fOutput->Add(fHistJetsZvsPt);
   
   if (fAnaType == kEMCAL) {
-    fHistJetsNEFvsPt = new TH2F("fHistJetsNEFvsPt", "fHistJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
+    fHistJetsNEFvsPt = new TH2F("fHistJetsNEFvsPt", "fHistJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
     fHistJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
     fHistJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
     fOutput->Add(fHistJetsNEFvsPt);
@@ -125,23 +121,18 @@ void AliJetResponseMaker::UserCreateOutputObjects()
   fHistMCJetPhiEta->GetYaxis()->SetTitle("#phi");
   fOutput->Add(fHistMCJetPhiEta);
   
-  fHistMCJetsPt = new TH1F("fHistMCJetsPt", "fHistMCJetsPt", fNbins, fMinPt, fMaxPt);
+  fHistMCJetsPt = new TH1F("fHistMCJetsPt", "fHistMCJetsPt", fNbins, fMinBinPt, fMaxBinPt);
   fHistMCJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
   fHistMCJetsPt->GetYaxis()->SetTitle("counts");
   fOutput->Add(fHistMCJetsPt);
   
-  fHistMCJetsPtNonBias = new TH1F("fHistMCJetsPtNonBias", "fHistMCJetsPtNonBias", fNbins, fMinPt, fMaxPt);
-  fHistMCJetsPtNonBias->GetXaxis()->SetTitle("p_{T} [GeV/c]");
-  fHistMCJetsPtNonBias->GetYaxis()->SetTitle("counts");
-  fOutput->Add(fHistMCJetsPtNonBias);
-  
-  fHistMCJetsZvsPt = new TH2F("fHistMCJetsZvsPt", "fHistMCJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
+  fHistMCJetsZvsPt = new TH2F("fHistMCJetsZvsPt", "fHistMCJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
   fHistMCJetsZvsPt->GetXaxis()->SetTitle("Z");
   fHistMCJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
   fOutput->Add(fHistMCJetsZvsPt);
   
   if (fAnaType == kEMCAL) {
-    fHistMCJetsNEFvsPt = new TH2F("fHistMCJetsNEFvsPt", "fHistMCJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
+    fHistMCJetsNEFvsPt = new TH2F("fHistMCJetsNEFvsPt", "fHistMCJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt);
     fHistMCJetsNEFvsPt->GetXaxis()->SetTitle("NEF");
     fHistMCJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
     fOutput->Add(fHistMCJetsNEFvsPt);
@@ -162,119 +153,143 @@ void AliJetResponseMaker::UserCreateOutputObjects()
   fHistClosestDeltaEta->GetYaxis()->SetTitle("counts");
   fOutput->Add(fHistClosestDeltaEta);
 
-  fHistClosestDeltaPt = new TH1F("fHistClosestDeltaPt", "fHistClosestDeltaPt", fNbins, -fMaxPt / 2, fMaxPt / 2);
+  fHistClosestDeltaPt = new TH1F("fHistClosestDeltaPt", "fHistClosestDeltaPt", fNbins, -fMaxBinPt / 2, fMaxBinPt / 2);
   fHistClosestDeltaPt->GetXaxis()->SetTitle("#Delta p_{T}");
   fHistClosestDeltaPt->GetYaxis()->SetTitle("counts");
   fOutput->Add(fHistClosestDeltaPt);
 
-  fHistPartvsDetecPt = new TH2F("fHistPartvsDetecPt", "fHistPartvsDetecPt", fNbins, fMinPt, fMaxPt, fNbins, fMinPt, fMaxPt);
-  fHistPartvsDetecPt->GetXaxis()->SetTitle("p_{T}^{det}");
-  fHistPartvsDetecPt->GetYaxis()->SetTitle("p_{T}^{rec}");
+  fHistNonMatchedMCJetPt = new TH1F("fHistNonMatchedMCJetPt", "fHistNonMatchedMCJetPt", fNbins, fMinBinPt, fMaxBinPt);
+  fHistNonMatchedMCJetPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+  fHistNonMatchedMCJetPt->GetYaxis()->SetTitle("counts");
+  fOutput->Add(fHistNonMatchedMCJetPt);
+
+  fHistNonMatchedJetPt = new TH1F("fHistNonMatchedJetPt", "fHistNonMatchedJetPt", fNbins, fMinBinPt, fMaxBinPt);
+  fHistNonMatchedJetPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+  fHistNonMatchedJetPt->GetYaxis()->SetTitle("counts");
+  fOutput->Add(fHistNonMatchedJetPt);
+
+  fHistPartvsDetecPt = new TH2F("fHistPartvsDetecPt", "fHistPartvsDetecPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt);
+  fHistPartvsDetecPt->GetXaxis()->SetTitle("p_{T}^{rec}");
+  fHistPartvsDetecPt->GetYaxis()->SetTitle("p_{T}^{gen}");
   fOutput->Add(fHistPartvsDetecPt);
 
-  PostData(1, fOutput); // Post data for ALL output slots >0 here, to get at least an empty histogram
+  fHistMissedMCJets = new TH1F("fHistMissedMCJets", "fHistMissedMCJets", fNbins, fMinBinPt, fMaxBinPt);
+  fHistMissedMCJets->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+  fHistMissedMCJets->GetYaxis()->SetTitle("counts");
+  fOutput->Add(fHistMissedMCJets);
+
+  PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram
 }
 
 //________________________________________________________________________
-void AliJetResponseMaker::FillHistograms()
+Bool_t AliJetResponseMaker::FillHistograms()
 {
   // Fill histograms.
 
   // Find the closest jets
-  DoJetLoop(fJets, fMCJets);
-  DoJetLoop(fMCJets, fJets);
+  DoJetLoop(fJets, fMCJets, kFALSE);
+  DoJetLoop(fMCJets, fJets, kTRUE);
 
-  Int_t nJets = fJets->GetEntriesFast();
+  const Int_t nMCJets = fMCJets->GetEntriesFast();
 
-  for (Int_t i = 0; i < nJets; i++) {
+  for (Int_t i = 0; i < nMCJets; i++) {
 
-    AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fJets->At(i));
+    AliEmcalJet* jet = static_cast<AliEmcalJet*>(fMCJets->At(i));
 
     if (!jet) {
       AliError(Form("Could not receive jet %d", i));
       continue;
     }  
 
-    if (!AcceptJet(jet))
+    if (!AcceptJet(jet, kTRUE, kFALSE))
       continue;
 
-    fHistJetsPtNonBias->Fill(jet->Pt());
-
-    if (jet->MaxTrackPt() < fPtBiasJetTrack && (fAnaType == kTPC || jet->MaxClusterPt() < fPtBiasJetClus))
-       continue;
+    if (jet->Pt() > fMaxBinPt)
+      continue;
 
-    if (jet->ClosestJet() && jet->ClosestJet()->ClosestJet() == jet && jet->ClosestJetDistance() < fMaxDistance) {    // Matched jet found
+    if (jet->ClosestJet() && jet->ClosestJet()->ClosestJet() == jet && 
+        jet->ClosestJetDistance() < fMaxDistance) {    // Matched jet found
       jet->SetMatchedToClosest();
       jet->ClosestJet()->SetMatchedToClosest();
-      fHistClosestDistance->Fill(jet->ClosestJetDistance());
-      Double_t deta = jet->Eta() - jet->MatchedJet()->Eta();
-      fHistClosestDeltaEta->Fill(deta);
-      Double_t dphi = jet->Phi() - jet->MatchedJet()->Phi();
-      fHistClosestDeltaPhi->Fill(dphi);
-      Double_t dpt = jet->Pt() - jet->MatchedJet()->Pt();
-      fHistClosestDeltaPt->Fill(dpt);
-      fHistPartvsDetecPt->Fill(jet->Pt(), jet->MatchedJet()->Pt());
+      if (jet->MatchedJet()->Pt() > fMaxBinPt) {
+       fHistMissedMCJets->Fill(jet->Pt());
+      }
+      else {
+       fHistClosestDistance->Fill(jet->ClosestJetDistance());
+       Double_t deta = jet->MatchedJet()->Eta() - jet->Eta();
+       fHistClosestDeltaEta->Fill(deta);
+       Double_t dphi = jet->MatchedJet()->Phi() - jet->Phi();
+       fHistClosestDeltaPhi->Fill(dphi);
+       Double_t dpt = jet->MatchedJet()->Pt() - jet->Pt();
+       fHistClosestDeltaPt->Fill(dpt);
+       fHistPartvsDetecPt->Fill(jet->MatchedJet()->Pt(), jet->Pt());
+      }
+    }
+    else {
+      fHistNonMatchedMCJetPt->Fill(jet->Pt());
+      fHistMissedMCJets->Fill(jet->Pt());
     }
 
-    fHistJetsPt->Fill(jet->Pt());
+    fHistMCJetsPt->Fill(jet->Pt());
 
-    fHistJetPhiEta->Fill(jet->Eta(), jet->Phi());
+    fHistMCJetPhiEta->Fill(jet->Eta(), jet->Phi());
 
     if (fAnaType == kEMCAL)
-      fHistJetsNEFvsPt->Fill(jet->NEF(), jet->Pt());
+      fHistMCJetsNEFvsPt->Fill(jet->NEF(), jet->Pt());
 
     for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
-      AliVParticle *track = jet->TrackAt(it, fTracks);
+      AliVParticle *track = jet->TrackAt(it, fMCTracks);
       if (track)
-       fHistJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt());
-    }
-
-    for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
-      AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
-      if (cluster) {
-       TLorentzVector nP;
-       cluster->GetMomentum(nP, fVertex);
-       fHistJetsZvsPt->Fill(nP.Pt() / jet->Pt(), jet->Pt());
-      }
+       fHistMCJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt());
     }
   }
 
Int_t nMCJets = fMCJets->GetEntriesFast();
 const Int_t nJets = fJets->GetEntriesFast();
 
-  for (Int_t i = 0; i < nMCJets; i++) {
+  for (Int_t i = 0; i < nJets; i++) {
 
-    AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fMCJets->At(i));
+    AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(i));
 
     if (!jet) {
       AliError(Form("Could not receive mc jet %d", i));
       continue;
     }  
-
+    
     if (!AcceptJet(jet))
       continue;
 
-    fHistMCJetsPtNonBias->Fill(jet->Pt());
-
-    if (jet->MaxTrackPt() < fPtBiasJetTrack)
-       continue;
+    if (!jet->MatchedJet()) {
+      fHistNonMatchedJetPt->Fill(jet->Pt());
+    }
 
-    fHistMCJetsPt->Fill(jet->Pt());
+    fHistJetsPt->Fill(jet->Pt());
 
-    fHistMCJetPhiEta->Fill(jet->Eta(), jet->Phi());
+    fHistJetPhiEta->Fill(jet->Eta(), jet->Phi());
 
     if (fAnaType == kEMCAL)
-      fHistMCJetsNEFvsPt->Fill(jet->NEF(), jet->Pt());
+      fHistJetsNEFvsPt->Fill(jet->NEF(), jet->Pt());
 
     for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
-      AliVParticle *track = jet->TrackAt(it, fMCTracks);
+      AliVParticle *track = jet->TrackAt(it, fTracks);
       if (track)
-       fHistMCJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt());
+       fHistJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt());
+    }
+
+    for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
+      AliVCluster *cluster = jet->ClusterAt(ic, fCaloClusters);
+      if (cluster) {
+       TLorentzVector nP;
+       cluster->GetMomentum(nP, fVertex);
+       fHistJetsZvsPt->Fill(nP.Pt() / jet->Pt(), jet->Pt());
+      }
     }
   }
+
+  return kTRUE;
 }
 
 //________________________________________________________________________
-void AliJetResponseMaker::DoJetLoop(TClonesArray *jets1, TClonesArray *jets2)
+void AliJetResponseMaker::DoJetLoop(TClonesArray *jets1, TClonesArray *jets2, Bool_t mc)
 {
   // Do the jet loop.
 
@@ -283,32 +298,26 @@ void AliJetResponseMaker::DoJetLoop(TClonesArray *jets1, TClonesArray *jets2)
 
   for (Int_t i = 0; i < nJets1; i++) {
 
-    AliEmcalJet* jet1 = dynamic_cast<AliEmcalJet*>(jets1->At(i));
+    AliEmcalJet* jet1 = static_cast<AliEmcalJet*>(jets1->At(i));
 
     if (!jet1) {
       AliError(Form("Could not receive jet %d", i));
       continue;
     }  
 
-    if (!AcceptJet(jet1))
+    if (!AcceptJet(jet1, kTRUE, mc))
       continue;
-    
-    if (jet1->MaxTrackPt() < fPtBiasJetTrack && (fAnaType == kTPC || jet1->IsMC() || jet1->MaxClusterPt() < fPtBiasJetClus))
-       continue;
 
     for (Int_t j = 0; j < nJets2; j++) {
       
-      AliEmcalJet* jet2 = dynamic_cast<AliEmcalJet*>(jets2->At(j));
+      AliEmcalJet* jet2 = static_cast<AliEmcalJet*>(jets2->At(j));
       
       if (!jet2) {
        AliError(Form("Could not receive jet %d", j));
        continue;
       }  
       
-      if (!AcceptJet(jet2))
-       continue;
-      
-      if (jet2->MaxTrackPt() < fPtBiasJetTrack && (fAnaType == kTPC || jet2->IsMC() || jet2->MaxClusterPt() < fPtBiasJetClus))
+      if (!AcceptJet(jet2, kTRUE, !mc))
        continue;
       
       Double_t deta = jet2->Eta() - jet1->Eta();
@@ -327,25 +336,43 @@ void AliJetResponseMaker::DoJetLoop(TClonesArray *jets1, TClonesArray *jets2)
 }
 
 //________________________________________________________________________
-void AliJetResponseMaker::RetrieveEventObjects()
+Bool_t AliJetResponseMaker::RetrieveEventObjects()
 {
   // Retrieve event objects.
 
-  AliAnalysisTaskEmcal::RetrieveEventObjects();
+  if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
+    return kFALSE;
   
-  if (!fMCJetsName.IsNull()) {
+  if (!fMCJetsName.IsNull() && !fMCJets) {
     fMCJets = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCJetsName));
     if (!fMCJets) {
-      AliWarning(Form("Could not retrieve MC jets %s!", fMCJetsName.Data())); 
+      AliError(Form("%s: Could not retrieve mc jets %s!", GetName(), fMCJetsName.Data()));
+      return kFALSE;
+    }
+    else if (!fMCJets->GetClass()->GetBaseClass("AliEmcalJet")) {
+      AliError(Form("%s: Collection %s does not contain AliEmcalJet objects!", GetName(), fMCJetsName.Data())); 
+      fMCJets = 0;
+      return kFALSE;
     }
   }
 
-  if (!fMCTracksName.IsNull()) {
+  if (!fMCTracksName.IsNull() && !fMCTracks) {
     fMCTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fMCTracksName));
-    if (!fMCJets) {
-      AliWarning(Form("Could not retrieve MC tracks %s!", fMCTracksName.Data())); 
+    if (!fMCTracks) {
+      AliError(Form("%s: Could not retrieve mc tracks %s!", GetName(), fMCTracksName.Data())); 
+      return kFALSE;
+    }
+    else {
+      TClass *cl = fMCTracks->GetClass();
+      if (!cl->GetBaseClass("AliVParticle") && !cl->GetBaseClass("AliEmcalParticle")) {
+       AliError(Form("%s: Collection %s does not contain AliVParticle nor AliEmcalParticle objects!", GetName(), fMCTracksName.Data())); 
+       fMCTracks = 0;
+       return kFALSE;
+      }
     }
   }
+
+  return kTRUE;
 }
 
 //________________________________________________________________________