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