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