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","",100,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.07","",100,0.,4000.,100,0.,5.,25,0.,50.);
 416     fh3JetDensityA4=new TH3F("Jet density vs multiplicity A>0.4","",100,0.,4000.,100,0.,5.,25,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]     = {100,6, 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
 489              //change binning in jet area
 490      Double_t *xPt6 = new Double_t[7];
 491      xPt6[0] = 0.;
 492      xPt6[1]=0.07;
 493      xPt6[2]=0.2;
 494      xPt6[3]=0.4;
 495      xPt6[4]=0.6;
 496      xPt6[5]=0.8;
 497      xPt6[6]=1;
 498     fHJetSpec->SetBinEdges(1,xPt6);
 499     delete [] xPt6;
 500
 501
 502
 503
 504
 505         fOutputList->Add(fHJetSpec);
 506
 507
 508         if(fRunAnaAzimuthalCorrelation)
 509           {
 510             fhTTPt = new TH2F("fhTTPt","Trigger track p_{T} vs centrality",10,0,100,100,0,100);
 511             fOutputList->Add(fhTTPt);
 512
 513             const Int_t dimCor = 5;
 514             const Int_t nBinsCor[dimCor]     = {50, 200, 100,              8,   10};
 515             const Double_t lowBinCor[dimCor] = {0,  -50, -0.5*TMath::Pi(), 0,   0};
 516             const Double_t hiBinCor[dimCor]  = {50, 150, 1.5*TMath::Pi(),  0.8, 100};
 517             fHJetPhiCorr = new THnSparseF("fHJetPhiCorr","TT p_{T} vs jet p_{T} vs dPhi vs area vs centrality",dimCor,nBinsCor,lowBinCor,hiBinCor);
 518             fOutputList->Add(fHJetPhiCorr);
 519           }
 520
 521
 522
 523
 524    // =========== Switch on Sumw2 for all histos ===========
 525    for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
 526       TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
 527       if (h1){
 528          h1->Sumw2();
 529          continue;
 530       }
 531     THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
 532       if (hn){
 533          hn->Sumw2();
 534       }
 535    }
 536    TH1::AddDirectory(oldStatus);
 537
 538    PostData(1, fOutputList);
 539 }
 540
 541 void AliAnalysisTaskJetCore::UserExec(Option_t *)
 542 {
 543
 544
 545    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 546       AliError("Jet branch name not set.");
 547       return;
 548    }
 549
 550    fESD=dynamic_cast<AliESDEvent*>(InputEvent());
 551    if (!fESD) {
 552       AliError("ESD not available");
 553       fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
 554    }
 555       fAODOut = dynamic_cast<AliAODEvent*>(AODEvent());
 556
 557        static AliAODEvent* aod = 0;
 558        // take all other information from the aod we take the tracks from
 559        if(!aod){
 560        if(!fESD)aod = fAODIn;
 561        else aod = fAODOut;}
 562
 563
 564
 565     if(fNonStdFile.Length()!=0){
 566     // case that we have an AOD extension we need can fetch the jets from the extended output
 567     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 568     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 569     if(!fAODExtension){
 570       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
 571     }}
 572
 573
 574    // -- event selection --
 575    fHistEvtSelection->Fill(1); // number of events before event selection
 576
 577    // physics selection
 578    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
 579    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 580          std::cout<<inputHandler->IsEventSelected()<<" "<<fOfflineTrgMask<<std::endl;
 581    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
 582       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
 583       fHistEvtSelection->Fill(2);
 584       PostData(1, fOutputList);
 585       return;
 586    }
 587
 588    // vertex selection
 589    if(!aod){
 590      if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__);
 591      fHistEvtSelection->Fill(3);
 592       PostData(1, fOutputList);
 593    }
 594    AliAODVertex* primVtx = aod->GetPrimaryVertex();
 595
 596    if(!primVtx){
 597      if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__);
 598      fHistEvtSelection->Fill(3);
 599      PostData(1, fOutputList);
 600      return;
 601    }
 602
 603    Int_t nTracksPrim = primVtx->GetNContributors();
 604    if ((nTracksPrim < fMinContribVtx) ||
 605          (primVtx->GetZ() < fVtxZMin) ||
 606          (primVtx->GetZ() > fVtxZMax) ){
 607       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
 608       fHistEvtSelection->Fill(3);
 609       PostData(1, fOutputList);
 610       return;
 611    }
 612
 613    // event class selection (from jet helper task)
 614    Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
 615    if(fDebug) Printf("Event class %d", eventClass);
 616    if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
 617       fHistEvtSelection->Fill(4);
 618       PostData(1, fOutputList);
 619       return;
 620    }
 621
 622    // centrality selection
 623    AliCentrality *cent = 0x0;
 624    Double_t centValue = 0.;
 625    if(fIsPbPb){
 626    if(fESD) {cent = fESD->GetCentrality();
 627      if(cent) centValue = cent->GetCentralityPercentile("V0M");}
 628    else     centValue=aod->GetHeader()->GetCentrality();
 629
 630    if(fDebug) printf("centrality: %f\n", centValue);
 631       if (centValue < fCentMin || centValue > fCentMax){
 632       fHistEvtSelection->Fill(4);
 633       PostData(1, fOutputList);
 634       return;
 635       }}
 636
 637
 638    fHistEvtSelection->Fill(0);
 639    // accepted events
 640    // -- end event selection --
 641
 642    // get background
 643    AliAODJetEventBackground* externalBackground = 0;
 644    if(fAODOut&&!externalBackground&&fBackgroundBranch.Length()){
 645       externalBackground =  (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
 646       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 647    }
 648    if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
 649      externalBackground =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
 650       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 651    }
 652
 653     if(fAODIn&&!externalBackground&&fBackgroundBranch.Length()){
 654       externalBackground =  (AliAODJetEventBackground*)(fAODIn->FindListObject(fBackgroundBranch.Data()));
 655       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 656     }
 657
 658    Float_t rho = 0;
 659
 660    if(fFlagRandom==0){
 661      if(externalBackground)rho = externalBackground->GetBackground(0);}
 662    if(fFlagRandom==1){
 663       if(externalBackground)rho = externalBackground->GetBackground(2);}
 664
 665    // fetch jets
 666    TClonesArray *aodJets[2];
 667    aodJets[0]=0;
 668    if(fAODOut&&!aodJets[0]){
 669    aodJets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data()));
 670    aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data()));  }
 671    if(fAODExtension && !aodJets[0]){
 672    aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
 673    aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));  }
 674      if(fAODIn&&!aodJets[0]){
 675    aodJets[0] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[0].Data()));
 676    aodJets[1] = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName[1].Data()));  }
 677
 678
 679    //Double_t ptsub[aodJets[0]->GetEntriesFast()];
 680    //Int_t inord[aodJets[0]->GetEntriesFast()];
 681    //for(Int_t n=0;n<aodJets[0]->GetEntriesFast();n++){
 682    //  ptsub[n]=0;
 683    //  inord[n]=0;}
 684
 685    TList ParticleList;
 686    Int_t nT = GetListOfTracks(&ParticleList);
 687      for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
 688       fListJets[iJetType]->Clear();
 689       if (!aodJets[iJetType]) continue;
 690       if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
 691       for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
 692          AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
 693          if (jet) fListJets[iJetType]->Add(jet);
 694          // if(iJetType==0){
 695          // ptsub[iJet]=jet->Pt()-rho*jet->EffectiveAreaCharged();}
 696       }}
 697
 698    Double_t etabig=0;
 699    Double_t ptbig=0;
 700    Double_t areabig=0;
 701    Double_t phibig=0.;
 702    Double_t etasmall=0;
 703    Double_t ptsmall=0;
 704    Double_t areasmall=0;
 705    Double_t phismall=0.;
 706
 707
 708
 709    Int_t iCount=0;
 710    Int_t trigJet=-1;
 711    Int_t trigBBTrack=-1;
 712    Int_t trigInTrack=-1;
 713    fRPAngle = aod->GetHeader()->GetEventplane();
 714
 715    AliVParticle *partback = (AliVParticle*)ParticleList.At(nT);
 716    if(!partback){
 717    PostData(1, fOutputList);
 718    return;}
 719
 720
 721    //for(Int_t tt=0;tt<ParticleList.GetEntries();tt++){
 722    //if(fFlagOnlyHardest!=0){if(tt!=nT) continue;}
 723    //AliVParticle *partback = (AliVParticle*)ParticleList.At(tt);
 724    Double_t accep=2.*TMath::Pi()*1.8;
 725    Int_t injet4=0;
 726    Int_t injet=0;
 727    if(fSemigoodCorrect){
 728    Double_t disthole=RelativePhi(partback->Phi(),fHolePos);
 729    if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-0.6){
 730    PostData(1, fOutputList);
 731    return;}
 732
 733    }
 734
 735    fh2Ntriggers->Fill(centValue,partback->Pt());
 736    Double_t phiBinT = RelativePhi(partback->Phi(),fRPAngle);
 737    if(centValue<20.) fh2RPTC20->Fill(TMath::Abs(phiBinT),partback->Pt());
 738    if(centValue<10.) fh2RPTC10->Fill(TMath::Abs(phiBinT),partback->Pt());
 739
 740
 741
 742    for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
 743            AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
 744            etabig  = jetbig->Eta();
 745            phibig  = jetbig->Phi();
 746            ptbig   = jetbig->Pt();
 747            if(ptbig==0) continue;
 748            Double_t phiBin = RelativePhi(phibig,fRPAngle);
 749            areabig = jetbig->EffectiveAreaCharged();
 750            Double_t ptcorr=ptbig-rho*areabig;
 751            if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
 752            if(areabig>=0.07) injet=injet+1;
 753            if(areabig>=0.4) injet4=injet4+1;
 754            Double_t dphi=RelativePhi(partback->Phi(),phibig);
 755
 756            if(fFlagEtaBkg==1){
 757            Double_t etadif= partback->Eta()-etabig;
 758            if(TMath::Abs(etadif)<=0.5){
 759
 760            if(centValue<20.) fh3JetTrackC20->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
 761            if(centValue>30. && centValue<60.) fh3JetTrackC3060->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));}}
 762            if(fFlagEtaBkg==0){
 763            if(centValue<20.) fh3JetTrackC20->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
 764            if(centValue>30. && centValue<60.) fh3JetTrackC3060->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));}
 765
 766
 767            if(fFlagJetHadron==0){
 768            if(fFlagPhiBkg==1) if((TMath::Abs(dphi)<TMath::Pi()/2.-0.1)||(TMath::Abs(dphi)>TMath::Pi()/2.+0.1)) continue;
 769            if(fFlagPhiBkg==0) if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;
 770            if(fFlagPhiBkg==2) if(TMath::Abs(dphi)<TMath::Pi()-0.7) continue;
 771            if(fFlagPhiBkg==3) if(TMath::Abs(dphi)<TMath::Pi()-0.5) continue;}
 772
 773            if(fFlagJetHadron!=0) if(TMath::Abs(dphi)>0.4) continue;
 774
 775
 776            if(centValue<10.) fh2RPJetsC10->Fill(TMath::Abs(phiBin), ptcorr);
 777            if(centValue<20.) fh2RPJetsC20->Fill(TMath::Abs(phiBin), ptcorr);
 778                    Double_t dismin=100.;
 779                    Double_t ptmax=-10.;
 780                    Int_t index1=-1;
 781                    Int_t index2=-1;
 782
 783                    Float_t phitt=partback->Phi();
 784                    if(phitt<0)phitt+=TMath::Pi()*2.;
 785                    Int_t phiBintt = GetPhiBin(phitt-fRPAngle);
 786
 787                    Double_t fillspec[] = {centValue,jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt(),phiBintt};
 788                   fHJetSpec->Fill(fillspec);
 789
 790
 791
 792                    if(ptcorr<=0) continue;
 793
 794                        AliAODTrack* leadtrack=0;
 795                        Int_t ippt=0;
 796                        Double_t ppt=-10;
 797                        if(fFlagJetHadron==0){
 798                          TRefArray *genTrackList = jetbig->GetRefTracks();
 799                          Int_t nTracksGenJet = genTrackList->GetEntriesFast();
 800                          AliAODTrack* genTrack;
 801                          for(Int_t ir=0; ir<nTracksGenJet; ++ir){
 802                            genTrack = (AliAODTrack*)(genTrackList->At(ir));
 803                            if(genTrack->Pt()>ppt){ppt=genTrack->Pt();
 804                              ippt=ir;}}
 805                          leadtrack=(AliAODTrack*)(genTrackList->At(ippt));
 806                          if(!leadtrack) continue;
 807                        }
 808
 809
 810
 811                        AliVParticle* leadtrackb=0;
 812                        if(fFlagJetHadron!=0){
 813                          Int_t nTb = GetHardestTrackBackToJet(jetbig);
 814                          leadtrackb = (AliVParticle*)ParticleList.At(nTb);
 815                          if(!leadtrackb) continue;
 816                        }
 817
 818
 819
 820
 821
 822                        //store one trigger info
 823                        if(iCount==0){
 824                          trigJet=i;
 825                          trigBBTrack=nT;
 826                          trigInTrack=ippt;
 827                          iCount=iCount+1;}
 828
 829
 830                        if(fCheckMethods){
 831                          for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
 832                            AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
 833                            etasmall  = jetsmall->Eta();
 834                            phismall = jetsmall->Phi();
 835                            ptsmall   = jetsmall->Pt();
 836                            areasmall = jetsmall->EffectiveAreaCharged();
 837                            Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
 838                            tmpDeltaR=TMath::Sqrt(tmpDeltaR);
 839                            //Fraction in the jet core
 840                            if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;
 841                              index2=j;}
 842                     if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
 843                       index1=j;}} //en of loop over R=0.2 jets
 844                 //method1:most concentric jet=core
 845                 if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));
 846                   if(centValue<10) fh2JetCoreMethod1C10->Fill(ptcorr,jetmethod1->Pt()/ptbig);
 847                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod1C20->Fill(ptcorr,jetmethod1->Pt()/ptbig);
 848                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod1C30->Fill(ptcorr,jetmethod1->Pt()/ptbig);
 849                   if(centValue>60) fh2JetCoreMethod1C60->Fill(ptcorr,jetmethod1->Pt()/ptbig); }
 850                 //method2:hardest contained jet=core
 851                 if(index2!=-1){
 852                   AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
 853                   if(centValue<10) fh2JetCoreMethod2C10->Fill(ptcorr,jetmethod2->Pt()/ptbig);
 854                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod2C20->Fill(ptcorr,jetmethod2->Pt()/ptbig);
 855                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptcorr,jetmethod2->Pt()/ptbig);
 856                   if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }}
 857                   if(centValue<10&&leadtrack) fh2Ntriggers2C10->Fill(leadtrack->Pt(),partback->Pt());
 858                   if(centValue<20&&leadtrack) fh2Ntriggers2C20->Fill(leadtrack->Pt(),partback->Pt());
 859          if(fDoEventMixing==0 && fFlagOnlyRecoil==0){
 860          for(int it = 0;it<ParticleList.GetEntries();++it){
 861           AliVParticle *part = (AliVParticle*)ParticleList.At(it);
 862           Double_t deltaR = jetbig->DeltaR(part);
 863           Double_t deltaEta = etabig-part->Eta();
 864
 865           Double_t deltaPhi=phibig-part->Phi();
 866           if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
 867           if(deltaPhi>3./2.*TMath::Pi()) deltaPhi-=2.*TMath::Pi();
 868           Double_t pTcont=0;
 869           if(fFlagJetHadron==0) pTcont=leadtrack->Pt();
 870           if(fFlagJetHadron!=0) pTcont=leadtrackb->Pt();
 871            Double_t jetEntries[8] = {centValue,ptcorr,part->Pt(),deltaR,deltaEta,deltaPhi,pTcont,partback->Pt()};
 872            fhnDeltaR->Fill(jetEntries);}
 873
 874
 875           }
 876
 877
 878           //end of track loop, we only do it if EM is switched off
 879
 880
 881
 882
 883
 884
 885
 886
 887
 888    }
 889    if(injet>0) fh3JetDensity->Fill(ParticleList.GetEntries(),injet/accep,partback->Pt());
 890    if(injet4>0)fh3JetDensityA4->Fill(ParticleList.GetEntries(),injet4/accep,partback->Pt());
 891           //end of jet loop
 892
 893    //}
 894
 895
 896           if(fDoEventMixing>0){
 897             //check before if the trigger exists
 898             // fTrigBuffer[i][0] = zvtx
 899             // fTrigBuffer[i][1] = phi
 900             // fTrigBuffer[i][2] = eta
 901             // fTrigBuffer[i][3] = pt_jet
 902             // fTrigBuffer[i][4] = pt_trig
 903             // fTrigBuffer[i][5]= centrality
 904             if(fTindex==10) fTindex=0;
 905             if(fTrigBuffer[fTindex][3]>0){
 906             if (TMath::Abs(fTrigBuffer[fTindex][0]-primVtx->GetZ()<2.)){
 907             if (TMath::Abs(fTrigBuffer[fTindex][5]-centValue<5)){
 908
 909                         for(int it = 0;it<nT;++it){
 910                         AliVParticle *part = (AliVParticle*)ParticleList.At(it);
 911                         Double_t DPhi = fTrigBuffer[fTindex][1] - part->Phi();
 912                         Double_t DEta = fTrigBuffer[fTindex][2] - part->Eta();
 913                         Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta);
 914                         if(DPhi<-0.5*TMath::Pi()) DPhi+=2.*TMath::Pi();
 915                         if(DPhi>3./2.*TMath::Pi()) DPhi-=2.*TMath::Pi();
 916                         Double_t triggerEntries[7] = {centValue,fTrigBuffer[fTindex][3],part->Pt(),DR,DEta,DPhi,fTrigBuffer[fTindex][4]};
 917                         fhnMixedEvents->Fill(triggerEntries);
 918                         }
 919                         fNevents=fNevents+1;
 920                         if(fNevents==10) fTindex=fTindex+1;
 921             }}}
 922
 923                if(fTindex==10&&fNevents==10) fCountAgain=0;
 924
 925                // Copy the triggers from the current event into the buffer.
 926                //again, only if the trigger exists:
 927                if(fCountAgain==0){
 928                 if(trigJet>-1){
 929                 AliAODJet* jetT = (AliAODJet*)(fListJets[0]->At(trigJet));                      AliVParticle *partT = (AliVParticle*)ParticleList.At(trigBBTrack);
 930                 fTrigBuffer[fTrigBufferIndex][0] = primVtx->GetZ();
 931                 fTrigBuffer[fTrigBufferIndex][1] = jetT->Phi();
 932                 fTrigBuffer[fTrigBufferIndex][2] = jetT->Eta();
 933                 fTrigBuffer[fTrigBufferIndex][3] = jetT->Pt()-rho*jetT->EffectiveAreaCharged();
 934                 fTrigBuffer[fTrigBufferIndex][4] = partT->Pt();
 935                 fTrigBuffer[fTrigBufferIndex][5] = centValue;
 936                 fTrigBufferIndex++;
 937                 if(fTrigBufferIndex==9) {fTrigBufferIndex=0;
 938                                          fCountAgain=1;}
 939                 }
 940                }
 941
 942           }
 943
 944           /////////////////////////////////////////////////////////////////////////////
 945           ////////////////////// Rongrong's analysis //////////////////////////////////
 946           if(fRunAnaAzimuthalCorrelation)
 947             {
 948               fhTTPt->Fill(centValue,partback->Pt());
 949               for(Int_t ij=0; ij<fListJets[0]->GetEntries(); ij++)
 950                 {
 951                   AliAODJet* jet = (AliAODJet*)(fListJets[0]->At(ij));
 952                   Double_t jetPt   = jet->Pt();
 953                   Double_t jetEta  = jet->Eta();
 954                   Double_t jetPhi  = jet->Phi();
 955                   if(jetPt==0) continue;
 956                   if((jetEta<fJetEtaMin)||(jetEta>fJetEtaMax)) continue;
 957                   Double_t jetArea = jet->EffectiveAreaCharged();
 958                   Double_t jetPtCorr=jetPt-rho*jetArea;
 959                   Double_t dPhi=jetPhi-partback->Phi();
 960                   if(dPhi>2*TMath::Pi()) dPhi -= 2*TMath::Pi();
 961                   if(dPhi<-2*TMath::Pi()) dPhi += 2*TMath::Pi();
 962                   if(dPhi<-0.5*TMath::Pi()) dPhi += 2*TMath::Pi();
 963                   if(dPhi>1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
 964
 965                   Double_t fill[] = {partback->Pt(),jetPtCorr,dPhi,jetArea,centValue};
 966                   fHJetPhiCorr->Fill(fill);
 967                 }
 968             }
 969           /////////////////////////////////////////////////////////////////////////////
 970           /////////////////////////////////////////////////////////////////////////////
 971
 972
 973      //////////////////ANGULAR STRUCTURE//////////////////////////////////////
 974
 975      //tracks up to R=0.8 distant from the jet axis
 976    //   if(fAngStructCloseTracks==1){
 977    //    TList CloseTrackList;
 978    //    Int_t nn=GetListOfTracksCloseToJet(&CloseTrackList,jetbig);
 979    //    Double_t difR=0.04;
 980    //    for(Int_t l=0;l<15;l++){
 981    //    Double_t rr=l*0.1+0.1;
 982    //     for(int it = 0;it<nn;++it){
 983    //         AliVParticle *part1 = (AliVParticle*)CloseTrackList.At(it);
 984    //         for(int itu=it+1;itu<CloseTrackList.GetEntries();itu++){
 985    //         AliVParticle *part2 = (AliVParticle*)CloseTrackList.At(itu);
 986    //         Double_t ptm=part1->Pt();
 987    //         Double_t ptn=part2->Pt();
 988    //         Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
 989    //                    Rnm=TMath::Sqrt(Rnm);
 990    //                    Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
 991    //                    Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
 992    //                    if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
 993    //                                                   down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
 994    //                    if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
 995    //                                                   down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
 996    //                    if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
 997    //                                                    down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
 998    //                    if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
 999    //                   down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}
1000    //   }
1001
1002    //   //only jet constituents
1003    //    if(fAngStructCloseTracks==2){
1004
1005    //    Double_t difR=0.04;
1006    //    for(Int_t l=0;l<15;l++){
1007    //    Double_t rr=l*0.1+0.1;
1008
1009
1010    //    AliAODTrack* part1;
1011    //    AliAODTrack* part2;
1012
1013    //        TRefArray *genTrackListb = jetbig->GetRefTracks();
1014    //        Int_t nTracksGenJetb = genTrackListb->GetEntriesFast();
1015
1016
1017
1018    //        for(Int_t it=0; it<nTracksGenJetb; ++it){
1019    //           part1 = (AliAODTrack*)(genTrackListb->At(it));
1020    //         for(Int_t itu=0; itu<nTracksGenJetb; ++itu){
1021    //           part2 = (AliAODTrack*)(genTrackListb->At(itu));
1022    //         Double_t ptm=part1->Pt();
1023    //         Double_t ptn=part2->Pt();
1024    //         Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
1025    //                    Rnm=TMath::Sqrt(Rnm);
1026    //                    Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
1027    //                    Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
1028    //                    if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
1029    //                                                   down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
1030    //                    if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
1031    //                                                   down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
1032    //                    if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
1033    //                                                    down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
1034    //                    if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
1035    //                   down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}}
1036    // }
1037      // //end loop over R=0.4 jets
1038      // if(fAngStructCloseTracks>0){
1039      // for(Int_t l=0;l<15;l++){
1040      // Double_t rr=l*0.1+0.1;
1041      //    if(down1[l]!=0){
1042      //         if(centValue<10.)fh2AngStructpt1C10->Fill(rr,rr*up1[l]/down1[l]);
1043      //    if(centValue>20. && centValue<40.) fh2AngStructpt1C20->Fill(rr,rr*up1[l]/down1[l]);
1044      //    if(centValue>30. && centValue<60.) fh2AngStructpt1C30->Fill(rr,rr*up1[l]/down1[l]);
1045      //    if(centValue>60.) fh2AngStructpt1C60->Fill(rr,rr*up1[l]/down1[l]);}
1046      //    if(down2[l]!=0){
1047      //         if(centValue<10.) fh2AngStructpt2C10->Fill(rr,rr*up2[l]/down2[l]);
1048      //    if(centValue>20. && centValue<40.) fh2AngStructpt2C20->Fill(rr,rr*up2[l]/down2[l]);
1049      //    if(centValue>30. && centValue<60.) fh2AngStructpt2C30->Fill(rr,rr*up2[l]/down2[l]);
1050      //    if(centValue>60.) fh2AngStructpt2C60->Fill(rr,rr*up2[l]/down2[l]);}
1051      //    if(down3[l]!=0){
1052      //         if(centValue<10.) fh2AngStructpt3C10->Fill(rr,rr*up3[l]/down3[l]);
1053      //    if(centValue>20. && centValue<40.) fh2AngStructpt3C20->Fill(rr,rr*up3[l]/down3[l]);
1054      //    if(centValue>30. && centValue<60.) fh2AngStructpt3C30->Fill(rr,rr*up3[l]/down3[l]);
1055      //    if(centValue>60.) fh2AngStructpt3C60->Fill(rr,rr*up3[l]/down3[l]);}
1056      //    if(down4[l]!=0){
1057      //         if(centValue<10.) fh2AngStructpt4C10->Fill(rr,rr*up4[l]/down4[l]);
1058      //    if(centValue>20. && centValue<40.) fh2AngStructpt4C20->Fill(rr,rr*up4[l]/down4[l]);
1059      //    if(centValue>30. && centValue<60.) fh2AngStructpt4C30->Fill(rr,rr*up4[l]/down4[l]);
1060      //    if(centValue>60.) fh2AngStructpt4C60->Fill(rr,rr*up4[l]/down4[l]);}}}
1061
1062
1063
1064
1065
1066
1067
1068    PostData(1, fOutputList);
1069 }
1070
1071 void AliAnalysisTaskJetCore::Terminate(const Option_t *)
1072 {
1073    // Draw result to the screen
1074    // Called once at the end of the query
1075
1076    if (!GetOutputData(1))
1077    return;
1078 }
1079
1080
1081
1082
1083
1084
1085 Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
1086
1087       Int_t iCount = 0;
1088       AliAODEvent *aod = 0;
1089
1090      if(!fESD)aod = fAODIn;
1091      else aod = fAODOut;
1092
1093      if(!aod)return 0;
1094
1095      Int_t index=-1;
1096      Double_t ptmax=-10;
1097
1098
1099
1100      for(int it = 0;it < aod->GetNumberOfTracks();++it){
1101       AliAODTrack *tr = aod->GetTrack(it);
1102       Bool_t bGood = false;
1103       if(fFilterType == 0)bGood = true;
1104       else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1105       else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1106      if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1107       if(bGood==false) continue;
1108       if(TMath::Abs(tr->Eta())>0.9)continue;
1109       if(tr->Pt()<0.15)continue;
1110       list->Add(tr);
1111       iCount++;
1112       if(fFilterType==2 && fFilterMaskBestPt>0){// only set the trigger track index for good quality tracks
1113         if(tr->TestFilterBit(fFilterMaskBestPt)){
1114           if(tr->Pt()>ptmax){
1115             ptmax=tr->Pt();
1116             index=iCount-1;
1117           }
1118         }
1119       }
1120       else{
1121         if(tr->Pt()>ptmax){
1122           ptmax=tr->Pt();
1123           index=iCount-1;
1124         }
1125       }
1126      }
1127
1128
1129     // else if (type == kTrackAODMCCharged) {
1130     // TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1131     // if(!tca)return iCount;
1132     // for(int it = 0;it < tca->GetEntriesFast();++it){
1133     //   AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1134     //   if(!part)continue;
1135     //   if(part->Pt()<0.15)continue;
1136     //   if(!part->IsPhysicalPrimary())continue;
1137     //   if(part->Charge()==0)continue;
1138     //   if(TMath::Abs(part->Eta())>0.9)continue;
1139     //   list->Add(part);
1140     //   iCount++;
1141     //   if(part->Pt()>ptmax){ ptmax=part->Pt();
1142     //  index=iCount-1;}}}
1143       return index;
1144
1145 }
1146
1147    Int_t  AliAnalysisTaskJetCore::GetHardestTrackBackToJet(AliAODJet *jetbig){
1148
1149     AliAODEvent *aod = 0;
1150     if(!fESD)aod = fAODIn;
1151     else aod = fAODOut;
1152     Int_t index=-1;
1153     Double_t ptmax=-10;
1154     Double_t dphi=0;
1155     Double_t dif=0;
1156     Int_t iCount=0;
1157     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1158       AliAODTrack *tr = aod->GetTrack(it);
1159       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1160       if(TMath::Abs(tr->Eta())>0.9)continue;
1161       if(tr->Pt()<0.15)continue;
1162       iCount=iCount+1;
1163       dphi=RelativePhi(tr->Phi(),jetbig->Phi());
1164       if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;
1165       if(tr->Pt()>ptmax){ ptmax=tr->Pt();
1166       index=iCount-1;
1167       dif=dphi;  }}
1168
1169       return index;
1170
1171    }
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181  Int_t  AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
1182
1183     Int_t iCount = 0;
1184       AliAODEvent *aod = 0;
1185      if(!fESD)aod = fAODIn;
1186      else aod = fAODOut;
1187
1188       for(int it = 0;it < aod->GetNumberOfTracks();++it){
1189       AliAODTrack *tr = aod->GetTrack(it);
1190       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1191       if(TMath::Abs(tr->Eta())>0.9)continue;
1192       if(tr->Pt()<0.15)continue;
1193       Double_t disR=jetbig->DeltaR(tr);
1194       if(disR>0.8)  continue;
1195       list->Add(tr);
1196       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
1197       iCount++;
1198     }
1199
1200    list->Sort();
1201    return iCount;
1202
1203 }
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215 Int_t AliAnalysisTaskJetCore::GetNInputTracks()
1216 {
1217
1218    Int_t nInputTracks = 0;
1219      AliAODEvent *aod = 0;
1220      if(!fESD)aod = fAODIn;
1221      else aod = fAODOut;
1222    TString jbname(fJetBranchName[1]);
1223    //needs complete event, use jets without background subtraction
1224    for(Int_t i=1; i<=3; ++i){
1225       if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
1226    }
1227    // use only HI event
1228    if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
1229    if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
1230
1231    if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
1232    TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(aod->FindListObject(jbname.Data()));
1233    if(!tmpAODjets){
1234       Printf("Jet branch %s not found", jbname.Data());
1235       Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
1236       return -1;
1237    }
1238
1239    for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
1240       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
1241       if(!jet) continue;
1242       TRefArray *trackList = jet->GetRefTracks();
1243       Int_t nTracks = trackList->GetEntriesFast();
1244       nInputTracks += nTracks;
1245       if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
1246    }
1247    if(fDebug) Printf("---> input tracks: %d", nInputTracks);
1248
1249    return nInputTracks;
1250 }
1251
1252
1253
1254 Double_t AliAnalysisTaskJetCore::RelativePhi(Double_t mphi,Double_t vphi){
1255
1256   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1257   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1258   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1259   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1260   double dphi = mphi-vphi;
1261   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1262   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1263   return dphi;//dphi in [-Pi, Pi]
1264 }
1265
1266 Int_t AliAnalysisTaskJetCore::GetPhiBin(Double_t phi)
1267 {
1268     Int_t phibin=-1;
1269     if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1270     Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1271     phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1272     if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1273     return phibin;
1274 }
1275
1276
1277
1278
1279 THnSparse* AliAnalysisTaskJetCore::NewTHnSparseF(const char* name, UInt_t entries)
1280 {
1281    // generate new THnSparseF, axes are defined in GetDimParams()
1282
1283    Int_t count = 0;
1284    UInt_t tmp = entries;
1285    while(tmp!=0){
1286       count++;
1287       tmp = tmp &~ -tmp;  // clear lowest bit
1288    }
1289
1290    TString hnTitle(name);
1291    const Int_t dim = count;
1292    Int_t nbins[dim];
1293    Double_t xmin[dim];
1294    Double_t xmax[dim];
1295
1296    Int_t i=0;
1297    Int_t c=0;
1298    while(c<dim && i<32){
1299       if(entries&(1<<i)){
1300
1301          TString label("");
1302          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1303          hnTitle += Form(";%s",label.Data());
1304          c++;
1305       }
1306
1307       i++;
1308    }
1309    hnTitle += ";";
1310
1311    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1312 }
1313
1314 void AliAnalysisTaskJetCore::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1315 {
1316    // stores label and binning of axis for THnSparse
1317
1318    const Double_t pi = TMath::Pi();
1319
1320    switch(iEntry){
1321
1322    case 0:
1323       label = "V0 centrality (%)";
1324
1325          nbins = 10;
1326          xmin = 0.;
1327          xmax = 100.;
1328          break;
1329
1330
1331    case 1:
1332       label = "corrected jet pt";
1333          nbins = 20;
1334          xmin = 0.;
1335          xmax = 200.;
1336           break;
1337
1338
1339    case 2:
1340       label = "track pT";
1341
1342          nbins = 9;
1343          xmin = 0.;
1344          xmax = 150;
1345          break;
1346
1347
1348     case 3:
1349       label = "deltaR";
1350       nbins = 15;
1351       xmin = 0.;
1352       xmax = 1.5;
1353       break;
1354
1355
1356
1357    case 4:
1358       label = "deltaEta";
1359       nbins = 8;
1360       xmin = -1.6;
1361       xmax = 1.6;
1362       break;
1363
1364
1365   case 5:
1366       label = "deltaPhi";
1367       nbins = 90;
1368       xmin = -0.5*pi;
1369       xmax = 1.5*pi;
1370       break;
1371
1372
1373
1374     case 6:
1375       label = "leading track";
1376       nbins = 13;
1377       xmin = 0;
1378       xmax = 50;
1379       break;
1380
1381      case 7:
1382
1383       label = "trigger track";
1384       nbins =10;
1385       xmin = 0;
1386       xmax = 50;
1387       break;
1388
1389
1390
1391
1392
1393
1394
1395
1396    }
1397
1398 }
1399