return kTRUE;
}
+//________________________________________________________________________
+Int_t AliParticleContainer::GetNAcceptedParticles()
+{
+
+ Int_t nPart = 0;
+
+ AliVParticle *vp = GetNextAcceptParticle(0);
+ if(vp) nPart = 1;
+ while (GetNextAcceptParticle())
+ nPart++;
+
+ return nPart;
+}
+
//________________________________________________________________________
void AliParticleContainer::SetClassName(const char *clname)
{
void GetMomentum(TLorentzVector &mom, Int_t i) const;
Bool_t AcceptParticle(AliVParticle *vp) const;
Int_t GetNParticles() const {return GetNEntries();}
+ Int_t GetNAcceptedParticles() ;
void SetClassName(const char *clname);
protected:
TString runPeriod(runperiod);
Bool_t includeNoITS = kFALSE;
runPeriod.ToLower();
- if (runPeriod == "lhc11h" || runPeriod == "lhc13b" || runPeriod == "lhc13c" || runPeriod == "lhc13d" || runPeriod == "lhc13e" || runPeriod == "lhc13f" || runPeriod == "lhc13g" || runPeriod == "lhc12g") {
+ if (runPeriod == "lhc11h" || runPeriod == "lhc13b" || runPeriod == "lhc13c" || runPeriod == "lhc13d" || runPeriod == "lhc13e" || runPeriod == "lhc13f" || runPeriod == "lhc13g" || runPeriod == "lhc12g" || runPeriod == "lhc10h") {
eTask->SetAODfilterBits(256,512); // hybrid tracks for LHC11h
eTask->SetMC(kFALSE);
+ if(runPeriod == "lhc10h")
+ includeNoITS = kTRUE;
} else if (runPeriod == "lhc12a15e" || runPeriod == "lhc13b4" || runPeriod == "lhc13b4_fix" || runPeriod == "lhc12a15f") {
eTask->SetAODfilterBits(256,512); // hybrid tracks for LHC12a15e, LHC13b4, and LHC12a15f
eTask->SetMC(kTRUE);
- } else if (runPeriod == "lhc11a") {
+ } else if (runPeriod == "lhc11a" || runPeriod == "lhc10hold") {
eTask->SetAODfilterBits(256,16); // hybrid tracks for LHC11a
eTask->SetMC(kFALSE);
includeNoITS = kTRUE;
EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalDiJetResponse.cxx
EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetHMEC.cxx
EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetHadCorQA.cxx
+ EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetMass.cxx
+ EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetMassBkg.cxx
EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetSpectra.cxx
EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetSpectraMECpA.cxx
+ EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetTagger.cxx
EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalJetTriggerQA.cxx
EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalTriggerInfoQA.cxx
EMCALJetTasks/UserTasks/AliAnalysisTaskFullpAJets.cxx
#include "AliRhoParameter.h"
#include "AliLog.h"
#include "AliJetContainer.h"
+#include "AliParticleContainer.h"
+#include "AliClusterContainer.h"
#include "AliAnalysisTaskEmcalJetSample.h"
//________________________________________________________________________
AliAnalysisTaskEmcalJetSample::AliAnalysisTaskEmcalJetSample() :
- AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetSample", kTRUE)
+ AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetSample", kTRUE),
+ fJetsCont(0),
+ fTracksCont(0),
+ fCaloClustersCont(0)
{
// Default constructor.
//________________________________________________________________________
AliAnalysisTaskEmcalJetSample::AliAnalysisTaskEmcalJetSample(const char *name) :
- AliAnalysisTaskEmcalJet(name, kTRUE)
+ AliAnalysisTaskEmcalJet(name, kTRUE),
+ fJetsCont(0),
+ fTracksCont(0),
+ fCaloClustersCont(0)
{
// Standard constructor.
AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
+ fJetsCont = GetJetContainer(0);
+ fTracksCont = fJetsCont->GetParticleContainer();
+ fCaloClustersCont = fJetsCont->GetClusterContainer();
+
TString histname;
for (Int_t i = 0; i < 4; i++) {
{
// Fill histograms.
- if (fTracks) {
- const Int_t ntracks = fTracks->GetEntriesFast();
- for (Int_t it = 0; it < ntracks; it++) {
- AliVTrack *track = static_cast<AliVTrack*>(fTracks->At(it));
-
- if (!track) {
- AliError(Form("Could not receive track %d", it));
- continue;
- }
-
- if (!AcceptTrack(track))
- continue;
-
+ if (fTracksCont) {
+ AliVParticle *track = fTracksCont->GetNextAcceptParticle(0);
+ while(track) {
fHistTracksPt[fCentBin]->Fill(track->Pt());
+
+ track = fTracksCont->GetNextAcceptParticle();
}
}
- if (fCaloClusters) {
- const Int_t nclusters = fCaloClusters->GetEntriesFast();
-
- for (Int_t ic = 0; ic < nclusters; ic++) {
- AliVCluster *cluster = static_cast<AliVCluster*>(fCaloClusters->At(ic));
-
- if (!cluster) {
- AliError(Form("Could not receive cluster %d", ic));
- continue;
- }
+ if (fCaloClustersCont) {
+ AliVCluster *cluster = fCaloClustersCont->GetNextAcceptCluster(0);
+ while(cluster) {
TLorentzVector nPart;
cluster->GetMomentum(nPart, fVertex);
fHistClustersPt[fCentBin]->Fill(nPart.Pt());
+
+ cluster = fCaloClustersCont->GetNextAcceptCluster();
}
}
- if (fJets) {
-
- fJets->Sort();
-
- const Int_t njets = fJets->GetEntriesFast();
- Bool_t leadJet = kFALSE;
- for (Int_t ij = 0; ij < njets; ij++) {
-
- AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(ij));
- if (!jet) {
- AliError(Form("Could not receive jet %d", ij));
- continue;
- }
-
- if (!AcceptJet(jet))
- continue;
-
- if (!leadJet) {
- fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
- leadJet = kTRUE;
- }
+ if (fJetsCont) {
+ AliEmcalJet *jet = fJetsCont->GetNextAcceptJet(0);
+ while(jet) {
fHistJetsPtArea[fCentBin]->Fill(jet->Pt(), jet->Area());
fHistJetsPhiEta[fCentBin]->Fill(jet->Eta(), jet->Phi());
Float_t ptLeading = GetLeadingHadronPt(jet);
fHistJetsPtLeadHad[fCentBin]->Fill(jet->Pt(), ptLeading);
- if (fRho) {
- Float_t corrPt = jet->Pt() - fRhoVal * jet->Area();
- fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());
- }
+ Float_t corrPt = jet->Pt() - fJetsCont->GetRhoVal() * jet->Area();
+ fHistJetsCorrPtArea[fCentBin]->Fill(corrPt, jet->Area());
+
+ jet = fJetsCont->GetNextAcceptJet();
}
+
+ jet = fJetsCont->GetLeadingJet();
+ fHistLeadingJetPt[fCentBin]->Fill(jet->Pt());
}
return kTRUE;
}
+//________________________________________________________________________
+void AliAnalysisTaskEmcalJetSample::ExecOnce() {
+
+ AliAnalysisTaskEmcalJet::ExecOnce();
+
+ if (fJetsCont && fJetsCont->GetArray() == 0) fJetsCont = 0;
+ if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
+ if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
+
+}
+
//________________________________________________________________________
Bool_t AliAnalysisTaskEmcalJetSample::Run()
{
-#ifndef ALIANALYSISTASKEMCALJETSAMPLEDEV_H
-#define ALIANALYSISTASKEMCALJETSAMPLEDEV_H
+#ifndef ALIANALYSISTASKEMCALJETSAMPLE_H
+#define ALIANALYSISTASKEMCALJETSAMPLE_H
// $Id$
class TH1;
class TH2;
+class AliJetContainer;
+class AliParticleContainer;
+class AliClusterContainer;
#include "AliAnalysisTaskEmcalJet.h"
void Terminate(Option_t *option);
protected:
+ void ExecOnce();
Bool_t FillHistograms() ;
Bool_t Run() ;
TH2 *fHistJetsPtLeadHad[4]; //!Jet pt vs. leading hadron
TH2 *fHistJetsCorrPtArea[4]; //!Jet pt - bkg vs. area
+ AliJetContainer *fJetsCont; //!Jets
+ AliParticleContainer *fTracksCont; //!Tracks
+ AliClusterContainer *fCaloClustersCont; //!Clusters
+
private:
AliAnalysisTaskEmcalJetSample(const AliAnalysisTaskEmcalJetSample&); // not implemented
AliAnalysisTaskEmcalJetSample &operator=(const AliAnalysisTaskEmcalJetSample&); // not implemented
- ClassDef(AliAnalysisTaskEmcalJetSample, 1) // jet sample analysis task
+ ClassDef(AliAnalysisTaskEmcalJetSample, 2) // jet sample analysis task
};
#endif
void SetJetPhiRange(Double_t pmi, Double_t pma) {fJetPhiMin = pmi; fJetPhiMax = pma; }
void SetGhostArea(Double_t gharea) { fGhostArea = gharea; }
void SetMinMCLabel(Int_t s) { fMinMCLabel = s ; }
- void SetRecombSheme(Int_t scheme) { fRecombScheme = scheme; }
+ void SetRecombScheme(Int_t scheme) { fRecombScheme = scheme; }
void SelectCollisionCandidates(UInt_t offlineTriggerMask = AliVEvent::kMB)
{
if(!fIsPSelSet)
-// $Id: AliAnalysisTaskCLQA.cxx 60694 2013-02-04 15:35:56Z morsch $
+// $Id$
//
// Constantin's Task
//
#ifndef ALIANALYSISTASKCLQA_H
#define ALIANALYSISTASKCLQA_H
-// $Id: AliAnalysisTaskCLQA.h 58847 2012-09-30 18:11:49Z loizides $
+// $Id$
class TClonesArray;
class TString;
--- /dev/null
+//
+// Jet mass analysis task.
+//
+// Author: M.Verweij
+
+#include <TClonesArray.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <THnSparse.h>
+#include <TList.h>
+#include <TLorentzVector.h>
+#include <TProfile.h>
+#include <TChain.h>
+#include <TSystem.h>
+#include <TFile.h>
+#include <TKey.h>
+
+#include "AliVCluster.h"
+#include "AliVTrack.h"
+#include "AliEmcalJet.h"
+#include "AliRhoParameter.h"
+#include "AliLog.h"
+#include "AliEmcalParticle.h"
+#include "AliMCEvent.h"
+#include "AliGenPythiaEventHeader.h"
+#include "AliAODMCHeader.h"
+#include "AliMCEvent.h"
+#include "AliAnalysisManager.h"
+#include "AliJetContainer.h"
+
+#include "AliAODEvent.h"
+
+#include "AliAnalysisTaskEmcalJetMass.h"
+
+ClassImp(AliAnalysisTaskEmcalJetMass)
+
+//________________________________________________________________________
+AliAnalysisTaskEmcalJetMass::AliAnalysisTaskEmcalJetMass() :
+ AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetMass", kTRUE),
+ fContainerBase(0),
+ fh2PtJet1VsLeadPtAllSel(0),
+ fh2PtJet1VsLeadPtTagged(0),
+ fh2PtVsMassJet1All(0),
+ fh2PtVsMassJet1Tagged(0),
+ fpPtVsMassJet1All(0),
+ fpPtVsMassJet1Tagged(0),
+ fh2MassVsAreaJet1All(0),
+ fh2MassVsAreaJet1Tagged(0)
+{
+ // Default constructor.
+
+ fh2PtJet1VsLeadPtAllSel = new TH2F*[fNcentBins];
+ fh2PtJet1VsLeadPtTagged = new TH2F*[fNcentBins];
+ fh2PtVsMassJet1All = new TH2F*[fNcentBins];
+ fh2PtVsMassJet1Tagged = new TH2F*[fNcentBins];
+ fpPtVsMassJet1All = new TProfile*[fNcentBins];
+ fpPtVsMassJet1Tagged = new TProfile*[fNcentBins];
+ fh2MassVsAreaJet1All = new TH2F*[fNcentBins];
+ fh2MassVsAreaJet1Tagged = new TH2F*[fNcentBins];
+
+ for (Int_t i = 0; i < fNcentBins; i++) {
+ fh2PtJet1VsLeadPtAllSel[i] = 0;
+ fh2PtJet1VsLeadPtTagged[i] = 0;
+ fh2PtVsMassJet1All[i] = 0;
+ fh2PtVsMassJet1Tagged[i] = 0;
+ fpPtVsMassJet1All[i] = 0;
+ fpPtVsMassJet1Tagged[i] = 0;
+ fh2MassVsAreaJet1All[i] = 0;
+ fh2MassVsAreaJet1Tagged[i] = 0;
+ }
+
+ SetMakeGeneralHistograms(kTRUE);
+
+}
+
+//________________________________________________________________________
+AliAnalysisTaskEmcalJetMass::AliAnalysisTaskEmcalJetMass(const char *name) :
+ AliAnalysisTaskEmcalJet(name, kTRUE),
+ fContainerBase(0),
+ fh2PtJet1VsLeadPtAllSel(0),
+ fh2PtJet1VsLeadPtTagged(0),
+ fh2PtVsMassJet1All(0),
+ fh2PtVsMassJet1Tagged(0),
+ fpPtVsMassJet1All(0),
+ fpPtVsMassJet1Tagged(0),
+ fh2MassVsAreaJet1All(0),
+ fh2MassVsAreaJet1Tagged(0)
+{
+ // Standard constructor.
+
+ fh2PtJet1VsLeadPtAllSel = new TH2F*[fNcentBins];
+ fh2PtJet1VsLeadPtTagged = new TH2F*[fNcentBins];
+ fh2PtVsMassJet1All = new TH2F*[fNcentBins];
+ fh2PtVsMassJet1Tagged = new TH2F*[fNcentBins];
+ fpPtVsMassJet1All = new TProfile*[fNcentBins];
+ fpPtVsMassJet1Tagged = new TProfile*[fNcentBins];
+ fh2MassVsAreaJet1All = new TH2F*[fNcentBins];
+ fh2MassVsAreaJet1Tagged = new TH2F*[fNcentBins];
+
+ for (Int_t i = 0; i < fNcentBins; i++) {
+ fh2PtJet1VsLeadPtAllSel[i] = 0;
+ fh2PtJet1VsLeadPtTagged[i] = 0;
+ fh2PtVsMassJet1All[i] = 0;
+ fh2PtVsMassJet1Tagged[i] = 0;
+ fpPtVsMassJet1All[i] = 0;
+ fpPtVsMassJet1Tagged[i] = 0;
+ fh2MassVsAreaJet1All[i] = 0;
+ fh2MassVsAreaJet1Tagged[i] = 0;
+ }
+
+ SetMakeGeneralHistograms(kTRUE);
+}
+
+//________________________________________________________________________
+AliAnalysisTaskEmcalJetMass::~AliAnalysisTaskEmcalJetMass()
+{
+ // Destructor.
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEmcalJetMass::UserCreateOutputObjects()
+{
+ // Create user output.
+
+ AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
+
+ Bool_t oldStatus = TH1::AddDirectoryStatus();
+ TH1::AddDirectory(kFALSE);
+
+ const Int_t nBinsPt = 250;
+ const Double_t minPt = -50.;
+ const Double_t maxPt = 200.;
+
+ const Int_t nBinsArea = 100;
+ const Double_t minArea = 0.;
+ const Double_t maxArea = 1.;
+
+ TString histName = "";
+ TString histTitle = "";
+ for (Int_t i = 0; i < fNcentBins; i++) {
+ histName = TString::Format("fh2PtJet1VsLeadPtAllSel_%d",i);
+ histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
+ fh2PtJet1VsLeadPtAllSel[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
+ fOutput->Add(fh2PtJet1VsLeadPtAllSel[i]);
+
+ histName = TString::Format("fh2PtJet1VsLeadPtTagged_%d",i);
+ histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
+ fh2PtJet1VsLeadPtTagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
+ fOutput->Add(fh2PtJet1VsLeadPtTagged[i]);
+
+ histName = TString::Format("fh2PtVsMassJet1All_%d",i);
+ histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
+ fh2PtVsMassJet1All[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
+ fOutput->Add(fh2PtVsMassJet1All[i]);
+
+ histName = TString::Format("fh2PtVsMassJet1Tagged_%d",i);
+ histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{M}_{jet1}",histName.Data());
+ fh2PtVsMassJet1Tagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
+ fOutput->Add(fh2PtVsMassJet1Tagged[i]);
+
+ histName = TString::Format("fpPtVsMassJet1All_%d",i);
+ histTitle = TString::Format("%s;#it{p}_{T,jet1};Avg #it{M}_{jet1}",histName.Data());
+ fpPtVsMassJet1All[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
+ fOutput->Add(fpPtVsMassJet1All[i]);
+
+ histName = TString::Format("fpPtVsMassJet1Tagged_%d",i);
+ histTitle = TString::Format("%s;#it{p}_{T,jet1};Avg #it{M}_{jet1}",histName.Data());
+ fpPtVsMassJet1Tagged[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
+ fOutput->Add(fpPtVsMassJet1Tagged[i]);
+
+ histName = TString::Format("fh2MassVsAreaJet1All_%d",i);
+ histTitle = TString::Format("%s;#it{M}_{jet1};#it{A}",histName.Data());
+ fh2MassVsAreaJet1All[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsArea,minArea,maxArea);
+ fOutput->Add(fh2MassVsAreaJet1All[i]);
+
+ histName = TString::Format("fh2MassVsAreaJet1Tagged_%d",i);
+ histTitle = TString::Format("%s;#it{M}_{jet1};#it{A}",histName.Data());
+ fh2MassVsAreaJet1Tagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsArea,minArea,maxArea);
+ fOutput->Add(fh2MassVsAreaJet1Tagged[i]);
+
+ }
+
+ // =========== Switch on Sumw2 for all histos ===========
+ for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
+ TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
+ if (h1){
+ h1->Sumw2();
+ continue;
+ }
+ THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
+ if(hn)hn->Sumw2();
+ }
+
+ TH1::AddDirectory(oldStatus);
+
+ PostData(1, fOutput); // Post data for ALL output slots > 0 here.
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskEmcalJetMass::Run()
+{
+ // Run analysis code here, if needed. It will be executed before FillHistograms().
+
+ return kTRUE;
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskEmcalJetMass::FillHistograms()
+{
+ // Fill histograms.
+
+ for(int i = 0; i < GetNJets(fContainerBase);++i) {
+ AliEmcalJet *jet1 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(i, fContainerBase));
+ if(!jet1) continue;
+
+ Double_t ptJet1 = jet1->Pt() - GetRhoVal(fContainerBase)*jet1->Area();
+ fh2PtJet1VsLeadPtAllSel[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
+ fh2PtVsMassJet1All[fCentBin]->Fill(ptJet1,jet1->M());
+ fpPtVsMassJet1All[fCentBin]->Fill(ptJet1,jet1->M());
+ fh2MassVsAreaJet1All[fCentBin]->Fill(jet1->M(),jet1->Area());
+
+ if(jet1->GetTagStatus()<1)
+ continue;
+
+ AliEmcalJet *jet2 = jet1->GetTaggedJet();
+ if(!jet2) continue;
+
+ fh2PtJet1VsLeadPtTagged[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
+ fh2PtVsMassJet1Tagged[fCentBin]->Fill(ptJet1,jet1->M());
+ fpPtVsMassJet1Tagged[fCentBin]->Fill(ptJet1,jet1->M());
+ fh2MassVsAreaJet1Tagged[fCentBin]->Fill(jet1->M(),jet1->Area());
+
+ }
+
+ return kTRUE;
+
+}
+
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskEmcalJetMass::RetrieveEventObjects() {
+ //
+ // retrieve event objects
+ //
+
+ if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
+ return kFALSE;
+
+ return kTRUE;
+
+}
+
+//_______________________________________________________________________
+void AliAnalysisTaskEmcalJetMass::Terminate(Option_t *)
+{
+ // Called once at the end of the analysis.
+}
+
--- /dev/null
+#ifndef ALIANALYSISTASKEMCALJETMASS_H
+#define ALIANALYSISTASKEMCALJETMASS_H
+
+class TH1;
+class TH2;
+class TH3;
+class TH3F;
+class THnSparse;
+class TClonesArray;
+class TArrayI;
+class AliAnalysisManager;
+class AliJetContainer;
+
+#include "AliAnalysisTaskEmcalJet.h"
+
+class AliAnalysisTaskEmcalJetMass : public AliAnalysisTaskEmcalJet {
+ public:
+
+ AliAnalysisTaskEmcalJetMass();
+ AliAnalysisTaskEmcalJetMass(const char *name);
+ virtual ~AliAnalysisTaskEmcalJetMass();
+
+ void UserCreateOutputObjects();
+ void Terminate(Option_t *option);
+
+ //Setters
+ void SetJetContainerBase(Int_t c) { fContainerBase = c;}
+
+ protected:
+ Bool_t RetrieveEventObjects();
+ Bool_t Run();
+ Bool_t FillHistograms();
+
+ Int_t fContainerBase; // jets to be tagged
+
+ private:
+ TH2F **fh2PtJet1VsLeadPtAllSel; //!all jets after std selection vs leading track pt
+ TH2F **fh2PtJet1VsLeadPtTagged; //!tagged jets vs leading track pt
+ TH2F **fh2PtVsMassJet1All; //!pT vs mass of all jets
+ TH2F **fh2PtVsMassJet1Tagged; //!pT vs mass of tagged jets
+ TProfile **fpPtVsMassJet1All; //!pT vs avg mass of all jets
+ TProfile **fpPtVsMassJet1Tagged; //!pT vs avg mass of tagged jets
+ TH2F **fh2MassVsAreaJet1All; //!mass vs area of all jets
+ TH2F **fh2MassVsAreaJet1Tagged; //!mass vs area of tagged jets
+
+ AliAnalysisTaskEmcalJetMass(const AliAnalysisTaskEmcalJetMass&); // not implemented
+ AliAnalysisTaskEmcalJetMass &operator=(const AliAnalysisTaskEmcalJetMass&); // not implemented
+
+ ClassDef(AliAnalysisTaskEmcalJetMass, 1)
+};
+#endif
+
--- /dev/null
+//
+// Jet mass background analysis task.
+//
+// Author: M.Verweij
+
+#include <TClonesArray.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <THnSparse.h>
+#include <TList.h>
+#include <TLorentzVector.h>
+#include <TProfile.h>
+#include <TChain.h>
+#include <TSystem.h>
+#include <TFile.h>
+#include <TKey.h>
+#include <TRandom3.h>
+
+#include "AliVCluster.h"
+#include "AliVTrack.h"
+#include "AliEmcalJet.h"
+#include "AliRhoParameter.h"
+#include "AliLog.h"
+#include "AliEmcalParticle.h"
+#include "AliMCEvent.h"
+#include "AliGenPythiaEventHeader.h"
+#include "AliAODMCHeader.h"
+#include "AliMCEvent.h"
+#include "AliAnalysisManager.h"
+#include "AliJetContainer.h"
+#include "AliClusterContainer.h"
+#include "AliParticleContainer.h"
+
+#include "AliAnalysisTaskEmcalJetMassBkg.h"
+
+ClassImp(AliAnalysisTaskEmcalJetMassBkg)
+
+//________________________________________________________________________
+AliAnalysisTaskEmcalJetMassBkg::AliAnalysisTaskEmcalJetMassBkg() :
+ AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetMassBkg", kTRUE),
+ fContainerBase(0),
+ fMinRC2LJ(-1),
+ fRCperEvent(10),
+ fConeRadius(0.2),
+ fConeMinEta(-0.9),
+ fConeMaxEta(0.9),
+ fConeMinPhi(0),
+ fConeMaxPhi(TMath::Pi()*2),
+ fJetsCont(0),
+ fTracksCont(0),
+ fCaloClustersCont(0),
+ fh2PtVsMassRC(0),
+ fpPtVsMassRC(0),
+ fh2PtVsMassRCExLJDPhi(0),
+ fpPtVsMassRCExLJ(0),
+ fh2PtVsMassPerpConeLJ(0),
+ fpPtVsMassPerpConeLJ(0),
+ fh2PtVsMassPerpConeTJ(0),
+ fpPtVsMassPerpConeTJ(0),
+ fh2CentVsMassRC(0),
+ fh2CentVsMassRCExLJ(0),
+ fh2CentVsMassPerpConeLJ(0),
+ fh2CentVsMassPerpConeTJ(0),
+ fh2MultVsMassRC(0),
+ fh2MultVsMassRCExLJ(0),
+ fh2MultVsMassPerpConeLJ(0),
+ fh2MultVsMassPerpConeTJ(0)
+{
+ // Default constructor.
+
+ fh2PtVsMassRC = new TH2F*[fNcentBins];
+ fpPtVsMassRC = new TProfile*[fNcentBins];
+ fh2PtVsMassRCExLJDPhi = new TH3F*[fNcentBins];
+ fpPtVsMassRCExLJ = new TProfile*[fNcentBins];
+ fh2PtVsMassPerpConeLJ = new TH2F*[fNcentBins];
+ fpPtVsMassPerpConeLJ = new TProfile*[fNcentBins];
+ fh2PtVsMassPerpConeTJ = new TH2F*[fNcentBins];
+ fpPtVsMassPerpConeTJ = new TProfile*[fNcentBins];
+
+ for (Int_t i = 0; i < fNcentBins; i++) {
+ fh2PtVsMassRC[i] = 0;
+ fpPtVsMassRC[i] = 0;
+ fh2PtVsMassRCExLJDPhi[i] = 0;
+ fpPtVsMassRCExLJ[i] = 0;
+ fh2PtVsMassPerpConeLJ[i] = 0;
+ fpPtVsMassPerpConeLJ[i] = 0;
+ fh2PtVsMassPerpConeTJ[i] = 0;
+ fpPtVsMassPerpConeTJ[i] = 0;
+ }
+
+ SetMakeGeneralHistograms(kTRUE);
+
+}
+
+//________________________________________________________________________
+AliAnalysisTaskEmcalJetMassBkg::AliAnalysisTaskEmcalJetMassBkg(const char *name) :
+ AliAnalysisTaskEmcalJet(name, kTRUE),
+ fContainerBase(0),
+ fMinRC2LJ(-1),
+ fRCperEvent(10),
+ fConeRadius(0.2),
+ fConeMinEta(-0.9),
+ fConeMaxEta(0.9),
+ fConeMinPhi(0),
+ fConeMaxPhi(TMath::Pi()*2),
+ fJetsCont(0),
+ fTracksCont(0),
+ fCaloClustersCont(0),
+ fh2PtVsMassRC(0),
+ fpPtVsMassRC(0),
+ fh2PtVsMassRCExLJDPhi(0),
+ fpPtVsMassRCExLJ(0),
+ fh2PtVsMassPerpConeLJ(0),
+ fpPtVsMassPerpConeLJ(0),
+ fh2PtVsMassPerpConeTJ(0),
+ fpPtVsMassPerpConeTJ(0),
+ fh2CentVsMassRC(0),
+ fh2CentVsMassRCExLJ(0),
+ fh2CentVsMassPerpConeLJ(0),
+ fh2CentVsMassPerpConeTJ(0),
+ fh2MultVsMassRC(0),
+ fh2MultVsMassRCExLJ(0),
+ fh2MultVsMassPerpConeLJ(0),
+ fh2MultVsMassPerpConeTJ(0)
+{
+ // Standard constructor.
+
+ fh2PtVsMassRC = new TH2F*[fNcentBins];
+ fpPtVsMassRC = new TProfile*[fNcentBins];
+ fh2PtVsMassRCExLJDPhi = new TH3F*[fNcentBins];
+ fpPtVsMassRCExLJ = new TProfile*[fNcentBins];
+ fh2PtVsMassPerpConeLJ = new TH2F*[fNcentBins];
+ fpPtVsMassPerpConeLJ = new TProfile*[fNcentBins];
+ fh2PtVsMassPerpConeTJ = new TH2F*[fNcentBins];
+ fpPtVsMassPerpConeTJ = new TProfile*[fNcentBins];
+
+ for (Int_t i = 0; i < fNcentBins; i++) {
+ fh2PtVsMassRC[i] = 0;
+ fpPtVsMassRC[i] = 0;
+ fh2PtVsMassRCExLJDPhi[i] = 0;
+ fpPtVsMassRCExLJ[i] = 0;
+ fh2PtVsMassPerpConeLJ[i] = 0;
+ fpPtVsMassPerpConeLJ[i] = 0;
+ fh2PtVsMassPerpConeTJ[i] = 0;
+ fpPtVsMassPerpConeTJ[i] = 0;
+ }
+
+ SetMakeGeneralHistograms(kTRUE);
+}
+
+//________________________________________________________________________
+AliAnalysisTaskEmcalJetMassBkg::~AliAnalysisTaskEmcalJetMassBkg()
+{
+ // Destructor.
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEmcalJetMassBkg::UserCreateOutputObjects()
+{
+ // Create user output.
+
+ AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
+
+ fJetsCont = GetJetContainer(fContainerBase);
+ fTracksCont = fJetsCont->GetParticleContainer();
+ fCaloClustersCont = fJetsCont->GetClusterContainer();
+
+ Bool_t oldStatus = TH1::AddDirectoryStatus();
+ TH1::AddDirectory(kFALSE);
+
+ const Int_t nBinsPt = 250;
+ const Double_t minPt = -50.;
+ const Double_t maxPt = 200.;
+
+ const Int_t nBinsCent = 100;
+ const Double_t minCent = 0.;
+ const Double_t maxCent = 100.;
+
+ const Int_t nBinsMult = 400;
+ const Double_t minMult = 0.;
+ const Double_t maxMult = 4000.;
+
+ fh2CentVsMassRC = new TH2F("fh2CentVsMassRC","fh2CentVsMassRC;cent;#it{M}_{RC}",nBinsCent,minCent,maxCent,nBinsPt,minPt,maxPt);
+ fOutput->Add(fh2CentVsMassRC);
+
+ fh2CentVsMassRCExLJ = new TH2F("fh2CentVsMassRCExLJ","fh2CentVsMassRCExLJ;cent;#it{M}_{RC}",nBinsCent,minCent,maxCent,nBinsPt,minPt,maxPt);
+ fOutput->Add(fh2CentVsMassRCExLJ);
+
+ fh2CentVsMassPerpConeLJ = new TH2F("fh2CentVsMassPerpConeLJ","fh2CentVsMassPerpConeLJ;cent;#it{M}_{RC}",nBinsCent,minCent,maxCent,nBinsPt,minPt,maxPt);
+ fOutput->Add(fh2CentVsMassPerpConeLJ);
+
+ fh2CentVsMassPerpConeTJ = new TH2F("fh2CentVsMassPerpConeTJ","fh2CentVsMassPerpConeTJ;cent;#it{M}_{RC}",nBinsCent,minCent,maxCent,nBinsPt,minPt,maxPt);
+ fOutput->Add(fh2CentVsMassPerpConeTJ);
+
+ fh2MultVsMassRC = new TH2F("fh2MultVsMassRC","fh2MultVsMassRC;#it{N}_{track};#it{M}_{RC}",nBinsMult,minMult,maxMult,nBinsPt,minPt,maxPt);
+ fOutput->Add(fh2MultVsMassRC);
+
+ fh2MultVsMassRCExLJ = new TH2F("fh2MultVsMassRCExLJ","fh2MultVsMassRCExLJ;#it{N}_{track};#it{M}_{RC}",nBinsMult,minMult,maxMult,nBinsPt,minPt,maxPt);
+ fOutput->Add(fh2MultVsMassRCExLJ);
+
+ fh2MultVsMassPerpConeLJ = new TH2F("fh2MultVsMassPerpConeLJ","fh2MultVsMassPerpConeLJ;#it{N}_{track};#it{M}_{RC}",nBinsMult,minMult,maxMult,nBinsPt,minPt,maxPt);
+ fOutput->Add(fh2MultVsMassPerpConeLJ);
+
+ fh2MultVsMassPerpConeTJ = new TH2F("fh2MultVsMassPerpConeTJ","fh2MultVsMassPerpConeTJ;#it{N}_{track};#it{M}_{RC}",nBinsMult,minMult,maxMult,nBinsPt,minPt,maxPt);
+ fOutput->Add(fh2MultVsMassPerpConeTJ);
+
+ TString histName = "";
+ TString histTitle = "";
+ for (Int_t i = 0; i < fNcentBins; i++) {
+ histName = TString::Format("fh2PtVsMassRC_%d",i);
+ histTitle = TString::Format("%s;#it{p}_{T,RC};#it{M}_{RC}",histName.Data());
+ fh2PtVsMassRC[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
+ fOutput->Add(fh2PtVsMassRC[i]);
+
+ histName = TString::Format("fpPtVsMassRC_%d",i);
+ histTitle = TString::Format("%s;#it{p}_{T,RC};Avg #it{M}_{RC}",histName.Data());
+ fpPtVsMassRC[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
+ fOutput->Add(fpPtVsMassRC[i]);
+
+ histName = TString::Format("fh2PtVsMassRCExLJDPhi_%d",i);
+ histTitle = TString::Format("%s;#it{p}_{T,RC};#it{M}_{RC}",histName.Data());
+ fh2PtVsMassRCExLJDPhi[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt,72,-0.5*TMath::Pi(),1.5*TMath::Pi());
+ fOutput->Add(fh2PtVsMassRCExLJDPhi[i]);
+
+ histName = TString::Format("fpPtVsMassRCExLJ_%d",i);
+ histTitle = TString::Format("%s;#it{p}_{T,RC};Avg #it{M}_{RC}",histName.Data());
+ fpPtVsMassRCExLJ[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
+ fOutput->Add(fpPtVsMassRCExLJ[i]);
+
+ histName = TString::Format("fh2PtVsMassPerpConeLJ_%d",i);
+ histTitle = TString::Format("%s;#it{p}_{T,PerpConeLJ};#it{M}_{PerpConeLJ}",histName.Data());
+ fh2PtVsMassPerpConeLJ[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
+ fOutput->Add(fh2PtVsMassPerpConeLJ[i]);
+
+ histName = TString::Format("fpPtVsMassPerpConeLJ_%d",i);
+ histTitle = TString::Format("%s;#it{p}_{T,RC};Avg #it{M}_{RC}",histName.Data());
+ fpPtVsMassPerpConeLJ[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
+ fOutput->Add(fpPtVsMassPerpConeLJ[i]);
+
+ histName = TString::Format("fh2PtVsMassPerpConeTJ_%d",i);
+ histTitle = TString::Format("%s;#it{p}_{T,PerpConeTJ};#it{M}_{PerpConeTJ}",histName.Data());
+ fh2PtVsMassPerpConeTJ[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
+ fOutput->Add(fh2PtVsMassPerpConeTJ[i]);
+
+ histName = TString::Format("fpPtVsMassPerpConeTJ_%d",i);
+ histTitle = TString::Format("%s;#it{p}_{T,RC};Avg #it{M}_{RC}",histName.Data());
+ fpPtVsMassPerpConeTJ[i] = new TProfile(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt);
+ fOutput->Add(fpPtVsMassPerpConeTJ[i]);
+ }
+
+ // =========== Switch on Sumw2 for all histos ===========
+ for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
+ TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
+ if (h1){
+ h1->Sumw2();
+ continue;
+ }
+ THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
+ if(hn)hn->Sumw2();
+ }
+
+ TH1::AddDirectory(oldStatus);
+
+ PostData(1, fOutput); // Post data for ALL output slots > 0 here.
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskEmcalJetMassBkg::Run()
+{
+ // Run analysis code here, if needed. It will be executed before FillHistograms().
+
+ return kTRUE;
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskEmcalJetMassBkg::FillHistograms()
+{
+ // Fill histograms.
+
+ const Float_t rcArea = fConeRadius * fConeRadius * TMath::Pi();
+ Double_t rho = GetRhoVal(fContainerBase);
+ Int_t trackMult = fTracksCont->GetNAcceptedParticles();
+
+ TLorentzVector lvRC(0.,0.,0.,0.);
+ Float_t RCpt = 0;
+ Float_t RCeta = 0;
+ Float_t RCphi = 0;
+ Float_t RCmass = 0.;
+ for (Int_t i = 0; i < fRCperEvent; i++) {
+ // Simple random cones
+ lvRC.SetPxPyPzE(0.,0.,0.,0.);
+ RCpt = 0;
+ RCeta = 0;
+ RCphi = 0;
+ GetRandomCone(lvRC,RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, 0);
+ RCmass = lvRC.M();
+ if (RCpt > 0) {
+ fh2PtVsMassRC[fCentBin]->Fill(RCpt - rho*rcArea,RCmass);
+ fpPtVsMassRC[fCentBin]->Fill(RCpt - rho*rcArea,RCmass);
+ fh2CentVsMassRC->Fill(fCent,RCmass);
+ fh2MultVsMassRC->Fill(trackMult,RCmass);
+ }
+
+ if (fJetsCont) {
+
+ // Random cones far from leading jet(s)
+ AliEmcalJet* jet = fJetsCont->GetLeadingJet("rho");
+ lvRC.SetPxPyPzE(0.,0.,0.,0.);
+ RCpt = 0;
+ RCeta = 0;
+ RCphi = 0;
+ GetRandomCone(lvRC,RCpt, RCeta, RCphi, fTracksCont, fCaloClustersCont, jet);
+ RCmass = lvRC.M();
+ if (RCpt > 0 && jet) {
+ Float_t dphi = RCphi - jet->Phi();
+ if (dphi > 1.5*TMath::Pi()) dphi -= TMath::Pi() * 2;
+ if (dphi < -0.5*TMath::Pi()) dphi += TMath::Pi() * 2;
+ fh2PtVsMassRCExLJDPhi[fCentBin]->Fill(RCpt - rho*rcArea,RCmass,dphi);
+ fpPtVsMassRCExLJ[fCentBin]->Fill(RCpt - rho*rcArea,RCmass);
+ fh2CentVsMassRCExLJ->Fill(fCent,RCmass);
+ fh2MultVsMassRCExLJ->Fill(trackMult,RCmass);
+ }
+ }
+ }//RC loop
+
+ if(fJetsCont) {
+ //cone perpendicular to leading jet
+ TLorentzVector lvPC(0.,0.,0.,0.);
+ Float_t PCpt = 0;
+ Float_t PCeta = 0;
+ Float_t PCphi = 0;
+ Float_t PCmass = 0.;
+ AliEmcalJet* jet = fJetsCont->GetLeadingJet("rho");
+ if(jet) {
+ GetPerpCone(lvPC,PCpt, PCeta, PCphi, fTracksCont, fCaloClustersCont, jet);
+ PCmass = lvPC.M();
+ if(PCpt>0.) {
+ fh2PtVsMassPerpConeLJ[fCentBin]->Fill(PCpt-rho*rcArea,PCmass);
+ fpPtVsMassPerpConeLJ[fCentBin]->Fill(PCpt-rho*rcArea,PCmass);
+ fh2CentVsMassPerpConeLJ->Fill(fCent,PCmass);
+ fh2MultVsMassPerpConeLJ->Fill(trackMult,PCmass);
+ }
+ }
+ //cone perpendicular to all tagged jets
+ for(int i = 0; i < fJetsCont->GetNJets();++i) {
+ jet = static_cast<AliEmcalJet*>(fJetsCont->GetAcceptJet(i));
+ if(!jet) continue;
+
+ if(jet->GetTagStatus()<1)
+ continue;
+
+ lvPC.SetPxPyPzE(0.,0.,0.,0.);
+ PCpt = 0;
+ PCeta = 0;
+ PCphi = 0;
+ GetPerpCone(lvPC,PCpt, PCeta, PCphi, fTracksCont, fCaloClustersCont, jet);
+ PCmass = lvPC.M();
+ if(PCpt>0.) {
+ fh2PtVsMassPerpConeTJ[fCentBin]->Fill(PCpt-rho*rcArea,PCmass);
+ fpPtVsMassPerpConeTJ[fCentBin]->Fill(PCpt-rho*rcArea,PCmass);
+ fh2CentVsMassPerpConeTJ->Fill(fCent,PCmass);
+ fh2MultVsMassPerpConeTJ->Fill(trackMult,PCmass);
+ }
+ }//jet loop
+ }
+
+ return kTRUE;
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEmcalJetMassBkg::GetRandomCone(TLorentzVector& lvRC,Float_t &pt, Float_t &eta, Float_t &phi,
+ AliParticleContainer* tracks, AliClusterContainer* clusters,
+ AliEmcalJet *jet) const
+{
+ // Get rigid cone.
+ lvRC.SetPxPyPzE(0.,0.,0.,0.);
+
+ eta = -999;
+ phi = -999;
+ pt = 0;
+
+ if (!tracks && !clusters)
+ return;
+
+ Float_t LJeta = 999;
+ Float_t LJphi = 999;
+
+ if (jet) {
+ LJeta = jet->Eta();
+ LJphi = jet->Phi();
+ }
+
+ Float_t maxEta = fConeMaxEta;
+ Float_t minEta = fConeMinEta;
+ Float_t maxPhi = fConeMaxPhi;
+ Float_t minPhi = fConeMinPhi;
+
+ if (maxPhi > TMath::Pi() * 2) maxPhi = TMath::Pi() * 2;
+ if (minPhi < 0) minPhi = 0;
+
+ Float_t dLJ = 0;
+ Int_t repeats = 0;
+ Bool_t reject = kTRUE;
+ do {
+ eta = gRandom->Rndm() * (maxEta - minEta) + minEta;
+ phi = gRandom->Rndm() * (maxPhi - minPhi) + minPhi;
+ dLJ = TMath::Sqrt((LJeta - eta) * (LJeta - eta) + (LJphi - phi) * (LJphi - phi));
+
+ repeats++;
+ } while (dLJ < fMinRC2LJ && repeats < 999 && reject);
+
+ if (repeats == 999) {
+ AliWarning(Form("%s: Could not get random cone!", GetName()));
+ return;
+ }
+
+ GetCone(lvRC,pt,eta,phi,tracks,clusters);
+
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEmcalJetMassBkg::GetCone(TLorentzVector& lvRC,Float_t &pt, Float_t eta, Float_t phi, AliParticleContainer* tracks, AliClusterContainer* clusters) const
+{
+
+ pt = 0.;
+ lvRC.SetPxPyPzE(0.,0.,0.,0.);
+
+ if (clusters) {
+ AliVCluster* cluster = clusters->GetNextAcceptCluster(0);
+ while (cluster) {
+ TLorentzVector nPart;
+ cluster->GetMomentum(nPart, const_cast<Double_t*>(fVertex));
+
+ Float_t cluseta = nPart.Eta();
+ Float_t clusphi = nPart.Phi();
+
+ if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi + 2 * TMath::Pi()))
+ clusphi += 2 * TMath::Pi();
+ if (TMath::Abs(clusphi - phi) > TMath::Abs(clusphi - phi - 2 * TMath::Pi()))
+ clusphi -= 2 * TMath::Pi();
+
+ Float_t d = TMath::Sqrt((cluseta - eta) * (cluseta - eta) + (clusphi - phi) * (clusphi - phi));
+ if (d <= fConeRadius) {
+ pt += nPart.Pt();
+ TLorentzVector lvcl(nPart.Px(),nPart.Py(),nPart.Pz(),nPart.E());
+ lvRC += lvcl;
+ }
+
+ cluster = clusters->GetNextAcceptCluster();
+ }
+ }
+
+ if (tracks) {
+ AliVParticle* track = tracks->GetNextAcceptParticle(0);
+ while(track) {
+ Float_t tracketa = track->Eta();
+ Float_t trackphi = track->Phi();
+
+ if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi + 2 * TMath::Pi()))
+ trackphi += 2 * TMath::Pi();
+ if (TMath::Abs(trackphi - phi) > TMath::Abs(trackphi - phi - 2 * TMath::Pi()))
+ trackphi -= 2 * TMath::Pi();
+
+ Float_t d = TMath::Sqrt((tracketa - eta) * (tracketa - eta) + (trackphi - phi) * (trackphi - phi));
+ if (d <= fConeRadius) {
+ pt += track->Pt();
+ TLorentzVector lvtr(track->Px(),track->Py(),track->Pz(),track->E());
+ lvRC += lvtr;
+ }
+
+ track = tracks->GetNextAcceptParticle();
+ }
+ }
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEmcalJetMassBkg::GetPerpCone(TLorentzVector& lvRC,Float_t &pt, Float_t &eta, Float_t &phi, AliParticleContainer* tracks, AliClusterContainer* clusters, AliEmcalJet *jet) const
+{
+ // Get rigid cone.
+ lvRC.SetPxPyPzE(0.,0.,0.,0.);
+
+ eta = -999;
+ phi = -999;
+ pt = 0;
+
+ if (!tracks && !clusters)
+ return;
+
+ if(!jet)
+ return;
+
+ Float_t LJeta = jet->Eta();
+ Float_t LJphi = jet->Phi();
+
+ eta = LJeta;
+ phi = LJphi + 0.5*TMath::Pi();
+ if(phi>TMath::TwoPi()) phi-=TMath::TwoPi();
+
+ GetCone(lvRC,pt,eta,phi,tracks,clusters);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEmcalJetMassBkg::SetConeEtaPhiEMCAL()
+{
+ // Set default cuts for full cones
+
+ SetConeEtaLimits(-0.7+fConeRadius,0.7-fConeRadius);
+ SetConePhiLimits(1.405+fConeRadius,3.135-fConeRadius);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEmcalJetMassBkg::SetConeEtaPhiTPC()
+{
+ // Set default cuts for charged cones
+
+ SetConeEtaLimits(-0.9+fConeRadius, 0.9-fConeRadius);
+ SetConePhiLimits(-10, 10);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEmcalJetMassBkg::ExecOnce() {
+
+ AliAnalysisTaskEmcalJet::ExecOnce();
+
+ if (fTracksCont && fTracksCont->GetArray() == 0) fTracksCont = 0;
+ if (fCaloClustersCont && fCaloClustersCont->GetArray() == 0) fCaloClustersCont = 0;
+
+ if (fRCperEvent < 0) {
+ Double_t area = (fConeMaxEta - fConeMinEta) * (fConeMaxPhi - fConeMinPhi);
+ Double_t rcArea = TMath::Pi() * fConeRadius * fConeRadius;
+ fRCperEvent = TMath::FloorNint(area / rcArea - 0.5);
+ if (fRCperEvent == 0)
+ fRCperEvent = 1;
+ }
+
+ if (fMinRC2LJ < 0)
+ fMinRC2LJ = fConeRadius * 1.5;
+
+ const Float_t maxDist = TMath::Max(fConeMaxPhi - fConeMinPhi, fConeMaxEta - fConeMinEta) / 2;
+ if (fMinRC2LJ > maxDist) {
+ AliWarning(Form("The parameter fMinRC2LJ = %f is too large for the considered acceptance. "
+ "Will use fMinRC2LJ = %f", fMinRC2LJ, maxDist));
+ fMinRC2LJ = maxDist;
+ }
+
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskEmcalJetMassBkg::RetrieveEventObjects() {
+ //
+ // retrieve event objects
+ //
+
+ if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
+ return kFALSE;
+
+ return kTRUE;
+
+}
+
+//_______________________________________________________________________
+void AliAnalysisTaskEmcalJetMassBkg::Terminate(Option_t *)
+{
+ // Called once at the end of the analysis.
+}
+
--- /dev/null
+#ifndef ALIANALYSISTASKEMCALJETMASSBKG_H
+#define ALIANALYSISTASKEMCALJETMASSBKG_H
+
+class TH1;
+class TH2;
+class TH3;
+class TH3F;
+class THnSparse;
+class TClonesArray;
+class TArrayI;
+class AliAnalysisManager;
+class AliJetContainer;
+
+#include "AliAnalysisTaskEmcalJet.h"
+
+class AliAnalysisTaskEmcalJetMassBkg : public AliAnalysisTaskEmcalJet {
+ public:
+
+ AliAnalysisTaskEmcalJetMassBkg();
+ AliAnalysisTaskEmcalJetMassBkg(const char *name);
+ virtual ~AliAnalysisTaskEmcalJetMassBkg();
+
+ void UserCreateOutputObjects();
+ void Terminate(Option_t *option);
+
+ //Setters
+ void SetJetContainerBase(Int_t c) { fContainerBase = c ; }
+
+ void SetJetMinRC2LJ(Float_t d) { fMinRC2LJ = d ; }
+ void SetRCperEvent(Int_t n) { fRCperEvent = n ; }
+ void SetConeRadius(Double_t r) { fConeRadius = r ; }
+ void SetConeEtaPhiEMCAL() ;
+ void SetConeEtaPhiTPC() ;
+ void SetConeEtaLimits(Float_t min, Float_t max) { fConeMinEta = min, fConeMaxEta = max ; }
+ void SetConePhiLimits(Float_t min, Float_t max) { fConeMinPhi = min, fConeMaxPhi = max ; }
+
+ protected:
+ void ExecOnce();
+ Bool_t RetrieveEventObjects();
+ Bool_t Run();
+ Bool_t FillHistograms();
+
+ void GetRandomCone(TLorentzVector& lvRC, Float_t &pt, Float_t &eta, Float_t &phi, AliParticleContainer* tracks, AliClusterContainer* clusters, AliEmcalJet *jet = 0) const;
+ void GetCone(TLorentzVector& lvRC,Float_t &pt, Float_t eta, Float_t phi, AliParticleContainer* tracks, AliClusterContainer* clusters) const;
+ void GetPerpCone(TLorentzVector& lvRC, Float_t &pt, Float_t &eta, Float_t &phi, AliParticleContainer* tracks, AliClusterContainer* clusters, AliEmcalJet *jet = 0) const;
+
+ Int_t fContainerBase; // jets to be tagged
+
+ Float_t fMinRC2LJ; // Minimum distance random cone to leading jet
+ Int_t fRCperEvent; // No. of random cones per event
+ Double_t fConeRadius; // Radius of the random cones
+ Float_t fConeMinEta; // Minimum eta of the random cones
+ Float_t fConeMaxEta; // Maximum eta of the random cones
+ Float_t fConeMinPhi; // Minimum phi of the random cones
+ Float_t fConeMaxPhi; // Maximum phi of the random cones
+
+ AliJetContainer *fJetsCont; //!Jets
+ AliParticleContainer *fTracksCont; //!Tracks
+ AliClusterContainer *fCaloClustersCont; //!Clusters
+
+ private:
+ TH2F **fh2PtVsMassRC; //!pT vs mass of RC
+ TProfile **fpPtVsMassRC; //!pT vs Avg mass of RC
+ TH3F **fh2PtVsMassRCExLJDPhi; //!pT vs mass of RC
+ TProfile **fpPtVsMassRCExLJ; //!pT vs Avg mass of RC excluding area around leading jet
+ TH2F **fh2PtVsMassPerpConeLJ; //!pT vs mass of cone perpendicular to leading jet
+ TProfile **fpPtVsMassPerpConeLJ; //!pT vs Avg mass of cone perpendicular to leading jet
+ TH2F **fh2PtVsMassPerpConeTJ; //!pT vs mass of cone perpendicular to all tagged jets
+ TProfile **fpPtVsMassPerpConeTJ; //!pT vs Avg mass of cone perpendicular to all tagged jets
+
+ TH2F *fh2CentVsMassRC; //!cent vs mass of RC
+ TH2F *fh2CentVsMassRCExLJ; //!cent vs mass of RC excluding area around leading jet
+ TH2F *fh2CentVsMassPerpConeLJ; //!cent vs mass of RC perpendicular to leading jet
+ TH2F *fh2CentVsMassPerpConeTJ; //!cent vs mass of RC perpendicular to tagged jets
+
+ TH2F *fh2MultVsMassRC; //!track multiplicity vs mass of RC
+ TH2F *fh2MultVsMassRCExLJ; //!track multiplicity vs mass of RC excluding area around leading jet
+ TH2F *fh2MultVsMassPerpConeLJ; //!track multiplicity vs mass of RC perpendicular to leading jet
+ TH2F *fh2MultVsMassPerpConeTJ; //!track multiplicity vs mass of RC perpendicular to tagged jets
+
+ AliAnalysisTaskEmcalJetMassBkg(const AliAnalysisTaskEmcalJetMassBkg&); // not implemented
+ AliAnalysisTaskEmcalJetMassBkg &operator=(const AliAnalysisTaskEmcalJetMassBkg&); // not implemented
+
+ ClassDef(AliAnalysisTaskEmcalJetMassBkg, 1)
+};
+#endif
+
-// $Id: AliAnalysisTaskEmcalJetSpectraMECpA.cxx 3010 2012-06-10 05:40:56Z loizides $
+// $Id$
//
// Jet spectrum task.
//
#ifndef AliAnalysisTaskEmcalJetSpectraMECpA_h
#define AliAnalysisTaskEmcalJetSpectraMECpA_h
-// $Id: AliAnalysisTaskEmcalJetSpectraMECpA.h 3010 2012-06-10 05:40:56Z loizides $
+// $Id$
class TH1F;
--- /dev/null
+//
+// Jet tagger analysis task.
+//
+// Author: M.Verweij
+
+#include <TClonesArray.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <THnSparse.h>
+#include <TList.h>
+#include <TLorentzVector.h>
+#include <TProfile.h>
+#include <TChain.h>
+#include <TSystem.h>
+#include <TFile.h>
+#include <TKey.h>
+
+#include "AliVCluster.h"
+#include "AliVTrack.h"
+#include "AliEmcalJet.h"
+#include "AliRhoParameter.h"
+#include "AliLog.h"
+#include "AliEmcalParticle.h"
+#include "AliMCEvent.h"
+#include "AliGenPythiaEventHeader.h"
+#include "AliAODMCHeader.h"
+#include "AliMCEvent.h"
+#include "AliAnalysisManager.h"
+#include "AliJetContainer.h"
+
+#include "AliAnalysisTaskEmcalJetTagger.h"
+
+ClassImp(AliAnalysisTaskEmcalJetTagger)
+
+//________________________________________________________________________
+AliAnalysisTaskEmcalJetTagger::AliAnalysisTaskEmcalJetTagger() :
+ AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalJetTagger", kTRUE),
+ fJetTaggingType(kTag),
+ fJetTaggingMethod(kGeo),
+ fContainerBase(0),
+ fContainerTag(1),
+ fMatchingDone(0),
+ fh3PtJet1VsDeltaEtaDeltaPhi(0),
+ fh2PtJet1VsDeltaR(0),
+ fh2PtJet1VsLeadPtAllSel(0),
+ fh2PtJet1VsLeadPtTagged(0),
+ fh2PtJet1VsPtJet2(0)
+{
+ // Default constructor.
+
+ fh3PtJet1VsDeltaEtaDeltaPhi = new TH3F*[fNcentBins];
+ fh2PtJet1VsDeltaR = new TH2F*[fNcentBins];
+ fh2PtJet1VsLeadPtAllSel = new TH2F*[fNcentBins];
+ fh2PtJet1VsLeadPtTagged = new TH2F*[fNcentBins];
+ fh2PtJet1VsPtJet2 = new TH2F*[fNcentBins];
+
+ for (Int_t i = 0; i < fNcentBins; i++) {
+ fh3PtJet1VsDeltaEtaDeltaPhi[i] = 0;
+ fh2PtJet1VsDeltaR[i] = 0;
+ fh2PtJet1VsLeadPtAllSel[i] = 0;
+ fh2PtJet1VsLeadPtTagged[i] = 0;
+ fh2PtJet1VsPtJet2[i] = 0;
+ }
+
+ SetMakeGeneralHistograms(kTRUE);
+
+}
+
+//________________________________________________________________________
+AliAnalysisTaskEmcalJetTagger::AliAnalysisTaskEmcalJetTagger(const char *name) :
+ AliAnalysisTaskEmcalJet(name, kTRUE),
+ fJetTaggingType(kTag),
+ fJetTaggingMethod(kGeo),
+ fContainerBase(0),
+ fContainerTag(1),
+ fMatchingDone(0),
+ fh3PtJet1VsDeltaEtaDeltaPhi(0),
+ fh2PtJet1VsDeltaR(0),
+ fh2PtJet1VsLeadPtAllSel(0),
+ fh2PtJet1VsLeadPtTagged(0),
+ fh2PtJet1VsPtJet2(0)
+{
+ // Standard constructor.
+
+ fh3PtJet1VsDeltaEtaDeltaPhi = new TH3F*[fNcentBins];
+ fh2PtJet1VsDeltaR = new TH2F*[fNcentBins];
+ fh2PtJet1VsLeadPtAllSel = new TH2F*[fNcentBins];
+ fh2PtJet1VsLeadPtTagged = new TH2F*[fNcentBins];
+ fh2PtJet1VsPtJet2 = new TH2F*[fNcentBins];
+
+ for (Int_t i = 0; i < fNcentBins; i++) {
+ fh3PtJet1VsDeltaEtaDeltaPhi[i] = 0;
+ fh2PtJet1VsDeltaR[i] = 0;
+ fh2PtJet1VsLeadPtAllSel[i] = 0;
+ fh2PtJet1VsLeadPtTagged[i] = 0;
+ fh2PtJet1VsPtJet2[i] = 0;
+ }
+
+ SetMakeGeneralHistograms(kTRUE);
+}
+
+//________________________________________________________________________
+AliAnalysisTaskEmcalJetTagger::~AliAnalysisTaskEmcalJetTagger()
+{
+ // Destructor.
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEmcalJetTagger::UserCreateOutputObjects()
+{
+ // Create user output.
+
+ AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
+
+ Bool_t oldStatus = TH1::AddDirectoryStatus();
+ TH1::AddDirectory(kFALSE);
+
+ const Int_t nBinsPt = 250;
+ const Int_t nBinsDPhi = 100;
+ const Int_t nBinsDEta = 100;
+ const Int_t nBinsDR = 50;
+
+ const Double_t minPt = -50.;
+ const Double_t maxPt = 200.;
+ const Double_t minDPhi = -0.5;
+ const Double_t maxDPhi = 0.5;
+ const Double_t minDEta = -0.5;
+ const Double_t maxDEta = 0.5;
+ const Double_t minDR = 0.;
+ const Double_t maxDR = 0.5;
+
+ TString histName = "";
+ TString histTitle = "";
+ for (Int_t i = 0; i < fNcentBins; i++) {
+
+ histName = TString::Format("fh3PtJet1VsDeltaEtaDeltaPhi_%d",i);
+ histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{#Delta#eta};#it{#Delta#varphi}",histName.Data());
+ fh3PtJet1VsDeltaEtaDeltaPhi[i] = new TH3F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDEta,minDEta,maxDEta,nBinsDPhi,minDPhi,maxDPhi);
+ fOutput->Add(fh3PtJet1VsDeltaEtaDeltaPhi[i]);
+
+ histName = TString::Format("fh2PtJet1VsDeltaR_%d",i);
+ histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{#Delta R}",histName.Data());
+ fh2PtJet1VsDeltaR[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsDR,minDR,maxDR);
+ fOutput->Add(fh2PtJet1VsDeltaR[i]);
+
+ histName = TString::Format("fh2PtJet1VsLeadPtAllSel_%d",i);
+ histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
+ fh2PtJet1VsLeadPtAllSel[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
+ fOutput->Add(fh2PtJet1VsLeadPtAllSel[i]);
+
+ histName = TString::Format("fh2PtJet1VsLeadPtTagged_%d",i);
+ histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,lead trk}",histName.Data());
+ fh2PtJet1VsLeadPtTagged[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,20,0.,20.);
+ fOutput->Add(fh2PtJet1VsLeadPtTagged[i]);
+
+ histName = TString::Format("fh2PtJet1VsPtJet2_%d",i);
+ histTitle = TString::Format("%s;#it{p}_{T,jet1};#it{p}_{T,jet2}",histName.Data());
+ fh2PtJet1VsPtJet2[i] = new TH2F(histName.Data(),histTitle.Data(),nBinsPt,minPt,maxPt,nBinsPt,minPt,maxPt);
+ fOutput->Add(fh2PtJet1VsPtJet2[i]);
+
+ }
+
+ // =========== Switch on Sumw2 for all histos ===========
+ for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
+ TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
+ if (h1){
+ h1->Sumw2();
+ continue;
+ }
+ THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
+ if(hn)hn->Sumw2();
+ }
+
+ TH1::AddDirectory(oldStatus);
+
+ PostData(1, fOutput); // Post data for ALL output slots > 0 here.
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskEmcalJetTagger::Run()
+{
+ // Run analysis code here, if needed. It will be executed before FillHistograms().
+
+ if(fJetTaggingMethod==kGeo)
+ MatchJetsGeo(fContainerBase,fContainerTag,0,0.3,2);
+
+ return kTRUE;
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskEmcalJetTagger::FillHistograms()
+{
+ // Fill histograms.
+
+ for(int i = 0; i < GetNJets(fContainerBase);++i) {
+ AliEmcalJet *jet1 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(i, fContainerBase));
+ if(!jet1) continue;
+
+ Double_t ptJet1 = jet1->Pt() - GetRhoVal(fContainerBase)*jet1->Area();
+ fh2PtJet1VsLeadPtAllSel[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
+
+ if(jet1->GetTagStatus()<1 && fJetTaggingType==kTag)
+ continue;
+
+ AliEmcalJet *jet2 = jet1->GetTaggedJet();
+ if(!jet2) continue;
+ Double_t ptJet2 = jet2->Pt() - GetRhoVal(fContainerTag)*jet2->Area();
+ fh2PtJet1VsLeadPtTagged[fCentBin]->Fill(ptJet1,jet1->MaxTrackPt());
+
+ fh2PtJet1VsPtJet2[fCentBin]->Fill(ptJet1,ptJet2);
+
+ Double_t dPhi = GetDeltaPhi(jet1->Phi(),jet2->Phi());
+ if(dPhi>TMath::Pi())
+ dPhi -= TMath::TwoPi();
+ if(dPhi<(-1.*TMath::Pi()))
+ dPhi += TMath::TwoPi();
+
+ fh3PtJet1VsDeltaEtaDeltaPhi[fCentBin]->Fill(ptJet1,jet1->Eta()-jet2->Eta(),dPhi);
+ fh2PtJet1VsDeltaR[fCentBin]->Fill(ptJet1,GetDeltaR(jet1,jet2));
+ }
+
+ return kTRUE;
+
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEmcalJetTagger::ResetTagging(const Int_t c) {
+
+ //Reset tagging of container c
+
+ for(int i = 0;i<GetNJets(c);i++){
+ AliEmcalJet *jet = static_cast<AliEmcalJet*>(GetJetFromArray(i, c));
+ if(fJetTaggingType==kClosest)
+ jet->ResetMatching();
+ else if(fJetTaggingType==kTag) {
+ jet->SetTaggedJet(0x0);
+ jet->SetTagStatus(-1);
+ }
+ }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskEmcalJetTagger::MatchJetsGeo(Int_t c1, Int_t c2,
+ Int_t iDebug, Float_t maxDist, Int_t type) {
+
+ //
+ // Match the full jets to the corresponding charged jets
+ // Translation of AliAnalysisHelperJetTasks::GetClosestJets to AliEmcalJet objects
+ // type: 0 = use acceptance cuts of container 1 = allow 0.1 one more for c2 in eta 2 = allow 0.1 more in eta and phi for c2
+
+ if(c1<0) c1 = fContainerBase;
+ if(c2<0) c2 = fContainerTag;
+
+ const Int_t nJets1 = GetNJets(c1);
+ const Int_t nJets2 = GetNJets(c2);
+
+ if(nJets1==0 || nJets2==0) return;
+
+ ResetTagging(c1);
+ ResetTagging(c2);
+
+ fMatchingDone = kFALSE;
+
+ TArrayI faMatchIndex1;
+ faMatchIndex1.Set(nJets2+1);
+ faMatchIndex1.Reset(-1);
+
+ TArrayI faMatchIndex2;
+ faMatchIndex2.Set(nJets1+1);
+ faMatchIndex2.Reset(-1);
+
+ static TArrayS iFlag(nJets1*nJets2);
+ if(iFlag.GetSize()<(nJets1*nJets2)){
+ iFlag.Set(nJets1*nJets1+1);
+ }
+ iFlag.Reset(0);
+
+ AliJetContainer *cont2 = GetJetContainer(c2);
+
+ // find the closest distance to the full jet
+ for(int i = 0;i<nJets1;i++){
+
+ AliEmcalJet *jet1 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(i, c1));
+ if(!jet1) continue;
+
+ Float_t dist = maxDist;
+
+ for(int j = 0;j <nJets2; j++){
+ AliEmcalJet *jet2 = 0x0;
+ if(type==0)
+ jet2 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(j, c2));
+ else {
+ jet2 = static_cast<AliEmcalJet*>(GetJetFromArray(j, c2));
+ if(!jet2) continue;
+ if(type>0) {
+ if(jet2->Eta()<(cont2->GetJetEtaMin()-0.1) || jet2->Eta()>(cont2->GetJetEtaMax()+0.1))
+ continue;
+ if(type==2) {
+ if(jet2->Phi()<(cont2->GetJetPhiMin()-0.1) || jet2->Phi()>(cont2->GetJetPhiMax()+0.1))
+ continue;
+ }
+ }
+ }
+ if(!jet2)
+ continue;
+
+ Double_t dR = GetDeltaR(jet1,jet2);
+ if(dR<dist && dR<maxDist){
+ faMatchIndex2[i]=j;
+ dist = dR;
+ }
+ }//j jet loop
+ if(faMatchIndex2[i]>=0) iFlag[i*nJets2+faMatchIndex2[i]]+=1;//j closest to i
+ if(iDebug>10) Printf("Full Distance (%d)--(%d) %3.3f flag[%d] = %d",i,faMatchIndex2[i],dist,i*nJets2+faMatchIndex2[i],iFlag[i*nJets2+faMatchIndex2[i]]);
+
+ }//i jet loop
+
+
+ // other way around
+ for(int j = 0;j<nJets2;j++){
+ AliEmcalJet *jet2 = 0x0;
+ if(type==0)
+ jet2 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(j, c2));
+ else {
+ jet2 = static_cast<AliEmcalJet*>(GetJetFromArray(j, c2));
+ if(!jet2) continue;;
+ if(type>0) {
+ if(jet2->Eta()<(cont2->GetJetEtaMin()-0.1) || jet2->Eta()>(cont2->GetJetEtaMax()+0.1))
+ continue;
+ if(type==2) {
+ if(jet2->Phi()<(cont2->GetJetPhiMin()-0.1) || jet2->Phi()>(cont2->GetJetPhiMax()+0.1))
+ continue;
+ }
+ }
+ }
+ if(!jet2)
+ continue;
+
+ Float_t dist = maxDist;
+ for(int i = 0;i<nJets1;i++){
+ AliEmcalJet *jet1 = static_cast<AliEmcalJet*>(GetAcceptJetFromArray(i, c1));
+ if(!jet1)
+ continue;
+
+ Double_t dR = GetDeltaR(jet1,jet2);
+ if(dR<dist && dR<maxDist){
+ faMatchIndex1[j]=i;
+ dist = dR;
+ }
+ }
+ if(faMatchIndex1[j]>=0) iFlag[faMatchIndex1[j]*nJets2+j]+=2;//i closest to j
+ if(iDebug>10) Printf("Other way Distance (%d)--(%d) %3.3f flag[%d] = %d",faMatchIndex1[j],j,dist,faMatchIndex1[j]*nJets2+j,iFlag[faMatchIndex1[j]*nJets2+j]);
+
+ }
+
+ // check for "true" correlations
+ for(int i = 0;i<nJets1;i++){
+ AliEmcalJet *jet1 = static_cast<AliEmcalJet*>(GetJetFromArray(i, c1));
+ for(int j = 0;j<nJets2;j++){
+ AliEmcalJet *jet2 = static_cast<AliEmcalJet*>(GetJetFromArray(j, c2));
+ if(iDebug>10) AliInfo(Form("%s: Flag[%d][%d] %d ",GetName(),i,j,iFlag[i*nJets2+j]));
+
+ // we have a uniqe correlation
+ if(iFlag[i*nJets2+j]==3){
+
+ Double_t dR = GetDeltaR(jet1,jet2);
+ if(iDebug>1) Printf("closest jets %d %d dR = %f",j,i,dR);
+
+ if(fJetTaggingType==kTag) {
+ jet1->SetTaggedJet(jet2);
+ jet1->SetTagStatus(1);
+
+ jet2->SetTaggedJet(jet1);
+ jet2->SetTagStatus(1);
+ }
+ else if(fJetTaggingType==kClosest) {
+ jet1->SetClosestJet(jet2,dR);
+ jet2->SetClosestJet(jet1,dR);
+ }
+ }
+ }
+ }
+
+ fMatchingDone = kTRUE;
+
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskEmcalJetTagger::GetDeltaR(const AliEmcalJet* jet1, const AliEmcalJet* jet2) const {
+ //
+ // Helper function to calculate the distance between two jets
+ //
+
+ Double_t dPhi = jet1->Phi() - jet2->Phi();
+ if(dPhi>TMath::Pi())
+ dPhi -= TMath::TwoPi();
+ if(dPhi<(-1.*TMath::Pi()))
+ dPhi += TMath::TwoPi();
+ Double_t dEta = jet1->Eta() - jet2->Eta();
+ Double_t dR = TMath::Sqrt(dPhi*dPhi+dEta*dEta);
+ return dR;
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskEmcalJetTagger::GetDeltaPhi(const AliEmcalJet* jet1, const AliEmcalJet* jet2) {
+ //
+ // Calculate azimuthal angle between the axises of the jets
+ //
+
+ return GetDeltaPhi(jet1->Phi(),jet2->Phi());
+
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskEmcalJetTagger::GetDeltaPhi(Double_t phi1,Double_t phi2) {
+ //
+ // Calculate azimuthal angle between the axises of the jets
+ //
+
+ Double_t dPhi = phi1-phi2;
+
+ if(dPhi <-0.5*TMath::Pi()) dPhi += TMath::TwoPi();
+ if(dPhi > 1.5*TMath::Pi()) dPhi -= TMath::TwoPi();
+
+ return dPhi;
+}
+
+//________________________________________________________________________
+Double_t AliAnalysisTaskEmcalJetTagger::GetFractionSharedPt(const AliEmcalJet *jet1, const AliEmcalJet *jet2) const
+{
+ //
+ // Get fraction of shared pT between matched full and charged jet
+ // Uses charged jet pT as baseline: fraction = \Sum_{const,full jet} pT,const,i / pT,jet,ch
+ //
+
+ Double_t fraction = 0.;
+ Double_t jetPt2 = jet2->Pt();
+
+ if(jetPt2>0) {
+
+ AliDebug(2,Form("%s: nConstituents: %d, ch: %d chne: %d ne: %d",GetName(),jet1->GetNumberOfConstituents(),jet2->GetNumberOfTracks(),jet1->GetNumberOfTracks(),jet1->GetNumberOfClusters()));
+
+ Double_t sumPt = 0.;
+ AliVParticle *vpf = 0x0;
+ Int_t iFound = 0;
+ for(Int_t icc=0; icc<jet2->GetNumberOfTracks(); icc++) {
+ Int_t idx = (Int_t)jet2->TrackAt(icc);
+ iFound = 0;
+ for(Int_t icf=0; icf<jet1->GetNumberOfTracks(); icf++) {
+ if(idx == jet1->TrackAt(icf) && iFound==0 ) {
+ iFound=1;
+ vpf = static_cast<AliVParticle*>(jet1->TrackAt(icf, fTracks));
+ if(!vpf) continue;
+ if(vpf->Charge()!=0)
+ sumPt += vpf->Pt();
+ continue;
+ }
+ }
+ }
+
+ fraction = sumPt/jetPt2;
+ }
+
+ return fraction;
+
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskEmcalJetTagger::RetrieveEventObjects() {
+ //
+ // retrieve event objects
+ //
+
+ if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
+ return kFALSE;
+
+ return kTRUE;
+
+}
+
+//_______________________________________________________________________
+void AliAnalysisTaskEmcalJetTagger::Terminate(Option_t *)
+{
+ // Called once at the end of the analysis.
+}
+
--- /dev/null
+#ifndef ALIANALYSISTASKEMCALJETTAGGER_H
+#define ALIANALYSISTASKEMCALJETTAGGER_H
+
+class TH1;
+class TH2;
+class TH3;
+class TH3F;
+class THnSparse;
+class TClonesArray;
+class TArrayI;
+class AliAnalysisManager;
+class AliJetContainer;
+
+#include "AliAnalysisTaskEmcalJet.h"
+
+class AliAnalysisTaskEmcalJetTagger : public AliAnalysisTaskEmcalJet {
+ public:
+ enum JetTaggingMethod {
+ kGeo = 0,
+ kFraction = 1
+ };
+
+ enum JetTaggingType {
+ kTag = 0,
+ kClosest = 1
+ };
+
+ AliAnalysisTaskEmcalJetTagger();
+ AliAnalysisTaskEmcalJetTagger(const char *name);
+ virtual ~AliAnalysisTaskEmcalJetTagger();
+
+ void UserCreateOutputObjects();
+ void Terminate(Option_t *option);
+
+ //Setters
+ void SetJetContainerBase(Int_t c) { fContainerBase = c;}
+ void SetJetContainerTag(Int_t c) { fContainerTag = c;}
+
+ void SetJetTaggingType(JetTaggingType t) { fJetTaggingType = t;}
+ void SetJetTaggingMethod(JetTaggingMethod m) { fJetTaggingMethod = m;}
+
+ protected:
+ Bool_t RetrieveEventObjects();
+ Bool_t Run();
+ Bool_t FillHistograms();
+
+
+ Double_t GetDeltaPhi(const AliEmcalJet* jet1, const AliEmcalJet* jet2);
+ Double_t GetDeltaPhi(Double_t phi1,Double_t phi2);
+ Double_t GetDeltaR(const AliEmcalJet* jet1, const AliEmcalJet* jet2) const;
+ Double_t GetFractionSharedPt(const AliEmcalJet *jet1, const AliEmcalJet *jet2) const;
+
+ void MatchJetsGeo(Int_t c1 = -1, Int_t c2 = -1, Int_t iDebug = 0, Float_t maxDist = 0.3, Int_t type = 2);
+ void ResetTagging(const Int_t c);
+
+ JetTaggingType fJetTaggingType; // jet matching type
+ JetTaggingMethod fJetTaggingMethod; // jet matching method
+
+ Int_t fContainerBase; // jets to be tagged
+ Int_t fContainerTag; // jets used for tagging
+
+ private:
+ Bool_t fMatchingDone; // flag to indicate if matching is done or not
+ TH3F **fh3PtJet1VsDeltaEtaDeltaPhi; //!pt jet 1 vs deta vs dphi
+ TH2F **fh2PtJet1VsDeltaR; //!pt jet 1 vs dR
+
+ TH2F **fh2PtJet1VsLeadPtAllSel; //!all jets after std selection
+ TH2F **fh2PtJet1VsLeadPtTagged; //!tagged jets
+ TH2F **fh2PtJet1VsPtJet2; //!pT of base jet vs tagged jet
+
+ AliAnalysisTaskEmcalJetTagger(const AliAnalysisTaskEmcalJetTagger&); // not implemented
+ AliAnalysisTaskEmcalJetTagger &operator=(const AliAnalysisTaskEmcalJetTagger&); // not implemented
+
+ ClassDef(AliAnalysisTaskEmcalJetTagger, 1)
+};
+#endif
+
if (!jet)
continue; //jet not selected
- Double_t jetPt = jet->Pt() - GetRhoVal(fContainerCharged)*jet->Area();;
+ Double_t jetPt = jet->Pt() - GetRhoVal(fContainerCharged)*jet->Area();
if(jetPt>ptLeadJet2) ptLeadJet2=jetPt;
fh3PtEtaPhiJetCharged->Fill(jetPt,jet->Eta(),jet->Phi());
fh3PtEtaAreaJetCharged->Fill(jetPt,jet->Eta(),jet->Area());
const Double_t minClPt = 0.30,
const Double_t ghostArea = 0.005,
const Double_t radius = 0.4,
+ const Int_t recombScheme = 1,
const char *tag = "Jet"
+
)
{
// Get the pointer to the existing analysis manager via the static access method.
else if (minClPt>=1.0)
sprintf(ETString,"ET%4.0f",minClPt*1000.0);
+ char recombSchemeString[200];
+ if(recombScheme==0)
+ sprintf(recombSchemeString,"%s","E_scheme");
+ else if(recombScheme==1)
+ sprintf(recombSchemeString,"%s","pt_scheme");
+ else if(recombScheme==2)
+ sprintf(recombSchemeString,"%s","pt2_scheme");
+ else if(recombScheme==3)
+ sprintf(recombSchemeString,"%s","Et_scheme");
+ else if(recombScheme==4)
+ sprintf(recombSchemeString,"%s","Et2_scheme");
+ else if(recombScheme==5)
+ sprintf(recombSchemeString,"%s","BIpt_scheme");
+ else if(recombScheme==6)
+ sprintf(recombSchemeString,"%s","BIpt2_scheme");
+ else if(recombScheme==99)
+ sprintf(recombSchemeString,"%s","ext_scheme");
+ else {
+ ::Error("AddTaskAliEmcalJet", "Recombination scheme not recognized.");
+ return NULL;
+ }
+
TString name;
if (*nTracks && *nClusters)
- name = TString(Form("%s_%s%s%s_%s_%s_%s_%s",
- tag,algoString,typeString,radiusString,nTracks,pTString,nClusters,ETString));
+ name = TString(Form("%s_%s%s%s_%s_%s_%s_%s_%s",
+ tag,algoString,typeString,radiusString,nTracks,pTString,nClusters,ETString,recombSchemeString));
else if (!*nClusters)
- name = TString(Form("%s_%s%s%s_%s_%s",
- tag,algoString,typeString,radiusString,nTracks,pTString));
+ name = TString(Form("%s_%s%s%s_%s_%s_%s",
+ tag,algoString,typeString,radiusString,nTracks,pTString,recombSchemeString));
else if (!*nTracks)
- name = TString(Form("%s_%s%s%s_%s_%s",
- tag,algoString,typeString,radiusString,nClusters,ETString));
+ name = TString(Form("%s_%s%s%s_%s_%s_%s",
+ tag,algoString,typeString,radiusString,nClusters,ETString,recombSchemeString));
AliEmcalJetTask* mgrTask = mgr->GetTask(name.Data());
if (mgrTask)
if ((type & (AliEmcalJetTask::kRX1Jet|AliEmcalJetTask::kRX2Jet|AliEmcalJetTask::kRX3Jet)) != 0)
jetTask->SetRadius(radius);
jetTask->SetGhostArea(ghostArea);
+ jetTask->SetRecombScheme(recombScheme);
//-------------------------------------------------------
// Final settings, pass to manager and set the containers
const Int_t type = 0,
const Double_t minTrPt = 0.15,
const Double_t minClPt = 0.30,
- const Double_t ghostArea = 0.01 ,
+ const Double_t ghostArea = 0.01,
+ const Int_t recombScheme = 1,
const char *tag = "Jet"
)
{
else
jetType |= AliEmcalJetTask::kRX1Jet;
- return AddTaskEmcalJet(jetType, nTracks, nClusters, minTrPt, minClPt, ghostArea, radius, tag);
+ return AddTaskEmcalJet(jetType, nTracks, nClusters, minTrPt, minClPt, ghostArea, radius, recombScheme, tag);
}
--- /dev/null
+enum AlgoType {kKT, kANTIKT};
+enum JetType {kFULLJETS, kCHARGEDJETS, kNEUTRALJETS};
+
+AliAnalysisTaskEmcalJetMass* AddTaskEmcalJetMass(const char * njetsBase,
+ const char * njetsTag,
+ const Double_t R,
+ const char * nrhoBase,
+ const char * nrhoTag,
+ const char * ntracks,
+ const char * nclusters,
+ const char *type,
+ const char *CentEst,
+ Int_t pSel,
+ TString trigClass = "",
+ TString kEmcalTriggers = "");
+
+AliAnalysisTaskEmcalJetMass* AddTaskEmcalJetMass(TString kTracksName = "PicoTracks",
+ TString kClusName = "caloClustersCorr",
+ Double_t R = 0.4,
+ Double_t ptminTrack = 0.15,
+ Double_t etminClus = 0.3,
+ Double_t ptminTag = 4.,
+ Int_t rhoType = 1,
+ const char *type = "EMCAL",
+ const char *CentEst = "V0M",
+ Int_t pSel = AliVEvent::kCentral | AliVEvent::kSemiCentral | AliVEvent::kMB,
+ TString trigClass = "",
+ TString kEmcalTriggers = "",
+ TString kPeriod = "LHC11h"
+ ) {
+
+ AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+ if (!mgr)
+ {
+ Error("AddTaskEmcalJetMass","No analysis manager found.");
+ return 0;
+ }
+ Bool_t ismc=kFALSE;
+ ismc = (mgr->GetMCtruthEventHandler())?kTRUE:kFALSE;
+
+ // Check the analysis type using the event handlers connected to the analysis manager.
+ //==============================================================================
+ if (!mgr->GetInputEventHandler())
+ {
+ ::Error("AddTaskEmcalJetMass", "This task requires an input event handler");
+ return NULL;
+ }
+
+ // #### Add necessary jet finder tasks
+ gROOT->LoadMacro("$ALICE_ROOT/PWGJE/EMCALJetTasks/macros/AddTaskEmcalJet.C");
+
+ AliEmcalJetTask* jetFinderTaskBase = 0x0;
+ if (strcmp(type,"TPC")==0)
+ jetFinderTaskBase = AddTaskEmcalJet(kTracksName, "", kANTIKT, R, kCHARGEDJETS, ptminTrack, etminClus,0.005,0);
+ else if (strcmp(type,"EMCAL")==0)
+ jetFinderTaskBase = AddTaskEmcalJet(kTracksName, kClusName, kANTIKT, R, kFULLJETS, ptminTrack, etminClus,0.005,0);
+
+ TString strJetsBase = jetFinderTaskBase->GetName();
+
+ AliAnalysisTaskRhoBase *rhoTaskBase;
+ TString rhoNameBase = "";
+ if(rhoType==1) {
+ rhoTaskBase = AttachRhoTask(kPeriod,kTracksName,kClusName,R,ptminTrack,etminClus);
+ rhoTaskBase->SetCentralityEstimator(CentEst);
+ rhoTaskBase->SelectCollisionCandidates(AliVEvent::kAny);
+
+ if (strcmp(type,"TPC")==0)
+ rhoNameBase = rhoTaskBase->GetOutRhoName();
+ if (strcmp(type,"EMCAL")==0)
+ rhoNameBase = rhoTaskBase->GetOutRhoScaledName();
+ }
+
+ //Configure jet tagger task
+ AliAnalysisTaskEmcalJetMass *task = AddTaskEmcalJetMass(jetFinderTaskBase->GetName(),
+ R,
+ rhoNameBase,
+ kTracksName.Data(),
+ kClusName.Data(),
+ type,
+ CentEst,
+ pSel,
+ trigClass,
+ kEmcalTriggers
+ );
+
+ return task;
+
+}
+
+AliAnalysisTaskEmcalJetMass* AddTaskEmcalJetMass(const char * njetsBase,
+ const Double_t R,
+ const char * nrhoBase,
+ const char * ntracks,
+ const char * nclusters,
+ const char * type,
+ const char * CentEst,
+ Int_t pSel,
+ TString trigClass,
+ TString kEmcalTriggers
+ ) {
+
+ AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+ if (!mgr)
+ {
+ Error("AddTaskEmcalJetMass","No analysis manager found.");
+ return 0;
+ }
+ Bool_t ismc=kFALSE;
+ ismc = (mgr->GetMCtruthEventHandler())?kTRUE:kFALSE;
+
+ // Check the analysis type using the event handlers connected to the analysis manager.
+ //==============================================================================
+ if (!mgr->GetInputEventHandler())
+ {
+ ::Error("AddTaskEmcalJetMass", "This task requires an input event handler");
+ return NULL;
+ }
+
+ TString wagonName = Form("JetMass_%s_TC%s",njetsBase,trigClass.Data());
+
+ //Configure jet tagger task
+ AliAnalysisTaskEmcalJetMass *task = new AliAnalysisTaskEmcalJetMass(wagonName.Data());
+
+ task->SetNCentBins(4);
+ //task->SetVzRange(-10.,10.);
+
+ AliParticleContainer *trackCont = task->AddParticleContainer(ntracks);
+ AliClusterContainer *clusterCont = task->AddClusterContainer(nclusters);
+
+ task->SetJetContainerBase(0);
+
+ TString strType(type);
+ AliJetContainer *jetContBase = task->AddJetContainer(njetsBase,strType,R);
+ if(jetContBase) {
+ jetContBase->SetRhoName(nrhoBase);
+ jetContBase->ConnectParticleContainer(trackCont);
+ jetContBase->ConnectClusterContainer(clusterCont);
+ jetContBase->SetZLeadingCut(0.98,0.98);
+ task->SetPercAreaCut(0.6, 0);
+ }
+
+ task->SetCaloTriggerPatchInfoName(kEmcalTriggers.Data());
+
+ task->SetCentralityEstimator(CentEst);
+
+ task->SelectCollisionCandidates(pSel);
+
+ task->SetUseAliAnaUtils(kFALSE);
+
+ mgr->AddTask(task);
+
+ //Connnect input
+ mgr->ConnectInput (task, 0, mgr->GetCommonInputContainer() );
+
+ //Connect output
+ TString contName(wagonName);
+ TString outputfile = Form("%s",AliAnalysisManager::GetCommonFileName());
+ AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contName.Data(), TList::Class(),AliAnalysisManager::kOutputContainer,outputfile);
+ mgr->ConnectOutput(task,1,coutput1);
+
+ return task;
+
+}
+
+
+AliAnalysisTaskRhoBase *AttachRhoTask(TString kPeriod ,// = "LHC13b",
+ TString kTracksName ,// = "PicoTracks",
+ TString kClusName ,// = "caloClustersCorr",
+ Double_t R ,// = 0.4,
+ Double_t ptminTrack ,// = 0.15,
+ Double_t etminClus // = 0.3
+ ) {
+
+ AliAnalysisTaskRhoBase *rhoTaskBase;
+
+ kPeriod.ToLower();
+
+ // Add kt jet finder and rho task in case we want background subtraction
+ AliEmcalJetTask *jetFinderKt;
+ AliEmcalJetTask *jetFinderAKt;
+ jetFinderKt = AddTaskEmcalJet(kTracksName, "", kKT, R, kCHARGEDJETS, ptminTrack, etminClus,0.005,0);
+ jetFinderAKt = AddTaskEmcalJet(kTracksName, "", kANTIKT, R, kCHARGEDJETS, ptminTrack, etminClus,0.005,0);
+
+ if(kPeriod.EqualTo("lhc13b") || kPeriod.EqualTo("lhc13c") || kPeriod.EqualTo("lhc13d") || kPeriod.EqualTo("lhc13e") || kPeriod.EqualTo("lhc13f")) {
+
+ jetFinderKt->SetMinJetPt(0.);
+
+ gROOT->LoadMacro("$ALICE_ROOT/PWGJE/EMCALJetTasks/macros/AddTaskRhoSparse.C");
+ TF1 *fScale = new TF1("fScale","1.42",0.,100.); //scale factor for pPb
+ AliAnalysisTaskRhoSparse *rhoTaskSparse = AddTaskRhoSparse(
+ jetFinderKt->GetName(),
+ jetFinderAKt->GetName(),
+ kTracksName,
+ kClusName,
+ Form("RhoSparseR%03d",(int)(100*R)),
+ R,
+ "TPC",
+ 0.01,
+ 0.15,
+ 0,
+ fScale,
+ 0,
+ kTRUE,
+ Form("RhoSparseR%03d",(int)(100*R)),
+ kTRUE
+ );
+ rhoTaskSparse->SetUseAliAnaUtils(kTRUE);
+
+ rhoTaskBase = dynamic_cast<AliAnalysisTaskRhoBase*>rhoTaskSparse;
+ }
+ else if(kPeriod.EqualTo("lhc10h") || kPeriod.EqualTo("lhc11h") ) {
+
+ gROOT->LoadMacro("$ALICE_ROOT/PWGJE/EMCALJetTasks/macros/AddTaskRho.C");
+
+ TF1* sfunc = new TF1("sfunc","[0]*x*x+[1]*x+[2]",-1,100);
+ sfunc->SetParameter(2,1.76458);
+ sfunc->SetParameter(1,-0.0111656);
+ sfunc->SetParameter(0,0.000107296);
+ AliAnalysisTaskRho *rhoTask = AddTaskRho(
+ jetFinderKt->GetName(),
+ kTracksName,
+ kClusName,
+ Form("RhoR%03dptmin%3.0f",(int)(100*R),ptminTrack*1000.0),
+ R,
+ "TPC",
+ 0.01,
+ 0,
+ sfunc,
+ 2,
+ kTRUE);
+ rhoTask->SetHistoBins(100,0,250);
+
+ rhoTaskBase = dynamic_cast<AliAnalysisTaskRhoBase*>rhoTask;
+
+ }
+
+ return rhoTaskBase;
+
+}
--- /dev/null
+enum AlgoType {kKT, kANTIKT};
+enum JetType {kFULLJETS, kCHARGEDJETS, kNEUTRALJETS};
+
+AliAnalysisTaskEmcalJetMassBkg* AddTaskEmcalJetMassBkg(const char * njetsBase,
+ const char * njetsTag,
+ const Double_t R,
+ const char * nrhoBase,
+ const char * nrhoTag,
+ const char * ntracks,
+ const char * nclusters
+ const char *type,
+ const char *CentEst,
+ Int_t pSel,
+ TString trigClass = "",
+ TString kEmcalTriggers = "");
+
+AliAnalysisTaskEmcalJetMassBkg* AddTaskEmcalJetMassBkg(TString kTracksName = "PicoTracks",
+ TString kClusName = "caloClustersCorr",
+ Double_t R = 0.4,
+ Double_t ptminTrack = 0.15,
+ Double_t etminClus = 0.3,
+ Double_t ptminTag = 4.,
+ Int_t rhoType = 1,
+ const char *type = "EMCAL",
+ const char *CentEst = "V0M",
+ Int_t pSel = AliVEvent::kCentral | AliVEvent::kSemiCentral | AliVEvent::kMB,
+ TString trigClass = "",
+ TString kEmcalTriggers = "",
+ TString kPeriod = "LHC11h"
+ ) {
+
+ AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+ if (!mgr)
+ {
+ Error("AddTaskEmcalJetMassBkg","No analysis manager found.");
+ return 0;
+ }
+ Bool_t ismc=kFALSE;
+ ismc = (mgr->GetMCtruthEventHandler())?kTRUE:kFALSE;
+
+ // Check the analysis type using the event handlers connected to the analysis manager.
+ //==============================================================================
+ if (!mgr->GetInputEventHandler())
+ {
+ ::Error("AddTaskEmcalJetMassBkg", "This task requires an input event handler");
+ return NULL;
+ }
+
+ // #### Add necessary jet finder tasks
+ gROOT->LoadMacro("$ALICE_ROOT/PWGJE/EMCALJetTasks/macros/AddTaskEmcalJet.C");
+
+ AliEmcalJetTask* jetFinderTaskBase = 0x0;
+ if (strcmp(type,"TPC")==0)
+ jetFinderTaskBase = AddTaskEmcalJet(kTracksName, "", kANTIKT, R, kCHARGEDJETS, ptminTrack, etminClus,0.005,0);
+ else if (strcmp(type,"EMCAL")==0)
+ jetFinderTaskBase = AddTaskEmcalJet(kTracksName, kClusName, kANTIKT, R, kFULLJETS, ptminTrack, etminClus,0.005,0);
+
+ TString strJetsBase = jetFinderTaskBase->GetName();
+
+ AliAnalysisTaskRhoBase *rhoTaskBase;
+ TString rhoNameBase = "";
+ if(rhoType==1) {
+ rhoTaskBase = AttachRhoTask(kPeriod,kTracksName,kClusName,R,ptminTrack,etminClus);
+ rhoTaskBase->SetCentralityEstimator(CentEst);
+ rhoTaskBase->SelectCollisionCandidates(AliVEvent::kAny);
+
+ if (strcmp(type,"TPC")==0)
+ rhoNameBase = rhoTaskBase->GetOutRhoName();
+ if (strcmp(type,"EMCAL")==0)
+ rhoNameBase = rhoTaskBase->GetOutRhoScaledName();
+ }
+
+ //Configure jet tagger task
+ AliAnalysisTaskEmcalJetMassBkg *task = AddTaskEmcalJetMassBkg(jetFinderTaskBase->GetName(),
+ R,
+ rhoNameBase,
+ kTracksName.Data(),
+ kClusName.Data(),
+ type,
+ CentEst,
+ pSel,
+ trigClass,
+ kEmcalTriggers
+ );
+
+ return task;
+
+}
+
+AliAnalysisTaskEmcalJetMassBkg* AddTaskEmcalJetMassBkg(const char * njetsBase,
+ const Double_t R,
+ const char * nrhoBase,
+ const char * ntracks,
+ const char * nclusters,
+ const char *type,
+ const char * CentEst,
+ Int_t pSel,
+ TString trigClass,
+ TString kEmcalTriggers
+ ) {
+
+ AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+ if (!mgr)
+ {
+ Error("AddTaskEmcalJetMassBkg","No analysis manager found.");
+ return 0;
+ }
+ Bool_t ismc=kFALSE;
+ ismc = (mgr->GetMCtruthEventHandler())?kTRUE:kFALSE;
+
+ // Check the analysis type using the event handlers connected to the analysis manager.
+ //==============================================================================
+ if (!mgr->GetInputEventHandler())
+ {
+ ::Error("AddTaskEmcalJetMassBkg", "This task requires an input event handler");
+ return NULL;
+ }
+
+ TString wagonName = Form("JetMassBkg_%s_TC%s",njetsBase,trigClass.Data());
+
+ //Configure jet tagger task
+ AliAnalysisTaskEmcalJetMassBkg *task = new AliAnalysisTaskEmcalJetMassBkg(wagonName.Data());
+
+ task->SetNCentBins(4);
+ task->SetConeRadius(R);
+
+ if (strcmp(type,"TPC")==0)
+ task->SetConeEtaPhiTPC();
+ else if (strcmp(type,"EMCAL")==0)
+ task->SetConeEtaPhiEMCAL();
+
+ AliParticleContainer *trackCont = task->AddParticleContainer(ntracks);
+ AliClusterContainer *clusterCont = task->AddClusterContainer(nclusters);
+
+ task->SetJetContainerBase(0);
+
+ TString strType(type);
+ AliJetContainer *jetContBase = task->AddJetContainer(njetsBase,strType,R);
+ if(jetContBase) {
+ jetContBase->SetRhoName(nrhoBase);
+ jetContBase->ConnectParticleContainer(trackCont);
+ jetContBase->ConnectClusterContainer(clusterCont);
+ jetContBase->SetZLeadingCut(0.98,0.98);
+ task->SetPercAreaCut(0.6, 0);
+ }
+
+ task->SetCaloTriggerPatchInfoName(kEmcalTriggers.Data());
+
+ task->SetCentralityEstimator(CentEst);
+
+ task->SelectCollisionCandidates(pSel);
+
+ task->SetUseAliAnaUtils(kFALSE);
+
+ mgr->AddTask(task);
+
+ //Connnect input
+ mgr->ConnectInput (task, 0, mgr->GetCommonInputContainer() );
+
+ //Connect output
+ TString contName(wagonName);
+ TString outputfile = Form("%s",AliAnalysisManager::GetCommonFileName());
+ AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contName.Data(), TList::Class(),AliAnalysisManager::kOutputContainer,outputfile);
+ mgr->ConnectOutput(task,1,coutput1);
+
+ return task;
+
+}
+
+
+AliAnalysisTaskRhoBase *AttachRhoTask(TString kPeriod ,// = "LHC13b",
+ TString kTracksName ,// = "PicoTracks",
+ TString kClusName ,// = "caloClustersCorr",
+ Double_t R ,// = 0.4,
+ Double_t ptminTrack ,// = 0.15,
+ Double_t etminClus // = 0.3
+ ) {
+
+ AliAnalysisTaskRhoBase *rhoTaskBase;
+
+ kPeriod.ToLower();
+
+ // Add kt jet finder and rho task in case we want background subtraction
+ AliEmcalJetTask *jetFinderKt;
+ AliEmcalJetTask *jetFinderAKt;
+ jetFinderKt = AddTaskEmcalJet(kTracksName, "", kKT, R, kCHARGEDJETS, ptminTrack, etminClus,0.005,0);
+ jetFinderAKt = AddTaskEmcalJet(kTracksName, "", kANTIKT, R, kCHARGEDJETS, ptminTrack, etminClus,0.005,0);
+
+ if(kPeriod.EqualTo("lhc13b") || kPeriod.EqualTo("lhc13c") || kPeriod.EqualTo("lhc13d") || kPeriod.EqualTo("lhc13e") || kPeriod.EqualTo("lhc13f")) {
+
+ jetFinderKt->SetMinJetPt(0.);
+
+ gROOT->LoadMacro("$ALICE_ROOT/PWGJE/EMCALJetTasks/macros/AddTaskRhoSparse.C");
+ TF1 *fScale = new TF1("fScale","1.42",0.,100.); //scale factor for pPb
+ AliAnalysisTaskRhoSparse *rhoTaskSparse = AddTaskRhoSparse(
+ jetFinderKt->GetName(),
+ jetFinderAKt->GetName(),
+ kTracksName,
+ kClusName,
+ Form("RhoSparseR%03d",(int)(100*R)),
+ R,
+ "TPC",
+ 0.01,
+ 0.15,
+ 0,
+ fScale,
+ 0,
+ kTRUE,
+ Form("RhoSparseR%03d",(int)(100*R)),
+ kTRUE
+ );
+ rhoTaskSparse->SetUseAliAnaUtils(kTRUE);
+
+ rhoTaskBase = dynamic_cast<AliAnalysisTaskRhoBase*>rhoTaskSparse;
+ }
+ else if(kPeriod.EqualTo("lhc10h") || kPeriod.EqualTo("lhc11h") ) {
+
+ gROOT->LoadMacro("$ALICE_ROOT/PWGJE/EMCALJetTasks/macros/AddTaskRho.C");
+
+ TF1* sfunc = new TF1("sfunc","[0]*x*x+[1]*x+[2]",-1,100);
+ sfunc->SetParameter(2,1.76458);
+ sfunc->SetParameter(1,-0.0111656);
+ sfunc->SetParameter(0,0.000107296);
+ AliAnalysisTaskRho *rhoTask = AddTaskRho(
+ jetFinderKt->GetName(),
+ kTracksName,
+ kClusName,
+ Form("RhoR%03dptmin%3.0f",(int)(100*R),ptminTrack*1000.0),
+ R,
+ "TPC",
+ 0.01,
+ 0,
+ sfunc,
+ 2,
+ kTRUE);
+ rhoTask->SetHistoBins(100,0,250);
+
+ rhoTaskBase = dynamic_cast<AliAnalysisTaskRhoBase*>rhoTask;
+
+ }
+
+ return rhoTaskBase;
+
+}
const char *nrho = "Rho",
Double_t jetradius = 0.2,
Double_t jetptcut = 1,
- Double_t jetareacut = 0.557,
- UInt_t type = AliAnalysisTaskEmcalJet::kEMCAL,
+ Double_t jetareacut = 0.6,
+ const char *type = "EMCAL",
Int_t leadhadtype = 0,
const char *taskname = "AliAnalysisTaskEmcalJetSample"
)
name += "_";
name += nrho;
}
- if (type == AliJetContainer::kTPC)
- name += "_TPC";
- else if (type == AliJetContainer::kEMCAL)
- name += "_EMCAL";
- else if (type == AliJetContainer::kUser)
- name += "_USER";
+ if (strcmp(type,"")) {
+ name += "_";
+ name += type;
+ }
Printf("name: %s",name.Data());
AliAnalysisTaskEmcalJetSample* jetTask = new AliAnalysisTaskEmcalJetSample(name);
- jetTask->SetAnaType(type);
- jetTask->SetTracksName(ntracks);
- jetTask->SetClusName(nclusters);
- jetTask->SetJetsName(njets);
- jetTask->SetRhoName(nrho);
- jetTask->SetJetRadius(jetradius);
- jetTask->SetJetPtCut(jetptcut);
- jetTask->SetPercAreaCut(jetareacut);
- jetTask->SetLeadingHadronType(leadhadtype);
+
+ AliParticleContainer *trackCont = task->AddParticleContainer(ntracks);
+ AliClusterContainer *clusterCont = task->AddClusterContainer(nclusters);
+
+ TString strType(type);
+ AliJetContainer *jetCont = jetTask->AddJetContainer(njets,strType,R);
+ if(jetCont) {
+ jetCont->SetRhoName(nrho);
+ jetCont->ConnectParticleContainer(trackCont);
+ jetCont->ConnectClusterContainer(clusterCont);
+ jetCont->SetZLeadingCut(0.98,0.98);
+ jetCont->SetPercAreaCut(0.6);
+ jetCont->SetJetPtCut(jetptcut);
+ jetCont->SetLeadingHadronType(leadhadtype);
+ }
//-------------------------------------------------------
// Final settings, pass to manager and set the containers
--- /dev/null
+enum AlgoType {kKT, kANTIKT};
+enum JetType {kFULLJETS, kCHARGEDJETS, kNEUTRALJETS};
+
+AliAnalysisTaskEmcalJetTagger* AddTaskEmcalJetTagger(const char * njetsBase,
+ const char * njetsTag,
+ const Double_t R,
+ const char * nrhoBase,
+ const char * nrhoTag,
+ const char * ntracks,
+ const char * nclusters,
+ const char *type,
+ const char *CentEst,
+ Int_t pSel,
+ TString trigClass = "",
+ TString kEmcalTriggers = "");
+
+AliAnalysisTaskEmcalJetTagger* AddTaskEmcalJetTagger(TString kTracksName = "PicoTracks",
+ TString kClusName = "caloClustersCorr",
+ Double_t R = 0.4,
+ Double_t ptminTrack = 0.15,
+ Double_t etminClus = 0.3,
+ Double_t ptminTag = 4.,
+ Int_t rhoType = 1,
+ const char *type = "EMCAL",
+ const char *CentEst = "V0M",
+ Int_t pSel = AliVEvent::kCentral | AliVEvent::kSemiCentral | AliVEvent::kMB,
+ TString trigClass = "",
+ TString kEmcalTriggers = "",
+ TString kPeriod = "LHC11h"
+ ) {
+
+ AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+ if (!mgr)
+ {
+ Error("AddTaskEmcalJetTagger","No analysis manager found.");
+ return 0;
+ }
+ Bool_t ismc=kFALSE;
+ ismc = (mgr->GetMCtruthEventHandler())?kTRUE:kFALSE;
+
+ // Check the analysis type using the event handlers connected to the analysis manager.
+ //==============================================================================
+ if (!mgr->GetInputEventHandler())
+ {
+ ::Error("AddTaskEmcalJetTagger", "This task requires an input event handler");
+ return NULL;
+ }
+
+ // #### Add necessary jet finder tasks
+ gROOT->LoadMacro("$ALICE_ROOT/PWGJE/EMCALJetTasks/macros/AddTaskEmcalJet.C");
+
+ AliEmcalJetTask* jetFinderTaskBase = 0x0;
+ if (strcmp(type,"TPC")==0)
+ jetFinderTaskBase = AddTaskEmcalJet(kTracksName, "", kANTIKT, R, kCHARGEDJETS, ptminTrack, etminClus,0.005,0);
+ else if (strcmp(type,"EMCAL")==0)
+ jetFinderTaskBase = AddTaskEmcalJet(kTracksName, kClusName, kANTIKT, R, kFULLJETS, ptminTrack, etminClus,0.005,0);
+
+ AliEmcalJetTask* jetFinderTaskTag = AddTaskEmcalJet(kTracksName, "", kANTIKT, R, kCHARGEDJETS, ptminTag, etminClus,0.005,0);
+
+ TString strJetsBase = jetFinderTaskBase->GetName();
+ TString strJetsTag = jetFinderTaskTag->GetName();
+
+ AliAnalysisTaskRhoBase *rhoTaskBase;
+ AliAnalysisTaskRhoBase *rhoTaskTag;
+ TString rhoNameBase = "";
+ TString rhoNameTag = "";
+ if(rhoType==1) {
+ rhoTaskBase = AttachRhoTaskTagger(kPeriod,kTracksName,kClusName,R,ptminTrack,etminClus);
+ rhoTaskBase->SetCentralityEstimator(CentEst);
+ rhoTaskBase->SelectCollisionCandidates(AliVEvent::kAny);
+
+ rhoTaskTag = AttachRhoTaskTagger(kPeriod,kTracksName,kClusName,R,ptminTag,0.);
+ rhoTaskTag->SetCentralityEstimator(CentEst);
+ rhoTaskTag->SelectCollisionCandidates(AliVEvent::kAny);
+
+ if (strcmp(type,"TPC")==0)
+ rhoNameBase = rhoTaskBase->GetOutRhoName();
+ if (strcmp(type,"EMCAL")==0)
+ rhoNameBase = rhoTaskBase->GetOutRhoScaledName();
+ rhoNameTag = rhoTaskTag->GetOutRhoName();
+ }
+
+ //Configure jet tagger task
+ AliAnalysisTaskEmcalJetTagger *task = AddTaskEmcalJetTagger(jetFinderTaskBase->GetName(),
+ jetFinderTaskTag->GetName(),
+ R,
+ rhoNameBase,
+ rhoNameTag,
+ kTracksName.Data(),
+ kClusName.Data(),
+ type,
+ CentEst,
+ pSel,
+ trigClass,
+ kEmcalTriggers
+ );
+
+ return task;
+
+}
+
+AliAnalysisTaskEmcalJetTagger* AddTaskEmcalJetTagger(const char * njetsBase,
+ const char * njetsTag,
+ const Double_t R,
+ const char * nrhoBase,
+ const char * nrhoTag,
+ const char * ntracks,
+ const char * nclusters,
+ const char * type,
+ const char * CentEst,
+ Int_t pSel,
+ TString trigClass,
+ TString kEmcalTriggers
+ ) {
+
+ AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+ if (!mgr)
+ {
+ Error("AddTaskEmcalJetTagger","No analysis manager found.");
+ return 0;
+ }
+ Bool_t ismc=kFALSE;
+ ismc = (mgr->GetMCtruthEventHandler())?kTRUE:kFALSE;
+
+ // Check the analysis type using the event handlers connected to the analysis manager.
+ //==============================================================================
+ if (!mgr->GetInputEventHandler())
+ {
+ ::Error("AddTaskEmcalJetTagger", "This task requires an input event handler");
+ return NULL;
+ }
+
+ TString wagonName = Form("JetTagger_%s_%s_TC%s",njetsBase,njetsTag,trigClass.Data());
+
+ //Configure jet tagger task
+ AliAnalysisTaskEmcalJetTagger *task = new AliAnalysisTaskEmcalJetTagger(wagonName.Data());
+
+ task->SetNCentBins(4);
+ //task->SetVzRange(-10.,10.);
+
+ AliParticleContainer *trackCont = task->AddParticleContainer(ntracks);
+ AliClusterContainer *clusterCont = task->AddClusterContainer(nclusters);
+
+ task->SetJetContainerBase(0);
+ task->SetJetContainerTag(1);
+
+ TString strType(type);
+ AliJetContainer *jetContBase = task->AddJetContainer(njetsBase,strType,R);
+ if(jetContBase) {
+ jetContBase->SetRhoName(nrhoBase);
+ jetContBase->ConnectParticleContainer(trackCont);
+ jetContBase->ConnectClusterContainer(clusterCont);
+ jetContBase->SetZLeadingCut(0.98,0.98);
+ }
+
+ AliJetContainer *jetContTag = task->AddJetContainer(njetsTag,"TPC",R);
+ if(jetContTag) {
+ jetContTag->SetRhoName(nrhoTag);
+ jetContTag->ConnectParticleContainer(trackCont);
+ jetContTag->ConnectClusterContainer(clusterCont);
+ }
+
+ for(Int_t i=0; i<2; i++) {
+ task->SetPercAreaCut(0.6, i);
+ }
+
+ task->SetCaloTriggerPatchInfoName(kEmcalTriggers.Data());
+
+ task->SetCentralityEstimator(CentEst);
+
+ task->SelectCollisionCandidates(pSel);
+
+ task->SetUseAliAnaUtils(kFALSE);
+
+ mgr->AddTask(task);
+
+ //Connnect input
+ mgr->ConnectInput (task, 0, mgr->GetCommonInputContainer() );
+
+ //Connect output
+ TString contName(wagonName);
+ TString outputfile = Form("%s",AliAnalysisManager::GetCommonFileName());
+ AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contName.Data(), TList::Class(),AliAnalysisManager::kOutputContainer,outputfile);
+ mgr->ConnectOutput(task,1,coutput1);
+
+ return task;
+
+}
+
+
+AliAnalysisTaskRhoBase *AttachRhoTaskTagger(TString kPeriod = "LHC13b",
+ TString kTracksName = "PicoTracks",
+ TString kClusName = "caloClustersCorr",
+ Double_t R = 0.4,
+ Double_t ptminTrack = 0.15,
+ Double_t etminClus = 0.3
+ ) {
+
+ AliAnalysisTaskRhoBase *rhoTaskBase;
+
+ kPeriod.ToLower();
+
+ // Add kt jet finder and rho task in case we want background subtraction
+ AliEmcalJetTask *jetFinderKt;
+ AliEmcalJetTask *jetFinderAKt;
+ jetFinderKt = AddTaskEmcalJet(kTracksName, "", kKT, R, kCHARGEDJETS, ptminTrack, etminClus,0.005,0);
+ jetFinderAKt = AddTaskEmcalJet(kTracksName, "", kANTIKT, R, kCHARGEDJETS, ptminTrack, etminClus,0.005,0);
+
+ if(kPeriod.EqualTo("lhc13b") || kPeriod.EqualTo("lhc13c") || kPeriod.EqualTo("lhc13d") || kPeriod.EqualTo("lhc13e") || kPeriod.EqualTo("lhc13f")) {
+
+
+ jetFinderKt->SetMinJetPt(0.);
+
+ gROOT->LoadMacro("$ALICE_ROOT/PWGJE/EMCALJetTasks/macros/AddTaskRhoSparse.C");
+ TF1 *fScale = new TF1("fScale","1.42",0.,100.); //scale factor for pPb
+ AliAnalysisTaskRhoSparse *rhoTaskSparse = AddTaskRhoSparse(
+ jetFinderKt->GetName(),
+ jetFinderAKt->GetName(),
+ kTracksName,
+ kClusName,
+ Form("RhoSparseR%03d",(int)(100*R)),
+ R,
+ "TPC",
+ 0.01,
+ 0.15,
+ 0,
+ fScale,
+ 0,
+ kTRUE,
+ Form("RhoSparseR%03d",(int)(100*R)),
+ kTRUE
+ );
+ rhoTaskSparse->SetUseAliAnaUtils(kTRUE);
+
+ rhoTaskBase = dynamic_cast<AliAnalysisTaskRhoBase*>rhoTaskSparse;
+ }
+ else if(kPeriod.EqualTo("lhc10h") || kPeriod.EqualTo("lhc11h") ) {
+
+ gROOT->LoadMacro("$ALICE_ROOT/PWGJE/EMCALJetTasks/macros/AddTaskRho.C");
+
+ TF1* sfunc = new TF1("sfunc","[0]*x*x+[1]*x+[2]",-1,100);
+ sfunc->SetParameter(2,1.76458);
+ sfunc->SetParameter(1,-0.0111656);
+ sfunc->SetParameter(0,0.000107296);
+ AliAnalysisTaskRho *rhoTask = AddTaskRho(
+ jetFinderKt->GetName(),
+ kTracksName,
+ kClusName,
+ Form("RhoR%03dptmin%3.0f",(int)(100*R),ptminTrack*1000.0),
+ R,
+ "TPC",
+ 0.01,
+ 0,
+ sfunc,
+ 2,
+ kTRUE);
+ rhoTask->SetHistoBins(100,0,250);
+
+ rhoTaskBase = dynamic_cast<AliAnalysisTaskRhoBase*>rhoTask;
+
+ }
+
+ return rhoTaskBase;
+
+}
if (useTender)
{
gROOT->LoadMacro("$ALICE_ROOT/PWG/EMCAL/macros/AddTaskEMCALTender.C");
- AliAnalysisTaskSE *tender = AddTaskEMCALTender(runPeriod, kFALSE, kTRUE, kTRUE, kTRUE, kTRUE, kTRUE, kFALSE, kTRUE, kTRUE, kTRUE,
- kFALSE, AliEMCALRecParam::kClusterizerNxN);
+ AliAnalysisTaskSE *tender = AddTaskEMCALTender(runPeriod, kTRUE, kTRUE, kTRUE, kTRUE, kTRUE, kFALSE, kTRUE, kTRUE, kTRUE,
+ AliEMCALRecoUtils::kBeamTestCorrected,kTRUE,0.1,0.05,AliEMCALRecParam::kClusterizerv2,
+ kFALSE,kFALSE,-1,1e6,1e6);
if (usedData != "AOD" && !useGrid) {
AliTender *alitender = dynamic_cast<AliTender*>(tender);
alitender->SetDefaultCDBStorage("local://$ALICE_ROOT/OCDB");
#pragma link C++ class AliAnalysisTaskEmcalDiJetResponse+;
#pragma link C++ class AliAnalysisTaskEmcalJetHMEC+;
#pragma link C++ class AliAnalysisTaskEmcalJetHadCorQA+;
+#pragma link C++ class AliAnalysisTaskEmcalJetMass+;
+#pragma link C++ class AliAnalysisTaskEmcalJetMassBkg+;
#pragma link C++ class AliAnalysisTaskEmcalJetSample+;
#pragma link C++ class AliAnalysisTaskEmcalJetSpectra+;
#pragma link C++ class AliAnalysisTaskEmcalJetSpectraMECpA;
+#pragma link C++ class AliAnalysisTaskEmcalJetTagger+;
#pragma link C++ class AliAnalysisTaskEmcalJetTriggerQA+;
#pragma link C++ class AliAnalysisTaskEmcalTriggerInfoQA+;
#pragma link C++ class AliAnalysisTaskFullpAJets+;
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff