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