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