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