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