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   else if(jbname.Contains("AODMCextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
 812   if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
 813   else if(jbname.Contains("AODMCextra")) jbname.ReplaceAll("AODextra","AOD");
 814
 815   if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
 816   TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
 817   if(!tmpAODjets){
 818     Printf("Jet branch %s not found", jbname.Data());
 819     Printf("AliAnalysisTaskJetResponseV2::GetNInputTracks FAILED");
 820     return -1;
 821   }
 822
 823   for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
 824     AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
 825     if(!jet) continue;
 826     TRefArray *trackList = jet->GetRefTracks();
 827     Int_t nTracks = trackList->GetEntriesFast();
 828     nInputTracks += nTracks;
 829     if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
 830   }
 831   if(fDebug) Printf("---> input tracks: %d", nInputTracks);
 832
 833   return nInputTracks;
 834 }
 835
 836 THnSparse* AliAnalysisTaskJetResponseV2::NewTHnSparseF(const char* name, UInt_t entries, UInt_t opt)
 837 {
 838   // generate new THnSparseF, axes are defined in GetDimParams()
 839
 840   Int_t count = 0;
 841   UInt_t tmp = entries;
 842   while(tmp!=0){
 843     count++;
 844     tmp = tmp &~ -tmp;  // clear lowest bit
 845   }
 846
 847   TString hnTitle(name);
 848   const Int_t dim = count;
 849   Int_t nbins[dim];
 850   Double_t xmin[dim];
 851   Double_t xmax[dim];
 852
 853   Int_t i=0;
 854   Int_t c=0;
 855   while(c<dim && i<32){
 856     if(entries&(1<<i)){
 857       Bool_t highres = opt&(1<<i);
 858       TString label("");
 859       GetDimParams(i, highres, label, nbins[c], xmin[c], xmax[c]);
 860       hnTitle += Form(";%s",label.Data());
 861       c++;
 862     }
 863
 864     i++;
 865   }
 866   hnTitle += ";";
 867
 868   return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
 869 }
 870
 871 void AliAnalysisTaskJetResponseV2::GetDimParams(Int_t iEntry, Bool_t hr, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
 872 {
 873   // stores label and binning of axis for THnSparse
 874
 875   const Double_t pi = TMath::Pi();
 876
 877   switch(iEntry){
 878
 879   case 0:
 880     label = "V0 centrality (%)";
 881     if(hr){
 882       nbins = 100;
 883       xmin = 0.;
 884       xmax = 100.;
 885     } else {
 886       nbins = 8;
 887       xmin = 0.;
 888       xmax = 80.;
 889     }
 890     break;
 891
 892
 893   case 1:
 894     label = "nb. of input tracks";
 895     if(fIsPbPb){
 896       if(hr){
 897         nbins = 400;
 898         xmin = 0.;
 899         xmax = 4000.;
 900       } else {
 901         nbins = 40;
 902         xmin = 0.;
 903         xmax = 4000.;
 904       }
 905     } else {
 906       nbins = 40;
 907       xmin = 0.;
 908       xmax = 400.;
 909     }
 910     break;
 911
 912
 913   case 2:
 914     label = "event plane #psi";
 915     if(hr){
 916       nbins = 30;
 917       xmin = 0.;
 918       xmax = pi;
 919     } else {
 920       nbins = 30;
 921       xmin = 0.;
 922       xmax = pi;
 923     }
 924     break;
 925
 926
 927   case 3:
 928     label = "event plane bin";
 929     nbins = 3;
 930     xmin = -.5;
 931     xmax = 2.5;
 932     break;
 933
 934
 935   case 4:
 936   case 5:
 937     if(iEntry==4)label = "#Delta#phi(RP-jet) (probe)";
 938     if(iEntry==5)label = "#Delta#phi(RP-jet) (rec)";
 939     nbins = 48;
 940     xmin = -pi;
 941     xmax =  pi;
 942     break;
 943
 944
 945   case 6:
 946   case 7:
 947   case 27:
 948     if(iEntry==6)label = "probe p_{T} (GeV/c)";
 949     if(iEntry==7 || iEntry==27)label = "rec p_{T} (GeV/c)";
 950     if(hr){
 951       nbins = 300;
 952       xmin = -50.;
 953       xmax = 250.;
 954     } else {
 955       nbins = 50;
 956       xmin = 0.;
 957       xmax = 250.;
 958     }
 959     break;
 960
 961
 962   case 8:
 963   case 9:
 964     if(iEntry==8)label = "probe #eta";
 965     if(iEntry==9)label = "rec #eta";
 966     if(hr){
 967       nbins = 56;
 968       xmin = -.7;
 969       xmax =  .7;
 970     } else {
 971       nbins = 28;
 972       xmin = -.7;
 973       xmax =  .7;
 974     }
 975     break;
 976
 977
 978   case 10:
 979   case 11:
 980     if(iEntry==10)label = "probe #phi";
 981     if(iEntry==11)label = "rec #phi";
 982     if(hr){
 983       nbins = 90;  // modulo 18 (sectors)
 984       xmin = 0.;
 985       xmax = 2*pi;
 986     } else {
 987       nbins = 25;
 988       xmin = 0.;
 989       xmax = 2*pi;
 990     }
 991     break;
 992
 993
 994   case 12:
 995   case 13:
 996   case 29:
 997     if(iEntry==12)label = "probe area";
 998     if(iEntry==13 || iEntry==29)label = "rec area";
 999     if(hr){
1000       nbins = 100;
1001       xmin = 0.;
1002       xmax = 1.;
1003     } else {
1004       nbins = 25;
1005       xmin = 0.;
1006       xmax = 1.;
1007     }
1008     break;
1009
1010   case 14:
1011   case 28:
1012     if(iEntry==14) label = "#delta p_{T}";
1013     if(iEntry==28) label = "#delta";
1014     if(hr){
1015       nbins = 241;
1016       xmin = -120.5;
1017       xmax =  120.5;
1018     } else {
1019       nbins = 101;
1020       xmin = -101.;
1021       xmax =  101.;
1022     }
1023     break;
1024
1025   case 15:
1026     label = "#Delta#eta";
1027     if(hr){
1028       nbins = 51;
1029       xmin = -1.02;
1030       xmax =  1.02;
1031     } else {
1032       nbins = 51;
1033       xmin = -1.02;
1034       xmax =  1.02;
1035     }
1036     break;
1037
1038
1039   case 16:
1040     label = "#Delta#phi";
1041     if(hr){
1042       nbins = 45;
1043       xmin = -pi;
1044       xmax =  pi;
1045     } else {
1046       nbins = 45;
1047       xmin = -pi;
1048       xmax =  pi;
1049     }
1050     break;
1051
1052
1053   case 17:
1054     label = "#DeltaR";
1055     if(hr){
1056       nbins = 50;
1057       xmin = 0.;
1058       xmax = 1.;
1059     } else {
1060       nbins = 50;
1061       xmin = 0.;
1062       xmax = 1.;
1063     }
1064     break;
1065
1066
1067   case 18:
1068     label = "#Deltaarea";
1069     if(hr){
1070       nbins = 81;
1071       xmin = -.81;
1072       xmax =  .81;
1073     } else {
1074       nbins = 33;
1075       xmin = -.825;
1076       xmax =  .825;
1077     }
1078     break;
1079
1080
1081   case 19:
1082   case 30:
1083     label = "fraction";
1084     if(hr){
1085       nbins = 52;
1086       xmin = 0.;
1087       xmax = 1.04;
1088     } else {
1089       nbins = 52;
1090       xmin = 0.;
1091       xmax = 1.04;
1092     }
1093     break;
1094
1095
1096   case 20:
1097     label = "distance to closest rec jet";
1098     if(hr){
1099       nbins = 51;
1100       xmin = -0.02;
1101       xmax =  1.;
1102     } else {
1103       nbins = 51;
1104       xmin = -0.02;
1105       xmax = 1.;
1106     }
1107     break;
1108
1109
1110   case 21:
1111     label = "p_{T} of closest rec jet";
1112     nbins = 100;
1113     xmin =  0.;
1114     xmax =  200.;
1115     break;
1116
1117   case 22:
1118     label = "#rho";
1119     nbins = 125;
1120     xmin  = 0.;
1121     xmax  = 250.;
1122     break;
1123
1124   case 23:
1125     label = "abs. correction of #rho for RP";
1126     nbins =  51;
1127     xmin  = -51.;
1128     xmax  =  51.;
1129     break;
1130
1131   case 24:
1132     label = "local #rho";
1133     nbins = 125;
1134     xmin  = 0.;
1135     xmax  = 250.;
1136     break;
1137
1138   case 25:
1139     label = "local #rho / #rho";
1140     nbins = 500;
1141     xmin  = 0.;
1142     xmax  = 5.;
1143     break;
1144
1145   case 26:
1146     label = "p_{T,hard} bin";
1147     nbins = 10;
1148     xmin  = -.5;
1149     xmax  = 9.5;
1150     break;
1151
1152   }
1153
1154 }
1155
1156 //____________________________________________________________________________
1157 Int_t AliAnalysisTaskJetResponseV2::GetPtHardBin(Double_t ptHard){
1158
1159   // returns pt hard bin (for LHC10e14, LHC11a1x, LHC11a2x simulations)
1160
1161   const Int_t nBins = 10;
1162   Double_t binLimits[nBins] = { 5., 11., 21., 36., 57., 84., 117., 156., 200., 249. }; // lower limits
1163
1164   Int_t bin = -1;
1165   while(bin<nBins-1 && binLimits[bin+1]<ptHard){
1166     bin++;
1167   }
1168
1169   return bin;
1170 }
1171
1172 //____________________________________________________________________________
1173 Double_t AliAnalysisTaskJetResponseV2::GetPt(AliAODJet *j, Int_t mode=0){
1174
1175   // returns jet pt, also negative pt after background subtraction if available
1176
1177   Double_t pt = 0.;
1178
1179   if(fKeepJets && mode>0){  // background subtracted pt, can be negative
1180     pt = j->GetPtSubtracted(0);
1181   }
1182   else{
1183     pt = j->Pt();  // if negative, pt=0.01
1184   }
1185
1186   return pt;
1187 }