PWG4/JetTasks/AliAnalysisTaskJetResponseV2.cxx

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