PWGJE/AliAnalysisTaskJetCore.cxx

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