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