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