From b16bb00196553a8ad7361d75b928d5acd27f7e4a Mon Sep 17 00:00:00 2001 From: loizides Date: Sat, 26 Jan 2013 23:11:01 +0000 Subject: [PATCH] Update from Salvatore for full jet embedding --- PWGJE/CMakelibPWGJEEMCALJetTasks.pkg | 1 + .../AliJetEmbeddingFromAODTask.cxx | 380 +++++++++++++++ .../AliJetEmbeddingFromAODTask.h | 75 +++ PWGJE/EMCALJetTasks/AliJetEmbeddingTask.cxx | 13 - PWGJE/EMCALJetTasks/AliJetEmbeddingTask.h | 3 - PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx | 443 ++++++++++++------ PWGJE/EMCALJetTasks/AliJetModelBaseTask.h | 43 +- PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx | 247 +++++----- PWGJE/EMCALJetTasks/AliJetResponseMaker.h | 88 ++-- .../UserTasks/AliAnalysisTaskSAQA.cxx | 13 +- .../macros/AddTaskJetEmbeddingFromAOD.C | 140 ++++++ .../macros/AddTaskJetRespPtHard.C | 3 +- .../macros/AddTaskJetResponseMaker.C | 2 + PWGJE/PWGJEEMCALJetTasksLinkDef.h | 1 + 14 files changed, 1132 insertions(+), 320 deletions(-) create mode 100644 PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.cxx create mode 100644 PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.h create mode 100644 PWGJE/EMCALJetTasks/macros/AddTaskJetEmbeddingFromAOD.C diff --git a/PWGJE/CMakelibPWGJEEMCALJetTasks.pkg b/PWGJE/CMakelibPWGJEEMCALJetTasks.pkg index bc4516d695b..3da6d6c5012 100644 --- a/PWGJE/CMakelibPWGJEEMCALJetTasks.pkg +++ b/PWGJE/CMakelibPWGJEEMCALJetTasks.pkg @@ -39,6 +39,7 @@ set ( SRCS EMCALJetTasks/AliHadCorrTask.cxx EMCALJetTasks/AliJetEmbeddingTask.cxx EMCALJetTasks/AliJetEmbeddingFromGenTask.cxx + EMCALJetTasks/AliJetEmbeddingFromAODTask.cxx EMCALJetTasks/AliJetModelBaseTask.cxx EMCALJetTasks/AliJetRandomizerTask.cxx EMCALJetTasks/AliJetResponseMaker.cxx diff --git a/PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.cxx b/PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.cxx new file mode 100644 index 00000000000..c1c86bab429 --- /dev/null +++ b/PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.cxx @@ -0,0 +1,380 @@ +// $Id: AliJetEmbeddingFromAODTask.cxx $ +// +// Jet embedding from AOD task. +// +// Author: S.Aiola, C.Loizides + +#include "AliJetEmbeddingFromAODTask.h" + +//#include + +#include +#include +#include +#include +#include +#include +#include + +#include "AliVEvent.h" +#include "AliAODTrack.h" +#include "AliESDtrack.h" +#include "AliPicoTrack.h" +#include "AliVTrack.h" +#include "AliVCluster.h" +#include "AliVCaloCells.h" +#include "AliCentrality.h" +#include "AliVHeader.h" +#include "AliVVertex.h" +#include "AliAODHeader.h" +#include "AliLog.h" +#include "AliInputEventHandler.h" + +ClassImp(AliJetEmbeddingFromAODTask) + +//________________________________________________________________________ +AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask() : + AliJetModelBaseTask("AliJetEmbeddingFromAODTask"), + fFileList(0), + fAODTreeName(), + fAODHeaderName(), + fAODVertexName(), + fAODTrackName(), + fAODClusName(), + fAODCellsName(), + fMinCentrality(0), + fMaxCentrality(10), + fTriggerMask(AliVEvent::kAny), + fZVertexCut(10), + fIncludeNoITS(kTRUE), + fTotalFiles(2200), + fEsdTreeMode(kFALSE), + fCurrentFileID(0), + fCurrentAODFileID(-1), + fCurrentAODFile(0), + fAODHeader(0), + fAODVertex(0), + fAODTracks(0), + fAODClusters(0), + fAODCaloCells(0), + fHistFileIDs(0) +{ + // Default constructor. + SetSuffix("AODEmbedding"); + fMarkMC = kFALSE; + fAODfilterBits[0] = -1; + fAODfilterBits[1] = -1; +} + +//________________________________________________________________________ +AliJetEmbeddingFromAODTask::AliJetEmbeddingFromAODTask(const char *name, Bool_t drawqa) : + AliJetModelBaseTask(name, drawqa), + fFileList(0), + fAODTreeName("aodTree"), + fAODHeaderName("header"), + fAODVertexName("vertices"), + fAODTrackName("tracks"), + fAODClusName("caloClusters"), + fAODCellsName("emcalCells"), + fMinCentrality(0), + fMaxCentrality(10), + fTriggerMask(AliVEvent::kAny), + fZVertexCut(10), + fIncludeNoITS(kTRUE), + fTotalFiles(2200), + fEsdTreeMode(kFALSE), + fCurrentFileID(0), + fCurrentAODFileID(-1), + fCurrentAODFile(0), + fAODHeader(0), + fAODVertex(0), + fAODTracks(0), + fAODClusters(0), + fAODCaloCells(0), + fHistFileIDs(0) +{ + // Standard constructor. + SetSuffix("AODEmbedding"); + fMarkMC = kFALSE; + fAODfilterBits[0] = -1; + fAODfilterBits[1] = -1; +} + +//________________________________________________________________________ +AliJetEmbeddingFromAODTask::~AliJetEmbeddingFromAODTask() +{ + // Destructor + + if (fCurrentAODFile) { + fCurrentAODFile->Close(); + delete fCurrentAODFile; + } +} + +//________________________________________________________________________ +void AliJetEmbeddingFromAODTask::UserCreateOutputObjects() +{ + if (!fQAhistos) + return; + + AliJetModelBaseTask::UserCreateOutputObjects(); + + fHistFileIDs = new TH2C("fHistFileIDs", "fHistFileIDs", fTotalFiles, -0.5, fTotalFiles - 0.5, fFileList->GetEntriesFast(), -0.5, fFileList->GetEntriesFast() -0.5); + fHistFileIDs->GetXaxis()->SetTitle("File no. (PYTHIA)"); + fHistFileIDs->GetYaxis()->SetTitle("File no. (Embedded AOD)"); + fOutput->Add(fHistFileIDs); + + PostData(1, fOutput); +} + +//________________________________________________________________________ +Bool_t AliJetEmbeddingFromAODTask::UserNotify() +{ + TString path(fInputHandler->GetTree()->GetCurrentFile()->GetPath()); + if (path.EndsWith("/")) + path.Remove(path.Length()-1); + path.Remove(path.Last('/')); + path.Remove(0, path.Last('/')+1); + if (!path.IsDec()) { + AliWarning(Form("Could not extract file number from path %s", path.Data())); + return kTRUE; + } + fCurrentFileID = path.Atoi(); + fCurrentAODFileID = fFileList->GetEntriesFast() * fCurrentFileID / fTotalFiles; + AliInfo(Form("Start embedding from file ID %d", fCurrentAODFileID)); + return kTRUE; +} + +//________________________________________________________________________ +Bool_t AliJetEmbeddingFromAODTask::ExecOnce() +{ + if (fAODTreeName.Contains("aod")) + fEsdTreeMode = kFALSE; + else + fEsdTreeMode = kTRUE; + + return AliJetModelBaseTask::ExecOnce(); +} + +//________________________________________________________________________ +Bool_t AliJetEmbeddingFromAODTask::OpenNextFile() +{ + if (fCurrentAODFile) { + fCurrentAODFile->Close(); + delete fCurrentAODFile; + fCurrentAODFile = 0; + } + + if (fCurrentAODFileID >= fFileList->GetEntriesFast()) + return kFALSE; + + TObjString *objFileName = static_cast(fFileList->At(fCurrentAODFileID)); + TString fileName = objFileName->GetString(); + + if (fileName.BeginsWith("alien://") && !gGrid) { + AliInfo("Trying to connect to AliEn ..."); + TGrid::Connect("alien://"); + } + + fCurrentAODFile = TFile::Open(fileName); + + if (!fCurrentAODFile || fCurrentAODFile->IsZombie()) { + return kFALSE; + } + + if (fQAhistos) + fHistFileIDs->Fill(fCurrentFileID, fCurrentAODFileID); + + fCurrentAODFileID++; + return kTRUE; +} + +//________________________________________________________________________ +Bool_t AliJetEmbeddingFromAODTask::GetNextEntry() +{ + static TTree *tree = 0; + static Int_t entry = 0; + + Int_t attempts = 0; + + do { + + if (!fCurrentAODFile || !tree || entry >= tree->GetEntries()) { + if (!OpenNextFile()) + return kFALSE; + + tree = static_cast(fCurrentAODFile->Get(fAODTreeName)); + if (!tree) + return kFALSE; + + if (!fAODHeaderName.IsNull()) + tree->SetBranchAddress(fAODHeaderName, &fAODHeader); + + if (!fAODVertexName.IsNull()) + tree->SetBranchAddress(fAODVertexName, &fAODVertex); + + if (!fAODTrackName.IsNull()) + tree->SetBranchAddress(fAODTrackName, &fAODTracks); + + if (!fAODClusName.IsNull()) + tree->SetBranchAddress(fAODClusName, &fAODClusters); + + if (!fAODCellsName.IsNull()) + tree->SetBranchAddress(fAODCellsName, &fAODCaloCells); + + entry = 0; + } + + tree->GetEntry(entry); + entry++; + attempts++; + + if (attempts == 1000) + AliWarning("After 1000 attempts no event has been accepted by the event selection (trigger, centrality...)!"); + + } while (!IsAODEventSelected() && tree); + + return (tree!=0); +} + +//________________________________________________________________________ +Bool_t AliJetEmbeddingFromAODTask::IsAODEventSelected() +{ + if (!fEsdTreeMode && fAODHeader) { + AliAODHeader *aodHeader = static_cast(fAODHeader); + + UInt_t offlineTrigger = aodHeader->GetOfflineTrigger(); + + if ((offlineTrigger & fTriggerMask) == 0) { + AliDebug(2, Form("Event rejected due to physics selection. Event trigger mask: %d, trigger mask selection: %d.", offlineTrigger, fTriggerMask)); + return kFALSE; + } + + AliCentrality *cent = aodHeader->GetCentralityP(); + Float_t centVal = cent->GetCentralityPercentile("V0M"); + if (centVal < fMinCentrality || centVal >= fMaxCentrality) { + AliDebug(2, Form("Event rejected due to centrality selection. Event centrality: %f, centrality range selection: %f to %f", centVal, fMinCentrality, fMaxCentrality)); + return kFALSE; + } + } + + if (fAODVertex) { + Double_t vert[3]={0}; + ((AliVVertex*)fAODVertex->At(0))->GetXYZ(vert); + if (TMath::Abs(vert[2]) > fZVertexCut) { + AliDebug(2,Form("Event rejected due to Z vertex selection. Event Z vertex: %f, Z vertex cut: %f", vert[2], fZVertexCut)); + return kFALSE; + } + } + + return kTRUE; +} + +//________________________________________________________________________ +void AliJetEmbeddingFromAODTask::Run() +{ + if (!GetNextEntry()) { + AliError("Unable to get AOD event to embed. Nothing will be embedded."); + return; + } + + if (fTracks && fOutTracks) { + + if (fCopyArray) + CopyTracks(); + + if (fAODTracks) { + AliDebug(2, Form("%d tracks will be processed for embedding.", fAODTracks->GetEntriesFast())); + for (Int_t i = 0; i < fAODTracks->GetEntriesFast(); i++) { + AliVTrack *track = static_cast(fAODTracks->At(i)); + if (!track) { + AliError(Form("Could not find track %d in branch %s of tree %s!", i, fAODTrackName.Data(), fAODTreeName.Data())); + continue; + } + + Int_t type = 0; + Bool_t isEmc = kFALSE; + if (!fEsdTreeMode) { + AliAODTrack *aodtrack = static_cast(track); + if (aodtrack->TestFilterBit(fAODfilterBits[0])) { + type = 0; + } + else if (aodtrack->TestFilterBit(fAODfilterBits[1])) { + if ((aodtrack->GetStatus()&AliESDtrack::kITSrefit)==0) { + if (fIncludeNoITS) + type = 2; + else + continue; + } + else { + type = 1; + } + } + else { /*not a good track*/ + continue; + } + + if (TMath::Abs(aodtrack->GetTrackEtaOnEMCal()) < 0.75 && + aodtrack->GetTrackPhiOnEMCal() > 70 * TMath::DegToRad() && + aodtrack->GetTrackPhiOnEMCal() < 190 * TMath::DegToRad()) + isEmc = kTRUE; + } + else { /*not AOD mode, let's see if it is PicoTrack*/ + AliPicoTrack *ptrack = dynamic_cast(track); + if (ptrack) { + type = ptrack->GetTrackType(); + isEmc = ptrack->IsEMCAL(); + } + } + + AliDebug(3, Form("Embedding track with pT = %f, eta = %f, phi = %f", track->Pt(), track->Eta(), track->Phi())); + AddTrack(track->Pt(), track->Eta(), track->Phi(), type, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), isEmc); + } + } + } + + if (fClusters && fOutClusters) { + + if (fCopyArray) + CopyClusters(); + + if (fAODClusters) { + for (Int_t i = 0; i < fAODClusters->GetEntriesFast(); i++) { + AliVCluster *clus = static_cast(fAODClusters->At(i)); + if (!clus) { + AliError(Form("Could not find cluster %d in branch %s of tree %s!", i, fAODClusName.Data(), fAODTreeName.Data())); + continue; + } + TLorentzVector vect; + Double_t vert[3] = {0,0,0}; + clus->GetMomentum(vect,vert); + AddCluster(clus->E(), vect.Eta(), vect.Phi()); + } + } + } + + if (fCaloCells && fOutCaloCells) { + + Int_t totalCells = fCaloCells->GetNumberOfCells(); + if (fAODCaloCells) + totalCells += fAODCaloCells->GetNumberOfCells(); + + SetNumberOfOutCells(totalCells); + + CopyCells(); + + if (fAODCaloCells) { + for (Short_t i = 0; i < fAODCaloCells->GetNumberOfCells(); i++) { + Short_t mclabel = -1; + Double_t efrac = 0.; + Double_t time = -1; + Short_t cellNum = -1; + Double_t amp = -1; + + fAODCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac); + AliDebug(3,Form("Adding cell with amplitude %f, absolute ID %d, time %f", amp, cellNum, time)); + AddCell(amp, cellNum, time); + } + } + } +} diff --git a/PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.h b/PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.h new file mode 100644 index 00000000000..b29db549d28 --- /dev/null +++ b/PWGJE/EMCALJetTasks/AliJetEmbeddingFromAODTask.h @@ -0,0 +1,75 @@ +#ifndef ALIJETEMBEDDINGFROMAODTASK_H +#define ALIJETEMBEDDINGFROMAODTASK_H + +// $Id: AliJetEmbeddingFromAODTask.h $ + +class TFile; +class TObjArray; +class TClonesArray; +class TString; +class AliVCaloCells; +class AliVHeader; +class TH2; + +#include "AliJetModelBaseTask.h" + +class AliJetEmbeddingFromAODTask : public AliJetModelBaseTask { + public: + AliJetEmbeddingFromAODTask(); + AliJetEmbeddingFromAODTask(const char *name, Bool_t drawqa=kFALSE); + virtual ~AliJetEmbeddingFromAODTask(); + + void UserCreateOutputObjects(); + Bool_t UserNotify(); + + void SetFileList(TObjArray *list) { fFileList = list ; } + void SetAODTreeName(const char *t) { fAODTreeName = t ; } + void SetAODHeaderName(const char *t) { fAODHeaderName = t ; } + void SetAODTracksName(const char *n) { fAODTrackName = n ; } + void SetAODClusName(const char *n) { fAODClusName = n ; } + void SetAODCellsName(const char *n) { fAODCellsName = n ; } + void SetCentralityRange(Double_t min, Double_t max) { fMinCentrality = min ; fMaxCentrality = max ; } + void SetTriggerMask(UInt_t mask) { fTriggerMask = mask ; } + void SetAODfilterBits(Int_t b0 = 0, Int_t b1 = 0) { fAODfilterBits[0] = b0 ; fAODfilterBits[1] = b1 ; } + void SetIncludeNoITS(Bool_t f) { fIncludeNoITS = f ; } + void SetTotalFiles(Int_t n) { fTotalFiles = n ; } + + protected: + Bool_t ExecOnce() ;// intialize task + void Run() ;// do jet model action + Bool_t OpenNextFile() ;// open next file in fFileList + Bool_t GetNextEntry() ;// get next entry in current tree + Bool_t IsAODEventSelected() ;// AOD event trigger/centrality selection + + TObjArray *fFileList ;// List of AOD files + TString fAODTreeName ;// Name of the tree in the AOD file + TString fAODHeaderName ;// Name of the header in the AOD tree + TString fAODVertexName ;// Name of the vertex branch in the AOD tree + TString fAODTrackName ;// Name of the track collection branch in the AOD tree + TString fAODClusName ;// Name of the cluster collection branch in the AOD tree + TString fAODCellsName ;// Name of the cell collection branch in the AOD tree + Double_t fMinCentrality ;// Minimum centrality + Double_t fMaxCentrality ;// Maximum centrality + UInt_t fTriggerMask ;// Trigger selection mask + Double_t fZVertexCut ;// Z vertex cut + Int_t fAODfilterBits[2] ;// AOD track filter bit map + Bool_t fIncludeNoITS ;// True = includes tracks with failed ITS refit + Int_t fTotalFiles ;// Total number of files per pt hard bin + Bool_t fEsdTreeMode ;//! True = embed from ESD (must be a skimmed ESD!) + Int_t fCurrentFileID ;//! Current file being processed (trough the event handler) + Int_t fCurrentAODFileID ;//! Current file ID + TFile *fCurrentAODFile ;//! Current open file + AliVHeader *fAODHeader ;//! AOD header + TClonesArray *fAODVertex ;//! AOD vertex + TClonesArray *fAODTracks ;//! AOD track collection + TClonesArray *fAODClusters ;//! AOD cluster collection + AliVCaloCells *fAODCaloCells ;//! AOD cell collection + TH2 *fHistFileIDs ;//! Current file ID vs. AOD file ID (to be embedded) + + private: + AliJetEmbeddingFromAODTask(const AliJetEmbeddingFromAODTask&); // not implemented + AliJetEmbeddingFromAODTask &operator=(const AliJetEmbeddingFromAODTask&); // not implemented + + ClassDef(AliJetEmbeddingFromAODTask, 1) // Jet embedding from AOD task +}; +#endif diff --git a/PWGJE/EMCALJetTasks/AliJetEmbeddingTask.cxx b/PWGJE/EMCALJetTasks/AliJetEmbeddingTask.cxx index 5fa78e2d365..59dc7054eab 100644 --- a/PWGJE/EMCALJetTasks/AliJetEmbeddingTask.cxx +++ b/PWGJE/EMCALJetTasks/AliJetEmbeddingTask.cxx @@ -6,19 +6,6 @@ #include "AliJetEmbeddingTask.h" -#include -#include -#include - -#include "AliAnalysisManager.h" -#include "AliVEvent.h" -#include "AliVCluster.h" -#include "AliEMCALDigit.h" -#include "AliEMCALRecPoint.h" -#include "AliPicoTrack.h" -#include "AliEMCALGeometry.h" -#include "AliLog.h" - ClassImp(AliJetEmbeddingTask) //________________________________________________________________________ diff --git a/PWGJE/EMCALJetTasks/AliJetEmbeddingTask.h b/PWGJE/EMCALJetTasks/AliJetEmbeddingTask.h index ab4f08d8893..b065d620343 100644 --- a/PWGJE/EMCALJetTasks/AliJetEmbeddingTask.h +++ b/PWGJE/EMCALJetTasks/AliJetEmbeddingTask.h @@ -3,9 +3,6 @@ // $Id$ -class TClonesArray; -class AliEMCALGeometry; - #include "AliJetModelBaseTask.h" class AliJetEmbeddingTask : public AliJetModelBaseTask { diff --git a/PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx b/PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx index af91d1becbb..d8b28978d59 100644 --- a/PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx +++ b/PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx @@ -10,17 +10,20 @@ #include #include #include +#include +#include "AliVEvent.h" #include "AliAODCaloCluster.h" -#include "AliAnalysisManager.h" +#include "AliESDCaloCluster.h" +#include "AliVCluster.h" #include "AliEMCALDigit.h" -#include "AliEMCALGeometry.h" #include "AliEMCALRecPoint.h" -#include "AliESDCaloCluster.h" -#include "AliLog.h" +#include "AliESDCaloCells.h" +#include "AliAODCaloCells.h" +#include "AliVCaloCells.h" #include "AliPicoTrack.h" -#include "AliVCluster.h" -#include "AliVEvent.h" +#include "AliEMCALGeometry.h" +#include "AliLog.h" ClassImp(AliJetModelBaseTask) @@ -32,6 +35,8 @@ AliJetModelBaseTask::AliJetModelBaseTask() : fOutTracksName(), fCaloName(), fOutCaloName(), + fCellsName(), + fOutCellsName(), fSuffix(), fEtaMin(-1), fEtaMax(1), @@ -41,27 +46,36 @@ AliJetModelBaseTask::AliJetModelBaseTask() : fPtMax(0), fCopyArray(kTRUE), fNClusters(0), + fNCells(0), fNTracks(0), fMarkMC(kTRUE), fPtSpectrum(0), + fQAhistos(kFALSE), fIsInit(0), fGeom(0), fClusters(0), fOutClusters(0), fTracks(0), - fOutTracks(0) + fOutTracks(0), + fCaloCells(0), + fOutCaloCells(0), + fAddedCells(0), + fEsdMode(kFALSE), + fOutput(0) { // Default constructor. } //________________________________________________________________________ -AliJetModelBaseTask::AliJetModelBaseTask(const char *name) : +AliJetModelBaseTask::AliJetModelBaseTask(const char *name, Bool_t drawqa) : AliAnalysisTaskSE(name), fGeomName(""), fTracksName("PicoTracks"), fOutTracksName("PicoTracksEmbedded"), fCaloName("CaloClustersCorr"), fOutCaloName("CaloClustersCorrEmbedded"), + fCellsName(""), + fOutCellsName(""), fSuffix("Processed"), fEtaMin(-1), fEtaMax(1), @@ -71,17 +85,28 @@ AliJetModelBaseTask::AliJetModelBaseTask(const char *name) : fPtMax(60), fCopyArray(kTRUE), fNClusters(0), + fNCells(0), fNTracks(1), fMarkMC(kTRUE), fPtSpectrum(0), + fQAhistos(drawqa), fIsInit(0), fGeom(0), fClusters(0), fOutClusters(0), fTracks(0), - fOutTracks(0) + fOutTracks(0), + fCaloCells(0), + fOutCaloCells(0), + fAddedCells(0), + fEsdMode(kFALSE), + fOutput(0) { // Standard constructor. + + if (fQAhistos) { + DefineOutput(1, TList::Class()); + } } //________________________________________________________________________ @@ -90,6 +115,20 @@ AliJetModelBaseTask::~AliJetModelBaseTask() // Destructor } +//________________________________________________________________________ +void AliJetModelBaseTask::UserCreateOutputObjects() +{ + // Create user output. + if (!fQAhistos) + return; + + OpenFile(1); + fOutput = new TList(); + fOutput->SetOwner(); + + PostData(1, fOutput); +} + //________________________________________________________________________ void AliJetModelBaseTask::UserExec(Option_t *) { @@ -108,9 +147,256 @@ void AliJetModelBaseTask::UserExec(Option_t *) fOutClusters->Delete(); } + if (fCaloCells) { + fAddedCells = 0; + if (!fCopyArray) + fCaloCells = static_cast(fCaloCells->Clone(Form("%s_old",fCaloCells->GetName()))); + } + Run(); } +//________________________________________________________________________ +Bool_t AliJetModelBaseTask::ExecOnce() +{ + // Init task. + + delete gRandom; + gRandom = new TRandom3(0); + + fEsdMode = InputEvent()->InheritsFrom("AliESDEvent"); + + if (fPtMax < fPtMin) { + AliWarning (Form("PtMax (%f) < PtMin (%f), setting PtMax = PtMin = %f", fPtMax, fPtMin, fPtMin)); + fPtMax = fPtMin; + } + + if (fEtaMax < fEtaMin) { + AliWarning (Form("EtaMax (%f) < EtaMin (%f), setting EtaMax = EtaMin = %f", fEtaMax, fEtaMin, fEtaMin)); + fEtaMax = fEtaMin; + } + + if (fPhiMax < fPhiMin) { + AliWarning (Form("PhiMax (%f) < PhiMin (%f), setting PhiMax = PhiMin = %f", fPhiMax, fPhiMin, fPhiMin)); + fPhiMax = fPhiMin; + } + + if (!fCaloCells && !fCellsName.IsNull()) { + fCaloCells = dynamic_cast(InputEvent()->FindListObject(fCellsName)); + if (!fCaloCells) { + AliError(Form("%s: Couldn't retrieve calo cells with name %s!", GetName(), fCellsName.Data())); + return kFALSE; + } + + if (!fCaloCells->InheritsFrom("AliVCaloCells")) { + AliError(Form("%s: Collection %s does not contain a AliVCaloCells object!", GetName(), fCellsName.Data())); + fCaloCells = 0; + return kFALSE; + } + + if (!fOutCaloCells) { + fOutCellsName = fCellsName; + if (fCopyArray) { + fOutCellsName += fSuffix; + + if (fEsdMode) + fOutCaloCells = new AliESDCaloCells(fOutCellsName,fOutCellsName); + else + fOutCaloCells = new AliAODCaloCells(fOutCellsName,fOutCellsName); + + if (InputEvent()->FindListObject(fOutCellsName)) { + AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCellsName.Data())); + return kFALSE; + } + else { + InputEvent()->AddObject(fOutCaloCells); + } + } + else { + fOutCaloCells = fCaloCells; + } + } + } + + if (fNTracks > 0 && !fTracks && !fTracksName.IsNull()) { + fTracks = dynamic_cast(InputEvent()->FindListObject(fTracksName)); + if (!fTracks) { + AliError(Form("%s: Couldn't retrieve tracks with name %s!", GetName(), fTracksName.Data())); + return kFALSE; + } + + if (!fTracks->GetClass()->GetBaseClass("AliPicoTrack")) { + AliError(Form("%s: Collection %s does not contain AliPicoTrack objects!", GetName(), fTracksName.Data())); + fTracks = 0; + return kFALSE; + } + + if (!fOutTracks) { + fOutTracksName = fTracksName; + if (fCopyArray) { + fOutTracksName += fSuffix; + fOutTracks = new TClonesArray("AliPicoTrack", fTracks->GetSize()); + fOutTracks->SetName(fOutTracksName); + if (InputEvent()->FindListObject(fOutTracksName)) { + AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutTracksName.Data())); + return kFALSE; + } + else { + InputEvent()->AddObject(fOutTracks); + } + } + else { + fOutTracks = fTracks; + } + } + } + + if (fNClusters > 0 && !fClusters && !fCaloName.IsNull()) { + fClusters = dynamic_cast(InputEvent()->FindListObject(fCaloName)); + + if (!fClusters) { + AliError(Form("%s: Couldn't retrieve clusters with name %s!", GetName(), fCaloName.Data())); + return kFALSE; + } + + if (!fClusters->GetClass()->GetBaseClass("AliVCluster")) { + AliError(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloName.Data())); + fClusters = 0; + return kFALSE; + } + + if (!fOutClusters) { + fOutCaloName = fCaloName; + if (fCopyArray) { + fOutCaloName += fSuffix; + fOutClusters = new TClonesArray(fClusters->GetClass()->GetName(), fClusters->GetSize()); + fOutClusters->SetName(fOutCaloName); + if (InputEvent()->FindListObject(fOutCaloName)) { + AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCaloName.Data())); + return kFALSE; + } + else { + InputEvent()->AddObject(fOutClusters); + } + } + else { + fOutClusters = fClusters; + } + } + } + +if (!fGeom && (fClusters || fCaloCells)) { + if (fGeomName.Length() > 0) { + fGeom = AliEMCALGeometry::GetInstance(fGeomName); + if (!fGeom) + AliError(Form("Could not get geometry with name %s!", fGeomName.Data())); + } else { + fGeom = AliEMCALGeometry::GetInstance(); + if (!fGeom) + AliError("Could not get default geometry!"); + } + } + + const Double_t EmcalMinEta = fGeom->GetArm1EtaMin(); + const Double_t EmcalMaxEta = fGeom->GetArm1EtaMax(); + const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad(); + const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad(); + + if (fEtaMax > EmcalMaxEta) fEtaMax = EmcalMaxEta; + if (fEtaMax < EmcalMinEta) fEtaMax = EmcalMinEta; + if (fEtaMin > EmcalMaxEta) fEtaMin = EmcalMaxEta; + if (fEtaMin < EmcalMinEta) fEtaMin = EmcalMinEta; + + if (fPhiMax > EmcalMaxPhi) fPhiMax = EmcalMaxPhi; + if (fPhiMax < EmcalMinPhi) fPhiMax = EmcalMinPhi; + if (fPhiMin > EmcalMaxPhi) fPhiMin = EmcalMaxPhi; + if (fPhiMin < EmcalMinPhi) fPhiMin = EmcalMinPhi; + + return kTRUE; +} + +//________________________________________________________________________ +Int_t AliJetModelBaseTask::SetNumberOfOutCells(Int_t n) +{ + if (fOutCaloCells->GetNumberOfCells() < n) { + fOutCaloCells->DeleteContainer(); + fOutCaloCells->CreateContainer(n); + } + else { + fOutCaloCells->SetNumberOfCells(n); + } + + fAddedCells = 0; + + return n; +} + +//________________________________________________________________________ +void AliJetModelBaseTask::CopyCells() +{ + for (Short_t i = 0; i < fCaloCells->GetNumberOfCells(); i++) { + Short_t mclabel = -1; + Double_t efrac = 0.; + Double_t time = -1; + Short_t cellNum = -1; + Double_t amp = -1; + + fCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac); + fOutCaloCells->SetCell(i, cellNum, amp, time, mclabel, efrac); + } + + fAddedCells = fCaloCells->GetNumberOfCells(); +} + +//________________________________________________________________________ +Int_t AliJetModelBaseTask::AddCell(Double_t e, Double_t eta, Double_t phi) +{ + // Add a cell to the event. + + Int_t absId = 0; + if (eta < -100 || phi < 0) { + GetRandomCell(eta, phi, absId); + } + else { + fGeom->EtaPhiFromIndex(absId, eta, phi); + } + + if (absId == -1) { + AliWarning(Form("Unable to embed cell in eta = %f, phi = %f!" + " Maybe the eta-phi range is not inside the EMCal acceptance (eta = [%f, %f], phi = [%f, %f])", + eta, phi, fEtaMin, fEtaMax, fPhiMin, fPhiMax)); + return 0; + } + + if (e < 0) { + Double_t pt = GetRandomPt(); + TLorentzVector nPart; + nPart.SetPtEtaPhiM(pt, eta, phi, 0); + e = nPart.E(); + } + + return AddCell(e, absId); +} + +//________________________________________________________________________ +Int_t AliJetModelBaseTask::AddCell(Double_t e, Int_t absId, Double_t time) +{ + // Add a cell to the event. + + Int_t label = fMarkMC ? 100 : -1; + + Bool_t r = fOutCaloCells->SetCell(fAddedCells, absId, e, time, label, 0); + + if (r) { + fAddedCells++; + return fAddedCells; + } + else { + return 0; + } +} + + //________________________________________________________________________ AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Double_t eta, Double_t phi) { @@ -183,7 +469,7 @@ AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Int_t absId) } //________________________________________________________________________ -AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t phi) +AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t phi, Byte_t type, Double_t etaemc, Double_t phiemc, Bool_t ise) { // Add a track to the event. @@ -196,16 +482,18 @@ AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t if (phi < 0) phi = GetRandomPhi(); - Int_t label = fMarkMC ? 100 : 0; + Int_t label = fMarkMC ? 100 : -1; AliPicoTrack *track = new ((*fOutTracks)[nTracks]) AliPicoTrack(pt, eta, phi, - 1, - label, // MC flag! - 0, - 0, - kFALSE); + 1, + label, // MC flag! + type, + etaemc, + phiemc, + ise); + return track; } @@ -214,11 +502,10 @@ void AliJetModelBaseTask::CopyClusters() { // Copy all the clusters in the new collection - Bool_t esdMode = (Bool_t)(fClusters->GetClass()->GetBaseClass("AliESDCaloCluster") != 0); const Int_t nClusters = fClusters->GetEntriesFast(); Int_t nCopiedClusters = 0; - if (esdMode) { + if (fEsdMode) { for (Int_t i = 0; i < nClusters; ++i) { AliESDCaloCluster *esdcluster = static_cast(fClusters->At(i)); if (!esdcluster || !esdcluster->IsEMCAL()) @@ -252,126 +539,6 @@ void AliJetModelBaseTask::CopyTracks() } } -//________________________________________________________________________ -Bool_t AliJetModelBaseTask::ExecOnce() -{ - // Init task. - - delete gRandom; - gRandom = new TRandom3(0); - - if (fPtMax < fPtMin) { - AliWarning (Form("PtMax (%f) < PtMin (%f), setting PtMax = PtMin = %f", fPtMax, fPtMin, fPtMin)); - fPtMax = fPtMin; - } - - if (fEtaMax < fEtaMin) { - AliWarning (Form("EtaMax (%f) < EtaMin (%f), setting EtaMax = EtaMin = %f", fEtaMax, fEtaMin, fEtaMin)); - fEtaMax = fEtaMin; - } - - if (fPhiMax < fPhiMin) { - AliWarning (Form("PhiMax (%f) < PhiMin (%f), setting PhiMax = PhiMin = %f", fPhiMax, fPhiMin, fPhiMin)); - fPhiMax = fPhiMin; - } - - if (fNTracks > 0 && !fTracks && !fTracksName.IsNull()) { - fTracks = dynamic_cast(InputEvent()->FindListObject(fTracksName)); - if (!fTracks) { - AliError(Form("%s: Couldn't retrieve tracks with name %s!", GetName(), fTracksName.Data())); - return kFALSE; - } - - if (!fTracks->GetClass()->GetBaseClass("AliPicoTrack")) { - AliError(Form("%s: Collection %s does not contain AliPicoTrack objects!", GetName(), fTracksName.Data())); - fTracks = 0; - return kFALSE; - } - - if (!fOutTracks) { - fOutTracksName = fTracksName; - if (fCopyArray) { - fOutTracksName += fSuffix; - fOutTracks = new TClonesArray("AliPicoTrack", fTracks->GetSize()); - fOutTracks->SetName(fOutTracksName); - if (InputEvent()->FindListObject(fOutTracksName)) { - AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutTracksName.Data())); - return kFALSE; - } - else { - InputEvent()->AddObject(fOutTracks); - } - } - else { - fOutTracks = fTracks; - } - } - } - - if (fNClusters > 0 && !fClusters && !fCaloName.IsNull()) { - fClusters = dynamic_cast(InputEvent()->FindListObject(fCaloName)); - - if (!fClusters) { - AliError(Form("%s: Couldn't retrieve clusters with name %s!", GetName(), fCaloName.Data())); - return kFALSE; - } - - if (!fClusters->GetClass()->GetBaseClass("AliVCluster")) { - AliError(Form("%s: Collection %s does not contain AliVCluster objects!", GetName(), fCaloName.Data())); - fClusters = 0; - return kFALSE; - } - - if (!fOutClusters) { - fOutCaloName = fCaloName; - if (fCopyArray) { - fOutCaloName += fSuffix; - fOutClusters = new TClonesArray(fClusters->GetClass()->GetName(), fClusters->GetSize()); - fOutClusters->SetName(fOutCaloName); - if (InputEvent()->FindListObject(fOutCaloName)) { - AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCaloName.Data())); - return kFALSE; - } - else { - InputEvent()->AddObject(fOutClusters); - } - } - else { - fOutClusters = fClusters; - } - } - - if (!fGeom) { - if (fGeomName.Length() > 0) { - fGeom = AliEMCALGeometry::GetInstance(fGeomName); - if (!fGeom) - AliError(Form("Could not get geometry with name %s!", fGeomName.Data())); - } else { - fGeom = AliEMCALGeometry::GetInstance(); - if (!fGeom) - AliError("Could not get default geometry!"); - } - } - - const Double_t EmcalMinEta = fGeom->GetArm1EtaMin(); - const Double_t EmcalMaxEta = fGeom->GetArm1EtaMax(); - const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad(); - const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad(); - - if (fEtaMax > EmcalMaxEta) fEtaMax = EmcalMaxEta; - if (fEtaMax < EmcalMinEta) fEtaMax = EmcalMinEta; - if (fEtaMin > EmcalMaxEta) fEtaMin = EmcalMaxEta; - if (fEtaMin < EmcalMinEta) fEtaMin = EmcalMinEta; - - if (fPhiMax > EmcalMaxPhi) fPhiMax = EmcalMaxPhi; - if (fPhiMax < EmcalMinPhi) fPhiMax = EmcalMinPhi; - if (fPhiMin > EmcalMaxPhi) fPhiMin = EmcalMaxPhi; - if (fPhiMin < EmcalMinPhi) fPhiMin = EmcalMinPhi; - } - - return kTRUE; -} - //________________________________________________________________________ void AliJetModelBaseTask::GetRandomCell(Double_t &eta, Double_t &phi, Int_t &absId) { diff --git a/PWGJE/EMCALJetTasks/AliJetModelBaseTask.h b/PWGJE/EMCALJetTasks/AliJetModelBaseTask.h index 523db00490d..f589e89ab47 100644 --- a/PWGJE/EMCALJetTasks/AliJetModelBaseTask.h +++ b/PWGJE/EMCALJetTasks/AliJetModelBaseTask.h @@ -7,6 +7,7 @@ class TClonesArray; class AliEMCALGeometry; class AliVCluster; class AliPicoTrack; +class AliVCaloCells; #include #include @@ -16,37 +17,48 @@ class AliPicoTrack; class AliJetModelBaseTask : public AliAnalysisTaskSE { public: AliJetModelBaseTask(); - AliJetModelBaseTask(const char *name); + AliJetModelBaseTask(const char *name, Bool_t drawqa=kFALSE); virtual ~AliJetModelBaseTask(); void UserExec(Option_t* /*option*/); + void UserCreateOutputObjects(); - void SetClusName(const char *n) { fCaloName = n; } - void SetCopyArray(Bool_t copy) { fCopyArray = copy; } void SetEtaRange(Float_t min, Float_t max) { fEtaMin = min; fEtaMax = max; } - void SetGeometryName(const char *n) { fGeomName = n; } - void SetMarkMC(Bool_t m) { fMarkMC = m; } - void SetNClusters(Int_t n) { fNClusters = n; } - void SetNTracks(Int_t n) { fNTracks = n; } 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(TH1 *f) { fPtSpectrum = f; } void SetPtSpectrum(TF1 *f) { fPtSpectrum = new TH1F("ptSpectrum","ptSpectrum",250,f->GetXmin(),f->GetXmax()); fPtSpectrum->Add(f); } - void SetSuffix(const char *s) { fSuffix = s; } + + void SetCopyArray(Bool_t copy) { fCopyArray = copy; } void SetTracksName(const char *n) { fTracksName = n; } + void SetClusName(const char *n) { fCaloName = n; } + void SetCellsName(const char *n) { fCellsName = n; } + void SetSuffix(const char *s) { fSuffix = s; } + + void SetGeometryName(const char *n) { fGeomName = n; } + void SetMarkMC(Bool_t m) { fMarkMC = m; } + virtual void SetNClusters(Int_t n) { fNClusters = n; } + virtual void SetNCells(Int_t n) { fNCells = n; } + virtual void SetNTracks(Int_t n) { fNTracks = n; } + protected: + Int_t SetNumberOfOutCells(Int_t n); // set the number of cells + Int_t AddCell(Double_t e = -1, Double_t eta = -999, Double_t phi = -1); // add a cell; if values are -1 generate random parameters + Int_t AddCell(Double_t e, Int_t absId, Double_t time = 0); // add a cell with given energy, position and times AliVCluster *AddCluster(Double_t e = -1, Double_t eta = -999, Double_t phi = -1); // add a cluster; if values are -1 generate random parameters AliVCluster *AddCluster(Double_t e, Int_t absId); // add a cluster with given energy and position - AliPicoTrack *AddTrack(Double_t pt = -1, Double_t eta = -999, Double_t phi = -1); // add a track; if values are -1 generate random parameters + AliPicoTrack *AddTrack(Double_t pt = -1, Double_t eta = -999, Double_t phi = -1, + Byte_t type=0, Double_t etaemc=0, Double_t phiemc=0, Bool_t ise=kFALSE); // add a track; if values are -1 generate random parameters + void CopyCells(); void CopyClusters(); void CopyTracks(); - virtual Bool_t ExecOnce(); void GetRandomCell(Double_t &eta, Double_t &phi, Int_t &absId); // generate a random cell in the calorimeter Double_t GetRandomEta(); // generate a random eta value in the given range Double_t GetRandomPhi(); // generate a random phi value in the given range Double_t GetRandomPt(); // generate a random pt value in the given range + virtual Bool_t ExecOnce(); // intialize task virtual void Run(); // do jet model action TString fGeomName; // EMCal geometry name @@ -54,6 +66,8 @@ class AliJetModelBaseTask : public AliAnalysisTaskSE { TString fOutTracksName; // name of output track collection TString fCaloName; // name of calo cluster collection TString fOutCaloName; // name of output cluster collection + TString fCellsName; // name of calo cells collection + TString fOutCellsName; // name of output cells collection TString fSuffix; // suffix to add in the name of new collections Float_t fEtaMin; // eta minimum value Float_t fEtaMax; // eta maximum value @@ -63,20 +77,27 @@ class AliJetModelBaseTask : public AliAnalysisTaskSE { Float_t fPtMax; // pt maximum value Bool_t fCopyArray; // whether or not the array will be copied to a new one before modelling Int_t fNClusters; // how many clusters are being processed + Int_t fNCells; // how many cells are being processed Int_t fNTracks; // how many tracks are being processed Bool_t fMarkMC; // whether or not mark new tracks/cluster as MC TH1 *fPtSpectrum; // pt spectrum parametrization to extract random pt values + Bool_t fQAhistos; // draw QA histograms Bool_t fIsInit; //=true if initialized AliEMCALGeometry *fGeom; //!pointer to EMCal geometry TClonesArray *fClusters; //!cluster collection TClonesArray *fOutClusters; //!output cluster collection TClonesArray *fTracks; //!track collection TClonesArray *fOutTracks; //!output track collection + AliVCaloCells *fCaloCells; //!cells collection + AliVCaloCells *fOutCaloCells; //!output cells collection + Int_t fAddedCells; //!number of added cells + Bool_t fEsdMode; //!ESD/AOD mode + TList *fOutput; //!output list for QA histograms private: AliJetModelBaseTask(const AliJetModelBaseTask&); // not implemented AliJetModelBaseTask &operator=(const AliJetModelBaseTask&); // not implemented - ClassDef(AliJetModelBaseTask, 4) // Jet modelling task + ClassDef(AliJetModelBaseTask, 5) // Jet modelling task }; #endif diff --git a/PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx b/PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx index e7955c50f27..aafdd5e4dfd 100644 --- a/PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx +++ b/PWGJE/EMCALJetTasks/AliJetResponseMaker.cxx @@ -6,20 +6,17 @@ #include "AliJetResponseMaker.h" -#include #include #include #include -#include +#include #include -#include "AliAnalysisManager.h" -#include "AliCentrality.h" -#include "AliFJWrapper.h" #include "AliVCluster.h" #include "AliVTrack.h" #include "AliEmcalJet.h" #include "AliGenPythiaEventHeader.h" +#include "AliLog.h" #include "AliMCEvent.h" ClassImp(AliJetResponseMaker) @@ -47,24 +44,27 @@ AliJetResponseMaker::AliJetResponseMaker() : fMCJets(0), fHistNTrials(0), fHistEvents(0), - fHistMCJetPhiEta(0), - fHistMCJetsPt(0), - fHistMCJetPhiEtaFiducial(0), - fHistMCJetsPtFiducial(0), + fHistMCJetsPhiEta(0), + fHistMCJetsPtArea(0), + fHistMCJetsPhiEtaFiducial(0), + fHistMCJetsPtAreaFiducial(0), fHistMCJetsNEFvsPt(0), fHistMCJetsZvsPt(0), - fHistJetPhiEta(0), - fHistJetsPt(0), + fHistJetsPhiEta(0), + fHistJetsPtArea(0), + fHistJetsCorrPtArea(0), fHistJetsNEFvsPt(0), fHistJetsZvsPt(0), - fHistClosestDistance(0), - fHistClosestDeltaPhi(0), - fHistClosestDeltaEta(0), - fHistClosestDeltaPt(0), - fHistNonMatchedMCJetPt(0), - fHistNonMatchedJetPt(0), + fHistMatchingLevelMCPt(0), + fHistClosestDeltaEtaPhiMCPt(0), + fHistClosestDeltaPtMCPt(0), + fHistClosestDeltaCorrPtMCPt(0), + fHistNonMatchedMCJetsPtArea(0), + fHistNonMatchedJetsPtArea(0), + fHistNonMatchedJetsCorrPtArea(0), fHistPartvsDetecPt(0), - fHistMissedMCJets(0) + fHistPartvsDetecCorrPt(0), + fHistMissedMCJetsPtArea(0) { // Default constructor. @@ -98,24 +98,27 @@ AliJetResponseMaker::AliJetResponseMaker(const char *name) : fMCJets(0), fHistNTrials(0), fHistEvents(0), - fHistMCJetPhiEta(0), - fHistMCJetsPt(0), - fHistMCJetPhiEtaFiducial(0), - fHistMCJetsPtFiducial(0), + fHistMCJetsPhiEta(0), + fHistMCJetsPtArea(0), + fHistMCJetsPhiEtaFiducial(0), + fHistMCJetsPtAreaFiducial(0), fHistMCJetsNEFvsPt(0), fHistMCJetsZvsPt(0), - fHistJetPhiEta(0), - fHistJetsPt(0), + fHistJetsPhiEta(0), + fHistJetsPtArea(0), + fHistJetsCorrPtArea(0), fHistJetsNEFvsPt(0), fHistJetsZvsPt(0), - fHistClosestDistance(0), - fHistClosestDeltaPhi(0), - fHistClosestDeltaEta(0), - fHistClosestDeltaPt(0), - fHistNonMatchedMCJetPt(0), - fHistNonMatchedJetPt(0), + fHistMatchingLevelMCPt(0), + fHistClosestDeltaEtaPhiMCPt(0), + fHistClosestDeltaPtMCPt(0), + fHistClosestDeltaCorrPtMCPt(0), + fHistNonMatchedMCJetsPtArea(0), + fHistNonMatchedJetsPtArea(0), + fHistNonMatchedJetsCorrPtArea(0), fHistPartvsDetecPt(0), - fHistMissedMCJets(0) + fHistPartvsDetecCorrPt(0), + fHistMissedMCJetsPtArea(0) { // Standard constructor. @@ -165,95 +168,121 @@ void AliJetResponseMaker::UserCreateOutputObjects() fHistEvents->GetXaxis()->SetBinLabel(i, Form("%d-%d",ptHardLo[i-1],ptHardHi[i-1])); } - fHistJetPhiEta = new TH2F("fHistJetPhiEta", "fHistJetPhiEta", 20, -2, 2, 32, 0, 6.4); - fHistJetPhiEta->GetXaxis()->SetTitle("#eta"); - fHistJetPhiEta->GetYaxis()->SetTitle("#phi"); - fOutput->Add(fHistJetPhiEta); + fHistJetsPhiEta = new TH2F("fHistJetsPhiEta", "fHistJetsPhiEta", 20, -2, 2, 32, 0, 6.4); + fHistJetsPhiEta->GetXaxis()->SetTitle("#eta"); + fHistJetsPhiEta->GetYaxis()->SetTitle("#phi"); + fOutput->Add(fHistJetsPhiEta); - fHistJetsPt = new TH1F("fHistJetsPt", "fHistJetsPt", fNbins, fMinBinPt, fMaxBinPt); - fHistJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]"); - fHistJetsPt->GetYaxis()->SetTitle("counts"); - fOutput->Add(fHistJetsPt); + fHistJetsPtArea = new TH2F("fHistJetsPtArea", "fHistJetsPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt); + fHistJetsPtArea->GetXaxis()->SetTitle("area"); + fHistJetsPtArea->GetYaxis()->SetTitle("p_{T}^{rec} (GeV/c)"); + fHistJetsPtArea->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistJetsPtArea); + + fHistJetsCorrPtArea = new TH2F("fHistJetsCorrPtArea", "fHistJetsCorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, (Int_t)(1.5*fNbins), -fMaxBinPt/2, fMaxBinPt); + fHistJetsCorrPtArea->GetXaxis()->SetTitle("area"); + fHistJetsCorrPtArea->GetYaxis()->SetTitle("p_{T}^{corr} (GeV/c)"); + fHistJetsCorrPtArea->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistJetsCorrPtArea); fHistJetsZvsPt = new TH2F("fHistJetsZvsPt", "fHistJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); fHistJetsZvsPt->GetXaxis()->SetTitle("Z"); - fHistJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]"); + fHistJetsZvsPt->GetYaxis()->SetTitle("p_{T}^{rec} (GeV/c)"); fOutput->Add(fHistJetsZvsPt); fHistJetsNEFvsPt = new TH2F("fHistJetsNEFvsPt", "fHistJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); fHistJetsNEFvsPt->GetXaxis()->SetTitle("NEF"); - fHistJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]"); + fHistJetsNEFvsPt->GetYaxis()->SetTitle("p_{T}^{rec} (GeV/c)"); fOutput->Add(fHistJetsNEFvsPt); - fHistMCJetPhiEta = new TH2F("fHistMCJetPhiEta", "fHistMCJetPhiEta", 20, -2, 2, 32, 0, 6.4); - fHistMCJetPhiEta->GetXaxis()->SetTitle("#eta"); - fHistMCJetPhiEta->GetYaxis()->SetTitle("#phi"); - fOutput->Add(fHistMCJetPhiEta); + fHistMCJetsPhiEta = new TH2F("fHistMCJetsPhiEta", "fHistMCJetsPhiEta", 20, -2, 2, 32, 0, 6.4); + fHistMCJetsPhiEta->GetXaxis()->SetTitle("#eta"); + fHistMCJetsPhiEta->GetYaxis()->SetTitle("#phi"); + fOutput->Add(fHistMCJetsPhiEta); - fHistMCJetsPt = new TH1F("fHistMCJetsPt", "fHistMCJetsPt", fNbins, fMinBinPt, fMaxBinPt); - fHistMCJetsPt->GetXaxis()->SetTitle("p_{T} [GeV/c]"); - fHistMCJetsPt->GetYaxis()->SetTitle("counts"); - fOutput->Add(fHistMCJetsPt); - - fHistMCJetPhiEtaFiducial = new TH2F("fHistMCJetPhiEtaFiducial", "fHistMCJetPhiEtaFiducial", 20, -2, 2, 32, 0, 6.4); - fHistMCJetPhiEtaFiducial->GetXaxis()->SetTitle("#eta"); - fHistMCJetPhiEtaFiducial->GetYaxis()->SetTitle("#phi"); - fOutput->Add(fHistMCJetPhiEtaFiducial); + fHistMCJetsPtArea = new TH2F("fHistMCJetsPtArea", "fHistMCJetsPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt); + fHistMCJetsPtArea->GetXaxis()->SetTitle("area"); + fHistMCJetsPtArea->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)"); + fHistMCJetsPtArea->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistMCJetsPtArea); + + fHistMCJetsPhiEtaFiducial = new TH2F("fHistMCJetsPhiEtaFiducial", "fHistMCJetsPhiEtaFiducial", 20, -2, 2, 32, 0, 6.4); + fHistMCJetsPhiEtaFiducial->GetXaxis()->SetTitle("#eta"); + fHistMCJetsPhiEtaFiducial->GetYaxis()->SetTitle("#phi"); + fOutput->Add(fHistMCJetsPhiEtaFiducial); - fHistMCJetsPtFiducial = new TH1F("fHistMCJetsPtFiducial", "fHistMCJetsPtFiducial", fNbins, fMinBinPt, fMaxBinPt); - fHistMCJetsPtFiducial->GetXaxis()->SetTitle("p_{T} [GeV/c]"); - fHistMCJetsPtFiducial->GetYaxis()->SetTitle("counts"); - fOutput->Add(fHistMCJetsPtFiducial); + fHistMCJetsPtAreaFiducial = new TH2F("fHistMCJetsPtAreaFiducial", "fHistMCJetsPtAreaFiducial", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt); + fHistMCJetsPtAreaFiducial->GetXaxis()->SetTitle("area"); + fHistMCJetsPtAreaFiducial->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)"); + fHistMCJetsPtAreaFiducial->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistMCJetsPtAreaFiducial); fHistMCJetsZvsPt = new TH2F("fHistMCJetsZvsPt", "fHistMCJetsZvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); fHistMCJetsZvsPt->GetXaxis()->SetTitle("Z"); - fHistMCJetsZvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]"); + fHistMCJetsZvsPt->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)"); fOutput->Add(fHistMCJetsZvsPt); fHistMCJetsNEFvsPt = new TH2F("fHistMCJetsNEFvsPt", "fHistMCJetsNEFvsPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); fHistMCJetsNEFvsPt->GetXaxis()->SetTitle("NEF"); - fHistMCJetsNEFvsPt->GetYaxis()->SetTitle("p_{T} [GeV/c]"); + fHistMCJetsNEFvsPt->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)"); fOutput->Add(fHistMCJetsNEFvsPt); - fHistClosestDistance = new TH1F("fHistClosestDistance", "fHistClosestDistance", 50, 0, fMaxDistance * 1.2); - fHistClosestDistance->GetXaxis()->SetTitle("d"); - fHistClosestDistance->GetYaxis()->SetTitle("counts"); - fOutput->Add(fHistClosestDistance); - - fHistClosestDeltaPhi = new TH1F("fHistClosestDeltaPhi", "fHistClosestDeltaPhi", 64, -1.6, 4.8); - fHistClosestDeltaPhi->GetXaxis()->SetTitle("#Delta#phi"); - fHistClosestDeltaPhi->GetYaxis()->SetTitle("counts"); - fOutput->Add(fHistClosestDeltaPhi); - - fHistClosestDeltaEta = new TH1F("fHistClosestDeltaEta", "fHistClosestDeltaEta", TMath::CeilNint(fJetMaxEta - fJetMinEta) * 20, fJetMinEta * 2, fJetMaxEta * 2); - fHistClosestDeltaEta->GetXaxis()->SetTitle("#Delta#eta"); - fHistClosestDeltaEta->GetYaxis()->SetTitle("counts"); - fOutput->Add(fHistClosestDeltaEta); - - fHistClosestDeltaPt = new TH1F("fHistClosestDeltaPt", "fHistClosestDeltaPt", fNbins, -fMaxBinPt / 2, fMaxBinPt / 2); - fHistClosestDeltaPt->GetXaxis()->SetTitle("#Delta p_{T}"); - fHistClosestDeltaPt->GetYaxis()->SetTitle("counts"); - fOutput->Add(fHistClosestDeltaPt); - - fHistNonMatchedMCJetPt = new TH1F("fHistNonMatchedMCJetPt", "fHistNonMatchedMCJetPt", fNbins, fMinBinPt, fMaxBinPt); - fHistNonMatchedMCJetPt->GetXaxis()->SetTitle("p_{T} [GeV/c]"); - fHistNonMatchedMCJetPt->GetYaxis()->SetTitle("counts"); - fOutput->Add(fHistNonMatchedMCJetPt); - - fHistNonMatchedJetPt = new TH1F("fHistNonMatchedJetPt", "fHistNonMatchedJetPt", fNbins, fMinBinPt, fMaxBinPt); - fHistNonMatchedJetPt->GetXaxis()->SetTitle("p_{T} [GeV/c]"); - fHistNonMatchedJetPt->GetYaxis()->SetTitle("counts"); - fOutput->Add(fHistNonMatchedJetPt); + fHistMatchingLevelMCPt = new TH2F("fHistMatchingLevelMCPt", "fHistMatchingLevelMCPt", fNbins, 0, 1.2, fNbins, fMinBinPt, fMaxBinPt); + fHistMatchingLevelMCPt->GetXaxis()->SetTitle("Matching level"); + fHistMatchingLevelMCPt->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)"); + fOutput->Add(fHistMatchingLevelMCPt); + + fHistClosestDeltaEtaPhiMCPt = new TH3F("fHistClosestDeltaEtaPhiMCPt", "fHistClosestDeltaEtaPhiMCPt", TMath::CeilNint(fJetMaxEta - fJetMinEta) * 20, fJetMinEta * 2, fJetMaxEta * 2, 64, -1.6, 4.8, fNbins, fMinBinPt, fMaxBinPt); + fHistClosestDeltaEtaPhiMCPt->GetXaxis()->SetTitle("#Delta#eta"); + fHistClosestDeltaEtaPhiMCPt->GetYaxis()->SetTitle("#Delta#phi"); + fHistClosestDeltaEtaPhiMCPt->GetZaxis()->SetTitle("p_{T}^{gen}"); + fOutput->Add(fHistClosestDeltaEtaPhiMCPt); + + fHistClosestDeltaPtMCPt = new TH2F("fHistClosestDeltaPtMCPt", "fHistClosestDeltaPtMCPt", fNbins, fMinBinPt, fMaxBinPt, (Int_t)(fNbins*1.5), -fMaxBinPt / 2, fMaxBinPt); + fHistClosestDeltaPtMCPt->GetXaxis()->SetTitle("p_{T}^{gen}"); + fHistClosestDeltaPtMCPt->GetYaxis()->SetTitle("#Deltap_{T}^{rec} (GeV/c)"); + fHistClosestDeltaPtMCPt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistClosestDeltaPtMCPt); + + fHistClosestDeltaCorrPtMCPt = new TH2F("fHistClosestDeltaCorrPtMCPt", "fHistClosestDeltaCorrPtMCPt", fNbins, fMinBinPt, fMaxBinPt, (Int_t)(fNbins*1.5), -fMaxBinPt / 2, fMaxBinPt); + fHistClosestDeltaCorrPtMCPt->GetXaxis()->SetTitle("p_{T}^{gen}"); + fHistClosestDeltaCorrPtMCPt->GetYaxis()->SetTitle("#Deltap_{T}^{corr} (GeV/c)"); + fHistClosestDeltaCorrPtMCPt->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistClosestDeltaCorrPtMCPt); + + fHistNonMatchedMCJetsPtArea = new TH2F("fHistNonMatchedMCJetsPtArea", "fHistNonMatchedMCJetsPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt); + fHistNonMatchedMCJetsPtArea->GetXaxis()->SetTitle("area"); + fHistNonMatchedMCJetsPtArea->GetYaxis()->SetTitle("p_{T}^{gen} (GeV/c)"); + fHistNonMatchedMCJetsPtArea->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistNonMatchedMCJetsPtArea); + + fHistNonMatchedJetsPtArea = new TH2F("fHistNonMatchedJetsPtArea", "fHistNonMatchedJetsPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt); + fHistNonMatchedJetsPtArea->GetXaxis()->SetTitle("area"); + fHistNonMatchedJetsPtArea->GetYaxis()->SetTitle("p_{T}^{rec} (GeV/c)"); + fHistNonMatchedJetsPtArea->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistNonMatchedJetsPtArea); + + fHistNonMatchedJetsCorrPtArea = new TH2F("fHistNonMatchedJetsCorrPtArea", "fHistNonMatchedJetsCorrPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, (Int_t)(1.5*fNbins), -fMaxBinPt/2, fMaxBinPt); + fHistNonMatchedJetsCorrPtArea->GetXaxis()->SetTitle("area"); + fHistNonMatchedJetsCorrPtArea->GetYaxis()->SetTitle("p_{T}^{corr} (GeV/c)"); + fHistNonMatchedJetsCorrPtArea->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistNonMatchedJetsCorrPtArea); fHistPartvsDetecPt = new TH2F("fHistPartvsDetecPt", "fHistPartvsDetecPt", fNbins, fMinBinPt, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt); fHistPartvsDetecPt->GetXaxis()->SetTitle("p_{T}^{rec}"); fHistPartvsDetecPt->GetYaxis()->SetTitle("p_{T}^{gen}"); fOutput->Add(fHistPartvsDetecPt); - fHistMissedMCJets = new TH1F("fHistMissedMCJets", "fHistMissedMCJets", fNbins, fMinBinPt, fMaxBinPt); - fHistMissedMCJets->GetXaxis()->SetTitle("p_{T} [GeV/c]"); - fHistMissedMCJets->GetYaxis()->SetTitle("counts"); - fOutput->Add(fHistMissedMCJets); + fHistPartvsDetecCorrPt = new TH2F("fHistPartvsDetecCorrPt", "fHistPartvsDetecCorrPt", (Int_t)(1.5*fNbins), -fMaxBinPt/2, fMaxBinPt, fNbins, fMinBinPt, fMaxBinPt); + fHistPartvsDetecCorrPt->GetXaxis()->SetTitle("p_{T}^{corr}"); + fHistPartvsDetecCorrPt->GetYaxis()->SetTitle("p_{T}^{gen}"); + fOutput->Add(fHistPartvsDetecCorrPt); + + fHistMissedMCJetsPtArea = new TH2F("fHistMissedMCJetsPtArea", "fHistMissedMCJetsPtArea", 40, 0, fJetRadius * fJetRadius * TMath::Pi() * 3, fNbins, fMinBinPt, fMaxBinPt); + fHistMissedMCJetsPtArea->GetXaxis()->SetTitle("area"); + fHistMissedMCJetsPtArea->GetYaxis()->SetTitle("p_{T} (GeV/c)"); + fHistMissedMCJetsPtArea->GetZaxis()->SetTitle("counts"); + fOutput->Add(fHistMissedMCJetsPtArea); PostData(1, fOutput); // Post data for ALL output slots > 0 here, to get at least an empty histogram } @@ -506,30 +535,30 @@ Bool_t AliJetResponseMaker::FillHistograms() if (!AcceptBiasJet(jet->MatchedJet()) || jet->MatchedJet()->MaxTrackPt() > fMaxTrackPt || jet->MatchedJet()->MaxClusterPt() > fMaxClusterPt || jet->MatchedJet()->Pt() > fMaxBinPt) { - fHistMissedMCJets->Fill(jet->Pt(), fEventWeight); + fHistMissedMCJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight); } else { - fHistClosestDistance->Fill(jet->ClosestJetDistance(), fEventWeight); + fHistMatchingLevelMCPt->Fill(jet->ClosestJetDistance(), jet->Pt(), fEventWeight); Double_t deta = jet->MatchedJet()->Eta() - jet->Eta(); - fHistClosestDeltaEta->Fill(deta, fEventWeight); - Double_t dphi = jet->MatchedJet()->Phi() - jet->Phi(); - fHistClosestDeltaPhi->Fill(dphi, fEventWeight); + fHistClosestDeltaEtaPhiMCPt->Fill(deta, dphi, jet->Pt(), fEventWeight); Double_t dpt = jet->MatchedJet()->Pt() - jet->Pt(); - fHistClosestDeltaPt->Fill(dpt, fEventWeight); + fHistClosestDeltaPtMCPt->Fill(jet->Pt(), dpt, fEventWeight); + fHistClosestDeltaCorrPtMCPt->Fill(jet->Pt(), dpt - fRhoVal * jet->MatchedJet()->Area(), fEventWeight); fHistPartvsDetecPt->Fill(jet->MatchedJet()->Pt(), jet->Pt(), fEventWeight); + fHistPartvsDetecCorrPt->Fill(jet->MatchedJet()->Pt() - fRhoVal * jet->MatchedJet()->Area(), jet->Pt(), fEventWeight); } } else { - fHistNonMatchedMCJetPt->Fill(jet->Pt(), fEventWeight); - fHistMissedMCJets->Fill(jet->Pt(), fEventWeight); + fHistNonMatchedMCJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight); + fHistMissedMCJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight); } - fHistMCJetsPt->Fill(jet->Pt(), fEventWeight); - fHistMCJetPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight); + fHistMCJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight); + fHistMCJetsPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight); fHistMCJetsNEFvsPt->Fill(jet->NEF(), jet->Pt(), fEventWeight); @@ -544,8 +573,8 @@ Bool_t AliJetResponseMaker::FillHistograms() if (jet->Eta() < fJetMinEta || jet->Eta() > fJetMaxEta || jet->Phi() < fJetMinPhi || jet->Phi() > fJetMaxPhi) continue; - fHistMCJetsPtFiducial->Fill(jet->Pt(), fEventWeight); - fHistMCJetPhiEtaFiducial->Fill(jet->Eta(), jet->Phi(), fEventWeight); + fHistMCJetsPtAreaFiducial->Fill(jet->Area(), jet->Pt(), fEventWeight); + fHistMCJetsPhiEtaFiducial->Fill(jet->Eta(), jet->Phi(), fEventWeight); } const Int_t nJets = fJets->GetEntriesFast(); @@ -569,12 +598,14 @@ Bool_t AliJetResponseMaker::FillHistograms() continue; if (!jet->MatchedJet()) { - fHistNonMatchedJetPt->Fill(jet->Pt(), fEventWeight); + fHistNonMatchedJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight); + fHistNonMatchedJetsCorrPtArea->Fill(jet->Area(), jet->Pt() - fRhoVal * jet->Area(), fEventWeight); } - fHistJetsPt->Fill(jet->Pt(), fEventWeight); + fHistJetsPtArea->Fill(jet->Area(), jet->Pt(), fEventWeight); + fHistJetsCorrPtArea->Fill(jet->Area(), jet->Pt() - fRhoVal * jet->Area(), fEventWeight); - fHistJetPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight); + fHistJetsPhiEta->Fill(jet->Eta(), jet->Phi(), fEventWeight); fHistJetsNEFvsPt->Fill(jet->NEF(), jet->Pt(), fEventWeight); diff --git a/PWGJE/EMCALJetTasks/AliJetResponseMaker.h b/PWGJE/EMCALJetTasks/AliJetResponseMaker.h index b25f0ceae5f..629ef7699bf 100644 --- a/PWGJE/EMCALJetTasks/AliJetResponseMaker.h +++ b/PWGJE/EMCALJetTasks/AliJetResponseMaker.h @@ -5,8 +5,9 @@ class AliGenPythiaEventHeader; class TClonesArray; -class TH1F; -class TH2F; +class TH1; +class TH2; +class TH3; #include "AliAnalysisTaskEmcalJet.h" @@ -36,55 +37,58 @@ class AliJetResponseMaker : public AliAnalysisTaskEmcalJet { Bool_t RetrieveEventObjects(); Bool_t Run(); - TString fMCTracksName; // name of MC particle collection - TString fMCJetsName; // name of MC jet collection - Double_t fMaxDistance; // maximum distance between matched particle and detector level jets - Bool_t fDoWeighting; // = true, weight using trials and given x section - Bool_t fEventWeightHist; // = true create event weight histogram - Bool_t fMCFiducial; // = true MC jets in fiducial acceptance - Double_t fMCminEta; // MC jets minimum eta - Double_t fMCmaxEta; // MC jets maximum eta - Double_t fMCminPhi; // MC jets minimum phi - Double_t fMCmaxPhi; // MC jets maximum phi - Int_t fSelectPtHardBin; // select one pt hard bin for analysis - Bool_t fDoMatching; // whether or not it should run the matching between MC and rec jets + TString fMCTracksName; // name of MC particle collection + TString fMCJetsName; // name of MC jet collection + Double_t fMaxDistance; // maximum distance between matched particle and detector level jets + Bool_t fDoWeighting; // = true, weight using trials and given x section + Bool_t fEventWeightHist; // = true create event weight histogram + Bool_t fMCFiducial; // = true MC jets in fiducial acceptance + Double_t fMCminEta; // MC jets minimum eta + Double_t fMCmaxEta; // MC jets maximum eta + Double_t fMCminPhi; // MC jets minimum phi + Double_t fMCmaxPhi; // MC jets maximum phi + Int_t fSelectPtHardBin; // select one pt hard bin for analysis + Bool_t fDoMatching; // whether or not it should run the matching between MC and rec jets - AliGenPythiaEventHeader *fPythiaHeader; //!event Pythia header - Double_t fEventWeight; //!event weight - Int_t fPtHardBin; //!event pt hard bin - Int_t fNTrials; //!event trials - TClonesArray *fMCTracks; //!MC particles - TClonesArray *fMCJets; //!MC jets + AliGenPythiaEventHeader *fPythiaHeader; //!event Pythia header + Double_t fEventWeight; //!event weight + Int_t fPtHardBin; //!event pt hard bin + Int_t fNTrials; //!event trials + TClonesArray *fMCTracks; //!MC particles + TClonesArray *fMCJets; //!MC jets // General histograms - TH1F *fHistNTrials; //!total number of trials per pt hard bin - TH1F *fHistEvents; //!total number of events per pt hard bin - TH1F *fHistEventWeight[11]; //!event weight + TH1 *fHistNTrials; //!total number of trials per pt hard bin + TH1 *fHistEvents; //!total number of events per pt hard bin + TH1 *fHistEventWeight[11]; //!event weight // Particle level jets - TH2F *fHistMCJetPhiEta; //!phi-eta distribution of jets - TH1F *fHistMCJetsPt; //!inclusive jet pt spectrum - TH2F *fHistMCJetPhiEtaFiducial; //!phi-eta distribution of jets in fiducial acceptance (plus lead hadron bias) - TH1F *fHistMCJetsPtFiducial; //!inclusive jet pt spectrum in fiducial acceptance (plus lead hadron bias) - TH2F *fHistMCJetsNEFvsPt; //!jet neutral energy fraction vs. jet pt - TH2F *fHistMCJetsZvsPt; //!constituent Pt over Jet Pt ratio vs. jet pt + TH2 *fHistMCJetsPhiEta; //!phi-eta distribution of jets + TH2 *fHistMCJetsPtArea; //!inclusive jet pt vs area histogram + TH2 *fHistMCJetsPhiEtaFiducial; //!phi-eta distribution of jets in fiducial acceptance (plus lead hadron bias) + TH2 *fHistMCJetsPtAreaFiducial; //!inclusive jet pt spectrum in fiducial acceptance (plus lead hadron bias) + TH2 *fHistMCJetsNEFvsPt; //!jet neutral energy fraction vs. jet pt + TH2 *fHistMCJetsZvsPt; //!constituent Pt over Jet Pt ratio vs. jet pt // Detector level jets - TH2F *fHistJetPhiEta; //!phi-eta distribution of jets - TH1F *fHistJetsPt; //!inclusive jet pt spectrum - TH2F *fHistJetsNEFvsPt; //!jet neutral energy fraction vs. jet pt - TH2F *fHistJetsZvsPt; //!constituent Pt over Jet Pt ratio vs. jet pt + TH2 *fHistJetsPhiEta; //!phi-eta distribution of jets + TH2 *fHistJetsPtArea; //!inclusive jet pt vs. area histogram + TH2 *fHistJetsCorrPtArea; //!inclusive jet pt vs. area histogram + TH2 *fHistJetsNEFvsPt; //!jet neutral energy fraction vs. jet pt + TH2 *fHistJetsZvsPt; //!constituent Pt over Jet Pt ratio vs. jet pt // Detector-particle level matching - TH1F *fHistClosestDistance; //!distance between closest particle to detector level jet - TH1F *fHistClosestDeltaPhi; //!delta phi between closest particle to detector level jet - TH1F *fHistClosestDeltaEta; //!delta eta between closest particle to detector level jet - TH1F *fHistClosestDeltaPt; //!delta pt between closest particle to detector level jet - TH1F *fHistNonMatchedMCJetPt; //!non-matched mc jet pt distribution - TH1F *fHistNonMatchedJetPt; //!non-matched jet pt distribution - TH2F *fHistPartvsDetecPt; //!particle vs detector level jet pt - TH1F *fHistMissedMCJets; //!mc jets not measured + TH2 *fHistMatchingLevelMCPt; //!matching level vs. particle level pt + TH3 *fHistClosestDeltaEtaPhiMCPt; //!delta eta-phi between closest particle level to detector level jet vs. particle level pt + TH2 *fHistClosestDeltaPtMCPt; //!delta pt between closest particle level to detector level jet vs. particle level pt + TH2 *fHistClosestDeltaCorrPtMCPt; //!delta pt between closest particle level to detector level jet vs. particle level pt + TH2 *fHistNonMatchedMCJetsPtArea; //!non-matched mc jet pt distribution + TH2 *fHistNonMatchedJetsPtArea; //!non-matched jet pt distribution + TH2 *fHistNonMatchedJetsCorrPtArea; //!non-matched jet pt distribution + TH2 *fHistPartvsDetecPt; //!particle vs detector level jet pt + TH2 *fHistPartvsDetecCorrPt; //!particle vs detector level jet pt + TH2 *fHistMissedMCJetsPtArea; //!mc jets not measured private: AliJetResponseMaker(const AliJetResponseMaker&); // not implemented AliJetResponseMaker &operator=(const AliJetResponseMaker&); // not implemented - ClassDef(AliJetResponseMaker, 8) // Jet response matrix producing task + ClassDef(AliJetResponseMaker, 9) // Jet response matrix producing task }; #endif diff --git a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAQA.cxx b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAQA.cxx index bfb5db76dae..bb282c0a169 100644 --- a/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAQA.cxx +++ b/PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSAQA.cxx @@ -25,6 +25,7 @@ #include "AliLog.h" #include "AliEMCALGeometry.h" #include "AliEMCALGeoParams.h" +#include "AliPicoTrack.h" #include "AliAnalysisTaskSAQA.h" @@ -498,12 +499,16 @@ Float_t AliAnalysisTaskSAQA::DoTrackLoop() } else { fHistTrPhiEtaPt[fCentBin][3]->Fill(track->Eta(), track->Phi(), track->Pt()); + Int_t type = 0; - Int_t label = track->GetLabel(); - if (label >= 0 && label < 3) - fHistTrPhiEtaPt[fCentBin][label]->Fill(track->Eta(), track->Phi(), track->Pt()); + AliPicoTrack* ptrack = dynamic_cast(track); + if (ptrack) + type = ptrack->GetTrackType(); + + if (type >= 0 && type < 3) + fHistTrPhiEtaPt[fCentBin][type]->Fill(track->Eta(), track->Phi(), track->Pt()); else - AliWarning(Form("%s: track label %d not recognized!", GetName(), label)); + AliWarning(Form("%s: track type %d not recognized!", GetName(), type)); } if (!vtrack) diff --git a/PWGJE/EMCALJetTasks/macros/AddTaskJetEmbeddingFromAOD.C b/PWGJE/EMCALJetTasks/macros/AddTaskJetEmbeddingFromAOD.C new file mode 100644 index 00000000000..1f0df3c39c7 --- /dev/null +++ b/PWGJE/EMCALJetTasks/macros/AddTaskJetEmbeddingFromAOD.C @@ -0,0 +1,140 @@ +// $Id: AddTaskJetEmbeddingFromAOD.C $ + +TObjArray* GenerateFileList(const char* list, Int_t nFiles); + +AliJetEmbeddingFromAODTask* AddTaskJetEmbeddingFromAOD( + const char *tracksName = "Tracks", + const char *cellsName = "EMCALCells", + const char *fileList = "files.txt", + const char *aodTreeName = "aodTree", + const char *aodTracksName = "tracks", + const char *aodCellsName = "emcalCells", + const char *runperiod = "", + Bool_t includeNoITS = kTRUE, + Double_t minCent = 0, + Double_t maxCent = 10, + UInt_t mask = AliVEvent::kAny, + const Int_t nTracks = 1234567890, + const Int_t nCells = 1234567890, + const Bool_t copyArray = kTRUE, + const Int_t nFiles = 1234567890, + const Bool_t makeQA = kFALSE, + const Double_t minPt = 0, + const Double_t maxPt = 1000, + const Double_t minEta = -0.9, + const Double_t maxEta = 0.9, + const Double_t minPhi = 0, + const Double_t maxPhi = TMath::Pi() * 2, + const char *taskName = "JetEmbeddingFromAODTask" +) +{ + // Get the pointer to the existing analysis manager via the static access method. + //============================================================================== + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) + { + ::Error("AddTaskJetEmbeddingFromAOD", "No analysis manager to connect to."); + return NULL; + } + + // Check the analysis type using the event handlers connected to the analysis manager. + //============================================================================== + if (!mgr->GetInputEventHandler()) + { + ::Error("AddTaskJetEmbeddingFromAOD", "This task requires an input event handler"); + return NULL; + } + + //------------------------------------------------------- + // Init the task and do settings + //------------------------------------------------------- + + AliJetEmbeddingFromAODTask *jetEmb = new AliJetEmbeddingFromAODTask(taskName,makeQA); + jetEmb->SetTracksName(tracksName); + jetEmb->SetCellsName(cellsName); + jetEmb->SetFileList(GenerateFileList(fileList, nFiles)); + jetEmb->SetAODTreeName(aodTreeName); + jetEmb->SetAODTracksName(aodTracksName); + jetEmb->SetAODCellsName(aodCellsName); + jetEmb->SetCentralityRange(minCent, maxCent); + jetEmb->SetTriggerMask(mask); + jetEmb->SetNCells(nCells); + jetEmb->SetNTracks(nTracks); + jetEmb->SetCopyArray(copyArray); + jetEmb->SetEtaRange(minEta, maxEta); + jetEmb->SetPhiRange(minPhi, maxPhi); + jetEmb->SetPtRange(minPt, maxPt); + + jetEmb->SetIncludeNoITS(includeNoITS); + TString runPeriod(runperiod); + runPeriod.ToLower(); + if (runPeriod == "lhc11h") { + jetEmb->SetAODfilterBits(256,512); // hybrid tracks for LHC11h + } + else if (runPeriod == "lhc11a" || runPeriod == "lhc12a15a" || runPeriod == "lhc12a15e") { + jetEmb->SetAODfilterBits(256,16); // hybrid tracks for LHC11a, LHC12a15a and LHC12a15e + } + else { + if (runPeriod.IsNull()) + ::Warning("Run period %s not known. It will use IsHybridGlobalConstrainedGlobal."); + } + + //------------------------------------------------------- + // Final settings, pass to manager and set the containers + //------------------------------------------------------- + + mgr->AddTask(jetEmb); + + // Create containers for input/output + mgr->ConnectInput(jetEmb, 0, mgr->GetCommonInputContainer()); + + if (makeQA) { + TString contName = taskName; + contName += "_histos"; + AliAnalysisDataContainer *outc = mgr->CreateContainer(contName, + TList::Class(), + AliAnalysisManager::kOutputContainer, + "AnalysisResults.root"); + mgr->ConnectOutput(jetEmb, 1, outc); + } + + return jetEmb; +} + +TObjArray* GenerateFileList(const char* list, Int_t nFiles) +{ + TObjArray *array = new TObjArray(9999); + + TString myList = list; + if (myList.Contains("alien:///")) { + TFile::Cp(myList,"file:./list.txt"); + myList = "./list.txt"; + } + + // Open the input stream + ifstream in; + in.open(myList.Data()); + + Int_t count = 0; + + // Read the input list of files and add them to the chain + TString line; + while (in.good()) { + if (nFiles != 1234567890) { + if (count >= nFiles) + break; + } + + in >> line; + + if (line.Length() == 0) + continue; + + TObjString *aodFile = new TObjString(line); + array->Add(aodFile); + + count++; + } + + return array; +} diff --git a/PWGJE/EMCALJetTasks/macros/AddTaskJetRespPtHard.C b/PWGJE/EMCALJetTasks/macros/AddTaskJetRespPtHard.C index 5123657c5d1..4adf5156eef 100644 --- a/PWGJE/EMCALJetTasks/macros/AddTaskJetRespPtHard.C +++ b/PWGJE/EMCALJetTasks/macros/AddTaskJetRespPtHard.C @@ -3,6 +3,7 @@ AliJetResponseMaker* AddTaskJetRespPtHard(const char *ntracks = "Tracks", const char *nclusters = "CaloClusters", const char *njets = "Jets", + const char *nrho = "Rho", const char *nmctracks = "MCParticles", const char *nmcjets = "MCJets", Double_t jetradius = 0.2, @@ -25,7 +26,7 @@ AliJetResponseMaker* AddTaskJetRespPtHard(const char *ntracks = "Trac AliJetResponseMaker *jetTask = new AliJetResponseMaker[maxPtBin - minPtBin + 1]; for (Int_t i = minPtBin; i <= maxPtBin; i++) { - AddTaskJetResponseMaker(ntracks, nclusters, njets, nmctracks, nmcjets, + AddTaskJetResponseMaker(ntracks, nclusters, njets, nrho, nmctracks, nmcjets, jetradius, jetptcut, jetareacut, jetBiasTrack, jetBiasClus, maxDistance, type, i, taskname, jetTask + i - minPtBin); jetTask[i - minPtBin].SetDoMatching(domatch); diff --git a/PWGJE/EMCALJetTasks/macros/AddTaskJetResponseMaker.C b/PWGJE/EMCALJetTasks/macros/AddTaskJetResponseMaker.C index 005f58da37a..d0724f6dafa 100644 --- a/PWGJE/EMCALJetTasks/macros/AddTaskJetResponseMaker.C +++ b/PWGJE/EMCALJetTasks/macros/AddTaskJetResponseMaker.C @@ -4,6 +4,7 @@ AliJetResponseMaker* AddTaskJetResponseMaker( const char *ntracks = "Tracks", const char *nclusters = "CaloClusters", const char *njets = "Jets", + const char *nrho = "Rho", const char *nmctracks = "MCParticles", const char *nmcjets = "MCJets", Double_t jetradius = 0.2, @@ -61,6 +62,7 @@ AliJetResponseMaker* AddTaskJetResponseMaker( jetTask->SetTracksName(ntracks); jetTask->SetClusName(nclusters); jetTask->SetJetsName(njets); + jetTask->SetRhoName(nrho); jetTask->SetMCJetsName(nmcjets); jetTask->SetMCTracksName(nmctracks); jetTask->SetJetRadius(jetradius); diff --git a/PWGJE/PWGJEEMCALJetTasksLinkDef.h b/PWGJE/PWGJEEMCALJetTasksLinkDef.h index 9cbd8884f02..210d81f127c 100644 --- a/PWGJE/PWGJEEMCALJetTasksLinkDef.h +++ b/PWGJE/PWGJEEMCALJetTasksLinkDef.h @@ -18,6 +18,7 @@ #pragma link C++ class AliHadCorrTask+; #pragma link C++ class AliJetEmbeddingTask+; #pragma link C++ class AliJetEmbeddingFromGenTask+; +#pragma link C++ class AliJetEmbeddingFromAODTask+; #pragma link C++ class AliJetModelBaseTask+; #pragma link C++ class AliJetRandomizerTask+; #pragma link C++ class AliJetResponseMaker+; -- 2.39.3