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