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