PWGJE/AliAnalysisTaskJetCore.cxx

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