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

#include "TChain.h"
#include "TTree.h"
#include "TMath.h"
#include "TH1F.h"
#include "TH2F.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(0),
  fEvtClassMax(10),
  fCentMin(0.),
  fCentMax(100.),
  fJetEtaMin(-.5),
  fJetEtaMax(.5),
  fJetDeltaEta(.4),
  fJetDeltaPhi(.4),
  fkNbranches(2),
  fkEvtClasses(10),
  fOutputList(0x0),
  fHistEvtSelection(0x0),
  fHistPtLeadingJet(new TH1F*[fkNbranches]),
  fHistEtaPhiLeadingJet(new TH2F*[fkNbranches]),
  fHistEtaPhiLeadingJetCut(new TH2F*[fkEvtClasses]),
  fHistDeltaEtaDeltaPhiLeadingJet(0x0),
  fHistDeltaEtaEtaLeadingJet(new TH2F*[fkEvtClasses]),
  fHistDeltaPtEtaLeadingJet(new TH2F*[fkEvtClasses]),
  fHistPtPtExtra(0x0),
  fHistPtResponse(new TH2F*[fkEvtClasses]),
  fHistPtSmearing(new TH2F*[fkEvtClasses]),
  fHistdR(new TH2F*[fkEvtClasses]),
  fHistArea(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(0),
  fEvtClassMax(10),
  fCentMin(0.),
  fCentMax(100.),
  fJetEtaMin(-.5),
  fJetEtaMax(.5),
  fJetDeltaEta(.4),
  fJetDeltaPhi(.4),
  fkNbranches(2),
  fkEvtClasses(10),
  fOutputList(0x0),
  fHistEvtSelection(0x0),
  fHistPtLeadingJet(new TH1F*[fkNbranches]),
  fHistEtaPhiLeadingJet(new TH2F*[fkNbranches]),
  fHistEtaPhiLeadingJetCut(new TH2F*[fkEvtClasses]),
  fHistDeltaEtaDeltaPhiLeadingJet(0x0),
  fHistDeltaEtaEtaLeadingJet(new TH2F*[fkEvtClasses]),
  fHistDeltaPtEtaLeadingJet(new TH2F*[fkEvtClasses]),
  fHistPtPtExtra(0x0),
  fHistPtResponse(new TH2F*[fkEvtClasses]),
  fHistPtSmearing(new TH2F*[fkEvtClasses]),
  fHistdR(new TH2F*[fkEvtClasses]),
  fHistArea(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

  fHistEvtSelection = new TH1F("fHistEvtSelection", "event selection", 5, -0.5, 5.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)");

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

    fHistEtaPhiLeadingJet[iBranch] = new TH2F(Form("fHistEtaPhiLeadingJet%i", iBranch),
					      Form("#eta - #phi of leading jet from branch %i;#eta;#phi", iBranch),
					      100, -2., 2., 100, 0., 2*TMath::Pi());

  }
  fHistDeltaEtaDeltaPhiLeadingJet = new TH2F("fHistDeltaEtaDeltaPhiLeadingJet",
					     "#Delta#eta - #Delta#phi of leading jet;#Delta#eta;#Delta#phi",
					     100, -4., 4., 100, -2.*TMath::Pi(), 2*TMath::Pi());

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

  for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
    fHistPtResponse[iEvtClass] = new TH2F(Form("pt_response%i",iEvtClass), Form("pt_response%i;p_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
					  300, 0., 300., 300, 0., 300.);
    fHistPtSmearing[iEvtClass] = new TH2F(Form("pt_smearing%i",iEvtClass), Form("pt_smearing%i;p_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
					  200, -50., 150., 300, 0., 300.);
    fHistdR[iEvtClass] = new TH2F(Form("hist_dR%i",iEvtClass), "", 200, -50.,150., 200, -1.,1.);
    fHistArea[iEvtClass] = new TH2F(Form("hist_Area%i",iEvtClass), "", 200, -50.,150., 100, 0.,1.);

    fHistEtaPhiLeadingJetCut[iEvtClass] = new TH2F(Form("fHistEtaPhiLeadingJetCut%i", iEvtClass),
					      Form("#eta - #phi of leading jet from event class %i;#eta;#phi", iEvtClass),
					      100, -2., 2., 100, 0., 2*TMath::Pi());
    fHistDeltaEtaEtaLeadingJet[iEvtClass] = new TH2F(Form("fHistDeltaEtaEtaLeadingJet%i", iEvtClass),
                                                Form("#eta - #Delta#eta of leading jet from event class %i;#eta;#Delta#eta", iEvtClass),
                                                100, -1., 1., 100, -.5, .5);
    fHistDeltaPtEtaLeadingJet[iEvtClass] = new TH2F(Form("fHistDeltaPtEtaLeadingJet%i", iEvtClass),
                                                Form("#eta - #Delta#eta of leading jet from event class %i;#eta;#Delta#p_{T}", iEvtClass),
                                                100, -1., 1., 200, -50., 150);
  }

  OpenFile(1);
  fOutputList = new TList;
  fOutputList->Add(fHistEvtSelection);
  for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
    fOutputList->Add(fHistPtLeadingJet[iBranch]);
    fOutputList->Add(fHistEtaPhiLeadingJet[iBranch]);
  }
  fOutputList->Add(fHistDeltaEtaDeltaPhiLeadingJet);
  fOutputList->Add(fHistPtPtExtra);
  for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
    fOutputList->Add(fHistPtResponse[iEvtClass]);
    fOutputList->Add(fHistPtSmearing[iEvtClass]);
    fOutputList->Add(fHistdR[iEvtClass]);
    fOutputList->Add(fHistArea[iEvtClass]);
    fOutputList->Add(fHistEtaPhiLeadingJetCut[iEvtClass]);
    fOutputList->Add(fHistDeltaEtaEtaLeadingJet[iEvtClass]);
    fOutputList->Add(fHistDeltaPtEtaLeadingJet[iEvtClass]);
  }
}

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

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

  // physics selection
  AliInputEventHandler* inputHandler = (AliInputEventHandler*)
    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
  if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
    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) ){
    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();
  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");
  printf("centrality: %f\n", centValue);
  if (centValue < fCentMin || centValue > fCentMax){
    fHistEvtSelection->Fill(5);
    PostData(1, fOutputList);
    return;
  }

  fHistEvtSelection->Fill(0); // accepted events


  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
  AliAODJet *leadingJet[2] = { 0x0, 0x0 };

  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();
    leadingJet[iJetType] = (AliAODJet*) fListJets[iJetType]->First();

    if (leadingJet[iJetType])
      fHistEtaPhiLeadingJet[iJetType]->Fill(leadingJet[iJetType]->Eta(), leadingJet[iJetType]->Phi());
  }

  // check if leading jets are close in eta-phi and compare their Pt
  if (leadingJet[0] && leadingJet[1]) {
    Float_t deltaEta = leadingJet[0]->Eta() - leadingJet[1]->Eta();
    Float_t deltaPhi = TVector2::Phi_mpi_pi(leadingJet[0]->Phi() - leadingJet[1]->Phi());

    if (eventClass > -1 && eventClass < fkEvtClasses){

       if(leadingJet[1]->Eta()<fJetEtaMax && leadingJet[1]->Eta()>fJetEtaMin){
          fHistEtaPhiLeadingJetCut[eventClass]->Fill(leadingJet[1]->Eta(), leadingJet[1]->Phi());
       }
    }

    // leading jets in "jet" acceptance
    if(leadingJet[0]->Eta()>fJetEtaMax || leadingJet[0]->Eta()<fJetEtaMin ||
       leadingJet[1]->Eta()>fJetEtaMax || leadingJet[1]->Eta()<fJetEtaMin){
       Printf("Jet not in eta acceptance.");
    }
    else{
       // check association of jets
       fHistDeltaEtaDeltaPhiLeadingJet->Fill(deltaEta, deltaPhi);

       if (TMath::Abs(deltaEta) > fJetDeltaEta || (TMath::Pi() - TMath::Abs(TMath::Abs(deltaPhi) - TMath::Pi())) > fJetDeltaPhi)
          printf("leading jets two far apart\n");
       else {
          Float_t pt0 = leadingJet[0]->Pt();
          Float_t pt1 = leadingJet[1]->Pt();
          Float_t dPt = pt1-pt0;
          Float_t jetarea = leadingJet[1]->EffectiveAreaCharged();

          fHistPtLeadingJet[0]->Fill(pt0);
          fHistPtLeadingJet[1]->Fill(pt1);

          fHistPtPtExtra->Fill(pt0, pt1);

          if (eventClass > -1 && eventClass < fkEvtClasses){
	      fHistPtResponse[eventClass]->Fill(pt1, pt0);
              fHistPtSmearing[eventClass]->Fill(dPt, pt0);

              Float_t dR = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi);
              fHistdR[eventClass]->Fill(dPt, dR);
              fHistArea[eventClass]->Fill(dPt, jetarea);

              fHistDeltaEtaEtaLeadingJet[eventClass]->Fill(leadingJet[0]->Eta(), deltaEta);
              fHistDeltaPtEtaLeadingJet[eventClass]->Fill(leadingJet[0]->Eta(), dPt);
          }
       }
    }
  }

  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"
	6	#include "TCanvas.h"
	7
	8	#include "AliLog.h"
	9
	10	#include "AliAnalysisTask.h"
	11	#include "AliAnalysisManager.h"
	12
	13	#include "AliVEvent.h"
	14	#include "AliESDEvent.h"
	15	#include "AliESDInputHandler.h"
	16	#include "AliCentrality.h"
	17	#include "AliAnalysisHelperJetTasks.h"
	18	#include "AliInputEventHandler.h"
	19
	20	#include "AliAODEvent.h"
	21	#include "AliAODJet.h"
	22
	23	#include "AliAnalysisTaskJetResponse.h"
	24
	25	ClassImp(AliAnalysisTaskJetResponse)
	26
	27	AliAnalysisTaskJetResponse::AliAnalysisTaskJetResponse() :
	28	AliAnalysisTaskSE(),
	29	fESD(0x0),
	30	fAOD(0x0),
	31	fOfflineTrgMask(AliVEvent::kAny),
	32	fMinContribVtx(1),
	33	fVtxZMin(-8.),
	34	fVtxZMax(8.),
	35	fEvtClassMin(0),
	36	fEvtClassMax(10),
	37	fCentMin(0.),
	38	fCentMax(100.),
	39	fJetEtaMin(-.5),
	40	fJetEtaMax(.5),
	41	fJetDeltaEta(.4),
	42	fJetDeltaPhi(.4),
	43	fkNbranches(2),
	44	fkEvtClasses(10),
	45	fOutputList(0x0),
	46	fHistEvtSelection(0x0),
	47	fHistPtLeadingJet(new TH1F*[fkNbranches]),
	48	fHistEtaPhiLeadingJet(new TH2F*[fkNbranches]),
a11972e9	49	fHistEtaPhiLeadingJetCut(new TH2F*[fkEvtClasses]),
2a436fa1	50	fHistDeltaEtaDeltaPhiLeadingJet(0x0),
a11972e9	51	fHistDeltaEtaEtaLeadingJet(new TH2F*[fkEvtClasses]),
a11972e9	52	fHistDeltaPtEtaLeadingJet(new TH2F*[fkEvtClasses]),
2a436fa1	53	fHistPtPtExtra(0x0),
2a436fa1	54	fHistPtResponse(new TH2F*[fkEvtClasses]),
a11972e9	55	fHistPtSmearing(new TH2F*[fkEvtClasses]),
	56	fHistdR(new TH2F*[fkEvtClasses]),
	57	fHistArea(new TH2F*[fkEvtClasses])
2a436fa1	58	{
	59	// default Constructor
	60
	61	fJetBranchName[0] = "";
	62	fJetBranchName[1] = "";
	63
	64	fListJets[0] = new TList;
	65	fListJets[1] = new TList;
	66	}
	67
	68	AliAnalysisTaskJetResponse::AliAnalysisTaskJetResponse(const char *name) :
	69	AliAnalysisTaskSE(name),
	70	fESD(0x0),
	71	fAOD(0x0),
	72	fOfflineTrgMask(AliVEvent::kAny),
	73	fMinContribVtx(1),
	74	fVtxZMin(-8.),
	75	fVtxZMax(8.),
	76	fEvtClassMin(0),
	77	fEvtClassMax(10),
	78	fCentMin(0.),
	79	fCentMax(100.),
	80	fJetEtaMin(-.5),
	81	fJetEtaMax(.5),
	82	fJetDeltaEta(.4),
	83	fJetDeltaPhi(.4),
	84	fkNbranches(2),
	85	fkEvtClasses(10),
	86	fOutputList(0x0),
	87	fHistEvtSelection(0x0),
	88	fHistPtLeadingJet(new TH1F*[fkNbranches]),
	89	fHistEtaPhiLeadingJet(new TH2F*[fkNbranches]),
a11972e9	90	fHistEtaPhiLeadingJetCut(new TH2F*[fkEvtClasses]),
2a436fa1	91	fHistDeltaEtaDeltaPhiLeadingJet(0x0),
a11972e9	92	fHistDeltaEtaEtaLeadingJet(new TH2F*[fkEvtClasses]),
a11972e9	93	fHistDeltaPtEtaLeadingJet(new TH2F*[fkEvtClasses]),
2a436fa1	94	fHistPtPtExtra(0x0),
2a436fa1	95	fHistPtResponse(new TH2F*[fkEvtClasses]),
a11972e9	96	fHistPtSmearing(new TH2F*[fkEvtClasses]),
	97	fHistdR(new TH2F*[fkEvtClasses]),
	98	fHistArea(new TH2F*[fkEvtClasses])
2a436fa1	99	{
	100	// Constructor
	101
	102	fJetBranchName[0] = "";
	103	fJetBranchName[1] = "";
	104
	105	fListJets[0] = new TList;
	106	fListJets[1] = new TList;
	107
	108	DefineOutput(1, TList::Class());
	109	}
	110
	111	AliAnalysisTaskJetResponse::~AliAnalysisTaskJetResponse()
	112	{
	113	delete fListJets[0];
	114	delete fListJets[1];
	115	}
	116
	117	void AliAnalysisTaskJetResponse::SetBranchNames(const TString &branch1, const TString &branch2)
	118	{
	119	fJetBranchName[0] = branch1;
	120	fJetBranchName[1] = branch2;
	121	}
	122
	123	void AliAnalysisTaskJetResponse::Init()
	124	{
	125
	126	// check for jet branches
	127	if(!strlen(fJetBranchName[0].Data()) \|\| !strlen(fJetBranchName[1].Data())){
	128	AliError("Jet branch name not set.");
	129	}
	130
	131	}
	132
	133	void AliAnalysisTaskJetResponse::UserCreateOutputObjects()
	134	{
	135	// Create histograms
	136	// Called once
	137
	138	fHistEvtSelection = new TH1F("fHistEvtSelection", "event selection", 5, -0.5, 5.5);
	139	fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
	140	fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
	141	fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
	142	fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
	143	fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
	144
	145	for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
	146	fHistPtLeadingJet[iBranch] = new TH1F(Form("fHistPtLeadingJet%i", iBranch),
	147	Form("p_{T} of leading jet from branch %i;p_{T} (GeV/c);counts", iBranch),
	148	300, 0., 300.);
	149	fHistPtLeadingJet[iBranch]->SetMarkerStyle(kFullCircle);
	150
	151	fHistEtaPhiLeadingJet[iBranch] = new TH2F(Form("fHistEtaPhiLeadingJet%i", iBranch),
	152	Form("#eta - #phi of leading jet from branch %i;#eta;#phi", iBranch),
	153	100, -2., 2., 100, 0., 2*TMath::Pi());
	154
	155	}
	156	fHistDeltaEtaDeltaPhiLeadingJet = new TH2F("fHistDeltaEtaDeltaPhiLeadingJet",
	157	"#Delta#eta - #Delta#phi of leading jet;#Delta#eta;#Delta#phi",
	158	100, -4., 4., 100, -2.TMath::Pi(), 2TMath::Pi());
	159
	160	fHistPtPtExtra = new TH2F("fHistPtPtExtra", "p_{T} response;p_{T} (GeV/c);p_{T} (GeV/c)",
	161	300, 0., 300., 300, 0., 300.);
	162
163	for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
164	fHistPtResponse[iEvtClass] = new TH2F(Form("pt_response%i",iEvtClass), Form("pt_response%i;p_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
165	300, 0., 300., 300, 0., 300.);
166	fHistPtSmearing[iEvtClass] = new TH2F(Form("pt_smearing%i",iEvtClass), Form("pt_smearing%i;p_{T} (GeV/c);p_{T} (GeV/c)",iEvtClass),
167	200, -50., 150., 300, 0., 300.);
a11972e9	168	fHistdR[iEvtClass] = new TH2F(Form("hist_dR%i",iEvtClass), "", 200, -50.,150., 200, -1.,1.);
	169	fHistArea[iEvtClass] = new TH2F(Form("hist_Area%i",iEvtClass), "", 200, -50.,150., 100, 0.,1.);
	170
	171	fHistEtaPhiLeadingJetCut[iEvtClass] = new TH2F(Form("fHistEtaPhiLeadingJetCut%i", iEvtClass),
	172	Form("#eta - #phi of leading jet from event class %i;#eta;#phi", iEvtClass),
	173	100, -2., 2., 100, 0., 2*TMath::Pi());
	174	fHistDeltaEtaEtaLeadingJet[iEvtClass] = new TH2F(Form("fHistDeltaEtaEtaLeadingJet%i", iEvtClass),
	175	Form("#eta - #Delta#eta of leading jet from event class %i;#eta;#Delta#eta", iEvtClass),
	176	100, -1., 1., 100, -.5, .5);
	177	fHistDeltaPtEtaLeadingJet[iEvtClass] = new TH2F(Form("fHistDeltaPtEtaLeadingJet%i", iEvtClass),
	178	Form("#eta - #Delta#eta of leading jet from event class %i;#eta;#Delta#p_{T}", iEvtClass),
	179	100, -1., 1., 200, -50., 150);
2a436fa1	180	}
	181
	182	OpenFile(1);
	183	fOutputList = new TList;
	184	fOutputList->Add(fHistEvtSelection);
	185	for (Int_t iBranch = 0; iBranch < fkNbranches; iBranch++) {
	186	fOutputList->Add(fHistPtLeadingJet[iBranch]);
	187	fOutputList->Add(fHistEtaPhiLeadingJet[iBranch]);
	188	}
	189	fOutputList->Add(fHistDeltaEtaDeltaPhiLeadingJet);
	190	fOutputList->Add(fHistPtPtExtra);
	191	for (Int_t iEvtClass =0; iEvtClass < fkEvtClasses; iEvtClass++) {
	192	fOutputList->Add(fHistPtResponse[iEvtClass]);
	193	fOutputList->Add(fHistPtSmearing[iEvtClass]);
a11972e9	194	fOutputList->Add(fHistdR[iEvtClass]);
	195	fOutputList->Add(fHistArea[iEvtClass]);
	196	fOutputList->Add(fHistEtaPhiLeadingJetCut[iEvtClass]);
	197	fOutputList->Add(fHistDeltaEtaEtaLeadingJet[iEvtClass]);
	198	fOutputList->Add(fHistDeltaPtEtaLeadingJet[iEvtClass]);
2a436fa1	199	}
	200	}
	201
	202	void AliAnalysisTaskJetResponse::UserExec(Option_t *)
	203	{
	204	// load events, apply event cuts, then compare leading jets
	205
	206	if(!strlen(fJetBranchName[0].Data()) \|\| !strlen(fJetBranchName[1].Data())){
	207	AliError("Jet branch name not set.");
	208	return;
	209	}
	210
	211	fESD=dynamic_cast<AliESDEvent*>(InputEvent());
	212	if (!fESD) {
	213	AliError("ESD not available");
	214	return;
	215	}
	216	fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
	217	if (!fAOD) {
	218	AliError("AOD not available");
	219	return;
	220	}
	221
	222	fHistEvtSelection->Fill(1); // number of events before event selection
	223
	224	// physics selection
	225	AliInputEventHandler* inputHandler = (AliInputEventHandler*)
	226	((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
	227	if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
	228	Printf(" Trigger Selection: event REJECTED ... ");
	229	fHistEvtSelection->Fill(2);
	230	PostData(1, fOutputList);
	231	return;
	232	}
	233
	234	// vertex selection
	235	AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
	236	Int_t nTracksPrim = primVtx->GetNContributors();
	237	if ((nTracksPrim < fMinContribVtx) \|\|
	238	(primVtx->GetZ() < fVtxZMin) \|\|
	239	(primVtx->GetZ() > fVtxZMax) ){
	240	Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
	241	fHistEvtSelection->Fill(3);
	242	PostData(1, fOutputList);
	243	return;
	244	}
	245
	246	// event class selection (from jet helper task)
	247	Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
	248	Printf("Event class %d", eventClass);
	249	if (eventClass < fEvtClassMin \|\| eventClass > fEvtClassMax){
	250	fHistEvtSelection->Fill(4);
	251	PostData(1, fOutputList);
	252	return;
	253	}
	254
	255	// centrality selection
	256	AliCentrality *cent = fESD->GetCentrality();
	257	Float_t centValue = cent->GetCentralityPercentile("TRK");
	258	printf("centrality: %f\n", centValue);
	259	if (centValue < fCentMin \|\| centValue > fCentMax){
	260	fHistEvtSelection->Fill(5);
	261	PostData(1, fOutputList);
	262	return;
263	}
264
265	fHistEvtSelection->Fill(0); // accepted events
266
267
268	TClonesArray *aodJets[2];
a11972e9	269	aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
a11972e9	270	aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE
2a436fa1	271	AliAODJet *leadingJet[2] = { 0x0, 0x0 };
	272
	273	for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
	274	fListJets[iJetType]->Clear();
	275	if (!aodJets[iJetType])
	276	continue;
	277	for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
	278	AliAODJet jet = dynamic_cast<AliAODJet>((*aodJets[iJetType])[iJet]);
	279	if (jet)
	280	fListJets[iJetType]->Add(jet);
	281	}
	282	fListJets[iJetType]->Sort();
	283	leadingJet[iJetType] = (AliAODJet*) fListJets[iJetType]->First();
	284
	285	if (leadingJet[iJetType])
	286	fHistEtaPhiLeadingJet[iJetType]->Fill(leadingJet[iJetType]->Eta(), leadingJet[iJetType]->Phi());
	287	}
	288
	289	// check if leading jets are close in eta-phi and compare their Pt
	290	if (leadingJet[0] && leadingJet[1]) {
a11972e9	291	Float_t deltaEta = leadingJet[0]->Eta() - leadingJet[1]->Eta();
	292	Float_t deltaPhi = TVector2::Phi_mpi_pi(leadingJet[0]->Phi() - leadingJet[1]->Phi());
	293
	294	if (eventClass > -1 && eventClass < fkEvtClasses){
	295
	296	if(leadingJet[1]->Eta()<fJetEtaMax && leadingJet[1]->Eta()>fJetEtaMin){
	297	fHistEtaPhiLeadingJetCut[eventClass]->Fill(leadingJet[1]->Eta(), leadingJet[1]->Phi());
	298	}
	299	}
2a436fa1	300
	301	// leading jets in "jet" acceptance
	302	if(leadingJet[0]->Eta()>fJetEtaMax \|\| leadingJet[0]->Eta()<fJetEtaMin \|\|
	303	leadingJet[1]->Eta()>fJetEtaMax \|\| leadingJet[1]->Eta()<fJetEtaMin){
	304	Printf("Jet not in eta acceptance.");
	305	}
	306	else{
	307	// check association of jets
2a436fa1	308	fHistDeltaEtaDeltaPhiLeadingJet->Fill(deltaEta, deltaPhi);
	309
	310	if (TMath::Abs(deltaEta) > fJetDeltaEta \|\| (TMath::Pi() - TMath::Abs(TMath::Abs(deltaPhi) - TMath::Pi())) > fJetDeltaPhi)
	311	printf("leading jets two far apart\n");
	312	else {
a11972e9	313	Float_t pt0 = leadingJet[0]->Pt();
	314	Float_t pt1 = leadingJet[1]->Pt();
	315	Float_t dPt = pt1-pt0;
	316	Float_t jetarea = leadingJet[1]->EffectiveAreaCharged();
	317
	318	fHistPtLeadingJet[0]->Fill(pt0);
	319	fHistPtLeadingJet[1]->Fill(pt1);
2a436fa1	320
a11972e9	321	fHistPtPtExtra->Fill(pt0, pt1);
2a436fa1	322
2a436fa1	323	if (eventClass > -1 && eventClass < fkEvtClasses){
a11972e9	324	fHistPtResponse[eventClass]->Fill(pt1, pt0);
	325	fHistPtSmearing[eventClass]->Fill(dPt, pt0);
	326
	327	Float_t dR = TMath::Sqrt(deltaEtadeltaEta+deltaPhideltaPhi);
	328	fHistdR[eventClass]->Fill(dPt, dR);
	329	fHistArea[eventClass]->Fill(dPt, jetarea);
	330
	331	fHistDeltaEtaEtaLeadingJet[eventClass]->Fill(leadingJet[0]->Eta(), deltaEta);
	332	fHistDeltaPtEtaLeadingJet[eventClass]->Fill(leadingJet[0]->Eta(), dPt);
2a436fa1	333	}
	334	}
	335	}
	336	}
	337
	338	PostData(1, fOutputList);
	339	}
	340
	341	void AliAnalysisTaskJetResponse::Terminate(const Option_t *)
	342	{
	343	// Draw result to the screen
	344	// Called once at the end of the query
	345
	346	if (!GetOutputData(1))
	347	return;
	348	}