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      if(!fESD)aod = fAODIn;
 976      else aod = fAODOut;
 977      Int_t index=-1;
 978      Double_t ptmax=-10;
 979     for(int it = 0;it < aod->GetNumberOfTracks();++it){
 980       AliAODTrack *tr = aod->GetTrack(it);
 981       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
 982       if(TMath::Abs(tr->Eta())>0.9)continue;
 983       if(tr->Pt()<0.15)continue;
 984       list->Add(tr);
 985       iCount++;
 986       if(tr->Pt()>ptmax){ ptmax=tr->Pt();
 987       index=iCount-1;}
 988
 989     }
 990
 991
 992   return index;
 993
 994 }
 995
 996    Int_t  AliAnalysisTaskJetCore::GetHardestTrackBackToJet(AliAODJet *jetbig){
 997
 998     AliAODEvent *aod = 0;
 999     if(!fESD)aod = fAODIn;
1000     else aod = fAODOut;
1001     Int_t index=-1;
1002     Double_t ptmax=-10;
1003     Double_t dphi=0;
1004     Double_t dif=0;
1005     Int_t iCount=0;
1006     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1007       AliAODTrack *tr = aod->GetTrack(it);
1008       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1009       if(TMath::Abs(tr->Eta())>0.9)continue;
1010       if(tr->Pt()<0.15)continue;
1011       iCount=iCount+1;
1012       dphi=RelativePhi(tr->Phi(),jetbig->Phi());
1013       if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;
1014       if(tr->Pt()>ptmax){ ptmax=tr->Pt();
1015       index=iCount-1;
1016       dif=dphi;  }}
1017
1018       return index;
1019
1020    }
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030  Int_t  AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
1031
1032     Int_t iCount = 0;
1033       AliAODEvent *aod = 0;
1034      if(!fESD)aod = fAODIn;
1035      else aod = fAODOut;
1036
1037       for(int it = 0;it < aod->GetNumberOfTracks();++it){
1038       AliAODTrack *tr = aod->GetTrack(it);
1039       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1040       if(TMath::Abs(tr->Eta())>0.9)continue;
1041       if(tr->Pt()<0.15)continue;
1042       Double_t disR=jetbig->DeltaR(tr);
1043       if(disR>0.8)  continue;
1044       list->Add(tr);
1045       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
1046       iCount++;
1047     }
1048
1049    list->Sort();
1050    return iCount;
1051
1052 }
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064 Int_t AliAnalysisTaskJetCore::GetNInputTracks()
1065 {
1066
1067    Int_t nInputTracks = 0;
1068      AliAODEvent *aod = 0;
1069      if(!fESD)aod = fAODIn;
1070      else aod = fAODOut;
1071    TString jbname(fJetBranchName[1]);
1072    //needs complete event, use jets without background subtraction
1073    for(Int_t i=1; i<=3; ++i){
1074       if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
1075    }
1076    // use only HI event
1077    if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
1078    if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
1079
1080    if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
1081    TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(aod->FindListObject(jbname.Data()));
1082    if(!tmpAODjets){
1083       Printf("Jet branch %s not found", jbname.Data());
1084       Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
1085       return -1;
1086    }
1087
1088    for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
1089       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
1090       if(!jet) continue;
1091       TRefArray *trackList = jet->GetRefTracks();
1092       Int_t nTracks = trackList->GetEntriesFast();
1093       nInputTracks += nTracks;
1094       if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
1095    }
1096    if(fDebug) Printf("---> input tracks: %d", nInputTracks);
1097
1098    return nInputTracks;
1099 }
1100
1101
1102
1103 Double_t AliAnalysisTaskJetCore::RelativePhi(Double_t mphi,Double_t vphi){
1104
1105   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1106   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1107   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1108   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1109   double dphi = mphi-vphi;
1110   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1111   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1112   return dphi;//dphi in [-Pi, Pi]
1113 }
1114
1115 Int_t AliAnalysisTaskJetCore::GetPhiBin(Double_t phi)
1116 {
1117     Int_t phibin=-1;
1118     if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1119     Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1120     phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1121     if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1122     return phibin;
1123 }
1124
1125
1126
1127
1128 THnSparse* AliAnalysisTaskJetCore::NewTHnSparseF(const char* name, UInt_t entries)
1129 {
1130    // generate new THnSparseF, axes are defined in GetDimParams()
1131
1132    Int_t count = 0;
1133    UInt_t tmp = entries;
1134    while(tmp!=0){
1135       count++;
1136       tmp = tmp &~ -tmp;  // clear lowest bit
1137    }
1138
1139    TString hnTitle(name);
1140    const Int_t dim = count;
1141    Int_t nbins[dim];
1142    Double_t xmin[dim];
1143    Double_t xmax[dim];
1144
1145    Int_t i=0;
1146    Int_t c=0;
1147    while(c<dim && i<32){
1148       if(entries&(1<<i)){
1149
1150          TString label("");
1151          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1152          hnTitle += Form(";%s",label.Data());
1153          c++;
1154       }
1155
1156       i++;
1157    }
1158    hnTitle += ";";
1159
1160    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1161 }
1162
1163 void AliAnalysisTaskJetCore::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1164 {
1165    // stores label and binning of axis for THnSparse
1166
1167    const Double_t pi = TMath::Pi();
1168
1169    switch(iEntry){
1170
1171    case 0:
1172       label = "V0 centrality (%)";
1173
1174          nbins = 10;
1175          xmin = 0.;
1176          xmax = 100.;
1177          break;
1178
1179
1180    case 1:
1181       label = "corrected jet pt";
1182          nbins = 20;
1183          xmin = 0.;
1184          xmax = 200.;
1185           break;
1186
1187
1188    case 2:
1189       label = "track pT";
1190
1191          nbins = 9;
1192          xmin = 0.;
1193          xmax = 150;
1194          break;
1195
1196
1197     case 3:
1198       label = "deltaR";
1199       nbins = 15;
1200       xmin = 0.;
1201       xmax = 1.5;
1202       break;
1203
1204
1205
1206    case 4:
1207       label = "deltaEta";
1208       nbins = 8;
1209       xmin = -1.6;
1210       xmax = 1.6;
1211       break;
1212
1213
1214   case 5:
1215       label = "deltaPhi";
1216       nbins = 90;
1217       xmin = -0.5*pi;
1218       xmax = 1.5*pi;
1219       break;
1220
1221
1222
1223     case 6:
1224       label = "leading track";
1225       nbins = 13;
1226       xmin = 0;
1227       xmax = 50;
1228       break;
1229
1230      case 7:
1231
1232       label = "trigger track";
1233       nbins =10;
1234       xmin = 0;
1235       xmax = 50;
1236       break;
1237
1238
1239
1240
1241
1242
1243
1244
1245    }
1246
1247 }
1248