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