EMCALJetTasks/AliRhoParameter.cxx
EMCALJetTasks/AliLocalRhoParameter.cxx
EMCALJetTasks/AliJetTriggerSelectionTask.cxx
+ EMCALJetTasks/AliStackPartonInfo.cxx
EMCALJetTasks/UserTasks/AliAnalysisTaskCLQA.cxx
EMCALJetTasks/UserTasks/AliAnalysisTaskChargedJetsPA.cxx
EMCALJetTasks/UserTasks/AliAnalysisTaskDcalDijetPerf.cxx
EMCALJetTasks/UserTasks/AliAnalysisTaskRhoMass.cxx
EMCALJetTasks/UserTasks/AliAnalysisTaskRhoMassBase.cxx
EMCALJetTasks/UserTasks/AliAnalysisTaskRhoMassScale.cxx
+ EMCALJetTasks/UserTasks/AliAnalysisTaskEmcalQGTagging.cxx
EMCALJetTasks/UserTasks/AliAnalysisTaskSAJF.cxx
EMCALJetTasks/UserTasks/AliAnalysisTaskSAQA.cxx
EMCALJetTasks/UserTasks/AliAnalysisTaskSOH.cxx
cont->SetRunNumber(InputEvent()->GetRunNumber());
cont->SetEMCALGeometry();
cont->SetArray(InputEvent());
- cont->LoadRho(InputEvent());
+ cont->LoadRho(InputEvent());
+ cont->LoadPartonsInfo(InputEvent());
}
//Get Jets, cuts and rho for first jet container
AliEmcalJet *jet_sub = new ((*fJetsSub)[jetCount])
AliEmcalJet(jets_sub[ijet].perp(), jets_sub[ijet].eta(), jets_sub[ijet].phi(), jets_sub[ijet].m());
jet_sub->SetLabel(ijet);
- // Fill constituent info
+ // Fill constituent info
std::vector<fastjet::PseudoJet> constituents_sub(fjw.GetJetConstituents(ijet));
jet_sub->SetNumberOfTracks(constituents_sub.size());
jet_sub->SetNumberOfClusters(constituents_sub.size());
} else {
AliError(Form("%s: No logical way to end up here.", GetName()));
continue;
- }
+ }
}
}
//______________________________________________________________________________________
#include "AliParticleContainer.h"
#include "AliClusterContainer.h"
#include "AliLocalRhoParameter.h"
+#include "AliStackPartonInfo.h"
#include "AliJetContainer.h"
fRhoName(),
fLocalRhoName(),
fRhoMassName(),
+ fPartonInfoName(),
fFlavourSelection(0),
fPtBiasJetTrack(0),
fPtBiasJetClus(0),
fRho(0),
fLocalRho(0),
fRhoMass(0),
+ fPartonsInfo(0),
fGeom(0),
fRunNumber(0)
{
fRhoName(),
fLocalRhoName(),
fRhoMassName(),
+ fPartonInfoName(),
fFlavourSelection(0),
fPtBiasJetTrack(0),
fPtBiasJetClus(0),
fRho(0),
fLocalRho(0),
fRhoMass(0),
+ fPartonsInfo(0),
fGeom(0),
fRunNumber(0)
{
}
}
}
+//________________________________________________________________________
+void AliJetContainer::LoadPartonsInfo(AliVEvent *event)
+{
+ // Load parton info
+
+ if (!fPartonInfoName.IsNull() && !fPartonsInfo) {
+ fPartonsInfo = dynamic_cast<AliStackPartonInfo*>(event->FindListObject(fPartonInfoName));
+ if (!fPartonsInfo) {
+ AliError(Form("%s: Could not retrieve parton infos! %s!", GetName(), fPartonInfoName.Data())); return;
+ }
+ }
+}
+
//________________________________________________________________________
class AliParticleContainer;
class AliClusterContainer;
class AliLocalRhoParameter;
+class AliStackPartonInfo;
#include "AliRhoParameter.h"
void LoadRho(AliVEvent *event);
void LoadLocalRho(AliVEvent *event);
void LoadRhoMass(AliVEvent *event);
+ void LoadPartonsInfo(AliVEvent *event);
void SetJetAcceptanceType(JetAcceptanceType type) { fJetAcceptanceType = type ; }
void PrintCuts();
virtual void SetRhoName(const char *n) { fRhoName = n ; }
virtual void SetLocalRhoName(const char *n) { fLocalRhoName = n ; }
virtual void SetRhoMassName(const char *n) { fRhoMassName = n ; }
+ virtual void SetPartonInfoName(const char *n) { fPartonInfoName = n; }
+
void ConnectParticleContainer(AliParticleContainer *c) { fParticleContainer = c ; }
void ConnectClusterContainer(AliClusterContainer *c) { fClusterContainer = c ; }
AliRhoParameter *GetRhoMassParameter() {return fRhoMass;}
Double_t GetRhoMassVal() const {if (fRhoMass) return fRhoMass->GetVal(); else return 0;}
const TString& GetRhoMassName() const {return fRhoMassName;}
+ const TString& GetPartonInfoName() const {return fPartonInfoName;}
+ AliStackPartonInfo *GetStackPartonInfo() const {return fPartonsInfo;}
Double_t GetJetPtCorr(Int_t i) const;
Double_t GetJetPtCorrLocal(Int_t i) const;
Float_t GetJetRadius() const {return fJetRadius;}
AliParticleContainer *GetParticleContainer() {return fParticleContainer;}
AliClusterContainer *GetClusterContainer() {return fClusterContainer;}
Double_t GetFractionSharedPt(AliEmcalJet *jet) const;
+
protected:
JetAcceptanceType fJetAcceptanceType; // acceptance type
TString fRhoName; // Name of rho object
TString fLocalRhoName; // Name of local rho object
TString fRhoMassName; // Name of rho mass object
+ TString fPartonInfoName; // Name of parton info object
Int_t fFlavourSelection; // selection on jet flavour
Float_t fPtBiasJetTrack; // select jets with a minimum pt track
Float_t fPtBiasJetClus; // select jets with a minimum pt cluster
AliRhoParameter *fRho; //! event rho for these jets
AliLocalRhoParameter *fLocalRho; //! event local rho for these jets
AliRhoParameter *fRhoMass; //! event rho mass for these jets
+ AliStackPartonInfo *fPartonsInfo; //! event parton info
AliEMCALGeometry *fGeom; //! emcal geometry
Int_t fRunNumber; //! run number
#include "AliVCluster.h"
#include "AliVEvent.h"
#include "AliGenPythiaEventHeader.h"
+#include "AliStackPartonInfo.h"
ClassImp(AliJetEmbeddingFromGenTask)
fHistPt = new TH1F("fHistpt","fHistPt;#it{p}_{T};N",100,0.,100.);
fOutput->Add(fHistPt);
- fHistEtaPhi = new TH2F("fHistEtaPhi","fHistEtaPhi;#eta;#varphi",100,-3.,3.,100.,0.,TMath::TwoPi());
+ fHistEtaPhi = new TH2F("fHistEtapHI","fHistEtaPhi;#eta;#varphi",100,-3.,3.,100.,0.,TMath::TwoPi());
fOutput->Add(fHistEtaPhi);
fHistTrials = new TH1F("fHistTrials", "fHistTrials", 1, 0, 1);
InputEvent()->AddObject(fOutTracks);
fNTracks = 0;
}
-
+
+ if (!(InputEvent()->FindListObject(fPartonInfoName))) {
+ fStackPartonInfo = new AliStackPartonInfo("PartonsInfo");
+ fStackPartonInfo->SetName(fPartonInfoName);
+ InputEvent()->AddObject(fStackPartonInfo);
+ }
return kTRUE;
}
stack->Reset();
fGen->Generate();
const Int_t nprim = stack->GetNprimary();
+ TParticle *part6 = stack->Particle(6);
+ TParticle *part7 = stack->Particle(7);
+
+ fStackPartonInfo->SetPartonFlag6(TMath::Abs(part6->GetPdgCode()));
+ fStackPartonInfo->SetPartonPt6(part6->Pt());
+ fStackPartonInfo->SetPartonEta6(part6->Eta());
+ fStackPartonInfo->SetPartonPhi6(part6->Phi());
+
+ fStackPartonInfo->SetPartonFlag7(TMath::Abs(part7->GetPdgCode()));
+ fStackPartonInfo->SetPartonPt7(part7->Pt());
+ fStackPartonInfo->SetPartonEta7(part7->Eta());
+ fStackPartonInfo->SetPartonPhi7(part7->Phi());
+
for (Int_t i=0;i<nprim;++i) {
if (!stack->IsPhysicalPrimary(i))
continue;
fOutMCParticlesMap(0),
fMCLabelShift(0),
fEsdMode(kFALSE),
- fOutput(0)
+ fOutput(0),
+ fStackPartonInfo(0x0)
{
// Default constructor.
fOutMCParticlesMap(0),
fMCLabelShift(0),
fEsdMode(kFALSE),
- fOutput(0)
+ fOutput(0),
+ fStackPartonInfo(0x0)
{
// Standard constructor.
if (fPtPhiEvPlDistribution || fAddV2)
fPsi = gRandom->Rndm() * TMath::Pi();
-
+
Run();
-
+
if (fCaloCells && !fCopyArray) {
delete fCaloCells;
fCaloCells = tempCaloCells;
class AliAODMCParticle;
class AliNamedArrayI;
class TF2;
+class AliStackPartonInfo;
#include <TH1F.h>
#include <TF1.h>
void SetEtaRange(Float_t min, Float_t max) { fEtaMin = min; fEtaMax = max; }
void SetPhiRange(Float_t min, Float_t max) { fPhiMin = min; fPhiMax = max; }
void SetPtRange(Float_t min, Float_t max) { fPtMin = min; fPtMax = max; }
- void SetPtSpectrum(TH1F *f) { fPtSpectrum = f; }
+ void SetPtSpectrum(TH1 *f) { fPtSpectrum = f; }
void SetPtSpectrum(TF1 *f) { fPtSpectrum = new TH1F("ptSpectrum","ptSpectrum",1000,f->GetXmin(),f->GetXmax());
fPtSpectrum->Add(f); }
void SetPtPhiEvPlDistribution(TF2 *f) { fPtPhiEvPlDistribution = f; }
- void SetDensitySpectrum(TH1F *f) { fDensitySpectrum = f; }
+ void SetDensitySpectrum(TH1 *f) { fDensitySpectrum = f; }
void SetDensitySpectrum(TF1 *f) { fDensitySpectrum = new TH1F("densitypectrum","densitypectrum",1000,f->GetXmin(),f->GetXmax());
fDensitySpectrum->Add(f); }
void SetDifferentialV2(TF1* f) { fDifferentialV2 = f; }
void SetClusName(const char *n) { fCaloName = n; }
void SetCellsName(const char *n) { fCellsName = n; }
void SetMCParticlesName(const char *n) { fMCParticlesName = n; }
+ void SetPartonInfoName(const char *n) { fPartonInfoName = n; }
void SetSuffix(const char *s) { fSuffix = s; }
void SetGeometryName(const char *n) { fGeomName = n; }
void SetMarkMC(Int_t m) { fMarkMC = m; }
TString fOutCellsName; // name of output cells collection
TString fMCParticlesName; // name of MC particle collection
TString fOutMCParticlesName; // name of output MC particle collection
+ TString fPartonInfoName; // name of partons info
Bool_t fIsMC; // whether the current event is MC or not
TString fSuffix; // suffix to add in the name of new collections
Float_t fEtaMin; // eta minimum value
Int_t fNCells; // how many cells are being processed
Int_t fNTracks; // how many tracks are being processed
Int_t fMarkMC; // which MC label is to be used (default=100)
- TH1F *fPtSpectrum; // pt spectrum to extract random pt values
+ TH1 *fPtSpectrum; // pt spectrum to extract random pt values
TF2 *fPtPhiEvPlDistribution; // pt vs. (phi-psi) distribution to extract random pt/phi values
- TH1F *fDensitySpectrum; // particle density spectrum to extract random density values
+ TH1 *fDensitySpectrum; // particle density spectrum to extract random density values
TF1 *fDifferentialV2; // v2 as function of pt
Bool_t fAddV2; // add v2 sampled from a tf1
Bool_t fFlowFluctuations; // introduce gaussian flow fluctuation
Int_t fMCLabelShift; //!MC label shift
Bool_t fEsdMode; //!ESD/AOD mode
TList *fOutput; //!output list for QA histograms
+ AliStackPartonInfo *fStackPartonInfo; //!Info on original partons:PDG,pt, eta, phi
private:
AliJetModelBaseTask(const AliJetModelBaseTask&); // not implemented
AliJetModelBaseTask &operator=(const AliJetModelBaseTask&); // not implemented
- ClassDef(AliJetModelBaseTask, 12) // Jet modelling task
+ ClassDef(AliJetModelBaseTask, 11) // Jet modelling task
};
#endif
--- /dev/null
+#include "AliStackPartonInfo.h"
+
+ClassImp(AliStackPartonInfo)
+
+AliStackPartonInfo::AliStackPartonInfo() :
+ TNamed(),
+ fPartonFlag6(0),
+ fPartonPt6(0),
+ fPartonEta6(0),
+ fPartonPhi6(0),
+ fPartonFlag7(0),
+ fPartonPt7(0),
+ fPartonEta7(0),
+ fPartonPhi7(0)
+{
+
+}
+
+//_______________________________________________
+
+AliStackPartonInfo::AliStackPartonInfo(const char* name) :
+ TNamed(),
+ fPartonFlag6(0),
+ fPartonPt6(0),
+ fPartonEta6(0),
+ fPartonPhi6(0),
+ fPartonFlag7(0),
+ fPartonPt7(0),
+ fPartonEta7(0),
+ fPartonPhi7(0)
+{
+
+}
+
--- /dev/null
+#ifndef ALISTACKPARTONINFO_H
+#define ALISTACKPARTONINFO_H
+
+#include <TMath.h>
+#include <TNamed.h>
+
+class AliStackPartonInfo : public TNamed{
+
+ public:
+ AliStackPartonInfo();
+ AliStackPartonInfo(const char* name);
+
+ void SetPartonFlag6(Int_t flag6) {fPartonFlag6 = flag6;}
+ void SetPartonPt6(Float_t pt6) {fPartonPt6 = pt6;}
+ void SetPartonEta6(Float_t eta6) {fPartonEta6 = eta6;}
+ void SetPartonPhi6(Float_t phi6) {fPartonPhi6 = phi6;}
+
+ void SetPartonFlag7(Int_t flag7) {fPartonFlag7 = flag7;}
+ void SetPartonPt7(Float_t pt7) {fPartonPt7 = pt7;}
+ void SetPartonEta7(Float_t eta7) {fPartonEta7 = eta7;}
+ void SetPartonPhi7(Float_t phi7) {fPartonPhi7 = phi7;}
+
+
+ Int_t GetPartonFlag6() {return fPartonFlag6;}
+ Float_t GetPartonPt6() {return fPartonPt6;}
+ Float_t GetPartonEta6() {return fPartonEta6;}
+ Float_t GetPartonPhi6() {return fPartonPhi6;}
+
+ Int_t GetPartonFlag7() {return fPartonFlag7;}
+ Float_t GetPartonPt7() {return fPartonPt7;}
+ Float_t GetPartonEta7() {return fPartonEta7;}
+ Float_t GetPartonPhi7() {return fPartonPhi7;}
+
+ private:
+ Int_t fPartonFlag6; //! parton flag
+ Float_t fPartonPt6; //! pT parton
+ Float_t fPartonEta6; //!eta parton
+ Float_t fPartonPhi6; //! phi parton
+
+ Int_t fPartonFlag7; //! parton flag
+ Float_t fPartonPt7; //! pT parton
+ Float_t fPartonEta7; //!eta parton
+ Float_t fPartonPhi7; //! phi parton
+
+ AliStackPartonInfo(const AliStackPartonInfo&);
+ AliStackPartonInfo& operator=(const AliStackPartonInfo&);
+
+ ClassDef(AliStackPartonInfo, 1);
+
+};
+#endif
--- /dev/null
+//
+// Jet QG tagging analysis task.
+//
+// Author: D. Caffarri, L. Cunqueiro
+
+#include <TClonesArray.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TH3F.h>
+#include <THnSparse.h>
+#include <TTree.h>
+#include <TList.h>
+#include <TLorentzVector.h>
+#include <TProfile.h>
+#include <TChain.h>
+#include <TSystem.h>
+#include <TFile.h>
+#include <TKey.h>
+#include "TMatrixD.h"
+#include "TMatrixDSym.h"
+#include "TMatrixDSymEigen.h"
+#include "TVector3.h"
+#include "TVector2.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 "AliParticleContainer.h"
+#include "AliStackPartonInfo.h"
+
+
+#include "AliAODEvent.h"
+
+#include "AliAnalysisTaskEmcalQGTagging.h"
+
+ClassImp(AliAnalysisTaskEmcalQGTagging)
+
+//________________________________________________________________________
+AliAnalysisTaskEmcalQGTagging::AliAnalysisTaskEmcalQGTagging() :
+ AliAnalysisTaskEmcalJet("AliAnalysisTaskEmcalQGTagging", kTRUE),
+ fContainer(0),
+ fMinFractionShared(0),
+ fJetShapeType(kRaw),
+ fShapesVar(0),
+ fIsMC(kFALSE),
+ fIsEmbedding(kFALSE),
+ fIsConstSub(kFALSE),
+ fRMatching(0.3),
+ fPhiJetCorr6(0x0),
+ fPhiJetCorr7(0x0),
+ fEtaJetCorr6(0x0),
+ fEtaJetCorr7(0x0),
+ fPtJetCorr(0x0),
+ fPtJet(0x0),
+ fTreeObservableTagging(0)
+{
+ SetMakeGeneralHistograms(kTRUE);
+}
+
+//________________________________________________________________________
+AliAnalysisTaskEmcalQGTagging::AliAnalysisTaskEmcalQGTagging(const char *name) :
+ AliAnalysisTaskEmcalJet(name, kTRUE),
+ fContainer(0),
+ fMinFractionShared(0),
+ fJetShapeType(kRaw),
+ fShapesVar(0),
+ fIsMC(kFALSE),
+ fIsEmbedding(kFALSE),
+ fIsConstSub(kFALSE),
+ fRMatching(0.3),
+ fPhiJetCorr6(0x0),
+ fPhiJetCorr7(0x0),
+ fEtaJetCorr6(0x0),
+ fEtaJetCorr7(0x0),
+ fPtJetCorr(0x0),
+ fPtJet(0x0),
+ fTreeObservableTagging(0)
+{
+ // Standard constructor.
+
+ SetMakeGeneralHistograms(kTRUE);
+
+ DefineOutput(1, TTree::Class());
+
+}
+
+//________________________________________________________________________
+AliAnalysisTaskEmcalQGTagging::~AliAnalysisTaskEmcalQGTagging()
+{
+ // Destructor.
+}
+
+//________________________________________________________________________
+ void AliAnalysisTaskEmcalQGTagging::UserCreateOutputObjects()
+{
+ // Create user output.
+
+ AliAnalysisTaskEmcalJet::UserCreateOutputObjects();
+
+ Bool_t oldStatus = TH1::AddDirectoryStatus();
+ TH1::AddDirectory(kFALSE);
+
+ fTreeObservableTagging = new TTree("fTreeJetShape", "fTreeJetShape");
+ Int_t nVar = 9;
+ fShapesVar = new Float_t [nVar];
+ TString *fShapesVarNames = new TString [nVar];
+
+ fShapesVarNames[0] = "partonCode";
+ fShapesVarNames[1] = "ptJet";
+ fShapesVarNames[2] = "ptDJet";
+ fShapesVarNames[3] = "mJet";
+ fShapesVarNames[4] = "nbOfConst";
+ fShapesVarNames[5] = "angularity";
+ fShapesVarNames[6] = "circularity";
+ fShapesVarNames[7] = "lesub";
+ fShapesVarNames[8] = "sigma2";
+
+ for(Int_t ivar=0; ivar < nVar; ivar++){
+ cout<<"looping over variables"<<endl;
+ fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data()));
+
+ //if( ivar == 4 ) fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/I", fShapesVarNames[ivar].Data()));
+
+ }
+
+ fPhiJetCorr6= new TH2F("fPhiJetCorr6", "fPhiJetCorr6", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
+ fOutput->Add(fPhiJetCorr6);
+ fEtaJetCorr6= new TH2F("fEtaJetCorr6", "fEtaJetCorr6", 50, -1.5, 1.5, 50, -1.5, 1.5);
+ fOutput->Add(fEtaJetCorr6);
+
+ fPhiJetCorr7= new TH2F("fPhiJetCorr7", "fPhiJetCorr7", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
+ fOutput->Add(fPhiJetCorr7);
+ fEtaJetCorr7= new TH2F("fEtaJetCorr7", "fEtaJetCorr7", 50, -1.5, 1.5, 50, -1.5, 1.5);
+ fOutput->Add(fEtaJetCorr7);
+
+ fPtJetCorr= new TH2F("fPtJetCorr", "fPtJetCorr", 100, 0, 200, 100, 0, 200);
+ fOutput->Add(fPtJetCorr);
+ fPtJet= new TH1F("fPtJet", "fPtJet", 100, 0, 200);
+ fOutput->Add(fPtJet);
+
+
+ fOutput->Add(fTreeObservableTagging);
+ TH1::AddDirectory(oldStatus);
+ PostData(1, fOutput); // Post data for ALL output slots > 0 here.
+
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskEmcalQGTagging::Run()
+{
+ // Run analysis code here, if needed. It will be executed before FillHistograms().
+
+ return kTRUE;
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskEmcalQGTagging::FillHistograms()
+{
+ // Fill histograms.
+ //cout<<"base container"<<endl;
+ AliEmcalJet* jet1 = NULL;
+ AliJetContainer *jetCont = GetJetContainer(0);
+
+ if(jetCont) {
+ jetCont->ResetCurrentID();
+ while((jet1 = jetCont->GetNextAcceptJet())) {
+ if (!jet1) continue;
+ fPtJet->Fill(jet1->Pt());
+ if (fIsMC) {
+ AliStackPartonInfo *partonsInfo = 0x0;
+ AliEmcalJet* jet2 = 0x0;
+ if (fIsEmbedding){
+ AliJetContainer *jetContTrue = GetJetContainer(1);
+ jet2 = jet1->ClosestJet();
+ if (!jet2) {
+ Printf("jet2 not exists, returning");
+ continue;
+ }
+
+ Double_t fraction = jetCont->GetFractionSharedPt(jet1);
+ if(fraction<fMinFractionShared) continue;
+ partonsInfo = (AliStackPartonInfo*) jetContTrue->GetStackPartonInfo();
+
+ }
+ else {
+ partonsInfo = (AliStackPartonInfo*) jetCont->GetStackPartonInfo();
+ jet2=jet1;
+ }
+
+ Double_t jp1=(jet2->Phi())-(partonsInfo->GetPartonPhi6());
+ Double_t detap1=(jet2->Eta())-(partonsInfo->GetPartonEta6());
+
+ if (jp1< -1*TMath::Pi()) jp1 = (-2*TMath::Pi())-jp1;
+ else if (jp1 > TMath::Pi()) jp1 = (2*TMath::Pi())-jp1;
+ Float_t dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
+ fEtaJetCorr6->Fill(jet2->Eta(), partonsInfo->GetPartonEta6());
+ fPhiJetCorr6->Fill(jet2->Phi(), partonsInfo->GetPartonPhi6());
+ if(dRp1 < fRMatching) {
+ fShapesVar[0] = partonsInfo->GetPartonFlag6();
+ fPtJetCorr ->Fill(partonsInfo->GetPartonPt6(), jet2->Pt());
+ }
+ else {
+ jp1=(jet2->Phi())-(partonsInfo->GetPartonPhi7());
+ detap1=(jet2->Eta())-(partonsInfo->GetPartonEta7());
+ if (jp1< -1*TMath::Pi()) jp1= (-2*TMath::Pi())-jp1;
+ else if (jp1 > TMath::Pi()) jp1 = (2*TMath::Pi())-jp1;
+ dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
+ fEtaJetCorr7->Fill(jet2->Eta(), partonsInfo->GetPartonEta7());
+ fPhiJetCorr7->Fill(jet2->Phi(), partonsInfo->GetPartonPhi7());
+ if(dRp1 < fRMatching) {
+ fShapesVar[0] = partonsInfo->GetPartonFlag7();
+ fPtJetCorr ->Fill(partonsInfo->GetPartonPt7(), jet2->Pt());
+ }
+ else continue;
+ }
+ }
+ else
+ fShapesVar[0] = 0.;
+
+
+ if ((fJetShapeType==AliAnalysisTaskEmcalQGTagging::kRaw && fIsConstSub==kFALSE) || (fJetShapeType==AliAnalysisTaskEmcalQGTagging::kDeriv)) fShapesVar[1] = jet1->Pt() - GetRhoVal(0)*jet1->Area();
+ else fShapesVar[1] = jet1->Pt();
+
+ fShapesVar[2] = GetJetpTD(jet1);
+ fShapesVar[3] = GetJetMass(jet1);
+ fShapesVar[4] = 1.*GetJetNumberOfConstituents(jet1);
+ fShapesVar[5] = GetJetAngularity(jet1);
+ fShapesVar[6] = GetJetCircularity(jet1);
+ fShapesVar[7] = GetJetLeSub(jet1);
+ fShapesVar[8] = GetSigma2(jet1);
+ fTreeObservableTagging->Fill();
+
+ }
+
+ }
+
+ return kTRUE;
+}
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::GetJetMass(AliEmcalJet *jet) {
+ //calc subtracted jet mass
+ if(fJetShapeType==kDeriv)
+ return jet->GetSecondOrderSubtracted();
+ else
+ return jet->M();
+}
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::Angularity(AliEmcalJet *jet){
+
+ AliJetContainer *jetCont = GetJetContainer(0);
+ if (!jet->GetNumberOfTracks())
+ return 0;
+ Double_t den=0.;
+ Double_t num = 0.;
+ AliVParticle *vp1 = 0x0;
+ for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
+ vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
+ Double_t dphi = vp1->Phi()-jet->Phi();
+ if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi();
+ if(dphi>TMath::Pi()) dphi-=TMath::TwoPi();
+ Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
+ Double_t dr = TMath::Sqrt(dr2);
+ num=num+vp1->Pt()*dr;
+ den=den+vp1->Pt();
+ }
+ return num/den;
+}
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::GetJetAngularity(AliEmcalJet *jet) {
+
+ if(fJetShapeType==kDeriv)
+ return jet->GetSecondOrderSubtractedAngularity();
+ else
+ return Angularity(jet);
+
+}
+
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::PTD(AliEmcalJet *jet){
+
+ AliJetContainer *jetCont = GetJetContainer(0);
+ if (!jet->GetNumberOfTracks())
+ return 0;
+ Double_t den=0.;
+ Double_t num = 0.;
+ AliVParticle *vp1 = 0x0;
+ for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
+ vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
+ num=num+vp1->Pt()*vp1->Pt();
+ den=den+vp1->Pt();
+ }
+ return TMath::Sqrt(num)/den;
+}
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::GetJetpTD(AliEmcalJet *jet) {
+ //calc subtracted jet mass
+ if(fJetShapeType==kDeriv)
+ return jet->GetSecondOrderSubtractedpTD();
+ else
+ return PTD(jet);
+
+}
+
+//_____________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::Circularity(AliEmcalJet *jet){
+
+ AliJetContainer *jetCont = GetJetContainer(0);
+ if (!jet->GetNumberOfTracks())
+ return 0;
+ Double_t mxx = 0.;
+ Double_t myy = 0.;
+ Double_t mxy = 0.;
+ int nc = 0;
+ Double_t sump2 = 0.;
+ Double_t pxjet=jet->Px();
+ Double_t pyjet=jet->Py();
+ Double_t pzjet=jet->Pz();
+
+
+ //2 general normalized vectors perpendicular to the jet
+ TVector3 ppJ1(pxjet, pyjet, pzjet);
+ TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
+ ppJ3.SetMag(1.);
+ TVector3 ppJ2(-pyjet, pxjet, 0);
+ ppJ2.SetMag(1.);
+ AliVParticle *vp1 = 0x0;
+ for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
+ vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
+
+
+ TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
+
+ //local frame
+ TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
+ TVector3 pPerp = pp - pLong;
+ //projection onto the two perpendicular vectors defined above
+
+ Float_t ppjX = pPerp.Dot(ppJ2);
+ Float_t ppjY = pPerp.Dot(ppJ3);
+ Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
+ if(ppjT<=0) return 0;
+
+ mxx += (ppjX * ppjX / ppjT);
+ myy += (ppjY * ppjY / ppjT);
+ mxy += (ppjX * ppjY / ppjT);
+ nc++;
+ sump2 += ppjT;}
+
+ if(nc<2) return 0;
+ if(sump2==0) return 0;
+ // Sphericity Matrix
+ Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
+ TMatrixDSym m0(2,ele);
+
+ // Find eigenvectors
+ TMatrixDSymEigen m(m0);
+ TVectorD eval(2);
+ TMatrixD evecm = m.GetEigenVectors();
+ eval = m.GetEigenValues();
+ // Largest eigenvector
+ int jev = 0;
+ // cout<<eval[0]<<" "<<eval[1]<<endl;
+ if (eval[0] < eval[1]) jev = 1;
+ TVectorD evec0(2);
+ // Principle axis
+ evec0 = TMatrixDColumn(evecm, jev);
+ Double_t compx=evec0[0];
+ Double_t compy=evec0[1];
+ TVector2 evec(compx, compy);
+ Double_t circ=0;
+ if(jev==1) circ=2*eval[0];
+ if(jev==0) circ=2*eval[1];
+
+ return circ;
+
+
+
+}
+
+
+
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::GetJetCircularity(AliEmcalJet *jet) {
+ //calc subtracted jet mass
+
+ if(fJetShapeType==kDeriv)
+ return jet->GetSecondOrderSubtractedCircularity();
+ else
+ return Circularity(jet);
+
+}
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::LeSub(AliEmcalJet *jet){
+
+ AliJetContainer *jetCont = GetJetContainer(0);
+ if (!jet->GetNumberOfTracks())
+ return 0;
+ Double_t den=0.;
+ Double_t num = 0.;
+ AliVParticle *vp1 = 0x0;
+ AliVParticle *vp2 = 0x0;
+ std::vector<int> ordindex;
+ ordindex=jet->SortConstituentsPt(jetCont->GetParticleContainer()->GetArray());
+
+ vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
+ vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
+
+ num=vp1->Pt();
+ den=vp2->Pt();
+
+return num-den;
+}
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::GetJetLeSub(AliEmcalJet *jet) {
+ //calc subtracted jet mass
+
+ if(fJetShapeType==kDeriv)
+ return jet->GetSecondOrderSubtractedLeSub();
+ else
+ return LeSub(jet);
+
+}
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::GetJetNumberOfConstituents(AliEmcalJet *jet) {
+ //calc subtracted jet mass
+
+ if(fJetShapeType==kDeriv)
+ return jet->GetSecondOrderSubtractedConstituent();
+ else
+ return jet->GetNumberOfTracks();
+
+}
+
+
+//______________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::Sigma2(AliEmcalJet *jet){
+
+ AliJetContainer *jetCont = GetJetContainer(0);
+ if (!jet->GetNumberOfTracks())
+ return 0;
+ Double_t mxx = 0.;
+ Double_t myy = 0.;
+ Double_t mxy = 0.;
+ int nc = 0;
+ Double_t sump2 = 0.;
+
+ AliVParticle *vp1 = 0x0;
+ for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
+ vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
+ Double_t ppt=vp1->Pt();
+ Double_t dphi = vp1->Phi()-jet->Phi();
+ if(dphi<-1.*TMath::Pi()) dphi+=TMath::TwoPi();
+ if(dphi>TMath::Pi()) dphi-=TMath::TwoPi();
+ Double_t deta = vp1->Eta()-jet->Eta();
+ mxx += ppt*ppt*deta*deta;
+ myy += ppt*ppt*dphi*dphi;
+ mxy -= ppt*ppt*deta*dphi;
+ nc++;
+ sump2 += ppt*ppt;
+
+ }
+ if(nc<2) return 0;
+ if(sump2==0) return 0;
+ // Sphericity Matrix
+ Double_t ele[4] = {mxx , mxy , mxy , myy };
+ TMatrixDSym m0(2,ele);
+
+ // Find eigenvectors
+ TMatrixDSymEigen m(m0);
+ TVectorD eval(2);
+ TMatrixD evecm = m.GetEigenVectors();
+ eval = m.GetEigenValues();
+ // Largest eigenvector
+ int jev = 0;
+ // cout<<eval[0]<<" "<<eval[1]<<endl;
+ if (eval[0] < eval[1]) jev = 1;
+ TVectorD evec0(2);
+ // Principle axis
+ evec0 = TMatrixDColumn(evecm, jev);
+ Double_t compx=evec0[0];
+ Double_t compy=evec0[1];
+ TVector2 evec(compx, compy);
+ Double_t sig=0;
+ if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
+ if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
+
+ return sig;
+
+}
+
+//________________________________________________________________________
+Float_t AliAnalysisTaskEmcalQGTagging::GetSigma2(AliEmcalJet *jet) {
+ //calc subtracted jet mass
+
+ if(fJetShapeType==kDeriv)
+ return jet->GetSecondOrderSubtractedSigma2();
+ else
+ return Sigma2(jet);
+
+}
+
+//________________________________________________________________________
+Bool_t AliAnalysisTaskEmcalQGTagging::RetrieveEventObjects() {
+ //
+ // retrieve event objects
+ //
+ if (!AliAnalysisTaskEmcalJet::RetrieveEventObjects())
+ return kFALSE;
+
+ return kTRUE;
+}
+
+//_______________________________________________________________________
+void AliAnalysisTaskEmcalQGTagging::Terminate(Option_t *)
+{
+ // Called once at the end of the analysis.
+
+ // fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(1));
+ // if (!fTreeObservableTagging){
+ // Printf("ERROR: fTreeObservableTagging not available");
+ // return;
+ // }
+
+}
+
--- /dev/null
+#ifndef ALIANALYSISTASKEMCALQGTAGGING_H
+#define ALIANALYSISTASKEMCALQGTAGGING_H
+
+class TH1;
+class TH2;
+class TH3;
+class TH3F;
+class TTree;
+class THnSparse;
+class TClonesArray;
+class TArrayI;
+class AliAnalysisManager;
+class AliJetContainer;
+
+#include "AliAnalysisTaskEmcalJet.h"
+
+
+
+class AliAnalysisTaskEmcalQGTagging : public AliAnalysisTaskEmcalJet {
+ public:
+ enum JetShapeType {
+ kRaw = 0, //mass form anti-kt 4-vector
+ kConstSub = 1, //constituent subtracted jetshape
+ kTrue = 2,
+ kDeriv = 3 //area based subtracted jet mass
+ };
+
+ AliAnalysisTaskEmcalQGTagging();
+ AliAnalysisTaskEmcalQGTagging(const char *name);
+ virtual ~AliAnalysisTaskEmcalQGTagging();
+
+ void UserCreateOutputObjects();
+ void Terminate(Option_t *option);
+
+ //Setters
+ void SetJetContainer(Int_t c) { fContainer = c ; }
+ void SetMinFractionShared(Double_t f) { fMinFractionShared = f ; }
+ void SetJetShapeType(JetShapeType t) { fJetShapeType = t ; }
+ void SetMCTask(Int_t f) { fIsMC = f ; }
+ void SetEmbeddingTask(Int_t f) { fIsEmbedding = f ; }
+ void SetIsConstSub(Bool_t f) { fIsConstSub = f ; }
+ void SetRMatching(Float_t f) { fRMatching = f ;}
+
+ protected:
+ Bool_t RetrieveEventObjects();
+ Bool_t Run();
+ Bool_t FillHistograms();
+
+ Float_t GetJetMass(AliEmcalJet *jet);
+ Float_t Angularity(AliEmcalJet *jet);
+ Float_t GetJetAngularity(AliEmcalJet *jet);
+ Float_t PTD(AliEmcalJet *jet);
+ Float_t GetJetpTD(AliEmcalJet *jet);
+ Float_t Circularity(AliEmcalJet *jet);
+ Float_t GetJetCircularity(AliEmcalJet *jet);
+ Float_t LeSub(AliEmcalJet *jet);
+ Float_t GetJetLeSub(AliEmcalJet *jet);
+ Float_t GetJetNumberOfConstituents(AliEmcalJet *jet);
+ Float_t GetSigma2(AliEmcalJet *jet);
+ Float_t Sigma2(AliEmcalJet *jet);
+
+
+ Int_t fContainer; // jets to be analyzed 0 for Base, 1 for subtracted.
+ Float_t fMinFractionShared; // only fill histos for jets if shared fraction larger than X
+ JetShapeType fJetShapeType; // jet mass type to be used
+ Float_t *fShapesVar; // jet shapes used for the tagging
+ Int_t fIsMC;
+ Int_t fIsEmbedding;
+ Bool_t fIsConstSub;
+ Float_t fRMatching;
+
+ TH2F *fPhiJetCorr6;
+ TH2F *fPhiJetCorr7;
+ TH2F *fEtaJetCorr6;
+ TH2F *fEtaJetCorr7;
+ TH2F *fPtJetCorr;
+ TH1F *fPtJet;
+
+
+ TTree *fTreeObservableTagging; //Tree with tagging variables subtracted MC or true MC or raw
+
+ private:
+ AliAnalysisTaskEmcalQGTagging(const AliAnalysisTaskEmcalQGTagging&); // not implemented
+ AliAnalysisTaskEmcalQGTagging &operator=(const AliAnalysisTaskEmcalQGTagging&); // not implemented
+
+ ClassDef(AliAnalysisTaskEmcalQGTagging, 1)
+};
+#endif
+
--- /dev/null
+AliAnalysisTaskEmcalQGTagging* AddTaskEmcalQGTagging(const char * njetsBase,
+ const char * njetsTrue,
+ const Double_t R,
+ const char * nrhoBase,
+ const char * ntracks,
+ const char * nclusters,
+ const char * ntracksTrue,
+ const char *type,
+ const char *CentEst,
+ Int_t pSel,
+ TString trigClass = "",
+ TString kEmcalTriggers = "",
+ TString tag = "",
+ AliAnalysisTaskEmcalQGTagging::JetShapeType jetShapeType, Int_t isEmbedding = 0) {
+
+ AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+ if (!mgr)
+ {
+ Error("AddTaskEmcalQGTagging","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("AddTaskEmcalQGTagging", "This task requires an input event handler");
+ return NULL;
+ }
+
+ TString wagonName = Form("JetQGTaggings_%s_TC%s%s",njetsBase,trigClass.Data(),tag.Data());
+
+ //Configure jet tagger task
+ AliAnalysisTaskEmcalQGTagging *task = new AliAnalysisTaskEmcalQGTagging(wagonName.Data());
+
+ //task->SetNCentBins(4);
+ task->SetJetShapeType(jetShapeType);
+ TString thename(njetsBase);
+ if(thename.Contains("Sub")) task->SetIsConstSub(kTRUE);
+ //task->SetVzRange(-10.,10.);
+
+ AliParticleContainer *trackCont = task->AddParticleContainer(ntracks);
+ AliParticleContainer *trackContTrue = task->AddParticleContainer(ntracksTrue);
+ AliClusterContainer *clusterCont = task->AddClusterContainer(nclusters);
+
+ task->SetJetContainer(0);
+
+ TString strType(type);
+ AliJetContainer *jetContBase = task->AddJetContainer(njetsBase,strType,R);
+ if(jetContBase) {
+ jetContBase->SetRhoName(nrhoBase);
+ jetContBase->ConnectParticleContainer(trackCont);
+ jetContBase->ConnectClusterContainer(clusterCont);
+ jetContBase->SetPercAreaCut(0.6);
+ if(jetShapeType == AliAnalysisTaskEmcalQGTagging::kTrue){ jetContBase->SetPartonInfoName("PartonsInfo");}
+ }
+
+
+ if(isEmbedding){
+ task->SetJetContainer(1);
+
+ AliJetContainer *jetContTrue = task->AddJetContainer(njetsTrue,strType,R);
+ if(jetContTrue) {
+ jetContTrue->SetRhoName(nrhoBase);
+ jetContTrue->ConnectParticleContainer(trackContTrue);
+ jetContTrue->SetPercAreaCut(0.6);
+ jetContTrue->SetPartonInfoName("PartonsInfo");
+ }
+ }
+ // task->SetJetContainer(1);
+
+ 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 contName1(wagonName);
+
+ TString outputfile = Form("%s",AliAnalysisManager::GetCommonFileName());
+ AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(contName1.Data(), TTree::Class(),AliAnalysisManager::kOutputContainer,outputfile);
+
+ mgr->ConnectOutput(task,1,coutput1);
+
+ return task;
+
+}
+
AliJetEmbeddingFromGenTask* AddTaskJetEmbeddingFromGen(
AliGenerator *genGen,
const char *tracksName = "GenParticles",
+ const char *partonInfoName = "PartonsInfo",
const char *taskName = "JetEmbeddingFromGenTask",
const Double_t minPt = 10,
const Double_t maxPt = 10,
if(TString(genGen->IsA()->GetName()).EqualTo("AliGenPythia")) genGen->AliGenPythia::SetEventListRange(0, 0);
jetEmb->SetGen(genGen);
jetEmb->SetTracksName(tracksName);
+ jetEmb->SetPartonInfoName(partonInfoName);
jetEmb->SetEtaRange(minEta, maxEta);
jetEmb->SetPhiRange(minPhi, maxPhi);
jetEmb->SetPtRange(minPt, maxPt);
#pragma link C++ class AliRhoParameter+;
#pragma link C++ class AliLocalRhoParameter+;
#pragma link C++ class AliJetTriggerSelectionTask+;
+#pragma link C++ class AliStackPartonInfo+;
// user tasks
#pragma link C++ class AliAnalysisTaskCLQA+;
#pragma link C++ class AliAnalysisTaskRhoMass+;
#pragma link C++ class AliAnalysisTaskRhoMassBase+;
#pragma link C++ class AliAnalysisTaskRhoMassScale+;
+#pragma link C++ class AliAnalysisTaskEmcalQGTagging+;
#pragma link C++ class AliAnalysisTaskSAJF+;
#pragma link C++ class AliAnalysisTaskSAQA+;
#pragma link C++ class AliAnalysisTaskSOH+;
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff