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
  34 #include "AliAnalysisTaskFlavourJetCorrelations.h"
  35 #include "AliAODMCHeader.h"
  36 #include "AliAODHandler.h"
  37 #include "AliAnalysisManager.h"
  38 #include "AliLog.h"
  39 #include "AliEmcalJet.h"
  40 #include "AliJetContainer.h"
  41 #include "AliAODRecoDecay.h"
  42 #include "AliAODRecoCascadeHF.h"
  43 #include "AliAODRecoDecayHF2Prong.h"
  44 #include "AliESDtrack.h"
  45 #include "AliAODMCParticle.h"
  46 #include "AliPicoTrack.h"
  47 #include "AliRDHFCutsD0toKpi.h"
  48 #include "AliRDHFCutsDStartoKpipi.h"
  49
  50 ClassImp(AliAnalysisTaskFlavourJetCorrelations)
  51
  52
  53 //_______________________________________________________________________________
  54
  55 AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations() :
  56 AliAnalysisTaskEmcalJet("",kFALSE),
  57 fUseMCInfo(kTRUE),
  58 fUseReco(kTRUE),
  59 fCandidateType(),
  60 fPDGmother(),
  61 fNProngs(),
  62 fPDGdaughters(),
  63 fBranchName(),
  64 fmyOutput(0),
  65 fCuts(0),
  66 fMinMass(),
  67 fMaxMass(),
  68 fJetArrName(0),
  69 fCandArrName(0),
  70 fLeadingJetOnly(kFALSE),
  71 fJetRadius(0),
  72 fCandidateArray(0),
  73 fSideBandArray(0)
  74 {
  75    //
  76    // Default ctor
  77    //
  78    //SetMakeGeneralHistograms(kTRUE);
  79
  80 }
  81
  82 //_______________________________________________________________________________
  83
  84 AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations(const Char_t* name, AliRDHFCuts* cuts,ECandidateType candtype) :
  85 AliAnalysisTaskEmcalJet(name,kFALSE),
  86 fUseMCInfo(kTRUE),
  87 fUseReco(kTRUE),
  88 fCandidateType(),
  89 fPDGmother(),
  90 fNProngs(),
  91 fPDGdaughters(),
  92 fBranchName(),
  93 fmyOutput(0),
  94 fCuts(0),
  95 fMinMass(),
  96 fMaxMass(),
  97 fJetArrName(0),
  98 fCandArrName(0),
  99 fLeadingJetOnly(kFALSE),
 100 fJetRadius(0),
 101 fCandidateArray(0),
 102 fSideBandArray(0)
 103 {
 104    //
 105    // Constructor. Initialization of Inputs and Outputs
 106    //
 107
 108    Info("AliAnalysisTaskFlavourJetCorrelations","Calling Constructor");
 109    fCuts=cuts;
 110    fCandidateType=candtype;
 111    const Int_t nptbins=fCuts->GetNPtBins();
 112    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};
 113
 114    switch(fCandidateType){
 115    case 0:
 116       fPDGmother=421;
 117       fNProngs=2;
 118       fPDGdaughters[0]=211;//pi
 119       fPDGdaughters[1]=321;//K
 120       fPDGdaughters[2]=0; //empty
 121       fPDGdaughters[3]=0; //empty
 122       fBranchName="D0toKpi";
 123       fCandArrName="D0";
 124       break;
 125    case 1:
 126       fPDGmother=413;
 127       fNProngs=3;
 128       fPDGdaughters[1]=211;//pi soft
 129       fPDGdaughters[0]=421;//D0
 130       fPDGdaughters[2]=211;//pi fromD0
 131       fPDGdaughters[3]=321; // K from D0
 132       fBranchName="Dstar";
 133       fCandArrName="Dstar";
 134
 135       if(nptbins<=13){
 136          for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]= defaultSigmaD013[ipt];
 137       } else {
 138          AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
 139       }
 140       break;
 141    default:
 142       printf("%d not accepted!!\n",fCandidateType);
 143       break;
 144    }
 145
 146    if(fCandidateType==kD0toKpi)SetMassLimits(0.15,fPDGmother);
 147    if(fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
 148
 149    DefineInput(1, TClonesArray::Class());
 150    DefineInput(2, TClonesArray::Class());
 151
 152    DefineOutput(1,TList::Class()); // histos
 153    DefineOutput(2,AliRDHFCuts::Class()); // my cuts
 154
 155 }
 156
 157 //_______________________________________________________________________________
 158
 159 AliAnalysisTaskFlavourJetCorrelations::~AliAnalysisTaskFlavourJetCorrelations() {
 160    //
 161    // destructor
 162    //
 163
 164    Info("~AliAnalysisTaskFlavourJetCorrelations","Calling Destructor");
 165
 166    delete fmyOutput;
 167    delete fCuts;
 168
 169 }
 170
 171 //_______________________________________________________________________________
 172
 173 void AliAnalysisTaskFlavourJetCorrelations::Init(){
 174    //
 175    // Initialization
 176    //
 177
 178    if(fDebug > 1) printf("AnalysisTaskRecoJetCorrelations::Init() \n");
 179    switch(fCandidateType){
 180    case 0:
 181       {
 182          AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
 183          copyfCuts->SetName("AnalysisCutsDzero");
 184          // Post the data
 185          PostData(2,copyfCuts);
 186       }
 187       break;
 188    case 1:
 189       {
 190          AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
 191          copyfCuts->SetName("AnalysisCutsDStar");
 192          // Post the cuts
 193          PostData(2,copyfCuts);
 194       }
 195       break;
 196    default:
 197       return;
 198    }
 199
 200    return;
 201 }
 202
 203 //_______________________________________________________________________________
 204
 205 void AliAnalysisTaskFlavourJetCorrelations::UserCreateOutputObjects() {
 206    // output
 207    Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 208    fmyOutput = new TList();
 209    fmyOutput->SetOwner();
 210    fmyOutput->SetName("pippo");
 211    // define histograms
 212    DefineHistoForAnalysis();
 213    PostData(1,fmyOutput);
 214
 215    return;
 216 }
 217
 218 //_______________________________________________________________________________
 219
 220 void AliAnalysisTaskFlavourJetCorrelations::UserExec(Option_t *)
 221 {
 222    // user exec
 223    if (!fInitialized){
 224       AliAnalysisTaskEmcalJet::ExecOnce();
 225    }
 226
 227    // Load the event
 228    AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 229
 230    TClonesArray *arrayDStartoD0pi=0;
 231    TClonesArray *trackArr = 0;
 232
 233    if (!aodEvent && AODEvent() && IsStandardAOD()) {
 234
 235       // In case there is an AOD handler writing a standard AOD, use the AOD
 236       // event in memory rather than the input (ESD) event.
 237       aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
 238
 239       // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 240       // have to taken from the AOD event hold by the AliAODExtension
 241       AliAODHandler* aodHandler = (AliAODHandler*)
 242       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 243       if(aodHandler->GetExtensions()) {
 244          AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 245          AliAODEvent *aodFromExt = ext->GetAOD();
 246          arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
 247       }
 248    } else {
 249       arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
 250    }
 251
 252    if (!arrayDStartoD0pi) {
 253       AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
 254       //  return;
 255    } else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
 256
 257    TClonesArray* mcArray = 0x0;
 258    if (fUseMCInfo) {
 259       mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 260       if (!mcArray) {
 261          printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
 262          return;
 263       }
 264    }
 265
 266    //retrieve jets
 267    trackArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("PicoTracks"));
 268    //clusArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClustersCorr"));
 269    //jetArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrName));
 270    fJetRadius=GetJetContainer(0)->GetJetRadius();
 271
 272    if(!trackArr){
 273       AliInfo(Form("Could not find the track array\n"));
 274       return;
 275    }
 276
 277
 278    fCandidateArray = dynamic_cast<TClonesArray*>(GetInputData(1));
 279    if (!fCandidateArray) return;
 280    if (fCandidateType==1) {
 281       fSideBandArray = dynamic_cast<TClonesArray*>(GetInputData(2));
 282       if (!fSideBandArray) return;
 283    }
 284    //Printf("ncandidates found %d",fCandidateArray->GetEntriesFast());
 285
 286    //Histograms
 287    TH1I* hstat=(TH1I*)fmyOutput->FindObject("hstat");
 288    TH1F* hPtJetTrks=(TH1F*)fmyOutput->FindObject("hPtJetTrks");
 289    TH1F* hPhiJetTrks=(TH1F*)fmyOutput->FindObject("hPhiJetTrks");
 290    TH1F* hEtaJetTrks=(TH1F*)fmyOutput->FindObject("hEtaJetTrks");
 291    TH1F* hEjetTrks=(TH1F*)fmyOutput->FindObject("hEjetTrks");
 292    TH1F* hPtJet=(TH1F*)fmyOutput->FindObject("hPtJet");
 293    TH1F* hPhiJet=(TH1F*)fmyOutput->FindObject("hPhiJet");
 294    TH1F* hEtaJet=(TH1F*)fmyOutput->FindObject("hEtaJet");
 295    TH1F* hEjet=(TH1F*)fmyOutput->FindObject("hEjet");
 296    TH1F* hNtrArr=(TH1F*)fmyOutput->FindObject("hNtrArr");
 297    TH1F* hNJetPerEv=(TH1F*)fmyOutput->FindObject("hNJetPerEv");
 298    TH1F* hdeltaRJetTracks=(TH1F*)fmyOutput->FindObject("hdeltaRJetTracks");
 299    TH1F* hNDPerEvNoJet=(TH1F*)fmyOutput->FindObject("hNDPerEvNoJet");
 300    TH1F* hptDPerEvNoJet=(TH1F*)fmyOutput->FindObject("hptDPerEvNoJet");
 301    TH1F* hNJetPerEvNoD=(TH1F*)fmyOutput->FindObject("hNJetPerEvNoD");
 302    TH1F* hPtJetPerEvNoD=(TH1F*)fmyOutput->FindObject("hPtJetPerEvNoD");
 303    THnSparseF* hnspDstandalone=(THnSparseF*)fmyOutput->FindObject("hsDstandalone");
 304
 305    hstat->Fill(0);
 306
 307    // fix for temporary bug in ESDfilter
 308    // the AODs with null vertex pointer didn't pass the PhysSel
 309    if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
 310
 311    //Event selection
 312    Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
 313    TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
 314    if(!iseventselected) return;
 315
 316    hstat->Fill(1);
 317
 318    //retrieve charm candidates selected
 319    Int_t candidates = fCandidateArray->GetEntriesFast();
 320
 321    //trigger on jets
 322
 323    Int_t njets=GetJetContainer()->GetNJets();
 324    if(njets == 0) {
 325       hstat->Fill(6, candidates);
 326       hNDPerEvNoJet->Fill(candidates);
 327       for(Int_t iD=0;iD<candidates;iD++){
 328          AliVParticle* cand=(AliVParticle*)fCandidateArray->At(iD);
 329          if(!cand) continue;
 330          hptDPerEvNoJet->Fill(cand->Pt());
 331
 332       }
 333       return;
 334
 335    }
 336
 337    //loop on candidates standalone (checking the candidates are there and their phi-eta distributions)
 338
 339    for(Int_t ic = 0; ic < candidates; ic++) {
 340
 341       // D* candidates
 342       AliAODRecoDecayHF* charm=0x0;
 343       AliAODRecoCascadeHF* dstar=0x0;
 344
 345
 346       charm=(AliAODRecoDecayHF*)fCandidateArray->At(ic);
 347       if(!charm) continue;
 348
 349       if(fCandidateType==kDstartoKpipi){
 350          dstar=(AliAODRecoCascadeHF*)fCandidateArray->At(ic);
 351          if(!dstar) continue;
 352       }
 353
 354       hstat->Fill(2);
 355
 356       Double_t candsparse[4]={charm->Eta(), charm->Phi(), charm->Pt(), 0};
 357
 358       if(fCandidateType==kDstartoKpipi) {
 359          Double_t deltamass= dstar->DeltaInvMass();
 360          candsparse[3]=deltamass;
 361          hnspDstandalone->Fill(candsparse);
 362       }
 363       if(fCandidateType==kD0toKpi){
 364          AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)charm;
 365          Int_t isselected=fCuts->IsSelected(dzero,AliRDHFCuts::kAll,aodEvent);
 366
 367          Double_t masses[2];
 368          Int_t pdgdaughtersD0[2]={211,321};//pi,K
 369          Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
 370
 371          masses[0]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
 372          masses[1]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
 373          if(isselected==1 || isselected==3) {
 374             candsparse[3]=masses[0];
 375             hnspDstandalone->Fill(candsparse);
 376          }
 377          if(isselected>=2){
 378             candsparse[3]=masses[1];
 379             hnspDstandalone->Fill(candsparse);
 380
 381          }
 382       }
 383    }
 384
 385    // we start with jets
 386    Double_t ejet   = 0;
 387    Double_t phiJet = 0;
 388    Double_t etaJet = 0;
 389    Double_t ptjet = 0;
 390    Double_t leadingJet =0;
 391
 392    Int_t ntrarr=trackArr->GetEntriesFast();
 393    hNtrArr->Fill(ntrarr);
 394
 395    for(Int_t i=0;i<ntrarr;i++){
 396       AliVTrack *jtrack=static_cast<AliVTrack*>(trackArr->At(i));
 397       //reusing histograms
 398       hPtJetTrks->Fill(jtrack->Pt());
 399       hPhiJetTrks->Fill(jtrack->Phi());
 400       hEtaJetTrks->Fill(jtrack->Eta());
 401       hEjetTrks->Fill(jtrack->E());
 402    }
 403
 404
 405    //option to use only the leading jet
 406    if(fLeadingJetOnly){
 407       for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) {
 408          AliEmcalJet* jetL = (AliEmcalJet*)GetJetFromArray(iJetsL);
 409          ptjet   = jetL->Pt();
 410          if(ptjet>leadingJet ) leadingJet = ptjet;
 411
 412       }
 413    }
 414
 415    Int_t cntjet=0;
 416    //loop over jets and consider only the leading jet in the event
 417    for (Int_t iJets = 0; iJets<njets; iJets++) {
 418       //Printf("Jet N %d",iJets);
 419       AliEmcalJet* jet = (AliEmcalJet*)GetJetFromArray(iJets);
 420       //Printf("Corr task Accept Jet");
 421       if(!AcceptJet(jet)) {
 422          hstat->Fill(5);
 423          continue;
 424       }
 425
 426       //jets variables
 427       ejet   = jet->E();
 428       phiJet = jet->Phi();
 429       etaJet = jet->Eta();
 430       ptjet = jet->Pt();
 431
 432       // choose the leading jet
 433       if(fLeadingJetOnly && (ejet<leadingJet)) continue;
 434       //used for normalization
 435       hstat->Fill(3);
 436       cntjet++;
 437       // fill energy, eta and phi of the jet
 438       hEjet   ->Fill(ejet);
 439       hPhiJet ->Fill(phiJet);
 440       hEtaJet ->Fill(etaJet);
 441       hPtJet  ->Fill(ptjet);
 442
 443       //loop on jet particles
 444       Int_t ntrjet=  jet->GetNumberOfTracks();
 445       for(Int_t itrk=0;itrk<ntrjet;itrk++){
 446
 447          AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,trackArr);
 448          hdeltaRJetTracks->Fill(DeltaR(jet,jetTrk));
 449
 450       }//end loop on jet tracks
 451
 452       if(candidates==0){
 453          hstat->Fill(7);
 454          hPtJetPerEvNoD->Fill(jet->Pt());
 455       }
 456
 457       //Printf("N candidates %d ", candidates);
 458       for(Int_t ic = 0; ic < candidates; ic++) {
 459
 460          // D* candidates
 461          AliVParticle* charm=0x0;
 462          charm=(AliVParticle*)fCandidateArray->At(ic);
 463          if(!charm) continue;
 464
 465
 466          FlagFlavour(charm, jet);
 467          if (jet->TestFlavourTag(AliEmcalJet::kDStar)) hstat->Fill(4);
 468
 469          FillHistogramsRecoJetCorr(charm, jet);
 470
 471       }
 472       //retrieve side band background candidates for Dstar
 473       Int_t nsbcand = 0; if(fSideBandArray) nsbcand=fSideBandArray->GetEntriesFast();
 474
 475       for(Int_t ib=0;ib<nsbcand;ib++){
 476          if(fCandidateType==kDstartoKpipi && !fUseMCInfo){
 477             AliAODRecoCascadeHF *sbcand=(AliAODRecoCascadeHF*)fSideBandArray->At(ib);
 478             if(!sbcand) continue;
 479             SideBandBackground(sbcand,jet);
 480          }
 481          if(fUseMCInfo){
 482             AliAODRecoDecayHF* charmbg = 0x0;
 483             charmbg=(AliAODRecoDecayHF*)fCandidateArray->At(ib);
 484             if(!charmbg) continue;
 485             MCBackground(charmbg,jet);
 486          }
 487       }
 488    } // end of jet loop
 489
 490    hNJetPerEv->Fill(cntjet);
 491    if(candidates==0) hNJetPerEvNoD->Fill(cntjet);
 492    PostData(1,fmyOutput);
 493
 494 }
 495
 496 //_______________________________________________________________________________
 497
 498 void AliAnalysisTaskFlavourJetCorrelations::Terminate(Option_t*)
 499 {
 500    // The Terminate() function is the last function to be called during
 501    // a query. It always runs on the client, it can be used to present
 502    // the results graphically or save the results to file.
 503
 504    Info("Terminate"," terminate");
 505    AliAnalysisTaskSE::Terminate();
 506
 507    fmyOutput = dynamic_cast<TList*> (GetOutputData(1));
 508    if (!fmyOutput) {
 509       printf("ERROR: fmyOutput not available\n");
 510       return;
 511    }
 512 }
 513
 514 //_______________________________________________________________________________
 515
 516 void  AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){
 517    Float_t mass=0;
 518    Int_t abspdg=TMath::Abs(pdg);
 519
 520    mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
 521    // compute the Delta mass for the D*
 522    if(fCandidateType==kDstartoKpipi){
 523       Float_t mass1=0;
 524       mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 525       mass = mass-mass1;
 526    }
 527
 528    fMinMass = mass-range;
 529    fMaxMass = mass+range;
 530
 531    AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
 532    if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
 533 }
 534
 535 //_______________________________________________________________________________
 536
 537 void  AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){
 538    if(uplimit>lowlimit)
 539    {
 540       fMinMass = lowlimit;
 541       fMaxMass = uplimit;
 542    }
 543    else{
 544       printf("Error! Lower limit larger than upper limit!\n");
 545       fMinMass = uplimit - uplimit*0.2;
 546       fMaxMass = uplimit;
 547    }
 548 }
 549
 550 //_______________________________________________________________________________
 551
 552 Bool_t AliAnalysisTaskFlavourJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){
 553    if(nptbins>30) {
 554       AliInfo("Maximum number of bins allowed is 30!");
 555       return kFALSE;
 556    }
 557    if(!width) return kFALSE;
 558    for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
 559    return kTRUE;
 560 }
 561
 562 //_______________________________________________________________________________
 563
 564 Double_t AliAnalysisTaskFlavourJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet) const{
 565    if(!part || !jet){
 566       printf("Particle or jet do not exist!\n");
 567       return -99;
 568    }
 569    Double_t p[3],pj[3];
 570    Bool_t okpp=part->PxPyPz(p);
 571    Bool_t okpj=jet->PxPyPz(pj);
 572    if(!okpp || !okpj){
 573       printf("Problems getting momenta\n");
 574       return -999;
 575    }
 576    Double_t pjet=jet->P();
 577    Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet*pjet);
 578    return z;
 579 }
 580
 581 //_______________________________________________________________________________
 582
 583 Bool_t  AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){
 584
 585    // Statistics
 586    TH1I* hstat=new TH1I("hstat","Statistics",8,-0.5,7.5);
 587    hstat->GetXaxis()->SetBinLabel(1,"N ev anal");
 588    hstat->GetXaxis()->SetBinLabel(2,"N ev sel");
 589    hstat->GetXaxis()->SetBinLabel(3,"N cand sel & jet");
 590    hstat->GetXaxis()->SetBinLabel(4,"N jets");
 591    hstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
 592    hstat->GetXaxis()->SetBinLabel(6,"N jet rej");
 593    hstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet");
 594    hstat->GetXaxis()->SetBinLabel(8,"N jets & !D");
 595    hstat->SetNdivisions(1);
 596    fmyOutput->Add(hstat);
 597
 598    const Int_t nbinsmass=200;
 599    const Int_t nbinsptjet=500;
 600    const Int_t nbinsptD=100;
 601    const Int_t nbinsz=100;
 602    const Int_t nbinsphi=200;
 603    const Int_t nbinseta=100;
 604
 605    const Float_t ptjetlims[2]={0.,200.};
 606    const Float_t ptDlims[2]={0.,50.};
 607    const Float_t zlims[2]={0.,1.2};
 608    const Float_t philims[2]={0.,6.3};
 609    const Float_t etalims[2]={-1.5,1.5};
 610
 611    if(fCandidateType==kDstartoKpipi)
 612    {
 613
 614       TH2F* hDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
 615       hDiffSideBand->SetStats(kTRUE);
 616       hDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
 617       hDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
 618       hDiffSideBand->Sumw2();
 619       fmyOutput->Add(hDiffSideBand);
 620
 621
 622       TH1F* hPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
 623       hPtPion->SetStats(kTRUE);
 624       hPtPion->GetXaxis()->SetTitle("GeV/c");
 625       hPtPion->GetYaxis()->SetTitle("Entries");
 626       hPtPion->Sumw2();
 627       fmyOutput->Add(hPtPion);
 628
 629    }
 630    // D related histograms
 631    TH1F *hNDPerEvNoJet=new TH1F("hNDPerEvNoJet","Number of candidates per event with no jets; N candidate/ev with no jet", 20, 0., 20.);
 632    hNDPerEvNoJet->Sumw2();
 633    fmyOutput->Add(hNDPerEvNoJet);
 634
 635    TH1F *hptDPerEvNoJet=new TH1F("hptDPerEvNoJet","pt distribution of candidates per events with no jets; p_{t}^{D} (GeV/c)",nbinsptD, ptDlims[0],ptDlims[1]);
 636    hptDPerEvNoJet->Sumw2();
 637    fmyOutput->Add(hptDPerEvNoJet);
 638
 639    const Int_t    nAxisD=4;
 640    const Int_t    nbinsSparseD[nAxisD]={nbinseta,nbinsphi,nbinsptD,nbinsmass};
 641    const Double_t minSparseD[nAxisD]  ={etalims[0],philims[0],ptDlims[0],fMinMass};
 642    const Double_t maxSparseD[nAxisD]  ={etalims[1],philims[1],ptDlims[1],fMaxMass};
 643    THnSparseF *hsDstandalone=new THnSparseF("hsDstandalone","#phi, #eta, p_{T}^{D}, and mass", nAxisD, nbinsSparseD, minSparseD, maxSparseD);
 644    hsDstandalone->Sumw2();
 645
 646    fmyOutput->Add(hsDstandalone);
 647
 648    // jet related fistograms
 649
 650    TH1F* hEjetTrks      = new TH1F("hEjetTrks",  "Jet tracks energy distribution;Energy (GeV)",500,0,200);
 651    hEjetTrks->Sumw2();
 652    TH1F* hPhiJetTrks    = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)",  nbinsphi,philims[0],philims[1]);
 653    hPhiJetTrks->Sumw2();
 654    TH1F* hEtaJetTrks    = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta",  nbinseta,etalims[0],etalims[1]);
 655    hEtaJetTrks->Sumw2();
 656    TH1F* hPtJetTrks     = new TH1F("hPtJetTrks",  "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
 657    hPtJetTrks->Sumw2();
 658
 659    TH1F* hEjet      = new TH1F("hEjet",  "Jet energy distribution;Energy (GeV)",500,0,200);
 660    hEjet->Sumw2();
 661    TH1F* hPhiJet    = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)",  nbinsphi,philims[0],philims[1]);
 662    hPhiJet->Sumw2();
 663    TH1F* hEtaJet    = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
 664    hEtaJet->Sumw2();
 665    TH1F* hPtJet      = new TH1F("hPtJet",  "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
 666    hPtJet->Sumw2();
 667
 668    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]);
 669    hPtJetWithD->Sumw2();
 670    //for the MC this histogram is filled with the real background
 671    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]);
 672    hPtJetWithDsb->Sumw2();
 673
 674    TH1F* hdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.);
 675    hdeltaRJetTracks->Sumw2();
 676
 677    TH1F* hNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000);
 678    hNtrArr->Sumw2();
 679
 680    TH1F *hNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);
 681    hNJetPerEv->Sumw2();
 682
 683    TH1F *hNJetPerEvNoD=new TH1F("hNJetPerEvNoD","Number of jets per event with no D; number of jets/ev with no D",10,-0.5,9.5);
 684    hNJetPerEvNoD->Sumw2();
 685
 686    TH1F *hPtJetPerEvNoD=new TH1F("hPtJetPerEvNoD","pt distribution of jets per event with no D; p_{T}^{jet} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
 687    hPtJetPerEvNoD->Sumw2();
 688
 689    fmyOutput->Add(hEjetTrks);
 690    fmyOutput->Add(hPhiJetTrks);
 691    fmyOutput->Add(hEtaJetTrks);
 692    fmyOutput->Add(hPtJetTrks);
 693    fmyOutput->Add(hEjet);
 694    fmyOutput->Add(hPhiJet);
 695    fmyOutput->Add(hEtaJet);
 696    fmyOutput->Add(hPtJet);
 697    fmyOutput->Add(hPtJetWithD);
 698    fmyOutput->Add(hPtJetWithDsb);
 699    fmyOutput->Add(hdeltaRJetTracks);
 700    fmyOutput->Add(hNtrArr);
 701    fmyOutput->Add(hNJetPerEv);
 702    fmyOutput->Add(hNJetPerEvNoD);
 703    fmyOutput->Add(hPtJetPerEvNoD);
 704
 705    TH1F* hDeltaRD=new TH1F("hDeltaRD","#Delta R distribution of D candidates selected;#Delta R",200, 0.,10.);
 706    hDeltaRD->Sumw2();
 707    fmyOutput->Add(hDeltaRD);
 708
 709    //background (side bands for the Dstar and like sign for D0)
 710    fJetRadius=GetJetContainer(0)->GetJetRadius();
 711    TH2F* hInvMassptD = new TH2F("hInvMassptD",Form("D (Delta R < %.1f) invariant mass distribution p_{T}^{j} > threshold",fJetRadius),nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);
 712    hInvMassptD->SetStats(kTRUE);
 713    hInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
 714    hInvMassptD->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
 715    hInvMassptD->Sumw2();
 716
 717    fmyOutput->Add(hInvMassptD);
 718
 719    if(fUseMCInfo){
 720       TH2F* hInvMassptDbg = 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]);
 721       hInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
 722       hInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
 723       hInvMassptDbg->Sumw2();
 724       fmyOutput->Add(hInvMassptDbg);
 725
 726    }
 727    Double_t pi=TMath::Pi(), philims2[2];
 728    philims2[0]=-pi*1./2.;
 729    philims2[1]=pi*3./2.;
 730    const Int_t nAxis=7;
 731    const Int_t nbinsSparse[nAxis]={nbinsz,nbinsphi,nbinsptjet,nbinsptD,nbinsmass,2, 2};
 732    const Double_t minSparse[nAxis]={zlims[0],philims2[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5, -0.5};
 733    const Double_t maxSparse[nAxis]={zlims[1],philims2[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5};
 734    THnSparseF *hsDphiz=new THnSparseF("hsDphiz","Z and #Delta#phi vs p_{T}^{jet}, p_{T}^{D}, mass. For SB or not and within the jet cone or not", nAxis, nbinsSparse, minSparse, maxSparse);
 735    hsDphiz->Sumw2();
 736
 737    fmyOutput->Add(hsDphiz);
 738
 739    PostData(1, fmyOutput);
 740
 741    return kTRUE;
 742 }
 743
 744 //_______________________________________________________________________________
 745
 746 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet){
 747
 748    Double_t ptD=candidate->Pt();
 749    Double_t ptjet=jet->Pt();
 750    Double_t deltaR=DeltaR(candidate,jet);
 751    Double_t deltaphi = jet->Phi()-candidate->Phi();
 752    if(deltaphi<=-(TMath::Pi())/2) deltaphi = deltaphi+2*(TMath::Pi());
 753    if(deltaphi>(3*(TMath::Pi()))/2) deltaphi = deltaphi-2*(TMath::Pi());
 754    Double_t z=Z(candidate,jet);
 755
 756    TH1F* hDeltaRD=(TH1F*)fmyOutput->FindObject("hDeltaRD");
 757    hDeltaRD->Fill(deltaR);
 758    if(fUseReco){
 759       if(fCandidateType==kD0toKpi) {
 760          AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;
 761          FillHistogramsD0JetCorr(dzero,deltaphi,z,ptD,ptjet,deltaR, AODEvent());
 762
 763       }
 764
 765       if(fCandidateType==kDstartoKpipi) {
 766          AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;
 767          FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,ptjet,deltaR);
 768
 769       }
 770    } else{ //generated level
 771       //AliInfo("Non reco");
 772       FillHistogramsMCGenDJetCorr(deltaphi,z,ptD,ptjet,deltaR);
 773    }
 774 }
 775
 776 //_______________________________________________________________________________
 777
 778 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj,Double_t deltaR, AliAODEvent* aodEvent){
 779
 780
 781    Double_t masses[2]={0.,0.};
 782    Int_t pdgdaughtersD0[2]={211,321};//pi,K
 783    Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
 784
 785    masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
 786    masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
 787
 788    TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");
 789    THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");
 790    Double_t point[7]={z,dPhi,ptj,ptD,masses[0],0, deltaR<fJetRadius ? 1 : 0};
 791
 792    Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
 793    if(isselected==1 || isselected==3) {
 794
 795       if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[0],ptD);
 796
 797       FillMassHistograms(masses[0], ptD, deltaR);
 798       hsDphiz->Fill(point,1.);
 799    }
 800    if(isselected>=2) {
 801       if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[1],ptD);
 802
 803       FillMassHistograms(masses[1], ptD, deltaR);
 804       point[4]=masses[1];
 805       hsDphiz->Fill(point,1.);
 806    }
 807
 808 }
 809
 810 //_______________________________________________________________________________
 811
 812 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj, Double_t deltaR){
 813   //dPhi and z not used at the moment,but will be (re)added
 814
 815    AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
 816    Double_t deltamass= dstar->DeltaInvMass();
 817    //Double_t massD0= dstar->InvMassD0();
 818
 819
 820    TH1F* hPtPion=(TH1F*)fmyOutput->FindObject("hPtPion");
 821    hPtPion->Fill(softpi->Pt());
 822
 823    TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");
 824    THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");
 825    Int_t isSB=0;//IsDzeroSideBand(dstar);
 826
 827    Double_t point[7]={z,dPhi,ptj,ptD,deltamass,isSB, deltaR<fJetRadius ? 1 : 0};
 828
 829    if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,deltamass,ptD);
 830
 831    FillMassHistograms(deltamass, ptD, deltaR);
 832    hsDphiz->Fill(point,1.);
 833
 834 }
 835
 836 //_______________________________________________________________________________
 837
 838 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t z,Double_t ptD,Double_t ptjet,Double_t deltaR){
 839
 840    Double_t pdgmass=0;
 841    TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");
 842    THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");
 843    Double_t point[7]={z,dPhi,ptjet,ptD,pdgmass,0, deltaR<fJetRadius ? 1 : 0};
 844
 845    if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 846    if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
 847    point[4]=pdgmass;
 848
 849    if(deltaR<fJetRadius) {
 850      hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet
 851      hsDphiz->Fill(point,1.);
 852    }
 853
 854 }
 855
 856 //_______________________________________________________________________________
 857
 858 void AliAnalysisTaskFlavourJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD, Double_t deltaR){
 859
 860    if(deltaR<fJetRadius) {
 861       TH2F* hInvMassptD=(TH2F*)fmyOutput->FindObject("hInvMassptD");
 862       hInvMassptD->Fill(mass,ptD);
 863    }
 864 }
 865
 866 //________________________________________________________________________________
 867
 868 void AliAnalysisTaskFlavourJetCorrelations::FlagFlavour(AliVParticle *charm, AliEmcalJet *jet){
 869    Double_t deltaR=DeltaR(charm, jet);
 870    AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
 871    if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;
 872    if (deltaR<fJetRadius) jet->AddFlavourTag(tag);
 873
 874    return;
 875
 876 }
 877
 878 //_______________________________________________________________________________
 879
 880 void AliAnalysisTaskFlavourJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){
 881
 882    //  D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
 883    // (expected detector resolution) on the left and right frm the D0 mass. Each band
 884    //  has a width of ~5 sigmas. Two band needed  for opening angle considerations
 885    TH2F* hDiffSideBand=(TH2F*)fmyOutput->FindObject("hDiffSideBand");
 886    TH3F* hPtJetWithDsb=(TH3F*)fmyOutput->FindObject("hPtJetWithDsb");
 887    THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");
 888
 889    Double_t deltaM=candDstar->DeltaInvMass();
 890    //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);
 891    Double_t z=Z(candDstar,jet);
 892    Double_t ptD=candDstar->Pt();
 893    Double_t ptjet=jet->Pt();
 894    Double_t dPhi=jet->Phi()-candDstar->Phi();
 895    Double_t deltaR=DeltaR(candDstar,jet);
 896    if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
 897    if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
 898
 899    Int_t isSideBand=1;//IsDzeroSideBand(candDstar);
 900    Double_t point[7]={z,dPhi,ptjet,ptD,deltaM,isSideBand, deltaR<fJetRadius ? 1 : 0};
 901    //fill the background histogram with the side bands or when looking at MC with the real background
 902    if(isSideBand==1){
 903       hDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background
 904       //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);
 905       hsDphiz->Fill(point,1.);
 906       if(deltaR<fJetRadius){  // evaluate in the near side
 907          //hzptDB->Fill(Z(candDstar,jet),deltaM,ptD);
 908          hPtJetWithDsb->Fill(ptjet,deltaM,ptD);
 909       }
 910    }
 911 }
 912
 913 //_______________________________________________________________________________
 914
 915 void AliAnalysisTaskFlavourJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg, AliEmcalJet *jet){
 916
 917    //need mass, deltaR, pt jet, ptD
 918    //two cases: D0 and Dstar
 919
 920    Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());
 921    if(!isselected) return;
 922
 923    Double_t ptD=candbg->Pt();
 924    Double_t ptjet=jet->Pt();
 925    Double_t deltaR=DeltaR(candbg,jet);
 926
 927    TH2F* hInvMassptDbg=(TH2F*)fmyOutput->FindObject("hInvMassptDbg");
 928    TH3F* hPtJetWithDsb=(TH3F*)fmyOutput->FindObject("hPtJetWithDsb");
 929
 930
 931    if(fCandidateType==kDstartoKpipi){
 932       AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;
 933       Double_t deltaM=dstarbg->DeltaInvMass();
 934       hInvMassptDbg->Fill(deltaM,ptD);
 935       if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,deltaM,ptD);
 936    }
 937
 938    if(fCandidateType==kD0toKpi){
 939       Double_t masses[2]={0.,0.};
 940       Int_t pdgdaughtersD0[2]={211,321};//pi,K
 941       Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
 942
 943       masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
 944       masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
 945
 946
 947       if(isselected==1 || isselected==3) {
 948          if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[0],ptD);
 949          hInvMassptDbg->Fill(masses[0],ptD);
 950       }
 951       if(isselected>=2) {
 952          if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[1],ptD);
 953          hInvMassptDbg->Fill(masses[1],ptD);
 954       }
 955
 956
 957    }
 958 }
 959
 960 //_______________________________________________________________________________
 961
 962 Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliVParticle *p1, AliVParticle *p2) const {
 963    //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
 964
 965    if(!p1 || !p2) return -1;
 966    Double_t phi1=p1->Phi(),eta1=p1->Eta();
 967    Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
 968
 969    Double_t dPhi=phi1-phi2;
 970    if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
 971    if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
 972
 973    Double_t dEta=eta1-eta2;
 974    Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
 975    return deltaR;
 976
 977 }
 978
 979 //_______________________________________________________________________________
 980
 981 Int_t AliAnalysisTaskFlavourJetCorrelations::IsDzeroSideBand(AliAODRecoCascadeHF *candDstar){
 982
 983    Double_t ptD=candDstar->Pt();
 984    Int_t bin = fCuts->PtBin(ptD);
 985    if (bin < 0){
 986       // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
 987       bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?
 988       return -1;  // out of bounds
 989    }
 990
 991    Double_t invM=candDstar->InvMassD0();
 992    Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 993
 994    Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
 995    Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
 996
 997    if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
 998    else return 0;
 999
1000 }