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) { //not used at the moment,uncomment return if you use
 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) || fUseMCInfo) {
 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     //If there's no background subtraction rhoval=0 and momentum is simply not took into account
 431     AliRhoParameter *rho = 0;
 432     Double_t rhoval = 0;
 433     TString sname("Rho");
 434     if (!sname.IsNull()) {
 435         rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(sname));
 436         if(rho) rhoval = rho->GetVal();
 437     }
 438
 439
 440    // we start with jets
 441    Double_t ejet   = 0;
 442    Double_t phiJet = 0;
 443    Double_t etaJet = 0;
 444    Double_t ptjet = 0;
 445    Double_t leadingJet =0;
 446    Double_t pointJ[6];
 447
 448    Int_t ntrarr=fTrackArr->GetEntriesFast();
 449    hNtrArr->Fill(ntrarr);
 450
 451    for(Int_t i=0;i<ntrarr;i++){
 452       AliVTrack *jtrack=static_cast<AliVTrack*>(fTrackArr->At(i));
 453       //reusing histograms
 454       hPtJetTrks->Fill(jtrack->Pt());
 455       hPhiJetTrks->Fill(jtrack->Phi());
 456       hEtaJetTrks->Fill(jtrack->Eta());
 457       hEjetTrks->Fill(jtrack->E());
 458    }
 459
 460
 461    //option to use only the leading jet
 462    if(fLeadingJetOnly){
 463       for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) {
 464          AliEmcalJet* jetL = (AliEmcalJet*)GetJetFromArray(iJetsL);
 465          ptjet   = jetL->Pt() - jetL->Area()*rhoval; //It takes into account the background subtraction
 466          if(ptjet>leadingJet ) leadingJet = ptjet;
 467
 468       }
 469    }
 470
 471    Int_t cntjet=0;
 472    //loop over jets and consider only the leading jet in the event
 473    for (Int_t iJets = 0; iJets<njets; iJets++) {
 474       fPmissing[0]=0;
 475       fPmissing[1]=0;
 476       fPmissing[2]=0;
 477
 478       //Printf("Jet N %d",iJets);
 479       AliEmcalJet* jet = (AliEmcalJet*)GetJetFromArray(iJets);
 480       //Printf("Corr task Accept Jet");
 481       if(!AcceptJet(jet)) {
 482          hstat->Fill(5);
 483          continue;
 484       }
 485
 486       //jets variables
 487       ejet   = jet->E();
 488       phiJet = jet->Phi();
 489       etaJet = jet->Eta();
 490       fPtJet = jet->Pt() - jet->Area()*rhoval; //It takes into account the background subtraction
 491       Double_t origPtJet=fPtJet;
 492
 493       pointJ[0] = phiJet;
 494       pointJ[1] = etaJet;
 495       pointJ[2] = ptjet - jet->Area()*rhoval; //It takes into account the background subtraction
 496       pointJ[3] = ejet;
 497       pointJ[4] = jet->GetNumberOfConstituents();
 498       pointJ[5] = jet->Area();
 499
 500       // choose the leading jet
 501       if(fLeadingJetOnly && (ejet<leadingJet)) continue;
 502       //used for normalization
 503       hstat->Fill(3);
 504       cntjet++;
 505       // fill energy, eta and phi of the jet
 506       hEjet   ->Fill(ejet);
 507       hPhiJet ->Fill(phiJet);
 508       hEtaJet ->Fill(etaJet);
 509       hPtJet  ->Fill(fPtJet);
 510       if(fJetOnlyMode && fSwitchOnSparses) hsJet->Fill(pointJ,1);
 511       //loop on jet particles
 512       Int_t ntrjet=  jet->GetNumberOfTracks();
 513       Double_t sumPtTracks=0,sumPzTracks=0;
 514       Int_t zg1jtrk=0;
 515       for(Int_t itrk=0;itrk<ntrjet;itrk++){
 516
 517          AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,fTrackArr);
 518          hdeltaRJetTracks->Fill(DeltaR(jet,jetTrk));
 519          Double_t ztrackj=Z(jetTrk,jet);
 520          hztracksinjet->Fill(ztrackj);
 521          sumPtTracks+=jetTrk->Pt();
 522          sumPzTracks+=jetTrk->Pz();
 523          if(ztrackj>1){
 524             zg1jtrk++;
 525          }
 526
 527          Double_t ztrackjTr=Z(jetTrk,jet,kTRUE);
 528          hztracksinjetT->Fill(ztrackjTr);
 529
 530       }//end loop on jet tracks
 531
 532       hNtrkjzNok->Fill(zg1jtrk);
 533       Double_t diffpt=origPtJet-sumPtTracks;
 534       Double_t diffpz=jet->Pz()-sumPzTracks;
 535       if(zg1jtrk>0){
 536          hDiffPtTrPtJzNok->Fill(diffpt);
 537          hDiffPzTrPtJzNok->Fill(diffpz);
 538
 539       }else{
 540          hDiffPtTrPtJzok->Fill(diffpt);
 541          hDiffPzTrPtJzok->Fill(diffpz);
 542       }
 543
 544       if(candidates==0){
 545
 546          hPtJetPerEvNoD->Fill(fPtJet);
 547       }
 548       if(!fJetOnlyMode) {
 549          //Printf("N candidates %d ", candidates);
 550          for(Int_t ic = 0; ic < candidates; ic++) {
 551             hstat->Fill(7);
 552             // D* candidates
 553             AliVParticle* charm=0x0;
 554             charm=(AliVParticle*)fCandidateArray->At(ic);
 555             if(!charm) continue;
 556             AliAODRecoDecayHF *charmdecay=(AliAODRecoDecayHF*) charm;
 557             fIsDInJet=IsDInJet(jet, charmdecay, kTRUE);
 558             if (fIsDInJet) FlagFlavour(jet);
 559             if (jet->TestFlavourTag(AliEmcalJet::kDStar) || jet->TestFlavourTag(AliEmcalJet::kD0)) hstat->Fill(4);
 560
 561             //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.
 562
 563             Double_t pjet[3];
 564             jet->PxPyPz(pjet);
 565              //It corrects the jet momentum if it was asked for jet background subtraction
 566              pjet[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
 567              pjet[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
 568              pjet[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
 569
 570             //It corrects the jet momentum due to D daughters out of the jet
 571             RecalculateMomentum(pjet,fPmissing);
 572             fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]); //recalculated jet pt
 573
 574
 575             //debugging histograms
 576             Double_t pmissing=TMath::Sqrt(fPmissing[0]*fPmissing[0]+fPmissing[1]*fPmissing[1]+fPmissing[2]*fPmissing[2]); //recalculated jet total momentum
 577             for(Int_t k=0;k<3;k++) ((TH1F*)fOutput->FindObject(Form("hMissP%d",k)))->Fill(fPmissing[k]);
 578
 579             ((TH1F*)fOutput->FindObject("hmissingp"))->Fill(pmissing);
 580             Double_t ptdiff=fPtJet-origPtJet;
 581             ((TH1F*)fOutput->FindObject("hDeltaPtJet"))->Fill(ptdiff);
 582             ((TH1F*)fOutput->FindObject("hRelDeltaPtJet"))->Fill(ptdiff/origPtJet);
 583
 584             FillHistogramsRecoJetCorr(charm, jet, aodEvent);
 585
 586          }//end loop on candidates
 587
 588          //retrieve side band background candidates for Dstar
 589          Int_t nsbcand = 0; if(fSideBandArray) nsbcand=fSideBandArray->GetEntriesFast();
 590
 591          for(Int_t ib=0;ib<nsbcand;ib++){
 592             if(fCandidateType==kDstartoKpipi && !fUseMCInfo){
 593                AliAODRecoCascadeHF *sbcand=(AliAODRecoCascadeHF*)fSideBandArray->At(ib);
 594                if(!sbcand) continue;
 595
 596                fIsDInJet=IsDInJet(jet, sbcand,kFALSE);
 597                Double_t pjet[3];
 598                jet->PxPyPz(pjet);
 599                 //It corrects the jet momentum if it was asked for jet background subtraction
 600                 pjet[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
 601                 pjet[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
 602                 pjet[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
 603
 604                //It corrects the jet momentum due to D daughters out of the jet
 605                RecalculateMomentum(pjet,fPmissing);
 606                fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]); //recalculated jet pt
 607
 608                SideBandBackground(sbcand,jet);
 609
 610             }
 611             if(fUseMCInfo){
 612
 613                AliAODRecoDecayHF* charmbg = 0x0;
 614                charmbg=(AliAODRecoDecayHF*)fSideBandArray->At(ib);
 615                if(!charmbg) continue;
 616                hstat->Fill(8);
 617                fIsDInJet=IsDInJet(jet, charmbg,kFALSE);
 618                if (fIsDInJet) {
 619                   FlagFlavour(jet); //this are backgroud HF jets, but flagged as signal at the moment. Can use the bkg flavour flag in the future. This info is not stored now a part in the jet
 620                   hstat->Fill(9);
 621                }
 622                Double_t pjet[3];
 623                jet->PxPyPz(pjet);
 624                //It corrects the jet momentum if it was asked for jet background subtraction
 625                pjet[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
 626                pjet[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
 627                pjet[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
 628
 629                //It corrects the jet momentum due to D daughters out of the jet
 630                RecalculateMomentum(pjet,fPmissing);
 631                fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]); //recalculated jet pt
 632
 633                MCBackground(charmbg,jet);
 634             }
 635          }
 636       }
 637    } // end of jet loop
 638
 639    hNJetPerEv->Fill(cntjet);
 640    if(candidates==0) hNJetPerEvNoD->Fill(cntjet);
 641    PostData(1,fOutput);
 642    return kTRUE;
 643 }
 644
 645 //_______________________________________________________________________________
 646
 647 void AliAnalysisTaskFlavourJetCorrelations::Terminate(Option_t*)
 648 {
 649    // The Terminate() function is the last function to be called during
 650    // a query. It always runs on the client, it can be used to present
 651    // the results graphically or save the results to file.
 652
 653    Info("Terminate"," terminate");
 654    AliAnalysisTaskSE::Terminate();
 655
 656    fOutput = dynamic_cast<TList*> (GetOutputData(1));
 657    if (!fOutput) {
 658       printf("ERROR: fOutput not available\n");
 659       return;
 660    }
 661 }
 662
 663 //_______________________________________________________________________________
 664
 665 void  AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){
 666    Float_t mass=0;
 667    Int_t abspdg=TMath::Abs(pdg);
 668
 669    mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
 670    // compute the Delta mass for the D*
 671    if(fCandidateType==kDstartoKpipi){
 672       Float_t mass1=0;
 673       mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 674       mass = mass-mass1;
 675    }
 676
 677    fMinMass = mass-range;
 678    fMaxMass = mass+range;
 679
 680    AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
 681    if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
 682 }
 683
 684 //_______________________________________________________________________________
 685
 686 void  AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){
 687    if(uplimit>lowlimit)
 688    {
 689       fMinMass = lowlimit;
 690       fMaxMass = uplimit;
 691    }
 692    else{
 693       printf("Error! Lower limit larger than upper limit!\n");
 694       fMinMass = uplimit - uplimit*0.2;
 695       fMaxMass = uplimit;
 696    }
 697 }
 698
 699 //_______________________________________________________________________________
 700
 701 Bool_t AliAnalysisTaskFlavourJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){
 702    if(nptbins>30) {
 703       AliInfo("Maximum number of bins allowed is 30!");
 704       return kFALSE;
 705    }
 706    if(!width) return kFALSE;
 707    for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
 708    return kTRUE;
 709 }
 710
 711 //_______________________________________________________________________________
 712
 713 Double_t AliAnalysisTaskFlavourJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet, Bool_t transverse) const{
 714    if(!part || !jet){
 715       printf("Particle or jet do not exist!\n");
 716       return -99;
 717    }
 718    Double_t p[3],pj[3];
 719    Bool_t okpp=part->PxPyPz(p);
 720    Bool_t okpj=jet->PxPyPz(pj);
 721
 722     //Background Subtraction
 723     //It corrects the each component of the jet momentum for Z calculation
 724     AliRhoParameter *rho = 0;
 725     Double_t rhoval = 0;
 726     TString sname("Rho");
 727     if (!sname.IsNull()) {
 728         rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(sname));
 729         if(rho){
 730             rhoval = rho->GetVal();
 731             pj[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
 732             pj[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
 733             pj[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
 734         }
 735     }
 736
 737
 738
 739    if(!okpp || !okpj){
 740       printf("Problems getting momenta\n");
 741       return -999;
 742    }
 743
 744    RecalculateMomentum(pj, fPmissing);
 745    if(transverse) return ZT(p,pj);
 746    else
 747    return Z(p,pj);
 748 }
 749
 750 //_______________________________________________________________________________
 751 Double_t AliAnalysisTaskFlavourJetCorrelations::Z(Double_t* p, Double_t *pj) const{
 752
 753    Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1]+pj[2]*pj[2];
 754    Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet2);
 755    return z;
 756 }
 757
 758
 759 //_______________________________________________________________________________
 760 Double_t AliAnalysisTaskFlavourJetCorrelations::ZT(Double_t* p, Double_t *pj) const{
 761
 762    Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1];
 763    Double_t z=(p[0]*pj[0]+p[1]*pj[1])/(pjet2);
 764    return z;
 765 }
 766
 767 //_______________________________________________________________________________
 768
 769 void AliAnalysisTaskFlavourJetCorrelations::RecalculateMomentum(Double_t* pj, const Double_t* pmissing) const {
 770
 771    pj[0]+=pmissing[0];
 772    pj[1]+=pmissing[1];
 773    pj[2]+=pmissing[2];
 774
 775 }
 776
 777 //_______________________________________________________________________________
 778
 779 Bool_t  AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){
 780
 781    // Statistics
 782    Int_t nbins=8;
 783    if(fUseMCInfo) nbins+=2;
 784
 785    TH1I* hstat=new TH1I("hstat","Statistics",nbins,-0.5,nbins-0.5);
 786    hstat->GetXaxis()->SetBinLabel(1,"N ev anal");
 787    hstat->GetXaxis()->SetBinLabel(2,"N ev sel");
 788    hstat->GetXaxis()->SetBinLabel(3,"N cand sel");
 789    hstat->GetXaxis()->SetBinLabel(4,"N jets");
 790    hstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
 791    hstat->GetXaxis()->SetBinLabel(6,"N jet rej");
 792    hstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet");
 793    hstat->GetXaxis()->SetBinLabel(8,"N jets & cand");
 794    if(fUseMCInfo) {
 795     hstat->GetXaxis()->SetBinLabel(3,"N Signal sel & jet");
 796     hstat->GetXaxis()->SetBinLabel(5,"N Signal in jet");
 797     hstat->GetXaxis()->SetBinLabel(9,"N Bkg sel & jet");
 798     hstat->GetXaxis()->SetBinLabel(10,"N Bkg in jet");
 799    }
 800    hstat->SetNdivisions(1);
 801    fOutput->Add(hstat);
 802
 803    const Int_t nbinsmass=200;
 804    const Int_t nbinsptjet=500;
 805    const Int_t nbinsptD=100;
 806    const Int_t nbinsz=100;
 807    const Int_t nbinsphi=200;
 808    const Int_t nbinseta=100;
 809
 810    //binning for THnSparse
 811    const Int_t nbinsSpsmass=50;
 812    const Int_t nbinsSpsptjet=100;
 813    const Int_t nbinsSpsptD=50;
 814    const Int_t nbinsSpsz=100;
 815    const Int_t nbinsSpsphi=100;
 816    const Int_t nbinsSpseta=60;
 817    const Int_t nbinsSpsContrib=100;
 818    const Int_t nbinsSpsA=100;
 819
 820    const Float_t ptjetlims[2]={0.,200.};
 821    const Float_t ptDlims[2]={0.,50.};
 822    const Float_t zlims[2]={0.,1.2};
 823    const Float_t philims[2]={0.,6.3};
 824    const Float_t etalims[2]={-1.5,1.5};
 825    const Int_t   nContriblims[2]={0,100};
 826    const Float_t arealims[2]={0.,2};
 827
 828    // jet related fistograms
 829
 830    TH1F* hEjetTrks      = new TH1F("hEjetTrks",  "Jet tracks energy distribution;Energy (GeV)",500,0,200);
 831    hEjetTrks->Sumw2();
 832    TH1F* hPhiJetTrks    = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)",  nbinsphi,philims[0],philims[1]);
 833    hPhiJetTrks->Sumw2();
 834    TH1F* hEtaJetTrks    = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta",  nbinseta,etalims[0],etalims[1]);
 835    hEtaJetTrks->Sumw2();
 836    TH1F* hPtJetTrks     = new TH1F("hPtJetTrks",  "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
 837    hPtJetTrks->Sumw2();
 838
 839    TH1F* hEjet      = new TH1F("hEjet",  "Jet energy distribution;Energy (GeV)",500,0,200);
 840    hEjet->Sumw2();
 841    TH1F* hPhiJet    = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)",  nbinsphi,philims[0],philims[1]);
 842    hPhiJet->Sumw2();
 843    TH1F* hEtaJet    = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
 844    hEtaJet->Sumw2();
 845    TH1F* hPtJet      = new TH1F("hPtJet",  "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
 846    hPtJet->Sumw2();
 847
 848    TH1F* hdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.);
 849    hdeltaRJetTracks->Sumw2();
 850
 851    TH1F* hNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000);
 852    hNtrArr->Sumw2();
 853
 854    TH1F *hNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);
 855    hNJetPerEv->Sumw2();
 856
 857
 858    fOutput->Add(hEjetTrks);
 859    fOutput->Add(hPhiJetTrks);
 860    fOutput->Add(hEtaJetTrks);
 861    fOutput->Add(hPtJetTrks);
 862    fOutput->Add(hEjet);
 863    fOutput->Add(hPhiJet);
 864    fOutput->Add(hEtaJet);
 865    fOutput->Add(hPtJet);
 866    fOutput->Add(hdeltaRJetTracks);
 867    fOutput->Add(hNtrArr);
 868    fOutput->Add(hNJetPerEv);
 869
 870    if(fJetOnlyMode && fSwitchOnSparses){
 871       //thnsparse for jets
 872       const Int_t nAxis=6;
 873       const Int_t nbinsSparse[nAxis]={nbinsSpsphi,nbinsSpseta, nbinsSpsptjet, nbinsSpsptjet,nbinsSpsContrib, nbinsSpsA};
 874       const Double_t minSparse[nAxis]={philims[0],etalims[0],ptjetlims[0],ptjetlims[0],static_cast<Double_t>(nContriblims[0]),arealims[0]};
 875       const Double_t
 876         maxSparse[nAxis]={philims[1],etalims[1],ptjetlims[1],ptjetlims[1],static_cast<Double_t>(nContriblims[1]),arealims[1]};
 877       THnSparseF *hsJet=new THnSparseF("hsJet","#phi, #eta, p_{T}^{jet}, E^{jet}, N contrib, Area", nAxis, nbinsSparse, minSparse, maxSparse);
 878       hsJet->Sumw2();
 879
 880       fOutput->Add(hsJet);
 881
 882    }
 883
 884    if(!fJetOnlyMode){
 885
 886       //debugging histograms
 887       TH1I* hControlDInJ=new TH1I("hControlDInJ","Checks D in Jet",8, -0.5,7.5);
 888       hControlDInJ->GetXaxis()->SetBinLabel(1,"DR In,1 daugh out");
 889       hControlDInJ->GetXaxis()->SetBinLabel(2,"DR In,2 daugh out");
 890       hControlDInJ->GetXaxis()->SetBinLabel(3,"DR In,3 daugh out");
 891       hControlDInJ->GetXaxis()->SetBinLabel(4,"Tot tracks non matched");
 892       hControlDInJ->GetXaxis()->SetBinLabel(5,"D0 daug missing");
 893       hControlDInJ->GetXaxis()->SetBinLabel(6,"soft pi missing");
 894       hControlDInJ->GetXaxis()->SetBinLabel(7,"IDprong diff GetDau");
 895       hControlDInJ->GetXaxis()->SetBinLabel(8,"still z>1 (cand)");
 896
 897       hControlDInJ->SetNdivisions(1);
 898       hControlDInJ->GetXaxis()->SetLabelSize(0.05);
 899       fOutput->Add(hControlDInJ);
 900
 901       TH1F *hmissingp=new TH1F("hmissingp", "Distribution of missing momentum (modulo);|p_{missing}|", 200, 0.,20.);
 902       fOutput->Add(hmissingp);
 903
 904       for(Int_t k=0;k<3;k++) {
 905          TH1F *hMissPi=new TH1F(Form("hMissP%d",k), Form("Missing p component {%d};p_{%d}",k,k), 400, -10.,10.);
 906          fOutput->Add(hMissPi);
 907       }
 908       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.);
 909
 910       fOutput->Add(hDeltaPtJet);
 911       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.);
 912       fOutput->Add(hRelDeltaPtJet);
 913
 914       TH1F* hzDinjet=new TH1F("hzDinjet","Z of candidates with daughters in jet;z",nbinsz,zlims[0],zlims[1]);
 915       fOutput->Add(hzDinjet);
 916       //frag func of tracks in the jet
 917       TH1F* hztracksinjet=new TH1F("hztracksinjet","Z of tracks in jet;z",nbinsz,zlims[0],zlims[1]);
 918       fOutput->Add(hztracksinjet);
 919
 920       //check jet and tracks
 921       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);
 922       fOutput->Add(hDiffPtTrPtJzok);
 923       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);
 924       fOutput->Add(hDiffPtTrPtJzNok);
 925       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);
 926       fOutput->Add(hDiffPzTrPtJzok);
 927       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);
 928       fOutput->Add(hDiffPzTrPtJzNok);
 929       TH1I* hNtrkjzNok=new TH1I("hNtrkjzNok", "Number of tracks in a jet with z>1;N^{tracks} (z>1)",5,0,5);
 930       fOutput->Add(hNtrkjzNok);
 931
 932       //calculate frag func with pt (simply ptD(or track)\cdot pt jet /ptjet^2)
 933       TH1F* hzDT=new TH1F("hzDT", Form("Z of D %s in jet in transverse components;(p_{T}^{D} dot p_{T}^{jet})/p_{T}^{jet}^{2} ", fUseMCInfo ? "(S+B)" : ""),nbinsz,zlims[0],zlims[1]);
 934       fOutput->Add(hzDT);
 935       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]);
 936       fOutput->Add(hztracksinjetT);
 937
 938       TH1I* hIDddaugh=new TH1I("hIDddaugh", "ID of daughters;ID", 2001,-1000,1000);
 939       fOutput->Add(hIDddaugh);
 940       TH1I* hIDddaughOut=new TH1I("hIDddaughOut", "ID of daughters not found in jet;ID", 2001,-1000,1000);
 941       fOutput->Add(hIDddaughOut);
 942       TH1I* hIDjetTracks=new TH1I("hIDjetTracks", "ID of jet tracks;ID", 2001,-1000,1000);
 943       fOutput->Add(hIDjetTracks);
 944
 945       TH1F* hDRdaughOut=new TH1F("hDRdaughOut", "#Delta R of daughters not belonging to the jet tracks (D in jet);#Delta R",200, 0.,2.);
 946       fOutput->Add(hDRdaughOut);
 947
 948
 949       if(fCandidateType==kDstartoKpipi)
 950       {
 951          if(fSwitchOnSB){
 952             TH2F* hDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
 953             hDiffSideBand->SetStats(kTRUE);
 954             hDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
 955             hDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
 956             hDiffSideBand->Sumw2();
 957             fOutput->Add(hDiffSideBand);
 958          }
 959
 960          TH1F* hPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
 961          hPtPion->SetStats(kTRUE);
 962          hPtPion->GetXaxis()->SetTitle("GeV/c");
 963          hPtPion->GetYaxis()->SetTitle("Entries");
 964          hPtPion->Sumw2();
 965          fOutput->Add(hPtPion);
 966
 967       }
 968       // D related histograms
 969       TH1F *hNDPerEvNoJet=new TH1F("hNDPerEvNoJet","Number of candidates per event with no jets; N candidate/ev with no jet", 20, 0., 20.);
 970       hNDPerEvNoJet->Sumw2();
 971       fOutput->Add(hNDPerEvNoJet);
 972
 973       TH1F *hptDPerEvNoJet=new TH1F("hptDPerEvNoJet","pt distribution of candidates per events with no jets; p_{t}^{D} (GeV/c)",nbinsptD, ptDlims[0],ptDlims[1]);
 974       hptDPerEvNoJet->Sumw2();
 975       fOutput->Add(hptDPerEvNoJet);
 976
 977       if(fSwitchOnSparses){
 978          const Int_t    nAxisD=4;
 979          const Int_t    nbinsSparseD[nAxisD]={nbinsSpseta,nbinsSpsphi,nbinsSpsptD,nbinsSpsmass};
 980          const Double_t minSparseD[nAxisD]  ={etalims[0],philims[0],ptDlims[0],fMinMass};
 981          const Double_t maxSparseD[nAxisD]  ={etalims[1],philims[1],ptDlims[1],fMaxMass};
 982          THnSparseF *hsDstandalone=new THnSparseF("hsDstandalone","#phi, #eta, p_{T}^{D}, and mass", nAxisD, nbinsSparseD, minSparseD, maxSparseD);
 983          hsDstandalone->Sumw2();
 984
 985          fOutput->Add(hsDstandalone);
 986       }
 987
 988       /*
 989       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]);
 990       hPtJetWithD->Sumw2();
 991       //for the MC this histogram is filled with the real background
 992       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]);
 993       hPtJetWithDsb->Sumw2();
 994
 995       fOutput->Add(hPtJetWithD);
 996       fOutput->Add(hPtJetWithDsb);
 997
 998       */
 999       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);
1000       hNJetPerEvNoD->Sumw2();
1001
1002       TH1F *hPtJetPerEvNoD=new TH1F("hPtJetPerEvNoD","pt distribution of jets per event with no D; #it{p}_{T}^{jet} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
1003       hPtJetPerEvNoD->Sumw2();
1004
1005       fOutput->Add(hNJetPerEvNoD);
1006       fOutput->Add(hPtJetPerEvNoD);
1007
1008       TH1F* hDeltaRD=new TH1F("hDeltaRD",Form("#Delta R distribution of D candidates %s selected;#Delta R", fUseMCInfo ? "(S+B)" : ""),200, 0.,10.);
1009       hDeltaRD->Sumw2();
1010       fOutput->Add(hDeltaRD);
1011
1012       TH3F* hDeltaRptDptj=new TH3F("hDeltaRptDptj",Form("#Delta R distribution of D candidates %s selected;#Delta R;#it{p}_{T}^{D} (GeV/c);#it{p}_{T}^{jet} (GeV/c)", fUseMCInfo ? "(S)" : ""),100, 0.,5.,nbinsptjet,ptjetlims[0],ptjetlims[1],nbinsptD, ptDlims[0],ptDlims[1]);
1013       hDeltaRptDptj->Sumw2();
1014       fOutput->Add(hDeltaRptDptj);
1015
1016       if(fUseMCInfo){
1017          TH3F* hDeltaRptDptjB=new TH3F("hDeltaRptDptjB",Form("#Delta R distribution of D candidates (B) selected;#Delta R;#it{p}_{T}^{D} (GeV/c);#it{p}_{T}^{jet} (GeV/c)"),100, 0.,5.,nbinsptjet,ptjetlims[0],ptjetlims[1],nbinsptD, ptDlims[0],ptDlims[1]);
1018          hDeltaRptDptjB->Sumw2();
1019          fOutput->Add(hDeltaRptDptjB);
1020       }
1021
1022       //background (side bands for the Dstar and like sign for D0)
1023       fJetRadius=GetJetContainer(0)->GetJetRadius();
1024       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]);
1025       hInvMassptD->SetStats(kTRUE);
1026       hInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
1027       hInvMassptD->GetYaxis()->SetTitle("#it{p}_{t}^{D} (GeV/c)");
1028       hInvMassptD->Sumw2();
1029
1030       fOutput->Add(hInvMassptD);
1031
1032       if(fUseMCInfo){
1033          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]);
1034          hInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
1035          hInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
1036          hInvMassptDbg->Sumw2();
1037          fOutput->Add(hInvMassptDbg);
1038
1039       }
1040
1041       if(fSwitchOnSparses){
1042          Double_t pi=TMath::Pi(), philims2[2];
1043          philims2[0]=-pi*1./2.;
1044          philims2[1]=pi*3./2.;
1045          THnSparseF *hsDphiz=0x0; //definition below according to the switches
1046
1047          if(fSwitchOnSB && fSwitchOnPhiAxis && fSwitchOnOutOfConeAxis){
1048             AliInfo("Creating a 9 axes container: This might cause large memory usage");
1049             const Int_t nAxis=9;
1050             const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsphi,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,2, 2, 2, 2};
1051             const Double_t minSparse[nAxis]={zlims[0],philims2[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5, -0.5,-0.5,-0.5};
1052             const Double_t maxSparse[nAxis]={zlims[1],philims2[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5, 1.5 , 1.5};
1053             fNAxesBigSparse=nAxis;
1054
1055             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);
1056          }
1057
1058          if(!fSwitchOnPhiAxis || !fSwitchOnOutOfConeAxis || !fSwitchOnSB){
1059             fSwitchOnPhiAxis=0;
1060             fSwitchOnOutOfConeAxis=0;
1061             fSwitchOnSB=0;
1062             if(fUseMCInfo){
1063                AliInfo("Creating a 7 axes container (MB background candidates)");
1064                const Int_t nAxis=7;
1065                const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,2, 2, 2};
1066                const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass, -0.5,-0.5,-0.5};
1067                const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5 , 1.5};
1068                fNAxesBigSparse=nAxis;
1069                hsDphiz=new THnSparseF("hsDphiz","Z vs p_{T}^{jet}, p_{T}^{D}, mass. Bkg?, D in EMCal acc?, jet in EMCal acc?", nAxis, nbinsSparse, minSparse, maxSparse);
1070
1071             } else{
1072                AliInfo("Creating a 6 axes container");
1073                const Int_t nAxis=6;
1074                const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass, 2, 2};
1075                const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5,-0.5};
1076                const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5};
1077                fNAxesBigSparse=nAxis;
1078
1079                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);
1080             }
1081          }
1082          if(!hsDphiz) AliFatal("No THnSparse created");
1083          hsDphiz->Sumw2();
1084
1085          fOutput->Add(hsDphiz);
1086       }
1087    }
1088    PostData(1, fOutput);
1089
1090    return kTRUE;
1091 }
1092
1093 //_______________________________________________________________________________
1094
1095 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet,  AliAODEvent* aodEvent){
1096
1097    Double_t ptD=candidate->Pt();
1098    Double_t deltaR=DeltaR(jet,candidate);
1099    Double_t phiD=candidate->Phi();
1100    Double_t deltaphi = jet->Phi()-phiD;
1101    if(deltaphi<=-(TMath::Pi())/2.) deltaphi = deltaphi+2.*(TMath::Pi());
1102    if(deltaphi>(3.*(TMath::Pi()))/2.) deltaphi = deltaphi-2.*(TMath::Pi());
1103    Double_t z=Z(candidate,jet);
1104    /*
1105    if(z>1) {
1106       ((TH1I*)fOutput->FindObject("hControlDInJ"))->Fill(7);
1107       Double_t pmissing=TMath::Sqrt(fPmissing[0]*fPmissing[0]+fPmissing[1]*fPmissing[1]+fPmissing[2]*fPmissing[2]);
1108
1109       for(Int_t k=0;k<3;k++) ((TH1F*)fOutput->FindObject(Form("hMissP%d",k)))->Fill(fPmissing[k]);
1110
1111       ((TH1F*)fOutput->FindObject("hmissingp"))->Fill(pmissing);
1112       Double_t ptdiff=fPtJet-jet->Pt();
1113       ((TH1F*)fOutput->FindObject("hDeltaPtJet"))->Fill(ptdiff);
1114       ((TH1F*)fOutput->FindObject("hRelDeltaPtJet"))->Fill(ptdiff/jet->Pt());
1115
1116
1117    }
1118    */
1119    if(fIsDInJet)((TH1F*)fOutput->FindObject("hzDT"))->Fill(Z(candidate,jet,kTRUE));
1120
1121    TH1F* hDeltaRD=(TH1F*)fOutput->FindObject("hDeltaRD");
1122    TH3F* hDeltaRptDptj=(TH3F*)fOutput->FindObject("hDeltaRptDptj");
1123
1124    hDeltaRD->Fill(deltaR);
1125    hDeltaRptDptj->Fill(deltaR,ptD,fPtJet);
1126
1127    Bool_t bDInEMCalAcc=InEMCalAcceptance(candidate);
1128    Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1129    if(fUseReco){
1130       if(fCandidateType==kD0toKpi) {
1131          AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;
1132
1133          FillHistogramsD0JetCorr(dzero,deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc, aodEvent);
1134
1135       }
1136
1137       if(fCandidateType==kDstartoKpipi) {
1138          AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;
1139          FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc);
1140
1141       }
1142    } else{ //generated level
1143       //AliInfo("Non reco");
1144       FillHistogramsMCGenDJetCorr(deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc);
1145    }
1146 }
1147
1148 //_______________________________________________________________________________
1149
1150 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){
1151
1152
1153    Double_t masses[2]={0.,0.};
1154    Int_t pdgdaughtersD0[2]={211,321};//pi,K
1155    Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1156
1157    masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1158    masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1159
1160    //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1161    THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
1162    Double_t *point=0x0;
1163    if(fNAxesBigSparse==9){
1164       point=new Double_t[9];
1165       point[0]=z;
1166       point[1]=dPhi;
1167       point[2]=ptj;
1168       point[3]=ptD;
1169       point[4]=masses[0];
1170       point[5]=0;
1171       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1172       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1173       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1174    }
1175    if(fNAxesBigSparse==6){
1176       point=new Double_t[6];
1177       point[0]=z;
1178       point[1]=ptj;
1179       point[2]=ptD;
1180       point[3]=masses[0];
1181       point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1182       point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1183 }
1184   if(fNAxesBigSparse==7){
1185       point=new Double_t[7];
1186       point[0]=z;
1187       point[1]=ptj;
1188       point[2]=ptD;
1189       point[3]=masses[0];
1190       point[4]=0;
1191       point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1192       point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1193
1194    }
1195
1196
1197    Printf("Candidate in FillHistogramsD0JetCorr IsA %s", (candidate->IsA())->GetName());
1198    Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
1199    if(isselected==1 || isselected==3) {
1200
1201       //if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[0],ptD);
1202
1203       FillMassHistograms(masses[0], ptD);
1204       if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) hsDphiz->Fill(point,1.);
1205    }
1206    if(isselected>=2) {
1207       //if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[1],ptD);
1208
1209       FillMassHistograms(masses[1], ptD);
1210       if(fNAxesBigSparse==9) point[4]=masses[1];
1211       if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[1];
1212       if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) hsDphiz->Fill(point,1.);
1213    }
1214
1215 }
1216
1217 //_______________________________________________________________________________
1218
1219 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi,  Double_t z, Double_t ptD, Double_t ptj, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc){
1220   //dPhi and z not used at the moment,but will be (re)added
1221
1222    AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
1223    Double_t deltamass= dstar->DeltaInvMass();
1224    //Double_t massD0= dstar->InvMassD0();
1225
1226
1227    TH1F* hPtPion=(TH1F*)fOutput->FindObject("hPtPion");
1228    hPtPion->Fill(softpi->Pt());
1229
1230    //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1231    THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
1232    Int_t isSB=0;//IsDzeroSideBand(dstar);
1233
1234    //Double_t point[]={z,dPhi,ptj,ptD,deltamass,static_cast<Double_t>(isSB), static_cast<Double_t>(fIsDInJet ? 1 : 0),bDInEMCalAcc,bJetInEMCalAcc};
1235    Double_t *point=0x0;
1236    if(fNAxesBigSparse==9){
1237       point=new Double_t[9];
1238       point[0]=z;
1239       point[1]=dPhi;
1240       point[2]=ptj;
1241       point[3]=ptD;
1242       point[4]=deltamass;
1243       point[5]=static_cast<Double_t>(isSB);
1244       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1245       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1246       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1247    }
1248    if(fNAxesBigSparse==6){
1249       point=new Double_t[6];
1250       point[0]=z;
1251       point[1]=ptj;
1252       point[2]=ptD;
1253       point[3]=deltamass;
1254       point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1255       point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1256    }
1257    if(fNAxesBigSparse==7){
1258       point=new Double_t[7];
1259       point[0]=z;
1260       point[1]=ptj;
1261       point[2]=ptD;
1262       point[3]=deltamass;
1263       point[4]=0;
1264       point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1265       point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1266    }
1267
1268    //if(fIsDInJet) hPtJetWithD->Fill(ptj,deltamass,ptD);
1269
1270    FillMassHistograms(deltamass, ptD);
1271    if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) hsDphiz->Fill(point,1.);
1272
1273 }
1274
1275 //_______________________________________________________________________________
1276
1277 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t z,Double_t ptD,Double_t ptjet, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc){
1278
1279    Double_t pdgmass=0;
1280    if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1281    if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
1282    //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1283    THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
1284    //Double_t point[9]={z,dPhi,ptjet,ptD,pdgmass,0, static_cast<Double_t>(fIsDInJet ? 1 : 0),bDInEMCalAcc,bJetInEMCalAcc};
1285
1286    Double_t *point=0x0;
1287    if(fNAxesBigSparse==9){
1288       point=new Double_t[9];
1289       point[0]=z;
1290       point[1]=dPhi;
1291       point[2]=ptjet;
1292       point[3]=ptD;
1293       point[4]=pdgmass;
1294       point[5]=0;
1295       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1296       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1297       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1298    }
1299    if(fNAxesBigSparse==6){
1300       point=new Double_t[6];
1301       point[0]=z;
1302       point[1]=ptjet;
1303       point[2]=ptD;
1304       point[3]=pdgmass;
1305       point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1306       point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1307    }
1308       if(fNAxesBigSparse==7){
1309       point=new Double_t[7];
1310       point[0]=z;
1311       point[1]=ptjet;
1312       point[2]=ptD;
1313       point[3]=pdgmass;
1314       point[4]=1;
1315       point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1316       point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1317    }
1318
1319
1320
1321    if(fNAxesBigSparse==9) point[4]=pdgmass;
1322    if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=pdgmass;
1323    if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) hsDphiz->Fill(point,1.);
1324    //if(fIsDInJet) {
1325    //  hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet
1326    //}
1327
1328 }
1329
1330 //_______________________________________________________________________________
1331
1332 void AliAnalysisTaskFlavourJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD){
1333
1334    if(fIsDInJet) {
1335       TH2F* hInvMassptD=(TH2F*)fOutput->FindObject("hInvMassptD");
1336       hInvMassptD->Fill(mass,ptD);
1337    }
1338 }
1339
1340 //________________________________________________________________________________
1341
1342 void AliAnalysisTaskFlavourJetCorrelations::FlagFlavour(AliEmcalJet *jet){
1343
1344    AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
1345    if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;
1346    if (fIsDInJet) jet->AddFlavourTag(tag);
1347
1348    return;
1349
1350 }
1351
1352 //_______________________________________________________________________________
1353
1354 void AliAnalysisTaskFlavourJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){
1355
1356    //  D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
1357    // (expected detector resolution) on the left and right frm the D0 mass. Each band
1358    //  has a width of ~5 sigmas. Two band needed  for opening angle considerations
1359    TH2F* hDiffSideBand=(TH2F*)fOutput->FindObject("hDiffSideBand");
1360    //TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
1361    THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
1362
1363    Bool_t bDInEMCalAcc=InEMCalAcceptance(candDstar);
1364    Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1365
1366    Double_t deltaM=candDstar->DeltaInvMass();
1367    //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);
1368    Double_t z=Z(candDstar,jet);
1369    Double_t ptD=candDstar->Pt();
1370
1371    Double_t dPhi=jet->Phi()-candDstar->Phi();
1372
1373    if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1374    if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1375
1376    Int_t isSideBand=1;//IsDzeroSideBand(candDstar);
1377    //if no SB analysis we should not enter here, so no need of checking the number of axes
1378    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)};
1379    //fill the background histogram with the side bands or when looking at MC with the real background
1380    if(isSideBand==1){
1381       hDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background
1382       //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);
1383       if(fSwitchOnSparses) hsDphiz->Fill(point,1.);
1384       //if(fIsDInJet){  // evaluate in the near side
1385       //         hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
1386       //}
1387    }
1388 }
1389
1390 //_______________________________________________________________________________
1391
1392 void AliAnalysisTaskFlavourJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg,AliEmcalJet* jet){
1393
1394    //need mass, deltaR, pt jet, ptD
1395    //two cases: D0 and Dstar
1396
1397    Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());
1398    if(!isselected) return;
1399
1400    Double_t ptD=candbg->Pt();
1401    Double_t deltaR=DeltaR(jet,candbg);
1402    Double_t phiD=candbg->Phi();
1403    Double_t deltaphi = jet->Phi()-phiD;
1404    if(deltaphi<=-(TMath::Pi())/2.) deltaphi = deltaphi+2.*(TMath::Pi());
1405    if(deltaphi>(3.*(TMath::Pi()))/2.) deltaphi = deltaphi-2.*(TMath::Pi());
1406    Double_t z=Z(candbg,jet);
1407
1408    if(fIsDInJet)((TH1F*)fOutput->FindObject("hzDT"))->Fill(Z(candbg,jet,kTRUE));
1409
1410    TH1F* hDeltaRD=(TH1F*)fOutput->FindObject("hDeltaRD");
1411    TH3F* hDeltaRptDptjB=(TH3F*)fOutput->FindObject("hDeltaRptDptjB");
1412
1413    hDeltaRD->Fill(deltaR);
1414    hDeltaRptDptjB->Fill(deltaR,ptD,fPtJet);
1415
1416    Bool_t bDInEMCalAcc=InEMCalAcceptance(candbg);
1417    Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1418
1419    TH2F* hInvMassptDbg=(TH2F*)fOutput->FindObject("hInvMassptDbg");
1420    //TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
1421
1422    THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
1423    Double_t *point=0x0;
1424    if(fNAxesBigSparse==9){
1425       point=new Double_t[9];
1426       point[0]=z;
1427       point[1]=deltaphi;
1428       point[2]=fPtJet;
1429       point[3]=ptD;
1430       point[4]=-999; //set below
1431       point[5]=1;
1432       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1433       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1434       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1435    }
1436
1437    if(fNAxesBigSparse==7){
1438       point=new Double_t[7];
1439       point[0]=z;
1440       point[1]=fPtJet;
1441       point[2]=ptD;
1442       point[3]=-999; //set below
1443       point[4]=1;
1444       point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1445       point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1446    }
1447
1448    if(fCandidateType==kDstartoKpipi){
1449       AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;
1450       Double_t deltaM=dstarbg->DeltaInvMass();
1451       hInvMassptDbg->Fill(deltaM,ptD);
1452       //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
1453       if(fNAxesBigSparse==9) point[4]=deltaM;
1454       if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=deltaM;
1455       if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) hsDphiz->Fill(point,1.);
1456    }
1457
1458    if(fCandidateType==kD0toKpi){
1459       Double_t masses[2]={0.,0.};
1460       Int_t pdgdaughtersD0[2]={211,321};//pi,K
1461       Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1462
1463       masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1464       masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1465
1466
1467       if(isselected==1 || isselected==3) {
1468          //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[0],ptD);
1469          hInvMassptDbg->Fill(masses[0],ptD);
1470          if(fNAxesBigSparse==9) point[4]=masses[0];
1471          if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[0];
1472          if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) hsDphiz->Fill(point,1.);
1473      }
1474       if(isselected>=2) {
1475          //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[1],ptD);
1476          hInvMassptDbg->Fill(masses[1],ptD);
1477          if(fNAxesBigSparse==9) point[4]=masses[1];
1478          if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[1];
1479          if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) hsDphiz->Fill(point,1.);
1480
1481       }
1482
1483
1484    }
1485 }
1486
1487 //_______________________________________________________________________________
1488
1489 Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliEmcalJet *p1, AliVParticle *p2) const {
1490    //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1491    //It recalculates the eta-phi values if it was asked for background subtraction of the jets
1492    if(!p1 || !p2) return -1;
1493
1494     Double_t phi1=p1->Phi(),eta1=p1->Eta();
1495
1496     //It subtracts the backgroud of jets if it was asked for it.
1497     TString sname("Rho");
1498     if (!sname.IsNull()) {
1499         AliRhoParameter *rho = 0;
1500         rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(sname));
1501         if(rho){
1502             Double_t pj[3];
1503             Bool_t okpj=p1->PxPyPz(pj);
1504             if(!okpj){
1505                 printf("Problems getting momenta\n");
1506                 return -999;
1507             }
1508             Double_t rhoval = rho->GetVal();
1509             pj[0] = p1->Px() - p1->Area()*(rhoval*TMath::Cos(p1->AreaPhi()));
1510             pj[1] = p1->Py() - p1->Area()*(rhoval*TMath::Sin(p1->AreaPhi()));
1511             pj[2] = p1->Pz() - p1->Area()*(rhoval*TMath::SinH(p1->AreaEta()));
1512             //Image of the function Arccos(px/pt) where pt = sqrt(px*px+py*py) is:
1513             //[0,pi]    if py > 0 and
1514             //[pi,2pi]  if py < 0
1515             phi1 = TMath::ACos(pj[0]/TMath::Sqrt(pj[0]*pj[0]+pj[1]*pj[1]));
1516             if(pj[1]<0) phi1 = 2*TMath::Pi() - phi1;
1517             eta1 = TMath::ASinH(pj[2]/TMath::Sqrt(pj[0]*pj[0]+pj[1]*pj[1]));
1518         }
1519     }
1520
1521    Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
1522
1523    Double_t dPhi=phi1-phi2;
1524    if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1525    if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1526
1527    Double_t dEta=eta1-eta2;
1528    Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
1529    return deltaR;
1530
1531 }
1532
1533 //_______________________________________________________________________________
1534
1535 Int_t AliAnalysisTaskFlavourJetCorrelations::IsDzeroSideBand(AliAODRecoCascadeHF *candDstar){
1536
1537    Double_t ptD=candDstar->Pt();
1538    Int_t bin = fCuts->PtBin(ptD);
1539    if (bin < 0){
1540       // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
1541       bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?
1542       return -1;  // out of bounds
1543    }
1544
1545    Double_t invM=candDstar->InvMassD0();
1546    Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1547
1548    Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
1549    Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
1550
1551    if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
1552    else return 0;
1553
1554 }
1555
1556 //_______________________________________________________________________________
1557
1558 Bool_t AliAnalysisTaskFlavourJetCorrelations::AreDaughtersInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Int_t*& daughOutOfJetID, AliAODTrack**& dtrks, Bool_t fillH){
1559    //returns true/false and the array of particles not found among jet constituents
1560
1561    TH1I* hControlDInJ=(TH1I*)fOutput->FindObject("hControlDInJ");
1562    TH1I* hIDddaugh   =(TH1I*)fOutput->FindObject("hIDddaugh");
1563    TH1I* hIDddaughOut=(TH1I*)fOutput->FindObject("hIDddaughOut");
1564    TH1I* hIDjetTracks=(TH1I*)fOutput->FindObject("hIDjetTracks");
1565
1566    Int_t daughtersID[3];
1567    Int_t ndaugh=0;
1568    Int_t dmatchedID[3]={0,0,0};
1569    Int_t countmatches=0;
1570    if (fCandidateType==kDstartoKpipi){
1571       AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)charm;
1572       AliAODRecoDecayHF2Prong* dzero=dstar->Get2Prong();
1573       daughtersID[0]=dzero->GetProngID(0);
1574       daughtersID[1]=dzero->GetProngID(1);
1575       daughtersID[2]=dstar->GetBachelor()->GetID();
1576       ndaugh=3;
1577
1578       dtrks=new AliAODTrack*[3];
1579       dtrks[0]=(AliAODTrack*)dzero->GetDaughter(0);
1580       dtrks[1]=(AliAODTrack*)dzero->GetDaughter(1);
1581       dtrks[2]=(AliAODTrack*)dstar->GetBachelor();
1582
1583       //check
1584       if(fillH){
1585          if(daughtersID[0]!=((AliAODTrack*)dtrks[0])->GetID() || daughtersID[1]!=((AliAODTrack*)dtrks[1])->GetID())  hControlDInJ->Fill(6);
1586
1587          hIDddaugh->Fill(daughtersID[0]);
1588          hIDddaugh->Fill(daughtersID[1]);
1589          hIDddaugh->Fill(daughtersID[2]);
1590
1591       }
1592       //Printf("ID daughters %d, %d, %d",daughtersID[0], daughtersID[1], daughtersID[2]);
1593    }
1594
1595    if (fCandidateType==kD0toKpi){
1596       daughtersID[0]=charm->GetProngID(0);
1597       daughtersID[1]=charm->GetProngID(1);
1598       ndaugh=2;
1599       if(fillH){
1600          hIDddaugh->Fill(daughtersID[0]);
1601          hIDddaugh->Fill(daughtersID[1]);
1602       }
1603       dtrks=new AliAODTrack*[2];
1604       dtrks[0]=(AliAODTrack*)charm->GetDaughter(0);
1605       dtrks[1]=(AliAODTrack*)charm->GetDaughter(1);
1606
1607    }
1608
1609    const Int_t cndaugh=ndaugh;
1610    daughOutOfJetID=new Int_t[cndaugh];
1611
1612    Int_t nchtrk=thejet->GetNumberOfTracks();
1613    for(Int_t itk=0;itk<nchtrk;itk++){
1614       AliVTrack* tkjet=((AliPicoTrack*)thejet->TrackAt(itk,fTrackArr))->GetTrack();
1615       Int_t idtkjet=tkjet->GetID();
1616       if(fillH) hIDjetTracks->Fill(idtkjet);
1617       for(Int_t id=0;id<ndaugh;id++){
1618          if(idtkjet==daughtersID[id]) {
1619             dmatchedID[id]=daughtersID[id]; //daughter which matches a track in the jet
1620             countmatches++;
1621
1622          }
1623
1624          if(countmatches==ndaugh) break;
1625       }
1626    }
1627    //reverse: include in the array the ID of daughters not matching
1628
1629    for(Int_t id=0;id<ndaugh;id++){
1630       if(dmatchedID[id]==0) {
1631          daughOutOfJetID[id]=daughtersID[id];
1632          if(fillH) hIDddaughOut->Fill(daughOutOfJetID[id]);
1633       }
1634       else daughOutOfJetID[id]=0;
1635    }
1636    if(fillH){
1637       if((ndaugh-countmatches) == 1) hControlDInJ->Fill(0);
1638       if((ndaugh-countmatches) == 2) hControlDInJ->Fill(1);
1639       if((ndaugh-countmatches) == 3) hControlDInJ->Fill(2);
1640    }
1641    if(ndaugh!=countmatches){
1642       return kFALSE;
1643    }
1644
1645    return kTRUE;
1646 }
1647
1648 //_______________________________________________________________________________
1649
1650 Bool_t AliAnalysisTaskFlavourJetCorrelations::IsDInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Bool_t fillH){
1651
1652    //check the conditions type:
1653    //type 0 : DeltaR < jet radius, don't check daughters (and don't correct pt jet)
1654    //type 1 : DeltaR < jet radius and check for all daughters among jet tracks, don't correct ptjet
1655    //type 2 (default) : DeltaR < jet radius and check for all daughters among jet tracks, if not present, correct the ptjet
1656
1657    TH1I* hControlDInJ=(TH1I*)fOutput->FindObject("hControlDInJ");
1658    TH1F* hDRdaughOut=(TH1F*)fOutput->FindObject("hDRdaughOut");
1659
1660    fPmissing[0]=0;
1661    fPmissing[1]=0;
1662    fPmissing[2]=0;
1663
1664    Bool_t testDeltaR=kFALSE;
1665    if(DeltaR(thejet,charm)<fJetRadius) testDeltaR=kTRUE;
1666
1667    Int_t* daughOutOfJet;
1668    AliAODTrack** charmDaugh;
1669    Bool_t testDaugh=AreDaughtersInJet(thejet, charm, daughOutOfJet,charmDaugh,fillH);
1670    if(testDaugh && testDeltaR) {
1671       //Z of candidates with daughters in the jet
1672       ((TH1F*)fOutput->FindObject("hzDinjet"))->Fill(Z(charm,thejet));
1673
1674    }
1675    if(!testDaugh && testDeltaR && fTypeDInJet==2){
1676
1677       Int_t ndaugh=3;
1678       if(fCandidateType==kD0toKpi) ndaugh=2;
1679       Int_t nOut=ndaugh;
1680
1681       for(Int_t id=0;id<ndaugh;id++){
1682          if(daughOutOfJet[id]!=0){ //non-matched daughter
1683             //get track and its momentum
1684             nOut--;
1685             if(fillH) {
1686                hControlDInJ->Fill(3);
1687                if(id<2) hControlDInJ->Fill(4);
1688                if(id==2)hControlDInJ->Fill(5);
1689                hDRdaughOut->Fill(DeltaR(thejet, charmDaugh[id]));
1690             }
1691             fPmissing[0]+=charmDaugh[id]->Px();
1692             fPmissing[1]+=charmDaugh[id]->Py();
1693             fPmissing[2]+=charmDaugh[id]->Pz();
1694          }
1695
1696       }
1697
1698       //now the D in within the jet
1699       testDaugh=kTRUE;
1700
1701    }
1702
1703    delete[] charmDaugh;
1704
1705    Bool_t result=0;
1706    switch(fTypeDInJet){
1707    case 0:
1708       result=testDeltaR;
1709       break;
1710    case 1:
1711       result=testDeltaR && testDaugh;
1712       break;
1713    case 2:
1714       result=testDeltaR && testDaugh;
1715       break;
1716    default:
1717       AliInfo("Selection type not specified, use 1");
1718       result=testDeltaR && testDaugh;
1719       break;
1720    }
1721  return result;
1722 }
1723
1724 Bool_t AliAnalysisTaskFlavourJetCorrelations::InEMCalAcceptance(AliVParticle *vpart){
1725    //check eta phi of a VParticle: return true if it is in the EMCal acceptance, false otherwise
1726
1727    Double_t phiEMCal[2]={1.405,3.135},etaEMCal[2]={-0.7,0.7};
1728    Bool_t binEMCal=kTRUE;
1729    Double_t phi=vpart->Phi(), eta=vpart->Eta();
1730    if(phi<phiEMCal[0] || phi>phiEMCal[1]) binEMCal=kFALSE;
1731    if(eta<etaEMCal[0] || eta>etaEMCal[1]) binEMCal=kFALSE;
1732    return binEMCal;
1733
1734
1735 }