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