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
 708    fh2Ntriggers->Fill(centValue,partback->Pt());
 709    Double_t phiBinT = RelativePhi(partback->Phi(),fRPAngle);
 710    if(centValue<20.) fh2RPTC20->Fill(TMath::Abs(phiBinT),partback->Pt());
 711    if(centValue<10.) fh2RPTC10->Fill(TMath::Abs(phiBinT),partback->Pt());
 712    Double_t accep=2.*TMath::Pi()*1.8;
 713    Int_t injet4=0;
 714    Int_t injet=0;
 715    if(fSemigoodCorrect){
 716    Double_t disthole=RelativePhi(partback->Phi(),fHolePos);
 717    if(TMath::Abs(disthole)+fHoleWidth>TMath::Pi()-0.6){
 718    PostData(1, fOutputList);
 719    return;}
 720
 721    }
 722
 723    for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
 724            AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
 725            etabig  = jetbig->Eta();
 726            phibig  = jetbig->Phi();
 727            ptbig   = jetbig->Pt();
 728            if(ptbig==0) continue;
 729            Double_t phiBin = RelativePhi(phibig,fRPAngle);
 730            areabig = jetbig->EffectiveAreaCharged();
 731            Double_t ptcorr=ptbig-rho*areabig;
 732            if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
 733            if(areabig>=0.07) injet=injet+1;
 734            if(areabig>=0.4) injet4=injet4+1;
 735            Double_t dphi=RelativePhi(partback->Phi(),phibig);
 736
 737            if(fFlagEtaBkg==1){
 738            Double_t etadif= partback->Eta()-etabig;
 739            if(TMath::Abs(etadif)<=0.5){
 740
 741            if(centValue<20.) fh3JetTrackC20->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
 742            if(centValue>30. && centValue<60.) fh3JetTrackC3060->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));}}
 743            if(fFlagEtaBkg==0){
 744            if(centValue<20.) fh3JetTrackC20->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));
 745            if(centValue>30. && centValue<60.) fh3JetTrackC3060->Fill(partback->Pt(),ptcorr,TMath::Abs(dphi));}
 746
 747
 748            if(fFlagJetHadron==0){
 749            if(fFlagPhiBkg!=0) if((TMath::Abs(dphi)<TMath::Pi()/2.-0.1)||(TMath::Abs(dphi)>TMath::Pi()/2.+0.1)) continue;
 750            if(fFlagPhiBkg==0) if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;}
 751
 752            if(fFlagJetHadron!=0) if(TMath::Abs(dphi)>0.4) continue;
 753
 754
 755            if(centValue<10.) fh2RPJetsC10->Fill(TMath::Abs(phiBin), ptcorr);
 756            if(centValue<20.) fh2RPJetsC20->Fill(TMath::Abs(phiBin), ptcorr);
 757                    Double_t dismin=100.;
 758                    Double_t ptmax=-10.;
 759                    Int_t index1=-1;
 760                    Int_t index2=-1;
 761
 762                    Float_t phitt=partback->Phi();
 763                    if(phitt<0)phitt+=TMath::Pi()*2.;
 764                    Int_t phiBintt = GetPhiBin(phitt-fRPAngle);
 765
 766                    Double_t fillspec[] = {centValue,jetbig->EffectiveAreaCharged(),ptcorr,partback->Pt(),phiBintt};
 767                   fHJetSpec->Fill(fillspec);
 768
 769
 770
 771                    if(ptcorr<=0) continue;
 772
 773                        AliAODTrack* leadtrack=0;
 774                        Int_t ippt=0;
 775                        Double_t ppt=-10;
 776                        if(fFlagJetHadron==0){
 777                        TRefArray *genTrackList = jetbig->GetRefTracks();
 778                        Int_t nTracksGenJet = genTrackList->GetEntriesFast();
 779                        AliAODTrack* genTrack;
 780                        for(Int_t ir=0; ir<nTracksGenJet; ++ir){
 781                        genTrack = (AliAODTrack*)(genTrackList->At(ir));
 782                        if(genTrack->Pt()>ppt){ppt=genTrack->Pt();
 783                        ippt=ir;}}
 784                        leadtrack=(AliAODTrack*)(genTrackList->At(ippt));
 785                        if(!leadtrack) continue;}
 786
 787                        AliVParticle* leadtrackb=0;
 788                        if(fFlagJetHadron!=0){
 789                          Int_t nTb = GetHardestTrackBackToJet(jetbig);
 790                          leadtrackb = (AliVParticle*)ParticleList.At(nTb);
 791                          if(!leadtrackb) continue;
 792                        }
 793
 794
 795
 796
 797
 798                         //store one trigger info
 799                         if(iCount==0){
 800                           trigJet=i;
 801                           trigBBTrack=nT;
 802                           trigInTrack=ippt;
 803                           iCount=iCount+1;}
 804
 805
 806                   if(fCheckMethods){
 807                   for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
 808                   AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
 809                   etasmall  = jetsmall->Eta();
 810                   phismall = jetsmall->Phi();
 811                   ptsmall   = jetsmall->Pt();
 812                   areasmall = jetsmall->EffectiveAreaCharged();
 813                   Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
 814                   tmpDeltaR=TMath::Sqrt(tmpDeltaR);
 815                      //Fraction in the jet core
 816                     if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;
 817                     index2=j;}
 818                     if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
 819                       index1=j;}} //en of loop over R=0.2 jets
 820                 //method1:most concentric jet=core
 821                 if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));
 822                   if(centValue<10) fh2JetCoreMethod1C10->Fill(ptcorr,jetmethod1->Pt()/ptbig);
 823                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod1C20->Fill(ptcorr,jetmethod1->Pt()/ptbig);
 824                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod1C30->Fill(ptcorr,jetmethod1->Pt()/ptbig);
 825                   if(centValue>60) fh2JetCoreMethod1C60->Fill(ptcorr,jetmethod1->Pt()/ptbig); }
 826                 //method2:hardest contained jet=core
 827                 if(index2!=-1){
 828                   AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
 829                   if(centValue<10) fh2JetCoreMethod2C10->Fill(ptcorr,jetmethod2->Pt()/ptbig);
 830                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod2C20->Fill(ptcorr,jetmethod2->Pt()/ptbig);
 831                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptcorr,jetmethod2->Pt()/ptbig);
 832                   if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }}
 833                   if(centValue<10) fh2Ntriggers2C10->Fill(leadtrack->Pt(),partback->Pt());
 834                   if(centValue<20) fh2Ntriggers2C20->Fill(leadtrack->Pt(),partback->Pt());
 835          if(fDoEventMixing==0 && fFlagOnlyRecoil==0){
 836          for(int it = 0;it<ParticleList.GetEntries();++it){
 837           AliVParticle *part = (AliVParticle*)ParticleList.At(it);
 838           Double_t deltaR = jetbig->DeltaR(part);
 839           Double_t deltaEta = etabig-part->Eta();
 840
 841           Double_t deltaPhi=phibig-part->Phi();
 842           if(deltaPhi<-0.5*TMath::Pi()) deltaPhi+=2.*TMath::Pi();
 843           if(deltaPhi>3./2.*TMath::Pi()) deltaPhi-=2.*TMath::Pi();
 844           Double_t pTcont=0;
 845           if(fFlagJetHadron==0) pTcont=leadtrack->Pt();
 846           if(fFlagJetHadron!=0) pTcont=leadtrackb->Pt();
 847            Double_t jetEntries[8] = {centValue,ptcorr,part->Pt(),deltaR,deltaEta,deltaPhi,pTcont,partback->Pt()};
 848            fhnDeltaR->Fill(jetEntries);}
 849
 850
 851           }
 852
 853
 854           //end of track loop, we only do it if EM is switched off
 855
 856
 857
 858
 859
 860
 861
 862
 863
 864    }
 865    if(injet>0) fh3JetDensity->Fill(ParticleList.GetEntries(),injet/accep,partback->Pt());
 866    if(injet4>0)fh3JetDensityA4->Fill(ParticleList.GetEntries(),injet4/accep,partback->Pt());
 867           //end of jet loop
 868
 869    //}
 870
 871
 872           if(fDoEventMixing>0){
 873             //check before if the trigger exists
 874             // fTrigBuffer[i][0] = zvtx
 875             // fTrigBuffer[i][1] = phi
 876             // fTrigBuffer[i][2] = eta
 877             // fTrigBuffer[i][3] = pt_jet
 878             // fTrigBuffer[i][4] = pt_trig
 879             // fTrigBuffer[i][5]= centrality
 880             if(fTindex==10) fTindex=0;
 881             if(fTrigBuffer[fTindex][3]>0){
 882             if (TMath::Abs(fTrigBuffer[fTindex][0]-primVtx->GetZ()<2.)){
 883             if (TMath::Abs(fTrigBuffer[fTindex][5]-centValue<5)){
 884
 885                         for(int it = 0;it<nT;++it){
 886                         AliVParticle *part = (AliVParticle*)ParticleList.At(it);
 887                         Double_t DPhi = fTrigBuffer[fTindex][1] - part->Phi();
 888                         Double_t DEta = fTrigBuffer[fTindex][2] - part->Eta();
 889                         Double_t DR=TMath::Sqrt(DPhi*DPhi+DEta*DEta);
 890                         if(DPhi<-0.5*TMath::Pi()) DPhi+=2.*TMath::Pi();
 891                         if(DPhi>3./2.*TMath::Pi()) DPhi-=2.*TMath::Pi();
 892                         Double_t triggerEntries[7] = {centValue,fTrigBuffer[fTindex][3],part->Pt(),DR,DEta,DPhi,fTrigBuffer[fTindex][4]};
 893                         fhnMixedEvents->Fill(triggerEntries);
 894                         }
 895                         fNevents=fNevents+1;
 896                         if(fNevents==10) fTindex=fTindex+1;
 897             }}}
 898
 899                if(fTindex==10&&fNevents==10) fCountAgain=0;
 900
 901                // Copy the triggers from the current event into the buffer.
 902                //again, only if the trigger exists:
 903                if(fCountAgain==0){
 904                 if(trigJet>-1){
 905                 AliAODJet* jetT = (AliAODJet*)(fListJets[0]->At(trigJet));                      AliVParticle *partT = (AliVParticle*)ParticleList.At(trigBBTrack);
 906                 fTrigBuffer[fTrigBufferIndex][0] = primVtx->GetZ();
 907                 fTrigBuffer[fTrigBufferIndex][1] = jetT->Phi();
 908                 fTrigBuffer[fTrigBufferIndex][2] = jetT->Eta();
 909                 fTrigBuffer[fTrigBufferIndex][3] = jetT->Pt()-rho*jetT->EffectiveAreaCharged();
 910                 fTrigBuffer[fTrigBufferIndex][4] = partT->Pt();
 911                 fTrigBuffer[fTrigBufferIndex][5] = centValue;
 912                 fTrigBufferIndex++;
 913                 if(fTrigBufferIndex==9) {fTrigBufferIndex=0;
 914                                          fCountAgain=1;}
 915                 }
 916                }
 917
 918           }
 919
 920           /////////////////////////////////////////////////////////////////////////////
 921           ////////////////////// Rongrong's analysis //////////////////////////////////
 922           if(fRunAnaAzimuthalCorrelation)
 923             {
 924               fhTTPt->Fill(centValue,partback->Pt());
 925               for(Int_t ij=0; ij<fListJets[0]->GetEntries(); ij++)
 926                 {
 927                   AliAODJet* jet = (AliAODJet*)(fListJets[0]->At(ij));
 928                   Double_t jetPt   = jet->Pt();
 929                   Double_t jetEta  = jet->Eta();
 930                   Double_t jetPhi  = jet->Phi();
 931                   if(jetPt==0) continue;
 932                   if((jetEta<fJetEtaMin)||(jetEta>fJetEtaMax)) continue;
 933                   Double_t jetArea = jet->EffectiveAreaCharged();
 934                   Double_t jetPtCorr=jetPt-rho*jetArea;
 935                   Double_t dPhi=jetPhi-partback->Phi();
 936                   if(dPhi>2*TMath::Pi()) dPhi -= 2*TMath::Pi();
 937                   if(dPhi<-2*TMath::Pi()) dPhi += 2*TMath::Pi();
 938                   if(dPhi<-0.5*TMath::Pi()) dPhi += 2*TMath::Pi();
 939                   if(dPhi>1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
 940
 941                   Double_t fill[] = {partback->Pt(),jetPtCorr,dPhi,jetArea,centValue};
 942                   fHJetPhiCorr->Fill(fill);
 943                 }
 944             }
 945           /////////////////////////////////////////////////////////////////////////////
 946           /////////////////////////////////////////////////////////////////////////////
 947
 948
 949      //////////////////ANGULAR STRUCTURE//////////////////////////////////////
 950
 951      //tracks up to R=0.8 distant from the jet axis
 952    //   if(fAngStructCloseTracks==1){
 953    //    TList CloseTrackList;
 954    //    Int_t nn=GetListOfTracksCloseToJet(&CloseTrackList,jetbig);
 955    //    Double_t difR=0.04;
 956    //    for(Int_t l=0;l<15;l++){
 957    //    Double_t rr=l*0.1+0.1;
 958    //     for(int it = 0;it<nn;++it){
 959    //         AliVParticle *part1 = (AliVParticle*)CloseTrackList.At(it);
 960    //         for(int itu=it+1;itu<CloseTrackList.GetEntries();itu++){
 961    //         AliVParticle *part2 = (AliVParticle*)CloseTrackList.At(itu);
 962    //         Double_t ptm=part1->Pt();
 963    //         Double_t ptn=part2->Pt();
 964    //         Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
 965    //                    Rnm=TMath::Sqrt(Rnm);
 966    //                    Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
 967    //                    Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
 968    //                    if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
 969    //                                                   down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
 970    //                    if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
 971    //                                                   down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
 972    //                    if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
 973    //                                                    down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
 974    //                    if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
 975    //                   down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}
 976    //   }
 977
 978    //   //only jet constituents
 979    //    if(fAngStructCloseTracks==2){
 980
 981    //    Double_t difR=0.04;
 982    //    for(Int_t l=0;l<15;l++){
 983    //    Double_t rr=l*0.1+0.1;
 984
 985
 986    //    AliAODTrack* part1;
 987    //    AliAODTrack* part2;
 988
 989    //        TRefArray *genTrackListb = jetbig->GetRefTracks();
 990    //        Int_t nTracksGenJetb = genTrackListb->GetEntriesFast();
 991
 992
 993
 994    //        for(Int_t it=0; it<nTracksGenJetb; ++it){
 995    //           part1 = (AliAODTrack*)(genTrackListb->At(it));
 996    //         for(Int_t itu=0; itu<nTracksGenJetb; ++itu){
 997    //           part2 = (AliAODTrack*)(genTrackListb->At(itu));
 998    //         Double_t ptm=part1->Pt();
 999    //         Double_t ptn=part2->Pt();
1000    //         Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
1001    //                    Rnm=TMath::Sqrt(Rnm);
1002    //                    Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
1003    //                    Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
1004    //                    if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
1005    //                                                   down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
1006    //                    if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
1007    //                                                   down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
1008    //                    if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
1009    //                                                    down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
1010    //                    if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
1011    //                   down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}}
1012    // }
1013      // //end loop over R=0.4 jets
1014      // if(fAngStructCloseTracks>0){
1015      // for(Int_t l=0;l<15;l++){
1016      // Double_t rr=l*0.1+0.1;
1017      //    if(down1[l]!=0){
1018      //         if(centValue<10.)fh2AngStructpt1C10->Fill(rr,rr*up1[l]/down1[l]);
1019      //    if(centValue>20. && centValue<40.) fh2AngStructpt1C20->Fill(rr,rr*up1[l]/down1[l]);
1020      //    if(centValue>30. && centValue<60.) fh2AngStructpt1C30->Fill(rr,rr*up1[l]/down1[l]);
1021      //    if(centValue>60.) fh2AngStructpt1C60->Fill(rr,rr*up1[l]/down1[l]);}
1022      //    if(down2[l]!=0){
1023      //         if(centValue<10.) fh2AngStructpt2C10->Fill(rr,rr*up2[l]/down2[l]);
1024      //    if(centValue>20. && centValue<40.) fh2AngStructpt2C20->Fill(rr,rr*up2[l]/down2[l]);
1025      //    if(centValue>30. && centValue<60.) fh2AngStructpt2C30->Fill(rr,rr*up2[l]/down2[l]);
1026      //    if(centValue>60.) fh2AngStructpt2C60->Fill(rr,rr*up2[l]/down2[l]);}
1027      //    if(down3[l]!=0){
1028      //         if(centValue<10.) fh2AngStructpt3C10->Fill(rr,rr*up3[l]/down3[l]);
1029      //    if(centValue>20. && centValue<40.) fh2AngStructpt3C20->Fill(rr,rr*up3[l]/down3[l]);
1030      //    if(centValue>30. && centValue<60.) fh2AngStructpt3C30->Fill(rr,rr*up3[l]/down3[l]);
1031      //    if(centValue>60.) fh2AngStructpt3C60->Fill(rr,rr*up3[l]/down3[l]);}
1032      //    if(down4[l]!=0){
1033      //         if(centValue<10.) fh2AngStructpt4C10->Fill(rr,rr*up4[l]/down4[l]);
1034      //    if(centValue>20. && centValue<40.) fh2AngStructpt4C20->Fill(rr,rr*up4[l]/down4[l]);
1035      //    if(centValue>30. && centValue<60.) fh2AngStructpt4C30->Fill(rr,rr*up4[l]/down4[l]);
1036      //    if(centValue>60.) fh2AngStructpt4C60->Fill(rr,rr*up4[l]/down4[l]);}}}
1037
1038
1039
1040
1041
1042
1043
1044    PostData(1, fOutputList);
1045 }
1046
1047 void AliAnalysisTaskJetCore::Terminate(const Option_t *)
1048 {
1049    // Draw result to the screen
1050    // Called once at the end of the query
1051
1052    if (!GetOutputData(1))
1053    return;
1054 }
1055
1056
1057
1058
1059
1060
1061 Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
1062
1063       Int_t iCount = 0;
1064       AliAODEvent *aod = 0;
1065
1066      if(!fESD)aod = fAODIn;
1067      else aod = fAODOut;
1068       Int_t index=-1;
1069      Double_t ptmax=-10;
1070
1071
1072
1073      for(int it = 0;it < aod->GetNumberOfTracks();++it){
1074       AliAODTrack *tr = aod->GetTrack(it);
1075       Bool_t bGood = false;
1076       if(fFilterType == 0)bGood = true;
1077       else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1078       else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1079      if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1080       if(bGood==false) continue;
1081       if(TMath::Abs(tr->Eta())>0.9)continue;
1082       if(tr->Pt()<0.15)continue;
1083       list->Add(tr);
1084       iCount++;
1085       if(fFilterType==2 && fFilterMaskBestPt>0){// only set the trigger track index for good quality tracks
1086         if(tr->TestFilterBit(fFilterMaskBestPt)){
1087           if(tr->Pt()>ptmax){
1088             ptmax=tr->Pt();
1089             index=iCount-1;
1090           }
1091         }
1092       }
1093       else{
1094         if(tr->Pt()>ptmax){
1095           ptmax=tr->Pt();
1096           index=iCount-1;
1097         }
1098       }
1099      }
1100
1101
1102     // else if (type == kTrackAODMCCharged) {
1103     // TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1104     // if(!tca)return iCount;
1105     // for(int it = 0;it < tca->GetEntriesFast();++it){
1106     //   AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1107     //   if(!part)continue;
1108     //   if(part->Pt()<0.15)continue;
1109     //   if(!part->IsPhysicalPrimary())continue;
1110     //   if(part->Charge()==0)continue;
1111     //   if(TMath::Abs(part->Eta())>0.9)continue;
1112     //   list->Add(part);
1113     //   iCount++;
1114     //   if(part->Pt()>ptmax){ ptmax=part->Pt();
1115     //  index=iCount-1;}}}
1116       return index;
1117
1118 }
1119
1120    Int_t  AliAnalysisTaskJetCore::GetHardestTrackBackToJet(AliAODJet *jetbig){
1121
1122     AliAODEvent *aod = 0;
1123     if(!fESD)aod = fAODIn;
1124     else aod = fAODOut;
1125     Int_t index=-1;
1126     Double_t ptmax=-10;
1127     Double_t dphi=0;
1128     Double_t dif=0;
1129     Int_t iCount=0;
1130     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1131       AliAODTrack *tr = aod->GetTrack(it);
1132       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1133       if(TMath::Abs(tr->Eta())>0.9)continue;
1134       if(tr->Pt()<0.15)continue;
1135       iCount=iCount+1;
1136       dphi=RelativePhi(tr->Phi(),jetbig->Phi());
1137       if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;
1138       if(tr->Pt()>ptmax){ ptmax=tr->Pt();
1139       index=iCount-1;
1140       dif=dphi;  }}
1141
1142       return index;
1143
1144    }
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154  Int_t  AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
1155
1156     Int_t iCount = 0;
1157       AliAODEvent *aod = 0;
1158      if(!fESD)aod = fAODIn;
1159      else aod = fAODOut;
1160
1161       for(int it = 0;it < aod->GetNumberOfTracks();++it){
1162       AliAODTrack *tr = aod->GetTrack(it);
1163       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1164       if(TMath::Abs(tr->Eta())>0.9)continue;
1165       if(tr->Pt()<0.15)continue;
1166       Double_t disR=jetbig->DeltaR(tr);
1167       if(disR>0.8)  continue;
1168       list->Add(tr);
1169       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
1170       iCount++;
1171     }
1172
1173    list->Sort();
1174    return iCount;
1175
1176 }
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188 Int_t AliAnalysisTaskJetCore::GetNInputTracks()
1189 {
1190
1191    Int_t nInputTracks = 0;
1192      AliAODEvent *aod = 0;
1193      if(!fESD)aod = fAODIn;
1194      else aod = fAODOut;
1195    TString jbname(fJetBranchName[1]);
1196    //needs complete event, use jets without background subtraction
1197    for(Int_t i=1; i<=3; ++i){
1198       if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
1199    }
1200    // use only HI event
1201    if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
1202    if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
1203
1204    if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
1205    TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(aod->FindListObject(jbname.Data()));
1206    if(!tmpAODjets){
1207       Printf("Jet branch %s not found", jbname.Data());
1208       Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
1209       return -1;
1210    }
1211
1212    for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
1213       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
1214       if(!jet) continue;
1215       TRefArray *trackList = jet->GetRefTracks();
1216       Int_t nTracks = trackList->GetEntriesFast();
1217       nInputTracks += nTracks;
1218       if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
1219    }
1220    if(fDebug) Printf("---> input tracks: %d", nInputTracks);
1221
1222    return nInputTracks;
1223 }
1224
1225
1226
1227 Double_t AliAnalysisTaskJetCore::RelativePhi(Double_t mphi,Double_t vphi){
1228
1229   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1230   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1231   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1232   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1233   double dphi = mphi-vphi;
1234   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1235   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1236   return dphi;//dphi in [-Pi, Pi]
1237 }
1238
1239 Int_t AliAnalysisTaskJetCore::GetPhiBin(Double_t phi)
1240 {
1241     Int_t phibin=-1;
1242     if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1243     Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1244     phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1245     if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1246     return phibin;
1247 }
1248
1249
1250
1251
1252 THnSparse* AliAnalysisTaskJetCore::NewTHnSparseF(const char* name, UInt_t entries)
1253 {
1254    // generate new THnSparseF, axes are defined in GetDimParams()
1255
1256    Int_t count = 0;
1257    UInt_t tmp = entries;
1258    while(tmp!=0){
1259       count++;
1260       tmp = tmp &~ -tmp;  // clear lowest bit
1261    }
1262
1263    TString hnTitle(name);
1264    const Int_t dim = count;
1265    Int_t nbins[dim];
1266    Double_t xmin[dim];
1267    Double_t xmax[dim];
1268
1269    Int_t i=0;
1270    Int_t c=0;
1271    while(c<dim && i<32){
1272       if(entries&(1<<i)){
1273
1274          TString label("");
1275          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1276          hnTitle += Form(";%s",label.Data());
1277          c++;
1278       }
1279
1280       i++;
1281    }
1282    hnTitle += ";";
1283
1284    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1285 }
1286
1287 void AliAnalysisTaskJetCore::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1288 {
1289    // stores label and binning of axis for THnSparse
1290
1291    const Double_t pi = TMath::Pi();
1292
1293    switch(iEntry){
1294
1295    case 0:
1296       label = "V0 centrality (%)";
1297
1298          nbins = 10;
1299          xmin = 0.;
1300          xmax = 100.;
1301          break;
1302
1303
1304    case 1:
1305       label = "corrected jet pt";
1306          nbins = 20;
1307          xmin = 0.;
1308          xmax = 200.;
1309           break;
1310
1311
1312    case 2:
1313       label = "track pT";
1314
1315          nbins = 9;
1316          xmin = 0.;
1317          xmax = 150;
1318          break;
1319
1320
1321     case 3:
1322       label = "deltaR";
1323       nbins = 15;
1324       xmin = 0.;
1325       xmax = 1.5;
1326       break;
1327
1328
1329
1330    case 4:
1331       label = "deltaEta";
1332       nbins = 8;
1333       xmin = -1.6;
1334       xmax = 1.6;
1335       break;
1336
1337
1338   case 5:
1339       label = "deltaPhi";
1340       nbins = 90;
1341       xmin = -0.5*pi;
1342       xmax = 1.5*pi;
1343       break;
1344
1345
1346
1347     case 6:
1348       label = "leading track";
1349       nbins = 13;
1350       xmin = 0;
1351       xmax = 50;
1352       break;
1353
1354      case 7:
1355
1356       label = "trigger track";
1357       nbins =10;
1358       xmin = 0;
1359       xmax = 50;
1360       break;
1361
1362
1363
1364
1365
1366
1367
1368
1369    }
1370
1371 }
1372