[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetResponse.cxx

#include "TChain.h"
#include "TTree.h"
#include "TMath.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TCanvas.h"

#include "AliLog.h"

#include "AliAnalysisTask.h"
#include "AliAnalysisManager.h"

#include "AliVEvent.h"
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliCentrality.h"
#include "AliAnalysisHelperJetTasks.h"
#include "AliInputEventHandler.h"

#include "AliAODEvent.h"
#include "AliAODJet.h"

#include "AliAnalysisTaskJetResponse.h"

ClassImp(AliAnalysisTaskJetResponse)

AliAnalysisTaskJetResponse::AliAnalysisTaskJetResponse() :
  AliAnalysisTaskSE(),
  fESD(0x0),
  fAOD(0x0),
  fOfflineTrgMask(AliVEvent::kAny),
  fMinContribVtx(1),
  fVtxZMin(-8.),
  fVtxZMax(8.),
  fEvtClassMin(1),
  fEvtClassMax(4),
  fCentMin(0.),
  fCentMax(100.),
  fJetEtaMin(-.5),
  fJetEtaMax(.5),
  fJetPtMin(20.),
  fJetPtFractionMin(0.5),
  fNMatchJets(4),
  //fJetDeltaEta(.4),
  //fJetDeltaPhi(.4),
  fkNbranches(2),
  fkEvtClasses(5),
  fOutputList(0x0),
  fHistEvtSelection(0x0),
  fHistEvtClass(0x0),
  fHistCentrality(0x0),
  fHistPtJet(new TH1F*[fkNbranches]),
  fHistEtaPhiJet(new TH2F*[fkNbranches]),
  fHistEtaPhiJetCut(new TH2F*[fkNbranches]),
  fHistDeltaEtaDeltaPhiJet(new TH2F*[fkEvtClasses]),
  fHistDeltaEtaDeltaPhiJetCut(new TH2F*[fkEvtClasses]),
  fHistDeltaEtaDeltaPhiJetNOMatching(new TH2F*[fkEvtClasses]),
  fHistDeltaEtaEtaJet(new TH2F*[fkEvtClasses]),
  fHistDeltaPtEtaJet(new TH2F*[fkEvtClasses]),
  fHistPtFraction(new TH2F*[fkEvtClasses]),
  fHistPtPtExtra(0x0),
  fHistPtResponse(new TH2F*[fkEvtClasses]),
  fHistPtSmearing(new TH2F*[fkEvtClasses]),
  fHistDeltaR(new TH2F*[fkEvtClasses]),
  fHistArea(new TH2F*[fkEvtClasses]),
  fHistDeltaArea(new TH2F*[fkEvtClasses])
{
  // default Constructor

  fJetBranchName[0] = "";
  fJetBranchName[1] = "";

  fListJets[0] = new TList;
  fListJets[1] = new TList;
}

AliAnalysisTaskJetResponse::AliAnalysisTaskJetResponse(const char *name) :
  AliAnalysisTaskSE(name),
  fESD(0x0),
  fAOD(0x0),
  fOfflineTrgMask(AliVEvent::kAny),
  fMinContribVtx(1),
  fVtxZMin(-8.),
  fVtxZMax(8.),
  fEvtClassMin(1),
  fEvtClassMax(4),
  fCentMin(0.),
  fCentMax(100.),
  fJetEtaMin(-.5),
  fJetEtaMax(.5),
  fJetPtMin(20.),
  fJetPtFractionMin(0.5),
  fNMatchJets(4),
  //fJetDeltaEta(.4),
  //fJetDeltaPhi(.4),
  fkNbranches(2),
  fkEvtClasses(5),
  fOutputList(0x0),
  fHistEvtSelection(0x0),
  fHistEvtClass(0x0),
  fHistCentrality(0x0),
  fHistPtJet(new TH1F*[fkNbranches]),
  fHistEtaPhiJet(new TH2F*[fkNbranches]),
  fHistEtaPhiJetCut(new TH2F*[fkNbranches]),
  fHistDeltaEtaDeltaPhiJet(new TH2F*[fkEvtClasses]),
  fHistDeltaEtaDeltaPhiJetCut(new TH2F*[fkEvtClasses]),
  fHistDeltaEtaDeltaPhiJetNOMatching(new TH2F*[fkEvtClasses]),
  fHistDeltaEtaEtaJet(new TH2F*[fkEvtClasses]),
  fHistDeltaPtEtaJet(new TH2F*[fkEvtClasses]),
  fHistPtFraction(new TH2F*[fkEvtClasses]),
  fHistPtPtExtra(0x0),
  fHistPtResponse(new TH2F*[fkEvtClasses]),
  fHistPtSmearing(new TH2F*[fkEvtClasses]),
  fHistDeltaR(new TH2F*[fkEvtClasses]),
  fHistArea(new TH2F*[fkEvtClasses]),
  fHistDeltaArea(new TH2F*[fkEvtClasses])
{
  // Constructor

  fJetBranchName[0] = "";
  fJetBranchName[1] = "";

  fListJets[0] = new TList;
  fListJets[1] = new TList;

  DefineOutput(1, TList::Class());
}

AliAnalysisTaskJetResponse::~AliAnalysisTaskJetResponse()
{
  delete fListJets[0];
  delete fListJets[1];
}

void AliAnalysisTaskJetResponse::SetBranchNames(const TString &branch1, const TString &branch2)
{
  fJetBranchName[0] = branch1;
  fJetBranchName[1] = branch2;
}

void AliAnalysisTaskJetResponse::Init()
{

   // check for jet branches
   if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
      AliError("Jet branch name not set.");
   }

}

void AliAnalysisTaskJetResponse::UserCreateOutputObjects()
{
  // Create histograms
  // Called once
  OpenFile(1);
  if(!fOutputList) fOutputList = new TList;
  fOutputList->SetOwner(kTRUE);

  Bool_t oldStatus = TH1::AddDirectoryStatus();
  TH1::AddDirectory(kFALSE);


  fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 5, -0.5, 4.5);
  fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
  fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
  fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
  fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
  fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
  
  fHistEvtClass   = new TH1I("fHistEvtClass", "event classes", fkEvtClasses, -0.5, fkEvtClasses-0.5);
  fHistCentrality = new TH1F("fHistCentrality", "event centrality", 100, 0., 100.);

  for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
    fHistPtJet[iBranch] = new TH1F(Form("fHistPtJet%i", iBranch),
					  Form("p_{T} of jets from branch %i;p_{T} (GeV/c);counts", iBranch),
					  250, 0., 250.);
    fHistPtJet[iBranch]->SetMarkerStyle(kFullCircle);

    fHistEtaPhiJet[iBranch]    = new TH2F(Form("fHistEtaPhiJet%i", iBranch),
					                 Form("#eta - #phi of jets from branch %i (before cut);#eta;#phi", iBranch),
					                 48, -1.2, 1.2, 100, 0., 2*TMath::Pi());
	fHistEtaPhiJetCut[iBranch] = new TH2F(Form("fHistEtaPhiJetCut%i", iBranch),
					                 Form("#eta - #phi of jets from branch %i (before cut);#eta;#phi", iBranch),
					                 48, -1.2, 1.2, 100, 0., 2*TMath::Pi());

  }
  

  fHistPtPtExtra = new TH2F("fHistPtPtExtra", "p_{T} response;p_{T} (GeV/c);p_{T} (GeV/c)",
			    250, 0., 250., 250, 0., 250.);

  for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
    fHistDeltaEtaDeltaPhiJet[iEvtClass] = new TH2F(Form("fHistDeltaEtaDeltaPhiJet%i",iEvtClass),
					                              "#Delta#eta - #Delta#phi of jet;#Delta#eta;#Delta#phi",
					                              101, -1.01, 1.01, 101, -1.01, 1.01);
    fHistDeltaEtaDeltaPhiJetCut[iEvtClass] = new TH2F(Form("fHistDeltaEtaDeltaPhiJetCut%i",iEvtClass),
					                              "#Delta#eta - #Delta#phi of jet;#Delta#eta;#Delta#phi",
					                              101, -1.01, 1.01, 101, -1.01, 1.01);												  
    fHistDeltaEtaDeltaPhiJetNOMatching[iEvtClass] = new TH2F("fHistDeltaEtaDeltaPhiJetNOMatching",
					                                         "#Delta#eta - #Delta#phi of jets which do not match;#Delta#eta;#Delta#phi",
					                                         100, -2., 2., 100, TMath::Pi(), TMath::Pi());
    fHistPtResponse[iEvtClass] = new TH2F(Form("pt_response%i",iEvtClass), Form("pt_response%i;p_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
					                      250, 0., 250., 250, 0., 250.);
    fHistPtSmearing[iEvtClass] = new TH2F(Form("pt_smearing%i",iEvtClass), Form("pt_smearing%i;#Deltap_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
					                      200, -50., 150., 250, 0., 250.);
    fHistDeltaR[iEvtClass]     = new TH2F(Form("hist_DeltaR%i",iEvtClass), "#DeltaR of matched jets;#Deltap_{T};#DeltaR", 200, -50.,150., 60, 0.,.6);
    fHistArea[iEvtClass]       = new TH2F(Form("hist_Area%i",iEvtClass), "jet area;#Deltap_{T};jet area", 200, -50.,150., 100, 0.,1.0);
    fHistDeltaArea[iEvtClass]  = new TH2F(Form("hist_DeltaArea%i",iEvtClass), "#DeltaArea of matched jets;#Deltap_{T};#Deltaarea", 200, -50., 150., 60, 0., .3); 
	
    fHistDeltaEtaEtaJet[iEvtClass] = new TH2F(Form("fHistDeltaEtaEtaJet%i", iEvtClass),
                                                Form("#eta - #Delta#eta of matched jets from event class %i;#eta;#Delta#eta", iEvtClass),
                                                60, -.6, .6, 100, -.5, .5);
    fHistDeltaPtEtaJet[iEvtClass] = new TH2F(Form("fHistDeltaPtEtaJet%i", iEvtClass),
                                                Form("#eta - #Deltap_{T} of matched jets from event class %i;#eta;#Deltap_{T}", iEvtClass),
                                                60, -.6, .6, 200, -50., 150);
    fHistPtFraction[iEvtClass] = new TH2F(Form("fHistPtFraction%i", iEvtClass),
	                                      Form("fraction from embedded jet in reconstructed jet;fraction (event class %i);p_{T}^{emb}", iEvtClass),
										  100, 0., 1., 250, 0., 250.);
  }

  
  fOutputList->Add(fHistEvtSelection);
  fOutputList->Add(fHistEvtClass);
  fOutputList->Add(fHistCentrality);

  for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
    fOutputList->Add(fHistPtJet[iBranch]);
    fOutputList->Add(fHistEtaPhiJet[iBranch]);
	fOutputList->Add(fHistEtaPhiJetCut[iBranch]);
  }
  
  fOutputList->Add(fHistPtPtExtra);
  
  for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
    fOutputList->Add(fHistDeltaEtaDeltaPhiJet[iEvtClass]);
    fOutputList->Add(fHistDeltaEtaDeltaPhiJetCut[iEvtClass]);
    fOutputList->Add(fHistDeltaEtaDeltaPhiJetNOMatching[iEvtClass]);
	fOutputList->Add(fHistDeltaEtaEtaJet[iEvtClass]);
    fOutputList->Add(fHistDeltaPtEtaJet[iEvtClass]);
	fOutputList->Add(fHistPtFraction[iEvtClass]);
	
    fOutputList->Add(fHistPtResponse[iEvtClass]);
    fOutputList->Add(fHistPtSmearing[iEvtClass]);
    fOutputList->Add(fHistDeltaR[iEvtClass]);
    fOutputList->Add(fHistArea[iEvtClass]);
    fOutputList->Add(fHistDeltaArea[iEvtClass]);
  }

  // =========== Switch on Sumw2 for all histos ===========
  for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
    TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
    if (h1){
      h1->Sumw2();
      continue;
    }
  }
  TH1::AddDirectory(oldStatus);

  PostData(1, fOutputList);
}

void AliAnalysisTaskJetResponse::UserExec(Option_t *)
{
  // load events, apply event cuts, then compare leading jets

  if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
      AliError("Jet branch name not set.");
      return;
   }

  fESD=dynamic_cast<AliESDEvent*>(InputEvent());
  if (!fESD) {
    AliError("ESD not available");
    return;
  }
  fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
  if (!fAOD) {
    AliError("AOD not available");
    return;
  }

  // -- event selection --
  fHistEvtSelection->Fill(1); // number of events before event selection

  // physics selection
  AliInputEventHandler* inputHandler = (AliInputEventHandler*)
    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
  if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
    if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
    fHistEvtSelection->Fill(2);
    PostData(1, fOutputList);
    return;
  }

  // vertex selection
  AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
  Int_t nTracksPrim = primVtx->GetNContributors();
  if ((nTracksPrim < fMinContribVtx) ||
      (primVtx->GetZ() < fVtxZMin) ||
      (primVtx->GetZ() > fVtxZMax) ){
    if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
    fHistEvtSelection->Fill(3);
    PostData(1, fOutputList);
    return;
  }

  // event class selection (from jet helper task)
  Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
  if(fDebug) Printf("Event class %d", eventClass);
  if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
    fHistEvtSelection->Fill(4);
    PostData(1, fOutputList);
    return;
  }

  // centrality selection
  AliCentrality *cent = fESD->GetCentrality();
  Float_t centValue = cent->GetCentralityPercentile("TRK");
  if(fDebug) printf("centrality: %f\n", centValue);
  if (centValue < fCentMin || centValue > fCentMax){
    fHistEvtSelection->Fill(4);
    PostData(1, fOutputList);
    return;
  }

  fHistEvtSelection->Fill(0); // accepted events
  fHistEvtClass->Fill(eventClass);
  fHistCentrality->Fill(centValue);
  // -- end event selection --

  // fetch jets
  TClonesArray *aodJets[2];
  aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
  aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE

  for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
    fListJets[iJetType]->Clear();
    if (!aodJets[iJetType]) continue;
	
    for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
      AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
      if (jet) fListJets[iJetType]->Add(jet);
    }
    fListJets[iJetType]->Sort();
  }
  
  // jet matching 
  static TArrayI aMatchIndex(fListJets[0]->GetEntries());
  static TArrayF aPtFraction(fListJets[0]->GetEntries());
  if(aMatchIndex.GetSize()<fListJets[0]->GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries());
  if(aPtFraction.GetSize()<fListJets[0]->GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries());
  
  // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction'
  AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
                                            fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()),
                                            aMatchIndex, aPtFraction, fDebug);
											
  // loop over matched jets
  Int_t ir = -1; // index of matched reconstruced jet
  Float_t fraction = 0.;
  AliAODJet *jet[2]  = { 0x0, 0x0 };
  Float_t jetEta[2]  = { 0., 0. };
  Float_t jetPhi[2]  = { 0., 0. };
  Float_t jetPt[2]   = { 0., 0. };
  Float_t jetArea[2] = { 0., 0. };
   
  for(Int_t ig=0; ig<fListJets[0]->GetEntries(); ++ig){
     ir = aMatchIndex[ig];
	 if(ir<0) continue;
	 fraction = aPtFraction[ig];
	 
	 // fetch jets
	 jet[0] = (AliAODJet*)(fListJets[0]->At(ig));
	 jet[1] = (AliAODJet*)(fListJets[1]->At(ir));
	 if(!jet[0] || !jet[1]) continue;
	 
	 for(Int_t i=0; i<fkNbranches; ++i){
	    jetEta[i]  = jet[i]->Eta();
		jetPhi[i]  = jet[i]->Phi();
		jetPt[i]   = jet[i]->Pt();
		jetArea[i] = jet[i]->EffectiveAreaCharged();
	 }
	 if(eventClass > -1 && eventClass < fkEvtClasses){
        fHistPtFraction[eventClass] -> Fill(fraction, jetPt[0]);
	 }
	 
	 if(fraction<fJetPtFractionMin) continue;
	 
	 // calculate parameters of associated jets
	 Float_t deltaPt    = jetPt[1]-jetPt[0];
	 Float_t deltaEta   = jetEta[1]-jetEta[0];
	 Float_t deltaPhi   = TVector2::Phi_mpi_pi(jetPhi[1]-jetPhi[0]);
	 Float_t deltaR     = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi);
	 Float_t deltaArea  = jetArea[1]-jetArea[0];
	 
	 // fill histograms before acceptance cut
	 for(Int_t i=0; i<fkNbranches; ++i){
	     fHistEtaPhiJet[i] -> Fill(jetEta[i], jetPhi[i]);
         }
	 if(eventClass > -1 && eventClass < fkEvtClasses){
             fHistDeltaEtaDeltaPhiJet[eventClass] -> Fill(deltaEta, deltaPhi);
	 }
	 
	 // jet acceptance + minimum pT check
	 if(jetEta[0]>fJetEtaMax || jetEta[0]<fJetEtaMin ||
	    jetEta[1]>fJetEtaMax || jetEta[1]<fJetEtaMin){
		
		if(fDebug){
     		Printf("Jet not in eta acceptance.");
			Printf("[0]: jet %d eta %.2f", ig, jetEta[0]);
			Printf("[1]: jet %d eta %.2f", ir, jetEta[1]);
		}
		continue;
     }
	 if(jetPt[1] < fJetPtMin){
	    if(fDebug) Printf("Jet %d (pT %.1f GeV/c) has less than required pT.", ir, jetPt[1]);
	    continue;
	 }
	 
	 
	 // fill histograms
	 for(Int_t i=0; i<fkNbranches; ++i){
	     fHistPtJet[i]        -> Fill(jetPt[i]);
             fHistEtaPhiJetCut[i] -> Fill(jetEta[i], jetPhi[i]);
	 }
	 
	 fHistPtPtExtra->Fill(jetPt[0], jetPt[1]);
	 
	 if(eventClass > -1 && eventClass < fkEvtClasses){
		 fHistDeltaEtaDeltaPhiJetCut[eventClass] -> Fill(deltaEta, deltaPhi);
		 
		 fHistPtResponse[eventClass]             -> Fill(jetPt[1], jetPt[0]);
		 fHistPtSmearing[eventClass]             -> Fill(deltaPt,  jetPt[0]);
		 
		 fHistDeltaPtEtaJet[eventClass]          -> Fill(jetEta[0], deltaPt);
                 fHistDeltaEtaEtaJet[eventClass]         -> Fill(jetEta[0], deltaEta);
		 
		 fHistDeltaR[eventClass]                 -> Fill(deltaPt, deltaR);
		 fHistArea[eventClass]                   -> Fill(deltaPt, jetArea[1]);
                 fHistDeltaArea[eventClass]              -> Fill(deltaPt, deltaArea);
		 
	 }
  }

  PostData(1, fOutputList);
}

void AliAnalysisTaskJetResponse::Terminate(const Option_t *)
{
  // Draw result to the screen
  // Called once at the end of the query

  if (!GetOutputData(1))
    return;
}
Commit	Line	Data
2a436fa1	1	#include "TChain.h"
	2	#include "TTree.h"
	3	#include "TMath.h"
	4	#include "TH1F.h"
	5	#include "TH2F.h"
46465e39	6	#include "TH3F.h"
2a436fa1	7	#include "TCanvas.h"
	8
	9	#include "AliLog.h"
	10
	11	#include "AliAnalysisTask.h"
	12	#include "AliAnalysisManager.h"
	13
	14	#include "AliVEvent.h"
	15	#include "AliESDEvent.h"
	16	#include "AliESDInputHandler.h"
	17	#include "AliCentrality.h"
	18	#include "AliAnalysisHelperJetTasks.h"
	19	#include "AliInputEventHandler.h"
	20
	21	#include "AliAODEvent.h"
	22	#include "AliAODJet.h"
	23
	24	#include "AliAnalysisTaskJetResponse.h"
	25
	26	ClassImp(AliAnalysisTaskJetResponse)
	27
	28	AliAnalysisTaskJetResponse::AliAnalysisTaskJetResponse() :
	29	AliAnalysisTaskSE(),
	30	fESD(0x0),
	31	fAOD(0x0),
	32	fOfflineTrgMask(AliVEvent::kAny),
	33	fMinContribVtx(1),
	34	fVtxZMin(-8.),
	35	fVtxZMax(8.),
46465e39	36	fEvtClassMin(1),
46465e39	37	fEvtClassMax(4),
2a436fa1	38	fCentMin(0.),
	39	fCentMax(100.),
	40	fJetEtaMin(-.5),
	41	fJetEtaMax(.5),
46465e39	42	fJetPtMin(20.),
	43	fJetPtFractionMin(0.5),
	44	fNMatchJets(4),
	45	//fJetDeltaEta(.4),
	46	//fJetDeltaPhi(.4),
2a436fa1	47	fkNbranches(2),
46465e39	48	fkEvtClasses(5),
2a436fa1	49	fOutputList(0x0),
2a436fa1	50	fHistEvtSelection(0x0),
46465e39	51	fHistEvtClass(0x0),
	52	fHistCentrality(0x0),
	53	fHistPtJet(new TH1F*[fkNbranches]),
	54	fHistEtaPhiJet(new TH2F*[fkNbranches]),
	55	fHistEtaPhiJetCut(new TH2F*[fkNbranches]),
	56	fHistDeltaEtaDeltaPhiJet(new TH2F*[fkEvtClasses]),
	57	fHistDeltaEtaDeltaPhiJetCut(new TH2F*[fkEvtClasses]),
	58	fHistDeltaEtaDeltaPhiJetNOMatching(new TH2F*[fkEvtClasses]),
	59	fHistDeltaEtaEtaJet(new TH2F*[fkEvtClasses]),
	60	fHistDeltaPtEtaJet(new TH2F*[fkEvtClasses]),
	61	fHistPtFraction(new TH2F*[fkEvtClasses]),
2a436fa1	62	fHistPtPtExtra(0x0),
2a436fa1	63	fHistPtResponse(new TH2F*[fkEvtClasses]),
a11972e9	64	fHistPtSmearing(new TH2F*[fkEvtClasses]),
46465e39	65	fHistDeltaR(new TH2F*[fkEvtClasses]),
	66	fHistArea(new TH2F*[fkEvtClasses]),
	67	fHistDeltaArea(new TH2F*[fkEvtClasses])
2a436fa1	68	{
	69	// default Constructor
	70
	71	fJetBranchName[0] = "";
	72	fJetBranchName[1] = "";
	73
	74	fListJets[0] = new TList;
	75	fListJets[1] = new TList;
	76	}
	77
	78	AliAnalysisTaskJetResponse::AliAnalysisTaskJetResponse(const char *name) :
	79	AliAnalysisTaskSE(name),
	80	fESD(0x0),
	81	fAOD(0x0),
	82	fOfflineTrgMask(AliVEvent::kAny),
	83	fMinContribVtx(1),
	84	fVtxZMin(-8.),
	85	fVtxZMax(8.),
46465e39	86	fEvtClassMin(1),
46465e39	87	fEvtClassMax(4),
2a436fa1	88	fCentMin(0.),
	89	fCentMax(100.),
	90	fJetEtaMin(-.5),
	91	fJetEtaMax(.5),
46465e39	92	fJetPtMin(20.),
	93	fJetPtFractionMin(0.5),
	94	fNMatchJets(4),
	95	//fJetDeltaEta(.4),
	96	//fJetDeltaPhi(.4),
2a436fa1	97	fkNbranches(2),
46465e39	98	fkEvtClasses(5),
2a436fa1	99	fOutputList(0x0),
2a436fa1	100	fHistEvtSelection(0x0),
46465e39	101	fHistEvtClass(0x0),
	102	fHistCentrality(0x0),
	103	fHistPtJet(new TH1F*[fkNbranches]),
	104	fHistEtaPhiJet(new TH2F*[fkNbranches]),
	105	fHistEtaPhiJetCut(new TH2F*[fkNbranches]),
	106	fHistDeltaEtaDeltaPhiJet(new TH2F*[fkEvtClasses]),
	107	fHistDeltaEtaDeltaPhiJetCut(new TH2F*[fkEvtClasses]),
	108	fHistDeltaEtaDeltaPhiJetNOMatching(new TH2F*[fkEvtClasses]),
	109	fHistDeltaEtaEtaJet(new TH2F*[fkEvtClasses]),
	110	fHistDeltaPtEtaJet(new TH2F*[fkEvtClasses]),
	111	fHistPtFraction(new TH2F*[fkEvtClasses]),
2a436fa1	112	fHistPtPtExtra(0x0),
2a436fa1	113	fHistPtResponse(new TH2F*[fkEvtClasses]),
a11972e9	114	fHistPtSmearing(new TH2F*[fkEvtClasses]),
46465e39	115	fHistDeltaR(new TH2F*[fkEvtClasses]),
	116	fHistArea(new TH2F*[fkEvtClasses]),
	117	fHistDeltaArea(new TH2F*[fkEvtClasses])
2a436fa1	118	{
	119	// Constructor
	120
	121	fJetBranchName[0] = "";
	122	fJetBranchName[1] = "";
	123
	124	fListJets[0] = new TList;
	125	fListJets[1] = new TList;
	126
	127	DefineOutput(1, TList::Class());
	128	}
	129
	130	AliAnalysisTaskJetResponse::~AliAnalysisTaskJetResponse()
	131	{
	132	delete fListJets[0];
	133	delete fListJets[1];
	134	}
	135
	136	void AliAnalysisTaskJetResponse::SetBranchNames(const TString &branch1, const TString &branch2)
	137	{
	138	fJetBranchName[0] = branch1;
	139	fJetBranchName[1] = branch2;
	140	}
	141
	142	void AliAnalysisTaskJetResponse::Init()
	143	{
	144
	145	// check for jet branches
	146	if(!strlen(fJetBranchName[0].Data()) \|\| !strlen(fJetBranchName[1].Data())){
	147	AliError("Jet branch name not set.");
	148	}
	149
	150	}
	151
	152	void AliAnalysisTaskJetResponse::UserCreateOutputObjects()
	153	{
	154	// Create histograms
	155	// Called once
46465e39	156	OpenFile(1);
	157	if(!fOutputList) fOutputList = new TList;
	158	fOutputList->SetOwner(kTRUE);
	159
	160	Bool_t oldStatus = TH1::AddDirectoryStatus();
	161	TH1::AddDirectory(kFALSE);
	162
2a436fa1	163
46465e39	164	fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 5, -0.5, 4.5);
2a436fa1	165	fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
	166	fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
	167	fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
	168	fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
	169	fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
46465e39	170
	171	fHistEvtClass = new TH1I("fHistEvtClass", "event classes", fkEvtClasses, -0.5, fkEvtClasses-0.5);
	172	fHistCentrality = new TH1F("fHistCentrality", "event centrality", 100, 0., 100.);
2a436fa1	173
2a436fa1	174	for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
46465e39	175	fHistPtJet[iBranch] = new TH1F(Form("fHistPtJet%i", iBranch),
	176	Form("p_{T} of jets from branch %i;p_{T} (GeV/c);counts", iBranch),
	177	250, 0., 250.);
	178	fHistPtJet[iBranch]->SetMarkerStyle(kFullCircle);
	179
	180	fHistEtaPhiJet[iBranch] = new TH2F(Form("fHistEtaPhiJet%i", iBranch),
	181	Form("#eta - #phi of jets from branch %i (before cut);#eta;#phi", iBranch),
	182	48, -1.2, 1.2, 100, 0., 2*TMath::Pi());
	183	fHistEtaPhiJetCut[iBranch] = new TH2F(Form("fHistEtaPhiJetCut%i", iBranch),
	184	Form("#eta - #phi of jets from branch %i (before cut);#eta;#phi", iBranch),
	185	48, -1.2, 1.2, 100, 0., 2*TMath::Pi());
2a436fa1	186
2a436fa1	187	}
46465e39	188
2a436fa1	189
2a436fa1	190	fHistPtPtExtra = new TH2F("fHistPtPtExtra", "p_{T} response;p_{T} (GeV/c);p_{T} (GeV/c)",
46465e39	191	250, 0., 250., 250, 0., 250.);
2a436fa1	192
2a436fa1	193	for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
46465e39	194	fHistDeltaEtaDeltaPhiJet[iEvtClass] = new TH2F(Form("fHistDeltaEtaDeltaPhiJet%i",iEvtClass),
	195	"#Delta#eta - #Delta#phi of jet;#Delta#eta;#Delta#phi",
	196	101, -1.01, 1.01, 101, -1.01, 1.01);
	197	fHistDeltaEtaDeltaPhiJetCut[iEvtClass] = new TH2F(Form("fHistDeltaEtaDeltaPhiJetCut%i",iEvtClass),
	198	"#Delta#eta - #Delta#phi of jet;#Delta#eta;#Delta#phi",
	199	101, -1.01, 1.01, 101, -1.01, 1.01);
	200	fHistDeltaEtaDeltaPhiJetNOMatching[iEvtClass] = new TH2F("fHistDeltaEtaDeltaPhiJetNOMatching",
	201	"#Delta#eta - #Delta#phi of jets which do not match;#Delta#eta;#Delta#phi",
	202	100, -2., 2., 100, TMath::Pi(), TMath::Pi());
2a436fa1	203	fHistPtResponse[iEvtClass] = new TH2F(Form("pt_response%i",iEvtClass), Form("pt_response%i;p_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
46465e39	204	250, 0., 250., 250, 0., 250.);
	205	fHistPtSmearing[iEvtClass] = new TH2F(Form("pt_smearing%i",iEvtClass), Form("pt_smearing%i;#Deltap_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
	206	200, -50., 150., 250, 0., 250.);
	207	fHistDeltaR[iEvtClass] = new TH2F(Form("hist_DeltaR%i",iEvtClass), "#DeltaR of matched jets;#Deltap_{T};#DeltaR", 200, -50.,150., 60, 0.,.6);
	208	fHistArea[iEvtClass] = new TH2F(Form("hist_Area%i",iEvtClass), "jet area;#Deltap_{T};jet area", 200, -50.,150., 100, 0.,1.0);
	209	fHistDeltaArea[iEvtClass] = new TH2F(Form("hist_DeltaArea%i",iEvtClass), "#DeltaArea of matched jets;#Deltap_{T};#Deltaarea", 200, -50., 150., 60, 0., .3);
	210
	211	fHistDeltaEtaEtaJet[iEvtClass] = new TH2F(Form("fHistDeltaEtaEtaJet%i", iEvtClass),
	212	Form("#eta - #Delta#eta of matched jets from event class %i;#eta;#Delta#eta", iEvtClass),
	213	60, -.6, .6, 100, -.5, .5);
	214	fHistDeltaPtEtaJet[iEvtClass] = new TH2F(Form("fHistDeltaPtEtaJet%i", iEvtClass),
	215	Form("#eta - #Deltap_{T} of matched jets from event class %i;#eta;#Deltap_{T}", iEvtClass),
	216	60, -.6, .6, 200, -50., 150);
	217	fHistPtFraction[iEvtClass] = new TH2F(Form("fHistPtFraction%i", iEvtClass),
	218	Form("fraction from embedded jet in reconstructed jet;fraction (event class %i);p_{T}^{emb}", iEvtClass),
	219	100, 0., 1., 250, 0., 250.);
2a436fa1	220	}
2a436fa1	221
46465e39	222
2a436fa1	223	fOutputList->Add(fHistEvtSelection);
46465e39	224	fOutputList->Add(fHistEvtClass);
	225	fOutputList->Add(fHistCentrality);
	226
2a436fa1	227	for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
46465e39	228	fOutputList->Add(fHistPtJet[iBranch]);
	229	fOutputList->Add(fHistEtaPhiJet[iBranch]);
	230	fOutputList->Add(fHistEtaPhiJetCut[iBranch]);
2a436fa1	231	}
46465e39	232
2a436fa1	233	fOutputList->Add(fHistPtPtExtra);
46465e39	234
2a436fa1	235	for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
46465e39	236	fOutputList->Add(fHistDeltaEtaDeltaPhiJet[iEvtClass]);
	237	fOutputList->Add(fHistDeltaEtaDeltaPhiJetCut[iEvtClass]);
	238	fOutputList->Add(fHistDeltaEtaDeltaPhiJetNOMatching[iEvtClass]);
	239	fOutputList->Add(fHistDeltaEtaEtaJet[iEvtClass]);
	240	fOutputList->Add(fHistDeltaPtEtaJet[iEvtClass]);
	241	fOutputList->Add(fHistPtFraction[iEvtClass]);
	242
2a436fa1	243	fOutputList->Add(fHistPtResponse[iEvtClass]);
2a436fa1	244	fOutputList->Add(fHistPtSmearing[iEvtClass]);
46465e39	245	fOutputList->Add(fHistDeltaR[iEvtClass]);
a11972e9	246	fOutputList->Add(fHistArea[iEvtClass]);
46465e39	247	fOutputList->Add(fHistDeltaArea[iEvtClass]);
2a436fa1	248	}
46465e39	249
	250	// =========== Switch on Sumw2 for all histos ===========
	251	for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
	252	TH1 h1 = dynamic_cast<TH1>(fOutputList->At(i));
	253	if (h1){
	254	h1->Sumw2();
	255	continue;
	256	}
	257	}
	258	TH1::AddDirectory(oldStatus);
	259
	260	PostData(1, fOutputList);
2a436fa1	261	}
	262
	263	void AliAnalysisTaskJetResponse::UserExec(Option_t *)
	264	{
	265	// load events, apply event cuts, then compare leading jets
	266
	267	if(!strlen(fJetBranchName[0].Data()) \|\| !strlen(fJetBranchName[1].Data())){
	268	AliError("Jet branch name not set.");
	269	return;
	270	}
	271
	272	fESD=dynamic_cast<AliESDEvent*>(InputEvent());
	273	if (!fESD) {
	274	AliError("ESD not available");
	275	return;
	276	}
	277	fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
	278	if (!fAOD) {
	279	AliError("AOD not available");
	280	return;
	281	}
	282
46465e39	283	// -- event selection --
2a436fa1	284	fHistEvtSelection->Fill(1); // number of events before event selection
	285
	286	// physics selection
	287	AliInputEventHandler* inputHandler = (AliInputEventHandler*)
	288	((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
	289	if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
46465e39	290	if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
2a436fa1	291	fHistEvtSelection->Fill(2);
	292	PostData(1, fOutputList);
	293	return;
	294	}
	295
	296	// vertex selection
	297	AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
	298	Int_t nTracksPrim = primVtx->GetNContributors();
	299	if ((nTracksPrim < fMinContribVtx) \|\|
	300	(primVtx->GetZ() < fVtxZMin) \|\|
	301	(primVtx->GetZ() > fVtxZMax) ){
46465e39	302	if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
2a436fa1	303	fHistEvtSelection->Fill(3);
	304	PostData(1, fOutputList);
	305	return;
	306	}
	307
	308	// event class selection (from jet helper task)
	309	Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
46465e39	310	if(fDebug) Printf("Event class %d", eventClass);
2a436fa1	311	if (eventClass < fEvtClassMin \|\| eventClass > fEvtClassMax){
	312	fHistEvtSelection->Fill(4);
	313	PostData(1, fOutputList);
	314	return;
	315	}
	316
	317	// centrality selection
	318	AliCentrality *cent = fESD->GetCentrality();
	319	Float_t centValue = cent->GetCentralityPercentile("TRK");
46465e39	320	if(fDebug) printf("centrality: %f\n", centValue);
2a436fa1	321	if (centValue < fCentMin \|\| centValue > fCentMax){
46465e39	322	fHistEvtSelection->Fill(4);
2a436fa1	323	PostData(1, fOutputList);
	324	return;
	325	}
	326
	327	fHistEvtSelection->Fill(0); // accepted events
46465e39	328	fHistEvtClass->Fill(eventClass);
	329	fHistCentrality->Fill(centValue);
	330	// -- end event selection --
2a436fa1	331
46465e39	332	// fetch jets
2a436fa1	333	TClonesArray *aodJets[2];
a11972e9	334	aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
a11972e9	335	aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE
2a436fa1	336
	337	for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
	338	fListJets[iJetType]->Clear();
46465e39	339	if (!aodJets[iJetType]) continue;
46465e39	340
2a436fa1	341	for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
2a436fa1	342	AliAODJet jet = dynamic_cast<AliAODJet>((*aodJets[iJetType])[iJet]);
46465e39	343	if (jet) fListJets[iJetType]->Add(jet);
2a436fa1	344	}
2a436fa1	345	fListJets[iJetType]->Sort();
2a436fa1	346	}
46465e39	347
	348	// jet matching
	349	static TArrayI aMatchIndex(fListJets[0]->GetEntries());
	350	static TArrayF aPtFraction(fListJets[0]->GetEntries());
	351	if(aMatchIndex.GetSize()<fListJets[0]->GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries());
	352	if(aPtFraction.GetSize()<fListJets[0]->GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries());
	353
	354	// stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction'
	355	AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
	356	fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()),
	357	aMatchIndex, aPtFraction, fDebug);
	358
	359	// loop over matched jets
	360	Int_t ir = -1; // index of matched reconstruced jet
	361	Float_t fraction = 0.;
	362	AliAODJet *jet[2] = { 0x0, 0x0 };
	363	Float_t jetEta[2] = { 0., 0. };
	364	Float_t jetPhi[2] = { 0., 0. };
	365	Float_t jetPt[2] = { 0., 0. };
	366	Float_t jetArea[2] = { 0., 0. };
	367
	368	for(Int_t ig=0; ig<fListJets[0]->GetEntries(); ++ig){
	369	ir = aMatchIndex[ig];
	370	if(ir<0) continue;
	371	fraction = aPtFraction[ig];
	372
	373	// fetch jets
	374	jet[0] = (AliAODJet*)(fListJets[0]->At(ig));
	375	jet[1] = (AliAODJet*)(fListJets[1]->At(ir));
	376	if(!jet[0] \|\| !jet[1]) continue;
	377
	378	for(Int_t i=0; i<fkNbranches; ++i){
	379	jetEta[i] = jet[i]->Eta();
	380	jetPhi[i] = jet[i]->Phi();
	381	jetPt[i] = jet[i]->Pt();
	382	jetArea[i] = jet[i]->EffectiveAreaCharged();
	383	}
	384	if(eventClass > -1 && eventClass < fkEvtClasses){
	385	fHistPtFraction[eventClass] -> Fill(fraction, jetPt[0]);
	386	}
	387
	388	if(fraction<fJetPtFractionMin) continue;
	389
	390	// calculate parameters of associated jets
	391	Float_t deltaPt = jetPt[1]-jetPt[0];
	392	Float_t deltaEta = jetEta[1]-jetEta[0];
	393	Float_t deltaPhi = TVector2::Phi_mpi_pi(jetPhi[1]-jetPhi[0]);
	394	Float_t deltaR = TMath::Sqrt(deltaEtadeltaEta+deltaPhideltaPhi);
	395	Float_t deltaArea = jetArea[1]-jetArea[0];
	396
	397	// fill histograms before acceptance cut
	398	for(Int_t i=0; i<fkNbranches; ++i){
	399	fHistEtaPhiJet[i] -> Fill(jetEta[i], jetPhi[i]);
	400	}
	401	if(eventClass > -1 && eventClass < fkEvtClasses){
	402	fHistDeltaEtaDeltaPhiJet[eventClass] -> Fill(deltaEta, deltaPhi);
	403	}
	404
	405	// jet acceptance + minimum pT check
	406	if(jetEta[0]>fJetEtaMax \|\| jetEta[0]<fJetEtaMin \|\|
	407	jetEta[1]>fJetEtaMax \|\| jetEta[1]<fJetEtaMin){
	408
	409	if(fDebug){
	410	Printf("Jet not in eta acceptance.");
411	Printf("[0]: jet %d eta %.2f", ig, jetEta[0]);
412	Printf("[1]: jet %d eta %.2f", ir, jetEta[1]);
413	}
414	continue;
415	}
416	if(jetPt[1] < fJetPtMin){
417	if(fDebug) Printf("Jet %d (pT %.1f GeV/c) has less than required pT.", ir, jetPt[1]);
418	continue;
419	}
420
421
422
423	// fill histograms
424	for(Int_t i=0; i<fkNbranches; ++i){
425	fHistPtJet[i] -> Fill(jetPt[i]);
426	fHistEtaPhiJetCut[i] -> Fill(jetEta[i], jetPhi[i]);
427	}
428
429	fHistPtPtExtra->Fill(jetPt[0], jetPt[1]);
430
431	if(eventClass > -1 && eventClass < fkEvtClasses){
432	fHistDeltaEtaDeltaPhiJetCut[eventClass] -> Fill(deltaEta, deltaPhi);
433
434	fHistPtResponse[eventClass] -> Fill(jetPt[1], jetPt[0]);
435	fHistPtSmearing[eventClass] -> Fill(deltaPt, jetPt[0]);
436
437	fHistDeltaPtEtaJet[eventClass] -> Fill(jetEta[0], deltaPt);
438	fHistDeltaEtaEtaJet[eventClass] -> Fill(jetEta[0], deltaEta);
439
440	fHistDeltaR[eventClass] -> Fill(deltaPt, deltaR);
441	fHistArea[eventClass] -> Fill(deltaPt, jetArea[1]);
442	fHistDeltaArea[eventClass] -> Fill(deltaPt, deltaArea);
443
444	}
2a436fa1	445	}
	446
	447	PostData(1, fOutputList);
	448	}
	449
	450	void AliAnalysisTaskJetResponse::Terminate(const Option_t *)
	451	{
	452	// Draw result to the screen
	453	// Called once at the end of the query
	454
	455	if (!GetOutputData(1))
	456	return;
	457	}
46465e39	458