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=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 = 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 = aod->GetTrack(it);
1170       Bool_t bGood = false;
1171       if(fFilterType == 0)bGood = true;
1172       else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1173       else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1174       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1175       if(fRequireITSRefit==1){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
1176       if(bGood==false) continue;
1177       if (fApplySharedClusterCut) {
1178            Double_t frac = Double_t(tr->GetTPCnclsS()) /Double_t(tr->GetTPCncls());
1179            if (frac > 0.4) continue;
1180       }
1181      if(TMath::Abs(tr->Eta())>0.9)continue;
1182       if(tr->Pt()<0.15)continue;
1183       list->Add(tr);
1184       iCount++;
1185       if(fFilterType==2 && fFilterMaskBestPt>0){// only set the trigger track index for good quality tracks
1186         if(tr->TestFilterBit(fFilterMaskBestPt)){
1187           if(tr->Pt()>ptmax){
1188             ptmax=tr->Pt();
1189             index=iCount-1;
1190           }
1191         }
1192       }
1193       else{
1194         if(tr->Pt()>ptmax){
1195           ptmax=tr->Pt();
1196           index=iCount-1;
1197         }
1198       }
1199      }
1200
1201
1202     // else if (type == kTrackAODMCCharged) {
1203     // TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1204     // if(!tca)return iCount;
1205     // for(int it = 0;it < tca->GetEntriesFast();++it){
1206     //   AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1207     //   if(!part)continue;
1208     //   if(part->Pt()<0.15)continue;
1209     //   if(!part->IsPhysicalPrimary())continue;
1210     //   if(part->Charge()==0)continue;
1211     //   if(TMath::Abs(part->Eta())>0.9)continue;
1212     //   list->Add(part);
1213     //   iCount++;
1214     //   if(part->Pt()>ptmax){ ptmax=part->Pt();
1215     //  index=iCount-1;}}}
1216       return index;
1217
1218 }
1219
1220
1221
1222 Int_t  AliAnalysisTaskJetCore::SelectTrigger(TList *list,Double_t minT,Double_t maxT,Int_t &number){
1223      Int_t iCount = 0;
1224      AliAODEvent *aod = 0;
1225      if(!fESD)aod = fAODIn;
1226      else aod = fAODOut;
1227      if(!aod)return 0;
1228      Int_t index=-1;
1229      Int_t triggers[100];
1230      for(Int_t cr=0;cr<100;cr++){triggers[cr]=-1;}
1231      Int_t im=0;
1232      for(int it = 0;it < aod->GetNumberOfTracks();++it){
1233       AliAODTrack *tr = aod->GetTrack(it);
1234       Bool_t bGood = false;
1235       if(fFilterType == 0)bGood = true;
1236       else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1237       else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1238       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1239        if(fRequireITSRefit==1){if((tr->GetStatus()&AliESDtrack::kITSrefit)==0)continue;}
1240       if(bGood==false) continue;
1241       if (fApplySharedClusterCut) {
1242            Double_t frac = Double_t(tr->GetTPCnclsS()) /Double_t(tr->GetTPCncls());
1243            if (frac > 0.4) continue;
1244       }
1245       if(TMath::Abs(tr->Eta())>0.9)continue;
1246       if(tr->Pt()<0.15)continue;
1247       list->Add(tr);
1248       iCount++;
1249
1250       if(tr->Pt()>=minT && tr->Pt()<maxT){
1251         triggers[im]=iCount-1;
1252         im=im+1;}
1253
1254      }
1255       number=im;
1256       Int_t rd=0;
1257       if(im==0) rd=0;
1258       if(im>0) rd=fRandom->Integer(im);
1259       index=triggers[rd];
1260
1261
1262
1263
1264
1265       return index;
1266
1267 }
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284    Int_t  AliAnalysisTaskJetCore::GetHardestTrackBackToJet(AliAODJet *jetbig){
1285
1286     AliAODEvent *aod = 0;
1287     if(!fESD)aod = fAODIn;
1288     else aod = fAODOut;
1289     Int_t index=-1;
1290     Double_t ptmax=-10;
1291     Double_t dphi=0;
1292     Double_t dif=0;
1293     Int_t iCount=0;
1294     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1295       AliAODTrack *tr = aod->GetTrack(it);
1296       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1297       if(TMath::Abs(tr->Eta())>0.9)continue;
1298       if(tr->Pt()<0.15)continue;
1299       iCount=iCount+1;
1300       dphi=RelativePhi(tr->Phi(),jetbig->Phi());
1301       if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;
1302       if(tr->Pt()>ptmax){ ptmax=tr->Pt();
1303       index=iCount-1;
1304       dif=dphi;  }}
1305
1306       return index;
1307
1308    }
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318  Int_t  AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
1319
1320     Int_t iCount = 0;
1321       AliAODEvent *aod = 0;
1322      if(!fESD)aod = fAODIn;
1323      else aod = fAODOut;
1324
1325       for(int it = 0;it < aod->GetNumberOfTracks();++it){
1326       AliAODTrack *tr = aod->GetTrack(it);
1327       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1328       if(TMath::Abs(tr->Eta())>0.9)continue;
1329       if(tr->Pt()<0.15)continue;
1330       Double_t disR=jetbig->DeltaR(tr);
1331       if(disR>0.8)  continue;
1332       list->Add(tr);
1333       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
1334       iCount++;
1335     }
1336
1337    list->Sort();
1338    return iCount;
1339
1340 }
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352 Int_t AliAnalysisTaskJetCore::GetNInputTracks()
1353 {
1354
1355    Int_t nInputTracks = 0;
1356      AliAODEvent *aod = 0;
1357      if(!fESD)aod = fAODIn;
1358      else aod = fAODOut;
1359    TString jbname(fJetBranchName[1]);
1360    //needs complete event, use jets without background subtraction
1361    for(Int_t i=1; i<=3; ++i){
1362       if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
1363    }
1364    // use only HI event
1365    if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
1366    if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
1367
1368    if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
1369    TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(aod->FindListObject(jbname.Data()));
1370    if(!tmpAODjets){
1371       Printf("Jet branch %s not found", jbname.Data());
1372       Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
1373       return -1;
1374    }
1375
1376    for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
1377       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
1378       if(!jet) continue;
1379       TRefArray *trackList = jet->GetRefTracks();
1380       Int_t nTracks = trackList->GetEntriesFast();
1381       nInputTracks += nTracks;
1382       if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
1383    }
1384    if(fDebug) Printf("---> input tracks: %d", nInputTracks);
1385
1386    return nInputTracks;
1387 }
1388
1389
1390
1391 Double_t AliAnalysisTaskJetCore::RelativePhi(Double_t mphi,Double_t vphi){
1392
1393   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1394   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1395   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1396   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1397   double dphi = mphi-vphi;
1398   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1399   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1400   return dphi;//dphi in [-Pi, Pi]
1401 }
1402
1403 Int_t AliAnalysisTaskJetCore::GetPhiBin(Double_t phi)
1404 {
1405     Int_t phibin=-1;
1406     if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1407     Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1408     phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1409     if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1410     return phibin;
1411 }
1412
1413
1414
1415
1416 THnSparse* AliAnalysisTaskJetCore::NewTHnSparseF(const char* name, UInt_t entries)
1417 {
1418    // generate new THnSparseF, axes are defined in GetDimParams()
1419
1420    Int_t count = 0;
1421    UInt_t tmp = entries;
1422    while(tmp!=0){
1423       count++;
1424       tmp = tmp &~ -tmp;  // clear lowest bit
1425    }
1426
1427    TString hnTitle(name);
1428    const Int_t dim = count;
1429    Int_t nbins[dim];
1430    Double_t xmin[dim];
1431    Double_t xmax[dim];
1432
1433    Int_t i=0;
1434    Int_t c=0;
1435    while(c<dim && i<32){
1436       if(entries&(1<<i)){
1437
1438          TString label("");
1439          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1440          hnTitle += Form(";%s",label.Data());
1441          c++;
1442       }
1443
1444       i++;
1445    }
1446    hnTitle += ";";
1447
1448    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1449 }
1450
1451 void AliAnalysisTaskJetCore::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1452 {
1453    // stores label and binning of axis for THnSparse
1454
1455    const Double_t pi = TMath::Pi();
1456
1457    switch(iEntry){
1458
1459    case 0:
1460       label = "V0 centrality (%)";
1461
1462          nbins = 10;
1463          xmin = 0.;
1464          xmax = 100.;
1465          break;
1466
1467
1468    case 1:
1469       label = "corrected jet pt";
1470          nbins = 20;
1471          xmin = 0.;
1472          xmax = 200.;
1473           break;
1474
1475
1476    case 2:
1477       label = "track pT";
1478
1479          nbins = 9;
1480          xmin = 0.;
1481          xmax = 150;
1482          break;
1483
1484
1485     case 3:
1486       label = "deltaR";
1487       nbins = 15;
1488       xmin = 0.;
1489       xmax = 1.5;
1490       break;
1491
1492
1493
1494    case 4:
1495       label = "deltaEta";
1496       nbins = 8;
1497       xmin = -1.6;
1498       xmax = 1.6;
1499       break;
1500
1501
1502   case 5:
1503       label = "deltaPhi";
1504       nbins = 90;
1505       xmin = -0.5*pi;
1506       xmax = 1.5*pi;
1507       break;
1508
1509
1510
1511     case 6:
1512       label = "leading track";
1513       nbins = 13;
1514       xmin = 0;
1515       xmax = 50;
1516       break;
1517
1518      case 7:
1519
1520       label = "trigger track";
1521       nbins =10;
1522       xmin = 0;
1523       xmax = 50;
1524       break;
1525
1526
1527
1528
1529
1530
1531
1532
1533    }
1534
1535 }
1536