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<<27 | 1<<14 | 1<<28 | 1<<12 | 1<<13 | 1<<29 | 1<<19 | 1<<30 | 1<<26;
 323     opt = 1<<6 | 1<<7 | 1<<27 | 1<<14 | 1<<28;
 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     fkNbranches=3;
 488   }
 489
 490   for (Int_t iJetType = 0; iJetType < fkNbranches; iJetType++) {
 491     fListJets[iJetType]->Clear();
 492     if (!aodJets[iJetType]) continue;
 493
 494     if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
 495
 496     for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
 497       AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
 498       if (jet) fListJets[iJetType]->Add(jet);
 499     }
 500     fListJets[iJetType]->Sort();
 501   }
 502
 503
 504   // jet matching
 505   static TArrayI aMatchIndex(fListJets[0]->GetEntries());
 506   static TArrayF aPtFraction(fListJets[0]->GetEntries());
 507   if(aMatchIndex.GetSize()<fListJets[0]->GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries());
 508   if(aPtFraction.GetSize()<fListJets[0]->GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries());
 509
 510   // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction'
 511   AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
 512                                             fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()),
 513                                             aMatchIndex, aPtFraction, fDebug, fMatchMaxDist, fIsPbPb?1:2);
 514   static TArrayI aMatchIndexv2(fListJets[0]->GetEntries());
 515   static TArrayF aPtFractionv2(fListJets[0]->GetEntries());
 516   if(aMatchIndexv2.GetSize()<fListJets[0]->GetEntries()) aMatchIndexv2.Set(fListJets[0]->GetEntries());
 517   if(aPtFractionv2.GetSize()<fListJets[0]->GetEntries()) aPtFractionv2.Set(fListJets[0]->GetEntries());
 518   if( strlen(fJetBranchName[2].Data()) ) {
 519     AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
 520                                               fListJets[2], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[2]->GetEntries()),
 521                                               aMatchIndexv2, aPtFractionv2, fDebug, fMatchMaxDist, fIsPbPb?1:2);
 522   }
 523
 524   // loop over matched jets
 525   Int_t ir = -1;  // index of matched reconstruced jet
 526   Int_t ir2 = -1; // index of matched reconstruced jet
 527   Float_t fraction = -1.;
 528   Float_t fraction2 = -1.;
 529   AliAODJet *jet[3]  = { 0x0, 0x0, 0x0 };
 530   Float_t jetEta[3]  = { -990., -990., -990. };
 531   Float_t jetPhi[3]  = { -990., -990., -990. };
 532   Float_t jetPt[3]   = { -990., -990., -990. };
 533   Float_t jetArea[3] = { -990., -990., -990. };
 534   Float_t rpJet[3]   = { -990., -990., -990. };
 535   Int_t rpBin = -990;
 536
 537   for(Int_t ig=0; ig<fListJets[0]->GetEntries(); ++ig){
 538     ir = aMatchIndex[ig];
 539     if(aMatchIndexv2.GetSize()>ig) ir2 = aMatchIndexv2[ig];
 540     else ir2 =0;
 541
 542     //fetch jets
 543     jet[0] = (AliAODJet*)(fListJets[0]->At(ig));
 544     if(ir>=0) jet[1] = (AliAODJet*)(fListJets[1]->At(ir));
 545     else      jet[1] = 0x0;
 546     if(ir2>=0) jet[2] = (AliAODJet*)(fListJets[2]->At(ir2));
 547     else      jet[2] = 0x0;
 548
 549     for(Int_t i=0; i<fkNbranches; ++i){
 550       if(!jet[i]){
 551         jetEta[i]  = -990;
 552         jetPhi[i]  = -990.;
 553         jetPt[i]   = -990.;
 554         jetArea[i] = -990.;
 555         rpJet[i]   = -990.;
 556         if(i==1) rpBin = -990;
 557       }
 558       else {
 559         jetEta[i]  = jet[i]->Eta();
 560         jetPhi[i]  = jet[i]->Phi();
 561         jetPt[i]   = GetPt(jet[i], i);
 562         jetArea[i] = jet[i]->EffectiveAreaCharged();
 563         rpJet[i]   = TVector2::Phi_mpi_pi(rp-jetPhi[i]);
 564         if(i==1) rpBin = AliAnalysisHelperJetTasks::GetPhiBin(TVector2::Phi_mpi_pi(rp-jetPhi[i]), 3);
 565       }
 566     }
 567     fraction  = aPtFraction[ig];
 568     if(aPtFractionv2.GetSize()>ig) fraction2 = aPtFractionv2[ig];
 569     else fraction2 = 0.;
 570
 571     // jet statistics
 572     fHistJetSelection->Fill(1); // all probe jets
 573     if(jet[0]) fh2JetSelection->Fill(1.,jetPt[0]); // all probe jets
 574
 575     if(ir<0){
 576       fHistJetSelection->Fill(2);
 577       if(jet[0]) fh2JetSelection->Fill(2.,jetPt[0]);
 578
 579       if(fbJetsMismatch1){
 580         if(!jet[0]) continue;
 581         Double_t jetEntriesMismatch1[7] = {
 582           (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 583           (Double_t)jetPt[0], (Double_t)jetEta[0], (Double_t)jetPhi[0], (Double_t)jetArea[0]
 584         };
 585         fhnJetsMismatch1->Fill(jetEntriesMismatch1);
 586       }
 587
 588       continue;
 589     }
 590
 591     if(!jet[0] || !jet[1]){
 592       fHistJetSelection->Fill(3);
 593       if(jet[0]) fh2JetSelection->Fill(3.,jetPt[0]);
 594       continue;
 595     }
 596
 597     // look for distance to next rec jet
 598     Float_t distNextJet = -0.01; // no neighbor
 599     Float_t ptNextJet = -1.;
 600     for(Int_t r=0; r<fListJets[1]->GetEntries(); ++r){
 601       if(r==ir) continue;
 602       Float_t tmpDeltaR = jet[1]->DeltaR((AliAODJet*)fListJets[1]->At(r));
 603       if(distNextJet<0. || distNextJet>tmpDeltaR){
 604         distNextJet = tmpDeltaR;
 605         if(fKeepJets) ptNextJet   = ((AliAODJet*)fListJets[1]->At(r))->GetPtSubtracted(0);
 606         else          ptNextJet   = ((AliAODJet*)fListJets[1]->At(r))->Pt();
 607       }
 608     }
 609
 610
 611
 612     //Float_t localRho = jetArea[1]>0. ? (jetPt[1]+rho*jetArea[1] - jetPt[0]) / jetArea[1] : 0.;
 613     //Float_t relRho = rho>0. ? localRho / rho : 0.;
 614
 615
 616     // calculate parameters of associated jets
 617     /* from Leticia, kT clusterizer
 618        Float_t par0[4] = { 0.00409,       0.01229,      0.05,         0.26 };
 619        Float_t par1[4] = { -2.97035e-03, -2.03182e-03, -1.25702e-03, -9.95107e-04 };
 620        Float_t par2[4] = { 1.02865e-01,   1.49039e-01,  1.53910e-01,  1.51109e-01 };
 621     */
 622     // own, from embedded tracks
 623     Float_t par0[4] = { 0.02841,       0.05039,      0.09092,      0.24089     };
 624     Float_t par1[4] = { -4.26725e-04, -1.15273e-03, -1.56827e-03, -3.08003e-03 };
 625     Float_t par2[4] = { 4.95415e-02,   9.79538e-02,  1.32814e-01,  1.71743e-01 };
 626
 627     Float_t rpCorr = 0.;
 628
 629     if(eventClass>0&&eventClass<4){
 630       rpCorr = par0[eventClass-1] + par1[eventClass-1] * 2*TMath::Cos(rpJet[1]) + par2[eventClass-1] * 2*TMath::Cos(2*rpJet[1]);
 631       rpCorr *= rho * jetArea[1];
 632     }
 633
 634     Float_t delta      = jetPt[2]-jetPt[1];
 635     Float_t deltaPt    = jetPt[1]-jetPt[0];
 636     Float_t deltaEta   = jetEta[1]-jetEta[0];
 637     Float_t deltaPhi   = TVector2::Phi_mpi_pi(jetPhi[1]-jetPhi[0]);
 638     Float_t deltaR     = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi);
 639     Float_t deltaArea  = jetArea[1]-jetArea[0];
 640
 641
 642     // fill thnsparse before acceptance cut
 643     if(fbJetsBeforeCut1){
 644       Double_t jetBeforeCutEntries1[10] = {
 645         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 646         (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1],
 647         (Double_t)jetPhi[0], (Double_t)jetPhi[1], (Double_t)fraction
 648       };
 649       fhnJetsBeforeCut1->Fill(jetBeforeCutEntries1);
 650     }
 651
 652     if(fbJetsBeforeCut2){
 653       Double_t jetBeforeCutEntries2[10] = {
 654         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 655         (Double_t)jetPt[0], (Double_t)jetPt[1],
 656         (Double_t)jetEta[0], (Double_t)jetEta[1],
 657         (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR
 658       };
 659       fhnJetsBeforeCut2->Fill(jetBeforeCutEntries2);
 660     }
 661
 662     // jet selection
 663     Bool_t jetAccepted = kTRUE;
 664     // minimum fraction required
 665     if(fraction<fJetPtFractionMin){
 666       fHistJetSelection->Fill(4);
 667       fh2JetSelection->Fill(4.,jetPt[0]);
 668       jetAccepted = kFALSE;
 669     }
 670
 671     if(jetAccepted){
 672       // jet acceptance + minimum pT check
 673       if(jetEta[0]>fJetEtaMax || jetEta[0]<fJetEtaMin ||
 674          jetEta[1]>fJetEtaMax || jetEta[1]<fJetEtaMin){
 675
 676         if(fDebug){
 677           Printf("Jet not in eta acceptance.");
 678           Printf("[0]: jet %d eta %.2f", ig, jetEta[0]);
 679           Printf("[1]: jet %d eta %.2f", ir, jetEta[1]);
 680         }
 681         fHistJetSelection->Fill(5);
 682         fh2JetSelection->Fill(5.,jetPt[0]);
 683         jetAccepted = kFALSE;
 684       }
 685     }
 686     if(jetAccepted){
 687       if(jetPt[0] < fJetPtMin){
 688         if(fDebug) Printf("Jet %d (pT %.1f GeV/c) has less than required pT.", ir, jetPt[0]);
 689         fHistJetSelection->Fill(6);
 690         fh2JetSelection->Fill(6.,jetPt[0]);
 691         jetAccepted = kFALSE;
 692       }
 693     }
 694
 695     if(jetAccepted){
 696       if(jet[1]->Trigger()&fJetTriggerExcludeMask){
 697         fHistJetSelection->Fill(7);
 698         fh2JetSelection->Fill(7.,jetPt[0]);
 699         jetAccepted = kFALSE;
 700       }
 701
 702     }
 703
 704     if(!jetAccepted){
 705       if(fbJetsMismatch2){
 706         Double_t jetEntriesMismatch2[10] = {
 707           (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 708           (Double_t)jetPt[0], (Double_t)jetPt[1],
 709           (Double_t)jetEta[0], (Double_t)jetEta[1],
 710           (Double_t)deltaEta, (Double_t)deltaR,
 711           (Double_t)fraction
 712         };
 713         fhnJetsMismatch2->Fill(jetEntriesMismatch2);
 714       }
 715       continue;
 716     }
 717
 718     // all accepted jets
 719     fHistJetSelection->Fill(0);
 720     fh2JetSelection->Fill(0.,jetPt[0]);
 721
 722     // fill thnsparse
 723     if(fbJetsRp){
 724       Double_t jetEntriesRp[11] = {
 725         (Double_t)centValue, (Double_t)nInputTracks, (Double_t) rp,
 726         (Double_t)rpBin, (Double_t)rpJet[0], (Double_t)rpJet[1],
 727         (Double_t)jetPt[0], (Double_t)deltaPt, (Double_t)rho, (Double_t)rpCorr,
 728         (Double_t)pthardbin
 729       };
 730       fhnJetsRp->Fill(jetEntriesRp);
 731     }
 732
 733     if(fbJetsDeltaPt){
 734       Double_t jetEntriesDeltaPt[8] = {
 735         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 736         (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)deltaPt, (Double_t)deltaArea,
 737         (Double_t)pthardbin
 738       };
 739       fhnJetsDeltaPt->Fill(jetEntriesDeltaPt);
 740     }
 741
 742     if(fbJetsEta){
 743       Double_t jetEntriesEta[11] = {
 744         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 745         (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetEta[0], (Double_t)jetEta[1],
 746         (Double_t)deltaPt, (Double_t)deltaEta, (Double_t)deltaR,
 747         (Double_t)pthardbin
 748       };
 749       fhnJetsEta->Fill(jetEntriesEta);
 750     }
 751
 752     if(fbJetsPhi){
 753       Double_t jetEntriesPhi[10] = {
 754         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 755         (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetPhi[0], (Double_t)jetPhi[1],
 756         (Double_t)deltaPt, (Double_t)deltaPhi,
 757         (Double_t)pthardbin
 758       };
 759       fhnJetsPhi->Fill(jetEntriesPhi);
 760     }
 761
 762     if(fbJetsArea){
 763       Double_t jetEntriesArea[14] = {
 764         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 765         (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetArea[0], (Double_t)jetArea[1],
 766         (Double_t)deltaPt, (Double_t)deltaR, (Double_t)deltaArea,
 767         (Double_t)fraction, (Double_t)distNextJet, (Double_t)ptNextJet,
 768         (Double_t)pthardbin
 769       };
 770       fhnJetsArea->Fill(jetEntriesArea);
 771     }
 772
 773     if(fbJets3Branches){
 774       Double_t jetEntries3Branches[13] = {
 775         (Double_t)centValue, (Double_t)nInputTracks,
 776         (Double_t)jetPt[0], (Double_t)jetPt[1],
 777         (Double_t)jetArea[0], (Double_t)jetArea[1],
 778         (Double_t)deltaPt, (Double_t)fraction, (Double_t)pthardbin,
 779         (Double_t)jetPt[2],(Double_t)delta,(Double_t)jetArea[2], (Double_t)fraction2
 780       };
 781       fhnJets3Branches->Fill(jetEntries3Branches);
 782     }
 783
 784   }
 785
 786   PostData(1, fOutputList);
 787 }
 788
 789 void AliAnalysisTaskJetResponseV2::Terminate(const Option_t *)
 790 {
 791   // Draw result to the screen
 792   // Called once at the end of the query
 793
 794   if (!GetOutputData(1))
 795     return;
 796 }
 797
 798 Int_t AliAnalysisTaskJetResponseV2::GetNInputTracks()
 799 {
 800   // number of tracks in the event, obtained via jet finder
 801
 802   Int_t nInputTracks = 0;
 803
 804   TString jbname(fJetBranchName[1]);
 805   //needs complete event, use jets without background subtraction
 806   for(Int_t i=1; i<=3; ++i){
 807     if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
 808   }
 809   // use only HI event
 810   if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
 811   if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
 812
 813   if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
 814   TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
 815   if(!tmpAODjets){
 816     Printf("Jet branch %s not found", jbname.Data());
 817     Printf("AliAnalysisTaskJetResponseV2::GetNInputTracks FAILED");
 818     return -1;
 819   }
 820
 821   for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
 822     AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
 823     if(!jet) continue;
 824     TRefArray *trackList = jet->GetRefTracks();
 825     Int_t nTracks = trackList->GetEntriesFast();
 826     nInputTracks += nTracks;
 827     if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
 828   }
 829   if(fDebug) Printf("---> input tracks: %d", nInputTracks);
 830
 831   return nInputTracks;
 832 }
 833
 834 THnSparse* AliAnalysisTaskJetResponseV2::NewTHnSparseF(const char* name, UInt_t entries, UInt_t opt)
 835 {
 836   // generate new THnSparseF, axes are defined in GetDimParams()
 837
 838   Int_t count = 0;
 839   UInt_t tmp = entries;
 840   while(tmp!=0){
 841     count++;
 842     tmp = tmp &~ -tmp;  // clear lowest bit
 843   }
 844
 845   TString hnTitle(name);
 846   const Int_t dim = count;
 847   Int_t nbins[dim];
 848   Double_t xmin[dim];
 849   Double_t xmax[dim];
 850
 851   Int_t i=0;
 852   Int_t c=0;
 853   while(c<dim && i<32){
 854     if(entries&(1<<i)){
 855       Bool_t highres = opt&(1<<i);
 856       TString label("");
 857       GetDimParams(i, highres, label, nbins[c], xmin[c], xmax[c]);
 858       hnTitle += Form(";%s",label.Data());
 859       c++;
 860     }
 861
 862     i++;
 863   }
 864   hnTitle += ";";
 865
 866   return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
 867 }
 868
 869 void AliAnalysisTaskJetResponseV2::GetDimParams(Int_t iEntry, Bool_t hr, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
 870 {
 871   // stores label and binning of axis for THnSparse
 872
 873   const Double_t pi = TMath::Pi();
 874
 875   switch(iEntry){
 876
 877   case 0:
 878     label = "V0 centrality (%)";
 879     if(hr){
 880       nbins = 100;
 881       xmin = 0.;
 882       xmax = 100.;
 883     } else {
 884       nbins = 8;
 885       xmin = 0.;
 886       xmax = 80.;
 887     }
 888     break;
 889
 890
 891   case 1:
 892     label = "nb. of input tracks";
 893     if(fIsPbPb){
 894       if(hr){
 895         nbins = 400;
 896         xmin = 0.;
 897         xmax = 4000.;
 898       } else {
 899         nbins = 40;
 900         xmin = 0.;
 901         xmax = 4000.;
 902       }
 903     } else {
 904       nbins = 40;
 905       xmin = 0.;
 906       xmax = 400.;
 907     }
 908     break;
 909
 910
 911   case 2:
 912     label = "event plane #psi";
 913     if(hr){
 914       nbins = 30;
 915       xmin = 0.;
 916       xmax = pi;
 917     } else {
 918       nbins = 30;
 919       xmin = 0.;
 920       xmax = pi;
 921     }
 922     break;
 923
 924
 925   case 3:
 926     label = "event plane bin";
 927     nbins = 3;
 928     xmin = -.5;
 929     xmax = 2.5;
 930     break;
 931
 932
 933   case 4:
 934   case 5:
 935     if(iEntry==4)label = "#Delta#phi(RP-jet) (probe)";
 936     if(iEntry==5)label = "#Delta#phi(RP-jet) (rec)";
 937     nbins = 48;
 938     xmin = -pi;
 939     xmax =  pi;
 940     break;
 941
 942
 943   case 6:
 944   case 7:
 945   case 27:
 946     if(iEntry==6)label = "probe p_{T} (GeV/c)";
 947     if(iEntry==7 || iEntry==27)label = "rec p_{T} (GeV/c)";
 948     if(hr){
 949       nbins = 300;
 950       xmin = -50.;
 951       xmax = 250.;
 952     } else {
 953       nbins = 50;
 954       xmin = 0.;
 955       xmax = 250.;
 956     }
 957     break;
 958
 959
 960   case 8:
 961   case 9:
 962     if(iEntry==8)label = "probe #eta";
 963     if(iEntry==9)label = "rec #eta";
 964     if(hr){
 965       nbins = 56;
 966       xmin = -.7;
 967       xmax =  .7;
 968     } else {
 969       nbins = 28;
 970       xmin = -.7;
 971       xmax =  .7;
 972     }
 973     break;
 974
 975
 976   case 10:
 977   case 11:
 978     if(iEntry==10)label = "probe #phi";
 979     if(iEntry==11)label = "rec #phi";
 980     if(hr){
 981       nbins = 90;  // modulo 18 (sectors)
 982       xmin = 0.;
 983       xmax = 2*pi;
 984     } else {
 985       nbins = 25;
 986       xmin = 0.;
 987       xmax = 2*pi;
 988     }
 989     break;
 990
 991
 992   case 12:
 993   case 13:
 994   case 29:
 995     if(iEntry==12)label = "probe area";
 996     if(iEntry==13 || iEntry==29)label = "rec area";
 997     if(hr){
 998       nbins = 100;
 999       xmin = 0.;
1000       xmax = 1.;
1001     } else {
1002       nbins = 25;
1003       xmin = 0.;
1004       xmax = 1.;
1005     }
1006     break;
1007
1008   case 14:
1009   case 28:
1010     if(iEntry==14) label = "#delta p_{T}";
1011     if(iEntry==28) label = "#delta";
1012     if(hr){
1013       nbins = 241;
1014       xmin = -120.5;
1015       xmax =  120.5;
1016     } else {
1017       nbins = 101;
1018       xmin = -101.;
1019       xmax =  101.;
1020     }
1021     break;
1022
1023   case 15:
1024     label = "#Delta#eta";
1025     if(hr){
1026       nbins = 51;
1027       xmin = -1.02;
1028       xmax =  1.02;
1029     } else {
1030       nbins = 51;
1031       xmin = -1.02;
1032       xmax =  1.02;
1033     }
1034     break;
1035
1036
1037   case 16:
1038     label = "#Delta#phi";
1039     if(hr){
1040       nbins = 45;
1041       xmin = -pi;
1042       xmax =  pi;
1043     } else {
1044       nbins = 45;
1045       xmin = -pi;
1046       xmax =  pi;
1047     }
1048     break;
1049
1050
1051   case 17:
1052     label = "#DeltaR";
1053     if(hr){
1054       nbins = 50;
1055       xmin = 0.;
1056       xmax = 1.;
1057     } else {
1058       nbins = 50;
1059       xmin = 0.;
1060       xmax = 1.;
1061     }
1062     break;
1063
1064
1065   case 18:
1066     label = "#Deltaarea";
1067     if(hr){
1068       nbins = 81;
1069       xmin = -.81;
1070       xmax =  .81;
1071     } else {
1072       nbins = 33;
1073       xmin = -.825;
1074       xmax =  .825;
1075     }
1076     break;
1077
1078
1079   case 19:
1080   case 30:
1081     label = "fraction";
1082     if(hr){
1083       nbins = 52;
1084       xmin = 0.;
1085       xmax = 1.04;
1086     } else {
1087       nbins = 52;
1088       xmin = 0.;
1089       xmax = 1.04;
1090     }
1091     break;
1092
1093
1094   case 20:
1095     label = "distance to closest rec jet";
1096     if(hr){
1097       nbins = 51;
1098       xmin = -0.02;
1099       xmax =  1.;
1100     } else {
1101       nbins = 51;
1102       xmin = -0.02;
1103       xmax = 1.;
1104     }
1105     break;
1106
1107
1108   case 21:
1109     label = "p_{T} of closest rec jet";
1110     nbins = 100;
1111     xmin =  0.;
1112     xmax =  200.;
1113     break;
1114
1115   case 22:
1116     label = "#rho";
1117     nbins = 125;
1118     xmin  = 0.;
1119     xmax  = 250.;
1120     break;
1121
1122   case 23:
1123     label = "abs. correction of #rho for RP";
1124     nbins =  51;
1125     xmin  = -51.;
1126     xmax  =  51.;
1127     break;
1128
1129   case 24:
1130     label = "local #rho";
1131     nbins = 125;
1132     xmin  = 0.;
1133     xmax  = 250.;
1134     break;
1135
1136   case 25:
1137     label = "local #rho / #rho";
1138     nbins = 500;
1139     xmin  = 0.;
1140     xmax  = 5.;
1141     break;
1142
1143   case 26:
1144     label = "p_{T,hard} bin";
1145     nbins = 10;
1146     xmin  = -.5;
1147     xmax  = 9.5;
1148     break;
1149
1150   }
1151
1152 }
1153
1154 //____________________________________________________________________________
1155 Int_t AliAnalysisTaskJetResponseV2::GetPtHardBin(Double_t ptHard){
1156
1157   // returns pt hard bin (for LHC10e14, LHC11a1x, LHC11a2x simulations)
1158
1159   const Int_t nBins = 10;
1160   Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
1161
1162   Int_t bin = -1;
1163   while(bin<nBins-1 && binLimits[bin+1]<ptHard){
1164     bin++;
1165   }
1166
1167   return bin;
1168 }
1169
1170 //____________________________________________________________________________
1171 Double_t AliAnalysisTaskJetResponseV2::GetPt(AliAODJet *j, Int_t mode=0){
1172
1173   // returns jet pt, also negative pt after background subtraction if available
1174
1175   Double_t pt = 0.;
1176
1177   if(fKeepJets && mode>0){  // background subtracted pt, can be negative
1178     pt = j->GetPtSubtracted(0);
1179   }
1180   else{
1181     pt = j->Pt();  // if negative, pt=0.01
1182   }
1183
1184   return pt;
1185 }