PWGJE/AliAnalysisTaskJetCore.cxx

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