PWGJE/AliAnalysisTaskJetCore.cxx

   1 // ******************************************
   2 // This task computes several jet observables like
   3 // the fraction of energy in inner and outer coronnas,
   4 // the distance from track to jet axis and a
   5 // correlation strength distribution of particles inside jets.
   6 // Author: lcunquei@cern.ch
   7 // *******************************************
   8
   9
  10 /**************************************************************************
  11  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  12  *                                                                        *
  13  * Author: The ALICE Off-line Project.                                    *
  14  * Contributors are mentioned in the code where appropriate.              *
  15  *                                                                        *
  16  * Permission to use, copy, modify and distribute this software and its   *
  17  * documentation strictly for non-commercial purposes is hereby granted   *
  18  * without fee, provided that the above copyright notice appears in all   *
  19  * copies and that both the copyright notice and this permission notice   *
  20  * appear in the supporting documentation. The authors make no claims     *
  21  * about the suitability of this software for any purpose. It is          *
  22  * provided "as is" without express or implied warranty.                  *
  23  **************************************************************************/
  24
  25
  26 #include "TChain.h"
  27 #include "TTree.h"
  28 #include "TMath.h"
  29 #include "TH1I.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 "AliAnalysisTaskFastEmbedding.h"
  49
  50 #include "AliAODEvent.h"
  51 #include "AliAODHandler.h"
  52 #include "AliAODJet.h"
  53
  54 #include "AliAnalysisTaskJetCore.h"
  55
  56 ClassImp(AliAnalysisTaskJetCore)
  57
  58 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore() :
  59 AliAnalysisTaskSE(),
  60 fESD(0x0),
  61 fAOD(0x0),
  62 fAODExtension(0x0),
  63 fBackgroundBranch(""),
  64 fNonStdFile(""),
  65 fIsPbPb(kTRUE),
  66 fOfflineTrgMask(AliVEvent::kAny),
  67 fMinContribVtx(1),
  68 fVtxZMin(-8.),
  69 fVtxZMax(8.),
  70 fEvtClassMin(0),
  71 fEvtClassMax(4),
  72 fRadioFrac(0.2),
  73 fMinDist(0.1),
  74 fCentMin(0.),
  75 fCentMax(100.),
  76 fNInputTracksMin(0),
  77 fNInputTracksMax(-1),
  78 fAngStructCloseTracks(0),
  79 fJetEtaMin(-.5),
  80 fJetEtaMax(.5),
  81 fJetPtMin(20.),
  82 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
  83 fJetPtFractionMin(0.5),
  84 fNMatchJets(4),
  85 fMatchMaxDist(0.8),
  86 fKeepJets(kFALSE),
  87 fkNbranches(2),
  88 fkEvtClasses(12),
  89 fOutputList(0x0),
  90 fbEvent(kTRUE),
  91 fHistEvtSelection(0x0),
  92 fHistJetSelection(0x0),
  93 fh2JetSelection(0x0),
  94 fh2JetCoreMethod1C10(0x0),
  95 fh2JetCoreMethod2C10(0x0),
  96 fh2JetCoreMethod3C10(0x0),
  97 fh2JetCoreMethod1C20(0x0),
  98 fh2JetCoreMethod2C20(0x0),
  99 fh2JetCoreMethod3C20(0x0),
 100 fh2JetCoreMethod1C30(0x0),
 101 fh2JetCoreMethod2C30(0x0),
 102 fh2JetCoreMethod3C30(0x0),
 103 fh2JetCoreMethod1C60(0x0),
 104 fh2JetCoreMethod2C60(0x0),
 105 fh2JetCoreMethod3C60(0x0),
 106 fh2SumPtInC10(0x0),
 107 fh2SumPtInC20(0x0),
 108 fh2SumPtInC30(0x0),
 109 fh2SumPtInC60(0x0),
 110 fh2SumPtOutC10(0x0),
 111 fh2SumPtOutC10b(0x0),
 112 fh2SumPtOutC20(0x0),
 113 fh2SumPtOutC30(0x0),
 114 fh2SumPtOutC60(0x0),
 115 fh2SumPtInC10bkg(0x0),
 116 fh2SumPtInC20bkg(0x0),
 117 fh2SumPtInC30bkg(0x0),
 118 fh2SumPtInC60bkg(0x0),
 119 fh2SumPtOutC10bkg(0x0),
 120 fh2SumPtOutC20bkg(0x0),
 121 fh2SumPtOutC30bkg(0x0),
 122 fh2SumPtOutC60bkg(0x0),
 123 fh2DeltaRC10pt1(0x0),
 124 fh2DeltaRC20pt1(0x0),
 125 fh2DeltaRC30pt1(0x0),
 126 fh2DeltaRC60pt1(0x0),
 127 fh2DeltaRC10pt2(0x0),
 128 fh2DeltaRC20pt2(0x0),
 129 fh2DeltaRC30pt2(0x0),
 130 fh2DeltaRC60pt2(0x0),
 131 fh2DeltaRC10pt3(0x0),
 132 fh2DeltaRC20pt3(0x0),
 133 fh2DeltaRC30pt3(0x0),
 134 fh2DeltaRC60pt3(0x0),
 135 fh2DeltaRC10pt4(0x0),
 136 fh2DeltaRC20pt4(0x0),
 137 fh2DeltaRC30pt4(0x0),
 138 fh2DeltaRC60pt4(0x0),
 139 fh2DeltaEtaC10pt1(0x0),
 140 fh2DeltaEtaC20pt1(0x0),
 141 fh2DeltaEtaC30pt1(0x0),
 142 fh2DeltaEtaC60pt1(0x0),
 143 fh2DeltaEtaC10pt2(0x0),
 144 fh2DeltaEtaC20pt2(0x0),
 145 fh2DeltaEtaC30pt2(0x0),
 146 fh2DeltaEtaC60pt2(0x0),
 147 fh2DeltaEtaC10pt3(0x0),
 148 fh2DeltaEtaC20pt3(0x0),
 149 fh2DeltaEtaC30pt3(0x0),
 150 fh2DeltaEtaC60pt3(0x0),
 151 fh2DeltaEtaC10pt4(0x0),
 152 fh2DeltaEtaC20pt4(0x0),
 153 fh2DeltaEtaC30pt4(0x0),
 154 fh2DeltaEtaC60pt4(0x0),
 155 fh2DeltaPhiC10pt1(0x0),
 156 fh2DeltaPhiC20pt1(0x0),
 157 fh2DeltaPhiC30pt1(0x0),
 158 fh2DeltaPhiC60pt1(0x0),
 159 fh2DeltaPhiC10pt2(0x0),
 160 fh2DeltaPhiC20pt2(0x0),
 161 fh2DeltaPhiC30pt2(0x0),
 162 fh2DeltaPhiC60pt2(0x0),
 163 fh2DeltaPhiC10pt3(0x0),
 164 fh2DeltaPhiC20pt3(0x0),
 165 fh2DeltaPhiC30pt3(0x0),
 166 fh2DeltaPhiC60pt3(0x0),
 167 fh2DeltaPhiC10pt4(0x0),
 168 fh2DeltaPhiC20pt4(0x0),
 169 fh2DeltaPhiC30pt4(0x0),
 170 fh2DeltaPhiC60pt4(0x0),
 171 fh2AngStructpt1C10(0x0),
 172 fh2AngStructpt2C10(0x0),
 173 fh2AngStructpt3C10(0x0),
 174 fh2AngStructpt4C10(0x0),
 175 fh2AngStructpt1C20(0x0),
 176 fh2AngStructpt2C20(0x0),
 177 fh2AngStructpt3C20(0x0),
 178 fh2AngStructpt4C20(0x0),
 179 fh2AngStructpt1C30(0x0),
 180 fh2AngStructpt2C30(0x0),
 181 fh2AngStructpt3C30(0x0),
 182 fh2AngStructpt4C30(0x0),
 183 fh2AngStructpt1C60(0x0),
 184 fh2AngStructpt2C60(0x0),
 185 fh2AngStructpt3C60(0x0),
 186 fh2AngStructpt4C60(0x0)
 187 {
 188    // default Constructor
 189
 190    fJetBranchName[0] = "";
 191    fJetBranchName[1] = "";
 192
 193    fListJets[0] = new TList;
 194    fListJets[1] = new TList;
 195 }
 196
 197 AliAnalysisTaskJetCore::AliAnalysisTaskJetCore(const char *name) :
 198 AliAnalysisTaskSE(name),
 199 fESD(0x0),
 200 fAOD(0x0),
 201 fAODExtension(0x0),
 202 fBackgroundBranch(""),
 203 fNonStdFile(""),
 204 fIsPbPb(kTRUE),
 205 fOfflineTrgMask(AliVEvent::kAny),
 206 fMinContribVtx(1),
 207 fVtxZMin(-8.),
 208 fVtxZMax(8.),
 209 fEvtClassMin(0),
 210 fEvtClassMax(4),
 211 fRadioFrac(0.2),
 212 fMinDist(0.1),
 213 fCentMin(0.),
 214 fCentMax(100.),
 215 fNInputTracksMin(0),
 216 fNInputTracksMax(-1),
 217 fAngStructCloseTracks(0),
 218 fJetEtaMin(-.5),
 219 fJetEtaMax(.5),
 220 fJetPtMin(20.),
 221 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
 222 fJetPtFractionMin(0.5),
 223 fNMatchJets(4),
 224 fMatchMaxDist(0.8),
 225 fKeepJets(kFALSE),
 226 fkNbranches(2),
 227 fkEvtClasses(12),
 228 fOutputList(0x0),
 229 fbEvent(kTRUE),
 230 fHistEvtSelection(0x0),
 231 fHistJetSelection(0x0),
 232 fh2JetSelection(0x0),
 233 fh2JetCoreMethod1C10(0x0),
 234 fh2JetCoreMethod2C10(0x0),
 235 fh2JetCoreMethod3C10(0x0),
 236 fh2JetCoreMethod1C20(0x0),
 237 fh2JetCoreMethod2C20(0x0),
 238 fh2JetCoreMethod3C20(0x0),
 239 fh2JetCoreMethod1C30(0x0),
 240 fh2JetCoreMethod2C30(0x0),
 241 fh2JetCoreMethod3C30(0x0),
 242 fh2JetCoreMethod1C60(0x0),
 243 fh2JetCoreMethod2C60(0x0),
 244 fh2JetCoreMethod3C60(0x0),
 245 fh2SumPtInC10(0x0),
 246 fh2SumPtInC20(0x0),
 247 fh2SumPtInC30(0x0),
 248 fh2SumPtInC60(0x0),
 249 fh2SumPtOutC10(0x0),
 250 fh2SumPtOutC10b(0x0),
 251 fh2SumPtOutC20(0x0),
 252 fh2SumPtOutC30(0x0),
 253 fh2SumPtOutC60(0x0),
 254 fh2SumPtInC10bkg(0x0),
 255 fh2SumPtInC20bkg(0x0),
 256 fh2SumPtInC30bkg(0x0),
 257 fh2SumPtInC60bkg(0x0),
 258 fh2SumPtOutC10bkg(0x0),
 259 fh2SumPtOutC20bkg(0x0),
 260 fh2SumPtOutC30bkg(0x0),
 261 fh2SumPtOutC60bkg(0x0),
 262 fh2DeltaRC10pt1(0x0),
 263 fh2DeltaRC20pt1(0x0),
 264 fh2DeltaRC30pt1(0x0),
 265 fh2DeltaRC60pt1(0x0),
 266 fh2DeltaRC10pt2(0x0),
 267 fh2DeltaRC20pt2(0x0),
 268 fh2DeltaRC30pt2(0x0),
 269 fh2DeltaRC60pt2(0x0),
 270 fh2DeltaRC10pt3(0x0),
 271 fh2DeltaRC20pt3(0x0),
 272 fh2DeltaRC30pt3(0x0),
 273 fh2DeltaRC60pt3(0x0),
 274 fh2DeltaRC10pt4(0x0),
 275 fh2DeltaRC20pt4(0x0),
 276 fh2DeltaRC30pt4(0x0),
 277 fh2DeltaRC60pt4(0x0),
 278 fh2DeltaEtaC10pt1(0x0),
 279 fh2DeltaEtaC20pt1(0x0),
 280 fh2DeltaEtaC30pt1(0x0),
 281 fh2DeltaEtaC60pt1(0x0),
 282 fh2DeltaEtaC10pt2(0x0),
 283 fh2DeltaEtaC20pt2(0x0),
 284 fh2DeltaEtaC30pt2(0x0),
 285 fh2DeltaEtaC60pt2(0x0),
 286 fh2DeltaEtaC10pt3(0x0),
 287 fh2DeltaEtaC20pt3(0x0),
 288 fh2DeltaEtaC30pt3(0x0),
 289 fh2DeltaEtaC60pt3(0x0),
 290 fh2DeltaEtaC10pt4(0x0),
 291 fh2DeltaEtaC20pt4(0x0),
 292 fh2DeltaEtaC30pt4(0x0),
 293 fh2DeltaEtaC60pt4(0x0),
 294 fh2DeltaPhiC10pt1(0x0),
 295 fh2DeltaPhiC20pt1(0x0),
 296 fh2DeltaPhiC30pt1(0x0),
 297 fh2DeltaPhiC60pt1(0x0),
 298 fh2DeltaPhiC10pt2(0x0),
 299 fh2DeltaPhiC20pt2(0x0),
 300 fh2DeltaPhiC30pt2(0x0),
 301 fh2DeltaPhiC60pt2(0x0),
 302 fh2DeltaPhiC10pt3(0x0),
 303 fh2DeltaPhiC20pt3(0x0),
 304 fh2DeltaPhiC30pt3(0x0),
 305 fh2DeltaPhiC60pt3(0x0),
 306 fh2DeltaPhiC10pt4(0x0),
 307 fh2DeltaPhiC20pt4(0x0),
 308 fh2DeltaPhiC30pt4(0x0),
 309 fh2DeltaPhiC60pt4(0x0),
 310 fh2AngStructpt1C10(0x0),
 311 fh2AngStructpt2C10(0x0),
 312 fh2AngStructpt3C10(0x0),
 313 fh2AngStructpt4C10(0x0),
 314 fh2AngStructpt1C20(0x0),
 315 fh2AngStructpt2C20(0x0),
 316 fh2AngStructpt3C20(0x0),
 317 fh2AngStructpt4C20(0x0),
 318 fh2AngStructpt1C30(0x0),
 319 fh2AngStructpt2C30(0x0),
 320 fh2AngStructpt3C30(0x0),
 321 fh2AngStructpt4C30(0x0),
 322 fh2AngStructpt1C60(0x0),
 323 fh2AngStructpt2C60(0x0),
 324 fh2AngStructpt3C60(0x0),
 325 fh2AngStructpt4C60(0x0)
 326
 327  {
 328    // Constructor
 329
 330    fJetBranchName[0] = "";
 331    fJetBranchName[1] = "";
 332
 333    fListJets[0] = new TList;
 334    fListJets[1] = new TList;
 335
 336    DefineOutput(1, TList::Class());
 337 }
 338
 339 AliAnalysisTaskJetCore::~AliAnalysisTaskJetCore()
 340 {
 341    delete fListJets[0];
 342    delete fListJets[1];
 343 }
 344
 345 void AliAnalysisTaskJetCore::SetBranchNames(const TString &branch1, const TString &branch2)
 346 {
 347    fJetBranchName[0] = branch1;
 348    fJetBranchName[1] = branch2;
 349 }
 350
 351 void AliAnalysisTaskJetCore::Init()
 352 {
 353
 354    // check for jet branches
 355    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 356       AliError("Jet branch name not set.");
 357    }
 358
 359 }
 360
 361 void AliAnalysisTaskJetCore::UserCreateOutputObjects()
 362 {
 363    // Create histograms
 364    // Called once
 365    OpenFile(1);
 366    if(!fOutputList) fOutputList = new TList;
 367    fOutputList->SetOwner(kTRUE);
 368
 369    Bool_t oldStatus = TH1::AddDirectoryStatus();
 370    TH1::AddDirectory(kFALSE);
 371
 372
 373    fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5);
 374    fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
 375    fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN");
 376    fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)");
 377    fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)");
 378    fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)");
 379    fHistEvtSelection->GetXaxis()->SetBinLabel(6,"multiplicity (rejected)");
 380
 381    fHistJetSelection = new TH1I("fHistJetSelection", "jet selection", 8, -0.5, 7.5);
 382    fHistJetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
 383    fHistJetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
 384    fHistJetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
 385    fHistJetSelection->GetXaxis()->SetBinLabel(4,"not in list");
 386    fHistJetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
 387    fHistJetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
 388    fHistJetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
 389    fHistJetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
 390
 391    fh2JetSelection = new TH2F("fh2JetSelection", "jet selection", 8, -0.5, 7.5,100,0.,200.);
 392    fh2JetSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
 393    fh2JetSelection->GetXaxis()->SetBinLabel(2,"probes IN");
 394    fh2JetSelection->GetXaxis()->SetBinLabel(3,"no matching jet");
 395    fh2JetSelection->GetXaxis()->SetBinLabel(4,"not in list");
 396    fh2JetSelection->GetXaxis()->SetBinLabel(5,"fraction cut");
 397    fh2JetSelection->GetXaxis()->SetBinLabel(6,"acceptance cut");
 398    fh2JetSelection->GetXaxis()->SetBinLabel(7,"p_{T} cut");
 399    fh2JetSelection->GetXaxis()->SetBinLabel(8,"trigger exclude mask");
 400
 401
 402    //UInt_t entries = 0; // bit coded, see GetDimParams() below
 403    //UInt_t opt = 0;  // bit coded, default (0) or high resolution (1)
 404
 405    //  Int_t bins[5]={10,10,150,150,15};
 406    //Double_t xmin[5]={0.,0.,0.,0.,0.};
 407    //Double_t xmax[5]={100.,6.5,150.,1.5,1.5};
 408    //fhnDeltaRjets = new THnSparseF("fhnDeltaRjets", "fhnDeltaRjets",5,bins,xmin,xmax);
 409
 410
 411
 412     fh2JetCoreMethod1C10 = new TH2F("JetCoreMethod1C10","",150, 0., 150.,100, 0., 1.5);
 413     fh2JetCoreMethod2C10 = new TH2F("JetCoreMethod2C10","",150, 0., 150.,100, 0., 1.5);
 414     fh2JetCoreMethod3C10 = new TH2F("JetCoreMethod3C10","",150, 0., 150.,100, 0., 1.5);
 415     fh2JetCoreMethod1C20 = new TH2F("JetCoreMethod1C20","",150, 0., 150.,100, 0., 1.5);
 416     fh2JetCoreMethod2C20 = new TH2F("JetCoreMethod2C20","",150, 0., 150.,100, 0., 1.5);
 417     fh2JetCoreMethod3C20 = new TH2F("JetCoreMethod3C20","",150, 0., 150.,100, 0., 1.5);
 418     fh2JetCoreMethod1C30 = new TH2F("JetCoreMethod1C30","",150, 0., 150.,100, 0., 1.5);
 419     fh2JetCoreMethod2C30 = new TH2F("JetCoreMethod2C30","",150, 0., 150.,100, 0., 1.5);
 420     fh2JetCoreMethod3C30 = new TH2F("JetCoreMethod3C30","",150, 0., 150.,100, 0., 1.5);
 421     fh2JetCoreMethod1C60 = new TH2F("JetCoreMethod1C60","",150, 0., 150.,100, 0., 1.5);
 422     fh2JetCoreMethod2C60 = new TH2F("JetCoreMethod2C60","",150, 0., 150.,100, 0., 1.5);
 423     fh2JetCoreMethod3C60 = new TH2F("JetCoreMethod3C60","",150, 0., 150.,100, 0., 1.5);
 424
 425     fh2SumPtInC10 = new TH2F("PtSumInC10","",150, 0., 150.,5000, 0., 50.);
 426     fh2SumPtInC20 = new TH2F("PtSumInC20","",150, 0., 150.,5000, 0., 50.);
 427     fh2SumPtInC30 = new TH2F("PtSumInC30","",150, 0., 150.,5000, 0., 50.);
 428     fh2SumPtInC60 = new TH2F("PtSumInC60","",150, 0., 150.,5000, 0., 50.);
 429     fh2SumPtOutC10 = new TH2F("PtSumOutC10","",150, 0., 150.,5000, 0., 50.);
 430     fh2SumPtOutC10b= new TH2F("PtSumOutC10b","",150, 0., 150.,500, 0., 50.);
 431     fh2SumPtOutC20 = new TH2F("PtSumOutC20","",150, 0., 150.,5000, 0., 50.);
 432     fh2SumPtOutC30 = new TH2F("PtSumOutC30","",150, 0., 150.,5000, 0., 50.);
 433     fh2SumPtOutC60 = new TH2F("PtSumOutC60","",150, 0., 150.,5000, 0., 50.);
 434
 435     fh2SumPtInC10bkg = new TH2F("PtSumInC10bkg","",150, 0., 150.,5000, 0., 50.);
 436     fh2SumPtInC20bkg = new TH2F("PtSumInC20bkg","",150, 0., 150.,5000, 0., 50.);
 437     fh2SumPtInC30bkg = new TH2F("PtSumInC30bkg","",150, 0., 150.,5000, 0., 50.);
 438     fh2SumPtInC60bkg = new TH2F("PtSumInC60bkg","",150, 0., 150.,5000, 0., 50.);
 439     fh2SumPtOutC10bkg = new TH2F("PtSumOutC10bkg","",150, 0., 150.,5000, 0., 50.);
 440     fh2SumPtOutC20bkg = new TH2F("PtSumOutC20bkg","",150, 0., 150.,5000, 0., 50.);
 441     fh2SumPtOutC30bkg = new TH2F("PtSumOutC30bkg","",150, 0., 150.,5000, 0., 50.);
 442     fh2SumPtOutC60bkg = new TH2F("PtSumOutC60bkg","",150, 0., 150.,5000, 0., 50.);
 443
 444
 445
 446     fh2DeltaRC10pt1 = new TH2F("DeltaRC10pt1","",150, 0., 150.,100,0.,1.5);
 447     fh2DeltaRC20pt1 = new TH2F("DeltaRC20pt1","",150, 0., 150.,100,0.,1.5);
 448     fh2DeltaRC30pt1 = new TH2F("DeltaRC30pt1","",150, 0., 150.,100,0.,1.5);
 449     fh2DeltaRC60pt1 = new TH2F("DeltaRC60pt1","",150, 0., 150.,100,0.,1.5);
 450     fh2DeltaRC10pt2 = new TH2F("DeltaRC10pt2","",150, 0., 150.,100,0.,1.5);
 451     fh2DeltaRC20pt2 = new TH2F("DeltaRC20pt2","",150, 0., 150.,100,0.,1.5);
 452     fh2DeltaRC30pt2 = new TH2F("DeltaRC30pt2","",150, 0., 150.,100,0.,1.5);
 453     fh2DeltaRC60pt2 = new TH2F("DeltaRC60pt2","",150, 0., 150.,100,0.,1.5);
 454     fh2DeltaRC10pt3 = new TH2F("DeltaRC10pt3","",150, 0., 150.,100,0.,1.5);
 455     fh2DeltaRC20pt3 = new TH2F("DeltaRC20pt3","",150, 0., 150.,100,0.,1.5);
 456     fh2DeltaRC30pt3 = new TH2F("DeltaRC30pt3","",150, 0., 150.,100,0.,1.5);
 457     fh2DeltaRC60pt3 = new TH2F("DeltaRC60pt3","",150, 0., 150.,100,0.,1.5);
 458     fh2DeltaRC10pt4 = new TH2F("DeltaRC10pt4","",150, 0., 150.,100,0.,1.5);
 459     fh2DeltaRC20pt4 = new TH2F("DeltaRC20pt4","",150, 0., 150.,100,0.,1.5);
 460     fh2DeltaRC30pt4 = new TH2F("DeltaRC30pt4","",150, 0., 150.,100,0.,1.5);
 461     fh2DeltaRC60pt4 = new TH2F("DeltaRC60pt4","",150, 0., 150.,100,0.,1.5);
 462
 463     fh2DeltaEtaC10pt1 = new TH2F("DeltaEtaC10pt1","",150, 0., 150.,100,-1.5,1.5);
 464     fh2DeltaEtaC20pt1 = new TH2F("DeltaEtaC20pt1","",150, 0., 150.,100,-1.5,1.5);
 465     fh2DeltaEtaC30pt1 = new TH2F("DeltaEtaC30pt1","",150, 0., 150.,100,-1.5,1.5);
 466     fh2DeltaEtaC60pt1 = new TH2F("DeltaEtaC60pt1","",150, 0., 150.,100,-1.5,1.5);
 467     fh2DeltaEtaC10pt2 = new TH2F("DeltaEtaC10pt2","",150, 0., 150.,100,-1.5,1.5);
 468     fh2DeltaEtaC20pt2 = new TH2F("DeltaEtaC20pt2","",150, 0., 150.,100,-1.5,1.5);
 469     fh2DeltaEtaC30pt2 = new TH2F("DeltaEtaC30pt2","",150, 0., 150.,100,-1.5,1.5);
 470     fh2DeltaEtaC60pt2 = new TH2F("DeltaEtaC60pt2","",150, 0., 150.,100,-1.5,1.5);
 471     fh2DeltaEtaC10pt3 = new TH2F("DeltaEtaC10pt3","",150, 0., 150.,100,-1.5,1.5);
 472     fh2DeltaEtaC20pt3 = new TH2F("DeltaEtaC20pt3","",150, 0., 150.,100,-1.5,1.5);
 473     fh2DeltaEtaC30pt3 = new TH2F("DeltaEtaC30pt3","",150, 0., 150.,100,-1.5,1.5);
 474     fh2DeltaEtaC60pt3 = new TH2F("DeltaEtaC60pt3","",150, 0., 150.,100,-1.5,1.5);
 475     fh2DeltaEtaC10pt4 = new TH2F("DeltaEtaC10pt4","",150, 0., 150.,100,-1.5,1.5);
 476     fh2DeltaEtaC20pt4 = new TH2F("DeltaEtaC20pt4","",150, 0., 150.,100,-1.5,1.5);
 477     fh2DeltaEtaC30pt4 = new TH2F("DeltaEtaC30pt4","",150, 0., 150.,100,-1.5,1.5);
 478     fh2DeltaEtaC60pt4 = new TH2F("DeltaEtaC60pt4","",150, 0., 150.,100,-1.5,1.5);
 479     fh2DeltaPhiC10pt1 = new TH2F("DeltaPhiC10pt1","",150, 0., 150.,100,-6.5,6.5);
 480     fh2DeltaPhiC20pt1 = new TH2F("DeltaPhiC20pt1","",150, 0., 150.,100,-6.5,6.5);
 481     fh2DeltaPhiC30pt1 = new TH2F("DeltaPhiC30pt1","",150, 0., 150.,100,-6.5,6.5);
 482     fh2DeltaPhiC60pt1 = new TH2F("DeltaPhiC60pt1","",150, 0., 150.,100,-6.5,6.5);
 483     fh2DeltaPhiC10pt2 = new TH2F("DeltaPhiC10pt2","",150, 0., 150.,100,-6.5,6.5);
 484     fh2DeltaPhiC20pt2 = new TH2F("DeltaPhiC20pt2","",150, 0., 150.,100,-6.5,6.5);
 485     fh2DeltaPhiC30pt2 = new TH2F("DeltaPhiC30pt2","",150, 0., 150.,100,-6.5,6.5);
 486     fh2DeltaPhiC60pt2 = new TH2F("DeltaPhiC60pt2","",150, 0., 150.,100,-6.5,6.5);
 487     fh2DeltaPhiC10pt3 = new TH2F("DeltaPhiC10pt3","",150, 0., 150.,100,-6.5,6.5);
 488     fh2DeltaPhiC20pt3 = new TH2F("DeltaPhiC20pt3","",150, 0., 150.,100,-6.5,6.5);
 489     fh2DeltaPhiC30pt3 = new TH2F("DeltaPhiC30pt3","",150, 0., 150.,100,-6.5,6.5);
 490     fh2DeltaPhiC60pt3 = new TH2F("DeltaPhiC60pt3","",150, 0., 150.,100,-6.5,6.5);
 491     fh2DeltaPhiC10pt4 = new TH2F("DeltaPhiC10pt4","",150, 0., 150.,100,-6.5,6.5);
 492     fh2DeltaPhiC20pt4 = new TH2F("DeltaPhiC20pt4","",150, 0., 150.,100,-6.5,6.5);
 493     fh2DeltaPhiC30pt4 = new TH2F("DeltaPhiC30pt4","",150, 0., 150.,100,-6.5,6.5);
 494     fh2DeltaPhiC60pt4 = new TH2F("DeltaPhiC60pt4","",150, 0., 150.,100,-6.5,6.5);
 495
 496     fh2AngStructpt1C10 = new TH2F("Ang struct pt1 C10","",15,0.,1.5,150,0.,5.);
 497     fh2AngStructpt2C10 = new TH2F("Ang struct pt2 C10","",15,0.,1.5,150,0.,5.);
 498     fh2AngStructpt3C10 = new TH2F("Ang struct pt3 C10","",15,0.,1.5,150,0.,5.);
 499     fh2AngStructpt4C10 = new TH2F("Ang struct pt4 C10","",15,0.,1.5,150,0.,5.);
 500     fh2AngStructpt1C20 = new TH2F("Ang struct pt1 C20","",15,0.,1.5,150,0.,5.);
 501     fh2AngStructpt2C20 = new TH2F("Ang struct pt2 C20","",15,0.,1.5,150,0.,5.);
 502     fh2AngStructpt3C20 = new TH2F("Ang struct pt3 C20","",15,0.,1.5,150,0.,5.);
 503     fh2AngStructpt4C20 = new TH2F("Ang struct pt4 C20","",15,0.,1.5,150,0.,5.);
 504     fh2AngStructpt1C30 = new TH2F("Ang struct pt1 C30","",15,0.,1.5,150,0.,5.);
 505     fh2AngStructpt2C30 = new TH2F("Ang struct pt2 C30","",15,0.,1.5,150,0.,5.);
 506     fh2AngStructpt3C30 = new TH2F("Ang struct pt3 C30","",15,0.,1.5,150,0.,5.);
 507     fh2AngStructpt4C30 = new TH2F("Ang struct pt4 C30","",15,0.,1.5,150,0.,5.);
 508     fh2AngStructpt1C60 = new TH2F("Ang struct pt1 C60","",15,0.,1.5,150,0.,5.);
 509     fh2AngStructpt2C60 = new TH2F("Ang struct pt2 C60","",15,0.,1.5,150,0.,5.);
 510     fh2AngStructpt3C60 = new TH2F("Ang struct pt3 C60","",15,0.,1.5,150,0.,5.);
 511     fh2AngStructpt4C60 = new TH2F("Ang struct pt4 C60","",15,0.,1.5,150,0.,5.);
 512
 513    fOutputList->Add(fHistEvtSelection);
 514    fOutputList->Add(fHistJetSelection);
 515    fOutputList->Add(fh2JetSelection);
 516       fOutputList->Add(fh2JetCoreMethod1C10);
 517       fOutputList->Add(fh2JetCoreMethod2C10);
 518       fOutputList->Add(fh2JetCoreMethod3C10);
 519       fOutputList->Add(fh2JetCoreMethod1C20);
 520       fOutputList->Add(fh2JetCoreMethod2C20);
 521       fOutputList->Add(fh2JetCoreMethod3C20);
 522       fOutputList->Add(fh2JetCoreMethod1C30);
 523       fOutputList->Add(fh2JetCoreMethod2C30);
 524       fOutputList->Add(fh2JetCoreMethod3C30);
 525       fOutputList->Add(fh2JetCoreMethod1C60);
 526       fOutputList->Add(fh2JetCoreMethod2C60);
 527       fOutputList->Add(fh2JetCoreMethod3C60);
 528       fOutputList->Add(fh2SumPtInC10);
 529       fOutputList->Add(fh2SumPtInC20);
 530       fOutputList->Add(fh2SumPtInC30);
 531       fOutputList->Add(fh2SumPtInC60);
 532       fOutputList->Add(fh2SumPtOutC10);
 533       fOutputList->Add(fh2SumPtOutC10b);
 534       fOutputList->Add(fh2SumPtOutC20);
 535       fOutputList->Add(fh2SumPtOutC30);
 536       fOutputList->Add(fh2SumPtOutC60);
 537
 538       fOutputList->Add(fh2SumPtInC10bkg);
 539       fOutputList->Add(fh2SumPtInC20bkg);
 540       fOutputList->Add(fh2SumPtInC30bkg);
 541       fOutputList->Add(fh2SumPtInC60bkg);
 542       fOutputList->Add(fh2SumPtOutC10bkg);
 543       fOutputList->Add(fh2SumPtOutC20bkg);
 544       fOutputList->Add(fh2SumPtOutC30bkg);
 545       fOutputList->Add(fh2SumPtOutC60bkg);
 546       fOutputList->Add(fh2DeltaRC10pt1);
 547       fOutputList->Add(fh2DeltaRC20pt1);
 548       fOutputList->Add(fh2DeltaRC30pt1);
 549       fOutputList->Add(fh2DeltaRC60pt1);
 550       fOutputList->Add(fh2DeltaRC10pt2);
 551       fOutputList->Add(fh2DeltaRC20pt2);
 552       fOutputList->Add(fh2DeltaRC30pt2);
 553       fOutputList->Add(fh2DeltaRC60pt2);
 554       fOutputList->Add(fh2DeltaRC10pt3);
 555       fOutputList->Add(fh2DeltaRC20pt3);
 556       fOutputList->Add(fh2DeltaRC30pt3);
 557       fOutputList->Add(fh2DeltaRC60pt3);
 558       fOutputList->Add(fh2DeltaRC10pt4);
 559       fOutputList->Add(fh2DeltaRC20pt4);
 560       fOutputList->Add(fh2DeltaRC30pt4);
 561       fOutputList->Add(fh2DeltaRC60pt4);
 562
 563       fOutputList->Add(fh2DeltaEtaC10pt1);
 564       fOutputList->Add(fh2DeltaEtaC20pt1);
 565       fOutputList->Add(fh2DeltaEtaC30pt1);
 566       fOutputList->Add(fh2DeltaEtaC60pt1);
 567       fOutputList->Add(fh2DeltaEtaC10pt2);
 568       fOutputList->Add(fh2DeltaEtaC20pt2);
 569       fOutputList->Add(fh2DeltaEtaC30pt2);
 570       fOutputList->Add(fh2DeltaEtaC60pt2);
 571       fOutputList->Add(fh2DeltaEtaC10pt3);
 572       fOutputList->Add(fh2DeltaEtaC20pt3);
 573       fOutputList->Add(fh2DeltaEtaC30pt3);
 574       fOutputList->Add(fh2DeltaEtaC60pt3);
 575       fOutputList->Add(fh2DeltaEtaC10pt4);
 576       fOutputList->Add(fh2DeltaEtaC20pt4);
 577       fOutputList->Add(fh2DeltaEtaC30pt4);
 578       fOutputList->Add(fh2DeltaEtaC60pt4);
 579       fOutputList->Add(fh2DeltaPhiC10pt1);
 580       fOutputList->Add(fh2DeltaPhiC20pt1);
 581       fOutputList->Add(fh2DeltaPhiC30pt1);
 582       fOutputList->Add(fh2DeltaPhiC60pt1);
 583       fOutputList->Add(fh2DeltaPhiC10pt2);
 584       fOutputList->Add(fh2DeltaPhiC20pt2);
 585       fOutputList->Add(fh2DeltaPhiC30pt2);
 586       fOutputList->Add(fh2DeltaPhiC60pt2);
 587       fOutputList->Add(fh2DeltaPhiC10pt3);
 588       fOutputList->Add(fh2DeltaPhiC20pt3);
 589       fOutputList->Add(fh2DeltaPhiC30pt3);
 590       fOutputList->Add(fh2DeltaPhiC60pt3);
 591       fOutputList->Add(fh2DeltaPhiC10pt4);
 592       fOutputList->Add(fh2DeltaPhiC20pt4);
 593       fOutputList->Add(fh2DeltaPhiC30pt4);
 594       fOutputList->Add(fh2DeltaPhiC60pt4);
 595
 596        fOutputList->Add(fh2AngStructpt1C10);
 597        fOutputList->Add(fh2AngStructpt2C10);
 598        fOutputList->Add(fh2AngStructpt3C10);
 599        fOutputList->Add(fh2AngStructpt4C10);
 600        fOutputList->Add(fh2AngStructpt1C20);
 601        fOutputList->Add(fh2AngStructpt2C20);
 602        fOutputList->Add(fh2AngStructpt3C20);
 603        fOutputList->Add(fh2AngStructpt4C20);
 604        fOutputList->Add(fh2AngStructpt1C30);
 605        fOutputList->Add(fh2AngStructpt2C30);
 606        fOutputList->Add(fh2AngStructpt3C30);
 607        fOutputList->Add(fh2AngStructpt4C30);
 608        fOutputList->Add(fh2AngStructpt1C60);
 609        fOutputList->Add(fh2AngStructpt2C60);
 610        fOutputList->Add(fh2AngStructpt3C60);
 611        fOutputList->Add(fh2AngStructpt4C60);
 612
 613
 614
 615
 616    // =========== Switch on Sumw2 for all histos ===========
 617    for (Int_t i=0; i<fOutputList->GetEntries(); ++i) {
 618       TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i));
 619       if (h1){
 620          h1->Sumw2();
 621          continue;
 622       }
 623     THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i));
 624       if (hn){
 625          hn->Sumw2();
 626       }
 627    }
 628    TH1::AddDirectory(oldStatus);
 629
 630    PostData(1, fOutputList);
 631 }
 632
 633 void AliAnalysisTaskJetCore::UserExec(Option_t *)
 634 {
 635
 636
 637    if(!strlen(fJetBranchName[0].Data()) || !strlen(fJetBranchName[1].Data())){
 638       AliError("Jet branch name not set.");
 639       return;
 640    }
 641
 642    fESD=dynamic_cast<AliESDEvent*>(InputEvent());
 643    if (!fESD) {
 644       AliError("ESD not available");
 645       fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
 646    } else {
 647       fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
 648    }
 649
 650     if(fNonStdFile.Length()!=0){
 651     // case that we have an AOD extension we need can fetch the jets from the extended output
 652     AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 653     fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
 654     if(!fAODExtension){
 655       if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
 656     }}
 657
 658
 659
 660
 661
 662    // -- event selection --
 663    fHistEvtSelection->Fill(1); // number of events before event selection
 664
 665    // physics selection
 666    AliInputEventHandler* inputHandler = (AliInputEventHandler*)
 667    ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
 668    if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){
 669       if(fDebug) Printf(" Trigger Selection: event REJECTED ... ");
 670       fHistEvtSelection->Fill(2);
 671       PostData(1, fOutputList);
 672       return;
 673    }
 674
 675    // vertex selection
 676    AliAODVertex* primVtx = fAOD->GetPrimaryVertex();
 677    Int_t nTracksPrim = primVtx->GetNContributors();
 678    if ((nTracksPrim < fMinContribVtx) ||
 679          (primVtx->GetZ() < fVtxZMin) ||
 680          (primVtx->GetZ() > fVtxZMax) ){
 681       if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
 682       fHistEvtSelection->Fill(3);
 683       PostData(1, fOutputList);
 684       return;
 685    }
 686
 687    // event class selection (from jet helper task)
 688    Int_t eventClass = AliAnalysisHelperJetTasks::EventClass();
 689    if(fDebug) Printf("Event class %d", eventClass);
 690    if (eventClass < fEvtClassMin || eventClass > fEvtClassMax){
 691       fHistEvtSelection->Fill(4);
 692       PostData(1, fOutputList);
 693       return;
 694    }
 695
 696    // centrality selection
 697    AliCentrality *cent = 0x0;
 698    Float_t centValue = 0.;
 699    if(fESD) cent = fESD->GetCentrality();
 700    if(cent) centValue = cent->GetCentralityPercentile("V0M");
 701    if(fDebug) printf("centrality: %f\n", centValue);
 702    if (centValue < fCentMin || centValue > fCentMax){
 703       fHistEvtSelection->Fill(4);
 704       PostData(1, fOutputList);
 705       return;
 706    }
 707
 708
 709    // multiplicity due to input tracks
 710    //Int_t nInputTracks = GetNInputTracks();
 711    //if (nInputTracks < fNInputTracksMin || (fNInputTracksMax > -1 && nInputTracks > fNInputTracksMax)){
 712    //   fHistEvtSelection->Fill(5);
 713    //   PostData(1, fOutputList);
 714    //   return;
 715    // }
 716
 717
 718    fHistEvtSelection->Fill(0);
 719    // accepted events
 720    // -- end event selection --
 721
 722    // get background
 723    AliAODJetEventBackground* externalBackground = 0;
 724    if(fAOD&&!externalBackground&&fBackgroundBranch.Length()){
 725       externalBackground =  (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
 726       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 727    }
 728    if(fAODExtension&&!externalBackground&&fBackgroundBranch.Length()){
 729      externalBackground =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
 730       if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
 731    }
 732
 733    Float_t rho = 0;
 734    if(externalBackground)rho = externalBackground->GetBackground(0);
 735
 736
 737    // fetch jets
 738    TClonesArray *aodJets[2];
 739    aodJets[0]=0;
 740    if(fAOD&&!aodJets[0]){
 741    aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); //
 742    aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data()));  }
 743    if(fAODExtension && !aodJets[0]){
 744      aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data())); //
 745      aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));  }
 746
 747    TList ParticleList;
 748    Int_t nT = GetListOfTracks(&ParticleList);
 749      for (Int_t iJetType = 0; iJetType < 2; iJetType++) {
 750       fListJets[iJetType]->Clear();
 751       if (!aodJets[iJetType]) continue;
 752
 753       if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
 754
 755       for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
 756          AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets[iJetType])[iJet]);
 757          if (jet) fListJets[iJetType]->Add(jet);
 758       }
 759       fListJets[iJetType]->Sort();
 760    }
 761
 762    Double_t etabig=0;
 763    Double_t ptbig=0;
 764    Double_t areabig=0;
 765    Double_t phibig=0.;
 766    Double_t etasmall=0;
 767    Double_t ptsmall=0;
 768    Double_t areasmall=0;
 769    Double_t distr=0.;
 770    Double_t phismall=0.;
 771
 772    Double_t up1[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 773    Double_t up2[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 774    Double_t up3[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 775    Double_t up4[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 776    Double_t down1[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 777    Double_t down2[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 778    Double_t down3[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 779    Double_t down4[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
 780
 781
 782
 783    for(Int_t i=0; i<fListJets[0]->GetEntries(); ++i){
 784            AliAODJet* jetbig = (AliAODJet*)(fListJets[0]->At(i));
 785            etabig  = jetbig->Eta();
 786            phibig  = jetbig->Phi();
 787            ptbig   = jetbig->Pt();
 788            if(ptbig==0) continue;
 789            areabig = jetbig->EffectiveAreaCharged();
 790            Double_t ptcorr=ptbig-rho*areabig;
 791            if(ptcorr<=0) continue;
 792            if((etabig<fJetEtaMin)||(etabig>fJetEtaMax)) continue;
 793                    Double_t dismin=100.;
 794                    Double_t ptmax=-10.;
 795                    Int_t index1=-1;
 796                    Int_t index2=-1;
 797                    Double_t fracin=0.;
 798                    Double_t sumPtIn=0.;
 799                    Double_t sumPtOut=0.;
 800            //compute sum of the pt of the tracks in a concentric cone
 801            TRefArray *genTrackList = jetbig->GetRefTracks();
 802            Int_t nTracksGenJet = genTrackList->GetEntriesFast();
 803            AliAODTrack* genTrack;
 804              for(Int_t ir=0; ir<nTracksGenJet; ++ir){
 805              genTrack = (AliAODTrack*)(genTrackList->At(ir));
 806              Float_t etr=genTrack->Eta();
 807              Float_t phir=genTrack->Phi();
 808              distr=(etr-etabig)*(etr-etabig)+(phir-phibig)*(phir-phibig);
 809              distr=TMath::Sqrt(distr);
 810              if(distr<=fRadioFrac){ fracin=fracin+genTrack->Pt();}}
 811
 812              if(centValue<10) fh2JetCoreMethod3C10->Fill(ptcorr,fracin/ptbig);
 813
 814              if((centValue>20)&&(centValue<40)) fh2JetCoreMethod3C20->Fill(ptcorr,fracin/ptbig);
 815
 816              if((centValue>30)&&(centValue<60)) fh2JetCoreMethod3C30->Fill(ptcorr,fracin/ptbig);
 817
 818              if(centValue>60) fh2JetCoreMethod3C60->Fill(ptcorr,fracin/ptbig);
 819
 820
 821
 822
 823                 for(Int_t j=0; j<fListJets[1]->GetEntries(); ++j){
 824                   AliAODJet* jetsmall = (AliAODJet*)(fListJets[1]->At(j));
 825                   etasmall  = jetsmall->Eta();
 826                   phismall = jetsmall->Phi();
 827                   ptsmall   = jetsmall->Pt();
 828                   areasmall = jetsmall->EffectiveAreaCharged();
 829                   Double_t tmpDeltaR=(phismall-phibig)*(phismall-phibig)+(etasmall-etabig)*(etasmall-etabig);
 830                   tmpDeltaR=TMath::Sqrt(tmpDeltaR);
 831                      //Fraction in the jet core
 832                     if((ptsmall>ptmax)&&(tmpDeltaR<=fRadioFrac)){ptmax=ptsmall;
 833                     index2=j;}
 834                     if(tmpDeltaR<=dismin){ dismin=tmpDeltaR;
 835                       index1=j;}} //en of loop over R=0.2 jets
 836                 //method1:most concentric jet=core
 837                 if(dismin<fMinDist){ AliAODJet* jetmethod1 = (AliAODJet*)(fListJets[1]->At(index1));
 838                   if(centValue<10) fh2JetCoreMethod1C10->Fill(ptcorr,jetmethod1->Pt()/ptbig);
 839                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod1C20->Fill(ptcorr,jetmethod1->Pt()/ptbig);
 840                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod1C30->Fill(ptcorr,jetmethod1->Pt()/ptbig);
 841                   if(centValue>60) fh2JetCoreMethod1C60->Fill(ptcorr,jetmethod1->Pt()/ptbig); }
 842                 //method2:hardest contained jet=core
 843                 if(index2!=-1){
 844                   AliAODJet* jetmethod2 = (AliAODJet*)(fListJets[1]->At(index2));
 845                   if(centValue<10) fh2JetCoreMethod2C10->Fill(ptcorr,jetmethod2->Pt()/ptbig);
 846                   if((centValue>20)&&(centValue<40)) fh2JetCoreMethod2C20->Fill(ptcorr,jetmethod2->Pt()/ptbig);
 847                   if((centValue>30)&&(centValue<60)) fh2JetCoreMethod2C30->Fill(ptcorr,jetmethod2->Pt()/ptbig);
 848                   if(centValue>60) fh2JetCoreMethod2C60->Fill(ptcorr,jetmethod2->Pt()/ptbig); }
 849
 850      Double_t R=fRadioFrac*2.;
 851
 852      for(int it = 0;it<nT;++it){
 853
 854           AliVParticle *part = (AliVParticle*)ParticleList.At(it);
 855           Float_t deltaR = jetbig->DeltaR(part);
 856           Float_t deltaEta = part->Eta()-etabig;
 857           Float_t deltaPhi = part->Phi()-phibig;
 858           if((deltaR>=0.4)&&(deltaR<=0.6))sumPtIn=sumPtIn+part->Pt();
 859           if((deltaR>=0.8)&&(deltaR<=1.))sumPtOut=sumPtOut+part->Pt();
 860           if(centValue<10.){
 861                if((ptcorr>=70.)&&(ptcorr<85.)) {fh2DeltaRC10pt1->Fill(part->Pt(),deltaR);
 862                if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R))fh2DeltaEtaC10pt1->Fill(part->Pt(),deltaEta);
 863                if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R))fh2DeltaPhiC10pt1->Fill(part->Pt(),deltaPhi);}
 864                if((ptcorr>=85.)&&(ptcorr<100.)) {fh2DeltaRC10pt2->Fill(part->Pt(),deltaR);
 865                if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC10pt2->Fill(part->Pt(),deltaEta);
 866                if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R)) fh2DeltaPhiC10pt2->Fill(part->Pt(),deltaPhi); }
 867                if((ptcorr>=100.)&&(ptcorr<120.)) {fh2DeltaRC10pt3->Fill(part->Pt(),deltaR);
 868                if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC10pt3->Fill(part->Pt(),deltaEta);
 869                if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R)) fh2DeltaPhiC10pt3->Fill(part->Pt(),deltaPhi);}
 870                if((ptcorr>=120.)&&(ptcorr<140.)) {fh2DeltaRC10pt4->Fill(part->Pt(),deltaR);
 871                if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC10pt4->Fill(part->Pt(),deltaEta);
 872                if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R))fh2DeltaPhiC10pt4->Fill(part->Pt(),deltaPhi); }}
 873
 874           if((centValue>20.)&&(centValue<40.)){
 875               if((ptcorr>=70.)&&(ptcorr<85.)) {fh2DeltaRC20pt1->Fill(part->Pt(),deltaR);
 876               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R))fh2DeltaEtaC20pt1->Fill(part->Pt(),deltaEta);
 877               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R))fh2DeltaPhiC20pt1->Fill(part->Pt(),deltaPhi);}
 878               if((ptcorr>=85.)&&(ptcorr<100.)) {fh2DeltaRC20pt2->Fill(part->Pt(),deltaR);
 879               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC20pt2->Fill(part->Pt(),deltaEta);
 880               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R)) fh2DeltaPhiC20pt2->Fill(part->Pt(),deltaPhi); }
 881               if((ptcorr>=100.)&&(ptcorr<120.)) {fh2DeltaRC20pt3->Fill(part->Pt(),deltaR);
 882               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC20pt3->Fill(part->Pt(),deltaEta);
 883               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R)) fh2DeltaPhiC20pt3->Fill(part->Pt(),deltaPhi);}
 884               if((ptcorr>=120.)&&(ptcorr<140.)) {fh2DeltaRC20pt4->Fill(part->Pt(),deltaR);
 885               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC20pt4->Fill(part->Pt(),deltaEta);
 886               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R))fh2DeltaPhiC20pt4->Fill(part->Pt(),deltaPhi); }}
 887
 888           if((centValue>30.)&&(centValue<60.)){
 889
 890               if((ptcorr>=70.)&&(ptcorr<85.)) {fh2DeltaRC30pt1->Fill(part->Pt(),deltaR);
 891               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R))fh2DeltaEtaC30pt1->Fill(part->Pt(),deltaEta);
 892               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R))fh2DeltaPhiC30pt1->Fill(part->Pt(),deltaPhi);}
 893               if((ptcorr>=85.)&&(ptcorr<100.)) {fh2DeltaRC30pt2->Fill(part->Pt(),deltaR);
 894               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC30pt2->Fill(part->Pt(),deltaEta);
 895               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R)) fh2DeltaPhiC30pt2->Fill(part->Pt(),deltaPhi); }
 896               if((ptcorr>=100.)&&(ptcorr<120.)) {fh2DeltaRC30pt3->Fill(part->Pt(),deltaR);
 897               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC30pt3->Fill(part->Pt(),deltaEta);
 898               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R)) fh2DeltaPhiC30pt3->Fill(part->Pt(),deltaPhi);}
 899               if((ptcorr>=120.)&&(ptcorr<140.)) {fh2DeltaRC30pt4->Fill(part->Pt(),deltaR);
 900               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC30pt4->Fill(part->Pt(),deltaEta);
 901               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R))fh2DeltaPhiC30pt4->Fill(part->Pt(),deltaPhi); }}
 902
 903
 904           if(centValue>60.){
 905               if((ptcorr>=70.)&&(ptcorr<85.)) {fh2DeltaRC60pt1->Fill(part->Pt(),deltaR);
 906               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R))fh2DeltaEtaC60pt1->Fill(part->Pt(),deltaEta);
 907               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R))fh2DeltaPhiC60pt1->Fill(part->Pt(),deltaPhi);}
 908               if((ptcorr>=85.)&&(ptcorr<100.)) {fh2DeltaRC60pt2->Fill(part->Pt(),deltaR);
 909               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC60pt2->Fill(part->Pt(),deltaEta);
 910               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R)) fh2DeltaPhiC60pt2->Fill(part->Pt(),deltaPhi); }
 911               if((ptcorr>=100.)&&(ptcorr<120.)) {fh2DeltaRC60pt3->Fill(part->Pt(),deltaR);
 912               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC60pt3->Fill(part->Pt(),deltaEta);
 913               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R)) fh2DeltaPhiC60pt3->Fill(part->Pt(),deltaPhi);}
 914               if((ptcorr>=120.)&&(ptcorr<140.)) {fh2DeltaRC60pt4->Fill(part->Pt(),deltaR);
 915               if((part->Phi()<=phibig+R)&&(part->Phi()>=phibig-R)) fh2DeltaEtaC60pt4->Fill(part->Pt(),deltaEta);
 916               if((part->Eta()<=etabig+R)&&(part->Eta()>=etabig-R))fh2DeltaPhiC60pt4->Fill(part->Pt(),deltaPhi);}}
 917
 918      }
 919      //end of track loop
 920      Double_t coronain=rho*TMath::Pi()*(1.-0.8*0.8);
 921      Double_t coronaout=rho*TMath::Pi()*(0.6*0.6-0.4*0.4);
 922      if(centValue<10.){
 923      fh2SumPtInC10bkg->Fill(ptcorr,coronain/ptbig);
 924      fh2SumPtOutC10bkg->Fill(ptcorr,coronaout/ptbig);
 925      fh2SumPtInC10->Fill(ptcorr,sumPtIn/ptbig);
 926      fh2SumPtOutC10->Fill(ptcorr,sumPtOut/ptbig);
 927      fh2SumPtOutC10b->Fill(ptcorr,sumPtOut/ptbig);
 928 }
 929      if((centValue>20.)&&(centValue<40.)){
 930      fh2SumPtInC20bkg->Fill(ptcorr,coronain/ptbig);
 931      fh2SumPtOutC20bkg->Fill(ptcorr,coronaout/ptbig);
 932      fh2SumPtInC20->Fill(ptcorr,sumPtIn/ptbig);
 933      fh2SumPtOutC20->Fill(ptcorr,sumPtOut/ptbig); }
 934      if((centValue>30.)&&(centValue<60.)){
 935      fh2SumPtInC30bkg->Fill(ptcorr,coronain/ptbig);
 936      fh2SumPtOutC30bkg->Fill(ptcorr,coronaout/ptbig);
 937      fh2SumPtInC30->Fill(ptcorr,sumPtIn/ptbig);
 938      fh2SumPtOutC30->Fill(ptcorr,sumPtOut/ptbig); }
 939      if(centValue>60.){
 940      fh2SumPtInC60bkg->Fill(ptcorr,coronain/ptbig);
 941      fh2SumPtOutC60bkg->Fill(ptcorr,coronaout/ptbig);
 942      fh2SumPtInC60->Fill(ptcorr,sumPtIn/ptbig);
 943      fh2SumPtOutC60->Fill(ptcorr,sumPtOut/ptbig);}
 944
 945      //////////////////ANGULAR STRUCTURE//////////////////////////////////////
 946
 947      //tracks up to R=0.8 distant from the jet axis
 948      if(fAngStructCloseTracks==1){
 949       TList CloseTrackList;
 950       Int_t nn=GetListOfTracksCloseToJet(&CloseTrackList,jetbig);
 951       Double_t difR=0.04;
 952       for(Int_t l=0;l<15;l++){
 953       Double_t rr=l*0.1+0.1;
 954        for(int it = 0;it<nn;++it){
 955            AliVParticle *part1 = (AliVParticle*)CloseTrackList.At(it);
 956            for(int itu=it+1;itu<CloseTrackList.GetEntries();itu++){
 957            AliVParticle *part2 = (AliVParticle*)CloseTrackList.At(itu);
 958            Double_t ptm=part1->Pt();
 959            Double_t ptn=part2->Pt();
 960            Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
 961                       Rnm=TMath::Sqrt(Rnm);
 962                       Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
 963                       Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
 964                       if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
 965                                                         down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
 966                       if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
 967                                                         down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
 968                       if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
 969                                                          down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
 970                       if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
 971                         down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}
 972      }
 973
 974      //only jet constituents
 975       if(fAngStructCloseTracks==0){
 976
 977       Double_t difR=0.04;
 978       for(Int_t l=0;l<15;l++){
 979       Double_t rr=l*0.1+0.1;
 980
 981
 982       AliAODTrack* part1;
 983       AliAODTrack* part2;
 984           for(Int_t it=0; it<nTracksGenJet; ++it){
 985              part1 = (AliAODTrack*)(genTrackList->At(it));
 986            for(Int_t itu=0; itu<nTracksGenJet; ++itu){
 987              part2 = (AliAODTrack*)(genTrackList->At(itu));
 988            Double_t ptm=part1->Pt();
 989            Double_t ptn=part2->Pt();
 990            Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
 991                       Rnm=TMath::Sqrt(Rnm);
 992                       Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
 993                       Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
 994                       if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
 995                                                         down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
 996                       if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
 997                                                         down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
 998                       if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
 999                                                          down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
1000                       if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
1001                         down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}}
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017    }
1018
1019
1020      //end loop over R=0.4 jets
1021
1022      for(Int_t l=0;l<15;l++){
1023      Double_t rr=l*0.1+0.1;
1024         if(down1[l]!=0){
1025         if(centValue<10.)fh2AngStructpt1C10->Fill(rr,rr*up1[l]/down1[l]);
1026         if(centValue>20. && centValue<40.) fh2AngStructpt1C20->Fill(rr,rr*up1[l]/down1[l]);
1027         if(centValue>30. && centValue<60.) fh2AngStructpt1C30->Fill(rr,rr*up1[l]/down1[l]);
1028         if(centValue>60.) fh2AngStructpt1C60->Fill(rr,rr*up1[l]/down1[l]);}
1029         if(down2[l]!=0){
1030         if(centValue<10.) fh2AngStructpt2C10->Fill(rr,rr*up2[l]/down2[l]);
1031         if(centValue>20. && centValue<40.) fh2AngStructpt2C20->Fill(rr,rr*up2[l]/down2[l]);
1032         if(centValue>30. && centValue<60.) fh2AngStructpt2C30->Fill(rr,rr*up2[l]/down2[l]);
1033         if(centValue>60.) fh2AngStructpt2C60->Fill(rr,rr*up2[l]/down2[l]);}
1034         if(down3[l]!=0){
1035         if(centValue<10.) fh2AngStructpt3C10->Fill(rr,rr*up3[l]/down3[l]);
1036         if(centValue>20. && centValue<40.) fh2AngStructpt3C20->Fill(rr,rr*up3[l]/down3[l]);
1037         if(centValue>30. && centValue<60.) fh2AngStructpt3C30->Fill(rr,rr*up3[l]/down3[l]);
1038         if(centValue>60.) fh2AngStructpt3C60->Fill(rr,rr*up3[l]/down3[l]);}
1039         if(down4[l]!=0){
1040         if(centValue<10.) fh2AngStructpt4C10->Fill(rr,rr*up4[l]/down4[l]);
1041         if(centValue>20. && centValue<40.) fh2AngStructpt4C20->Fill(rr,rr*up4[l]/down4[l]);
1042         if(centValue>30. && centValue<60.) fh2AngStructpt4C30->Fill(rr,rr*up4[l]/down4[l]);
1043         if(centValue>60.) fh2AngStructpt4C60->Fill(rr,rr*up4[l]/down4[l]);}}
1044
1045
1046
1047
1048
1049
1050
1051    PostData(1, fOutputList);
1052 }
1053
1054 void AliAnalysisTaskJetCore::Terminate(const Option_t *)
1055 {
1056    // Draw result to the screen
1057    // Called once at the end of the query
1058
1059    if (!GetOutputData(1))
1060    return;
1061 }
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073 Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
1074
1075     Int_t iCount = 0;
1076
1077   Int_t fFilterMask=272;
1078     for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
1079       AliAODTrack *tr = fAOD->GetTrack(it);
1080       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1081       if(TMath::Abs(tr->Eta())>0.9)continue;
1082       if(tr->Pt()<0.15)continue;
1083       list->Add(tr);
1084       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
1085       iCount++;
1086     }
1087
1088    list->Sort();
1089   return iCount;
1090
1091 }
1092
1093  Int_t  AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
1094
1095     Int_t iCount = 0;
1096
1097   Int_t fFilterMask=272;
1098     for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
1099       AliAODTrack *tr = fAOD->GetTrack(it);
1100       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1101       if(TMath::Abs(tr->Eta())>0.9)continue;
1102       if(tr->Pt()<0.15)continue;
1103       Double_t disR=jetbig->DeltaR(tr);
1104       if(disR>0.8)  continue;
1105       list->Add(tr);
1106       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
1107       iCount++;
1108     }
1109
1110    list->Sort();
1111    return iCount;
1112
1113 }
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125 Int_t AliAnalysisTaskJetCore::GetNInputTracks()
1126 {
1127
1128    Int_t nInputTracks = 0;
1129
1130    TString jbname(fJetBranchName[1]);
1131    //needs complete event, use jets without background subtraction
1132    for(Int_t i=1; i<=3; ++i){
1133       if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
1134    }
1135    // use only HI event
1136    if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
1137    if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
1138
1139    if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
1140    TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
1141    if(!tmpAODjets){
1142       Printf("Jet branch %s not found", jbname.Data());
1143       Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
1144       return -1;
1145    }
1146
1147    for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
1148       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
1149       if(!jet) continue;
1150       TRefArray *trackList = jet->GetRefTracks();
1151       Int_t nTracks = trackList->GetEntriesFast();
1152       nInputTracks += nTracks;
1153       if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
1154    }
1155    if(fDebug) Printf("---> input tracks: %d", nInputTracks);
1156
1157    return nInputTracks;
1158 }
1159
1160
1161
1162 Double_t AliAnalysisTaskJetCore::RelativePhi(Double_t mphi,Double_t vphi){
1163
1164   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1165   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1166   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1167   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1168   double dphi = mphi-vphi;
1169   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1170   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1171   return dphi;//dphi in [-Pi, Pi]
1172 }
1173
1174
1175