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