PWGJE/FlavourJetTasks/AliAnalysisTaskFlavourJetCorrelations.cxx

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