PWGJE/FlavourJetTasks/AliAnalysisTaskFlavourJetCorrelations.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15 //
  16 //
  17 //  Analysis Taks for reconstructed particle correlation
  18 //  (first implementation done for D mesons) with jets
  19 //  (use the so called Emcal framework)
  20 //
  21 //-----------------------------------------------------------------------
  22 // Authors:
  23 // C. Bianchin (Utrecht University) chiara.bianchin@cern.ch
  24 // A. Grelli (Utrecht University) a.grelli@uu.nl
  25 // X. Zhang (LBNL)  XMZhang@lbl.gov
  26 //-----------------------------------------------------------------------
  27
  28 #include <TDatabasePDG.h>
  29 #include <TParticle.h>
  30 #include "TROOT.h"
  31 #include <TH3F.h>
  32 #include <THnSparse.h>
  33 #include <TSystem.h>
  34 #include <TObjectTable.h>
  35
  36 #include "AliAnalysisTaskFlavourJetCorrelations.h"
  37 #include "AliAODMCHeader.h"
  38 #include "AliAODHandler.h"
  39 #include "AliAnalysisManager.h"
  40 #include "AliLog.h"
  41 #include "AliEmcalJet.h"
  42 #include "AliJetContainer.h"
  43 #include "AliAODRecoDecay.h"
  44 #include "AliAODRecoCascadeHF.h"
  45 #include "AliAODRecoDecayHF2Prong.h"
  46 #include "AliESDtrack.h"
  47 #include "AliAODMCParticle.h"
  48 #include "AliPicoTrack.h"
  49 #include "AliRDHFCutsD0toKpi.h"
  50 #include "AliRDHFCutsDStartoKpipi.h"
  51 #include "AliRhoParameter.h"
  52 #include "AliParticleContainer.h"
  53
  54 ClassImp(AliAnalysisTaskFlavourJetCorrelations)
  55
  56
  57 //_______________________________________________________________________________
  58
  59 AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations() :
  60 AliAnalysisTaskEmcalJet("",kTRUE),
  61 fUseMCInfo(kTRUE),
  62 fUseReco(kTRUE),
  63 fCandidateType(),
  64 fPDGmother(),
  65 fNProngs(),
  66 fPDGdaughters(),
  67 fBranchName(),
  68 fCuts(0),
  69 fMinMass(),
  70 fMaxMass(),
  71 fJetArrName(0),
  72 fTrackArrName(0),
  73 fCandArrName(0),
  74 fLeadingJetOnly(kFALSE),
  75 fJetRadius(0),
  76 fCandidateArray(0),
  77 fSideBandArray(0),
  78 fJetOnlyMode(0),
  79 fPmissing(),
  80 fPtJet(0),
  81 fIsDInJet(0),
  82 fTypeDInJet(2),
  83 fTrackArr(0),
  84 fSwitchOnSB(0),
  85 fSwitchOnPhiAxis(0),
  86 fSwitchOnOutOfConeAxis(0),
  87 fSwitchOnSparses(1),
  88 fNAxesBigSparse(9),
  89 fJetCont(0),
  90 fTrackCont(0),
  91 fClusterCont(0),
  92 fhstat(),
  93 fhPtJetTrks(),
  94 fhPhiJetTrks(),
  95 fhEtaJetTrks(),
  96 fhEjetTrks(),
  97 fhPtJet(),
  98 fhPhiJet(),
  99 fhEtaJet(),
 100 fhEjet(),
 101 fhNtrArr(),
 102 fhNJetPerEv(),
 103 fhdeltaRJetTracks(),
 104 fhsJet(),
 105 fhNDPerEvNoJet(),
 106 fhptDPerEvNoJet(),
 107 fhNJetPerEvNoD(),
 108 fhPtJetPerEvNoD(),
 109 fhsDstandalone(),
 110 fhInvMassptD(),
 111 fhDiffSideBand(),
 112 fhInvMassptDbg(),
 113 fhPtPion(),
 114 fhztracksinjet(),
 115 fhDiffPtTrPtJzNok(),
 116 fhDiffPtTrPtJzok(),
 117 fhDiffPzTrPtJzok(),
 118 fhDiffPzTrPtJzNok(),
 119 fhNtrkjzNok(),
 120 fhztracksinjetT(),
 121 fhControlDInJ(),
 122 fhIDddaugh(),
 123 fhIDddaughOut(),
 124 fhIDjetTracks(),
 125 fhDRdaughOut(),
 126 fhzDinjet(),
 127 fhmissingp(),
 128 fhMissPi(),
 129 fhDeltaPtJet(),
 130 fhRelDeltaPtJet(),
 131 fhzDT(),
 132 fhDeltaRD(),
 133 fhDeltaRptDptj(),
 134 fhDeltaRptDptjB(),
 135 fhsDphiz()
 136 {
 137    //
 138    // Default ctor
 139    //
 140    //SetMakeGeneralHistograms(kTRUE)(),
 141
 142 }
 143
 144 //_______________________________________________________________________________
 145
 146 AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations(const Char_t* name, AliRDHFCuts* cuts,ECandidateType candtype, Bool_t jetOnly) :
 147 AliAnalysisTaskEmcalJet(name,kTRUE),
 148 fUseMCInfo(kTRUE),
 149 fUseReco(kTRUE),
 150 fCandidateType(),
 151 fPDGmother(),
 152 fNProngs(),
 153 fPDGdaughters(),
 154 fBranchName(),
 155 fCuts(0),
 156 fMinMass(),
 157 fMaxMass(),
 158 fJetArrName(0),
 159 fTrackArrName(0),
 160 fCandArrName(0),
 161 fLeadingJetOnly(kFALSE),
 162 fJetRadius(0),
 163 fCandidateArray(0),
 164 fSideBandArray(0),
 165 fJetOnlyMode(jetOnly),
 166 fPmissing(),
 167 fPtJet(0),
 168 fIsDInJet(0),
 169 fTypeDInJet(2),
 170 fTrackArr(0),
 171 fSwitchOnSB(0),
 172 fSwitchOnPhiAxis(0),
 173 fSwitchOnOutOfConeAxis(0),
 174 fSwitchOnSparses(1),
 175 fNAxesBigSparse(9),
 176 fJetCont(0),
 177 fTrackCont(0),
 178 fClusterCont(0),
 179 fhstat(),
 180 fhPtJetTrks(),
 181 fhPhiJetTrks(),
 182 fhEtaJetTrks(),
 183 fhEjetTrks(),
 184 fhPtJet(),
 185 fhPhiJet(),
 186 fhEtaJet(),
 187 fhEjet(),
 188 fhNtrArr(),
 189 fhNJetPerEv(),
 190 fhdeltaRJetTracks(),
 191 fhsJet(),
 192 fhNDPerEvNoJet(),
 193 fhptDPerEvNoJet(),
 194 fhNJetPerEvNoD(),
 195 fhPtJetPerEvNoD(),
 196 fhsDstandalone(),
 197 fhInvMassptD(),
 198 fhDiffSideBand(),
 199 fhInvMassptDbg(),
 200 fhPtPion(),
 201 fhztracksinjet(),
 202 fhDiffPtTrPtJzNok(),
 203 fhDiffPtTrPtJzok(),
 204 fhDiffPzTrPtJzok(),
 205 fhDiffPzTrPtJzNok(),
 206 fhNtrkjzNok(),
 207 fhztracksinjetT(),
 208 fhControlDInJ(),
 209 fhIDddaugh(),
 210 fhIDddaughOut(),
 211 fhIDjetTracks(),
 212 fhDRdaughOut(),
 213 fhzDinjet(),
 214 fhmissingp(),
 215 fhMissPi(),
 216 fhDeltaPtJet(),
 217 fhRelDeltaPtJet(),
 218 fhzDT(),
 219 fhDeltaRD(),
 220 fhDeltaRptDptj(),
 221 fhDeltaRptDptjB(),
 222 fhsDphiz()
 223 {
 224    //
 225    // Constructor. Initialization of Inputs and Outputs
 226    //
 227
 228    Info("AliAnalysisTaskFlavourJetCorrelations","Calling Constructor");
 229    fCuts=cuts;
 230    fCandidateType=candtype;
 231    const Int_t nptbins=fCuts->GetNPtBins();
 232    Float_t defaultSigmaD013[13]={0.012, 0.012, 0.012, 0.015, 0.015,0.018,0.018,0.020,0.020,0.030,0.030,0.037,0.040};
 233
 234    switch(fCandidateType){
 235    case 0:
 236       fPDGmother=421;
 237       fNProngs=2;
 238       fPDGdaughters[0]=211;//pi
 239       fPDGdaughters[1]=321;//K
 240       fPDGdaughters[2]=0; //empty
 241       fPDGdaughters[3]=0; //empty
 242       fBranchName="D0toKpi";
 243       fCandArrName="D0";
 244       break;
 245    case 1:
 246       fPDGmother=413;
 247       fNProngs=3;
 248       fPDGdaughters[1]=211;//pi soft
 249       fPDGdaughters[0]=421;//D0
 250       fPDGdaughters[2]=211;//pi fromD0
 251       fPDGdaughters[3]=321; // K from D0
 252       fBranchName="Dstar";
 253       fCandArrName="Dstar";
 254
 255       if(nptbins<=13){
 256          for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]= defaultSigmaD013[ipt];
 257       } else {
 258          AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
 259       }
 260       break;
 261    default:
 262       printf("%d not accepted!!\n",fCandidateType);
 263       break;
 264    }
 265
 266    if(fCandidateType==kD0toKpi)SetMassLimits(0.15,fPDGmother);
 267    if(fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
 268
 269    if(fJetOnlyMode){
 270       DefineOutput(1,TList::Class()); //histos with jet info
 271       DefineOutput(2,AliRDHFCuts::Class()); // my cuts
 272    }
 273    else{
 274       DefineInput(1, TClonesArray::Class());
 275       DefineInput(2, TClonesArray::Class());
 276
 277       DefineOutput(1,TList::Class()); // histos
 278       DefineOutput(2,AliRDHFCuts::Class()); // my cuts
 279    }
 280 }
 281
 282 //_______________________________________________________________________________
 283
 284 AliAnalysisTaskFlavourJetCorrelations::~AliAnalysisTaskFlavourJetCorrelations() {
 285    //
 286    // destructor
 287    //
 288
 289    Info("~AliAnalysisTaskFlavourJetCorrelations","Calling Destructor");
 290
 291    delete fCuts;
 292
 293 }
 294
 295 //_______________________________________________________________________________
 296
 297 void AliAnalysisTaskFlavourJetCorrelations::Init(){
 298    //
 299    // Initialization
 300    //
 301
 302    if(fDebug > 1) printf("AnalysisTaskRecoJetCorrelations::Init() \n");
 303
 304    switch(fCandidateType){
 305    case 0:
 306       {
 307          AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
 308          copyfCuts->SetName("AnalysisCutsDzero");
 309          // Post the data
 310          PostData(2,copyfCuts);
 311       }
 312       break;
 313    case 1:
 314       {
 315          AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
 316          copyfCuts->SetName("AnalysisCutsDStar");
 317          // Post the cuts
 318          PostData(2,copyfCuts);
 319       }
 320       break;
 321    default:
 322       return;
 323    }
 324
 325    return;
 326 }
 327
 328 //_______________________________________________________________________________
 329
 330 void AliAnalysisTaskFlavourJetCorrelations::UserCreateOutputObjects() {
 331    // output
 332    Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 333    AliAnalysisTaskEmcal::UserCreateOutputObjects();
 334
 335    fJetCont = GetJetContainer(0);
 336    if(fJetCont){
 337       fTrackCont =   fJetCont->GetParticleContainer();
 338       fClusterCont = fJetCont->GetClusterContainer();
 339    }
 340
 341
 342    // define histograms
 343    // the TList fOutput is already defined in  AliAnalysisTaskEmcal::UserCreateOutputObjects()
 344    DefineHistoForAnalysis();
 345    PostData(1,fOutput);
 346
 347    return;
 348 }
 349
 350 //_______________________________________________________________________________
 351
 352 Bool_t AliAnalysisTaskFlavourJetCorrelations::Run()
 353 {
 354    // user exec from AliAnalysisTaskEmcal is used
 355
 356    // Load the event
 357    AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 358
 359    TClonesArray *arrayDStartoD0pi=0;
 360
 361    if (!aodEvent && AODEvent() && IsStandardAOD()) {
 362
 363       // In case there is an AOD handler writing a standard AOD, use the AOD
 364       // event in memory rather than the input (ESD) event.
 365       aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
 366
 367       // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 368       // have to taken from the AOD event hold by the AliAODExtension
 369       AliAODHandler* aodHandler = (AliAODHandler*)
 370       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 371       if(aodHandler->GetExtensions()) {
 372          AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 373          AliAODEvent *aodFromExt = ext->GetAOD();
 374          arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
 375       }
 376    } else {
 377       arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
 378    }
 379
 380    if (!arrayDStartoD0pi) {
 381       AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
 382       //  return;
 383    } else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
 384
 385    TClonesArray* mcArray = 0x0;
 386    if (fUseMCInfo) { //not used at the moment,uncomment return if you use
 387       mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 388       if (!mcArray) {
 389          printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
 390          //return kFALSE;
 391       }
 392    }
 393
 394    //retrieve jets
 395    //this is a duplication of fTrackCont, but is is used in the loop of line 598 and changing it needs a thorough test
 396    fTrackArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTrackArrName));
 397    //clusArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClustersCorr"));
 398    //fJetArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrName));
 399    //fJetContainer=GetJetContainer(0);
 400    //if(!fJetContainer) {
 401    //   AliError("Jet Container 0 not found");
 402    //   return kFALSE;
 403    //}
 404    fJetRadius=fJetCont->GetJetRadius();
 405
 406    if(!fTrackArr){
 407       AliInfo(Form("Could not find the track array\n"));
 408       return kFALSE;
 409    }
 410
 411
 412    fCandidateArray = dynamic_cast<TClonesArray*>(GetInputData(1));
 413    if (!fCandidateArray) return kFALSE;
 414    if ((fCandidateType==1 && fSwitchOnSB) || fUseMCInfo) {
 415       fSideBandArray = dynamic_cast<TClonesArray*>(GetInputData(2));
 416       if (!fSideBandArray) return kFALSE;
 417    }
 418    //Printf("ncandidates found %d",fCandidateArray->GetEntriesFast());
 419
 420    fhstat->Fill(0);
 421
 422    // fix for temporary bug in ESDfilter
 423    // the AODs with null vertex pointer didn't pass the PhysSel
 424    if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return kFALSE;
 425
 426    //Event selection
 427    Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
 428    TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
 429    if(!iseventselected) return kFALSE;
 430
 431    fhstat->Fill(1);
 432
 433    //retrieve charm candidates selected
 434    Int_t candidates = 0;
 435    Int_t njets=fJetCont->GetNJets();
 436    //Printf("N jets in this event %d",njets);
 437    if(!fJetOnlyMode) {
 438       candidates = fCandidateArray->GetEntriesFast();
 439
 440    //trigger on jets
 441    if(njets == 0) {
 442       fhstat->Fill(6, candidates);
 443       fhNDPerEvNoJet->Fill(candidates);
 444       for(Int_t iD=0;iD<candidates;iD++){
 445          AliVParticle* cand=(AliVParticle*)fCandidateArray->At(iD);
 446          if(!cand) continue;
 447          fhptDPerEvNoJet->Fill(cand->Pt());
 448
 449       }
 450       return kFALSE;
 451
 452    }
 453
 454    //loop on candidates standalone (checking the candidates are there and their phi-eta distributions)
 455
 456    for(Int_t ic = 0; ic < candidates; ic++) {
 457
 458       // D* candidates
 459       AliAODRecoDecayHF* charm=0x0;
 460       AliAODRecoCascadeHF* dstar=0x0;
 461
 462
 463       charm=(AliAODRecoDecayHF*)fCandidateArray->At(ic);
 464       if(!charm) continue;
 465
 466       if(fCandidateType==kDstartoKpipi){
 467          dstar=(AliAODRecoCascadeHF*)fCandidateArray->At(ic);
 468          if(!dstar) continue;
 469       }
 470
 471       fhstat->Fill(2);
 472
 473       Double_t candsparse[4]={charm->Eta(), charm->Phi(), charm->Pt(), 0};
 474
 475       if(fCandidateType==kDstartoKpipi) {
 476          if(fUseReco){
 477             Double_t deltamass= dstar->DeltaInvMass();
 478             candsparse[3]=deltamass;
 479          } else candsparse[3] = 0.145; //for generated
 480          if(fSwitchOnSparses) fhsDstandalone->Fill(candsparse);
 481       }
 482       if(fCandidateType==kD0toKpi){
 483          if(fUseReco){
 484             AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)charm;
 485             Int_t isselected=fCuts->IsSelected(dzero,AliRDHFCuts::kAll,aodEvent);
 486             //not a further selection,just to get the value of isselected for the D0
 487             Double_t masses[2];
 488             Int_t pdgdaughtersD0[2]={211,321};//pi,K
 489             Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
 490
 491             masses[0]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
 492             masses[1]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
 493             if(isselected==1 || isselected==3) {
 494                candsparse[3]=masses[0];
 495                if(fSwitchOnSparses) fhsDstandalone->Fill(candsparse);
 496             }
 497             if(isselected>=2){
 498                candsparse[3]=masses[1];
 499                if(fSwitchOnSparses) fhsDstandalone->Fill(candsparse);
 500
 501             }
 502          } else { //generated
 503             Int_t pdg=((AliAODMCParticle*)charm)->GetPdgCode();
 504             candsparse[3]=TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
 505             if(fSwitchOnSparses) fhsDstandalone->Fill(candsparse);
 506          }
 507       }
 508    }
 509    }
 510
 511     //Background Subtraction for jets
 512     //If there's no background subtraction rhoval=0 and momentum is simply not took into account
 513     AliRhoParameter *rho = 0;
 514     Double_t rhoval = 0;
 515     TString sname("Rho");
 516     if (!sname.IsNull()) {
 517         rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(sname));
 518         if(rho) rhoval = rho->GetVal();
 519     }
 520
 521
 522    // we start with jets
 523    Double_t ejet   = 0;
 524    Double_t phiJet = 0;
 525    Double_t etaJet = 0;
 526    Double_t ptjet = 0;
 527    Double_t leadingJet =0;
 528    Double_t pointJ[6];
 529
 530    Int_t ntrarr=fTrackCont->GetNParticles();
 531    fhNtrArr->Fill(ntrarr);
 532
 533    for(Int_t i=0;i<ntrarr;i++){
 534       AliVTrack *jtrack=static_cast<AliVTrack*>(fTrackCont->GetParticle(i));
 535       //reusing histograms
 536       if(!jtrack) continue;
 537       fhPtJetTrks->Fill(jtrack->Pt());
 538       fhPhiJetTrks->Fill(jtrack->Phi());
 539       fhEtaJetTrks->Fill(jtrack->Eta());
 540       fhEjetTrks->Fill(jtrack->E());
 541    }
 542
 543
 544    //option to use only the leading jet
 545    if(fLeadingJetOnly){
 546       for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) {
 547          AliEmcalJet* jetL = (AliEmcalJet*)fJetCont->GetJet(iJetsL);
 548          ptjet   = jetL->Pt() - jetL->Area()*rhoval; //It takes into account the background subtraction
 549          if(ptjet>leadingJet ) leadingJet = ptjet;
 550
 551       }
 552    }
 553
 554    Int_t cntjet=0;
 555    //loop over jets and consider only the leading jet in the event
 556    for (Int_t iJets = 0; iJets<njets; iJets++) {
 557       fPmissing[0]=0;
 558       fPmissing[1]=0;
 559       fPmissing[2]=0;
 560
 561       //Printf("Jet N %d",iJets);
 562       AliEmcalJet* jet = (AliEmcalJet*)fJetCont->GetJet(iJets);
 563       //Printf("Corr task Accept Jet");
 564       if(!AcceptJet(jet)) {
 565          fhstat->Fill(5);
 566          continue;
 567       }
 568
 569       //jets variables
 570       ejet   = jet->E();
 571       phiJet = jet->Phi();
 572       etaJet = jet->Eta();
 573       fPtJet = jet->Pt() - jet->Area()*rhoval; //It takes into account the background subtraction
 574       Double_t origPtJet=fPtJet;
 575
 576       pointJ[0] = phiJet;
 577       pointJ[1] = etaJet;
 578       pointJ[2] = ptjet - jet->Area()*rhoval; //It takes into account the background subtraction
 579       pointJ[3] = ejet;
 580       pointJ[4] = jet->GetNumberOfConstituents();
 581       pointJ[5] = jet->Area();
 582
 583       // choose the leading jet
 584       if(fLeadingJetOnly && (ejet<leadingJet)) continue;
 585       //used for normalization
 586       fhstat->Fill(3);
 587       cntjet++;
 588       // fill energy, eta and phi of the jet
 589       fhEjet   ->Fill(ejet);
 590       fhPhiJet ->Fill(phiJet);
 591       fhEtaJet ->Fill(etaJet);
 592       fhPtJet  ->Fill(fPtJet);
 593       if(fJetOnlyMode && fSwitchOnSparses) fhsJet->Fill(pointJ,1);
 594       //loop on jet particles
 595       Int_t ntrjet=  jet->GetNumberOfTracks();
 596       Double_t sumPtTracks=0,sumPzTracks=0;
 597       Int_t zg1jtrk=0;
 598       for(Int_t itrk=0;itrk<ntrjet;itrk++){
 599
 600          AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,fTrackArr);
 601          fhdeltaRJetTracks->Fill(DeltaR(jet,jetTrk));
 602          Double_t ztrackj=Z(jetTrk,jet);
 603          fhztracksinjet->Fill(ztrackj);
 604          sumPtTracks+=jetTrk->Pt();
 605          sumPzTracks+=jetTrk->Pz();
 606          if(ztrackj>1){
 607             zg1jtrk++;
 608          }
 609
 610          Double_t ztrackjTr=Z(jetTrk,jet,kTRUE);
 611          fhztracksinjetT->Fill(ztrackjTr);
 612
 613       }//end loop on jet tracks
 614
 615       fhNtrkjzNok->Fill(zg1jtrk);
 616       Double_t diffpt=origPtJet-sumPtTracks;
 617       Double_t diffpz=jet->Pz()-sumPzTracks;
 618       if(zg1jtrk>0){
 619          fhDiffPtTrPtJzNok->Fill(diffpt);
 620          fhDiffPzTrPtJzNok->Fill(diffpz);
 621
 622       }else{
 623          fhDiffPtTrPtJzok->Fill(diffpt);
 624          fhDiffPzTrPtJzok->Fill(diffpz);
 625       }
 626
 627       if(candidates==0){
 628
 629          fhPtJetPerEvNoD->Fill(fPtJet);
 630       }
 631       if(!fJetOnlyMode) {
 632          //Printf("N candidates %d ", candidates);
 633          for(Int_t ic = 0; ic < candidates; ic++) {
 634             fhstat->Fill(7);
 635             // D* candidates
 636             AliVParticle* charm=0x0;
 637             charm=(AliVParticle*)fCandidateArray->At(ic);
 638             if(!charm) continue;
 639             AliAODRecoDecayHF *charmdecay=(AliAODRecoDecayHF*) charm;
 640             fIsDInJet=IsDInJet(jet, charmdecay, kTRUE);
 641             if (fIsDInJet) FlagFlavour(jet);
 642             if (jet->TestFlavourTag(AliEmcalJet::kDStar) || jet->TestFlavourTag(AliEmcalJet::kD0)) fhstat->Fill(4);
 643
 644             //Note: the z component of the jet momentum comes from the eta-phi direction of the jet particles, it is not calculated from the z component of the tracks since, as default, the scheme used for jet reco is the pt-scheme which sums the scalar component, not the vectors. Addind the D daughter momentum component by componet as done here is not 100% correct, but the difference is small, for fairly collimated particles.
 645
 646             Double_t pjet[3];
 647             jet->PxPyPz(pjet);
 648              //It corrects the jet momentum if it was asked for jet background subtraction
 649              pjet[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
 650              pjet[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
 651              pjet[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
 652
 653             //It corrects the jet momentum due to D daughters out of the jet
 654             RecalculateMomentum(pjet,fPmissing);
 655             fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]); //recalculated jet pt
 656
 657
 658             //debugging histograms
 659             Double_t pmissing=TMath::Sqrt(fPmissing[0]*fPmissing[0]+fPmissing[1]*fPmissing[1]+fPmissing[2]*fPmissing[2]); //recalculated jet total momentum
 660             for(Int_t k=0;k<3;k++) fhMissPi[k]->Fill(fPmissing[k]);
 661
 662             fhmissingp->Fill(pmissing);
 663             Double_t ptdiff=fPtJet-origPtJet;
 664             fhDeltaPtJet->Fill(ptdiff);
 665             fhRelDeltaPtJet->Fill(ptdiff/origPtJet);
 666
 667             FillHistogramsRecoJetCorr(charm, jet, aodEvent);
 668
 669          }//end loop on candidates
 670
 671          //retrieve side band background candidates for Dstar
 672          Int_t nsbcand = 0; if(fSideBandArray) nsbcand=fSideBandArray->GetEntriesFast();
 673
 674          for(Int_t ib=0;ib<nsbcand;ib++){
 675             if(fCandidateType==kDstartoKpipi && !fUseMCInfo){
 676                AliAODRecoCascadeHF *sbcand=(AliAODRecoCascadeHF*)fSideBandArray->At(ib);
 677                if(!sbcand) continue;
 678
 679                fIsDInJet=IsDInJet(jet, sbcand,kFALSE);
 680                Double_t pjet[3];
 681                jet->PxPyPz(pjet);
 682                 //It corrects the jet momentum if it was asked for jet background subtraction
 683                 pjet[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
 684                 pjet[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
 685                 pjet[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
 686
 687                //It corrects the jet momentum due to D daughters out of the jet
 688                RecalculateMomentum(pjet,fPmissing);
 689                fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]); //recalculated jet pt
 690
 691                SideBandBackground(sbcand,jet);
 692
 693             }
 694             if(fUseMCInfo){
 695
 696                AliAODRecoDecayHF* charmbg = 0x0;
 697                charmbg=(AliAODRecoDecayHF*)fSideBandArray->At(ib);
 698                if(!charmbg) continue;
 699                fhstat->Fill(8);
 700                fIsDInJet=IsDInJet(jet, charmbg,kFALSE);
 701                if (fIsDInJet) {
 702                   FlagFlavour(jet); //this are backgroud HF jets, but flagged as signal at the moment. Can use the bkg flavour flag in the future. This info is not stored now a part in the jet
 703                   fhstat->Fill(9);
 704                }
 705                Double_t pjet[3];
 706                jet->PxPyPz(pjet);
 707                //It corrects the jet momentum if it was asked for jet background subtraction
 708                pjet[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
 709                pjet[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
 710                pjet[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
 711
 712                //It corrects the jet momentum due to D daughters out of the jet
 713                RecalculateMomentum(pjet,fPmissing);
 714                fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]); //recalculated jet pt
 715
 716                MCBackground(charmbg,jet);
 717             }
 718          }
 719       }
 720    } // end of jet loop
 721
 722    fhNJetPerEv->Fill(cntjet);
 723    if(candidates==0) fhNJetPerEvNoD->Fill(cntjet);
 724    PostData(1,fOutput);
 725    return kTRUE;
 726 }
 727
 728 //_______________________________________________________________________________
 729
 730 void AliAnalysisTaskFlavourJetCorrelations::Terminate(Option_t*)
 731 {
 732    // The Terminate() function is the last function to be called during
 733    // a query. It always runs on the client, it can be used to present
 734    // the results graphically or save the results to file.
 735
 736    Info("Terminate"," terminate");
 737    AliAnalysisTaskSE::Terminate();
 738
 739    fOutput = dynamic_cast<TList*> (GetOutputData(1));
 740    if (!fOutput) {
 741       printf("ERROR: fOutput not available\n");
 742       return;
 743    }
 744 }
 745
 746 //_______________________________________________________________________________
 747
 748 void  AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){
 749    Float_t mass=0;
 750    Int_t abspdg=TMath::Abs(pdg);
 751
 752    mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
 753    // compute the Delta mass for the D*
 754    if(fCandidateType==kDstartoKpipi){
 755       Float_t mass1=0;
 756       mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 757       mass = mass-mass1;
 758    }
 759
 760    fMinMass = mass-range;
 761    fMaxMass = mass+range;
 762
 763    AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
 764    if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
 765 }
 766
 767 //_______________________________________________________________________________
 768
 769 void  AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){
 770    if(uplimit>lowlimit)
 771    {
 772       fMinMass = lowlimit;
 773       fMaxMass = uplimit;
 774    }
 775    else{
 776       printf("Error! Lower limit larger than upper limit!\n");
 777       fMinMass = uplimit - uplimit*0.2;
 778       fMaxMass = uplimit;
 779    }
 780 }
 781
 782 //_______________________________________________________________________________
 783
 784 Bool_t AliAnalysisTaskFlavourJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){
 785    if(nptbins>30) {
 786       AliInfo("Maximum number of bins allowed is 30!");
 787       return kFALSE;
 788    }
 789    if(!width) return kFALSE;
 790    for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
 791    return kTRUE;
 792 }
 793
 794 //_______________________________________________________________________________
 795
 796 Double_t AliAnalysisTaskFlavourJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet, Bool_t transverse) const{
 797    if(!part || !jet){
 798       printf("Particle or jet do not exist!\n");
 799       return -99;
 800    }
 801    Double_t p[3],pj[3];
 802    Bool_t okpp=part->PxPyPz(p);
 803    Bool_t okpj=jet->PxPyPz(pj);
 804
 805     //Background Subtraction
 806     //It corrects the each component of the jet momentum for Z calculation
 807     AliRhoParameter *rho = 0;
 808     Double_t rhoval = 0;
 809     TString sname("Rho");
 810     if (!sname.IsNull()) {
 811         rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(sname));
 812         if(rho){
 813             rhoval = rho->GetVal();
 814             pj[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
 815             pj[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
 816             pj[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
 817         }
 818     }
 819
 820
 821
 822    if(!okpp || !okpj){
 823       printf("Problems getting momenta\n");
 824       return -999;
 825    }
 826
 827    RecalculateMomentum(pj, fPmissing);
 828    if(transverse) return ZT(p,pj);
 829    else
 830    return Z(p,pj);
 831 }
 832
 833 //_______________________________________________________________________________
 834 Double_t AliAnalysisTaskFlavourJetCorrelations::Z(Double_t* p, Double_t *pj) const{
 835
 836    Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1]+pj[2]*pj[2];
 837    Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet2);
 838    return z;
 839 }
 840
 841
 842 //_______________________________________________________________________________
 843 Double_t AliAnalysisTaskFlavourJetCorrelations::ZT(Double_t* p, Double_t *pj) const{
 844
 845    Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1];
 846    Double_t z=(p[0]*pj[0]+p[1]*pj[1])/(pjet2);
 847    return z;
 848 }
 849
 850 //_______________________________________________________________________________
 851
 852 void AliAnalysisTaskFlavourJetCorrelations::RecalculateMomentum(Double_t* pj, const Double_t* pmissing) const {
 853
 854    pj[0]+=pmissing[0];
 855    pj[1]+=pmissing[1];
 856    pj[2]+=pmissing[2];
 857
 858 }
 859
 860 //_______________________________________________________________________________
 861
 862 Bool_t  AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){
 863
 864    // Statistics
 865    Int_t nbins=8;
 866    if(fUseMCInfo) nbins+=2;
 867
 868    fhstat=new TH1I("hstat","Statistics",nbins,-0.5,nbins-0.5);
 869    fhstat->GetXaxis()->SetBinLabel(1,"N ev anal");
 870    fhstat->GetXaxis()->SetBinLabel(2,"N ev sel");
 871    fhstat->GetXaxis()->SetBinLabel(3,"N cand sel");
 872    fhstat->GetXaxis()->SetBinLabel(4,"N jets");
 873    fhstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
 874    fhstat->GetXaxis()->SetBinLabel(6,"N jet rej");
 875    fhstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet");
 876    fhstat->GetXaxis()->SetBinLabel(8,"N jets & cand");
 877    if(fUseMCInfo) {
 878     fhstat->GetXaxis()->SetBinLabel(3,"N Signal sel & jet");
 879     fhstat->GetXaxis()->SetBinLabel(5,"N Signal in jet");
 880     fhstat->GetXaxis()->SetBinLabel(9,"N Bkg sel & jet");
 881     fhstat->GetXaxis()->SetBinLabel(10,"N Bkg in jet");
 882    }
 883    fhstat->SetNdivisions(1);
 884    fOutput->Add(fhstat);
 885
 886    const Int_t nbinsmass=200;
 887    const Int_t nbinsptjet=500;
 888    const Int_t nbinsptD=100;
 889    const Int_t nbinsz=100;
 890    const Int_t nbinsphi=200;
 891    const Int_t nbinseta=100;
 892
 893    //binning for THnSparse
 894    const Int_t nbinsSpsmass=50;
 895    const Int_t nbinsSpsptjet=100;
 896    const Int_t nbinsSpsptD=50;
 897    const Int_t nbinsSpsz=100;
 898    const Int_t nbinsSpsphi=100;
 899    const Int_t nbinsSpseta=60;
 900    const Int_t nbinsSpsContrib=100;
 901    const Int_t nbinsSpsA=100;
 902
 903    const Float_t ptjetlims[2]={0.,200.};
 904    const Float_t ptDlims[2]={0.,50.};
 905    const Float_t zlims[2]={0.,1.2};
 906    const Float_t philims[2]={0.,6.3};
 907    const Float_t etalims[2]={-1.5,1.5};
 908    const Int_t   nContriblims[2]={0,100};
 909    const Float_t arealims[2]={0.,2};
 910
 911    // jet related fistograms
 912
 913    fhEjetTrks      = new TH1F("hEjetTrks",  "Jet tracks energy distribution;Energy (GeV)",500,0,200);
 914    fhEjetTrks->Sumw2();
 915    fhPhiJetTrks    = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)",  nbinsphi,philims[0],philims[1]);
 916    fhPhiJetTrks->Sumw2();
 917    fhEtaJetTrks    = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta",  nbinseta,etalims[0],etalims[1]);
 918    fhEtaJetTrks->Sumw2();
 919    fhPtJetTrks     = new TH1F("hPtJetTrks",  "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
 920    fhPtJetTrks->Sumw2();
 921
 922    fhEjet      = new TH1F("hEjet",  "Jet energy distribution;Energy (GeV)",500,0,200);
 923    fhEjet->Sumw2();
 924    fhPhiJet    = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)",  nbinsphi,philims[0],philims[1]);
 925    fhPhiJet->Sumw2();
 926    fhEtaJet    = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
 927    fhEtaJet->Sumw2();
 928    fhPtJet      = new TH1F("hPtJet",  "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
 929    fhPtJet->Sumw2();
 930
 931    fhdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.);
 932    fhdeltaRJetTracks->Sumw2();
 933
 934    fhNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000);
 935    fhNtrArr->Sumw2();
 936
 937    fhNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);
 938    fhNJetPerEv->Sumw2();
 939
 940
 941    fOutput->Add(fhEjetTrks);
 942    fOutput->Add(fhPhiJetTrks);
 943    fOutput->Add(fhEtaJetTrks);
 944    fOutput->Add(fhPtJetTrks);
 945    fOutput->Add(fhEjet);
 946    fOutput->Add(fhPhiJet);
 947    fOutput->Add(fhEtaJet);
 948    fOutput->Add(fhPtJet);
 949    fOutput->Add(fhdeltaRJetTracks);
 950    fOutput->Add(fhNtrArr);
 951    fOutput->Add(fhNJetPerEv);
 952
 953    if(fJetOnlyMode && fSwitchOnSparses){
 954       //thnsparse for jets
 955       const Int_t nAxis=6;
 956       const Int_t nbinsSparse[nAxis]={nbinsSpsphi,nbinsSpseta, nbinsSpsptjet, nbinsSpsptjet,nbinsSpsContrib, nbinsSpsA};
 957       const Double_t minSparse[nAxis]={philims[0],etalims[0],ptjetlims[0],ptjetlims[0],static_cast<Double_t>(nContriblims[0]),arealims[0]};
 958       const Double_t
 959         maxSparse[nAxis]={philims[1],etalims[1],ptjetlims[1],ptjetlims[1],static_cast<Double_t>(nContriblims[1]),arealims[1]};
 960       fhsJet=new THnSparseF("hsJet","#phi, #eta, p_{T}^{jet}, E^{jet}, N contrib, Area", nAxis, nbinsSparse, minSparse, maxSparse);
 961       fhsJet->Sumw2();
 962
 963       fOutput->Add(fhsJet);
 964
 965    }
 966
 967    if(!fJetOnlyMode){
 968
 969       //debugging histograms
 970       fhControlDInJ=new TH1I("hControlDInJ","Checks D in Jet",8, -0.5,7.5);
 971       fhControlDInJ->GetXaxis()->SetBinLabel(1,"DR In,1 daugh out");
 972       fhControlDInJ->GetXaxis()->SetBinLabel(2,"DR In,2 daugh out");
 973       fhControlDInJ->GetXaxis()->SetBinLabel(3,"DR In,3 daugh out");
 974       fhControlDInJ->GetXaxis()->SetBinLabel(4,"Tot tracks non matched");
 975       fhControlDInJ->GetXaxis()->SetBinLabel(5,"D0 daug missing");
 976       fhControlDInJ->GetXaxis()->SetBinLabel(6,"soft pi missing");
 977       fhControlDInJ->GetXaxis()->SetBinLabel(7,"IDprong diff GetDau");
 978       fhControlDInJ->GetXaxis()->SetBinLabel(8,"still z>1 (cand)");
 979
 980       fhControlDInJ->SetNdivisions(1);
 981       fhControlDInJ->GetXaxis()->SetLabelSize(0.05);
 982       fOutput->Add(fhControlDInJ);
 983
 984       fhmissingp=new TH1F("hmissingp", "Distribution of missing momentum (modulo);|p_{missing}|", 200, 0.,20.);
 985       fOutput->Add(fhmissingp);
 986
 987       fhMissPi=new TH1F*[3];
 988       for(Int_t k=0;k<3;k++) {
 989          fhMissPi[k]=new TH1F(Form("hMissP%d",k), Form("Missing p component {%d};p_{%d}",k,k), 400, -10.,10.);
 990          fOutput->Add(fhMissPi[k]);
 991       }
 992       fhDeltaPtJet=new TH1F("hDeltaPtJet", "Difference between the jet pt before and after correction;p_{T}^{jet,recalc}-p_{T}^{jet,orig}", 200, 0.,20.);
 993
 994       fOutput->Add(fhDeltaPtJet);
 995       fhRelDeltaPtJet=new TH1F("hRelDeltaPtJet", "Difference between the jet pt before and after correction/ original pt;(p_{T}^{jet,recalc}-p_{T}^{jet,orig})/p_{T}^{jet,orig}", 200, 0.,1.);
 996       fOutput->Add(fhRelDeltaPtJet);
 997
 998       fhzDinjet=new TH1F("hzDinjet","Z of candidates with daughters in jet;z",nbinsz,zlims[0],zlims[1]);
 999       fOutput->Add(fhzDinjet);
1000       //frag func of tracks in the jet
1001       fhztracksinjet=new TH1F("hztracksinjet","Z of tracks in jet;z",nbinsz,zlims[0],zlims[1]);
1002       fOutput->Add(fhztracksinjet);
1003
1004       //check jet and tracks
1005       fhDiffPtTrPtJzok=new TH1F("hDiffPtTrPtJzok","Sum p_{T}^{trks} - p_{T}^{jet}, for z<1;(#Sigma p_{T}^{trks}) - p_{T}^{jet}", 100,-0.2,0.2);
1006       fOutput->Add(fhDiffPtTrPtJzok);
1007       fhDiffPtTrPtJzNok=new TH1F("hDiffPtTrPtJzNok","Sum p_{T}^{trks} - p_{T}^{jet}, for z>1;(#Sigma p_{T}^{trks}) - p_{T}^{jet}", 100,-0.2,0.2);
1008       fOutput->Add(fhDiffPtTrPtJzNok);
1009       fhDiffPzTrPtJzok=new TH1F("hDiffPzTrPtJzok","Sum p_{z}^{trks} - p_{z}^{jet}, for z<1;(#Sigma p_{z}^{trks}) - p_{z}^{jet}", 100,-0.2,0.2);
1010       fOutput->Add(fhDiffPzTrPtJzok);
1011       fhDiffPzTrPtJzNok=new TH1F("hDiffPzTrPtJzNok","Sum p_{z}^{trks} - p_{z}^{jet}, for z>1;(#Sigma p_{z}^{trks}) - p_{z}^{jet}", 100,-0.2,0.2);
1012       fOutput->Add(fhDiffPzTrPtJzNok);
1013       fhNtrkjzNok=new TH1I("hNtrkjzNok", "Number of tracks in a jet with z>1;N^{tracks} (z>1)",5,0,5);
1014       fOutput->Add(fhNtrkjzNok);
1015
1016       //calculate frag func with pt (simply ptD(or track)\cdot pt jet /ptjet^2)
1017       fhzDT=new TH1F("hzDT", Form("Z of D %s in jet in transverse components;(p_{T}^{D} dot p_{T}^{jet})/p_{T}^{jet}^{2} ", fUseMCInfo ? "(S+B)" : ""),nbinsz,zlims[0],zlims[1]);
1018       fOutput->Add(fhzDT);
1019       fhztracksinjetT=new TH1F("hztracksinjetT", "Z of jet tracks in transverse components;(p_{T}^{trks} dot p_{T}^{jet})/p_{T}^{jet}^{2}",nbinsz,zlims[0],zlims[1]);
1020       fOutput->Add(fhztracksinjetT);
1021
1022       fhIDddaugh=new TH1I("hIDddaugh", "ID of daughters;ID", 2001,-1000,1000);
1023       fOutput->Add(fhIDddaugh);
1024       fhIDddaughOut=new TH1I("hIDddaughOut", "ID of daughters not found in jet;ID", 2001,-1000,1000);
1025       fOutput->Add(fhIDddaughOut);
1026       fhIDjetTracks=new TH1I("hIDjetTracks", "ID of jet tracks;ID", 2001,-1000,1000);
1027       fOutput->Add(fhIDjetTracks);
1028
1029       fhDRdaughOut=new TH1F("hDRdaughOut", "#Delta R of daughters not belonging to the jet tracks (D in jet);#Delta R",200, 0.,2.);
1030       fOutput->Add(fhDRdaughOut);
1031
1032
1033       if(fCandidateType==kDstartoKpipi)
1034       {
1035          if(fSwitchOnSB){
1036             fhDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
1037             fhDiffSideBand->SetStats(kTRUE);
1038             fhDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
1039             fhDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
1040             fhDiffSideBand->Sumw2();
1041             fOutput->Add(fhDiffSideBand);
1042          }
1043
1044          fhPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
1045          fhPtPion->SetStats(kTRUE);
1046          fhPtPion->GetXaxis()->SetTitle("GeV/c");
1047          fhPtPion->GetYaxis()->SetTitle("Entries");
1048          fhPtPion->Sumw2();
1049          fOutput->Add(fhPtPion);
1050
1051       }
1052       // D related histograms
1053       fhNDPerEvNoJet=new TH1F("hNDPerEvNoJet","Number of candidates per event with no jets; N candidate/ev with no jet", 20, 0., 20.);
1054       fhNDPerEvNoJet->Sumw2();
1055       fOutput->Add(fhNDPerEvNoJet);
1056
1057       fhptDPerEvNoJet=new TH1F("hptDPerEvNoJet","pt distribution of candidates per events with no jets; p_{t}^{D} (GeV/c)",nbinsptD, ptDlims[0],ptDlims[1]);
1058       fhptDPerEvNoJet->Sumw2();
1059       fOutput->Add(fhptDPerEvNoJet);
1060
1061       if(fSwitchOnSparses){
1062          const Int_t    nAxisD=4;
1063          const Int_t    nbinsSparseD[nAxisD]={nbinsSpseta,nbinsSpsphi,nbinsSpsptD,nbinsSpsmass};
1064          const Double_t minSparseD[nAxisD]  ={etalims[0],philims[0],ptDlims[0],fMinMass};
1065          const Double_t maxSparseD[nAxisD]  ={etalims[1],philims[1],ptDlims[1],fMaxMass};
1066          fhsDstandalone=new THnSparseF("hsDstandalone","#phi, #eta, p_{T}^{D}, and mass", nAxisD, nbinsSparseD, minSparseD, maxSparseD);
1067          fhsDstandalone->Sumw2();
1068
1069          fOutput->Add(fhsDstandalone);
1070       }
1071
1072       /*
1073       TH3F* hPtJetWithD=new TH3F("hPtJetWithD","D-Jet Pt distribution; p_{T} (GeV/c);delta mass (GeV/c^{2})",nbinsptjet,ptjetlims[0],ptjetlims[1],nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
1074       hPtJetWithD->Sumw2();
1075       //for the MC this histogram is filled with the real background
1076       TH3F* hPtJetWithDsb=new TH3F("hPtJetWithDsb","D(background)-Jet Pt distribution; p_{T} (GeV/c);delta mass (GeV/c^{2});p_{T}^{D} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1],nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
1077       hPtJetWithDsb->Sumw2();
1078
1079       fOutput->Add(hPtJetWithD);
1080       fOutput->Add(hPtJetWithDsb);
1081
1082       */
1083       fhNJetPerEvNoD=new TH1F("hNJetPerEvNoD","Number of jets per event with no D; number of jets/ev with no D",10,-0.5,9.5);
1084       fhNJetPerEvNoD->Sumw2();
1085
1086       fhPtJetPerEvNoD=new TH1F("hPtJetPerEvNoD","pt distribution of jets per event with no D; #it{p}_{T}^{jet} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
1087       fhPtJetPerEvNoD->Sumw2();
1088
1089       fOutput->Add(fhNJetPerEvNoD);
1090       fOutput->Add(fhPtJetPerEvNoD);
1091
1092       fhDeltaRD=new TH1F("hDeltaRD",Form("#Delta R distribution of D candidates %s selected;#Delta R", fUseMCInfo ? "(S+B)" : ""),200, 0.,10.);
1093       fhDeltaRD->Sumw2();
1094       fOutput->Add(fhDeltaRD);
1095
1096       fhDeltaRptDptj=new TH3F("hDeltaRptDptj",Form("#Delta R distribution of D candidates %s selected;#Delta R;#it{p}_{T}^{D} (GeV/c);#it{p}_{T}^{jet} (GeV/c)", fUseMCInfo ? "(S)" : ""),100, 0.,5.,nbinsptD, ptDlims[0],ptDlims[1],nbinsptjet,ptjetlims[0],ptjetlims[1]);
1097       fhDeltaRptDptj->Sumw2();
1098       fOutput->Add(fhDeltaRptDptj);
1099
1100       if(fUseMCInfo){
1101          fhDeltaRptDptjB=new TH3F("hDeltaRptDptjB",Form("#Delta R distribution of D candidates (B) selected;#Delta R;#it{p}_{T}^{D} (GeV/c);#it{p}_{T}^{jet} (GeV/c)"),100, 0.,5.,nbinsptD, ptDlims[0],ptDlims[1],nbinsptjet,ptjetlims[0],ptjetlims[1]);
1102          fhDeltaRptDptjB->Sumw2();
1103          fOutput->Add(fhDeltaRptDptjB);
1104       }
1105
1106       //background (side bands for the Dstar and like sign for D0)
1107       AliJetContainer *jetCont=GetJetContainer(0);
1108       if(!jetCont){
1109          Printf("Container 0 not found, try with name %s", fJetArrName.Data());
1110          jetCont=GetJetContainer(fJetArrName);
1111          if(!jetCont) Printf("NOT FOUND AGAIN");
1112          return kFALSE;
1113       }
1114       Printf("CONTAINER NAME IS %s", jetCont->GetArrayName().Data());
1115       //fJetRadius=jetCont->GetJetRadius();
1116       fhInvMassptD = new TH2F("hInvMassptD",Form("D (Delta R < %.1f) invariant mass distribution p_{T}^{j} > threshold",fJetRadius),nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);
1117       fhInvMassptD->SetStats(kTRUE);
1118       fhInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
1119       fhInvMassptD->GetYaxis()->SetTitle("#it{p}_{t}^{D} (GeV/c)");
1120       fhInvMassptD->Sumw2();
1121
1122       fOutput->Add(fhInvMassptD);
1123
1124       if(fUseMCInfo){
1125          fhInvMassptDbg = new TH2F("hInvMassptDbg",Form("Background D (Delta R < %.1f) invariant mass distribution p_{T}^{j} > threshold",fJetRadius),nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);
1126          fhInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
1127          fhInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
1128          fhInvMassptDbg->Sumw2();
1129          fOutput->Add(fhInvMassptDbg);
1130
1131       }
1132
1133       if(fSwitchOnSparses){
1134          Double_t pi=TMath::Pi(), philims2[2];
1135          philims2[0]=-pi*1./2.;
1136          philims2[1]=pi*3./2.;
1137          fhsDphiz=0x0; //definition below according to the switches
1138
1139          if(fSwitchOnSB && fSwitchOnPhiAxis && fSwitchOnOutOfConeAxis){
1140             AliInfo("Creating a 9 axes container: This might cause large memory usage");
1141             const Int_t nAxis=9;
1142             const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsphi,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,2, 2, 2, 2};
1143             const Double_t minSparse[nAxis]={zlims[0],philims2[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5, -0.5,-0.5,-0.5};
1144             const Double_t maxSparse[nAxis]={zlims[1],philims2[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5, 1.5 , 1.5};
1145             fNAxesBigSparse=nAxis;
1146
1147             fhsDphiz=new THnSparseF("hsDphiz","Z and #Delta#phi vs p_{T}^{jet}, p_{T}^{D}, mass. SB? D within the jet cone?, D in EMCal acc?, jet in EMCal acc?", nAxis, nbinsSparse, minSparse, maxSparse);
1148          }
1149
1150          if(!fSwitchOnPhiAxis || !fSwitchOnOutOfConeAxis || !fSwitchOnSB){
1151             fSwitchOnPhiAxis=0;
1152             fSwitchOnOutOfConeAxis=0;
1153             fSwitchOnSB=0;
1154             if(fUseMCInfo){
1155                AliInfo("Creating a 7 axes container (MB background candidates)");
1156                const Int_t nAxis=7;
1157                const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,2, 2, 2};
1158                const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass, -0.5,-0.5,-0.5};
1159                const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5 , 1.5};
1160                fNAxesBigSparse=nAxis;
1161                fhsDphiz=new THnSparseF("hsDphiz","Z vs p_{T}^{jet}, p_{T}^{D}, mass. Bkg?, D in EMCal acc?, jet in EMCal acc?", nAxis, nbinsSparse, minSparse, maxSparse);
1162
1163             } else{
1164                AliInfo("Creating a 6 axes container");
1165                const Int_t nAxis=6;
1166                const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass, 2, 2};
1167                const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5,-0.5};
1168                const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5};
1169                fNAxesBigSparse=nAxis;
1170
1171                fhsDphiz=new THnSparseF("hsDphiz","Z vs p_{T}^{jet}, p_{T}^{D}, mass., D in EMCal acc?, jet in EMCal acc?", nAxis, nbinsSparse, minSparse, maxSparse);
1172             }
1173          }
1174          if(!fhsDphiz) AliFatal("No THnSparse created");
1175          fhsDphiz->Sumw2();
1176
1177          fOutput->Add(fhsDphiz);
1178       }
1179    }
1180    PostData(1, fOutput);
1181
1182    return kTRUE;
1183 }
1184
1185 //_______________________________________________________________________________
1186
1187 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet,  AliAODEvent* aodEvent){
1188
1189    Double_t ptD=candidate->Pt();
1190    Double_t deltaR=DeltaR(jet,candidate);
1191    Double_t phiD=candidate->Phi();
1192    Double_t deltaphi = jet->Phi()-phiD;
1193    if(deltaphi<=-(TMath::Pi())/2.) deltaphi = deltaphi+2.*(TMath::Pi());
1194    if(deltaphi>(3.*(TMath::Pi()))/2.) deltaphi = deltaphi-2.*(TMath::Pi());
1195    Double_t z=Z(candidate,jet);
1196    /*
1197    if(z>1) {
1198       ((TH1I*)fOutput->FindObject("hControlDInJ"))->Fill(7);
1199       Double_t pmissing=TMath::Sqrt(fPmissing[0]*fPmissing[0]+fPmissing[1]*fPmissing[1]+fPmissing[2]*fPmissing[2]);
1200
1201       for(Int_t k=0;k<3;k++) ((TH1F*)fOutput->FindObject(Form("hMissP%d",k)))->Fill(fPmissing[k]);
1202
1203       ((TH1F*)fOutput->FindObject("hmissingp"))->Fill(pmissing);
1204       Double_t ptdiff=fPtJet-jet->Pt();
1205       ((TH1F*)fOutput->FindObject("hDeltaPtJet"))->Fill(ptdiff);
1206       ((TH1F*)fOutput->FindObject("hRelDeltaPtJet"))->Fill(ptdiff/jet->Pt());
1207
1208
1209    }
1210    */
1211    if(fIsDInJet)fhzDT->Fill(Z(candidate,jet,kTRUE));
1212
1213    fhDeltaRD->Fill(deltaR);
1214    fhDeltaRptDptj->Fill(deltaR,ptD,fPtJet);
1215
1216    Bool_t bDInEMCalAcc=InEMCalAcceptance(candidate);
1217    Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1218    if(fUseReco){
1219       if(fCandidateType==kD0toKpi) {
1220          AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;
1221
1222          FillHistogramsD0JetCorr(dzero,deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc, aodEvent);
1223
1224       }
1225
1226       if(fCandidateType==kDstartoKpipi) {
1227          AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;
1228          FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc);
1229
1230       }
1231    } else{ //generated level
1232       //AliInfo("Non reco");
1233       FillHistogramsMCGenDJetCorr(deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc);
1234    }
1235 }
1236
1237 //_______________________________________________________________________________
1238
1239 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc, AliAODEvent* aodEvent){
1240
1241
1242    Double_t masses[2]={0.,0.};
1243    Int_t pdgdaughtersD0[2]={211,321};//pi,K
1244    Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1245
1246    masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1247    masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1248
1249    //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1250    Double_t *point=0x0;
1251    if(fNAxesBigSparse==9){
1252       point=new Double_t[9];
1253       point[0]=z;
1254       point[1]=dPhi;
1255       point[2]=ptj;
1256       point[3]=ptD;
1257       point[4]=masses[0];
1258       point[5]=0;
1259       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1260       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1261       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1262    }
1263    if(fNAxesBigSparse==6){
1264       point=new Double_t[6];
1265       point[0]=z;
1266       point[1]=ptj;
1267       point[2]=ptD;
1268       point[3]=masses[0];
1269       point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1270       point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1271 }
1272   if(fNAxesBigSparse==7){
1273       point=new Double_t[7];
1274       point[0]=z;
1275       point[1]=ptj;
1276       point[2]=ptD;
1277       point[3]=masses[0];
1278       point[4]=0;
1279       point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1280       point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1281
1282    }
1283
1284
1285    Printf("Candidate in FillHistogramsD0JetCorr IsA %s", (candidate->IsA())->GetName());
1286    Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
1287    if(isselected==1 || isselected==3) {
1288
1289       //if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[0],ptD);
1290
1291       FillMassHistograms(masses[0], ptD);
1292       if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1293    }
1294    if(isselected>=2) {
1295       //if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[1],ptD);
1296
1297       FillMassHistograms(masses[1], ptD);
1298       if(fNAxesBigSparse==9) point[4]=masses[1];
1299       if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[1];
1300       if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1301    }
1302
1303 }
1304
1305 //_______________________________________________________________________________
1306
1307 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi,  Double_t z, Double_t ptD, Double_t ptj, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc){
1308   //dPhi and z not used at the moment,but will be (re)added
1309
1310    AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
1311    Double_t deltamass= dstar->DeltaInvMass();
1312    //Double_t massD0= dstar->InvMassD0();
1313
1314
1315    fhPtPion->Fill(softpi->Pt());
1316
1317    //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1318    Int_t isSB=0;//IsDzeroSideBand(dstar);
1319
1320    //Double_t point[]={z,dPhi,ptj,ptD,deltamass,static_cast<Double_t>(isSB), static_cast<Double_t>(fIsDInJet ? 1 : 0),bDInEMCalAcc,bJetInEMCalAcc};
1321    Double_t *point=0x0;
1322    if(fNAxesBigSparse==9){
1323       point=new Double_t[9];
1324       point[0]=z;
1325       point[1]=dPhi;
1326       point[2]=ptj;
1327       point[3]=ptD;
1328       point[4]=deltamass;
1329       point[5]=static_cast<Double_t>(isSB);
1330       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1331       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1332       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1333    }
1334    if(fNAxesBigSparse==6){
1335       point=new Double_t[6];
1336       point[0]=z;
1337       point[1]=ptj;
1338       point[2]=ptD;
1339       point[3]=deltamass;
1340       point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1341       point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1342    }
1343    if(fNAxesBigSparse==7){
1344       point=new Double_t[7];
1345       point[0]=z;
1346       point[1]=ptj;
1347       point[2]=ptD;
1348       point[3]=deltamass;
1349       point[4]=0;
1350       point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1351       point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1352    }
1353
1354    //if(fIsDInJet) hPtJetWithD->Fill(ptj,deltamass,ptD);
1355
1356    FillMassHistograms(deltamass, ptD);
1357    if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1358
1359 }
1360
1361 //_______________________________________________________________________________
1362
1363 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t z,Double_t ptD,Double_t ptjet, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc){
1364
1365    Double_t pdgmass=0;
1366    if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1367    if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
1368    //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1369
1370    //Double_t point[9]={z,dPhi,ptjet,ptD,pdgmass,0, static_cast<Double_t>(fIsDInJet ? 1 : 0),bDInEMCalAcc,bJetInEMCalAcc};
1371
1372    Double_t *point=0x0;
1373    if(fNAxesBigSparse==9){
1374       point=new Double_t[9];
1375       point[0]=z;
1376       point[1]=dPhi;
1377       point[2]=ptjet;
1378       point[3]=ptD;
1379       point[4]=pdgmass;
1380       point[5]=0;
1381       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1382       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1383       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1384    }
1385    if(fNAxesBigSparse==6){
1386       point=new Double_t[6];
1387       point[0]=z;
1388       point[1]=ptjet;
1389       point[2]=ptD;
1390       point[3]=pdgmass;
1391       point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1392       point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1393    }
1394       if(fNAxesBigSparse==7){
1395       point=new Double_t[7];
1396       point[0]=z;
1397       point[1]=ptjet;
1398       point[2]=ptD;
1399       point[3]=pdgmass;
1400       point[4]=1;
1401       point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1402       point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1403    }
1404
1405
1406
1407    if(fNAxesBigSparse==9) point[4]=pdgmass;
1408    if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=pdgmass;
1409    if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1410    //if(fIsDInJet) {
1411    //  hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet
1412    //}
1413
1414 }
1415
1416 //_______________________________________________________________________________
1417
1418 void AliAnalysisTaskFlavourJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD){
1419
1420    if(fIsDInJet) {
1421       fhInvMassptD->Fill(mass,ptD);
1422    }
1423 }
1424
1425 //________________________________________________________________________________
1426
1427 void AliAnalysisTaskFlavourJetCorrelations::FlagFlavour(AliEmcalJet *jet){
1428
1429    AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
1430    if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;
1431    if (fIsDInJet) jet->AddFlavourTag(tag);
1432
1433    return;
1434
1435 }
1436
1437 //_______________________________________________________________________________
1438
1439 void AliAnalysisTaskFlavourJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){
1440
1441    //  D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
1442    // (expected detector resolution) on the left and right frm the D0 mass. Each band
1443    //  has a width of ~5 sigmas. Two band needed  for opening angle considerations
1444
1445    //TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
1446
1447    Bool_t bDInEMCalAcc=InEMCalAcceptance(candDstar);
1448    Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1449
1450    Double_t deltaM=candDstar->DeltaInvMass();
1451    //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);
1452    Double_t z=Z(candDstar,jet);
1453    Double_t ptD=candDstar->Pt();
1454
1455    Double_t dPhi=jet->Phi()-candDstar->Phi();
1456
1457    if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1458    if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1459
1460    Int_t isSideBand=1;//IsDzeroSideBand(candDstar);
1461    //if no SB analysis we should not enter here, so no need of checking the number of axes
1462    Double_t point[9]={z,dPhi,fPtJet,ptD,deltaM,static_cast<Double_t>(isSideBand), static_cast<Double_t>(fIsDInJet ? 1 : 0),static_cast<Double_t>(bDInEMCalAcc? 1 : 0),static_cast<Double_t>(bJetInEMCalAcc? 1 : 0)};
1463    //fill the background histogram with the side bands or when looking at MC with the real background
1464    if(isSideBand==1){
1465       fhDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background
1466       //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);
1467       if(fSwitchOnSparses) fhsDphiz->Fill(point,1.);
1468       //if(fIsDInJet){  // evaluate in the near side
1469       //         hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
1470       //}
1471    }
1472 }
1473
1474 //_______________________________________________________________________________
1475
1476 void AliAnalysisTaskFlavourJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg,AliEmcalJet* jet){
1477
1478    //need mass, deltaR, pt jet, ptD
1479    //two cases: D0 and Dstar
1480
1481    Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());
1482    if(!isselected) return;
1483
1484    Double_t ptD=candbg->Pt();
1485    Double_t deltaR=DeltaR(jet,candbg);
1486    Double_t phiD=candbg->Phi();
1487    Double_t deltaphi = jet->Phi()-phiD;
1488    if(deltaphi<=-(TMath::Pi())/2.) deltaphi = deltaphi+2.*(TMath::Pi());
1489    if(deltaphi>(3.*(TMath::Pi()))/2.) deltaphi = deltaphi-2.*(TMath::Pi());
1490    Double_t z=Z(candbg,jet);
1491
1492    if(fIsDInJet) fhzDT->Fill(Z(candbg,jet,kTRUE));
1493
1494
1495
1496
1497    fhDeltaRD->Fill(deltaR);
1498    fhDeltaRptDptjB->Fill(deltaR,ptD,fPtJet);
1499
1500    Bool_t bDInEMCalAcc=InEMCalAcceptance(candbg);
1501    Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1502
1503    //TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
1504
1505    Double_t *point=0x0;
1506    if(fNAxesBigSparse==9){
1507       point=new Double_t[9];
1508       point[0]=z;
1509       point[1]=deltaphi;
1510       point[2]=fPtJet;
1511       point[3]=ptD;
1512       point[4]=-999; //set below
1513       point[5]=1;
1514       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1515       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1516       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1517    }
1518
1519    if(fNAxesBigSparse==7){
1520       point=new Double_t[7];
1521       point[0]=z;
1522       point[1]=fPtJet;
1523       point[2]=ptD;
1524       point[3]=-999; //set below
1525       point[4]=1;
1526       point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1527       point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1528    }
1529
1530    if(fCandidateType==kDstartoKpipi){
1531       AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;
1532       Double_t deltaM=dstarbg->DeltaInvMass();
1533       fhInvMassptDbg->Fill(deltaM,ptD);
1534       //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
1535       if(fNAxesBigSparse==9) point[4]=deltaM;
1536       if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=deltaM;
1537       if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1538    }
1539
1540    if(fCandidateType==kD0toKpi){
1541       Double_t masses[2]={0.,0.};
1542       Int_t pdgdaughtersD0[2]={211,321};//pi,K
1543       Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1544
1545       masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1546       masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1547
1548
1549       if(isselected==1 || isselected==3) {
1550          //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[0],ptD);
1551          fhInvMassptDbg->Fill(masses[0],ptD);
1552          if(fNAxesBigSparse==9) point[4]=masses[0];
1553          if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[0];
1554          if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1555      }
1556       if(isselected>=2) {
1557          //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[1],ptD);
1558          fhInvMassptDbg->Fill(masses[1],ptD);
1559          if(fNAxesBigSparse==9) point[4]=masses[1];
1560          if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[1];
1561          if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1562
1563       }
1564
1565
1566    }
1567 }
1568
1569 //_______________________________________________________________________________
1570
1571 Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliEmcalJet *p1, AliVParticle *p2) const {
1572    //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1573    //It recalculates the eta-phi values if it was asked for background subtraction of the jets
1574    if(!p1 || !p2) return -1;
1575
1576     Double_t phi1=p1->Phi(),eta1=p1->Eta();
1577
1578     //It subtracts the backgroud of jets if it was asked for it.
1579     TString sname("Rho");
1580     if (!sname.IsNull()) {
1581         AliRhoParameter *rho = 0;
1582         rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(sname));
1583         if(rho){
1584             Double_t pj[3];
1585             Bool_t okpj=p1->PxPyPz(pj);
1586             if(!okpj){
1587                 printf("Problems getting momenta\n");
1588                 return -999;
1589             }
1590             Double_t rhoval = rho->GetVal();
1591             pj[0] = p1->Px() - p1->Area()*(rhoval*TMath::Cos(p1->AreaPhi()));
1592             pj[1] = p1->Py() - p1->Area()*(rhoval*TMath::Sin(p1->AreaPhi()));
1593             pj[2] = p1->Pz() - p1->Area()*(rhoval*TMath::SinH(p1->AreaEta()));
1594             //Image of the function Arccos(px/pt) where pt = sqrt(px*px+py*py) is:
1595             //[0,pi]    if py > 0 and
1596             //[pi,2pi]  if py < 0
1597             phi1 = TMath::ACos(pj[0]/TMath::Sqrt(pj[0]*pj[0]+pj[1]*pj[1]));
1598             if(pj[1]<0) phi1 = 2*TMath::Pi() - phi1;
1599             eta1 = TMath::ASinH(pj[2]/TMath::Sqrt(pj[0]*pj[0]+pj[1]*pj[1]));
1600         }
1601     }
1602
1603    Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
1604
1605    Double_t dPhi=phi1-phi2;
1606    if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1607    if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1608
1609    Double_t dEta=eta1-eta2;
1610    Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
1611    return deltaR;
1612
1613 }
1614
1615 //_______________________________________________________________________________
1616
1617 Int_t AliAnalysisTaskFlavourJetCorrelations::IsDzeroSideBand(AliAODRecoCascadeHF *candDstar){
1618
1619    Double_t ptD=candDstar->Pt();
1620    Int_t bin = fCuts->PtBin(ptD);
1621    if (bin < 0){
1622       // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
1623       bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?
1624       return -1;  // out of bounds
1625    }
1626
1627    Double_t invM=candDstar->InvMassD0();
1628    Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1629
1630    Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
1631    Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
1632
1633    if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
1634    else return 0;
1635
1636 }
1637
1638 //_______________________________________________________________________________
1639
1640 Bool_t AliAnalysisTaskFlavourJetCorrelations::AreDaughtersInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Int_t*& daughOutOfJetID, AliAODTrack**& dtrks, Bool_t fillH){
1641    //returns true/false and the array of particles not found among jet constituents
1642
1643
1644    Int_t daughtersID[3];
1645    Int_t ndaugh=0;
1646    Int_t dmatchedID[3]={0,0,0};
1647    Int_t countmatches=0;
1648    if (fCandidateType==kDstartoKpipi){
1649       AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)charm;
1650       AliAODRecoDecayHF2Prong* dzero=dstar->Get2Prong();
1651       daughtersID[0]=dzero->GetProngID(0);
1652       daughtersID[1]=dzero->GetProngID(1);
1653       daughtersID[2]=dstar->GetBachelor()->GetID();
1654       ndaugh=3;
1655
1656       dtrks=new AliAODTrack*[3];
1657       dtrks[0]=(AliAODTrack*)dzero->GetDaughter(0);
1658       dtrks[1]=(AliAODTrack*)dzero->GetDaughter(1);
1659       dtrks[2]=(AliAODTrack*)dstar->GetBachelor();
1660
1661       //check
1662       if(fillH){
1663          if(daughtersID[0]!=((AliAODTrack*)dtrks[0])->GetID() || daughtersID[1]!=((AliAODTrack*)dtrks[1])->GetID())  fhControlDInJ->Fill(6);
1664
1665          fhIDddaugh->Fill(daughtersID[0]);
1666          fhIDddaugh->Fill(daughtersID[1]);
1667          fhIDddaugh->Fill(daughtersID[2]);
1668
1669       }
1670       //Printf("ID daughters %d, %d, %d",daughtersID[0], daughtersID[1], daughtersID[2]);
1671    }
1672
1673    if (fCandidateType==kD0toKpi){
1674       daughtersID[0]=charm->GetProngID(0);
1675       daughtersID[1]=charm->GetProngID(1);
1676       ndaugh=2;
1677       if(fillH){
1678          fhIDddaugh->Fill(daughtersID[0]);
1679          fhIDddaugh->Fill(daughtersID[1]);
1680       }
1681       dtrks=new AliAODTrack*[2];
1682       dtrks[0]=(AliAODTrack*)charm->GetDaughter(0);
1683       dtrks[1]=(AliAODTrack*)charm->GetDaughter(1);
1684
1685    }
1686
1687    const Int_t cndaugh=ndaugh;
1688    daughOutOfJetID=new Int_t[cndaugh];
1689
1690    Int_t nchtrk=thejet->GetNumberOfTracks();
1691    for(Int_t itk=0;itk<nchtrk;itk++){
1692       AliVTrack* tkjet=((AliPicoTrack*)thejet->TrackAt(itk,fTrackArr))->GetTrack();
1693       Int_t idtkjet=tkjet->GetID();
1694       if(fillH) fhIDjetTracks->Fill(idtkjet);
1695       for(Int_t id=0;id<ndaugh;id++){
1696          if(idtkjet==daughtersID[id]) {
1697             dmatchedID[id]=daughtersID[id]; //daughter which matches a track in the jet
1698             countmatches++;
1699
1700          }
1701
1702          if(countmatches==ndaugh) break;
1703       }
1704    }
1705    //reverse: include in the array the ID of daughters not matching
1706
1707    for(Int_t id=0;id<ndaugh;id++){
1708       if(dmatchedID[id]==0) {
1709          daughOutOfJetID[id]=daughtersID[id];
1710          if(fillH) fhIDddaughOut->Fill(daughOutOfJetID[id]);
1711       }
1712       else daughOutOfJetID[id]=0;
1713    }
1714    if(fillH){
1715       if((ndaugh-countmatches) == 1) fhControlDInJ->Fill(0);
1716       if((ndaugh-countmatches) == 2) fhControlDInJ->Fill(1);
1717       if((ndaugh-countmatches) == 3) fhControlDInJ->Fill(2);
1718    }
1719    if(ndaugh!=countmatches){
1720       return kFALSE;
1721    }
1722
1723    return kTRUE;
1724 }
1725
1726 //_______________________________________________________________________________
1727
1728 Bool_t AliAnalysisTaskFlavourJetCorrelations::IsDInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Bool_t fillH){
1729
1730    //check the conditions type:
1731    //type 0 : DeltaR < jet radius, don't check daughters (and don't correct pt jet)
1732    //type 1 : DeltaR < jet radius and check for all daughters among jet tracks, don't correct ptjet
1733    //type 2 (default) : DeltaR < jet radius and check for all daughters among jet tracks, if not present, correct the ptjet
1734
1735    fPmissing[0]=0;
1736    fPmissing[1]=0;
1737    fPmissing[2]=0;
1738
1739    Bool_t testDeltaR=kFALSE;
1740    if(DeltaR(thejet,charm)<fJetRadius) testDeltaR=kTRUE;
1741
1742    Int_t* daughOutOfJet;
1743    AliAODTrack** charmDaugh;
1744    Bool_t testDaugh=AreDaughtersInJet(thejet, charm, daughOutOfJet,charmDaugh,fillH);
1745    if(testDaugh && testDeltaR) {
1746       //Z of candidates with daughters in the jet
1747       fhzDinjet->Fill(Z(charm,thejet));
1748
1749    }
1750    if(!testDaugh && testDeltaR && fTypeDInJet==2){
1751
1752       Int_t ndaugh=3;
1753       if(fCandidateType==kD0toKpi) ndaugh=2;
1754       Int_t nOut=ndaugh;
1755
1756       for(Int_t id=0;id<ndaugh;id++){
1757          if(daughOutOfJet[id]!=0){ //non-matched daughter
1758             //get track and its momentum
1759             nOut--;
1760             if(fillH) {
1761                fhControlDInJ->Fill(3);
1762                if(id<2) fhControlDInJ->Fill(4);
1763                if(id==2)fhControlDInJ->Fill(5);
1764                fhDRdaughOut->Fill(DeltaR(thejet, charmDaugh[id]));
1765             }
1766             fPmissing[0]+=charmDaugh[id]->Px();
1767             fPmissing[1]+=charmDaugh[id]->Py();
1768             fPmissing[2]+=charmDaugh[id]->Pz();
1769          }
1770
1771       }
1772
1773       //now the D in within the jet
1774       testDaugh=kTRUE;
1775
1776    }
1777
1778    delete[] charmDaugh;
1779
1780    Bool_t result=0;
1781    switch(fTypeDInJet){
1782    case 0:
1783       result=testDeltaR;
1784       break;
1785    case 1:
1786       result=testDeltaR && testDaugh;
1787       break;
1788    case 2:
1789       result=testDeltaR && testDaugh;
1790       break;
1791    default:
1792       AliInfo("Selection type not specified, use 1");
1793       result=testDeltaR && testDaugh;
1794       break;
1795    }
1796  return result;
1797 }
1798
1799 Bool_t AliAnalysisTaskFlavourJetCorrelations::InEMCalAcceptance(AliVParticle *vpart){
1800    //check eta phi of a VParticle: return true if it is in the EMCal acceptance, false otherwise
1801
1802    Double_t phiEMCal[2]={1.405,3.135},etaEMCal[2]={-0.7,0.7};
1803    Bool_t binEMCal=kTRUE;
1804    Double_t phi=vpart->Phi(), eta=vpart->Eta();
1805    if(phi<phiEMCal[0] || phi>phiEMCal[1]) binEMCal=kFALSE;
1806    if(eta<etaEMCal[0] || eta>etaEMCal[1]) binEMCal=kFALSE;
1807    return binEMCal;
1808
1809
1810 }