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