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