PWGJE/AliAnalysisTaskJetCore.cxx

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