adding task for subtracting background after jet finding, used for all clustering...
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 20 Dec 2010 07:08:29 +0000 (07:08 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 20 Dec 2010 07:08:29 +0000 (07:08 +0000)
JETAN/AliAnalysisTaskJetBackgroundSubtract.cxx [new file with mode: 0644]
JETAN/AliAnalysisTaskJetBackgroundSubtract.h [new file with mode: 0644]
JETAN/AliAnalysisTaskJetCluster.cxx
JETAN/CMakelibFASTJETAN.pkg
JETAN/FASTJETANLinkDef.h
JETAN/libFASTJETAN.pkg

diff --git a/JETAN/AliAnalysisTaskJetBackgroundSubtract.cxx b/JETAN/AliAnalysisTaskJetBackgroundSubtract.cxx
new file mode 100644 (file)
index 0000000..1775f4d
--- /dev/null
@@ -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 <TROOT.h>
+#include <TH1F.h>
+#include <THnSparse.h>
+#include <TSystem.h>
+#include <TInterpreter.h>
+#include <TList.h>
+#include <TLorentzVector.h>
+#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<AliAODHandler*>(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<AliAODEvent*>(InputEvent());
+
+
+  // collect the jet branches 
+
+  if(!fInJetArrayList)fInJetArrayList =new TList();
+  if(!fOutJetArrayList)fOutJetArrayList =new TList();
+
+  for(int iJB = 0;iJB<fJBArray->GetEntries();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; i<fHistList->GetEntries(); ++i) {
+    TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
+    if (h1){
+      h1->Sumw2();
+      continue;
+    }
+    THnSparse *hn = dynamic_cast<THnSparse*>(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;iJB<fInJetArrayList->GetEntries();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;iJB<fOutJetArrayList->GetEntries();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 (file)
index 0000000..ca50720
--- /dev/null
@@ -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
index 5ba6a3e..0c0f1dd 100644 (file)
@@ -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<fastjet::PseudoJet> 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){       
index 1322701..865f78f 100644 (file)
@@ -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})
 
index 4497b75..a334785 100644 (file)
@@ -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
index 3a0e627..0b5ab56 100644 (file)
@@ -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