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