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=jetPhi-partback->Phi();
 923                   if(dPhi>2*TMath::Pi()) dPhi -= 2*TMath::Pi();
 924                   if(dPhi<-2*TMath::Pi()) dPhi += 2*TMath::Pi();
 925                   if(dPhi<-0.5*TMath::Pi()) dPhi += 2*TMath::Pi();
 926                   if(dPhi>1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
 927
 928                   Double_t fill[] = {partback->Pt(),jetPtCorr,dPhi,jetArea,centValue};
 929                   fHJetPhiCorr->Fill(fill);
 930                 }
 931             }
 932           /////////////////////////////////////////////////////////////////////////////
 933           /////////////////////////////////////////////////////////////////////////////
 934
 935
 936      //////////////////ANGULAR STRUCTURE//////////////////////////////////////
 937
 938      //tracks up to R=0.8 distant from the jet axis
 939    //   if(fAngStructCloseTracks==1){
 940    //    TList CloseTrackList;
 941    //    Int_t nn=GetListOfTracksCloseToJet(&CloseTrackList,jetbig);
 942    //    Double_t difR=0.04;
 943    //    for(Int_t l=0;l<15;l++){
 944    //    Double_t rr=l*0.1+0.1;
 945    //     for(int it = 0;it<nn;++it){
 946    //         AliVParticle *part1 = (AliVParticle*)CloseTrackList.At(it);
 947    //         for(int itu=it+1;itu<CloseTrackList.GetEntries();itu++){
 948    //         AliVParticle *part2 = (AliVParticle*)CloseTrackList.At(itu);
 949    //         Double_t ptm=part1->Pt();
 950    //         Double_t ptn=part2->Pt();
 951    //         Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
 952    //                    Rnm=TMath::Sqrt(Rnm);
 953    //                    Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
 954    //                    Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
 955    //                    if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
 956    //                                                   down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
 957    //                    if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
 958    //                                                   down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
 959    //                    if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
 960    //                                                    down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
 961    //                    if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
 962    //                   down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}
 963    //   }
 964
 965    //   //only jet constituents
 966    //    if(fAngStructCloseTracks==2){
 967
 968    //    Double_t difR=0.04;
 969    //    for(Int_t l=0;l<15;l++){
 970    //    Double_t rr=l*0.1+0.1;
 971
 972
 973    //    AliAODTrack* part1;
 974    //    AliAODTrack* part2;
 975
 976    //        TRefArray *genTrackListb = jetbig->GetRefTracks();
 977    //        Int_t nTracksGenJetb = genTrackListb->GetEntriesFast();
 978
 979
 980
 981    //        for(Int_t it=0; it<nTracksGenJetb; ++it){
 982    //           part1 = (AliAODTrack*)(genTrackListb->At(it));
 983    //         for(Int_t itu=0; itu<nTracksGenJetb; ++itu){
 984    //           part2 = (AliAODTrack*)(genTrackListb->At(itu));
 985    //         Double_t ptm=part1->Pt();
 986    //         Double_t ptn=part2->Pt();
 987    //         Double_t Rnm = (part1->Eta()-part2->Eta())*(part1->Eta()-part2->Eta())+(part1->Phi()-part2->Phi())*(part1->Phi()-part2->Phi());
 988    //                    Rnm=TMath::Sqrt(Rnm);
 989    //                    Double_t deltag=(1./(TMath::Sqrt(2*TMath::Pi())*difR))*TMath::Exp(-1.*(rr-Rnm)*(rr-Rnm)/(2.*difR*difR));
 990    //                    Double_t stepf=0.5*(1.+TMath::Erf((rr-Rnm)/(TMath::Sqrt(2.)*difR)));
 991    //                    if((ptcorr<85.) && (ptcorr>=70.)){up1[l]=up1[l]+ptm*ptn*Rnm*Rnm*deltag;
 992    //                                                   down1[l]=down1[l]+ptm*ptn*Rnm*Rnm*stepf;}
 993    //                    if((ptcorr<100.) && (ptcorr>=85.)){up2[l]=up2[l]+ptm*ptn*Rnm*Rnm*deltag;
 994    //                                                   down2[l]=down2[l]+ptm*ptn*Rnm*Rnm*stepf;}
 995    //                    if((ptcorr<120.) && (ptcorr>=100.)){up3[l]=up3[l]+ptm*ptn*Rnm*Rnm*deltag;
 996    //                                                    down3[l]=down3[l]+ptm*ptn*Rnm*Rnm*stepf;}
 997    //                    if((ptcorr<140.) && (ptcorr>=120.)){up4[l]=up4[l]+ptm*ptn*Rnm*Rnm*deltag;
 998    //                   down4[l]=down4[l]+ptm*ptn*Rnm*Rnm*stepf;}}}}}
 999    // }
1000      // //end loop over R=0.4 jets
1001      // if(fAngStructCloseTracks>0){
1002      // for(Int_t l=0;l<15;l++){
1003      // Double_t rr=l*0.1+0.1;
1004      //    if(down1[l]!=0){
1005      //         if(centValue<10.)fh2AngStructpt1C10->Fill(rr,rr*up1[l]/down1[l]);
1006      //    if(centValue>20. && centValue<40.) fh2AngStructpt1C20->Fill(rr,rr*up1[l]/down1[l]);
1007      //    if(centValue>30. && centValue<60.) fh2AngStructpt1C30->Fill(rr,rr*up1[l]/down1[l]);
1008      //    if(centValue>60.) fh2AngStructpt1C60->Fill(rr,rr*up1[l]/down1[l]);}
1009      //    if(down2[l]!=0){
1010      //         if(centValue<10.) fh2AngStructpt2C10->Fill(rr,rr*up2[l]/down2[l]);
1011      //    if(centValue>20. && centValue<40.) fh2AngStructpt2C20->Fill(rr,rr*up2[l]/down2[l]);
1012      //    if(centValue>30. && centValue<60.) fh2AngStructpt2C30->Fill(rr,rr*up2[l]/down2[l]);
1013      //    if(centValue>60.) fh2AngStructpt2C60->Fill(rr,rr*up2[l]/down2[l]);}
1014      //    if(down3[l]!=0){
1015      //         if(centValue<10.) fh2AngStructpt3C10->Fill(rr,rr*up3[l]/down3[l]);
1016      //    if(centValue>20. && centValue<40.) fh2AngStructpt3C20->Fill(rr,rr*up3[l]/down3[l]);
1017      //    if(centValue>30. && centValue<60.) fh2AngStructpt3C30->Fill(rr,rr*up3[l]/down3[l]);
1018      //    if(centValue>60.) fh2AngStructpt3C60->Fill(rr,rr*up3[l]/down3[l]);}
1019      //    if(down4[l]!=0){
1020      //         if(centValue<10.) fh2AngStructpt4C10->Fill(rr,rr*up4[l]/down4[l]);
1021      //    if(centValue>20. && centValue<40.) fh2AngStructpt4C20->Fill(rr,rr*up4[l]/down4[l]);
1022      //    if(centValue>30. && centValue<60.) fh2AngStructpt4C30->Fill(rr,rr*up4[l]/down4[l]);
1023      //    if(centValue>60.) fh2AngStructpt4C60->Fill(rr,rr*up4[l]/down4[l]);}}}
1024
1025
1026
1027
1028
1029
1030
1031    PostData(1, fOutputList);
1032 }
1033
1034 void AliAnalysisTaskJetCore::Terminate(const Option_t *)
1035 {
1036    // Draw result to the screen
1037    // Called once at the end of the query
1038
1039    if (!GetOutputData(1))
1040    return;
1041 }
1042
1043
1044
1045
1046
1047
1048 Int_t  AliAnalysisTaskJetCore::GetListOfTracks(TList *list){
1049
1050       Int_t iCount = 0;
1051       AliAODEvent *aod = 0;
1052
1053      if(!fESD)aod = fAODIn;
1054      else aod = fAODOut;
1055       Int_t index=-1;
1056      Double_t ptmax=-10;
1057
1058
1059
1060      for(int it = 0;it < aod->GetNumberOfTracks();++it){
1061       AliAODTrack *tr = aod->GetTrack(it);
1062       Bool_t bGood = false;
1063       if(fFilterType == 0)bGood = true;
1064       else if(fFilterType == 1)bGood = tr->IsHybridTPCConstrainedGlobal();
1065       else if(fFilterType == 2)bGood = tr->IsHybridGlobalConstrainedGlobal();
1066      if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1067       if(bGood==false) continue;
1068       if(TMath::Abs(tr->Eta())>0.9)continue;
1069       if(tr->Pt()<0.15)continue;
1070       list->Add(tr);
1071       iCount++;
1072       if(fFilterType==2 && fFilterMaskBestPt>0){// only set the trigger track index for good quality tracks
1073         if(tr->TestFilterBit(fFilterMaskBestPt)){
1074           if(tr->Pt()>ptmax){
1075             ptmax=tr->Pt();
1076             index=iCount-1;
1077           }
1078         }
1079       }
1080       else{
1081         if(tr->Pt()>ptmax){
1082           ptmax=tr->Pt();
1083           index=iCount-1;
1084         }
1085       }
1086      }
1087
1088
1089     // else if (type == kTrackAODMCCharged) {
1090     // TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1091     // if(!tca)return iCount;
1092     // for(int it = 0;it < tca->GetEntriesFast();++it){
1093     //   AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1094     //   if(!part)continue;
1095     //   if(part->Pt()<0.15)continue;
1096     //   if(!part->IsPhysicalPrimary())continue;
1097     //   if(part->Charge()==0)continue;
1098     //   if(TMath::Abs(part->Eta())>0.9)continue;
1099     //   list->Add(part);
1100     //   iCount++;
1101     //   if(part->Pt()>ptmax){ ptmax=part->Pt();
1102     //  index=iCount-1;}}}
1103       return index;
1104
1105 }
1106
1107    Int_t  AliAnalysisTaskJetCore::GetHardestTrackBackToJet(AliAODJet *jetbig){
1108
1109     AliAODEvent *aod = 0;
1110     if(!fESD)aod = fAODIn;
1111     else aod = fAODOut;
1112     Int_t index=-1;
1113     Double_t ptmax=-10;
1114     Double_t dphi=0;
1115     Double_t dif=0;
1116     Int_t iCount=0;
1117     for(int it = 0;it < aod->GetNumberOfTracks();++it){
1118       AliAODTrack *tr = aod->GetTrack(it);
1119       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1120       if(TMath::Abs(tr->Eta())>0.9)continue;
1121       if(tr->Pt()<0.15)continue;
1122       iCount=iCount+1;
1123       dphi=RelativePhi(tr->Phi(),jetbig->Phi());
1124       if(TMath::Abs(dphi)<TMath::Pi()-0.6) continue;
1125       if(tr->Pt()>ptmax){ ptmax=tr->Pt();
1126       index=iCount-1;
1127       dif=dphi;  }}
1128
1129       return index;
1130
1131    }
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141  Int_t  AliAnalysisTaskJetCore::GetListOfTracksCloseToJet(TList *list,AliAODJet *jetbig){
1142
1143     Int_t iCount = 0;
1144       AliAODEvent *aod = 0;
1145      if(!fESD)aod = fAODIn;
1146      else aod = fAODOut;
1147
1148       for(int it = 0;it < aod->GetNumberOfTracks();++it){
1149       AliAODTrack *tr = aod->GetTrack(it);
1150       if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1151       if(TMath::Abs(tr->Eta())>0.9)continue;
1152       if(tr->Pt()<0.15)continue;
1153       Double_t disR=jetbig->DeltaR(tr);
1154       if(disR>0.8)  continue;
1155       list->Add(tr);
1156       //cout<<fAOD->GetNumberOfTracks()<<" "<<tr->Pt()<<endl;
1157       iCount++;
1158     }
1159
1160    list->Sort();
1161    return iCount;
1162
1163 }
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175 Int_t AliAnalysisTaskJetCore::GetNInputTracks()
1176 {
1177
1178    Int_t nInputTracks = 0;
1179      AliAODEvent *aod = 0;
1180      if(!fESD)aod = fAODIn;
1181      else aod = fAODOut;
1182    TString jbname(fJetBranchName[1]);
1183    //needs complete event, use jets without background subtraction
1184    for(Int_t i=1; i<=3; ++i){
1185       if(jbname.Contains(Form("B%d",i))) jbname.ReplaceAll(Form("B%d",i),"B0");
1186    }
1187    // use only HI event
1188    if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
1189    if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
1190
1191    if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
1192    TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(aod->FindListObject(jbname.Data()));
1193    if(!tmpAODjets){
1194       Printf("Jet branch %s not found", jbname.Data());
1195       Printf("AliAnalysisTaskJetCore::GetNInputTracks FAILED");
1196       return -1;
1197    }
1198
1199    for (Int_t iJet=0; iJet<tmpAODjets->GetEntriesFast(); iJet++){
1200       AliAODJet *jet = dynamic_cast<AliAODJet*>((*tmpAODjets)[iJet]);
1201       if(!jet) continue;
1202       TRefArray *trackList = jet->GetRefTracks();
1203       Int_t nTracks = trackList->GetEntriesFast();
1204       nInputTracks += nTracks;
1205       if(fDebug) Printf("#jet%d: %d tracks", iJet, nTracks);
1206    }
1207    if(fDebug) Printf("---> input tracks: %d", nInputTracks);
1208
1209    return nInputTracks;
1210 }
1211
1212
1213
1214 Double_t AliAnalysisTaskJetCore::RelativePhi(Double_t mphi,Double_t vphi){
1215
1216   if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1217   else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1218   if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1219   else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1220   double dphi = mphi-vphi;
1221   if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1222   else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1223   return dphi;//dphi in [-Pi, Pi]
1224 }
1225
1226 Int_t AliAnalysisTaskJetCore::GetPhiBin(Double_t phi)
1227 {
1228     Int_t phibin=-1;
1229     if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1230     Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1231     phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1232     if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1233     return phibin;
1234 }
1235
1236
1237
1238
1239 THnSparse* AliAnalysisTaskJetCore::NewTHnSparseF(const char* name, UInt_t entries)
1240 {
1241    // generate new THnSparseF, axes are defined in GetDimParams()
1242
1243    Int_t count = 0;
1244    UInt_t tmp = entries;
1245    while(tmp!=0){
1246       count++;
1247       tmp = tmp &~ -tmp;  // clear lowest bit
1248    }
1249
1250    TString hnTitle(name);
1251    const Int_t dim = count;
1252    Int_t nbins[dim];
1253    Double_t xmin[dim];
1254    Double_t xmax[dim];
1255
1256    Int_t i=0;
1257    Int_t c=0;
1258    while(c<dim && i<32){
1259       if(entries&(1<<i)){
1260
1261          TString label("");
1262          GetDimParams(i, label, nbins[c], xmin[c], xmax[c]);
1263          hnTitle += Form(";%s",label.Data());
1264          c++;
1265       }
1266
1267       i++;
1268    }
1269    hnTitle += ";";
1270
1271    return new THnSparseF(name, hnTitle.Data(), dim, nbins, xmin, xmax);
1272 }
1273
1274 void AliAnalysisTaskJetCore::GetDimParams(Int_t iEntry, TString &label, Int_t &nbins, Double_t &xmin, Double_t &xmax)
1275 {
1276    // stores label and binning of axis for THnSparse
1277
1278    const Double_t pi = TMath::Pi();
1279
1280    switch(iEntry){
1281
1282    case 0:
1283       label = "V0 centrality (%)";
1284
1285          nbins = 10;
1286          xmin = 0.;
1287          xmax = 100.;
1288          break;
1289
1290
1291    case 1:
1292       label = "corrected jet pt";
1293          nbins = 20;
1294          xmin = 0.;
1295          xmax = 200.;
1296           break;
1297
1298
1299    case 2:
1300       label = "track pT";
1301
1302          nbins = 9;
1303          xmin = 0.;
1304          xmax = 150;
1305          break;
1306
1307
1308     case 3:
1309       label = "deltaR";
1310       nbins = 15;
1311       xmin = 0.;
1312       xmax = 1.5;
1313       break;
1314
1315
1316
1317    case 4:
1318       label = "deltaEta";
1319       nbins = 8;
1320       xmin = -1.6;
1321       xmax = 1.6;
1322       break;
1323
1324
1325   case 5:
1326       label = "deltaPhi";
1327       nbins = 90;
1328       xmin = -0.5*pi;
1329       xmax = 1.5*pi;
1330       break;
1331
1332
1333
1334     case 6:
1335       label = "leading track";
1336       nbins = 13;
1337       xmin = 0;
1338       xmax = 50;
1339       break;
1340
1341      case 7:
1342
1343       label = "trigger track";
1344       nbins =10;
1345       xmin = 0;
1346       xmax = 50;
1347       break;
1348
1349
1350
1351
1352
1353
1354
1355
1356    }
1357
1358 }
1359