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