PWGJE/AliAnalysisTaskJetCore.cxx

   1
   2 // ******************************************
   3 // This task computes several jet observables like
   4 // the fraction of energy in inner and outer coronnas,
   5 // jet-track correlations,triggered jet shapes and
   6 // correlation strength distribution of particles inside jets.
   7 // Author: lcunquei@cern.ch
   8 // *******************************************
   9
  10
  11 /**************************************************************************
  12  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  13  *                                                                        *
  14  * Author: The ALICE Off-line Project.                                    *
  15  * Contributors are mentioned in the code where appropriate.              *
  16  *                                                                        *
  17  * Permission to use, copy, modify and distribute this software and its   *
  18  * documentation strictly for non-commercial purposes is hereby granted   *
  19  * without fee, provided that the above copyright notice appears in all   *
  20  * copies and that both the copyright notice and this permission notice   *
  21  * appear in the supporting documentation. The authors make no claims     *
  22  * about the suitability of this software for any purpose. It is          *
  23  * provided "as is" without express or implied warranty.                  *
  24  **************************************************************************/
  25
  26
  27 #include "TChain.h"
  28 #include "TTree.h"
  29 #include "TMath.h"
  30 #include "TH1F.h"
  31 #include "TH2F.h"
  32 #include "TH3F.h"
  33 #include "THnSparse.h"
  34 #include "TCanvas.h"
  35
  36 #include "AliLog.h"
  37
  38 #include "AliAnalysisTask.h"
  39 #include "AliAnalysisManager.h"
  40
  41 #include "AliVEvent.h"
  42 #include "AliESDEvent.h"
  43 #include "AliESDInputHandler.h"
  44 #include "AliCentrality.h"
  45 #include "AliAnalysisHelperJetTasks.h"
  46 #include "AliInputEventHandler.h"
  47 #include "AliAODJetEventBackground.h"
  48 #include "AliAnalysisTaskFastEmbedding.h"
  49 #include "AliAODEvent.h"
  50 #include "AliAODHandler.h"
  51 #include "AliAODJet.h"
  52
  53 #include "AliAnalysisTaskJetCore.h"
  54
  55 ClassImp(AliAnalysisTaskJetCore)
  56
  57 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore() :
  58 AliAnalysisTaskSE(),
  59 fESD(0x0),
  60 fAODIn(0x0),
  61 fAODOut(0x0),
  62 fAODExtension(0x0),
  63 fBackgroundBranch(""),
  64 fNonStdFile(""),
  65 fIsPbPb(kTRUE),
  66 fOfflineTrgMask(AliVEvent::kAny),
  67 fMinContribVtx(1),
  68 fVtxZMin(-10.),
  69 fVtxZMax(10.),
  70 fEvtClassMin(0),
  71 fEvtClassMax(4),
  72 fFilterMask(0),
  73 fRadioFrac(0.2),
  74 fMinDist(0.1),
  75 fCentMin(0.),
  76 fCentMax(100.),
  77 fNInputTracksMin(0),
  78 fNInputTracksMax(-1),
  79 fAngStructCloseTracks(0),
  80 fCheckMethods(0),
  81 fDoEventMixing(0),
  82 fJetEtaMin(-.5),
  83 fJetEtaMax(.5),
  84 fNevents(0x0),
  85 fTindex(0x0),
  86 fTrigBufferIndex(0x0),
  87 fJetPtMin(20.),
  88 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
  89 fJetPtFractionMin(0.5),
  90 fNMatchJets(4),
  91 fMatchMaxDist(0.8),
  92 fKeepJets(kFALSE),
  93 fkNbranches(2),
  94 fkEvtClasses(12),
  95 fOutputList(0x0),
  96 fbEvent(kTRUE),
  97 fHistEvtSelection(0x0),
  98 fhnDeltaR(0x0),
  99 fhnMixedEvents(0x0),
 100 fh2JetCoreMethod1C10(0x0),
 101 fh2JetCoreMethod2C10(0x0),
 102 fh2JetCoreMethod1C20(0x0),
 103 fh2JetCoreMethod2C20(0x0),
 104 fh2JetCoreMethod1C30(0x0),
 105 fh2JetCoreMethod2C30(0x0),
 106 fh2JetCoreMethod1C60(0x0),
 107 fh2JetCoreMethod2C60(0x0),
 108 fh2JetTrackC10(0x0),
 109 fh2JetTrackC20(0x0),
 110 fh2AngStructpt1C10(0x0),
 111 fh2AngStructpt2C10(0x0),
 112 fh2AngStructpt3C10(0x0),
 113 fh2AngStructpt4C10(0x0),
 114 fh2AngStructpt1C20(0x0),
 115 fh2AngStructpt2C20(0x0),
 116 fh2AngStructpt3C20(0x0),
 117 fh2AngStructpt4C20(0x0),
 118 fh2AngStructpt1C30(0x0),
 119 fh2AngStructpt2C30(0x0),
 120 fh2AngStructpt3C30(0x0),
 121 fh2AngStructpt4C30(0x0),
 122 fh2AngStructpt1C60(0x0),
 123 fh2AngStructpt2C60(0x0),
 124 fh2AngStructpt3C60(0x0),
 125 fh2AngStructpt4C60(0x0),
 126 fh2JetsumHT3R2a(0x0),
 127 fh2JetsumHT3R2ap(0x0),
 128 fh2JetsumHT3R4a(0x0),
 129 fh2JetsumHT3R4ap(0x0),
 130 fh2JetsumHT3R6a(0x0),
 131 fh2JetsumHT3R6ap(0x0),
 132 fh2JetsumHT3R8a(0x0),
 133 fh2JetsumHT3R8ap(0x0),
 134 fh2JetsumHT3R10a(0x0),
 135 fh2JetsumHT3R10ap(0x0),
 136 fh2JetsumHT3R2aa(0x0),
 137 fh2JetsumHT3R2aap(0x0),
 138 fh2JetsumHT3R4aa(0x0),
 139 fh2JetsumHT3R4aap(0x0),
 140 fh2JetsumHT3R6aa(0x0),
 141 fh2JetsumHT3R6aap(0x0),
 142 fh2JetsumHT3R8aa(0x0),
 143 fh2JetsumHT3R8aap(0x0),
 144 fh2JetsumHT3R10aa(0x0),
 145 fh2JetsumHT3R10aap(0x0),
 146 fh2JetsumHT3R2aaa(0x0),
 147 fh2JetsumHT3R2aaap(0x0),
 148 fh2JetsumHT3R4aaa(0x0),
 149 fh2JetsumHT3R4aaap(0x0),
 150 fh2JetsumHT3R6aaa(0x0),
 151 fh2JetsumHT3R6aaap(0x0),
 152 fh2JetsumHT3R8aaa(0x0),
 153 fh2JetsumHT3R8aaap(0x0),
 154 fh2JetsumHT3R10aaa(0x0),
 155 fh2JetsumHT3R10aaap(0x0),
 156 fh2JetsumHT3R2b(0x0),
 157 fh2JetsumHT3R2bp(0x0),
 158 fh2JetsumHT3R4b(0x0),
 159 fh2JetsumHT3R4bp(0x0),
 160 fh2JetsumHT3R6b(0x0),
 161 fh2JetsumHT3R6bp(0x0),
 162 fh2JetsumHT3R8b(0x0),
 163 fh2JetsumHT3R8bp(0x0),
 164 fh2JetsumHT3R10b(0x0),
 165 fh2JetsumHT3R10bp(0x0),
 166 fh2JetsumHT3R2bb(0x0),
 167 fh2JetsumHT3R2bbp(0x0),
 168 fh2JetsumHT3R4bb(0x0),
 169 fh2JetsumHT3R4bbp(0x0),
 170 fh2JetsumHT3R6bb(0x0),
 171 fh2JetsumHT3R6bbp(0x0),
 172 fh2JetsumHT3R8bb(0x0),
 173 fh2JetsumHT3R8bbp(0x0),
 174 fh2JetsumHT3R10bb(0x0),
 175 fh2JetsumHT3R10bbp(0x0),
 176 fh2JetsumHT3R2bbb(0x0),
 177 fh2JetsumHT3R2bbbp(0x0),
 178 fh2JetsumHT3R4bbb(0x0),
 179 fh2JetsumHT3R4bbbp(0x0),
 180 fh2JetsumHT3R6bbb(0x0),
 181 fh2JetsumHT3R6bbbp(0x0),
 182 fh2JetsumHT3R8bbb(0x0),
 183 fh2JetsumHT3R8bbbp(0x0),
 184 fh2JetsumHT3R10bbb(0x0),
 185 fh2JetsumHT3R10bbbp(0x0),
 186 fh2Ntriggers(0x0),
 187 fh2JetDensity(0x0),
 188 fh2JetDensityA4(0x0),
 189 fh3spectriggeredC10(0x0),
 190 fh3spectriggeredC20(0x0),
 191 fh3spectriggeredC3060(0x0),
 192 fh3specbiased(0x0),
 193 fh3spectot(0x0),
 194 fh3spectotb(0x0)
 195
 196 {
 197    // default Constructor
 198
 199
 200  // Trigger buffer.
 201    for(Int_t i=0; i<10; i++) {
 202                 for(Int_t j=0; j<7; j++) {
 203                         fTrigBuffer[i][j]=0;
 204                 }
 205         }
 206
 207
 208
 209
 210
 211    fJetBranchName[0] = "";
 212    fJetBranchName[1] = "";
 213
 214    fListJets[0] = new TList;
 215    fListJets[1] = new TList;
 216 }
 217
 218 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore(const char *name) :
 219 AliAnalysisTaskSE(name),
 220 fESD(0x0),
 221 fAODIn(0x0),
 222 fAODOut(0x0),
 223 fAODExtension(0x0),
 224 fBackgroundBranch(""),
 225 fNonStdFile(""),
 226 fIsPbPb(kTRUE),
 227 fOfflineTrgMask(AliVEvent::kAny),
 228 fMinContribVtx(1),
 229 fVtxZMin(-10.),
 230 fVtxZMax(10.),
 231 fEvtClassMin(0),
 232 fEvtClassMax(4),
 233 fFilterMask(0),
 234 fRadioFrac(0.2),
 235 fMinDist(0.1),
 236 fCentMin(0.),
 237 fCentMax(100.),
 238 fNInputTracksMin(0),
 239 fNInputTracksMax(-1),
 240 fAngStructCloseTracks(0),
 241 fCheckMethods(0),
 242 fDoEventMixing(0),
 243 fJetEtaMin(-.5),
 244 fJetEtaMax(.5),
 245 fNevents(0x0),
 246 fTindex(0x0),
 247 fTrigBufferIndex(0x0),
 248 fJetPtMin(20.),
 249 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
 250 fJetPtFractionMin(0.5),
 251 fNMatchJets(4),
 252 fMatchMaxDist(0.8),
 253 fKeepJets(kFALSE),
 254 fkNbranches(2),
 255 fkEvtClasses(12),
 256 fOutputList(0x0),
 257 fbEvent(kTRUE),
 258 fHistEvtSelection(0x0),
 259 fhnDeltaR(0x0),
 260 fhnMixedEvents(0x0),
 261 fh2JetCoreMethod1C10(0x0),
 262 fh2JetCoreMethod2C10(0x0),
 263 fh2JetCoreMethod1C20(0x0),
 264 fh2JetCoreMethod2C20(0x0),
 265 fh2JetCoreMethod1C30(0x0),
 266 fh2JetCoreMethod2C30(0x0),
 267 fh2JetCoreMethod1C60(0x0),
 268 fh2JetCoreMethod2C60(0x0),
 269 fh2JetTrackC10(0x0),
 270 fh2JetTrackC20(0x0),
 271 fh2AngStructpt1C10(0x0),
 272 fh2AngStructpt2C10(0x0),
 273 fh2AngStructpt3C10(0x0),
 274 fh2AngStructpt4C10(0x0),
 275 fh2AngStructpt1C20(0x0),
 276 fh2AngStructpt2C20(0x0),
 277 fh2AngStructpt3C20(0x0),
 278 fh2AngStructpt4C20(0x0),
 279 fh2AngStructpt1C30(0x0),
 280 fh2AngStructpt2C30(0x0),
 281 fh2AngStructpt3C30(0x0),
 282 fh2AngStructpt4C30(0x0),
 283 fh2AngStructpt1C60(0x0),
 284 fh2AngStructpt2C60(0x0),
 285 fh2AngStructpt3C60(0x0),
 286 fh2AngStructpt4C60(0x0),
 287 fh2JetsumHT3R2a(0x0),
 288 fh2JetsumHT3R2ap(0x0),
 289 fh2JetsumHT3R4a(0x0),
 290 fh2JetsumHT3R4ap(0x0),
 291 fh2JetsumHT3R6a(0x0),
 292 fh2JetsumHT3R6ap(0x0),
 293 fh2JetsumHT3R8a(0x0),
 294 fh2JetsumHT3R8ap(0x0),
 295 fh2JetsumHT3R10a(0x0),
 296 fh2JetsumHT3R10ap(0x0),
 297 fh2JetsumHT3R2aa(0x0),
 298 fh2JetsumHT3R2aap(0x0),
 299 fh2JetsumHT3R4aa(0x0),
 300 fh2JetsumHT3R4aap(0x0),
 301 fh2JetsumHT3R6aa(0x0),
 302 fh2JetsumHT3R6aap(0x0),
 303 fh2JetsumHT3R8aa(0x0),
 304 fh2JetsumHT3R8aap(0x0),
 305 fh2JetsumHT3R10aa(0x0),
 306 fh2JetsumHT3R10aap(0x0),
 307 fh2JetsumHT3R2aaa(0x0),
 308 fh2JetsumHT3R2aaap(0x0),
 309 fh2JetsumHT3R4aaa(0x0),
 310 fh2JetsumHT3R4aaap(0x0),
 311 fh2JetsumHT3R6aaa(0x0),
 312 fh2JetsumHT3R6aaap(0x0),
 313 fh2JetsumHT3R8aaa(0x0),
 314 fh2JetsumHT3R8aaap(0x0),
 315 fh2JetsumHT3R10aaa(0x0),
 316 fh2JetsumHT3R10aaap(0x0),
 317 fh2JetsumHT3R2b(0x0),
 318 fh2JetsumHT3R2bp(0x0),
 319 fh2JetsumHT3R4b(0x0),
 320 fh2JetsumHT3R4bp(0x0),
 321 fh2JetsumHT3R6b(0x0),
 322 fh2JetsumHT3R6bp(0x0),
 323 fh2JetsumHT3R8b(0x0),
 324 fh2JetsumHT3R8bp(0x0),
 325 fh2JetsumHT3R10b(0x0),
 326 fh2JetsumHT3R10bp(0x0),
 327 fh2JetsumHT3R2bb(0x0),
 328 fh2JetsumHT3R2bbp(0x0),
 329 fh2JetsumHT3R4bb(0x0),
 330 fh2JetsumHT3R4bbp(0x0),
 331 fh2JetsumHT3R6bb(0x0),
 332 fh2JetsumHT3R6bbp(0x0),
 333 fh2JetsumHT3R8bb(0x0),
 334 fh2JetsumHT3R8bbp(0x0),
 335 fh2JetsumHT3R10bb(0x0),
 336 fh2JetsumHT3R10bbp(0x0),
 337 fh2JetsumHT3R2bbb(0x0),
 338 fh2JetsumHT3R2bbbp(0x0),
 339 fh2JetsumHT3R4bbb(0x0),
 340 fh2JetsumHT3R4bbbp(0x0),
 341 fh2JetsumHT3R6bbb(0x0),
 342 fh2JetsumHT3R6bbbp(0x0),
 343 fh2JetsumHT3R8bbb(0x0),
 344 fh2JetsumHT3R8bbbp(0x0),
 345 fh2JetsumHT3R10bbb(0x0),
 346 fh2JetsumHT3R10bbbp(0x0),
 347 fh2Ntriggers(0x0),
 348 fh2JetDensity(0x0),
 349 fh2JetDensityA4(0x0),
 350 fh3spectriggeredC10(0x0),
 351 fh3spectriggeredC20(0x0),
 352 fh3spectriggeredC3060(0x0),
 353 fh3specbiased(0x0),
 354 fh3spectot(0x0),
 355 fh3spectotb(0x0)
 356  {
 357    // Constructor
 358
 359
 360     for(Int_t i=0; i<10; i++) {
 361        for(Int_t j=0; j<7; j++) {
 362             fTrigBuffer[i][j]=0;
 363                 }
 364     }
 365
 366
 367
 368    fJetBranchName[0] = "";
 369    fJetBranchName[1] = "";
 370
 371    fListJets[0] = new TList;
 372    fListJets[1] = new TList;
 373
 374    DefineOutput(1, TList::Class());
 375 }
 376
 377 AliAnalysisTaskJetCore::~AliAnalysisTaskJetCore()
 378 {
 379    delete fListJets[0];
 380    delete fListJets[1];
 381 }
 382
 383 void AliAnalysisTaskJetCore::SetBranchNames(const TString &branch1, const TString &branch2)
 384 {
 385    fJetBranchName[0] = branch1;
 386    fJetBranchName[1] = branch2;
 387 }
 388
 389 void AliAnalysisTaskJetCore::Init()
 390 {
 391
 392    // check for jet branches
 393    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 394       AliError("Jet branch name not set.");
 395    }
 396
 397 }
 398
 399 void AliAnalysisTaskJetCore::UserCreateOutputObjects()
 400 {
 401    // Create histograms
 402    // Called once
 403    OpenFile(1);
 404    if(!fOutputList) fOutputList = new TList;
 405    fOutputList->SetOwner(kTRUE);
 406
 407    Bool_t oldStatus = TH1::AddDirectoryStatus();
 408    TH1::AddDirectory(kFALSE);
 409
 410
 411    fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
 412    fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
 413    fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
 414    fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
 415    fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
 416    fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
 417    fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
 418
 419      UInt_t entries = 0; // bit coded, see GetDimParams() below
 420      entries = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 |1<<7;
 421      fhnDeltaR = NewTHnSparseF("fhnDeltaR", entries);
 422
 423      //change binning in pTtrack
 424      Double_t *xPt3 = new Double_t[10];
 425      xPt3[0] = 0.;
 426      for(int i = 1; i<=9;i++){
 427       if(xPt3[i-1]<2)xPt3[i] = xPt3[i-1] + 0.4; // 1 - 5
 428       else if(xPt3[i-1]<11)xPt3[i] = xPt3[i-1] + 3; // 5 - 12
 429       else xPt3[i] = xPt3[i-1] + 150.; // 18
 430      }
 431     fhnDeltaR->SetBinEdges(2,xPt3);
 432     delete [] xPt3;
 433
 434     //change binning in HTI
 435      Double_t *xPt4 = new Double_t[14];
 436      xPt4[0] = 0.;
 437      for(int i = 1; i<=13;i++){
 438       if(xPt4[i-1]<10)xPt4[i] = xPt4[i-1] + 1; // 1 - 10
 439       else if(xPt4[i-1]<20)xPt4[i] = xPt4[i-1] + 5; // 10 - 12
 440       else xPt4[i] = xPt4[i-1] + 30.; // 13
 441      }
 442     fhnDeltaR->SetBinEdges(6,xPt4);
 443     delete [] xPt4;
 444
 445
 446
 447
 448
 449
 450
 451      if(fDoEventMixing){
 452      UInt_t cifras = 0; // bit coded, see GetDimParams() below
 453      cifras = 1<<0 | 1<<1 | 1<<2 | 1<<3 | 1<<4 | 1<<5 | 1<<6 |1<<7;
 454      fhnMixedEvents = NewTHnSparseF("fhnMixedEvents", cifras);}
 455
 456     if(fCheckMethods){
 457
 458     fh2JetCoreMethod1C10 = new TH2F("JetCoreMethod1C10","",150, 0., 150.,100, 0., 1.5);
 459     fh2JetCoreMethod2C10 = new TH2F("JetCoreMethod2C10","",150, 0., 150.,100, 0., 1.5);
 460     fh2JetCoreMethod1C20 = new TH2F("JetCoreMethod1C20","",150, 0., 150.,100, 0., 1.5);
 461     fh2JetCoreMethod2C20 = new TH2F("JetCoreMethod2C20","",150, 0., 150.,100, 0., 1.5);
 462     fh2JetCoreMethod1C30 = new TH2F("JetCoreMethod1C30","",150, 0., 150.,100, 0., 1.5);
 463     fh2JetCoreMethod2C30 = new TH2F("JetCoreMethod2C30","",150, 0., 150.,100, 0., 1.5);
 464     fh2JetCoreMethod1C60 = new TH2F("JetCoreMethod1C60","",150, 0., 150.,100, 0., 1.5);
 465     fh2JetCoreMethod2C60 = new TH2F("JetCoreMethod2C60","",150, 0., 150.,100, 0., 1.5);}
 466
 467     fh2JetTrackC10=new TH2F("JetTrackC10","",150,0.,150.,32,0,3.2);
 468     fh2JetTrackC20=new TH2F("JetTrackC20","",150,0.,150.,32,0,3.2);
 469
 470     if(fAngStructCloseTracks>0){
 471     fh2AngStructpt1C10 = new TH2F("Ang struct pt1 C10","",15,0.,1.5,150,0.,10.);
 472     fh2AngStructpt2C10 = new TH2F("Ang struct pt2 C10","",15,0.,1.5,150,0.,10.);
 473     fh2AngStructpt3C10 = new TH2F("Ang struct pt3 C10","",15,0.,1.5,150,0.,10.);
 474     fh2AngStructpt4C10 = new TH2F("Ang struct pt4 C10","",15,0.,1.5,150,0.,10.);
 475     fh2AngStructpt1C20 = new TH2F("Ang struct pt1 C20","",15,0.,1.5,150,0.,10.);
 476     fh2AngStructpt2C20 = new TH2F("Ang struct pt2 C20","",15,0.,1.5,150,0.,10.);
 477     fh2AngStructpt3C20 = new TH2F("Ang struct pt3 C20","",15,0.,1.5,150,0.,10.);
 478     fh2AngStructpt4C20 = new TH2F("Ang struct pt4 C20","",15,0.,1.5,150,0.,10.);
 479     fh2AngStructpt1C30 = new TH2F("Ang struct pt1 C30","",15,0.,1.5,150,0.,10.);
 480     fh2AngStructpt2C30 = new TH2F("Ang struct pt2 C30","",15,0.,1.5,150,0.,10.);
 481     fh2AngStructpt3C30 = new TH2F("Ang struct pt3 C30","",15,0.,1.5,150,0.,10.);
 482     fh2AngStructpt4C30 = new TH2F("Ang struct pt4 C30","",15,0.,1.5,150,0.,10.);
 483     fh2AngStructpt1C60 = new TH2F("Ang struct pt1 C60","",15,0.,1.5,150,0.,10.);
 484     fh2AngStructpt2C60 = new TH2F("Ang struct pt2 C60","",15,0.,1.5,150,0.,10.);
 485     fh2AngStructpt3C60 = new TH2F("Ang struct pt3 C60","",15,0.,1.5,150,0.,10.);
 486     fh2AngStructpt4C60 = new TH2F("Ang struct pt4 C60","",15,0.,1.5,150,0.,10.); }
 487
 488     fh2JetsumHT3R2a = new TH2F("Pt sum R02 HT0 TT10","",20,0.,200.,100,0.,10.);
 489     fh2JetsumHT3R2ap = new TH2F("Pt sum R02 HT0 TT10 p","",20,0.,200.,100,0.,10.);
 490     fh2JetsumHT3R4a = new TH2F("Pt sum R04 HT0 TT10","",20,0.,200.,100,0.,10.);
 491     fh2JetsumHT3R4ap = new TH2F("Pt sum R04 HT0 TT10 p","",20,0.,200.,100,0.,10.);
 492     fh2JetsumHT3R6a = new TH2F("Pt sum R06 HT0 TT10","",20,0.,200.,100,0.,10.);
 493     fh2JetsumHT3R6ap = new TH2F("Pt sum R06 HT0 TT10 p","",20,0.,200.,100,0.,10.);
 494     fh2JetsumHT3R8a = new TH2F("Pt sum R08 HT0 TT10","",20,0.,200.,100,0.,10.);
 495     fh2JetsumHT3R8ap = new TH2F("Pt sum R08 HT0 TT10 p","",20,0.,200.,100,0.,10.);
 496     fh2JetsumHT3R10a = new TH2F("Pt sum R10 HT0 TT10","",20,0.,200.,100,0.,10.);
 497     fh2JetsumHT3R10ap = new TH2F("Pt sum R10 HT0 TT10 p","",20,0.,200.,100,0.,10.);
 498
 499     fh2JetsumHT3R2aa = new TH2F("Pt sum R02 HT0 TT20","",20,0.,200.,100,0.,10.);
 500     fh2JetsumHT3R2aap = new TH2F("Pt sum R02 HT0 TT20 p","",20,0.,200.,100,0.,10.);
 501     fh2JetsumHT3R4aa = new TH2F("Pt sum R04 HT0 TT20","",20,0.,200.,100,0.,10.);
 502     fh2JetsumHT3R4aap = new TH2F("Pt sum R04 HT0 TT20 p","",20,0.,200.,100,0.,10.);
 503     fh2JetsumHT3R6aa = new TH2F("Pt sum R06 HT0 TT20","",20,0.,200.,100,0.,10.);
 504     fh2JetsumHT3R6aap = new TH2F("Pt sum R06 HT0 TT20 p","",20,0.,200.,100,0.,10.);
 505     fh2JetsumHT3R8aa = new TH2F("Pt sum R08 HT0 TT20","",20,0.,200.,100,0.,10.);
 506     fh2JetsumHT3R8aap = new TH2F("Pt sum R08 HT0 TT20 p","",20,0.,200.,100,0.,10.);
 507     fh2JetsumHT3R10aa = new TH2F("Pt sum R10 HT0 TT20","",20,0.,200.,100,0.,10.);
 508     fh2JetsumHT3R10aap = new TH2F("Pt sum R10 HT0 TT20 p","",20,0.,200.,100,0.,10.);
 509
 510     fh2JetsumHT3R2aaa = new TH2F("Pt sum R02 HT0 TT0","",20,0.,200.,100,0.,10.);
 511     fh2JetsumHT3R2aaap = new TH2F("Pt sum R02 HT0 TT0 p","",20,0.,200.,100,0.,10.);
 512     fh2JetsumHT3R4aaa = new TH2F("Pt sum R04 HT0 TT0","",20,0.,200.,100,0.,10.);
 513     fh2JetsumHT3R4aaap = new TH2F("Pt sum R04 HT0 TT0 p","",20,0.,200.,100,0.,10.);
 514     fh2JetsumHT3R6aaa = new TH2F("Pt sum R06 HT0 TT0","",20,0.,200.,100,0.,10.);
 515     fh2JetsumHT3R6aaap = new TH2F("Pt sum R06 HT0 TT0 p","",20,0.,200.,100,0.,10.);
 516     fh2JetsumHT3R8aaa = new TH2F("Pt sum R08 HT0 TT0","",20,0.,200.,100,0.,10.);
 517     fh2JetsumHT3R8aaap = new TH2F("Pt sum R08 HT0 TT0 p","",20,0.,200.,100,0.,10.);
 518     fh2JetsumHT3R10aaa = new TH2F("Pt sum R10 HT0 TT0","",20,0.,200.,100,0.,10.);
 519     fh2JetsumHT3R10aaap = new TH2F("Pt sum R10 HT0 TT0 p","",20,0.,200.,100,0.,10.);
 520
 521     fh2JetsumHT3R2b = new TH2F("Pt sum R02 HT6 TT10","",20,0.,200.,100,0.,10.);
 522     fh2JetsumHT3R2bp = new TH2F("Pt sum R02 HT6 TT10 p","",20,0.,200.,100,0.,10.);
 523     fh2JetsumHT3R4b = new TH2F("Pt sum R04 HT6 TT10","",20,0.,200.,100,0.,10.);
 524     fh2JetsumHT3R4bp = new TH2F("Pt sum R04 HT6 TT10 p","",20,0.,200.,100,0.,10.);
 525     fh2JetsumHT3R6b = new TH2F("Pt sum R06 HT6 TT10","",20,0.,200.,100,0.,10.);
 526     fh2JetsumHT3R6bp = new TH2F("Pt sum R06 HT6 TT10 p","",20,0.,200.,100,0.,10.);
 527     fh2JetsumHT3R8b = new TH2F("Pt sum R08 HT6 TT10","",20,0.,200.,100,0.,10.);
 528     fh2JetsumHT3R8bp = new TH2F("Pt sum R08 HT6 TT10 p","",20,0.,200.,100,0.,10.);
 529     fh2JetsumHT3R10b = new TH2F("Pt sum R10 HT6 TT10","",20,0.,200.,100,0.,10.);
 530     fh2JetsumHT3R10bp = new TH2F("Pt sum R10 HT6 TT10 p","",20,0.,200.,100,0.,10.);
 531
 532     fh2JetsumHT3R2bb = new TH2F("Pt sum R02 HT6 TT20","",20,0.,200.,100,0.,10.);
 533     fh2JetsumHT3R2bbp = new TH2F("Pt sum R02 HT6 TT20 p","",20,0.,200.,100,0.,10.);
 534     fh2JetsumHT3R4bb = new TH2F("Pt sum R04 HT6 TT20","",20,0.,200.,100,0.,10.);
 535     fh2JetsumHT3R4bbp = new TH2F("Pt sum R04 HT6 TT20 p","",20,0.,200.,100,0.,10.);
 536     fh2JetsumHT3R6bb = new TH2F("Pt sum R06 HT6 TT20","",20,0.,200.,100,0.,10.);
 537     fh2JetsumHT3R6bbp = new TH2F("Pt sum R06 HT6 TT20 p","",20,0.,200.,100,0.,10.);
 538     fh2JetsumHT3R8bb = new TH2F("Pt sum R08 HT6 TT20","",20,0.,200.,100,0.,10.);
 539     fh2JetsumHT3R8bbp = new TH2F("Pt sum R08 HT6 TT20 p","",20,0.,200.,100,0.,10.);
 540     fh2JetsumHT3R10bb = new TH2F("Pt sum R10 HT6 TT20","",20,0.,200.,100,0.,10.);
 541     fh2JetsumHT3R10bbp = new TH2F("Pt sum R10 HT6 TT20 p","",20,0.,200.,100,0.,10.);
 542
 543     fh2JetsumHT3R2bbb = new TH2F("Pt sum R02 HT6 TT0","",20,0.,200.,100,0.,10.);
 544     fh2JetsumHT3R2bbbp = new TH2F("Pt sum R02 HT6 TT0 p","",20,0.,200.,100,0.,10.);
 545     fh2JetsumHT3R4bbb = new TH2F("Pt sum R04 HT6 TT0","",20,0.,200.,100,0.,10.);
 546     fh2JetsumHT3R4bbbp = new TH2F("Pt sum R04 HT6 TT0 p","",20,0.,200.,100,0.,10.);
 547     fh2JetsumHT3R6bbb = new TH2F("Pt sum R06 HT6 TT0","",20,0.,200.,100,0.,10.);
 548     fh2JetsumHT3R6bbbp = new TH2F("Pt sum R06 HT6 TT0 p","",20,0.,200.,100,0.,10.);
 549     fh2JetsumHT3R8bbb = new TH2F("Pt sum R08 HT6 TT0","",20,0.,200.,100,0.,10.);
 550     fh2JetsumHT3R8bbbp = new TH2F("Pt sum R08 HT6 TT0 p","",20,0.,200.,100,0.,10.);
 551     fh2JetsumHT3R10bbb = new TH2F("Pt sum R10 HT6 TT0","",20,0.,200.,100,0.,10.);
 552     fh2JetsumHT3R10bbbp = new TH2F("Pt sum R10 HT6 TT0 p","",20,0.,200.,100,0.,10.);
 553
 554
 555
 556     fh2Ntriggers=new TH2F("# of triggers","",10,0.,100.,50,0.,50.);
 557     fh2JetDensity=new TH2F("Jet density vs centrality A>0.4","",100,0.,4000.,100,0.,5.);
 558     fh2JetDensityA4=new TH2F("Jet density vs multiplicity A>0.4","",100,0.,4000.,100,0.,5.);
 559     fh3spectriggeredC10 = new TH3F("Triggered spectrumC10","",5,0.,1.,140,-80.,200.,50,0.,50.);
 560     fh3spectriggeredC20 = new TH3F("Triggered spectrumC20","",5,0.,1.,140,-80.,200.,50,0.,50.);
 561     fh3spectriggeredC3060 = new TH3F("Triggered spectrumC3060","",5,0.,1.,140,-80.,200.,50,0.,50.);
 562
 563     fh3specbiased = new TH3F("Biased spectrum","",10,0,100,50,0.,200.,50,0.,50.);
 564     fh3spectot = new TH3F("Total spectrum 0-10","",50,0.,200.,50,0.,50.,50,0.,50.);
 565     fh3spectotb = new TH3F("Total spectrum 30-60","",50,0.,200.,50,0.,50.,50,0.,50.);
 566    fOutputList->Add(fHistEvtSelection);
 567
 568    fOutputList->Add(fhnDeltaR);
 569
 570    fOutputList->Add(fhnMixedEvents);
 571
 572
 573
 574       if(fCheckMethods){
 575
 576       fOutputList->Add(fh2JetCoreMethod1C10);
 577       fOutputList->Add(fh2JetCoreMethod2C10);
 578       fOutputList->Add(fh2JetCoreMethod1C20);
 579       fOutputList->Add(fh2JetCoreMethod2C20);
 580       fOutputList->Add(fh2JetCoreMethod1C30);
 581       fOutputList->Add(fh2JetCoreMethod2C30);
 582       fOutputList->Add(fh2JetCoreMethod1C60);
 583       fOutputList->Add(fh2JetCoreMethod2C60);}
 584
 585       fOutputList->Add(fh2JetTrackC10);
 586       fOutputList->Add(fh2JetTrackC20);
 587
 588
 589
 590
 591         if(fAngStructCloseTracks>0){
 592        fOutputList->Add(fh2AngStructpt1C10);
 593        fOutputList->Add(fh2AngStructpt2C10);
 594        fOutputList->Add(fh2AngStructpt3C10);
 595        fOutputList->Add(fh2AngStructpt4C10);
 596        fOutputList->Add(fh2AngStructpt1C20);
 597        fOutputList->Add(fh2AngStructpt2C20);
 598        fOutputList->Add(fh2AngStructpt3C20);
 599        fOutputList->Add(fh2AngStructpt4C20);
 600        fOutputList->Add(fh2AngStructpt1C30);
 601        fOutputList->Add(fh2AngStructpt2C30);
 602        fOutputList->Add(fh2AngStructpt3C30);
 603        fOutputList->Add(fh2AngStructpt4C30);
 604        fOutputList->Add(fh2AngStructpt1C60);
 605        fOutputList->Add(fh2AngStructpt2C60);
 606        fOutputList->Add(fh2AngStructpt3C60);
 607        fOutputList->Add(fh2AngStructpt4C60);}
 608
 609
 610         fOutputList->Add(fh2JetsumHT3R2a);
 611         fOutputList->Add(fh2JetsumHT3R2ap);
 612         fOutputList->Add(fh2JetsumHT3R4a);
 613         fOutputList->Add(fh2JetsumHT3R4ap);
 614         fOutputList->Add(fh2JetsumHT3R6a);
 615         fOutputList->Add(fh2JetsumHT3R6ap);
 616         fOutputList->Add(fh2JetsumHT3R8a);
 617         fOutputList->Add(fh2JetsumHT3R8ap);
 618         fOutputList->Add(fh2JetsumHT3R10a);
 619         fOutputList->Add(fh2JetsumHT3R10ap);
 620         fOutputList->Add(fh2JetsumHT3R2aa);
 621         fOutputList->Add(fh2JetsumHT3R2aap);
 622         fOutputList->Add(fh2JetsumHT3R4aa);
 623         fOutputList->Add(fh2JetsumHT3R4aap);
 624         fOutputList->Add(fh2JetsumHT3R6aa);
 625         fOutputList->Add(fh2JetsumHT3R6aap);
 626         fOutputList->Add(fh2JetsumHT3R8aa);
 627         fOutputList->Add(fh2JetsumHT3R8aap);
 628         fOutputList->Add(fh2JetsumHT3R10aa);
 629         fOutputList->Add(fh2JetsumHT3R10aap);
 630         fOutputList->Add(fh2JetsumHT3R2aaa);
 631         fOutputList->Add(fh2JetsumHT3R2aaap);
 632         fOutputList->Add(fh2JetsumHT3R4aaa);
 633         fOutputList->Add(fh2JetsumHT3R4aaap);
 634         fOutputList->Add(fh2JetsumHT3R6aaa);
 635         fOutputList->Add(fh2JetsumHT3R6aaap);
 636         fOutputList->Add(fh2JetsumHT3R8aaa);
 637         fOutputList->Add(fh2JetsumHT3R8aaap);
 638         fOutputList->Add(fh2JetsumHT3R10aaa);
 639         fOutputList->Add(fh2JetsumHT3R10aaap);
 640
 641         fOutputList->Add(fh2JetsumHT3R2b);
 642         fOutputList->Add(fh2JetsumHT3R2bp);
 643         fOutputList->Add(fh2JetsumHT3R4b);
 644         fOutputList->Add(fh2JetsumHT3R4bp);
 645         fOutputList->Add(fh2JetsumHT3R6b);
 646         fOutputList->Add(fh2JetsumHT3R6bp);
 647         fOutputList->Add(fh2JetsumHT3R8b);
 648         fOutputList->Add(fh2JetsumHT3R8bp);
 649         fOutputList->Add(fh2JetsumHT3R10b);
 650         fOutputList->Add(fh2JetsumHT3R10bp);
 651         fOutputList->Add(fh2JetsumHT3R2bb);
 652         fOutputList->Add(fh2JetsumHT3R2bbp);
 653         fOutputList->Add(fh2JetsumHT3R4bb);
 654         fOutputList->Add(fh2JetsumHT3R4bbp);
 655         fOutputList->Add(fh2JetsumHT3R6bb);
 656         fOutputList->Add(fh2JetsumHT3R6bbp);
 657         fOutputList->Add(fh2JetsumHT3R8bb);
 658         fOutputList->Add(fh2JetsumHT3R8bbp);
 659         fOutputList->Add(fh2JetsumHT3R10bb);
 660         fOutputList->Add(fh2JetsumHT3R10bbp);
 661         fOutputList->Add(fh2JetsumHT3R2bbb);
 662         fOutputList->Add(fh2JetsumHT3R2bbbp);
 663         fOutputList->Add(fh2JetsumHT3R4bbb);
 664         fOutputList->Add(fh2JetsumHT3R4bbbp);
 665         fOutputList->Add(fh2JetsumHT3R6bbb);
 666         fOutputList->Add(fh2JetsumHT3R6bbbp);
 667         fOutputList->Add(fh2JetsumHT3R8bbb);
 668         fOutputList->Add(fh2JetsumHT3R8bbbp);
 669         fOutputList->Add(fh2JetsumHT3R10bbb);
 670         fOutputList->Add(fh2JetsumHT3R10bbbp);
 671
 672
 673         fOutputList->Add(fh2Ntriggers);
 674         fOutputList->Add(fh2JetDensity);
 675         fOutputList->Add(fh2JetDensityA4);
 676
 677        fOutputList->Add(fh3spectriggeredC10);
 678        fOutputList->Add(fh3spectriggeredC20);
 679        fOutputList->Add(fh3spectriggeredC3060);
 680
 681        fOutputList->Add(fh3specbiased);
 682        fOutputList->Add(fh3spectot);
 683        fOutputList->Add(fh3spectotb);
 684    // =========== Switch on Sumw2 for all histos ===========
 685    for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
 686       TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
 687       if (h1){
 688          h1->Sumw2();
 689          continue;
 690       }
 691     THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
 692       if (hn){
 693          hn->Sumw2();
 694       }
 695    }
 696    TH1::AddDirectory(oldStatus);
 697
 698    PostData(1, fOutputList);
 699 }
 700
 701 void AliAnalysisTaskJetCore::UserExec(Option_t *)
 702 {
 703
 704
 705    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 706       AliError("Jet branch name not set.");
 707       return;
 708    }
 709
 710    fESD=dynamic_cast<AliESDEvent*>(InputEvent());
 711    if (!fESD) {
 712       AliError("ESD not available");
 713       fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 714    }
 715       fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
 716
 717        static AliAODEvent* aod = 0;
 718        // take all other information from the aod we take the tracks from
 719        if(!aod){
 720        if(!fESD)aod = fAODIn;
 721        else aod = fAODOut;}
 722
 723
 724
 725     if(fNonStdFile.Length()!=0){
 726     // case that we have an AOD extension we need can fetch the jets from the extended output
 727     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 728     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 729     if(!fAODExtension){
 730       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
 731     }}
 732
 733
 734
 735
 736
 737    // -- event selection --
 738    fHistEvtSelection->Fill(1); // number of events before event selection
 739
 740    // physics selection
 741    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
 742    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 743    cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<endl;
 744    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
 745       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
 746       fHistEvtSelection->Fill(2);
 747       PostData(1, fOutputList);
 748       return;
 749    }
 750
 751    // vertex selection
 752    if(!aod){
 753      if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
 754      fHistEvtSelection->Fill(3);
 755       PostData(1, fOutputList);
 756    }
 757    AliAODVertex* primVtx = aod->GetPrimaryVertex();
 758
 759    if(!primVtx){
 760      if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
 761      fHistEvtSelection->Fill(3);
 762      PostData(1, fOutputList);
 763      return;
 764    }
 765
 766    Int_t nTracksPrim = primVtx->GetNContributors();
 767    if ((nTracksPrim < fMinContribVtx) ||
 768          (primVtx->GetZ() < fVtxZMin) ||
 769          (primVtx->GetZ() > fVtxZMax) ){
 770       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
 771       fHistEvtSelection->Fill(3);
 772       PostData(1, fOutputList);
 773       return;
 774    }
 775
 776    // event class selection (from jet helper task)
 777    Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
 778    if(fDebug) Printf("Event class %d", eventClass);
 779    if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
 780       fHistEvtSelection->Fill(4);
 781       PostData(1, fOutputList);
 782       return;
 783    }
 784
 785    // centrality selection
 786    AliCentrality *cent = 0x0;
 787    Double_t centValue = 0.;
 788    if(fIsPbPb){
 789    if(fESD) {cent = fESD->GetCentrality();
 790      if(cent) centValue = cent->GetCentralityPercentile("V0M");}
 791    else     centValue=aod->GetHeader()->GetCentrality();
 792
 793    if(fDebug) printf("centrality: %f\n", centValue);
 794       if (centValue < fCentMin || centValue > fCentMax){
 795       fHistEvtSelection->Fill(4);
 796       PostData(1, fOutputList);
 797       return;
 798       }}
 799
 800
 801    fHistEvtSelection->Fill(0);
 802    // accepted events
 803    // -- end event selection --
 804
 805    // get background
 806    AliAODJetEventBackground* externalBackground = 0;
 807    if(fAODOut&&!externalBackground&&fBackgroundBranch.Length()){
 808       externalBackground =  (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
 809       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 810    }
 811    if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
 812      externalBackground =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
 813       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 814    }
 815
 816     if(fAODIn&&!externalBackground&&fBackgroundBranch.Length()){
 817       externalBackground =  (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data()));
 818       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 819     }
 820
 821    Float_t rho = 0;
 822    if(externalBackground)rho = externalBackground->GetBackground(0);
 823
 824
 825    // fetch jets
 826    TClonesArray *aodJets[2];
 827    aodJets[0]=0;
 828    if(fAODOut&&!aodJets[0]){
 829    aodJets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data()));
 830    aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data()));  }
 831    if(fAODExtension && !aodJets[0]){
 832    aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
 833    aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));  }
 834      if(fAODIn&&!aodJets[0]){
 835    aodJets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data()));
 836    aodJets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data()));  }
 837
 838
 839    //Double_t ptsub[aodJets[0]->GetEntriesFast()];
 840    //Int_t inord[aodJets[0]->GetEntriesFast()];
 841    //for(Int_t n=0;n<aodJets[0]->GetEntriesFast();n++){
 842    //  ptsub[n]=0;
 843    //  inord[n]=0;}
 844
 845    TList ParticleList;
 846    Int_t nT = GetListOfTracks(&ParticleList);
 847      for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
 848       fListJets[iJetType]->Clear();
 849       if (!aodJets[iJetType]) continue;
 850
 851       if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
 852
 853
 854       for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
 855          AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
 856          if (jet) fListJets[iJetType]->Add(jet);
 857          // if(iJetType==0){
 858          // ptsub[iJet]=jet->Pt()-rho*jet->EffectiveAreaCharged();}
 859       }}
 860
 861    Double_t etabig=0;
 862    Double_t ptbig=0;
 863    Double_t areabig=0;
 864    Double_t phibig=0.;
 865    Double_t etasmall=0;
 866    Double_t ptsmall=0;
 867    Double_t areasmall=0;
 868    Double_t phismall=0.;
 869
 870
 871    // Double_t up1[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 872    // Double_t up2[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 873    // Double_t up3[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 874    // Double_t up4[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 875    // Double_t down1[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 876    // Double_t down2[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 877    // Double_t down3[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 878    // Double_t down4[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 879    Int_t iCount=0;
 880    Int_t trigJet=-1;
 881    Int_t trigBBTrack=-1;
 882    Int_t trigInTrack=-1;
 883
 884
 885    AliVParticle *partback = (AliVParticle*)ParticleList.At(nT);
 886    if(!partback){
 887    PostData(1, fOutputList);
 888    return;}
 889    fh2Ntriggers->Fill(centValue,partback->Pt());
 890    Double_t accep=2.*TMath::Pi()*1.8;
 891    Int_t injet4=0;
 892    Int_t injet=0;
 893    for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
 894            AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
 895            etabig  = jetbig->Eta();
 896            phibig  = jetbig->Phi();
 897            ptbig   = jetbig->Pt();
 898            if(ptbig==0) continue;
 899            areabig = jetbig->EffectiveAreaCharged();
 900            Double_t ptcorr=ptbig-rho*areabig;
 901            if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
 902
 903            if(areabig>=0.2) injet=injet+1;
 904            if(areabig>=0.4) injet4=injet4+1;
 905            Double_t dphi=RelativePhi(partback->Phi(),phibig);
 906            if(centValue<10.) fh2JetTrackC10->Fill(ptcorr,TMath::Abs(dphi));
 907            if(centValue<20.) fh2JetTrackC20->Fill(ptcorr,TMath::Abs(dphi));
 908
 909            if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;
 910                    Double_t dismin=100.;
 911                    Double_t ptmax=-10.;
 912                    Int_t index1=-1;
 913                    Int_t index2=-1;
 914            if(centValue<10.)  fh3spectriggeredC10->Fill(jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt());
 915            if(centValue<20.)  fh3spectriggeredC20->Fill(jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt());
 916            if(centValue>30. && centValue<60.)  fh3spectriggeredC3060->Fill(jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt());
 917
 918                    if(ptcorr<=0) continue;
 919                        AliAODTrack* leadtrack;
 920                        Int_t ippt=0;
 921                        Double_t ppt=-10;
 922                        TRefArray *genTrackList = jetbig->GetRefTracks();
 923                        Int_t nTracksGenJet = genTrackList->GetEntriesFast();
 924                        AliAODTrack* genTrack;
 925                        for(Int_t ir=0; ir<nTracksGenJet; ++ir){
 926                        genTrack = (AliAODTrack*)(genTrackList->At(ir));
 927                        if(genTrack->Pt()>ppt){ppt=genTrack->Pt();
 928                          ippt=ir;}}
 929                         leadtrack=(AliAODTrack*)(genTrackList->At(ippt));
 930                         if(!leadtrack) continue;
 931                         fh3specbiased->Fill(centValue,ptcorr,leadtrack->Pt());
 932                         if(centValue<10)fh3spectot->Fill(ptcorr,leadtrack->Pt(),partback->Pt());
 933                         if(centValue>30. && centValue<60.)fh3spectotb->Fill(ptcorr,leadtrack->Pt(),partback->Pt());
 934                         //store one trigger info
 935                         if((partback->Pt()>10.)&&(iCount==0)){
 936                           trigJet=i;
 937                           trigBBTrack=nT;
 938                           trigInTrack=ippt;
 939                           iCount=iCount+1;}
 940
 941
 942                   if(fCheckMethods){
 943                   for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
 944                   AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
 945                   etasmall  = jetsmall->Eta();
 946                   phismall = jetsmall->Phi();
 947                   ptsmall   = jetsmall->Pt();
 948                   areasmall = jetsmall->EffectiveAreaCharged();
 949                   Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
 950                   tmpDeltaR=TMath::Sqrt(tmpDeltaR);
 951                      //Fraction in the jet core
 952                     if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;
 953                     index2=j;}
 954                     if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
 955                       index1=j;}} //en of loop over R=0.2 jets
 956                 //method1:most concentric jet=core
 957                 if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));
 958                   if(centValue<10) fh2JetCoreMethod1C10->Fill(ptcorr,jetmethod1->Pt()/ptbig);
 959                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod1C20->Fill(ptcorr,jetmethod1->Pt()/ptbig);
 960                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod1C30->Fill(ptcorr,jetmethod1->Pt()/ptbig);
 961                   if(centValue>60) fh2JetCoreMethod1C60->Fill(ptcorr,jetmethod1->Pt()/ptbig); }
 962                 //method2:hardest contained jet=core
 963                 if(index2!=-1){
 964                   AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
 965                   if(centValue<10) fh2JetCoreMethod2C10->Fill(ptcorr,jetmethod2->Pt()/ptbig);
 966                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod2C20->Fill(ptcorr,jetmethod2->Pt()/ptbig);
 967                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptcorr,jetmethod2->Pt()/ptbig);
 968                   if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }}
 969                   Double_t sumpt2a=0.;
 970                   Double_t sumpt2b=0.;
 971                   Double_t sumpt4a=0.;
 972                   Double_t sumpt4b=0.;
 973                   Double_t sumpt6a=0.;
 974                   Double_t sumpt6b=0.;
 975                   Double_t sumpt8a=0.;
 976                   Double_t sumpt8b=0.;
 977                   Double_t sumpt10a=0.;
 978                   Double_t sumpt10b=0.;
 979                   Double_t sumpt2aaa=0.;
 980                   Double_t sumpt2bbb=0.;
 981                   Double_t sumpt4aaa=0.;
 982                   Double_t sumpt4bbb=0.;
 983                   Double_t sumpt6aaa=0.;
 984                   Double_t sumpt6bbb=0.;
 985                   Double_t sumpt8aaa=0.;
 986                   Double_t sumpt8bbb=0.;
 987                   Double_t sumpt10aaa=0.;
 988                   Double_t sumpt10bbb=0.;
 989                   Double_t sumpt2ap=0.;
 990                   Double_t sumpt2bp=0.;
 991                   Double_t sumpt4ap=0.;
 992                   Double_t sumpt4bp=0.;
 993                   Double_t sumpt6ap=0.;
 994                   Double_t sumpt6bp=0.;
 995                   Double_t sumpt8ap=0.;
 996                   Double_t sumpt8bp=0.;
 997                   Double_t sumpt10ap=0.;
 998                   Double_t sumpt10bp=0.;
 999                   Double_t sumpt2aaap=0.;
1000                   Double_t sumpt2bbbp=0.;
1001                   Double_t sumpt4aaap=0.;
1002                   Double_t sumpt4bbbp=0.;
1003                   Double_t sumpt6aaap=0.;
1004                   Double_t sumpt6bbbp=0.;
1005                   Double_t sumpt8aaap=0.;
1006                   Double_t sumpt8bbbp=0.;
1007                   Double_t sumpt10aaap=0.;
1008                   Double_t sumpt10bbbp=0.;
1009
1010
1011
1012
1013
1014                   for(int it = 0;it<ParticleList.GetEntries();++it){
1015           AliVParticle *part = (AliVParticle*)ParticleList.At(it);
1016           Double_t deltaR = jetbig->DeltaR(part);
1017           Double_t deltaEta = etabig-part->Eta();
1018           if(partback->Pt()>14.){
1019
1020           if(centValue<10){
1021             //for one more centrality
1022                if(deltaR<0.2){if(leadtrack->Pt()>0.)sumpt2a=sumpt2a+part->Pt();
1023                               if(leadtrack->Pt()>6) sumpt2b=sumpt2b+part->Pt();}
1024                if(deltaR>=0.2 && deltaR<0.4){if(leadtrack->Pt()>0.)sumpt4a=sumpt4a+part->Pt();
1025                                                if(leadtrack->Pt()>6.)sumpt4b=sumpt4b+part->Pt();}
1026
1027                if(deltaR>=0.4 && deltaR<0.6){if(leadtrack->Pt()>0.)sumpt6a=sumpt6a+part->Pt();
1028                                              if(leadtrack->Pt()>6.) sumpt6b=sumpt6b+part->Pt();}
1029                if(deltaR>=0.6 && deltaR<0.8){if(leadtrack->Pt()>0.)sumpt8a=sumpt8a+part->Pt();
1030                                              if(leadtrack->Pt()>6.) sumpt8b=sumpt8b+part->Pt();}
1031                if(deltaR>=0.8 && deltaR<1.2){if(leadtrack->Pt()>0.)sumpt10a=sumpt10a+part->Pt();
1032                  if(leadtrack->Pt()>6.) sumpt10b=sumpt10b+part->Pt();}}
1033
1034                if(centValue>30. && centValue<60.){
1035             //for one more centrality
1036                if(deltaR<0.2){if(leadtrack->Pt()>0.)sumpt2ap=sumpt2ap+part->Pt();
1037                               if(leadtrack->Pt()>6) sumpt2bp=sumpt2bp+part->Pt();}
1038                if(deltaR>=0.2 && deltaR<0.4){if(leadtrack->Pt()>0.)sumpt4ap=sumpt4ap+part->Pt();
1039                                                if(leadtrack->Pt()>6.)sumpt4bp=sumpt4bp+part->Pt();}
1040
1041                if(deltaR>=0.4 && deltaR<0.6){if(leadtrack->Pt()>0.)sumpt6ap=sumpt6ap+part->Pt();
1042                                              if(leadtrack->Pt()>6.) sumpt6bp=sumpt6bp+part->Pt();}
1043                if(deltaR>=0.6 && deltaR<0.8){if(leadtrack->Pt()>0.)sumpt8ap=sumpt8ap+part->Pt();
1044                                              if(leadtrack->Pt()>6.) sumpt8bp=sumpt8bp+part->Pt();}
1045                if(deltaR>=0.8 && deltaR<1.2){if(leadtrack->Pt()>0.)sumpt10ap=sumpt10ap+part->Pt();
1046                  if(leadtrack->Pt()>6.) sumpt10bp=sumpt10bp+part->Pt();}}}
1047
1048
1049
1050
1051              if(partback->Pt()<1.){
1052
1053               if(centValue<10){
1054             //for one more centrality
1055                if(deltaR<0.2){if(leadtrack->Pt()>0.)sumpt2aaa=sumpt2aaa+part->Pt();
1056                               if(leadtrack->Pt()>6.) sumpt2bbb=sumpt2bbb+part->Pt();}
1057                if(deltaR>=0.2 && deltaR<0.4){if(leadtrack->Pt()>0.)sumpt4aaa=sumpt4aaa+part->Pt();
1058                                                if(leadtrack->Pt()>6.)sumpt4bbb=sumpt4bbb+part->Pt();}
1059
1060                if(deltaR>=0.4 && deltaR<0.6){if(leadtrack->Pt()>0.)sumpt6aaa=sumpt6aaa+part->Pt();
1061                                              if(leadtrack->Pt()>6.) sumpt6bbb=sumpt6bbb+part->Pt();}
1062                if(deltaR>=0.6 && deltaR<0.8){if(leadtrack->Pt()>0.)sumpt8aaa=sumpt8aaa+part->Pt();
1063                                              if(leadtrack->Pt()>6.) sumpt8bbb=sumpt8bbb+part->Pt();}
1064                if(deltaR>=0.8 && deltaR<1.2){if(leadtrack->Pt()>0.)sumpt10aaa=sumpt10aaa+part->Pt();
1065                  if(leadtrack->Pt()>6.) sumpt10bbb=sumpt10bbb+part->Pt();}}
1066
1067                if(centValue>30. && centValue<60.){
1068             //for one more centrality
1069                if(deltaR<0.2){if(leadtrack->Pt()>0.)sumpt2aaap=sumpt2aaap+part->Pt();
1070                               if(leadtrack->Pt()>6) sumpt2bbbp=sumpt2bbbp+part->Pt();}
1071                if(deltaR>=0.2 && deltaR<0.4){if(leadtrack->Pt()>0.)sumpt4aaap=sumpt4aaap+part->Pt();
1072                                                if(leadtrack->Pt()>6.)sumpt4bbbp=sumpt4bbbp+part->Pt();}
1073
1074                if(deltaR>=0.4 && deltaR<0.6){if(leadtrack->Pt()>0.)sumpt6aaap=sumpt6aaap+part->Pt();
1075                                              if(leadtrack->Pt()>6.) sumpt6bbbp=sumpt6bbbp+part->Pt();}
1076                if(deltaR>=0.6 && deltaR<0.8){if(leadtrack->Pt()>0.)sumpt8aaap=sumpt8aaap+part->Pt();
1077                                              if(leadtrack->Pt()>6.) sumpt8bbbp=sumpt8bbbp+part->Pt();}
1078                if(deltaR>=0.8 && deltaR<1.2){if(leadtrack->Pt()>0.)sumpt10aaap=sumpt10aaap+part->Pt();
1079                  if(leadtrack->Pt()>6.) sumpt10bbbp=sumpt10bbbp+part->Pt();}}}
1080
1081
1082
1083
1084
1085
1086           Double_t deltaPhi=phibig-part->Phi();
1087           if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
1088           if(deltaPhi>3./2.*TMath::Pi()) deltaPhi-=2.*TMath::Pi();
1089           Double_t jetEntries[8] = {centValue,ptcorr,part->Pt(),deltaR,deltaEta,deltaPhi,leadtrack->Pt(),partback->Pt()};                     fhnDeltaR->Fill(jetEntries);
1090           }
1091           //end of track loop
1092             Double_t rhoin2=rho*TMath::Pi()*0.2*0.2;
1093             Double_t rhoin4=rho*TMath::Pi()*(0.4*0.4-0.2*0.2);
1094             Double_t rhoin6=rho*TMath::Pi()*(0.6*0.6-0.4*0.4);
1095             Double_t rhoin8=rho*TMath::Pi()*(0.8*0.8-0.6*0.6);
1096             Double_t rhoin10=rho*TMath::Pi()*(1.2*1.2-0.8*0.8);
1097
1098
1099           if(rho!=0){
1100
1101           if(partback->Pt()>14.){
1102             if(centValue<10.){
1103             if(leadtrack->Pt()>0.){
1104           fh2JetsumHT3R2a->Fill(ptcorr,sumpt2a/rhoin2);
1105           fh2JetsumHT3R4a->Fill(ptcorr,sumpt4a/rhoin4);
1106           fh2JetsumHT3R6a->Fill(ptcorr,sumpt6a/rhoin6);
1107           fh2JetsumHT3R8a->Fill(ptcorr,sumpt8a/rhoin8);
1108           fh2JetsumHT3R10a->Fill(ptcorr,sumpt10a/rhoin10);}
1109             if(leadtrack->Pt()>6.){
1110           fh2JetsumHT3R2b->Fill(ptcorr,sumpt2b/rhoin2);
1111           fh2JetsumHT3R4b->Fill(ptcorr,sumpt4b/rhoin4);
1112           fh2JetsumHT3R6b->Fill(ptcorr,sumpt6b/rhoin6);
1113           fh2JetsumHT3R8b->Fill(ptcorr,sumpt8b/rhoin8);
1114           fh2JetsumHT3R10b->Fill(ptcorr,sumpt10b/rhoin10);}}
1115
1116                     if(centValue>30 && centValue<60.){
1117              if(leadtrack->Pt()>0.){
1118           fh2JetsumHT3R2ap->Fill(ptcorr,sumpt2ap/rhoin2);
1119           fh2JetsumHT3R4ap->Fill(ptcorr,sumpt4ap/rhoin4);
1120           fh2JetsumHT3R6ap->Fill(ptcorr,sumpt6ap/rhoin6);
1121           fh2JetsumHT3R8ap->Fill(ptcorr,sumpt8ap/rhoin8);
1122           fh2JetsumHT3R10ap->Fill(ptcorr,sumpt10ap/rhoin10);}
1123             if(leadtrack->Pt()>6.){
1124           fh2JetsumHT3R2bp->Fill(ptcorr,sumpt2bp/rhoin2);
1125           fh2JetsumHT3R4bp->Fill(ptcorr,sumpt4bp/rhoin4);
1126           fh2JetsumHT3R6bp->Fill(ptcorr,sumpt6bp/rhoin6);
1127           fh2JetsumHT3R8bp->Fill(ptcorr,sumpt8bp/rhoin8);
1128           fh2JetsumHT3R10bp->Fill(ptcorr,sumpt10bp/rhoin10);}}}
1129
1130
1131
1132                         if(partback->Pt()<1.){
1133             if(centValue<10.){
1134               if(leadtrack->Pt()>0.){
1135           fh2JetsumHT3R2aaa->Fill(ptcorr,sumpt2aaa/rhoin2);
1136           fh2JetsumHT3R4aaa->Fill(ptcorr,sumpt4aaa/rhoin4);
1137           fh2JetsumHT3R6aaa->Fill(ptcorr,sumpt6aaa/rhoin6);
1138           fh2JetsumHT3R8aaa->Fill(ptcorr,sumpt8aaa/rhoin8);
1139           fh2JetsumHT3R10aaa->Fill(ptcorr,sumpt10aaa/rhoin10);}
1140             if(leadtrack->Pt()>6.){
1141           fh2JetsumHT3R2bbb->Fill(ptcorr,sumpt2bbb/rhoin2);
1142           fh2JetsumHT3R4bbb->Fill(ptcorr,sumpt4bbb/rhoin4);
1143           fh2JetsumHT3R6bbb->Fill(ptcorr,sumpt6bbb/rhoin6);
1144           fh2JetsumHT3R8bbb->Fill(ptcorr,sumpt8bbb/rhoin8);
1145           fh2JetsumHT3R10bbb->Fill(ptcorr,sumpt10bbb/rhoin10);}}
1146
1147                     if(centValue>30 && centValue<60.){
1148              if(leadtrack->Pt()>0.){
1149           fh2JetsumHT3R2aaap->Fill(ptcorr,sumpt2aaap/rhoin2);
1150           fh2JetsumHT3R4aaap->Fill(ptcorr,sumpt4aaap/rhoin4);
1151           fh2JetsumHT3R6aaap->Fill(ptcorr,sumpt6aaap/rhoin6);
1152           fh2JetsumHT3R8aaap->Fill(ptcorr,sumpt8aaap/rhoin8);
1153           fh2JetsumHT3R10aaap->Fill(ptcorr,sumpt10aaap/rhoin10);}
1154             if(leadtrack->Pt()>6.){
1155           fh2JetsumHT3R2bbbp->Fill(ptcorr,sumpt2bbbp/rhoin2);
1156           fh2JetsumHT3R4bbbp->Fill(ptcorr,sumpt4bbbp/rhoin4);
1157           fh2JetsumHT3R6bbbp->Fill(ptcorr,sumpt6bbbp/rhoin6);
1158           fh2JetsumHT3R8bbbp->Fill(ptcorr,sumpt8bbbp/rhoin8);
1159           fh2JetsumHT3R10bbbp->Fill(ptcorr,sumpt10bbbp/rhoin10);}}}
1160
1161           }
1162
1163
1164
1165
1166
1167
1168
1169
1170    }
1171    if(injet>0) fh2JetDensity->Fill(ParticleList.GetEntries(),injet/accep);
1172    if(injet4>0)fh2JetDensityA4->Fill(ParticleList.GetEntries(),injet4/accep);
1173           //end of jet loop
1174
1175
1176
1177
1178           if(fDoEventMixing){
1179             //check before if the trigger exists
1180             // fTrigBuffer[i][0] = zvtx
1181             // fTrigBuffer[i][1] = phi
1182             // fTrigBuffer[i][2] = eta
1183             // fTrigBuffer[i][3] = pt_jet
1184             // fTrigBuffer[i][4] = pt_trig
1185             // fTrigBuffer[i][5]= pt_track_in
1186             // fTrigBuffer[i][6]= centrality
1187             if(fTindex==11) fTindex=0;
1188             if(fTrigBuffer[fTindex][3]>0){
1189             if (TMath::Abs(fTrigBuffer[fTindex][0]-primVtx->GetZ()<2.)){
1190             if (TMath::Abs(fTrigBuffer[fTindex][6]-centValue<10)){
1191
1192                         for(int it = 0;it<nT;++it){
1193                         AliVParticle *part = (AliVParticle*)ParticleList.At(it);
1194                         Double_t DPhi = fTrigBuffer[fTindex][1] - part->Phi();
1195                         Double_t DEta = fTrigBuffer[fTindex][2] - part->Eta();
1196                         Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta);
1197                         if(DPhi<-0.5*TMath::Pi()) DPhi+=2.*TMath::Pi();
1198                         if(DPhi>3./2.*TMath::Pi()) DPhi-=2.*TMath::Pi();
1199                         Double_t triggerEntries[8] = {centValue,fTrigBuffer[fTindex][3],part->Pt(),DR,DEta,DPhi,fTrigBuffer[fTindex][4],fTrigBuffer[fTindex][5]};
1200                         fhnMixedEvents->Fill(triggerEntries);
1201                         }
1202                         fNevents=fNevents+1;
1203                         if(fNevents==9) {fTindex=fTindex+1;
1204                         fNevents=0;}
1205             }}}
1206
1207
1208                // Copy the triggers from the current event into the buffer.
1209                //again, only if the trigger exists:
1210                 if(trigJet>-1){
1211                 AliAODJet* jetT = (AliAODJet*)(fListJets[0]->At(trigJet));
1212                 AliVParticle *partL = (AliVParticle*)ParticleList.At(trigInTrack);
1213                 AliVParticle *partT = (AliVParticle*)ParticleList.At(trigBBTrack);
1214                 fTrigBuffer[fTrigBufferIndex][0] = primVtx->GetZ();
1215                 fTrigBuffer[fTrigBufferIndex][1] = jetT->Phi();
1216                 fTrigBuffer[fTrigBufferIndex][2] = jetT->Eta();
1217                 fTrigBuffer[fTrigBufferIndex][3] = jetT->Pt()-rho*jetT->EffectiveAreaCharged();
1218                 fTrigBuffer[fTrigBufferIndex][4] = partT->Pt();
1219                 fTrigBuffer[fTrigBufferIndex][5] = partL->Pt();
1220                 fTrigBuffer[fTrigBufferIndex][6] = centValue;
1221                 fTrigBufferIndex++;
1222                 if(fTrigBufferIndex==9) fTrigBufferIndex=0;
1223                 }
1224           }
1225
1226
1227
1228
1229
1230      //////////////////ANGULAR STRUCTURE//////////////////////////////////////
1231
1232      //tracks up to R=0.8 distant from the jet axis
1233    //   if(fAngStructCloseTracks==1){
1234    //    TList CloseTrackList;
1235    //    Int_t nn=GetListOfTracksCloseToJet(&CloseTrackList,jetbig);
1236    //    Double_t difR=0.04;
1237    //    for(Int_t l=0;l<15;l++){
1238    //    Double_t rr=l*0.1+0.1;
1239    //     for(int it = 0;it<nn;++it){
1240    //         AliVParticle *part1 = (AliVParticle*)CloseTrackList.At(it);
1241    //         for(int itu=it+1;itu<CloseTrackList.GetEntries();itu++){
1242    //         AliVParticle *part2 = (AliVParticle*)CloseTrackList.At(itu);
1243    //         Double_t ptm=part1->Pt();
1244    //         Double_t ptn=part2->Pt();
1245    //         Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
1246    //                    Rnm=TMath::Sqrt(Rnm);
1247    //                    Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
1248    //                    Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
1249    //                    if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
1250    //                                                   down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
1251    //                    if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
1252    //                                                   down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
1253    //                    if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
1254    //                                                    down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
1255    //                    if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
1256    //                   down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}
1257    //   }
1258
1259    //   //only jet constituents
1260    //    if(fAngStructCloseTracks==2){
1261
1262    //    Double_t difR=0.04;
1263    //    for(Int_t l=0;l<15;l++){
1264    //    Double_t rr=l*0.1+0.1;
1265
1266
1267    //    AliAODTrack* part1;
1268    //    AliAODTrack* part2;
1269
1270    //        TRefArray *genTrackListb = jetbig->GetRefTracks();
1271    //        Int_t nTracksGenJetb = genTrackListb->GetEntriesFast();
1272
1273
1274
1275    //        for(Int_t it=0; it<nTracksGenJetb; ++it){
1276    //           part1 = (AliAODTrack*)(genTrackListb->At(it));
1277    //         for(Int_t itu=0; itu<nTracksGenJetb; ++itu){
1278    //           part2 = (AliAODTrack*)(genTrackListb->At(itu));
1279    //         Double_t ptm=part1->Pt();
1280    //         Double_t ptn=part2->Pt();
1281    //         Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
1282    //                    Rnm=TMath::Sqrt(Rnm);
1283    //                    Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
1284    //                    Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
1285    //                    if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
1286    //                                                   down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
1287    //                    if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
1288    //                                                   down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
1289    //                    if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
1290    //                                                    down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
1291    //                    if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
1292    //                   down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}}
1293    // }
1294      // //end loop over R=0.4 jets
1295      // if(fAngStructCloseTracks>0){
1296      // for(Int_t l=0;l<15;l++){
1297      // Double_t rr=l*0.1+0.1;
1298      //    if(down1[l]!=0){
1299      //         if(centValue<10.)fh2AngStructpt1C10->Fill(rr,rr*up1[l]/down1[l]);
1300      //    if(centValue>20. && centValue<40.) fh2AngStructpt1C20->Fill(rr,rr*up1[l]/down1[l]);
1301      //    if(centValue>30. && centValue<60.) fh2AngStructpt1C30->Fill(rr,rr*up1[l]/down1[l]);
1302      //    if(centValue>60.) fh2AngStructpt1C60->Fill(rr,rr*up1[l]/down1[l]);}
1303      //    if(down2[l]!=0){
1304      //         if(centValue<10.) fh2AngStructpt2C10->Fill(rr,rr*up2[l]/down2[l]);
1305      //    if(centValue>20. && centValue<40.) fh2AngStructpt2C20->Fill(rr,rr*up2[l]/down2[l]);
1306      //    if(centValue>30. && centValue<60.) fh2AngStructpt2C30->Fill(rr,rr*up2[l]/down2[l]);
1307      //    if(centValue>60.) fh2AngStructpt2C60->Fill(rr,rr*up2[l]/down2[l]);}
1308      //    if(down3[l]!=0){
1309      //         if(centValue<10.) fh2AngStructpt3C10->Fill(rr,rr*up3[l]/down3[l]);
1310      //    if(centValue>20. && centValue<40.) fh2AngStructpt3C20->Fill(rr,rr*up3[l]/down3[l]);
1311      //    if(centValue>30. && centValue<60.) fh2AngStructpt3C30->Fill(rr,rr*up3[l]/down3[l]);
1312      //    if(centValue>60.) fh2AngStructpt3C60->Fill(rr,rr*up3[l]/down3[l]);}
1313      //    if(down4[l]!=0){
1314      //         if(centValue<10.) fh2AngStructpt4C10->Fill(rr,rr*up4[l]/down4[l]);
1315      //    if(centValue>20. && centValue<40.) fh2AngStructpt4C20->Fill(rr,rr*up4[l]/down4[l]);
1316      //    if(centValue>30. && centValue<60.) fh2AngStructpt4C30->Fill(rr,rr*up4[l]/down4[l]);
1317      //    if(centValue>60.) fh2AngStructpt4C60->Fill(rr,rr*up4[l]/down4[l]);}}}
1318
1319
1320
1321
1322
1323
1324
1325    PostData(1, fOutputList);
1326 }
1327
1328 void AliAnalysisTaskJetCore::Terminate(const Option_t *)
1329 {
1330    // Draw result to the screen
1331    // Called once at the end of the query
1332
1333    if (!GetOutputData(1))
1334    return;
1335 }
1336
1337
1338
1339
1340
1341
1342 Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
1343
1344      Int_t iCount = 0;
1345      AliAODEvent *aod = 0;
1346      if(!fESD)aod = fAODIn;
1347      else aod = fAODOut;
1348      Int_t index=-1;
1349      Double_t ptmax=-10;
1350     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1351       AliAODTrack *tr = aod->GetTrack(it);
1352       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1353       if(TMath::Abs(tr->Eta())>0.9)continue;
1354       if(tr->Pt()<0.15)continue;
1355       list->Add(tr);
1356       iCount++;
1357       if(tr->Pt()>ptmax){ ptmax=tr->Pt();
1358       index=iCount-1;}
1359
1360     }
1361
1362
1363   return index;
1364
1365 }
1366
1367    Int_t  AliAnalysisTaskJetCore::GetHardestTrackBackToJet(AliAODJet *jetbig){
1368
1369     AliAODEvent *aod = 0;
1370     if(!fESD)aod = fAODIn;
1371     else aod = fAODOut;
1372     Int_t index=-1;
1373     Double_t ptmax=-10;
1374     Double_t dphi=0;
1375     Double_t dif=0;
1376     Int_t iCount=0;
1377     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1378       AliAODTrack *tr = aod->GetTrack(it);
1379       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1380       if(TMath::Abs(tr->Eta())>0.9)continue;
1381       if(tr->Pt()<0.15)continue;
1382       iCount=iCount+1;
1383       dphi=RelativePhi(tr->Phi(),jetbig->Phi());
1384       if(TMath::Abs(dphi)<TMath::Pi()-0.2) continue;
1385       if(tr->Pt()>ptmax){ ptmax=tr->Pt();
1386       index=iCount-1;
1387       dif=dphi;  }}
1388
1389       return index;
1390
1391    }
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401  Int_t  AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
1402
1403     Int_t iCount = 0;
1404       AliAODEvent *aod = 0;
1405      if(!fESD)aod = fAODIn;
1406      else aod = fAODOut;
1407
1408     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1409       AliAODTrack *tr = aod->GetTrack(it);
1410       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1411       if(TMath::Abs(tr->Eta())>0.9)continue;
1412       if(tr->Pt()<0.15)continue;
1413       Double_t disR=jetbig->DeltaR(tr);
1414       if(disR>0.8)  continue;
1415       list->Add(tr);
1416       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
1417       iCount++;
1418     }
1419
1420    list->Sort();
1421    return iCount;
1422
1423 }
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435 Int_t AliAnalysisTaskJetCore::GetNInputTracks()
1436 {
1437
1438    Int_t nInputTracks = 0;
1439      AliAODEvent *aod = 0;
1440      if(!fESD)aod = fAODIn;
1441      else aod = fAODOut;
1442    TString jbname(fJetBranchName[1]);
1443    //needs complete event, use jets without background subtraction
1444    for(Int_t i=1; i<=3; ++i){
1445       if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
1446    }
1447    // use only HI event
1448    if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
1449    if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
1450
1451    if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
1452    TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(aod->FindListObject(jbname.Data()));
1453    if(!tmpAODjets){
1454       Printf("Jet branch %s not found", jbname.Data());
1455       Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
1456       return -1;
1457    }
1458
1459    for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
1460       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
1461       if(!jet) continue;
1462       TRefArray *trackList = jet->GetRefTracks();
1463       Int_t nTracks = trackList->GetEntriesFast();
1464       nInputTracks += nTracks;
1465       if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
1466    }
1467    if(fDebug) Printf("---> input tracks: %d", nInputTracks);
1468
1469    return nInputTracks;
1470 }
1471
1472
1473
1474 Double_t AliAnalysisTaskJetCore::RelativePhi(Double_t mphi,Double_t vphi){
1475
1476   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1477   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1478   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1479   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1480   double dphi = mphi-vphi;
1481   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1482   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1483   return dphi;//dphi in [-Pi, Pi]
1484 }
1485
1486
1487
1488 THnSparse* AliAnalysisTaskJetCore::NewTHnSparseF(const char* name, UInt_t entries)
1489 {
1490    // generate new THnSparseF, axes are defined in GetDimParams()
1491
1492    Int_t count = 0;
1493    UInt_t tmp = entries;
1494    while(tmp!=0){
1495       count++;
1496       tmp = tmp &~ -tmp;  // clear lowest bit
1497    }
1498
1499    TString hnTitle(name);
1500    const Int_t dim = count;
1501    Int_t nbins[dim];
1502    Double_t xmin[dim];
1503    Double_t xmax[dim];
1504
1505    Int_t i=0;
1506    Int_t c=0;
1507    while(c<dim && i<32){
1508       if(entries&(1<<i)){
1509
1510          TString label("");
1511          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1512          hnTitle += Form(";%s",label.Data());
1513          c++;
1514       }
1515
1516       i++;
1517    }
1518    hnTitle += ";";
1519
1520    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1521 }
1522
1523 void AliAnalysisTaskJetCore::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1524 {
1525    // stores label and binning of axis for THnSparse
1526
1527    const Double_t pi = TMath::Pi();
1528
1529    switch(iEntry){
1530
1531    case 0:
1532       label = "V0 centrality (%)";
1533
1534          nbins = 10;
1535          xmin = 0.;
1536          xmax = 100.;
1537          break;
1538
1539
1540    case 1:
1541       label = "corrected jet pt";
1542          nbins = 20;
1543          xmin = 0.;
1544          xmax = 200.;
1545           break;
1546
1547
1548    case 2:
1549       label = "track pT";
1550
1551          nbins = 9;
1552          xmin = 0.;
1553          xmax = 150;
1554          break;
1555
1556
1557     case 3:
1558       label = "deltaR";
1559       nbins = 15;
1560       xmin = 0.;
1561       xmax = 1.5;
1562       break;
1563
1564
1565
1566    case 4:
1567       label = "deltaEta";
1568       nbins = 8;
1569       xmin = -1.6;
1570       xmax = 1.6;
1571       break;
1572
1573
1574   case 5:
1575       label = "deltaPhi";
1576       nbins = 90;
1577       xmin = -0.5*pi;
1578       xmax = 1.5*pi;
1579       break;
1580
1581
1582
1583     case 6:
1584       label = "leading track";
1585       nbins = 13;
1586       xmin = 0;
1587       xmax = 50;
1588       break;
1589
1590      case 7:
1591
1592       label = "trigger track";
1593       nbins =10;
1594       xmin = 0;
1595       xmax = 50;
1596       break;
1597
1598
1599
1600
1601
1602
1603
1604
1605    }
1606
1607 }
1608