PWGJE/AliAnalysisTaskJetResponseV2.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //
  17 // task compares jets in two or three branches,
  18 // written for analysis of jet embedding in HI events
  19 //
  20
  21
  22 #include "TChain.h"
  23 #include "TTree.h"
  24 #include "TMath.h"
  25 #include "TH1F.h"
  26 #include "TH1I.h"
  27 #include "TH2F.h"
  28 #include "TH3F.h"
  29 #include "THnSparse.h"
  30 #include "TCanvas.h"
  31
  32 #include "AliLog.h"
  33
  34 #include "AliAnalysisTask.h"
  35 #include "AliAnalysisManager.h"
  36
  37 #include "AliVEvent.h"
  38 #include "AliESDEvent.h"
  39 #include "AliESDInputHandler.h"
  40 #include "AliCentrality.h"
  41 #include "AliAnalysisHelperJetTasks.h"
  42 #include "AliInputEventHandler.h"
  43 #include "AliAODJetEventBackground.h"
  44 #include "AliAnalysisTaskFastEmbedding.h"
  45
  46 #include "AliAODEvent.h"
  47 #include "AliAODJet.h"
  48 #include "AliAODHandler.h"
  49
  50 #include "AliAnalysisTaskJetResponseV2.h"
  51
  52 using std::cout;
  53 using std::endl;
  54
  55 ClassImp(AliAnalysisTaskJetResponseV2)
  56
  57 AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2() :
  58 AliAnalysisTaskSE(),
  59   fESD(0x0),
  60   fAOD(0x0),
  61   fAODOut(0x0),
  62   fAODExtension(0x0),
  63   fNonStdFile(""),
  64   fBackgroundBranch(""),
  65   fIsPbPb(kTRUE),
  66   fOfflineTrgMask(AliVEvent::kAny),
  67   fMinContribVtx(1),
  68   fVtxZMin(-8.),
  69   fVtxZMax(8.),
  70   fEvtClassMin(0),
  71   fEvtClassMax(4),
  72   fCentMin(0.),
  73   fCentMax(100.),
  74   fNInputTracksMin(0),
  75   fNInputTracksMax(-1),
  76   fJetEtaMin(-.5),
  77   fJetEtaMax(.5),
  78   fJetPtMin(20.),
  79   fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
  80   fJetPtFractionMin(0.5),
  81   fNMatchJets(4),
  82   fMatchMaxDist(0.8),
  83   fKeepJets(kFALSE),
  84   fkNbranches(2),
  85   fkEvtClasses(12),
  86   fOutputList(0x0),
  87   fbEvent(kTRUE),
  88   fbJetsMismatch1(kTRUE),
  89   fbJetsMismatch2(kTRUE),
  90   fbJetsRp(kTRUE),
  91   fbJetsDeltaPt(kTRUE),
  92   fbJetsEta(kTRUE),
  93   fbJetsPhi(kTRUE),
  94   fbJetsArea(kTRUE),
  95   fbJets3Branches(kFALSE),
  96   fbJetsBeforeCut1(kTRUE),
  97   fbJetsBeforeCut2(kTRUE),
  98   fHistEvtSelection(0x0),
  99   fHistJetSelection(0x0),
 100   fh2JetSelection(0x0),
 101   fhnEvent(0x0),
 102   fhnJetsMismatch1(0x0),
 103   fhnJetsMismatch2(0x0),
 104   fhnJetsRp(0x0),
 105   fhnJetsDeltaPt(0x0),
 106   fhnJetsEta(0x0),
 107   fhnJetsPhi(0x0),
 108   fhnJetsArea(0x0),
 109   fhnJets3Branches(0x0),
 110   fhnJetsBeforeCut1(0x0),
 111   fhnJetsBeforeCut2(0x0)
 112 {
 113   // default Constructor
 114
 115   fJetBranchName[0] = "";
 116   fJetBranchName[1] = "";
 117   fJetBranchName[2] = "";
 118
 119   fListJets[0] = new TList;
 120   fListJets[1] = new TList;
 121   fListJets[2] = new TList;
 122 }
 123
 124 AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2(const char *name) :
 125   AliAnalysisTaskSE(name),
 126   fESD(0x0),
 127   fAOD(0x0),
 128   fAODOut(0x0),
 129   fAODExtension(0x0),
 130   fNonStdFile(""),
 131   fBackgroundBranch(""),
 132   fIsPbPb(kTRUE),
 133   fOfflineTrgMask(AliVEvent::kAny),
 134   fMinContribVtx(1),
 135   fVtxZMin(-8.),
 136   fVtxZMax(8.),
 137   fEvtClassMin(0),
 138   fEvtClassMax(4),
 139   fCentMin(0.),
 140   fCentMax(100.),
 141   fNInputTracksMin(0),
 142   fNInputTracksMax(-1),
 143   fJetEtaMin(-.5),
 144   fJetEtaMax(.5),
 145   fJetPtMin(20.),
 146   fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
 147   fJetPtFractionMin(0.5),
 148   fNMatchJets(4),
 149   fMatchMaxDist(0.8),
 150   fKeepJets(kFALSE),
 151   fkNbranches(2),
 152   fkEvtClasses(12),
 153   fOutputList(0x0),
 154   fbEvent(kTRUE),
 155   fbJetsMismatch1(kTRUE),
 156   fbJetsMismatch2(kTRUE),
 157   fbJetsRp(kTRUE),
 158   fbJetsDeltaPt(kTRUE),
 159   fbJetsEta(kTRUE),
 160   fbJetsPhi(kTRUE),
 161   fbJetsArea(kTRUE),
 162   fbJets3Branches(kFALSE),
 163   fbJetsBeforeCut1(kTRUE),
 164   fbJetsBeforeCut2(kTRUE),
 165   fHistEvtSelection(0x0),
 166   fHistJetSelection(0x0),
 167   fh2JetSelection(0x0),
 168   fhnEvent(0x0),
 169   fhnJetsMismatch1(0x0),
 170   fhnJetsMismatch2(0x0),
 171   fhnJetsRp(0x0),
 172   fhnJetsDeltaPt(0x0),
 173   fhnJetsEta(0x0),
 174   fhnJetsPhi(0x0),
 175   fhnJetsArea(0x0),
 176   fhnJets3Branches(0x0),
 177   fhnJetsBeforeCut1(0x0),
 178   fhnJetsBeforeCut2(0x0)
 179 {
 180   // Constructor
 181
 182   fJetBranchName[0] = "";
 183   fJetBranchName[1] = "";
 184   fJetBranchName[2] = "";
 185
 186   fListJets[0] = new TList;
 187   fListJets[1] = new TList;
 188   fListJets[2] = new TList;
 189
 190   DefineOutput(1, TList::Class());
 191 }
 192
 193 AliAnalysisTaskJetResponseV2::~AliAnalysisTaskJetResponseV2()
 194 {
 195   delete fListJets[0];
 196   delete fListJets[1];
 197   delete fListJets[2];
 198 }
 199
 200 void AliAnalysisTaskJetResponseV2::SetBranchNames(const TString &branch1, const TString &branch2, const TString &branch3)
 201 {
 202   fJetBranchName[0] = branch1;
 203   fJetBranchName[1] = branch2;
 204   fJetBranchName[2] = branch3;
 205
 206   if(strlen(fJetBranchName[2].Data()) ) {
 207     fkNbranches = 3;
 208     fbJets3Branches = kTRUE;
 209   }
 210 }
 211
 212 void AliAnalysisTaskJetResponseV2::Init()
 213 {
 214
 215   // check for jet branches
 216   if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 217     AliError("Jet branch name not set.");
 218   }
 219
 220 }
 221
 222 void AliAnalysisTaskJetResponseV2::UserCreateOutputObjects()
 223 {
 224   // Create histograms
 225   // Called once
 226   OpenFile(1);
 227   if(!fOutputList) fOutputList = new TList;
 228   fOutputList->SetOwner(kTRUE);
 229
 230   Bool_t oldStatus = TH1::AddDirectoryStatus();
 231   TH1::AddDirectory(kFALSE);
 232
 233
 234   fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
 235   fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
 236   fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
 237   fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
 238   fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
 239   fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
 240   fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
 241
 242   fHistJetSelection = new TH1I("fHistJetSelection", "jet selection", 8, -0.5, 7.5);
 243   fHistJetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
 244   fHistJetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
 245   fHistJetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
 246   fHistJetSelection->GetXaxis()->SetBinLabel(4,"not in list");
 247   fHistJetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
 248   fHistJetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
 249   fHistJetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
 250   fHistJetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
 251
 252   fh2JetSelection = new TH2F("fh2JetSelection", "jet selection", 8, -0.5, 7.5,100,0.,200.);
 253   fh2JetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
 254   fh2JetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
 255   fh2JetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
 256   fh2JetSelection->GetXaxis()->SetBinLabel(4,"not in list");
 257   fh2JetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
 258   fh2JetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
 259   fh2JetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
 260   fh2JetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
 261
 262
 263   UInt_t entries = 0; // bit coded, see GetDimParams() below
 264   UInt_t opt = 0;  // bit coded, default (0) or high resolution (1)
 265
 266   if(fbEvent){
 267     entries = 1<<0 | 1<<1 | 1<<2 | 1<<26;  // cent : nInpTrks : rp psi : pT hard bin
 268     opt = 1<<0 | 1<<1; // centrality and nInpTrks in high resolution
 269     fhnEvent = NewTHnSparseF("fhnEvent", entries, opt);
 270   }
 271
 272   if(fbJetsMismatch1){ // real mismatch (no related rec jet found)
 273     // cent : nInpTrks : rp bins : probe pt : probe eta : probe phi : probe area : pT hard bin
 274     entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<8 | 1<<10 | 1<<12 | 1<<26;
 275     opt =  1<<6 | 1<<8 | 1<<10;
 276     fhnJetsMismatch1 = NewTHnSparseF("fhnJetsMismatch1", entries, opt);
 277   }
 278
 279   if(fbJetsMismatch2){  // acceptance + fraction cut
 280     // cent : nInpTrks : rp bins : jetPt(2x) : jetEta(2x) : deltaEta : deltaR : fraction : pT hard bin
 281     entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<15 | 1<<17 | 1<<19 | 1<<26;
 282     opt = 1<<6 | 1<<7 | 1<<8 | 1<<9;
 283     fhnJetsMismatch2 = NewTHnSparseF("fhnJetsMismatch2", entries, opt);
 284
 285   }
 286
 287
 288   if(fbJetsRp){
 289     // cent : nInpTrks : rp bins : rp wrt jet(2x) : probe pT : probe area: deltaPt : rp phi : rho : correction for RP : local rho
 290     /*
 291       entries = 1<<0 | 1<<1 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<12 | 1<<14 | 1<<2 | 1<<22 | 1<<23 | 1<<24 | 1<<25;
 292       opt = 1<<4 | 1<<5;*/
 293     // cent : nInpTrks : rp bins : rp wrt jet(2x) : probe pT : deltaPt : rp phi : rho : correction for RP | pT hard bin
 294     entries = 1<<0 | 1<<1 | 1<<3 | 1<<4 | 1<<5 | 1<<6 | 1<<14 | 1<<2 | 1<<22 | 1<<23 | 1<<26;
 295     opt = 1<<4 | 1<<5;
 296     fhnJetsRp = NewTHnSparseF("fhnJetsRp", entries, opt);
 297   }
 298
 299   // cent : nInpTrks : rp bins:  deltaPt : jetPt(2x) : deltaArea : pT hard bin (hr delta pt)
 300   if(fbJetsDeltaPt){
 301     entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<18 | 1<<26;
 302     opt = 1<<1 | 1<<14 | 1<<6 | 1<<7 | 1<<18;
 303     fhnJetsDeltaPt = NewTHnSparseF("fhnJetsDeltaPt", entries, opt);
 304   }
 305
 306   // cent : nInpTrks : rp bins : deltaPt : jetPt(2x) : deltaR : deltaEta : jetEta(2x) : pT hard bin (hr for eta)
 307   if(fbJetsEta){
 308     entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<17 | 1<<15 | 1<<8 | 1<<9 | 1<<26;
 309     opt = 1<<15 | 1<<8 | 1<<9;
 310     fhnJetsEta = NewTHnSparseF("fhnJetsEta", entries, opt);
 311   }
 312
 313   // cent : nInpTrks : rp bins : jetPt(2x) : jetPhi(2x) : deltaPt : deltaPhi : pT hard bin
 314   if(fbJetsPhi){
 315     entries = 1<<0 | 1<<1 | 1<<3 | 1<<6 | 1<<7 | 1<<10 | 1<<11 | 1<<14 | 1<<16 | 1<<26;
 316     opt = 1<<10 | 1<<11;
 317     fhnJetsPhi = NewTHnSparseF("fhnJetsPhi", entries, opt);
 318   }
 319
 320   // cent : nInpTrks : rp bins : deltaArea : jetArea(2x) : deltaR : fraction : distance next rec jet : pT next jet : deltaPt : jetPt(2x) : pT hard bin (hr for area)
 321   if(fbJetsArea){
 322     entries = 1<<0 | 1<<1 | 1<<3 | 1<<18 | 1<<12 | 1<<13 | 1<<17 | 1<<19 | 1<<20 | 1<<21 | 1<<14 | 1<<6 | 1<<7 | 1<<26;
 323     opt = 1<<18 | 1<<12 | 1<<13;
 324     fhnJetsArea = NewTHnSparseF("fhnJetsArea", entries, opt);
 325   }
 326
 327   // cent : nInpTrks : jetPt(3x) : deltaPt : delta : jetArea(3x) : fraction(2x) : deltaR(1x) : pT hard bin
 328   if(fbJets3Branches){
 329     entries = 1<<0 | 1<<1 | 1<<6 | 1<<7 | 1<<27 | 1<<14 | 1<<28 | 1<<12 | 1<<13 | 1<<29 | 1<<19 | 1<<30 | 1<<17 | 1<<26;
 330     opt = 1<<6 | 1<<7 | 1<<27 | 1<<14 | 1<<28;
 331     fhnJets3Branches = NewTHnSparseF("fhnJets3Branches", entries, opt);
 332   }
 333
 334
 335
 336
 337   //before cut
 338
 339   // cent : nInpTrks : rp bins : fraction : jetPt(2x) : jetEta(2x) : jetPhi(2x) (low resolution) (with fraction, eta, phi, pt cuts possible)
 340   if(fbJetsBeforeCut1){
 341     entries = 1<<0 | 1<<1 | 1<<3 | 1<<19 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<10 | 1<<11;
 342     opt = 0;
 343     fhnJetsBeforeCut1 = NewTHnSparseF("fhnJetsBeforeCut1", entries, opt);
 344   }
 345
 346   // cent : nInpTrks : rp bins : deltaPt : jetPt(2x) : deltaR : deltaEta : jetEta(2x) (low resolution)
 347   if(fbJetsBeforeCut2){
 348     entries = 1<<0 | 1<<1 | 1<<3 | 1<<14 | 1<<6 | 1<<7 | 1<<17 | 1<<15 | 1<<8 | 1<<9;
 349     opt = 0;
 350     fhnJetsBeforeCut2 = NewTHnSparseF("fhnJetsBeforeCut2", entries, opt);
 351   }
 352
 353   fOutputList->Add(fHistEvtSelection);
 354   fOutputList->Add(fHistJetSelection);
 355   fOutputList->Add(fh2JetSelection);
 356   if(fbEvent)           fOutputList->Add(fhnEvent);
 357   if(fbJetsMismatch1)   fOutputList->Add(fhnJetsMismatch1);
 358   if(fbJetsMismatch2)   fOutputList->Add(fhnJetsMismatch2);
 359   if(fbJetsRp)          fOutputList->Add(fhnJetsRp);
 360   if(fbJetsDeltaPt)     fOutputList->Add(fhnJetsDeltaPt);
 361   if(fbJetsEta)         fOutputList->Add(fhnJetsEta);
 362   if(fbJetsPhi)         fOutputList->Add(fhnJetsPhi);
 363   if(fbJetsArea)        fOutputList->Add(fhnJetsArea);
 364   if(fbJets3Branches)   fOutputList->Add(fhnJets3Branches);
 365   if(fbJetsBeforeCut1)  fOutputList->Add(fhnJetsBeforeCut1);
 366   if(fbJetsBeforeCut2)  fOutputList->Add(fhnJetsBeforeCut2);
 367
 368   // =========== Switch on Sumw2 for all histos ===========
 369   for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
 370     TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
 371     if (h1){
 372       h1->Sumw2();
 373       continue;
 374     }
 375     THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
 376     if (hn){
 377       hn->Sumw2();
 378     }
 379   }
 380   TH1::AddDirectory(oldStatus);
 381
 382   PostData(1, fOutputList);
 383 }
 384
 385 void AliAnalysisTaskJetResponseV2::UserExec(Option_t *)
 386 {
 387   // load events, apply event cuts, then compare leading jets
 388
 389   if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 390     AliError("Jet branch name not set.");
 391     return;
 392   }
 393
 394   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
 395   if (!fESD) {
 396     AliError("ESD not available");
 397
 398     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 399     //  assume that the AOD is in the general output...
 400     fAODOut  = AODEvent();
 401   } else {
 402     fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
 403   }
 404   if (!fAOD) {
 405     AliError("AOD not available");
 406     return;
 407   }
 408
 409   if(fNonStdFile.Length()!=0){
 410     // case that we have an AOD extension we need can fetch the jets from the extended output
 411     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 412     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 413     if(!fAODExtension){
 414       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
 415     }
 416   }
 417
 418   // -- event selection --
 419   fHistEvtSelection->Fill(1); // number of events before event selection
 420
 421   // physics selection
 422   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
 423     ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 424   if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
 425     if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
 426     fHistEvtSelection->Fill(2);
 427     PostData(1, fOutputList);
 428     return;
 429   }
 430
 431   // vertex selection
 432   AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
 433   Int_t nTracksPrim = primVtx->GetNContributors();
 434   if ((nTracksPrim < fMinContribVtx) ||
 435       (primVtx->GetZ() < fVtxZMin) ||
 436       (primVtx->GetZ() > fVtxZMax) ){
 437     if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
 438     fHistEvtSelection->Fill(3);
 439     PostData(1, fOutputList);
 440     return;
 441   }
 442
 443   // event class selection (from jet helper task)
 444   Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
 445   if(fDebug) Printf("Event class %d", eventClass);
 446   if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
 447     fHistEvtSelection->Fill(4);
 448     PostData(1, fOutputList);
 449     return;
 450   }
 451
 452   // centrality selection
 453   AliCentrality *cent = 0x0;
 454   Float_t centValue = 0.;
 455   if(fESD) cent = fESD->GetCentrality();
 456   if(cent) centValue = cent->GetCentralityPercentile("V0M");
 457   if(fDebug) printf("centrality: %f\n", centValue);
 458   if (centValue < fCentMin || centValue > fCentMax){
 459     fHistEvtSelection->Fill(4);
 460     PostData(1, fOutputList);
 461     return;
 462   }
 463
 464
 465   // multiplicity due to input tracks
 466   Int_t nInputTracks = GetNInputTracks();
 467
 468   if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
 469     fHistEvtSelection->Fill(5);
 470     PostData(1, fOutputList);
 471     return;
 472   }
 473
 474
 475   fHistEvtSelection->Fill(0); // accepted events
 476   // -- end event selection --
 477
 478   // pt hard
 479   Double_t pthard = AliAnalysisTaskFastEmbedding::GetPtHard();
 480   Int_t pthardbin = GetPtHardBin(pthard);
 481
 482   // reaction plane
 483   Double_t rp = AliAnalysisHelperJetTasks::ReactionPlane(kFALSE);
 484   if(fbEvent){
 485     Double_t eventEntries[3] = { (Double_t)centValue, (Double_t)nInputTracks, rp };
 486     fhnEvent->Fill(eventEntries);
 487   }
 488
 489
 490   // get background
 491   AliAODJetEventBackground* externalBackground = 0;
 492   if(!externalBackground&&fBackgroundBranch.Length()){
 493     externalBackground =  (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
 494     if(!externalBackground && fAODOut) externalBackground =  (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
 495     if(!externalBackground && fAODExtension)  externalBackground =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
 496     //if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 497   }
 498   Float_t rho = 0;
 499   if(externalBackground) rho = externalBackground->GetBackground(0);
 500
 501
 502   // fetch jets
 503   TClonesArray *aodJets[3];
 504   aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
 505   if(!aodJets[0] && fAODOut) aodJets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
 506   if(!aodJets[0] && fAODExtension) aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
 507   aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE version1
 508   if(!aodJets[1] && fAODOut) aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet
 509   if(!aodJets[1] && fAODExtension) aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));
 510   if( strlen(fJetBranchName[2].Data()) ) {
 511     aodJets[2] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[2].Data())); // in general: embedded jet + UE version2
 512     if(!aodJets[2] && fAODOut) aodJets[2] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[2].Data())); // in general: embedded jet
 513     if(!aodJets[2] && fAODExtension) aodJets[2] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[2].Data()));
 514     if(aodJets[2]) fkNbranches=3;
 515     if(fDebug>10) printf("3rd branch: %s\n",fJetBranchName[2].Data());
 516   }
 517
 518   if(fDebug>10)  printf("fkNbranches %d\n",fkNbranches);
 519   for (Int_t iJetType = 0; iJetType < fkNbranches; iJetType++) {
 520     fListJets[iJetType]->Clear();
 521     if (!aodJets[iJetType]) {
 522         if(fDebug) Printf("%s: no jet branch\n",fJetBranchName[iJetType].Data());
 523       continue;
 524     }
 525     if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
 526
 527     for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
 528       AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
 529       if (jet) fListJets[iJetType]->Add(jet);
 530     }
 531     fListJets[iJetType]->Sort();
 532   }
 533
 534
 535   // jet matching
 536   static TArrayI aMatchIndex(fListJets[0]->GetEntries());
 537   static TArrayF aPtFraction(fListJets[0]->GetEntries());
 538   if(aMatchIndex.GetSize()<fListJets[0]->GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries());
 539   if(aPtFraction.GetSize()<fListJets[0]->GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries());
 540
 541   // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction'
 542   AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
 543                                             fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()),
 544                                             aMatchIndex, aPtFraction, fDebug, fMatchMaxDist, fIsPbPb?1:2);
 545   static TArrayI aMatchIndexv2(fListJets[0]->GetEntries());
 546   static TArrayF aPtFractionv2(fListJets[0]->GetEntries());
 547   if(aMatchIndexv2.GetSize()<fListJets[0]->GetEntries()) aMatchIndexv2.Set(fListJets[0]->GetEntries());
 548   if(aPtFractionv2.GetSize()<fListJets[0]->GetEntries()) aPtFractionv2.Set(fListJets[0]->GetEntries());
 549   if( strlen(fJetBranchName[2].Data()) ) {
 550     AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
 551                                               fListJets[2], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[2]->GetEntries()),
 552                                               aMatchIndexv2, aPtFractionv2, fDebug, fMatchMaxDist, fIsPbPb?1:2);
 553   }
 554
 555   // loop over matched jets
 556   Int_t ir = -1;  // index of matched reconstruced jet
 557   Int_t ir2 = -1; // index of matched reconstruced jet
 558   Float_t fraction = -1.;
 559   Float_t fraction2 = -1.;
 560   AliAODJet *jet[3]  = { 0x0, 0x0, 0x0 };
 561   Float_t jetEta[3]  = { -990., -990., -990. };
 562   Float_t jetPhi[3]  = { -990., -990., -990. };
 563   Float_t jetPt[3]   = { -990., -990., -990. };
 564   Float_t jetArea[3] = { -990., -990., -990. };
 565   Float_t rpJet[3]   = { -990., -990., -990. };
 566   Int_t rpBin = -990;
 567
 568   for(Int_t ig=0; ig<fListJets[0]->GetEntries(); ++ig){
 569     ir = aMatchIndex[ig];
 570     if(aMatchIndexv2.GetSize()>ig) ir2 = aMatchIndexv2[ig];
 571     else ir2 =0;
 572
 573     //fetch jets
 574     jet[0] = (AliAODJet*)(fListJets[0]->At(ig));
 575     if(ir>=0) jet[1] = (AliAODJet*)(fListJets[1]->At(ir));
 576     else      jet[1] = 0x0;
 577     if(ir2>=0) jet[2] = (AliAODJet*)(fListJets[2]->At(ir2));
 578     else      jet[2] = 0x0;
 579
 580     for(Int_t i=0; i<fkNbranches; ++i){
 581       if(!jet[i]){
 582         jetEta[i]  = -990;
 583         jetPhi[i]  = -990.;
 584         jetPt[i]   = -990.;
 585         jetArea[i] = -990.;
 586         rpJet[i]   = -990.;
 587         if(i==1) rpBin = -990;
 588       }
 589       else {
 590         jetEta[i]  = jet[i]->Eta();
 591         jetPhi[i]  = jet[i]->Phi();
 592         jetPt[i]   = GetPt(jet[i], i);
 593         jetArea[i] = jet[i]->EffectiveAreaCharged();
 594         rpJet[i]   = TVector2::Phi_mpi_pi(rp-jetPhi[i]);
 595         if(i==1) rpBin = AliAnalysisHelperJetTasks::GetPhiBin(TVector2::Phi_mpi_pi(rp-jetPhi[i]), 3);
 596       }
 597     }
 598     fraction  = aPtFraction[ig];
 599     if(aPtFractionv2.GetSize()>ig) fraction2 = aPtFractionv2[ig];
 600     else fraction2 = 0.;
 601
 602     // jet statistics
 603     fHistJetSelection->Fill(1); // all probe jets
 604     if(jet[0]) fh2JetSelection->Fill(1.,jetPt[0]); // all probe jets
 605
 606     if(ir<0){
 607       fHistJetSelection->Fill(2);
 608       if(jet[0]) fh2JetSelection->Fill(2.,jetPt[0]);
 609
 610       if(fbJetsMismatch1){
 611         if(!jet[0]) continue;
 612         Double_t jetEntriesMismatch1[7] = {
 613           (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 614           (Double_t)jetPt[0], (Double_t)jetEta[0], (Double_t)jetPhi[0], (Double_t)jetArea[0]
 615         };
 616         fhnJetsMismatch1->Fill(jetEntriesMismatch1);
 617       }
 618
 619       continue;
 620     }
 621
 622     if(!jet[0] || !jet[1]){
 623       fHistJetSelection->Fill(3);
 624       if(jet[0]) fh2JetSelection->Fill(3.,jetPt[0]);
 625       continue;
 626     }
 627
 628     // look for distance to next rec jet
 629     Float_t distNextJet = -0.01; // no neighbor
 630     Float_t ptNextJet = -1.;
 631     for(Int_t r=0; r<fListJets[1]->GetEntries(); ++r){
 632       if(r==ir) continue;
 633       Float_t tmpDeltaR = jet[1]->DeltaR((AliAODJet*)fListJets[1]->At(r));
 634       if(distNextJet<0. || distNextJet>tmpDeltaR){
 635         distNextJet = tmpDeltaR;
 636         if(fKeepJets) ptNextJet   = ((AliAODJet*)fListJets[1]->At(r))->GetPtSubtracted(0);
 637         else          ptNextJet   = ((AliAODJet*)fListJets[1]->At(r))->Pt();
 638       }
 639     }
 640
 641
 642
 643     //Float_t localRho = jetArea[1]>0. ? (jetPt[1]+rho*jetArea[1] - jetPt[0]) / jetArea[1] : 0.;
 644     //Float_t relRho = rho>0. ? localRho / rho : 0.;
 645
 646
 647     // calculate parameters of associated jets
 648     /* from Leticia, kT clusterizer
 649        Float_t par0[4] = { 0.00409,       0.01229,      0.05,         0.26 };
 650        Float_t par1[4] = { -2.97035e-03, -2.03182e-03, -1.25702e-03, -9.95107e-04 };
 651        Float_t par2[4] = { 1.02865e-01,   1.49039e-01,  1.53910e-01,  1.51109e-01 };
 652     */
 653     // own, from embedded tracks
 654     Float_t par0[4] = { 0.02841,       0.05039,      0.09092,      0.24089     };
 655     Float_t par1[4] = { -4.26725e-04, -1.15273e-03, -1.56827e-03, -3.08003e-03 };
 656     Float_t par2[4] = { 4.95415e-02,   9.79538e-02,  1.32814e-01,  1.71743e-01 };
 657
 658     Float_t rpCorr = 0.;
 659
 660     if(eventClass>0&&eventClass<4){
 661       rpCorr = par0[eventClass-1] + par1[eventClass-1] * 2*TMath::Cos(rpJet[1]) + par2[eventClass-1] * 2*TMath::Cos(2*rpJet[1]);
 662       rpCorr *= rho * jetArea[1];
 663     }
 664
 665     Float_t delta      = jetPt[2]-jetPt[1];
 666     Float_t deltaPt    = jetPt[1]-jetPt[0];
 667     Float_t deltaEta   = jetEta[1]-jetEta[0];
 668     Float_t deltaPhi   = TVector2::Phi_mpi_pi(jetPhi[1]-jetPhi[0]);
 669     Float_t deltaR     = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi);
 670     Float_t deltaArea  = jetArea[1]-jetArea[0];
 671
 672
 673     // fill thnsparse before acceptance cut
 674     if(fbJetsBeforeCut1){
 675       Double_t jetBeforeCutEntries1[10] = {
 676         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 677         (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1],
 678         (Double_t)jetPhi[0], (Double_t)jetPhi[1], (Double_t)fraction
 679       };
 680       fhnJetsBeforeCut1->Fill(jetBeforeCutEntries1);
 681     }
 682
 683     if(fbJetsBeforeCut2){
 684       Double_t jetBeforeCutEntries2[10] = {
 685         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 686         (Double_t)jetPt[0], (Double_t)jetPt[1],
 687         (Double_t)jetEta[0], (Double_t)jetEta[1],
 688         (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR
 689       };
 690       fhnJetsBeforeCut2->Fill(jetBeforeCutEntries2);
 691     }
 692
 693     // jet selection
 694     Bool_t jetAccepted = kTRUE;
 695     // minimum fraction required
 696     if(fraction<fJetPtFractionMin){
 697       fHistJetSelection->Fill(4);
 698       fh2JetSelection->Fill(4.,jetPt[0]);
 699       jetAccepted = kFALSE;
 700     }
 701
 702     if(jetAccepted){
 703       // jet acceptance + minimum pT check
 704       if(jetEta[0]>fJetEtaMax || jetEta[0]<fJetEtaMin ||
 705          jetEta[1]>fJetEtaMax || jetEta[1]<fJetEtaMin){
 706
 707         if(fDebug){
 708           Printf("Jet not in eta acceptance.");
 709           Printf("[0]: jet %d eta %.2f", ig, jetEta[0]);
 710           Printf("[1]: jet %d eta %.2f", ir, jetEta[1]);
 711         }
 712         fHistJetSelection->Fill(5);
 713         fh2JetSelection->Fill(5.,jetPt[0]);
 714         jetAccepted = kFALSE;
 715       }
 716     }
 717     if(jetAccepted){
 718       if(jetPt[0] < fJetPtMin){
 719         if(fDebug) Printf("Jet %d (pT %.1f GeV/c) has less than required pT.", ir, jetPt[0]);
 720         fHistJetSelection->Fill(6);
 721         fh2JetSelection->Fill(6.,jetPt[0]);
 722         jetAccepted = kFALSE;
 723       }
 724     }
 725
 726     if(jetAccepted){
 727       if(jet[1]->Trigger()&fJetTriggerExcludeMask){
 728         fHistJetSelection->Fill(7);
 729         fh2JetSelection->Fill(7.,jetPt[0]);
 730         jetAccepted = kFALSE;
 731       }
 732
 733     }
 734
 735     if(!jetAccepted){
 736       if(fbJetsMismatch2){
 737         Double_t jetEntriesMismatch2[10] = {
 738           (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 739           (Double_t)jetPt[0], (Double_t)jetPt[1],
 740           (Double_t)jetEta[0], (Double_t)jetEta[1],
 741           (Double_t)deltaEta, (Double_t)deltaR,
 742           (Double_t)fraction
 743         };
 744         fhnJetsMismatch2->Fill(jetEntriesMismatch2);
 745       }
 746       continue;
 747     }
 748
 749     // all accepted jets
 750     fHistJetSelection->Fill(0);
 751     fh2JetSelection->Fill(0.,jetPt[0]);
 752
 753     // fill thnsparse
 754     if(fbJetsRp){
 755       Double_t jetEntriesRp[11] = {
 756         (Double_t)centValue, (Double_t)nInputTracks, (Double_t) rp,
 757         (Double_t)rpBin, (Double_t)rpJet[0], (Double_t)rpJet[1],
 758         (Double_t)jetPt[0], (Double_t)deltaPt, (Double_t)rho, (Double_t)rpCorr,
 759         (Double_t)pthardbin
 760       };
 761       fhnJetsRp->Fill(jetEntriesRp);
 762     }
 763
 764     if(fbJetsDeltaPt){
 765       Double_t jetEntriesDeltaPt[8] = {
 766         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 767         (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)deltaPt, (Double_t)deltaArea,
 768         (Double_t)pthardbin
 769       };
 770       fhnJetsDeltaPt->Fill(jetEntriesDeltaPt);
 771     }
 772
 773     if(fbJetsEta){
 774       Double_t jetEntriesEta[11] = {
 775         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 776         (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1],
 777         (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR,
 778         (Double_t)pthardbin
 779       };
 780       fhnJetsEta->Fill(jetEntriesEta);
 781     }
 782
 783     if(fbJetsPhi){
 784       Double_t jetEntriesPhi[10] = {
 785         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 786         (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetPhi[0], (Double_t)jetPhi[1],
 787         (Double_t)deltaPt, (Double_t)deltaPhi,
 788         (Double_t)pthardbin
 789       };
 790       fhnJetsPhi->Fill(jetEntriesPhi);
 791     }
 792
 793     if(fbJetsArea){
 794       Double_t jetEntriesArea[14] = {
 795         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 796         (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetArea[0], (Double_t)jetArea[1],
 797         (Double_t)deltaPt, (Double_t)deltaR, (Double_t)deltaArea,
 798         (Double_t)fraction, (Double_t)distNextJet, (Double_t)ptNextJet,
 799         (Double_t)pthardbin
 800       };
 801       fhnJetsArea->Fill(jetEntriesArea);
 802     }
 803
 804     if(fbJets3Branches){
 805       Double_t jetEntries3Branches[14] = {
 806         (Double_t)centValue, (Double_t)nInputTracks,
 807         (Double_t)jetPt[0], (Double_t)jetPt[1],
 808         (Double_t)jetArea[0], (Double_t)jetArea[1],
 809         (Double_t)deltaPt, (Double_t)fraction, (Double_t)pthardbin,
 810         (Double_t)jetPt[2],(Double_t)delta,(Double_t)jetArea[2], (Double_t)fraction2, (Double_t)deltaR
 811       };
 812       fhnJets3Branches->Fill(jetEntries3Branches);
 813     }
 814
 815   }
 816
 817   PostData(1, fOutputList);
 818 }
 819
 820 void AliAnalysisTaskJetResponseV2::Terminate(const Option_t *)
 821 {
 822   // Draw result to the screen
 823   // Called once at the end of the query
 824
 825   if (!GetOutputData(1))
 826     return;
 827 }
 828
 829 Int_t AliAnalysisTaskJetResponseV2::GetNInputTracks()
 830 {
 831   // number of tracks in the event, obtained via jet finder
 832
 833   Int_t nInputTracks = 0;
 834
 835   TString jbname(fJetBranchName[1]);
 836   //needs complete event, use jets without background subtraction
 837   for(Int_t i=1; i<=3; ++i){
 838     if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
 839   }
 840   // use only HI event
 841   if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
 842   else if(jbname.Contains("AODMCextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
 843   if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
 844   else if(jbname.Contains("AODMCextra")) jbname.ReplaceAll("AODextra","AOD");
 845
 846   if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
 847   TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
 848   if(!tmpAODjets && fAODOut) tmpAODjets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(jbname.Data()));
 849   if(!tmpAODjets && fAODExtension) tmpAODjets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(jbname.Data()));
 850   if(!tmpAODjets){
 851     Printf("Jet branch %s not found", jbname.Data());
 852     Printf("AliAnalysisTaskJetResponseV2::GetNInputTracks FAILED");
 853     return -1;
 854   }
 855
 856   for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
 857     AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
 858     if(!jet) continue;
 859     TRefArray *trackList = jet->GetRefTracks();
 860     Int_t nTracks = trackList->GetEntriesFast();
 861     nInputTracks += nTracks;
 862     if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
 863   }
 864   if(fDebug) Printf("---> input tracks: %d", nInputTracks);
 865
 866   return nInputTracks;
 867 }
 868
 869 THnSparse* AliAnalysisTaskJetResponseV2::NewTHnSparseF(const char* name, UInt_t entries, UInt_t opt)
 870 {
 871   // generate new THnSparseF, axes are defined in GetDimParams()
 872
 873   Int_t count = 0;
 874   UInt_t tmp = entries;
 875   while(tmp!=0){
 876     count++;
 877     tmp = tmp &~ -tmp;  // clear lowest bit
 878   }
 879
 880   TString hnTitle(name);
 881   const Int_t dim = count;
 882   Int_t nbins[dim];
 883   Double_t xmin[dim];
 884   Double_t xmax[dim];
 885
 886   Int_t i=0;
 887   Int_t c=0;
 888   while(c<dim && i<32){
 889     if(entries&(1<<i)){
 890       Bool_t highres = opt&(1<<i);
 891       TString label("");
 892       GetDimParams(i, highres, label, nbins[c], xmin[c], xmax[c]);
 893       hnTitle += Form(";%s",label.Data());
 894       c++;
 895     }
 896
 897     i++;
 898   }
 899   hnTitle += ";";
 900
 901   return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
 902 }
 903
 904 void AliAnalysisTaskJetResponseV2::GetDimParams(Int_t iEntry, Bool_t hr, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
 905 {
 906   // stores label and binning of axis for THnSparse
 907
 908   const Double_t pi = TMath::Pi();
 909
 910   switch(iEntry){
 911
 912   case 0:
 913     label = "V0 centrality (%)";
 914     if(hr){
 915       nbins = 100;
 916       xmin = 0.;
 917       xmax = 100.;
 918     } else {
 919       nbins = 8;
 920       xmin = 0.;
 921       xmax = 80.;
 922     }
 923     break;
 924
 925
 926   case 1:
 927     label = "nb. of input tracks";
 928     if(fIsPbPb){
 929       if(hr){
 930         nbins = 400;
 931         xmin = 0.;
 932         xmax = 4000.;
 933       } else {
 934         nbins = 40;
 935         xmin = 0.;
 936         xmax = 4000.;
 937       }
 938     } else {
 939       nbins = 40;
 940       xmin = 0.;
 941       xmax = 400.;
 942     }
 943     break;
 944
 945
 946   case 2:
 947     label = "event plane #psi";
 948     if(hr){
 949       nbins = 30;
 950       xmin = 0.;
 951       xmax = pi;
 952     } else {
 953       nbins = 30;
 954       xmin = 0.;
 955       xmax = pi;
 956     }
 957     break;
 958
 959
 960   case 3:
 961     label = "event plane bin";
 962     nbins = 3;
 963     xmin = -.5;
 964     xmax = 2.5;
 965     break;
 966
 967
 968   case 4:
 969   case 5:
 970     if(iEntry==4)label = "#Delta#phi(RP-jet) (probe)";
 971     if(iEntry==5)label = "#Delta#phi(RP-jet) (rec)";
 972     nbins = 48;
 973     xmin = -pi;
 974     xmax =  pi;
 975     break;
 976
 977
 978   case 6:
 979   case 7:
 980   case 27:
 981     if(iEntry==6)label = "probe p_{T} (GeV/c)";
 982     if(iEntry==7 || iEntry==27)label = "rec p_{T} (GeV/c)";
 983     if(hr){
 984       nbins = 300;
 985       xmin = -50.;
 986       xmax = 250.;
 987     } else {
 988       nbins = 50;
 989       xmin = 0.;
 990       xmax = 250.;
 991     }
 992     break;
 993
 994
 995   case 8:
 996   case 9:
 997     if(iEntry==8)label = "probe #eta";
 998     if(iEntry==9)label = "rec #eta";
 999     if(hr){
1000       nbins = 56;
1001       xmin = -.7;
1002       xmax =  .7;
1003     } else {
1004       nbins = 28;
1005       xmin = -.7;
1006       xmax =  .7;
1007     }
1008     break;
1009
1010
1011   case 10:
1012   case 11:
1013     if(iEntry==10)label = "probe #phi";
1014     if(iEntry==11)label = "rec #phi";
1015     if(hr){
1016       nbins = 90;  // modulo 18 (sectors)
1017       xmin = 0.;
1018       xmax = 2*pi;
1019     } else {
1020       nbins = 25;
1021       xmin = 0.;
1022       xmax = 2*pi;
1023     }
1024     break;
1025
1026
1027   case 12:
1028   case 13:
1029   case 29:
1030     if(iEntry==12)label = "probe area";
1031     if(iEntry==13 || iEntry==29)label = "rec area";
1032     if(hr){
1033       nbins = 100;
1034       xmin = 0.;
1035       xmax = 1.;
1036     } else {
1037       nbins = 25;
1038       xmin = 0.;
1039       xmax = 1.;
1040     }
1041     break;
1042
1043   case 14:
1044   case 28:
1045     if(iEntry==14) label = "#delta p_{T}";
1046     if(iEntry==28) label = "#delta";
1047     if(hr){
1048       nbins = 241;
1049       xmin = -120.5;
1050       xmax =  120.5;
1051     } else {
1052       nbins = 101;
1053       xmin = -101.;
1054       xmax =  101.;
1055     }
1056     break;
1057
1058   case 15:
1059     label = "#Delta#eta";
1060     if(hr){
1061       nbins = 51;
1062       xmin = -1.02;
1063       xmax =  1.02;
1064     } else {
1065       nbins = 51;
1066       xmin = -1.02;
1067       xmax =  1.02;
1068     }
1069     break;
1070
1071
1072   case 16:
1073     label = "#Delta#phi";
1074     if(hr){
1075       nbins = 45;
1076       xmin = -pi;
1077       xmax =  pi;
1078     } else {
1079       nbins = 45;
1080       xmin = -pi;
1081       xmax =  pi;
1082     }
1083     break;
1084
1085
1086   case 17:
1087     label = "#DeltaR";
1088     if(hr){
1089       nbins = 50;
1090       xmin = 0.;
1091       xmax = 1.;
1092     } else {
1093       nbins = 50;
1094       xmin = 0.;
1095       xmax = 1.;
1096     }
1097     break;
1098
1099
1100   case 18:
1101     label = "#Deltaarea";
1102     if(hr){
1103       nbins = 81;
1104       xmin = -.81;
1105       xmax =  .81;
1106     } else {
1107       nbins = 33;
1108       xmin = -.825;
1109       xmax =  .825;
1110     }
1111     break;
1112
1113
1114   case 19:
1115   case 30:
1116     label = "fraction";
1117     if(hr){
1118       nbins = 52;
1119       xmin = 0.;
1120       xmax = 1.04;
1121     } else {
1122       nbins = 52;
1123       xmin = 0.;
1124       xmax = 1.04;
1125     }
1126     break;
1127
1128
1129   case 20:
1130     label = "distance to closest rec jet";
1131     if(hr){
1132       nbins = 51;
1133       xmin = -0.02;
1134       xmax =  1.;
1135     } else {
1136       nbins = 51;
1137       xmin = -0.02;
1138       xmax = 1.;
1139     }
1140     break;
1141
1142
1143   case 21:
1144     label = "p_{T} of closest rec jet";
1145     nbins = 100;
1146     xmin =  0.;
1147     xmax =  200.;
1148     break;
1149
1150   case 22:
1151     label = "#rho";
1152     nbins = 125;
1153     xmin  = 0.;
1154     xmax  = 250.;
1155     break;
1156
1157   case 23:
1158     label = "abs. correction of #rho for RP";
1159     nbins =  51;
1160     xmin  = -51.;
1161     xmax  =  51.;
1162     break;
1163
1164   case 24:
1165     label = "local #rho";
1166     nbins = 125;
1167     xmin  = 0.;
1168     xmax  = 250.;
1169     break;
1170
1171   case 25:
1172     label = "local #rho / #rho";
1173     nbins = 500;
1174     xmin  = 0.;
1175     xmax  = 5.;
1176     break;
1177
1178   case 26:
1179     label = "p_{T,hard} bin";
1180     nbins = 10;
1181     xmin  = -.5;
1182     xmax  = 9.5;
1183     break;
1184
1185   }
1186
1187 }
1188
1189 //____________________________________________________________________________
1190 Int_t AliAnalysisTaskJetResponseV2::GetPtHardBin(Double_t ptHard){
1191
1192   // returns pt hard bin (for LHC10e14, LHC11a1x, LHC11a2x simulations)
1193
1194   const Int_t nBins = 10;
1195   Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
1196
1197   Int_t bin = -1;
1198   while(bin<nBins-1 && binLimits[bin+1]<ptHard){
1199     bin++;
1200   }
1201
1202   return bin;
1203 }
1204
1205 //____________________________________________________________________________
1206 Double_t AliAnalysisTaskJetResponseV2::GetPt(AliAODJet *j, Int_t mode=0){
1207
1208   // returns jet pt, also negative pt after background subtraction if available
1209
1210   Double_t pt = 0.;
1211
1212   if(fKeepJets && mode>0){  // background subtracted pt, can be negative
1213     pt = j->GetPtSubtracted(0);
1214   }
1215   else{
1216     pt = j->Pt();  // if negative, pt=0.01
1217   }
1218
1219   return pt;
1220 }