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