[u/mrichter/AliRoot.git] / JETAN / AliAnalysisTaskJetBackgroundSubtract.cxx

// **************************************
// 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 <TH2F.h>
#include <THnSparse.h>
#include <TSystem.h>
#include <TInterpreter.h>
#include <TList.h>
#include <TLorentzVector.h>
#include <TRefArray.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(kRhoRecalc),
  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(kRhoRecalc),
  fInJetArrayList(0x0),
  fOutJetArrayList(0x0),
  fHistList(0x0)  
{
 DefineOutput(1, TList::Class());  
}


Bool_t AliAnalysisTaskJetBackgroundSubtract::Notify()
{
  //
  fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());

  ResetOutJets();

  // Now we also have the Input Event Available! Fillvthe pointers in the list

  fInJetArrayList->Clear();
  fOutJetArrayList->Clear();

  for(int iJB = 0;iJB<fJBArray->GetEntries();iJB++){
    TObjString *ostr = (TObjString*)fJBArray->At(iJB);
 
  
    TClonesArray* jarray = 0;      
    if(!jarray&&fAODOut){
      jarray = (TClonesArray*)(fAODOut->FindListObject(ostr->GetString().Data()));
    }
    if(!jarray&&fAODExtension){
      jarray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(ostr->GetString().Data()));
    }
    if(!jarray&&fAODIn){
      jarray = (TClonesArray*)(fAODIn->FindListObject(ostr->GetString().Data()));
    }

    if(!jarray){
      if(fDebug){
	Printf("%s:%d Input jet branch %s not found",(char*)__FILE__,__LINE__,ostr->GetString().Data());
      }
      continue;
    }
    
    TString newName(ostr->GetString().Data());
    newName.ReplaceAll(fReplaceString1.Data(),Form(fReplaceString2.Data(),fSubtraction));
    TClonesArray* jarrayOut = 0;      
    if(!jarrayOut&&fAODOut){
      jarrayOut = (TClonesArray*)(fAODOut->FindListObject(newName.Data()));
    }
    if(!jarrayOut&&fAODExtension){
      jarrayOut = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(newName.Data()));
    }

    if(!jarrayOut){
      if(fDebug){
	Printf("%s:%d Output jet branch %s not found",(char*)__FILE__,__LINE__,newName.Data());
	PrintAODContents();
      }
      continue;
    }
    if(jarrayOut&&jarray){
      fOutJetArrayList->Add(jarrayOut);
      fInJetArrayList->Add(jarray);
    }
  }
  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();

  // usually we do not have the input already here

  if(!fInJetArrayList)fInJetArrayList =new TList();
  if(!fOutJetArrayList)fOutJetArrayList =new TList();

  for(int iJB = 0;iJB<fJBArray->GetEntries();iJB++){
    TObjString *ostr = (TObjString*)fJBArray->At(iJB);
    TString newName(ostr->GetString().Data());
    if(!newName.Contains(fReplaceString1.Data())){
      Printf("%s:%d cannot replace string %s in %s",(char*)__FILE__,__LINE__,fReplaceString1.Data(),
	     newName.Data());
      continue;
    }

    // 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(),
	     ostr->GetString().Data());
    }
    tca->SetName(newName.Data());
    AddAODBranch("TClonesArray",&tca,fNonStdFile.Data());
  }
  

  if(!fHistList)fHistList = new TList();
  fHistList->SetOwner();

  for(int iJB = 0;iJB<fJBArray->GetEntries();iJB++){
    TObjString *ostr = (TObjString*)fJBArray->At(iJB);
    TString oldName(ostr->GetString().Data()); 
    TString newName(ostr->GetString().Data()); 
    if(!newName.Contains(fReplaceString1.Data())){
      Printf("%s:%d cannot replace string %s in %s",(char*)__FILE__,__LINE__,fReplaceString1.Data(),
	     newName.Data());
      continue;
    }
    newName.ReplaceAll(fReplaceString1.Data(),Form(fReplaceString2.Data(),fSubtraction));
    TH2F *hTmp = new TH2F(Form("h2PtInPtOut_%d",iJB),Form(";%s p_{T}; %s p_{T}",oldName.Data(),newName.Data()),200,0,200.,200,0.,200.);
    fHistList->Add(hTmp);
  }

  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*/)
{

  if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::UserExec() \n");
  ResetOutJets();
  if(fBackgroundBranch.Length()==0||fJBArray->GetEntries()==0){
    if(fDebug)Printf("%s:%d No background subtraction done",(char*)__FILE__,__LINE__);
    PostData(1,fHistList);
  }
  if(fJBArray->GetEntries()!=fInJetArrayList->GetEntries()){
    if(fDebug)Printf("%s:%d different Array  sizes %d %d %d",(char*)__FILE__,__LINE__,fJBArray->GetEntries(),fInJetArrayList->GetEntries(),fOutJetArrayList->GetEntries());
    PostData(1,fHistList);
  }


  AliAODJetEventBackground*  evBkg = 0;
  TClonesArray*              bkgClusters = 0;
  TClonesArray*              bkgClustersRC = 0;
  TString bkgClusterName(fBackgroundBranch.Data());
  bkgClusterName.ReplaceAll(Form("%s_",AliAODJetEventBackground::StdBranchName()),"");
  TString bkgClusterRCName(Form("%s%s",bkgClusterName.Data(),"RandomCone")); 

  if(!evBkg&&!bkgClusters&&!bkgClustersRC&&fAODOut){
    evBkg = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
    bkgClusters = (TClonesArray*)(fAODOut->FindListObject(bkgClusterName.Data()));
    bkgClustersRC = (TClonesArray*)(fAODOut->FindListObject(bkgClusterRCName.Data()));

    if(fDebug&&bkgClusters)Printf("%s:%d Background cluster branch %s found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
    if(fDebug&&bkgClustersRC)Printf("%s:%d Background cluster RC branch %s found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
    if(fDebug&&evBkg)Printf("%s:%d Backgroundbranch %s found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
  }
  if(!evBkg&&!bkgClusters&&!bkgClustersRC&&fAODExtension){
    evBkg = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
    bkgClusters = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(bkgClusterName.Data()));
    bkgClustersRC = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(bkgClusterRCName.Data()));
    if(fDebug&&bkgClusters)Printf("%s:%d Background cluster branch %s found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
    if(fDebug&&bkgClustersRC)Printf("%s:%d Background cluster RC branch %s found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());

    if(fDebug&&evBkg)Printf("%s:%d Backgroundbranch %s found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
  }

  if(!evBkg&&!bkgClusters&&!bkgClustersRC&&fAODIn){
    evBkg = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data()));
    bkgClusters = (TClonesArray*)(fAODIn->FindListObject(bkgClusterName.Data()));
    bkgClustersRC = (TClonesArray*)(fAODIn->FindListObject(bkgClusterRCName.Data()));

    if(fDebug&&bkgClusters)Printf("%s:%d Background cluster branch %s found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
    if(fDebug&&bkgClustersRC)Printf("%s:%d Background cluster RC branch %s found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
    if(fDebug&&evBkg)Printf("%s:%d Backgroundbranch %s found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
  }

  if(!evBkg&&(fSubtraction==kArea||fSubtraction==kRhoRecalc)){
    if(fDebug){
      Printf("%s:%d Backgroundbranch %s not found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
      PrintAODContents();
    }
    PostData(1,fHistList);
    return;
  }

  if(!bkgClusters&&(fSubtraction==kRhoRecalc)){
    if(fDebug){
      Printf("%s:%d Background cluster branch %s not found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
      PrintAODContents();
    }
    PostData(1,fHistList);
    return;
  }

  if(!bkgClustersRC&&(fSubtraction==kRhoRC)){
    if(fDebug){
      Printf("%s:%d Background cluster RC branch %s not found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
      PrintAODContents();
    }
    PostData(1,fHistList);
    return;
  }
  // LOOP over all jet branches and subtract the background

   Float_t rho = 0;
   Double_t meanarea = 0;
   if(fSubtraction==kArea)rho =  evBkg->GetBackground(2);
   if(fSubtraction==kRhoRecalc){
     meanarea=evBkg->GetMeanarea(2);
     rho =RecalcRho(bkgClusters,meanarea);
   }
   if(fSubtraction==kRhoRC)rho =RhoRC(bkgClustersRC);
   
  for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
    TClonesArray* jarray = (TClonesArray*)fInJetArrayList->At(iJB);
    TClonesArray* jarrayOut = (TClonesArray*)fOutJetArrayList->At(iJB);
    
    if(!jarray||!jarrayOut){
    Printf("%s:%d Array not found %d: %p %p",(char*)__FILE__,__LINE__,iJB,jarray->GetName(),jarrayOut->GetName());
    continue;
      }
    TH2F* h2PtInOut = (TH2F*)fHistList->FindObject(Form("h2PtInPtOut_%d",iJB));
    // 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){	
	Double_t background = rho * jet->EffectiveAreaCharged();
	Float_t ptSub = jet->Pt() - background;	
	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);
	  if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),0.1);
	}
	else{
	  bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
	  if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
	}
	// add background estimates to the new jet object
	// allows to recover old p_T and rho...
	tmpNewJet.SetBgEnergy(background,0);
      }// kAREA
      else if(fSubtraction==kRhoRecalc){
 	Double_t background = rho * jet->EffectiveAreaCharged();
	Float_t ptSub = jet->Pt() - background;	
        if(fDebug>2){
	  Printf("%s:%d Jet %d %3.3f %3.3f %3.3f %3.3f",(char*)__FILE__,__LINE__,i,jet->Pt(),ptSub,background,rho);}
	if(ptSub<0){
	  // optionally rescale it and keep??
	  bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
	  if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),0.1);
	}
	else{
	  bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
	  if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
	}
	// add background estimates to the new jet object
	// allows to recover old p_T and rho...
	tmpNewJet.SetBgEnergy(background,0);

      }//kRhoRecalc
       else if(fSubtraction==kRhoRC){
	Double_t background = rho * jet->EffectiveAreaCharged();
	Float_t ptSub = jet->Pt() - background;	
	if(fDebug>2){	Printf("%s:%d Jet %d %3.3f %3.3f %3.3f %3.3f",(char*)__FILE__,__LINE__,i,jet->Pt(),ptSub,background,rho);}
	if(ptSub<0){
	  // optionally rescale it and keep??
	  bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
	  if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),0.1);
	}
	else{
	  bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
	  if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
	}
	// add background estimates to the new jet object
	// allows to recover old p_T and rho...
	tmpNewJet.SetBgEnergy(background,0);

       }//kRhoRC

      if(bAdd){
        AliAODJet *newJet = new ((*jarrayOut)[nOut++]) AliAODJet(tmpNewJet);
	// what about track references, clear for now...
	newJet->GetRefTracks()->Clear();
      }


    }


    // 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;
}

Double_t AliAnalysisTaskJetBackgroundSubtract::RecalcRho(TClonesArray* bkgClusters,Double_t meanarea){
  
       Double_t ptarea=0.;
       Int_t count=0;
       Double_t rho=0.; 
       const Double_t Rlimit2=0.8*0.8;  //2*jet radius.
       TClonesArray* jarray=0;
       
        for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
        TObjString *ostr = (TObjString*)fInJetArrayList->At(iJB);
        TString jetref=ostr->GetString().Data();
        if(jetref.Contains("ANTIKT04")){ 
	  jarray = (TClonesArray*)fInJetArrayList->At(iJB);}}
	if(jarray->GetEntries()>=2){ 
	  AliAODJet *first = (AliAODJet*)(jarray->At(0)); 
	  AliAODJet *second= (AliAODJet*)(jarray->At(1)); 
	  for(Int_t k=0;k<bkgClusters->GetEntriesFast();k++){
	    AliAODJet *clus = (AliAODJet*)(bkgClusters->At(k));
	    if(TMath::Abs(clus->Eta())>0.5) continue;
	    if((clus->EffectiveAreaCharged())<0.1*meanarea) continue; 
	    Double_t distance1=(first->Eta()-clus->Eta())*(first->Eta()-clus->Eta())+
	      (first->Phi()-clus->Phi())*(first->Phi()-clus->Phi());
	    Double_t distance2= (second->Eta()-clus->Eta())*(second->Eta()-clus->Eta())+
	      (second->Phi()-clus->Phi())*(second->Phi()-clus->Phi());
	    if((distance1<Rlimit2)||(distance2<Rlimit2)) continue;    
	    ptarea=ptarea+clus->Pt()/clus->EffectiveAreaCharged(); 
	    count=count+1;}
	  if(count!=0) rho=ptarea/count; 
	}        
	return rho;
}

   Double_t AliAnalysisTaskJetBackgroundSubtract::RhoRC(TClonesArray* bkgClustersRC){
  
       Double_t ptarea=0.;
       Int_t count=0;
       Double_t rho=0.; 
       const Double_t Rlimit2=0.8*0.8;  //2*jet radius.
       TClonesArray* jarray=0;
        for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
        TObjString *ostr = (TObjString*)fInJetArrayList->At(iJB);
        TString jetref=ostr->GetString().Data();
        if(jetref.Contains("ANTIKT04")){ 
	 jarray = (TClonesArray*)fInJetArrayList->At(iJB);}}
         if(jarray->GetEntries()>=2){ 
	   AliAODJet *first = (AliAODJet*)(jarray->At(0)); 
	   AliAODJet *second=(AliAODJet*)(jarray->At(1)); 
         for(Int_t k=0;k<bkgClustersRC->GetEntriesFast();k++){
	   AliAODJet *clus = (AliAODJet*)(bkgClustersRC->At(k));
	   if(TMath::Abs(clus->Eta())>0.5) continue;
	   Double_t distance1=(first->Eta()-clus->Eta())*(first->Eta()-clus->Eta())+
	     (first->Phi()-clus->Phi())*(first->Phi()-clus->Phi());
	   Double_t distance2= (second->Eta()-clus->Eta())*(second->Eta()-clus->Eta())+
	     (second->Phi()-clus->Phi())*(second->Phi()-clus->Phi());
	   if((distance1<Rlimit2)||(distance2<Rlimit2)) continue;    
	   ptarea=ptarea+clus->Pt()/clus->EffectiveAreaCharged(); 
	   count=count+1;}
         if(count!=0) rho=ptarea/count;  }
         return rho;
}


void AliAnalysisTaskJetBackgroundSubtract::ResetOutJets(){
  if(!fOutJetArrayList)return;
  for(int iJB = 0;iJB<fOutJetArrayList->GetEntries();iJB++){
    TClonesArray* jarray = (TClonesArray*)fOutJetArrayList->At(iJB);
    if(jarray)jarray->Delete();
  }
}


void AliAnalysisTaskJetBackgroundSubtract::PrintAODContents(){
  if(fAODIn){
    Printf("%s:%d >>>>>> Input",(char*)__FILE__,__LINE__);
    fAODIn->Print();
  }
  if(fAODExtension){
    Printf("%s:%d >>>>>> Extenstion",(char*)__FILE__,__LINE__);
    fAODExtension->GetAOD()->Print();
  }
  if(fAODOut){
    Printf("%s:%d >>>>>> Output",(char*)__FILE__,__LINE__);
    fAODOut->Print();
  }
}
Commit	Line	Data
fb10c4b8	1	// **************************************
	2	// Task used for the correction of determiantion of reconstructed jet spectra
	3	// Compares input (gen) and output (rec) jets
	4	// *******************************************
	5
	6
	7	/**************************************************************************
	8	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	9	* *
	10	* Author: The ALICE Off-line Project. *
	11	* Contributors are mentioned in the code where appropriate. *
	12	* *
	13	* Permission to use, copy, modify and distribute this software and its *
	14	* documentation strictly for non-commercial purposes is hereby granted *
	15	* without fee, provided that the above copyright notice appears in all *
	16	* copies and that both the copyright notice and this permission notice *
	17	* appear in the supporting documentation. The authors make no claims *
	18	* about the suitability of this software for any purpose. It is *
	19	* provided "as is" without express or implied warranty. *
	20	**************************************************************************/
	21
	22
	23	#include <TROOT.h>
	24	#include <TH1F.h>
d9c0c425	25	#include <TH2F.h>
fb10c4b8	26	#include <THnSparse.h>
	27	#include <TSystem.h>
	28	#include <TInterpreter.h>
	29	#include <TList.h>
	30	#include <TLorentzVector.h>
d9c0c425	31	#include <TRefArray.h>
fb10c4b8	32	#include "TDatabasePDG.h"
	33
	34	#include "AliAnalysisTaskJetBackgroundSubtract.h"
	35	#include "AliAnalysisManager.h"
	36	#include "AliAODHandler.h"
	37	#include "AliAODTrack.h"
	38	#include "AliAODJet.h"
	39	#include "AliAODEvent.h"
	40	#include "AliInputEventHandler.h"
	41	#include "AliAODJetEventBackground.h"
	42
	43
	44	ClassImp(AliAnalysisTaskJetBackgroundSubtract)
	45
	46	AliAnalysisTaskJetBackgroundSubtract::~AliAnalysisTaskJetBackgroundSubtract(){
	47	delete fJBArray;
	48	delete fOutJetArrayList;
	49	delete fInJetArrayList;
	50	}
	51
	52	AliAnalysisTaskJetBackgroundSubtract::AliAnalysisTaskJetBackgroundSubtract():
	53	AliAnalysisTaskSE(),
	54	fAODOut(0x0),
	55	fAODIn(0x0),
	56	fAODExtension(0x0),
	57	fJBArray(new TObjArray()),
	58	fBackgroundBranch(""),
	59	fNonStdFile(""),
	60	fReplaceString1("B0"),
	61	fReplaceString2("B%d"),
a18eedbb	62	fSubtraction(kRhoRecalc),
fb10c4b8	63	fInJetArrayList(0x0),
	64	fOutJetArrayList(0x0),
	65	fHistList(0x0)
	66	{
	67
	68	}
	69
	70	AliAnalysisTaskJetBackgroundSubtract::AliAnalysisTaskJetBackgroundSubtract(const char* name):
	71
	72	AliAnalysisTaskSE(name),
	73	fAODOut(0x0),
	74	fAODIn(0x0),
	75	fAODExtension(0x0),
	76	fJBArray(new TObjArray()),
	77	fBackgroundBranch(""),
	78	fNonStdFile(""),
	79	fReplaceString1("B0"),
	80	fReplaceString2("B%d"),
a18eedbb	81	fSubtraction(kRhoRecalc),
fb10c4b8	82	fInJetArrayList(0x0),
	83	fOutJetArrayList(0x0),
	84	fHistList(0x0)
	85	{
	86	DefineOutput(1, TList::Class());
	87	}
	88
	89
	90
	91	Bool_t AliAnalysisTaskJetBackgroundSubtract::Notify()
	92	{
	93	//
1c21a0e7	94	fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
	95
	96	ResetOutJets();
	97
	98	// Now we also have the Input Event Available! Fillvthe pointers in the list
	99
	100	fInJetArrayList->Clear();
	101	fOutJetArrayList->Clear();
	102
	103	for(int iJB = 0;iJB<fJBArray->GetEntries();iJB++){
	104	TObjString ostr = (TObjString)fJBArray->At(iJB);
a18eedbb	105
a18eedbb	106
1c21a0e7	107	TClonesArray* jarray = 0;
	108	if(!jarray&&fAODOut){
	109	jarray = (TClonesArray*)(fAODOut->FindListObject(ostr->GetString().Data()));
	110	}
	111	if(!jarray&&fAODExtension){
	112	jarray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(ostr->GetString().Data()));
	113	}
	114	if(!jarray&&fAODIn){
	115	jarray = (TClonesArray*)(fAODIn->FindListObject(ostr->GetString().Data()));
	116	}
	117
	118	if(!jarray){
	119	if(fDebug){
	120	Printf("%s:%d Input jet branch %s not found",(char*)__FILE__,__LINE__,ostr->GetString().Data());
	121	}
	122	continue;
	123	}
	124
	125	TString newName(ostr->GetString().Data());
	126	newName.ReplaceAll(fReplaceString1.Data(),Form(fReplaceString2.Data(),fSubtraction));
	127	TClonesArray* jarrayOut = 0;
	128	if(!jarrayOut&&fAODOut){
	129	jarrayOut = (TClonesArray*)(fAODOut->FindListObject(newName.Data()));
	130	}
	131	if(!jarrayOut&&fAODExtension){
	132	jarrayOut = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(newName.Data()));
	133	}
	134
	135	if(!jarrayOut){
	136	if(fDebug){
	137	Printf("%s:%d Output jet branch %s not found",(char*)__FILE__,__LINE__,newName.Data());
	138	PrintAODContents();
	139	}
	140	continue;
	141	}
	142	if(jarrayOut&&jarray){
	143	fOutJetArrayList->Add(jarrayOut);
	144	fInJetArrayList->Add(jarray);
	145	}
	146	}
fb10c4b8	147	return kTRUE;
	148	}
	149
	150	void AliAnalysisTaskJetBackgroundSubtract::UserCreateOutputObjects()
	151	{
	152
	153	//
	154	// Create the output container
	155	//
	156	// Connect the AOD
	157
fb10c4b8	158	if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::UserCreateOutputObjects() \n");
fb10c4b8	159	if(fNonStdFile.Length()!=0){
a18eedbb	160
fb10c4b8	161	// case that we have an AOD extension we need to fetch the jets from the extended output
	162	// we identifay the extension aod event by looking for the branchname
	163	AliAODHandler aodH = dynamic_cast<AliAODHandler>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
	164	fAODExtension = aodH->GetExtension(fNonStdFile.Data());
	165
	166	if(!fAODExtension){
	167	if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
	168	}
	169	}
	170	fAODOut = AODEvent();
fb10c4b8	171
1c21a0e7	172	// usually we do not have the input already here
fb10c4b8	173
	174	if(!fInJetArrayList)fInJetArrayList =new TList();
	175	if(!fOutJetArrayList)fOutJetArrayList =new TList();
	176
	177	for(int iJB = 0;iJB<fJBArray->GetEntries();iJB++){
fb10c4b8	178	TObjString ostr = (TObjString)fJBArray->At(iJB);
1c21a0e7	179	TString newName(ostr->GetString().Data());
fb10c4b8	180	if(!newName.Contains(fReplaceString1.Data())){
fb10c4b8	181	Printf("%s:%d cannot replace string %s in %s",(char*)__FILE__,__LINE__,fReplaceString1.Data(),
1c21a0e7	182	newName.Data());
fb10c4b8	183	continue;
	184	}
	185
fb10c4b8	186	// add a new branch to the output for the background subtracted jets take the names from
	187	// the input jets and replace the background flag names
	188	TClonesArray *tca = new TClonesArray("AliAODJet", 0);
fb10c4b8	189	newName.ReplaceAll(fReplaceString1.Data(),Form(fReplaceString2.Data(),fSubtraction));
	190	if(fDebug){
	191	Printf("%s:%d created branch \n %s from \n %s",(char*)__FILE__,__LINE__,newName.Data(),
1c21a0e7	192	ostr->GetString().Data());
fb10c4b8	193	}
	194	tca->SetName(newName.Data());
	195	AddAODBranch("TClonesArray",&tca,fNonStdFile.Data());
fb10c4b8	196	}
	197
	198
fb10c4b8	199	if(!fHistList)fHistList = new TList();
	200	fHistList->SetOwner();
	201
1c21a0e7	202	for(int iJB = 0;iJB<fJBArray->GetEntries();iJB++){
	203	TObjString ostr = (TObjString)fJBArray->At(iJB);
	204	TString oldName(ostr->GetString().Data());
	205	TString newName(ostr->GetString().Data());
	206	if(!newName.Contains(fReplaceString1.Data())){
	207	Printf("%s:%d cannot replace string %s in %s",(char*)__FILE__,__LINE__,fReplaceString1.Data(),
	208	newName.Data());
	209	continue;
	210	}
	211	newName.ReplaceAll(fReplaceString1.Data(),Form(fReplaceString2.Data(),fSubtraction));
	212	TH2F *hTmp = new TH2F(Form("h2PtInPtOut_%d",iJB),Form(";%s p_{T}; %s p_{T}",oldName.Data(),newName.Data()),200,0,200.,200,0.,200.);
d9c0c425	213	fHistList->Add(hTmp);
	214	}
	215
fb10c4b8	216	Bool_t oldStatus = TH1::AddDirectoryStatus();
	217	TH1::AddDirectory(kFALSE);
	218
	219	//
	220	// Histogram booking, add som control histograms here
	221	//
	222
	223	// =========== Switch on Sumw2 for all histos ===========
	224	for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
	225	TH1 h1 = dynamic_cast<TH1>(fHistList->At(i));
	226	if (h1){
	227	h1->Sumw2();
	228	continue;
	229	}
	230	THnSparse hn = dynamic_cast<THnSparse>(fHistList->At(i));
	231	if(hn)hn->Sumw2();
	232	}
	233	TH1::AddDirectory(oldStatus);
	234
	235	if(fBackgroundBranch.Length()==0)
	236	AliError(Form("%s:%d No BackgroundBranch defined",(char*)__FILE__,__LINE__));
	237	if(fJBArray->GetEntries()==0)
	238	AliError(Form("%s:%d No Jet Branches defined defined",(char*)__FILE__,__LINE__));
	239	}
	240
	241	void AliAnalysisTaskJetBackgroundSubtract::Init()
	242	{
	243	//
	244	// Initialization
	245	//
	246	if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::Init() \n");
	247	}
	248
	249	void AliAnalysisTaskJetBackgroundSubtract::UserExec(Option_t /option*/)
	250	{
	251
1c21a0e7	252	if (fDebug > 1) printf("AnalysisTaskJetBackgroundSubtract::UserExec() \n");
fb10c4b8	253	ResetOutJets();
fb10c4b8	254	if(fBackgroundBranch.Length()==0\|\|fJBArray->GetEntries()==0){
1c21a0e7	255	if(fDebug)Printf("%s:%d No background subtraction done",(char*)__FILE__,__LINE__);
fb10c4b8	256	PostData(1,fHistList);
fb10c4b8	257	}
1c21a0e7	258	if(fJBArray->GetEntries()!=fInJetArrayList->GetEntries()){
	259	if(fDebug)Printf("%s:%d different Array sizes %d %d %d",(char*)__FILE__,__LINE__,fJBArray->GetEntries(),fInJetArrayList->GetEntries(),fOutJetArrayList->GetEntries());
	260	PostData(1,fHistList);
	261	}
	262
fb10c4b8	263
fb10c4b8	264
a18eedbb	265	AliAODJetEventBackground* evBkg = 0;
	266	TClonesArray* bkgClusters = 0;
	267	TClonesArray* bkgClustersRC = 0;
	268	TString bkgClusterName(fBackgroundBranch.Data());
1c21a0e7	269	bkgClusterName.ReplaceAll(Form("%s_",AliAODJetEventBackground::StdBranchName()),"");
a18eedbb	270	TString bkgClusterRCName(Form("%s%s",bkgClusterName.Data(),"RandomCone"));
	271
	272	if(!evBkg&&!bkgClusters&&!bkgClustersRC&&fAODOut){
fb10c4b8	273	evBkg = (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
d9c0c425	274	bkgClusters = (TClonesArray*)(fAODOut->FindListObject(bkgClusterName.Data()));
a18eedbb	275	bkgClustersRC = (TClonesArray*)(fAODOut->FindListObject(bkgClusterRCName.Data()));
a18eedbb	276
1c21a0e7	277	if(fDebug&&bkgClusters)Printf("%s:%d Background cluster branch %s found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
a18eedbb	278	if(fDebug&&bkgClustersRC)Printf("%s:%d Background cluster RC branch %s found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
1c21a0e7	279	if(fDebug&&evBkg)Printf("%s:%d Backgroundbranch %s found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
fb10c4b8	280	}
a18eedbb	281	if(!evBkg&&!bkgClusters&&!bkgClustersRC&&fAODExtension){
fb10c4b8	282	evBkg = (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
d9c0c425	283	bkgClusters = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(bkgClusterName.Data()));
a18eedbb	284	bkgClustersRC = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(bkgClusterRCName.Data()));
1c21a0e7	285	if(fDebug&&bkgClusters)Printf("%s:%d Background cluster branch %s found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
a18eedbb	286	if(fDebug&&bkgClustersRC)Printf("%s:%d Background cluster RC branch %s found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
a18eedbb	287
1c21a0e7	288	if(fDebug&&evBkg)Printf("%s:%d Backgroundbranch %s found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
fb10c4b8	289	}
1c21a0e7	290
a18eedbb	291	if(!evBkg&&!bkgClusters&&!bkgClustersRC&&fAODIn){
fb10c4b8	292	evBkg = (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data()));
1c21a0e7	293	bkgClusters = (TClonesArray*)(fAODIn->FindListObject(bkgClusterName.Data()));
a18eedbb	294	bkgClustersRC = (TClonesArray*)(fAODIn->FindListObject(bkgClusterRCName.Data()));
a18eedbb	295
1c21a0e7	296	if(fDebug&&bkgClusters)Printf("%s:%d Background cluster branch %s found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
a18eedbb	297	if(fDebug&&bkgClustersRC)Printf("%s:%d Background cluster RC branch %s found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
1c21a0e7	298	if(fDebug&&evBkg)Printf("%s:%d Backgroundbranch %s found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
fb10c4b8	299	}
fb10c4b8	300
8a56597f	301	if(!evBkg&&(fSubtraction==kArea\|\|fSubtraction==kRhoRecalc)){
1c21a0e7	302	if(fDebug){
	303	Printf("%s:%d Backgroundbranch %s not found",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());
	304	PrintAODContents();
	305	}
fb10c4b8	306	PostData(1,fHistList);
d9c0c425	307	return;
	308	}
	309
8a56597f	310	if(!bkgClusters&&(fSubtraction==kRhoRecalc)){
1c21a0e7	311	if(fDebug){
	312	Printf("%s:%d Background cluster branch %s not found",(char*)__FILE__,__LINE__,bkgClusterName.Data());
	313	PrintAODContents();
	314	}
d9c0c425	315	PostData(1,fHistList);
d9c0c425	316	return;
fb10c4b8	317	}
fb10c4b8	318
8a56597f	319	if(!bkgClustersRC&&(fSubtraction==kRhoRC)){
a18eedbb	320	if(fDebug){
	321	Printf("%s:%d Background cluster RC branch %s not found",(char*)__FILE__,__LINE__,bkgClusterRCName.Data());
	322	PrintAODContents();
	323	}
	324	PostData(1,fHistList);
	325	return;
	326	}
fb10c4b8	327	// LOOP over all jet branches and subtract the background
fb10c4b8	328
a18eedbb	329	Float_t rho = 0;
	330	Double_t meanarea = 0;
	331	if(fSubtraction==kArea)rho = evBkg->GetBackground(2);
	332	if(fSubtraction==kRhoRecalc){
	333	meanarea=evBkg->GetMeanarea(2);
	334	rho =RecalcRho(bkgClusters,meanarea);
	335	}
	336	if(fSubtraction==kRhoRC)rho =RhoRC(bkgClustersRC);
	337
fb10c4b8	338	for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
	339	TClonesArray* jarray = (TClonesArray*)fInJetArrayList->At(iJB);
	340	TClonesArray* jarrayOut = (TClonesArray*)fOutJetArrayList->At(iJB);
1c21a0e7	341
1c21a0e7	342	if(!jarray\|\|!jarrayOut){
a18eedbb	343	Printf("%s:%d Array not found %d: %p %p",(char*)__FILE__,__LINE__,iJB,jarray->GetName(),jarrayOut->GetName());
	344	continue;
	345	}
d9c0c425	346	TH2F* h2PtInOut = (TH2F*)fHistList->FindObject(Form("h2PtInPtOut_%d",iJB));
fb10c4b8	347	// loop over all jets
fb10c4b8	348	Int_t nOut = 0;
a18eedbb	349
a18eedbb	350
fb10c4b8	351	for(int i = 0;i < jarray->GetEntriesFast();i++){
	352	AliAODJet jet = (AliAODJet)jarray->At(i);
	353	AliAODJet tmpNewJet(*jet);
	354	Bool_t bAdd = false;
a18eedbb	355
a18eedbb	356
fb10c4b8	357	if(fSubtraction==kArea){
d9c0c425	358	Double_t background = rho * jet->EffectiveAreaCharged();
d9c0c425	359	Float_t ptSub = jet->Pt() - background;
fb10c4b8	360	if(fDebug>2){
	361	Printf("%s:%d Jet %d %3.3f %3.3f",(char*)__FILE__,__LINE__,i,jet->Pt(),ptSub);
	362	}
	363	if(ptSub<0){
	364	// optionally rescale it and keep??
	365	bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
d9c0c425	366	if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),0.1);
fb10c4b8	367	}
	368	else{
	369	bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
d9c0c425	370	if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
fb10c4b8	371	}
d9c0c425	372	// add background estimates to the new jet object
	373	// allows to recover old p_T and rho...
	374	tmpNewJet.SetBgEnergy(background,0);
fb10c4b8	375	}// kAREA
a18eedbb	376	else if(fSubtraction==kRhoRecalc){
	377	Double_t background = rho * jet->EffectiveAreaCharged();
	378	Float_t ptSub = jet->Pt() - background;
	379	if(fDebug>2){
	380	Printf("%s:%d Jet %d %3.3f %3.3f %3.3f %3.3f",(char*)__FILE__,__LINE__,i,jet->Pt(),ptSub,background,rho);}
	381	if(ptSub<0){
	382	// optionally rescale it and keep??
	383	bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
	384	if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),0.1);
	385	}
	386	else{
	387	bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
	388	if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
	389	}
	390	// add background estimates to the new jet object
	391	// allows to recover old p_T and rho...
	392	tmpNewJet.SetBgEnergy(background,0);
	393
	394	}//kRhoRecalc
	395	else if(fSubtraction==kRhoRC){
	396	Double_t background = rho * jet->EffectiveAreaCharged();
	397	Float_t ptSub = jet->Pt() - background;
	398	if(fDebug>2){ Printf("%s:%d Jet %d %3.3f %3.3f %3.3f %3.3f",(char*)__FILE__,__LINE__,i,jet->Pt(),ptSub,background,rho);}
	399	if(ptSub<0){
	400	// optionally rescale it and keep??
	401	bAdd = RescaleJetMomentum(&tmpNewJet,0.1);
	402	if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),0.1);
	403	}
	404	else{
	405	bAdd = RescaleJetMomentum(&tmpNewJet,ptSub);
	406	if(h2PtInOut)h2PtInOut->Fill(jet->Pt(),ptSub);
	407	}
	408	// add background estimates to the new jet object
	409	// allows to recover old p_T and rho...
	410	tmpNewJet.SetBgEnergy(background,0);
	411
	412	}//kRhoRC
fb10c4b8	413
fb10c4b8	414	if(bAdd){
d9c0c425	415	AliAODJet newJet = new ((jarrayOut)[nOut++]) AliAODJet(tmpNewJet);
	416	// what about track references, clear for now...
	417	newJet->GetRefTracks()->Clear();
fb10c4b8	418	}
fb10c4b8	419
a18eedbb	420
	421
	422
	423
	424
	425	}
fb10c4b8	426
	427
	428
	429	// subtract the background
	430
	431
	432	// remove jets??
	433
	434	// sort jets...
	435
	436	}
	437	PostData(1, fHistList);
	438	}
	439
	440	void AliAnalysisTaskJetBackgroundSubtract::Terminate(Option_t /option*/)
	441	{
	442	// Terminate analysis
	443	//
	444	if (fDebug > 1) printf("AnalysisJetBackgroundSubtract: Terminate() \n");
	445	}
	446
	447	Bool_t AliAnalysisTaskJetBackgroundSubtract::RescaleJetMomentum(AliAODJet *jet,Float_t pT){
	448	// keep the direction and the jet mass
	449	if(pT<=0)return kFALSE;
	450	Double_t pTold = jet->Pt();
	451	Double_t scale = pT/pTold;
	452	Double_t mass = jet->M();
	453	Double_t pNew = jet->P() * scale;
	454	jet->SetPxPyPzE(scalejet->Px(),scalejet->Py(),scalejet->Pz(),TMath::Sqrt(massmass+pNew*pNew));
	455	return kTRUE;
	456	}
	457
a18eedbb	458	Double_t AliAnalysisTaskJetBackgroundSubtract::RecalcRho(TClonesArray* bkgClusters,Double_t meanarea){
	459
	460	Double_t ptarea=0.;
	461	Int_t count=0;
	462	Double_t rho=0.;
	463	const Double_t Rlimit2=0.80.8; //2jet radius.
	464	TClonesArray* jarray=0;
	465
	466	for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
	467	TObjString ostr = (TObjString)fInJetArrayList->At(iJB);
	468	TString jetref=ostr->GetString().Data();
	469	if(jetref.Contains("ANTIKT04")){
	470	jarray = (TClonesArray*)fInJetArrayList->At(iJB);}}
	471	if(jarray->GetEntries()>=2){
	472	AliAODJet first = (AliAODJet)(jarray->At(0));
	473	AliAODJet second= (AliAODJet)(jarray->At(1));
	474	for(Int_t k=0;k<bkgClusters->GetEntriesFast();k++){
	475	AliAODJet clus = (AliAODJet)(bkgClusters->At(k));
	476	if(TMath::Abs(clus->Eta())>0.5) continue;
	477	if((clus->EffectiveAreaCharged())<0.1*meanarea) continue;
	478	Double_t distance1=(first->Eta()-clus->Eta())*(first->Eta()-clus->Eta())+
	479	(first->Phi()-clus->Phi())*(first->Phi()-clus->Phi());
	480	Double_t distance2= (second->Eta()-clus->Eta())*(second->Eta()-clus->Eta())+
	481	(second->Phi()-clus->Phi())*(second->Phi()-clus->Phi());
	482	if((distance1<Rlimit2)\|\|(distance2<Rlimit2)) continue;
	483	ptarea=ptarea+clus->Pt()/clus->EffectiveAreaCharged();
	484	count=count+1;}
	485	if(count!=0) rho=ptarea/count;
	486	}
	487	return rho;
	488	}
	489
	490	Double_t AliAnalysisTaskJetBackgroundSubtract::RhoRC(TClonesArray* bkgClustersRC){
	491
	492	Double_t ptarea=0.;
	493	Int_t count=0;
	494	Double_t rho=0.;
	495	const Double_t Rlimit2=0.80.8; //2jet radius.
	496	TClonesArray* jarray=0;
	497	for(int iJB = 0;iJB<fInJetArrayList->GetEntries();iJB++){
	498	TObjString ostr = (TObjString)fInJetArrayList->At(iJB);
	499	TString jetref=ostr->GetString().Data();
	500	if(jetref.Contains("ANTIKT04")){
	501	jarray = (TClonesArray*)fInJetArrayList->At(iJB);}}
	502	if(jarray->GetEntries()>=2){
	503	AliAODJet first = (AliAODJet)(jarray->At(0));
	504	AliAODJet second=(AliAODJet)(jarray->At(1));
	505	for(Int_t k=0;k<bkgClustersRC->GetEntriesFast();k++){
	506	AliAODJet clus = (AliAODJet)(bkgClustersRC->At(k));
	507	if(TMath::Abs(clus->Eta())>0.5) continue;
	508	Double_t distance1=(first->Eta()-clus->Eta())*(first->Eta()-clus->Eta())+
	509	(first->Phi()-clus->Phi())*(first->Phi()-clus->Phi());
	510	Double_t distance2= (second->Eta()-clus->Eta())*(second->Eta()-clus->Eta())+
	511	(second->Phi()-clus->Phi())*(second->Phi()-clus->Phi());
	512	if((distance1<Rlimit2)\|\|(distance2<Rlimit2)) continue;
	513	ptarea=ptarea+clus->Pt()/clus->EffectiveAreaCharged();
	514	count=count+1;}
	515	if(count!=0) rho=ptarea/count; }
	516	return rho;
	517	}
	518
	519
	520
	521
522
523
524
525
526
fb10c4b8	527	void AliAnalysisTaskJetBackgroundSubtract::ResetOutJets(){
	528	if(!fOutJetArrayList)return;
	529	for(int iJB = 0;iJB<fOutJetArrayList->GetEntries();iJB++){
	530	TClonesArray* jarray = (TClonesArray*)fOutJetArrayList->At(iJB);
	531	if(jarray)jarray->Delete();
	532	}
	533	}
1c21a0e7	534
	535
	536	void AliAnalysisTaskJetBackgroundSubtract::PrintAODContents(){
	537	if(fAODIn){
	538	Printf("%s:%d >>>>>> Input",(char*)__FILE__,__LINE__);
	539	fAODIn->Print();
	540	}
	541	if(fAODExtension){
	542	Printf("%s:%d >>>>>> Extenstion",(char*)__FILE__,__LINE__);
	543	fAODExtension->GetAOD()->Print();
	544	}
	545	if(fAODOut){
	546	Printf("%s:%d >>>>>> Output",(char*)__FILE__,__LINE__);
	547	fAODOut->Print();
	548	}
	549	}