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