From fb10c4b8e46003aa13c169e9886ceaa35fa3b470 Mon Sep 17 00:00:00 2001 From: kleinb Date: Mon, 20 Dec 2010 07:08:29 +0000 Subject: [PATCH] adding task for subtracting background after jet finding, used for all clustering jet finders --- .../AliAnalysisTaskJetBackgroundSubtract.cxx | 302 ++++++++++++++++++ JETAN/AliAnalysisTaskJetBackgroundSubtract.h | 89 ++++++ JETAN/AliAnalysisTaskJetCluster.cxx | 15 +- JETAN/CMakelibFASTJETAN.pkg | 2 +- JETAN/FASTJETANLinkDef.h | 1 + JETAN/libFASTJETAN.pkg | 2 +- 6 files changed, 403 insertions(+), 8 deletions(-) create mode 100644 JETAN/AliAnalysisTaskJetBackgroundSubtract.cxx create mode 100644 JETAN/AliAnalysisTaskJetBackgroundSubtract.h diff --git a/JETAN/AliAnalysisTaskJetBackgroundSubtract.cxx b/JETAN/AliAnalysisTaskJetBackgroundSubtract.cxx new file mode 100644 index 00000000000..1775f4ddd80 --- /dev/null +++ b/JETAN/AliAnalysisTaskJetBackgroundSubtract.cxx @@ -0,0 +1,302 @@ +// ************************************** +// Task used for the correction of determiantion of reconstructed jet spectra +// Compares input (gen) and output (rec) jets +// ******************************************* + + +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + + +#include +#include +#include +#include +#include +#include +#include +#include "TDatabasePDG.h" + +#include "AliAnalysisTaskJetBackgroundSubtract.h" +#include "AliAnalysisManager.h" +#include "AliAODHandler.h" +#include "AliAODTrack.h" +#include "AliAODJet.h" +#include "AliAODEvent.h" +#include "AliInputEventHandler.h" +#include "AliAODJetEventBackground.h" + + +ClassImp(AliAnalysisTaskJetBackgroundSubtract) + +AliAnalysisTaskJetBackgroundSubtract::~AliAnalysisTaskJetBackgroundSubtract(){ + delete fJBArray; + delete fOutJetArrayList; + delete fInJetArrayList; +} + +AliAnalysisTaskJetBackgroundSubtract::AliAnalysisTaskJetBackgroundSubtract(): + AliAnalysisTaskSE(), + fAODOut(0x0), + fAODIn(0x0), + fAODExtension(0x0), + fJBArray(new TObjArray()), + fBackgroundBranch(""), + fNonStdFile(""), + fReplaceString1("B0"), + fReplaceString2("B%d"), + fSubtraction(kArea), + fInJetArrayList(0x0), + fOutJetArrayList(0x0), + fHistList(0x0) +{ + +} + +AliAnalysisTaskJetBackgroundSubtract::AliAnalysisTaskJetBackgroundSubtract(const char* name): + + AliAnalysisTaskSE(name), + fAODOut(0x0), + fAODIn(0x0), + fAODExtension(0x0), + fJBArray(new TObjArray()), + fBackgroundBranch(""), + fNonStdFile(""), + fReplaceString1("B0"), + fReplaceString2("B%d"), + fSubtraction(kArea), + fInJetArrayList(0x0), + fOutJetArrayList(0x0), + fHistList(0x0) +{ + DefineOutput(1, TList::Class()); +} + + + +Bool_t AliAnalysisTaskJetBackgroundSubtract::Notify() +{ + // + return kTRUE; +} + +void AliAnalysisTaskJetBackgroundSubtract::UserCreateOutputObjects() +{ + + // + // Create the output container + // + // Connect the AOD + + + if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::UserCreateOutputObjects() \n"); + if(fNonStdFile.Length()!=0){ + // + // case that we have an AOD extension we need to fetch the jets from the extended output + // we identifay the extension aod event by looking for the branchname + AliAODHandler *aodH = dynamic_cast(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()); + fAODExtension = aodH->GetExtension(fNonStdFile.Data()); + + if(!fAODExtension){ + if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data()); + } + } + fAODOut = AODEvent(); + fAODIn = dynamic_cast(InputEvent()); + + + // collect the jet branches + + if(!fInJetArrayList)fInJetArrayList =new TList(); + if(!fOutJetArrayList)fOutJetArrayList =new TList(); + + for(int iJB = 0;iJBGetEntries();iJB++){ + TClonesArray* jarray = 0; + TObjString *ostr = (TObjString*)fJBArray->At(iJB); + if(!jarray&&fAODOut){ + jarray = (TClonesArray*)(fAODOut->FindListObject(ostr->GetString().Data())); + } + if(!jarray&&fAODExtension){ + jarray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(ostr->GetString().Data())); + } + if(!jarray){ + if(fDebug)Printf("%s:%d jet branch %s not found",(char*)__FILE__,__LINE__,ostr->GetString().Data()); + continue; + } + + TString newName(jarray->GetName()); + if(!newName.Contains(fReplaceString1.Data())){ + Printf("%s:%d cannot replace string %s in %s",(char*)__FILE__,__LINE__,fReplaceString1.Data(), + jarray->GetName()); + continue; + } + + fInJetArrayList->Add(jarray); + + // add a new branch to the output for the background subtracted jets take the names from + // the input jets and replace the background flag names + TClonesArray *tca = new TClonesArray("AliAODJet", 0); + + newName.ReplaceAll(fReplaceString1.Data(),Form(fReplaceString2.Data(),fSubtraction)); + if(fDebug){ + Printf("%s:%d created branch \n %s from \n %s",(char*)__FILE__,__LINE__,newName.Data(), + jarray->GetName()); + } + tca->SetName(newName.Data()); + AddAODBranch("TClonesArray",&tca,fNonStdFile.Data()); + fOutJetArrayList->Add(tca); + } + + + + + if(!fHistList)fHistList = new TList(); + fHistList->SetOwner(); + + Bool_t oldStatus = TH1::AddDirectoryStatus(); + TH1::AddDirectory(kFALSE); + + // + // Histogram booking, add som control histograms here + // + + // =========== Switch on Sumw2 for all histos =========== + for (Int_t i=0; iGetEntries(); ++i) { + TH1 *h1 = dynamic_cast(fHistList->At(i)); + if (h1){ + h1->Sumw2(); + continue; + } + THnSparse *hn = dynamic_cast(fHistList->At(i)); + if(hn)hn->Sumw2(); + } + TH1::AddDirectory(oldStatus); + + if(fBackgroundBranch.Length()==0) + AliError(Form("%s:%d No BackgroundBranch defined",(char*)__FILE__,__LINE__)); + if(fJBArray->GetEntries()==0) + AliError(Form("%s:%d No Jet Branches defined defined",(char*)__FILE__,__LINE__)); +} + +void AliAnalysisTaskJetBackgroundSubtract::Init() +{ + // + // Initialization + // + if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::Init() \n"); +} + +void AliAnalysisTaskJetBackgroundSubtract::UserExec(Option_t */*option*/) +{ + + ResetOutJets(); + if(fBackgroundBranch.Length()==0||fJBArray->GetEntries()==0){ + PostData(1,fHistList); + } + if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::UserExec() \n"); + + + static AliAODJetEventBackground* evBkg = 0; + if(!evBkg&&fAODOut){ + evBkg = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data())); + } + if(!evBkg&&fAODExtension){ + evBkg = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data())); + } + if(!evBkg&&fAODIn){ + evBkg = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data())); + } + + if(!evBkg){ + if(fDebug)Printf("%s:%d Backroundbranch %s not found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data()); + PostData(1,fHistList); + } + + + // LOOP over all jet branches and subtract the background + + Float_t rho = 0; + if(fSubtraction==kArea)rho = evBkg->GetBackground(2); + + for(int iJB = 0;iJBGetEntries();iJB++){ + TClonesArray* jarray = (TClonesArray*)fInJetArrayList->At(iJB); + TClonesArray* jarrayOut = (TClonesArray*)fOutJetArrayList->At(iJB); + if(!jarray||!jarrayOut)continue; + + // loop over all jets + Int_t nOut = 0; + for(int i = 0;i < jarray->GetEntriesFast();i++){ + AliAODJet *jet = (AliAODJet*)jarray->At(i); + AliAODJet tmpNewJet(*jet); + Bool_t bAdd = false; + if(fSubtraction==kArea){ + Float_t ptSub = jet->Pt() - rho * jet->EffectiveAreaCharged(); + if(fDebug>2){ + Printf("%s:%d Jet %d %3.3f %3.3f",(char*)__FILE__,__LINE__,i,jet->Pt(),ptSub); + } + if(ptSub<0){ + // optionally rescale it and keep?? + bAdd = RescaleJetMomentum(&tmpNewJet,0.1); + } + else{ + bAdd = RescaleJetMomentum(&tmpNewJet,ptSub); + } + }// kAREA + + if(bAdd){ + new ((*jarrayOut)[nOut++]) AliAODJet(tmpNewJet); + // optionally add background estimates to the jet object? CKB + } + + } + + + + // subtract the background + + + // remove jets?? + + // sort jets... + + } + PostData(1, fHistList); +} + +void AliAnalysisTaskJetBackgroundSubtract::Terminate(Option_t */*option*/) +{ + // Terminate analysis + // + if (fDebug > 1) printf("AnalysisJetBackgroundSubtract: Terminate() \n"); +} + +Bool_t AliAnalysisTaskJetBackgroundSubtract::RescaleJetMomentum(AliAODJet *jet,Float_t pT){ + // keep the direction and the jet mass + if(pT<=0)return kFALSE; + Double_t pTold = jet->Pt(); + Double_t scale = pT/pTold; + Double_t mass = jet->M(); + Double_t pNew = jet->P() * scale; + jet->SetPxPyPzE(scale*jet->Px(),scale*jet->Py(),scale*jet->Pz(),TMath::Sqrt(mass*mass+pNew*pNew)); + return kTRUE; +} + +void AliAnalysisTaskJetBackgroundSubtract::ResetOutJets(){ + if(!fOutJetArrayList)return; + for(int iJB = 0;iJBGetEntries();iJB++){ + TClonesArray* jarray = (TClonesArray*)fOutJetArrayList->At(iJB); + if(jarray)jarray->Delete(); + } +} diff --git a/JETAN/AliAnalysisTaskJetBackgroundSubtract.h b/JETAN/AliAnalysisTaskJetBackgroundSubtract.h new file mode 100644 index 00000000000..ca507206b23 --- /dev/null +++ b/JETAN/AliAnalysisTaskJetBackgroundSubtract.h @@ -0,0 +1,89 @@ +#ifndef ALIANALYSISTASKJETBACKGROUNDSUBTRACT_H +#define ALIANALYSISTASKJETBACKGROUNDSUBTRACT_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +// ************************************** +// task used for background subtration of all jets found with clustering algos +// ******************************************* + +#include "AliAnalysisTaskSE.h" +#include "TObjString.h" +#include "TString.h" +#include "TObjArray.h" + + +//////////////// +class AliJetHeader; +class AliESDEvent; +class AliAODEvent; +class AliAODExtension; +class AliAODJet; +class AliAODJetEventBackground; +class AliJetFinder; +class TList; +class TChain; +class TH2F; +class TH1F; +class TH3F; +class TProfile; +class TRandom3; +class TRefArray; + + +class AliAnalysisTaskJetBackgroundSubtract : public AliAnalysisTaskSE +{ + public: + AliAnalysisTaskJetBackgroundSubtract(); + AliAnalysisTaskJetBackgroundSubtract(const char* name); + virtual ~AliAnalysisTaskJetBackgroundSubtract(); + // Implementation of interface methods + virtual void UserCreateOutputObjects(); + virtual void Init(); + virtual void LocalInit() { Init(); } + virtual void UserExec(Option_t *option); + virtual void Terminate(Option_t *option); + virtual Bool_t Notify(); + + // Task specific methods... + virtual void AddJetBranch(const char* c){fJBArray->Add(new TObjString(c));} + virtual void SetSubtractionMethod(Int_t i){fSubtraction = i;} + virtual Int_t GetSubtractionMethod(){return fSubtraction;} + virtual void SetBackgroundBranch(char* c){fBackgroundBranch = c;} + virtual void SetNonStdFile(char* c){fNonStdFile = c;} + virtual void SetToReplace(char* c){fReplaceString1 = c;} + const char* GetToReplace(){return fReplaceString1.Data();} + virtual void SetReplacementMask(char* c){fReplaceString2 = c;} + const char* GetReplacementMask(){fReplaceString2.Data();} + + enum {kNoSubtract = 0,kArea,k4Area}; + + private: + + + + AliAnalysisTaskJetBackgroundSubtract(const AliAnalysisTaskJetBackgroundSubtract&); + AliAnalysisTaskJetBackgroundSubtract& operator=(const AliAnalysisTaskJetBackgroundSubtract&); + Bool_t RescaleJetMomentum(AliAODJet *jet,Float_t pT); + void ResetOutJets(); + + + AliAODEvent *fAODOut; // ! where we take the jets from and they are modified + AliAODEvent *fAODIn; // ! where we may take the background from, only in case we do not find it in the output + AliAODExtension *fAODExtension; // ! where we take the jets from can be input or output AOD + TObjArray *fJBArray; // Array that stores the name of all jet branches to be subtracted + TString fBackgroundBranch; // name of the branch used for background subtraction + // + TString fNonStdFile; // The optional name of the output file the non-std brnach is written to + TString fReplaceString1; // To construct the new output name + TString fReplaceString2; // To construct the new output name + Int_t fSubtraction; // Parameter for subtraction mode + TList *fInJetArrayList; //! transient list to make ease the handling of input jets + TList *fOutJetArrayList; //! transient list to make ease the reset of output jets + TList *fHistList; //! the histograms output list + + ClassDef(AliAnalysisTaskJetBackgroundSubtract, 1) +}; + +#endif diff --git a/JETAN/AliAnalysisTaskJetCluster.cxx b/JETAN/AliAnalysisTaskJetCluster.cxx index 5ba6a3efc31..0c0f1dd37d0 100644 --- a/JETAN/AliAnalysisTaskJetCluster.cxx +++ b/JETAN/AliAnalysisTaskJetCluster.cxx @@ -904,6 +904,14 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) aodOutJet = 0; nAodOutTracks = 0; Float_t tmpPt = tmpRec.Pt(); + + if(tmpPt>fJetOutputMinPt&&jarray){// cut on the non-background subtracted... + aodOutJet = new ((*jarray)[nAodOutJets++]) AliAODJet(tmpRec); + Double_t area1 = clustSeq.area(sortedJets[j]); + aodOutJet->SetEffArea(area1,0); + } + + Float_t tmpPtBack = 0; if(externalBackground){ // carefull has to be filled in a task before @@ -916,6 +924,7 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) fh1PtJetsRecIn->Fill(tmpPt); // Fill Spectra with constituents vector constituents = clustSeq.constituents(sortedJets[j]); + fh1NConstRec->Fill(constituents.size()); fh2PtNch->Fill(nCh,tmpPt); fh2PtNchN->Fill(nCh,tmpPt,constituents.size()); @@ -923,12 +932,6 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/) // loop over constiutents and fill spectrum // Add the jet information and the track references to the output AOD - - if(tmpPt>fJetOutputMinPt&&jarray){ - aodOutJet = new ((*jarray)[nAodOutJets++]) AliAODJet(tmpRec); - Double_t area1 = clustSeq.area(sortedJets[j]); - aodOutJet->SetEffArea(area1,0); - } if(fUseBackgroundCalc){ diff --git a/JETAN/CMakelibFASTJETAN.pkg b/JETAN/CMakelibFASTJETAN.pkg index 1322701fc5e..865f78fccf4 100644 --- a/JETAN/CMakelibFASTJETAN.pkg +++ b/JETAN/CMakelibFASTJETAN.pkg @@ -25,7 +25,7 @@ # SHLIBS - Shared Libraries and objects for linking (Executables only) # #--------------------------------------------------------------------------------# -set ( SRCS AliFastJetFinder.cxx AliFastJetHeaderV1.cxx AliFastJetInput.cxx AliJetBkg.cxx AliSISConeJetFinder.cxx AliSISConeJetHeader.cxx AliAnalysisTaskJetCluster.cxx) +set ( SRCS AliFastJetFinder.cxx AliFastJetHeaderV1.cxx AliFastJetInput.cxx AliJetBkg.cxx AliSISConeJetFinder.cxx AliSISConeJetHeader.cxx AliAnalysisTaskJetCluster.cxx AliAnalysisTaskJetBackgroundSubtract.cxx) set(FASTJET_ROOT $ENV{FASTJET_ROOT}) diff --git a/JETAN/FASTJETANLinkDef.h b/JETAN/FASTJETANLinkDef.h index 4497b75bf88..a3347851dc5 100644 --- a/JETAN/FASTJETANLinkDef.h +++ b/JETAN/FASTJETANLinkDef.h @@ -11,4 +11,5 @@ #pragma link C++ class AliSISConeJetFinder+; #pragma link C++ class AliSISConeJetHeader+; #pragma link C++ class AliAnalysisTaskJetCluster+; +#pragma link C++ class AliAnalysisTaskJetBackgroundSubtract+; #endif diff --git a/JETAN/libFASTJETAN.pkg b/JETAN/libFASTJETAN.pkg index 3a0e6277b80..0b5ab56a638 100644 --- a/JETAN/libFASTJETAN.pkg +++ b/JETAN/libFASTJETAN.pkg @@ -1,7 +1,7 @@ #-*-Mode: Makefile-*- SRCS = AliFastJetFinder.cxx AliFastJetHeaderV1.cxx AliFastJetInput.cxx AliJetBkg.cxx\ - AliSISConeJetFinder.cxx AliSISConeJetHeader.cxx AliAnalysisTaskJetCluster.cxx + AliSISConeJetFinder.cxx AliSISConeJetHeader.cxx AliAnalysisTaskJetCluster.cxx AliAnalysisTaskJetBackgroundSubtract.cxx ifneq ($(FASTJET_ROOT),) EDEFINE =-isystem$(FASTJET_ROOT)/include -- 2.39.3