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