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
  22 #include "AliAODEvent.h"
  23 #include "AliAODJet.h"
  24
  25 #include "AliAnalysisTaskJetResponseV2.h"
  26
  27 ClassImp(AliAnalysisTaskJetResponseV2)
  28
  29 AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2() :
  30   AliAnalysisTaskSE(),
  31   fESD(0x0),
  32   fAOD(0x0),
  33   fOfflineTrgMask(AliVEvent::kAny),
  34   fMinContribVtx(1),
  35   fVtxZMin(-8.),
  36   fVtxZMax(8.),
  37   fEvtClassMin(1),
  38   fEvtClassMax(4),
  39   fCentMin(0.),
  40   fCentMax(100.),
  41   fNInputTracksMin(0),
  42   fNInputTracksMax(-1),
  43   fReactionPlaneBin(-1),
  44   fJetEtaMin(-.5),
  45   fJetEtaMax(.5),
  46   fJetPtMin(20.),
  47   fJetPtFractionMin(0.5),
  48   fNMatchJets(4),
  49   fMatchMaxDist(0.8),
  50   fkNbranches(2),
  51   fkEvtClasses(12),
  52   fOutputList(0x0),
  53   fHistEvtSelection(0x0),
  54   fhnEvent(0x0),
  55   fhnJetsRp(0x0),
  56   fhnJetsDeltaPt(0x0),
  57   fhnJetsEta(0x0),
  58   fhnJetsArea(0x0),
  59   fhnJetsBeforeCut1(0x0),
  60   fhnJetsBeforeCut2(0x0)
  61 {
  62   // default Constructor
  63
  64   fJetBranchName[0] = "";
  65   fJetBranchName[1] = "";
  66
  67   fListJets[0] = new TList;
  68   fListJets[1] = new TList;
  69 }
  70
  71 AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2(const char *name) :
  72   AliAnalysisTaskSE(name),
  73   fESD(0x0),
  74   fAOD(0x0),
  75   fOfflineTrgMask(AliVEvent::kAny),
  76   fMinContribVtx(1),
  77   fVtxZMin(-8.),
  78   fVtxZMax(8.),
  79   fEvtClassMin(1),
  80   fEvtClassMax(4),
  81   fCentMin(0.),
  82   fCentMax(100.),
  83   fNInputTracksMin(0),
  84   fNInputTracksMax(-1),
  85   fReactionPlaneBin(-1),
  86   fJetEtaMin(-.5),
  87   fJetEtaMax(.5),
  88   fJetPtMin(20.),
  89   fJetPtFractionMin(0.5),
  90   fNMatchJets(4),
  91   fMatchMaxDist(0.8),
  92   fkNbranches(2),
  93   fkEvtClasses(12),
  94   fOutputList(0x0),
  95   fHistEvtSelection(0x0),
  96   fhnEvent(0x0),
  97   fhnJetsRp(0x0),
  98   fhnJetsDeltaPt(0x0),
  99   fhnJetsEta(0x0),
 100   fhnJetsArea(0x0),
 101   fhnJetsBeforeCut1(0x0),
 102   fhnJetsBeforeCut2(0x0)
 103 {
 104   // Constructor
 105
 106   fJetBranchName[0] = "";
 107   fJetBranchName[1] = "";
 108
 109   fListJets[0] = new TList;
 110   fListJets[1] = new TList;
 111
 112   DefineOutput(1, TList::Class());
 113 }
 114
 115 AliAnalysisTaskJetResponseV2::~AliAnalysisTaskJetResponseV2()
 116 {
 117   delete fListJets[0];
 118   delete fListJets[1];
 119 }
 120
 121 void AliAnalysisTaskJetResponseV2::SetBranchNames(const TString &branch1, const TString &branch2)
 122 {
 123   fJetBranchName[0] = branch1;
 124   fJetBranchName[1] = branch2;
 125 }
 126
 127 void AliAnalysisTaskJetResponseV2::Init()
 128 {
 129
 130    // check for jet branches
 131    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 132       AliError("Jet branch name not set.");
 133    }
 134
 135 }
 136
 137 void AliAnalysisTaskJetResponseV2::UserCreateOutputObjects()
 138 {
 139   // Create histograms
 140   // Called once
 141   OpenFile(1);
 142   if(!fOutputList) fOutputList = new TList;
 143   fOutputList->SetOwner(kTRUE);
 144
 145   Bool_t oldStatus = TH1::AddDirectoryStatus();
 146   TH1::AddDirectory(kFALSE);
 147
 148
 149   fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
 150   fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
 151   fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
 152   fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
 153   fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
 154   fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
 155   fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
 156
 157   Float_t pi = TMath::Pi();
 158
 159   const Int_t dim1 = 3;
 160   // cent : nInpTrks : rp
 161   Int_t    nbins1[dim1] = { 100,   400,  30  };
 162   Double_t xmin1[dim1]  = {   0.,    0.,  0. };
 163   Double_t xmax1[dim1]  = { 100., 4000., pi  };
 164
 165   TString hnTitle("variables per event;centrality;nb. of input tracks;reaction plane #phi");
 166
 167   fhnEvent = new THnSparseF("fhnEvent", hnTitle.Data(), dim1, nbins1, xmin1, xmax1);
 168
 169   /* original thnsparse
 170   const Int_t dim2 = 18;
 171   // cent : nInpTrks : rp bins : rp wrt jet : fraction :
 172   // jetPt(2x) : jetEta(2x) : jetPhi (2x) : jetArea(2x)
 173   // deltaPt : deltaEta : deltaPhi : deltaR : deltaArea
 174
 175   Int_t    nbins2[dim2] = {  16,   400,    3,  48,    52, 125,  125,    56,   56,   25,   25, 100, 100,    201,   101,    51,  50,   81 };
 176   Double_t xmin2[dim2]  = {   0.,    0., -.5, -pi,  0.  ,   0.,   0., -0.7, -0.7,   0.,   0., 0. , 0. , -100.5, -1.01, -1.02,  0., -.81 };
 177   Double_t xmax2[dim2]  = {  80., 4000., 2.5,  pi,  1.04, 250., 250.,  0.7,  0.7, 2*pi, 2*pi, 1.0, 1.0,  100.5,  1.01,  1.02,  1.,  .81 };
 178
 179   fhnJets = new THnSparseF("fhnJets", "variables per jet", dim2, nbins2, xmin2, xmax2);
 180   */
 181
 182
 183   const Int_t dim2 = 6;
 184   // cent : nInpTrks : rp bins : rp wrt jet : probe pT
 185   Int_t    nbins2[dim2] = {   8 ,    40,   3,  48,  48,  50 };
 186   Double_t xmin2[dim2]  = {   0.,    0., -.5, -pi, -pi,  0. };
 187   Double_t xmax2[dim2]  = {  80., 4000., 2.5,  pi,  pi, 250.};
 188   hnTitle = "variables per jet;centrality;nb. of input tracks; reaction plane bin;#Delta#phi(RP-jet);probe p_{T}";
 189   fhnJetsRp = new THnSparseF("fhnJetsRp", hnTitle.Data(), dim2, nbins2, xmin2, xmax2);
 190
 191
 192   const Int_t dim3 = 6;
 193   // cent : nInpTrks : rp bins:  deltaPt : jetPt(2x) (hr delta pt)
 194   Int_t    nbins3[dim3] = {  16,   400,    3,    241,  250,  250 };
 195   Double_t xmin3[dim3]  = {   0.,    0., -.5, -120.5,   0.,   0. };
 196   Double_t xmax3[dim3]  = {  80., 4000., 2.5,  120.5, 250., 250. };
 197   hnTitle = "variables per jet;centrality;nb. of input tracks; reaction plane bin;#delta p_{T};probe p_{T};rec p_{T}";
 198   fhnJetsDeltaPt = new THnSparseF("fhnJetsDeltaPt", hnTitle.Data(), dim3, nbins3, xmin3, xmax3);
 199
 200   const Int_t dim4 = 10;
 201   // cent : nInpTrks : rp bins : deltaPt : jetPt(2x) : deltaR : deltaEta : jetEta(2x) (hr for eta)
 202   Int_t    nbins4[dim4] = {   8 ,   40 ,   3,  101 ,  50 ,  50 , 50,    51,   56,   56 };
 203   Double_t xmin4[dim4]  = {   0.,    0., -.5, -101.,   0.,   0., 0., -1.02, -0.7, -0.7 };
 204   Double_t xmax4[dim4]  = {  80., 4000., 2.5,  101., 250., 250., 1.,  1.02,  0.7,  0.7 };
 205   hnTitle = "variables per jet;centrality;nb. of input tracks; reaction plane bin;#delta p_{T};probe p_{T};rec p_{T};#DeltaR;#Delta#eta;#eta(probe);#eta(rec)";
 206   fhnJetsEta = new THnSparseF("fhnJetsEta", hnTitle.Data(), dim4, nbins4, xmin4, xmax4);
 207
 208   const Int_t dim5 = 13;
 209   // cent : nInpTrks : rp bins : deltaArea : jetArea(2x) : deltaR : fraction : distance next rec jet : pT next jet : deltaPt : jetPt(2x) (hr for area)
 210   hnTitle = "variables per jet;centrality;nb. of input tracks; reaction plane bin;#Deltaarea;probe area;rec area;#DeltaR;fraction;distance to closest rec jet;p_{T} of closest rec jet;#delta p_{T};probe p_{T};rec p_{T}";
 211   Int_t    nbins5[dim5] = {   8 ,   40 ,   3,   81, 100, 100, 50,   52,   51 , 100 ,  101 ,  50 ,  50  };
 212   Double_t xmin5[dim5]  = {   0.,    0., -.5, -.81,  0.,  0., 0., 0.  , -0.02,   0., -101.,   0.,   0. };
 213   Double_t xmax5[dim5]  = {  80., 4000., 2.5,  .81,  1.,  1., 1., 1.04,  1.  , 200.,  101., 250., 250. };
 214   fhnJetsArea = new THnSparseF("fhnJetsArea", hnTitle.Data(), dim5, nbins5, xmin5, xmax5);
 215
 216
 217   //before cut
 218   const Int_t dim6 = 10;
 219   // cent : nInpTrks : rp bins : fraction : jetPt(2x) : jetEta(2x) : jetPhi(2x) (low resolution) (with fraction, eta, phi, pt cuts possible)
 220   Int_t    nbins6[dim6] = {  8 ,   40 ,   3, 52 ,  50 ,  50 ,   28,   28,   25,   25 };
 221   Double_t xmin6[dim6]  = {  0.,    0., -.5,  0.,   0.,   0., -0.7, -0.7,   0.,   0. };
 222   Double_t xmax6[dim6]  = { 80., 4000., 2.5, 1.04, 250., 250.,  0.7,  0.7, 2*pi, 2*pi };
 223   hnTitle = "variables before cut;centrality;nb. of input tracks; reaction plane bin;fraction;probe p_{T};rec p_{T};probe #eta; rec #eta;probe #phi;rec #phi";
 224   fhnJetsBeforeCut1 = new THnSparseF("fhnJetsBeforeCut1", hnTitle.Data(), dim6, nbins6, xmin6, xmax6);
 225
 226   const Int_t dim7 = 10;
 227   // cent : nInpTrks : rp bins : deltaPt : jetPt(2x) : deltaR : deltaEta : jetEta(2x) (low resolution)
 228   Int_t    nbins7[dim7] = {   8 ,   40 ,   3,  101 ,  50 ,  50 , 50,    51,   28,   28 };
 229   Double_t xmin7[dim7]  = {   0.,    0., -.5, -101.,   0.,   0., 0., -1.02, -0.7, -0.7 };
 230   Double_t xmax7[dim7]  = {  80., 4000., 2.5,  101., 250., 250., 1.,  1.02,  0.7,  0.7 };
 231   hnTitle = "variables before cut;centrality;nb. of input tracks; reaction plane bin;#delta p_{T};probe p_{T};rec p_{T};#Delta R;#Delta #eta;probe #eta;rec #eta";
 232   fhnJetsBeforeCut2 = new THnSparseF("fhnJetsBeforeCut2", hnTitle.Data(), dim7, nbins7, xmin7, xmax7);
 233
 234
 235
 236   fOutputList->Add(fHistEvtSelection);
 237   fOutputList->Add(fhnEvent);
 238   fOutputList->Add(fhnJetsRp);
 239   fOutputList->Add(fhnJetsDeltaPt);
 240   fOutputList->Add(fhnJetsEta);
 241   fOutputList->Add(fhnJetsArea);
 242   fOutputList->Add(fhnJetsBeforeCut1);
 243   fOutputList->Add(fhnJetsBeforeCut2);
 244
 245   // =========== Switch on Sumw2 for all histos ===========
 246   for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
 247     TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
 248     if (h1){
 249       h1->Sumw2();
 250       continue;
 251     }
 252         THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
 253     if (hn){
 254       hn->Sumw2();
 255     }
 256   }
 257   TH1::AddDirectory(oldStatus);
 258
 259   PostData(1, fOutputList);
 260 }
 261
 262 void AliAnalysisTaskJetResponseV2::UserExec(Option_t *)
 263 {
 264   // load events, apply event cuts, then compare leading jets
 265
 266   if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 267       AliError("Jet branch name not set.");
 268       return;
 269    }
 270
 271   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
 272   if (!fESD) {
 273     AliError("ESD not available");
 274     return;
 275   }
 276   fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
 277   if (!fAOD) {
 278     AliError("AOD not available");
 279     return;
 280   }
 281
 282   // -- event selection --
 283   fHistEvtSelection->Fill(1); // number of events before event selection
 284
 285   // physics selection
 286   AliInputEventHandler* inputHandler = (AliInputEventHandler*)
 287     ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 288   if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
 289     if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
 290     fHistEvtSelection->Fill(2);
 291     PostData(1, fOutputList);
 292     return;
 293   }
 294
 295   // vertex selection
 296   AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
 297   Int_t nTracksPrim = primVtx->GetNContributors();
 298   if ((nTracksPrim < fMinContribVtx) ||
 299       (primVtx->GetZ() < fVtxZMin) ||
 300       (primVtx->GetZ() > fVtxZMax) ){
 301     if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
 302     fHistEvtSelection->Fill(3);
 303     PostData(1, fOutputList);
 304     return;
 305   }
 306
 307   // event class selection (from jet helper task)
 308   Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
 309   if(fDebug) Printf("Event class %d", eventClass);
 310   if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
 311     fHistEvtSelection->Fill(4);
 312     PostData(1, fOutputList);
 313     return;
 314   }
 315
 316   // centrality selection
 317   AliCentrality *cent = fESD->GetCentrality();
 318   Float_t centValue = cent->GetCentralityPercentile("V0M");
 319   if(fDebug) printf("centrality: %f\n", centValue);
 320   if (centValue < fCentMin || centValue > fCentMax){
 321     fHistEvtSelection->Fill(4);
 322     PostData(1, fOutputList);
 323     return;
 324   }
 325
 326
 327   // multiplicity due to input tracks
 328   Int_t nInputTracks = GetNInputTracks();
 329
 330   if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
 331        fHistEvtSelection->Fill(5);
 332        PostData(1, fOutputList);
 333        return;
 334   }
 335
 336
 337   fHistEvtSelection->Fill(0); // accepted events
 338   // -- end event selection --
 339
 340   Double_t rp = AliAnalysisHelperJetTasks::ReactionPlane(kFALSE);
 341   Double_t eventEntries[3] = { (Double_t)centValue, (Double_t)nInputTracks, rp };
 342   fhnEvent->Fill(eventEntries);
 343
 344
 345   // fetch jets
 346   TClonesArray *aodJets[2];
 347   aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
 348   aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE
 349
 350   for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
 351     fListJets[iJetType]->Clear();
 352     if (!aodJets[iJetType]) continue;
 353
 354     for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
 355       AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
 356       if (jet) fListJets[iJetType]->Add(jet);
 357     }
 358     fListJets[iJetType]->Sort();
 359   }
 360
 361
 362   // jet matching
 363   static TArrayI aMatchIndex(fListJets[0]->GetEntries());
 364   static TArrayF aPtFraction(fListJets[0]->GetEntries());
 365   if(aMatchIndex.GetSize()<fListJets[0]->GetEntries()) aMatchIndex.Set(fListJets[0]->GetEntries());
 366   if(aPtFraction.GetSize()<fListJets[0]->GetEntries()) aPtFraction.Set(fListJets[0]->GetEntries());
 367
 368   // stores matched jets in 'aMatchIndex' and fraction of pT in 'aPtFraction'
 369   AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()),
 370                                             fListJets[1], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[1]->GetEntries()),
 371                                             aMatchIndex, aPtFraction, fDebug, fMatchMaxDist);
 372
 373   // loop over matched jets
 374   Int_t ir = -1; // index of matched reconstruced jet
 375   Float_t fraction = -1.;
 376   AliAODJet *jet[2]  = { 0x0, 0x0 };
 377   Float_t jetEta[2]  = { -990., -990. };
 378   Float_t jetPhi[2]  = { -990., -990. };
 379   Float_t jetPt[2]   = { -990., -990. };
 380   Float_t jetArea[2] = { -990., -990. };
 381   Float_t rpJet[2]   = { -990., -990. };
 382
 383   for(Int_t ig=0; ig<fListJets[0]->GetEntries(); ++ig){
 384      ir = aMatchIndex[ig];
 385          if(ir<0) continue;
 386          fraction = aPtFraction[ig];
 387
 388          // fetch jets
 389          jet[0] = (AliAODJet*)(fListJets[0]->At(ig));
 390          jet[1] = (AliAODJet*)(fListJets[1]->At(ir));
 391          if(!jet[0] || !jet[1]) continue;
 392
 393          // look for distance to next rec jet
 394          Float_t distNextJet = -0.01; // no neighbor
 395          Float_t ptNextJet = -1.;
 396          for(Int_t r=0; r<fListJets[1]->GetEntries(); ++r){
 397             if(r==ir) continue;
 398                 Float_t tmpDeltaR = jet[1]->DeltaR((AliAODJet*)fListJets[1]->At(r));
 399                 if(distNextJet<0. || distNextJet>tmpDeltaR){
 400                 distNextJet = tmpDeltaR;
 401                     ptNextJet   = ((AliAODJet*)fListJets[1]->At(r))->Pt();
 402                 }
 403          }
 404
 405          // get parameters
 406          for(Int_t i=0; i<fkNbranches; ++i){
 407             jetEta[i]  = jet[i]->Eta();
 408                 jetPhi[i]  = jet[i]->Phi();
 409                 jetPt[i]   = jet[i]->Pt();
 410                 jetArea[i] = jet[i]->EffectiveAreaCharged();
 411                 rpJet[i]   = TVector2::Phi_mpi_pi(rp-jetPhi[i]);
 412          }
 413          Int_t rpBin = AliAnalysisHelperJetTasks::GetPhiBin(TVector2::Phi_mpi_pi(rp-jetPhi[1]), 3);
 414
 415
 416          // calculate parameters of associated jets
 417          Float_t deltaPt    = jetPt[1]-jetPt[0];
 418          Float_t deltaEta   = jetEta[1]-jetEta[0];
 419          Float_t deltaPhi   = TVector2::Phi_mpi_pi(jetPhi[1]-jetPhi[0]);
 420          Float_t deltaR     = TMath::Sqrt(deltaEta*deltaEta+deltaPhi*deltaPhi);
 421          Float_t deltaArea  = jetArea[1]-jetArea[0];
 422
 423
 424          // fill thnsparse before acceptance cut
 425          Double_t jetBeforeCutEntries1[10] = {
 426             (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 427                 (Double_t)fraction,     (Double_t)jetPt[0], (Double_t)jetPt[1],
 428                 (Double_t)jetEta[0], (Double_t)jetEta[1], (Double_t)jetPhi[0], (Double_t)jetPhi[1] };
 429
 430          Double_t jetBeforeCutEntries2[10] = {
 431         (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 432                 (Double_t)deltaPt, (Double_t)jetPt[0], (Double_t)jetPt[1],
 433                 (Double_t)deltaR, (Double_t)deltaEta,
 434                 (Double_t)jetEta[0], (Double_t)jetEta[1] };
 435
 436          fhnJetsBeforeCut1->Fill(jetBeforeCutEntries1);
 437          fhnJetsBeforeCut2->Fill(jetBeforeCutEntries2);
 438
 439
 440          // minimum fraction required
 441          if(fraction<fJetPtFractionMin) continue;
 442
 443          // jet acceptance + minimum pT check
 444          if(jetEta[0]>fJetEtaMax || jetEta[0]<fJetEtaMin ||
 445             jetEta[1]>fJetEtaMax || jetEta[1]<fJetEtaMin){
 446
 447                 if(fDebug){
 448                 Printf("Jet not in eta acceptance.");
 449                         Printf("[0]: jet %d eta %.2f", ig, jetEta[0]);
 450                         Printf("[1]: jet %d eta %.2f", ir, jetEta[1]);
 451                 }
 452                 continue;
 453      }
 454          if(jetPt[1] < fJetPtMin){
 455             if(fDebug) Printf("Jet %d (pT %.1f GeV/c) has less than required pT.", ir, jetPt[1]);
 456             continue;
 457          }
 458
 459          // fill thnsparse
 460          Double_t jetEntriesRp[6] = {
 461           (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 462           (Double_t)rpJet[0], (Double_t)rpJet[1], (Double_t)jetPt[0]
 463           };
 464      fhnJetsRp->Fill(jetEntriesRp);
 465
 466          Double_t jetEntriesDeltaPt[6] = {
 467            (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 468                    (Double_t)deltaPt, (Double_t)jetPt[0], (Double_t)jetPt[1]
 469                    };
 470      fhnJetsDeltaPt->Fill(jetEntriesDeltaPt);
 471
 472          Double_t jetEntriesEta[10] = {
 473            (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 474                    (Double_t)deltaPt, (Double_t)jetPt[0], (Double_t)jetPt[1],
 475                    (Double_t)deltaR, (Double_t)deltaEta, (Double_t)jetEta[0], (Double_t)jetEta[1]
 476                    };
 477      fhnJetsEta->Fill(jetEntriesEta);
 478
 479          Double_t jetEntriesArea[13] = {
 480            (Double_t)centValue, (Double_t)nInputTracks, (Double_t)rpBin,
 481            (Double_t)deltaArea, (Double_t)jetArea[0], (Double_t)jetArea[1],
 482                    (Double_t)deltaR, (Double_t)fraction, (Double_t)distNextJet, (Double_t)ptNextJet,
 483                    (Double_t)deltaPt, (Double_t)jetPt[0], (Double_t)jetPt[1]
 484                    };
 485      fhnJetsArea->Fill(jetEntriesArea);
 486
 487   }
 488
 489   PostData(1, fOutputList);
 490 }
 491
 492 void AliAnalysisTaskJetResponseV2::Terminate(const Option_t *)
 493 {
 494   // Draw result to the screen
 495   // Called once at the end of the query
 496
 497   if (!GetOutputData(1))
 498     return;
 499 }
 500
 501 Int_t AliAnalysisTaskJetResponseV2::GetNInputTracks()
 502 {
 503
 504   Int_t nInputTracks = 0;
 505
 506   TString jbname(fJetBranchName[1]);
 507   //needs complete event, use jets without background subtraction
 508   for(Int_t i=1; i<=3; ++i){
 509     if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
 510   }
 511   // use only HI event
 512   if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
 513   if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
 514
 515   if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
 516   TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
 517   if(!tmpAODjets){
 518     Printf("Jet branch %s not found", jbname.Data());
 519         Printf("AliAnalysisTaskJetResponseV2::GetNInputTracks FAILED");
 520         return -1;
 521   }
 522
 523   for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
 524       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
 525           if(!jet) continue;
 526           TRefArray *trackList = jet->GetRefTracks();
 527           Int_t nTracks = trackList->GetEntriesFast();
 528           nInputTracks += nTracks;
 529           if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
 530   }
 531   if(fDebug) Printf("---> input tracks: %d", nInputTracks);
 532
 533   return nInputTracks;
 534 }