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