PWGJE/AliAnalysisTaskJetCore.cxx

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