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