response maker
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 May 2012 19:35:29 +0000 (19:35 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 May 2012 19:35:29 +0000 (19:35 +0000)
PWGGA/EMCALJetTasks/AliJetResponseMaker.cxx [new file with mode: 0644]
PWGGA/EMCALJetTasks/AliJetResponseMaker.h [new file with mode: 0644]

diff --git a/PWGGA/EMCALJetTasks/AliJetResponseMaker.cxx b/PWGGA/EMCALJetTasks/AliJetResponseMaker.cxx
new file mode 100644 (file)
index 0000000..a8e8a4e
--- /dev/null
@@ -0,0 +1,249 @@
+// $Id: AliJetResponseMaker.cxx  $
+//
+// Emcal jet response matrix maker task.
+//
+// Author: S. Aiola
+
+#include "AliJetResponseMaker.h"
+
+#include <TChain.h>
+#include <TClonesArray.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TList.h>
+#include <TLorentzVector.h>
+
+#include "AliAnalysisManager.h"
+#include "AliCentrality.h"
+#include "AliFJWrapper.h"
+#include "AliVCluster.h"
+#include "AliVTrack.h"
+#include "AliEmcalJet.h"
+#include "AliMCEvent.h"
+
+ClassImp(AliJetResponseMaker)
+
+//________________________________________________________________________
+AliJetResponseMaker::AliJetResponseMaker() : 
+  AliAnalysisTaskEmcal("AliJetResponseMaker"),
+  fMCPartName("Tracks"),
+  fMCJetsName("mcJets"),
+  fMCParts(0),
+  fMCJets(0),
+  fHistMCJetPhiEta(0),
+  fHistMCJetsPt(0),
+  fHistMCJetsPtTrack(0),
+  fHistMCJetsPtClus(0),
+  fHistMCJetsPtNonBias(0),
+  fHistMCLeadingJetPt(0),
+  fHistMCJetsNEFvsPt(0),
+  fHistMCJetsZvsPt(0),
+  fHistJetPhiEta(0),
+  fHistJetsPt(0),
+  fHistJetsPtTrack(0),
+  fHistJetsPtClus(0),
+  fHistJetsPtNonBias(0),
+  fHistLeadingJetPt(0),
+  fHistJetsNEFvsPt(0),
+  fHistJetsZvsPt(0)
+{
+  // Default constructor.
+}
+
+//________________________________________________________________________
+AliJetResponseMaker::AliJetResponseMaker(const char *name) : 
+  AliAnalysisTaskEmcal(name),
+  fMCPartName("Tracks"),
+  fMCJetsName("mcJets"),
+  fMCParts(0),
+  fMCJets(0),
+  fHistMCJetPhiEta(0),
+  fHistMCJetsPt(0),
+  fHistMCJetsPtTrack(0),
+  fHistMCJetsPtClus(0),
+  fHistMCJetsPtNonBias(0),
+  fHistMCLeadingJetPt(0),
+  fHistMCJetsNEFvsPt(0),
+  fHistMCJetsZvsPt(0),
+  fHistJetPhiEta(0),
+  fHistJetsPt(0),
+  fHistJetsPtTrack(0),
+  fHistJetsPtClus(0),
+  fHistJetsPtNonBias(0),
+  fHistLeadingJetPt(0),
+  fHistJetsNEFvsPt(0),
+  fHistJetsZvsPt(0)
+{
+  // Standard constructor.
+
+  fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
+}
+
+//________________________________________________________________________
+AliJetResponseMaker::~AliJetResponseMaker()
+{
+  // Destructor
+}
+
+//________________________________________________________________________
+void AliJetResponseMaker::UserCreateOutputObjects()
+{
+  // Create user objects.
+
+  OpenFile(1);
+  fOutput = new TList();
+  fOutput->SetOwner();
+
+  fHistJetsPt = new TH1F("fHistJetsPt", "fHistJetsPt", fNbins, fMinPt, fMaxPt);
+  fHistJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+  fHistJetsPt->GetYaxis()->SetTitle("counts");
+  fOutput->Add(fHistJetsPt);
+  
+  if (fAnaType == kEMCAL) {
+    fHistJetsPtClus = new TH1F("fHistJetsPtClus", "fHistJetsPtClus", fNbins, fMinPt, fMaxPt);
+    fHistJetsPtClus->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+    fHistJetsPtClus->GetYaxis()->SetTitle("counts");
+    fOutput->Add(fHistJetsPtClus);
+  }
+  
+  fHistJetsPtTrack = new TH1F("fHistJetsPtTrack", "fHistJetsPtTrack", fNbins, fMinPt, fMaxPt);
+  fHistJetsPtTrack->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+  fHistJetsPtTrack->GetYaxis()->SetTitle("counts");
+  fOutput->Add(fHistJetsPtTrack);
+  
+  fHistJetsPtNonBias = new TH1F("fHistJetsPtNonBias", "fHistJetsPtNonBias", fNbins, fMinPt, fMaxPt);
+  fHistJetsPtNonBias->GetXaxis()->SetTitle("p_{T} [GeV/c]");
+  fHistJetsPtNonBias->GetYaxis()->SetTitle("counts");
+  fOutput->Add(fHistJetsPtNonBias);
+  
+  fHistLeadingJetPt = new TH1F("fHistLeadingJetPt", "fHistLeadingJetPt", fNbins, fMinPt, fMaxPt);
+  fHistLeadingJetPt->GetXaxis()->SetTitle("p_{T} [GeV]");
+  fHistLeadingJetPt->GetYaxis()->SetTitle("counts");
+  fOutput->Add(fHistLeadingJetPt);
+  
+  fHistJetsZvsPt = new TH2F("fHistJetsZvsPt", "fHistJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinPt, fMaxPt);
+  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->GetXaxis()->SetTitle("NEF");
+    fHistJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
+    fOutput->Add(fHistJetsNEFvsPt);
+  }
+}
+
+//________________________________________________________________________
+void AliJetResponseMaker::FillHistograms()
+{
+  // Fill histograms.
+
+  Int_t maxJetIndex  = -1;
+  Int_t maxMCJetIndex  = -1;
+
+  DoJetLoop(maxJetIndex, fJets, fTracks, fCaloClusters); 
+  if (maxJetIndex < 0) 
+    return;
+  AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(fJets->At(maxJetIndex));
+  if (!jet) 
+    return;
+
+  DoJetLoop(maxMCJetIndex, fMCJets, fMCParts);  
+  if (maxMCJetIndex < 0) 
+    return;
+  AliEmcalJet* mcjet = dynamic_cast<AliEmcalJet*>(fMCJets->At(maxMCJetIndex));
+  if (!mcjet) 
+    return;
+
+  fHistLeadingJetPt->Fill(jet->Pt());
+  fHistMCLeadingJetPt->Fill(mcjet->Pt());
+}
+
+//________________________________________________________________________
+void AliJetResponseMaker::DoJetLoop(Int_t &maxJetIndex, TClonesArray *jets, TClonesArray *tracks, TClonesArray *clusters)
+{
+  // Do the jet loop.
+
+  if (!fJets)
+    return;
+
+  Int_t njets = jets->GetEntriesFast();
+
+  Float_t maxJetPt = 0;
+  for (Int_t ij = 0; ij < njets; ij++) {
+
+    AliEmcalJet* jet = dynamic_cast<AliEmcalJet*>(jets->At(ij));
+
+    if (!jet) {
+      AliError(Form("Could not receive jet %d", ij));
+      continue;
+    }  
+
+    if (!AcceptJet(jet))
+      continue;
+
+    fHistJetsPtNonBias->Fill(jet->Pt());
+
+    if (jet->MaxTrackPt() > fPtBiasJetTrack)
+      fHistJetsPtTrack->Fill(jet->Pt());
+    
+    if (fAnaType == kEMCAL && jet->MaxClusterPt() > fPtBiasJetClus) 
+      fHistJetsPtClus->Fill(jet->Pt());
+    
+    if (jet->MaxTrackPt() < fPtBiasJetTrack && (fAnaType == kTPC || jet->MaxClusterPt() < fPtBiasJetClus))
+       continue;
+
+    if (maxJetPt < jet->Pt()) {
+      maxJetPt = jet->Pt();
+      maxJetIndex = ij;
+    }
+
+    fHistJetsPt->Fill(jet->Pt());
+
+    fHistJetPhiEta->Fill(jet->Eta(), jet->Phi());
+
+    if (fAnaType == kEMCAL)
+      fHistJetsNEFvsPt->Fill(jet->NEF(), jet->Pt());
+
+    for (Int_t it = 0; it < jet->GetNumberOfTracks(); it++) {
+      AliVParticle *track = jet->TrackAt(it, tracks);
+      if (track)
+       fHistJetsZvsPt->Fill(track->Pt() / jet->Pt(), jet->Pt());
+    }
+
+    if (fAnaType != kEMCAL || !clusters)
+      continue;
+
+    for (Int_t ic = 0; ic < jet->GetNumberOfClusters(); ic++) {
+      AliVCluster *cluster = jet->ClusterAt(ic, clusters);
+
+      if (cluster) {
+       TLorentzVector nPart;
+       cluster->GetMomentum(nPart, fVertex);
+       fHistJetsZvsPt->Fill(nPart.Et() / jet->Pt(), jet->Pt());
+      }
+    }
+  }
+}
+
+//________________________________________________________________________
+void AliJetResponseMaker::RetrieveEventObjects()
+{
+  // Retrieve event objects.
+
+  AliAnalysisTaskEmcal::RetrieveEventObjects();
+  
+  if (strcmp(fMCJetsName,"")) {
+    fMCJets = dynamic_cast<TClonesArray*>(MCEvent()->FindListObject(fMCJetsName));
+    if (!fMCJets) {
+      AliWarning(Form("Could not retrieve MC jets %s!", fMCJetsName.Data())); 
+    }
+  }
+}
+
+//________________________________________________________________________
+void AliJetResponseMaker::Terminate(Option_t *) 
+{
+  // Called once at the end of the analysis.
+}
diff --git a/PWGGA/EMCALJetTasks/AliJetResponseMaker.h b/PWGGA/EMCALJetTasks/AliJetResponseMaker.h
new file mode 100644 (file)
index 0000000..5e387af
--- /dev/null
@@ -0,0 +1,59 @@
+#ifndef ALIJETRESPONSEMAKER_H
+#define ALIJETRESPONSEMAKER_H
+
+// $Id: AliJetResponseMaker.h $
+
+class TClonesArray;
+class TH1F;
+class TH2F;
+
+#include "AliAnalysisTaskEmcal.h"
+
+class AliJetResponseMaker : public AliAnalysisTaskEmcal {
+ public:
+  AliJetResponseMaker();
+  AliJetResponseMaker(const char *name);
+  virtual ~AliJetResponseMaker();
+
+  void         UserCreateOutputObjects();
+  void         Terminate(Option_t *option);
+  void         SetMCPartName(const char *n)     { fMCPartName    = n;  }
+  void         SetMCJetsName(const char *n)     { fMCJetsName    = n;  }
+
+ protected:
+  void                        DoJetLoop(Int_t &maxJetIndex, TClonesArray *jets, TClonesArray *tracks, TClonesArray *clusters = 0);
+  void                        FillHistograms();
+  AliEmcalJet                *GetMCJet(const Int_t i) const;
+  Int_t                       GetNumberOfMCJets()     const;
+  void                        RetrieveEventObjects();
+
+  TString                     fMCPartName;                // name of MC particle collection
+  TString                     fMCJetsName;                // name of MC jet collection
+  TClonesArray               *fMCParts;                   //!MC particles
+  TClonesArray               *fMCJets;                    //!MC jets
+  // Particle level jets
+  TH2F                       *fHistMCJetPhiEta;           //!phi-eta distribution of jets
+  TH1F                       *fHistMCJetsPt;              //!inclusive jet pt spectrum
+  TH1F                       *fHistMCJetsPtTrack;         //!inclusive jet pt spectrum track biased
+  TH1F                       *fHistMCJetsPtClus;          //!inclusive jet pt spectrum cluster biased
+  TH1F                       *fHistMCJetsPtNonBias;       //!non biased inclusive jet pt spectrum
+  TH1F                       *fHistMCLeadingJetPt;        //!leading jet pt spectrum
+  TH2F                       *fHistMCJetsNEFvsPt;         //!jet neutral energy fraction vs. jet pt
+  TH2F                       *fHistMCJetsZvsPt;           //!constituent Pt over Jet Pt ratio vs. jet pt
+  // Detector level jets
+  TH2F                       *fHistJetPhiEta;             //!phi-eta distribution of jets
+  TH1F                       *fHistJetsPt;                //!inclusive jet pt spectrum
+  TH1F                       *fHistJetsPtTrack;           //!inclusive jet pt spectrum track biased
+  TH1F                       *fHistJetsPtClus;            //!inclusive jet pt spectrum cluster biased
+  TH1F                       *fHistJetsPtNonBias;         //!non biased inclusive jet pt spectrum
+  TH1F                       *fHistLeadingJetPt;          //!leading jet pt spectrum
+  TH2F                       *fHistJetsNEFvsPt;           //!jet neutral energy fraction vs. jet pt
+  TH2F                       *fHistJetsZvsPt;             //!constituent Pt over Jet Pt ratio vs. jet pt
+
+ private:
+  AliJetResponseMaker(const AliJetResponseMaker&);            // not implemented
+  AliJetResponseMaker &operator=(const AliJetResponseMaker&); // not implemented
+
+  ClassDef(AliJetResponseMaker, 1) // Jet response matrix producing task
+};
+#endif