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