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    hstat->Fill(0);
 317
 318    // fix for temporary bug in ESDfilter
 319    // the AODs with null vertex pointer didn't pass the PhysSel
 320    if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return kFALSE;
 321
 322    //Event selection
 323    Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
 324    TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
 325    if(!iseventselected) return kFALSE;
 326
 327    hstat->Fill(1);
 328
 329    //retrieve charm candidates selected
 330    Int_t candidates = 0;
 331    Int_t njets=GetJetContainer()->GetNJets();
 332
 333    if(!fJetOnlyMode) {
 334       candidates = fCandidateArray->GetEntriesFast();
 335
 336    //trigger on jets
 337    if(njets == 0) {
 338       hstat->Fill(6, candidates);
 339       hNDPerEvNoJet->Fill(candidates);
 340       for(Int_t iD=0;iD<candidates;iD++){
 341          AliVParticle* cand=(AliVParticle*)fCandidateArray->At(iD);
 342          if(!cand) continue;
 343          hptDPerEvNoJet->Fill(cand->Pt());
 344
 345       }
 346       return kFALSE;
 347
 348    }
 349
 350    //loop on candidates standalone (checking the candidates are there and their phi-eta distributions)
 351
 352    for(Int_t ic = 0; ic < candidates; ic++) {
 353
 354       // D* candidates
 355       AliAODRecoDecayHF* charm=0x0;
 356       AliAODRecoCascadeHF* dstar=0x0;
 357
 358
 359       charm=(AliAODRecoDecayHF*)fCandidateArray->At(ic);
 360       if(!charm) continue;
 361
 362       if(fCandidateType==kDstartoKpipi){
 363          dstar=(AliAODRecoCascadeHF*)fCandidateArray->At(ic);
 364          if(!dstar) continue;
 365       }
 366
 367       hstat->Fill(2);
 368
 369       Double_t candsparse[4]={charm->Eta(), charm->Phi(), charm->Pt(), 0};
 370
 371       if(fCandidateType==kDstartoKpipi) {
 372          if(fUseReco){
 373             Double_t deltamass= dstar->DeltaInvMass();
 374             candsparse[3]=deltamass;
 375          } else candsparse[3] = 0.145; //for generated
 376          hnspDstandalone->Fill(candsparse);
 377       }
 378       if(fCandidateType==kD0toKpi){
 379          if(fUseReco){
 380             AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)charm;
 381             Int_t isselected=fCuts->IsSelected(dzero,AliRDHFCuts::kAll,aodEvent);
 382             //not a further selection,just to get the value of isselected for the D0
 383             Double_t masses[2];
 384             Int_t pdgdaughtersD0[2]={211,321};//pi,K
 385             Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
 386
 387             masses[0]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
 388             masses[1]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
 389             if(isselected==1 || isselected==3) {
 390                candsparse[3]=masses[0];
 391                hnspDstandalone->Fill(candsparse);
 392             }
 393             if(isselected>=2){
 394                candsparse[3]=masses[1];
 395                hnspDstandalone->Fill(candsparse);
 396
 397             }
 398          } else { //generated
 399             Int_t pdg=((AliAODMCParticle*)charm)->GetPdgCode();
 400             candsparse[3]=TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
 401             hnspDstandalone->Fill(candsparse);
 402          }
 403       }
 404    }
 405    }
 406
 407    // we start with jets
 408    Double_t ejet   = 0;
 409    Double_t phiJet = 0;
 410    Double_t etaJet = 0;
 411    Double_t ptjet = 0;
 412    Double_t leadingJet =0;
 413    Double_t pointJ[6];
 414
 415    Int_t ntrarr=fTrackArr->GetEntriesFast();
 416    hNtrArr->Fill(ntrarr);
 417
 418    for(Int_t i=0;i<ntrarr;i++){
 419       AliVTrack *jtrack=static_cast<AliVTrack*>(fTrackArr->At(i));
 420       //reusing histograms
 421       hPtJetTrks->Fill(jtrack->Pt());
 422       hPhiJetTrks->Fill(jtrack->Phi());
 423       hEtaJetTrks->Fill(jtrack->Eta());
 424       hEjetTrks->Fill(jtrack->E());
 425    }
 426
 427
 428    //option to use only the leading jet
 429    if(fLeadingJetOnly){
 430       for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) {
 431          AliEmcalJet* jetL = (AliEmcalJet*)GetJetFromArray(iJetsL);
 432          ptjet   = jetL->Pt();
 433          if(ptjet>leadingJet ) leadingJet = ptjet;
 434
 435       }
 436    }
 437
 438    Int_t cntjet=0;
 439    //loop over jets and consider only the leading jet in the event
 440    for (Int_t iJets = 0; iJets<njets; iJets++) {
 441       //Printf("Jet N %d",iJets);
 442       AliEmcalJet* jet = (AliEmcalJet*)GetJetFromArray(iJets);
 443       //Printf("Corr task Accept Jet");
 444       if(!AcceptJet(jet)) {
 445          hstat->Fill(5);
 446          continue;
 447       }
 448
 449       //jets variables
 450       ejet   = jet->E();
 451       phiJet = jet->Phi();
 452       etaJet = jet->Eta();
 453       fPtJet = jet->Pt();
 454       Double_t origPtJet=fPtJet;
 455       pointJ[0] = phiJet;
 456       pointJ[1] = etaJet;
 457       pointJ[2] = ptjet;
 458       pointJ[3] = ejet;
 459       pointJ[4] = jet->GetNumberOfConstituents();
 460       pointJ[5] = jet->Area();
 461
 462       // choose the leading jet
 463       if(fLeadingJetOnly && (ejet<leadingJet)) continue;
 464       //used for normalization
 465       hstat->Fill(3);
 466       cntjet++;
 467       // fill energy, eta and phi of the jet
 468       hEjet   ->Fill(ejet);
 469       hPhiJet ->Fill(phiJet);
 470       hEtaJet ->Fill(etaJet);
 471       hPtJet  ->Fill(fPtJet);
 472       if(fJetOnlyMode) hsJet->Fill(pointJ,1);
 473       //loop on jet particles
 474       Int_t ntrjet=  jet->GetNumberOfTracks();
 475       for(Int_t itrk=0;itrk<ntrjet;itrk++){
 476
 477          AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,fTrackArr);
 478          hdeltaRJetTracks->Fill(DeltaR(jet,jetTrk));
 479
 480       }//end loop on jet tracks
 481
 482       if(candidates==0){
 483          hstat->Fill(7);
 484          hPtJetPerEvNoD->Fill(fPtJet);
 485       }
 486       if(!fJetOnlyMode) {
 487          //Printf("N candidates %d ", candidates);
 488          for(Int_t ic = 0; ic < candidates; ic++) {
 489
 490             // D* candidates
 491             AliVParticle* charm=0x0;
 492             charm=(AliVParticle*)fCandidateArray->At(ic);
 493             if(!charm) continue;
 494             AliAODRecoDecayHF *charmdecay=(AliAODRecoDecayHF*) charm;
 495             fIsDInJet=IsDInJet(jet, charmdecay, kTRUE);
 496             if (fIsDInJet) FlagFlavour(jet);
 497             if (jet->TestFlavourTag(AliEmcalJet::kDStar)) hstat->Fill(4);
 498
 499             //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.
 500
 501             Double_t pjet[3];
 502             jet->PxPyPz(pjet);
 503             RecalculateMomentum(pjet,fPmissing);
 504             fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]);
 505
 506             Double_t pmissing=TMath::Sqrt(fPmissing[0]*fPmissing[0]+fPmissing[1]*fPmissing[1]+fPmissing[2]*fPmissing[2]);
 507
 508             ((TH1F*)fOutput->FindObject("hmissingp"))->Fill(pmissing);
 509             Double_t ptdiff=fPtJet-origPtJet;
 510             ((TH1F*)fOutput->FindObject("hDeltaPtJet"))->Fill(ptdiff);
 511             ((TH1F*)fOutput->FindObject("hRelDeltaPtJet"))->Fill(ptdiff/origPtJet);
 512
 513             FillHistogramsRecoJetCorr(charm, jet, aodEvent);
 514
 515          }//end loop on candidates
 516
 517          //retrieve side band background candidates for Dstar
 518          Int_t nsbcand = 0; if(fSideBandArray) nsbcand=fSideBandArray->GetEntriesFast();
 519
 520          for(Int_t ib=0;ib<nsbcand;ib++){
 521             if(fCandidateType==kDstartoKpipi && !fUseMCInfo){
 522                AliAODRecoCascadeHF *sbcand=(AliAODRecoCascadeHF*)fSideBandArray->At(ib);
 523                if(!sbcand) continue;
 524
 525                fIsDInJet=IsDInJet(jet, sbcand,kFALSE);
 526                Double_t pjet[3];
 527                jet->PxPyPz(pjet);
 528                RecalculateMomentum(pjet,fPmissing);
 529                fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]);
 530
 531                SideBandBackground(sbcand,jet);
 532
 533             }
 534             if(fUseMCInfo){
 535                AliAODRecoDecayHF* charmbg = 0x0;
 536                charmbg=(AliAODRecoDecayHF*)fCandidateArray->At(ib);
 537                if(!charmbg) continue;
 538                fIsDInJet=IsDInJet(jet, charmbg,kFALSE);
 539
 540                Double_t pjet[3];
 541                jet->PxPyPz(pjet);
 542                RecalculateMomentum(pjet,fPmissing);
 543                fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]);
 544
 545                MCBackground(charmbg);
 546             }
 547          }
 548       }
 549    } // end of jet loop
 550
 551    hNJetPerEv->Fill(cntjet);
 552    if(candidates==0) hNJetPerEvNoD->Fill(cntjet);
 553    PostData(1,fOutput);
 554    return kTRUE;
 555 }
 556
 557 //_______________________________________________________________________________
 558
 559 void AliAnalysisTaskFlavourJetCorrelations::Terminate(Option_t*)
 560 {
 561    // The Terminate() function is the last function to be called during
 562    // a query. It always runs on the client, it can be used to present
 563    // the results graphically or save the results to file.
 564
 565    Info("Terminate"," terminate");
 566    AliAnalysisTaskSE::Terminate();
 567
 568    fOutput = dynamic_cast<TList*> (GetOutputData(1));
 569    if (!fOutput) {
 570       printf("ERROR: fOutput not available\n");
 571       return;
 572    }
 573 }
 574
 575 //_______________________________________________________________________________
 576
 577 void  AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){
 578    Float_t mass=0;
 579    Int_t abspdg=TMath::Abs(pdg);
 580
 581    mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
 582    // compute the Delta mass for the D*
 583    if(fCandidateType==kDstartoKpipi){
 584       Float_t mass1=0;
 585       mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 586       mass = mass-mass1;
 587    }
 588
 589    fMinMass = mass-range;
 590    fMaxMass = mass+range;
 591
 592    AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
 593    if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
 594 }
 595
 596 //_______________________________________________________________________________
 597
 598 void  AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){
 599    if(uplimit>lowlimit)
 600    {
 601       fMinMass = lowlimit;
 602       fMaxMass = uplimit;
 603    }
 604    else{
 605       printf("Error! Lower limit larger than upper limit!\n");
 606       fMinMass = uplimit - uplimit*0.2;
 607       fMaxMass = uplimit;
 608    }
 609 }
 610
 611 //_______________________________________________________________________________
 612
 613 Bool_t AliAnalysisTaskFlavourJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){
 614    if(nptbins>30) {
 615       AliInfo("Maximum number of bins allowed is 30!");
 616       return kFALSE;
 617    }
 618    if(!width) return kFALSE;
 619    for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
 620    return kTRUE;
 621 }
 622
 623 //_______________________________________________________________________________
 624
 625 Double_t AliAnalysisTaskFlavourJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet) const{
 626    if(!part || !jet){
 627       printf("Particle or jet do not exist!\n");
 628       return -99;
 629    }
 630    Double_t p[3],pj[3];
 631    Bool_t okpp=part->PxPyPz(p);
 632    Bool_t okpj=jet->PxPyPz(pj);
 633    if(!okpp || !okpj){
 634       printf("Problems getting momenta\n");
 635       return -999;
 636    }
 637
 638    RecalculateMomentum(pj, fPmissing);
 639    Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1]+pj[2]*pj[2];
 640    Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet2);
 641    return z;
 642 }
 643
 644 //_______________________________________________________________________________
 645
 646 void AliAnalysisTaskFlavourJetCorrelations::RecalculateMomentum(Double_t* pj, const Double_t* pmissing) const {
 647
 648    pj[0]+=pmissing[0];
 649    pj[1]+=pmissing[1];
 650    pj[2]+=pmissing[2];
 651
 652 }
 653
 654 //_______________________________________________________________________________
 655
 656 Bool_t  AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){
 657
 658    // Statistics
 659    TH1I* hstat=new TH1I("hstat","Statistics",8,-0.5,7.5);
 660    hstat->GetXaxis()->SetBinLabel(1,"N ev anal");
 661    hstat->GetXaxis()->SetBinLabel(2,"N ev sel");
 662    hstat->GetXaxis()->SetBinLabel(3,"N cand sel & jet");
 663    hstat->GetXaxis()->SetBinLabel(4,"N jets");
 664    hstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
 665    hstat->GetXaxis()->SetBinLabel(6,"N jet rej");
 666    hstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet");
 667    hstat->GetXaxis()->SetBinLabel(8,"N jets & !D");
 668    hstat->SetNdivisions(1);
 669    fOutput->Add(hstat);
 670
 671    const Int_t nbinsmass=200;
 672    const Int_t nbinsptjet=500;
 673    const Int_t nbinsptD=100;
 674    const Int_t nbinsz=100;
 675    const Int_t nbinsphi=200;
 676    const Int_t nbinseta=100;
 677    const Int_t nbinsContrib=100;
 678    const Int_t nbinsA=100;
 679
 680    const Float_t ptjetlims[2]={0.,200.};
 681    const Float_t ptDlims[2]={0.,50.};
 682    const Float_t zlims[2]={0.,1.2};
 683    const Float_t philims[2]={0.,6.3};
 684    const Float_t etalims[2]={-1.5,1.5};
 685    const Int_t   nContriblims[2]={0,100};
 686    const Float_t arealims[2]={0.,2};
 687
 688    // jet related fistograms
 689
 690    TH1F* hEjetTrks      = new TH1F("hEjetTrks",  "Jet tracks energy distribution;Energy (GeV)",500,0,200);
 691    hEjetTrks->Sumw2();
 692    TH1F* hPhiJetTrks    = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)",  nbinsphi,philims[0],philims[1]);
 693    hPhiJetTrks->Sumw2();
 694    TH1F* hEtaJetTrks    = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta",  nbinseta,etalims[0],etalims[1]);
 695    hEtaJetTrks->Sumw2();
 696    TH1F* hPtJetTrks     = new TH1F("hPtJetTrks",  "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
 697    hPtJetTrks->Sumw2();
 698
 699    TH1F* hEjet      = new TH1F("hEjet",  "Jet energy distribution;Energy (GeV)",500,0,200);
 700    hEjet->Sumw2();
 701    TH1F* hPhiJet    = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)",  nbinsphi,philims[0],philims[1]);
 702    hPhiJet->Sumw2();
 703    TH1F* hEtaJet    = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
 704    hEtaJet->Sumw2();
 705    TH1F* hPtJet      = new TH1F("hPtJet",  "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
 706    hPtJet->Sumw2();
 707
 708    TH1F* hdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.);
 709    hdeltaRJetTracks->Sumw2();
 710
 711    TH1F* hNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000);
 712    hNtrArr->Sumw2();
 713
 714    TH1F *hNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);
 715    hNJetPerEv->Sumw2();
 716
 717
 718    fOutput->Add(hEjetTrks);
 719    fOutput->Add(hPhiJetTrks);
 720    fOutput->Add(hEtaJetTrks);
 721    fOutput->Add(hPtJetTrks);
 722    fOutput->Add(hEjet);
 723    fOutput->Add(hPhiJet);
 724    fOutput->Add(hEtaJet);
 725    fOutput->Add(hPtJet);
 726    fOutput->Add(hdeltaRJetTracks);
 727    fOutput->Add(hNtrArr);
 728    fOutput->Add(hNJetPerEv);
 729
 730    if(fJetOnlyMode){
 731       //thnsparse for jets
 732       const Int_t nAxis=6;
 733       const Int_t nbinsSparse[nAxis]={nbinsphi,nbinseta, nbinsptjet, nbinsptjet,nbinsContrib, nbinsA};
 734       const Double_t minSparse[nAxis]={philims[0],etalims[0],ptjetlims[0],ptjetlims[0],nContriblims[0],arealims[0]};
 735       const Double_t
 736       maxSparse[nAxis]={philims[1],etalims[1],ptjetlims[1],ptjetlims[1],nContriblims[1],arealims[1]};
 737       THnSparseF *hsJet=new THnSparseF("hsJet","#phi, #eta, p_{T}^{jet}, E^{jet}, N contrib, Area", nAxis, nbinsSparse, minSparse, maxSparse);
 738       hsJet->Sumw2();
 739
 740       fOutput->Add(hsJet);
 741
 742    }
 743
 744    if(!fJetOnlyMode){
 745    TH1I* hControlDInJ=new TH1I("hControlDInJ","Checks D in Jet",8, -0.5,7.5);
 746    hControlDInJ->GetXaxis()->SetBinLabel(1,"DR In,1 daugh out");
 747    hControlDInJ->GetXaxis()->SetBinLabel(2,"DR In,2 daugh out");
 748    hControlDInJ->GetXaxis()->SetBinLabel(3,"DR In,3 daugh out");
 749    hControlDInJ->GetXaxis()->SetBinLabel(4,"Tot tracks non matched");
 750    hControlDInJ->GetXaxis()->SetBinLabel(5,"D0 daug missing");
 751    hControlDInJ->GetXaxis()->SetBinLabel(6,"soft pi missing");
 752    hControlDInJ->GetXaxis()->SetBinLabel(7,"IDprong diff GetDau");
 753    hControlDInJ->GetXaxis()->SetBinLabel(8,"still z>1 (cand)");
 754
 755    hControlDInJ->SetNdivisions(1);
 756    hControlDInJ->GetXaxis()->SetLabelSize(0.05);
 757    fOutput->Add(hControlDInJ);
 758
 759    TH1F *hmissingp=new TH1F("hmissingp", "Distribution of missing momentum (modulo)", 50, 0.,30.);
 760    fOutput->Add(hmissingp);
 761    TH1F *hDeltaPtJet=new TH1F("hDeltaPtJet", "Difference between the jet pt before and after correction", 50, 0.,30.);
 762    fOutput->Add(hDeltaPtJet);
 763    TH1F *hRelDeltaPtJet=new TH1F("hRelDeltaPtJet", "Difference between the jet pt before and after correction/ original pt", 200, 0.,2.);
 764    fOutput->Add(hRelDeltaPtJet);
 765
 766    TH1I* hIDddaugh=new TH1I("hIDddaugh", "ID of daughters", 2001,-1000,1000);
 767    fOutput->Add(hIDddaugh);
 768    TH1I* hIDddaughOut=new TH1I("hIDddaughOut", "ID of daughters not found in jet", 2001,-1000,1000);
 769    fOutput->Add(hIDddaughOut);
 770    TH1I* hIDjetTracks=new TH1I("hIDjetTracks", "ID of jet tracks", 2001,-1000,1000);
 771    fOutput->Add(hIDjetTracks);
 772
 773       if(fCandidateType==kDstartoKpipi)
 774       {
 775
 776          TH2F* hDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
 777          hDiffSideBand->SetStats(kTRUE);
 778          hDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
 779          hDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
 780          hDiffSideBand->Sumw2();
 781          fOutput->Add(hDiffSideBand);
 782
 783
 784          TH1F* hPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
 785          hPtPion->SetStats(kTRUE);
 786          hPtPion->GetXaxis()->SetTitle("GeV/c");
 787          hPtPion->GetYaxis()->SetTitle("Entries");
 788          hPtPion->Sumw2();
 789          fOutput->Add(hPtPion);
 790
 791       }
 792       // D related histograms
 793       TH1F *hNDPerEvNoJet=new TH1F("hNDPerEvNoJet","Number of candidates per event with no jets; N candidate/ev with no jet", 20, 0., 20.);
 794       hNDPerEvNoJet->Sumw2();
 795       fOutput->Add(hNDPerEvNoJet);
 796
 797       TH1F *hptDPerEvNoJet=new TH1F("hptDPerEvNoJet","pt distribution of candidates per events with no jets; p_{t}^{D} (GeV/c)",nbinsptD, ptDlims[0],ptDlims[1]);
 798       hptDPerEvNoJet->Sumw2();
 799       fOutput->Add(hptDPerEvNoJet);
 800
 801       const Int_t    nAxisD=4;
 802       const Int_t    nbinsSparseD[nAxisD]={nbinseta,nbinsphi,nbinsptD,nbinsmass};
 803       const Double_t minSparseD[nAxisD]  ={etalims[0],philims[0],ptDlims[0],fMinMass};
 804       const Double_t maxSparseD[nAxisD]  ={etalims[1],philims[1],ptDlims[1],fMaxMass};
 805       THnSparseF *hsDstandalone=new THnSparseF("hsDstandalone","#phi, #eta, p_{T}^{D}, and mass", nAxisD, nbinsSparseD, minSparseD, maxSparseD);
 806       hsDstandalone->Sumw2();
 807
 808       fOutput->Add(hsDstandalone);
 809
 810       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]);
 811       hPtJetWithD->Sumw2();
 812       //for the MC this histogram is filled with the real background
 813       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]);
 814       hPtJetWithDsb->Sumw2();
 815
 816       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);
 817       hNJetPerEvNoD->Sumw2();
 818
 819       TH1F *hPtJetPerEvNoD=new TH1F("hPtJetPerEvNoD","pt distribution of jets per event with no D; p_{T}^{jet} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
 820       hPtJetPerEvNoD->Sumw2();
 821
 822       fOutput->Add(hPtJetWithD);
 823       fOutput->Add(hPtJetWithDsb);
 824       fOutput->Add(hNJetPerEvNoD);
 825       fOutput->Add(hPtJetPerEvNoD);
 826
 827       TH1F* hDeltaRD=new TH1F("hDeltaRD","#Delta R distribution of D candidates selected;#Delta R",200, 0.,10.);
 828       hDeltaRD->Sumw2();
 829       fOutput->Add(hDeltaRD);
 830
 831       //background (side bands for the Dstar and like sign for D0)
 832       fJetRadius=GetJetContainer(0)->GetJetRadius();
 833       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]);
 834       hInvMassptD->SetStats(kTRUE);
 835       hInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
 836       hInvMassptD->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
 837       hInvMassptD->Sumw2();
 838
 839       fOutput->Add(hInvMassptD);
 840
 841       if(fUseMCInfo){
 842          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]);
 843          hInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
 844          hInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
 845          hInvMassptDbg->Sumw2();
 846          fOutput->Add(hInvMassptDbg);
 847
 848       }
 849       Double_t pi=TMath::Pi(), philims2[2];
 850       philims2[0]=-pi*1./2.;
 851       philims2[1]=pi*3./2.;
 852       const Int_t nAxis=7;
 853       const Int_t nbinsSparse[nAxis]={nbinsz,nbinsphi,nbinsptjet,nbinsptD,nbinsmass,2, 2};
 854       const Double_t minSparse[nAxis]={zlims[0],philims2[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5, -0.5};
 855       const Double_t maxSparse[nAxis]={zlims[1],philims2[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5};
 856       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);
 857       hsDphiz->Sumw2();
 858
 859       fOutput->Add(hsDphiz);
 860    }
 861    PostData(1, fOutput);
 862
 863    return kTRUE;
 864 }
 865
 866 //_______________________________________________________________________________
 867
 868 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet,  AliAODEvent* aodEvent){
 869
 870    Double_t ptD=candidate->Pt();
 871    Double_t deltaR=DeltaR(candidate,jet);
 872    Double_t phiD=candidate->Phi();
 873    Double_t deltaphi = jet->Phi()-phiD;
 874    if(deltaphi<=-(TMath::Pi())/2.) deltaphi = deltaphi+2.*(TMath::Pi());
 875    if(deltaphi>(3.*(TMath::Pi()))/2.) deltaphi = deltaphi-2.*(TMath::Pi());
 876    Double_t z=Z(candidate,jet);
 877    if(z>1) ((TH1I*)fOutput->FindObject("hControlDInJ"))->Fill(7);
 878
 879    TH1F* hDeltaRD=(TH1F*)fOutput->FindObject("hDeltaRD");
 880    hDeltaRD->Fill(deltaR);
 881    if(fUseReco){
 882       if(fCandidateType==kD0toKpi) {
 883          AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;
 884
 885          FillHistogramsD0JetCorr(dzero,deltaphi,phiD,z,ptD,fPtJet, aodEvent);
 886
 887       }
 888
 889       if(fCandidateType==kDstartoKpipi) {
 890          AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;
 891          FillHistogramsDstarJetCorr(dstar,deltaphi,phiD,z,ptD,fPtJet);
 892
 893       }
 894    } else{ //generated level
 895       //AliInfo("Non reco");
 896       FillHistogramsMCGenDJetCorr(deltaphi,phiD,z,ptD,fPtJet);
 897    }
 898 }
 899
 900 //_______________________________________________________________________________
 901
 902 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t phiD, Double_t z, Double_t ptD, Double_t ptj, AliAODEvent* aodEvent){
 903
 904
 905    Double_t masses[2]={0.,0.};
 906    Int_t pdgdaughtersD0[2]={211,321};//pi,K
 907    Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
 908
 909    masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
 910    masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
 911
 912    TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
 913    THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
 914    Double_t point[8]={z,dPhi,ptj,ptD,masses[0],0, fIsDInJet ? 1 : 0,phiD};
 915    Printf("Candidate in FillHistogramsD0JetCorr IsA %s", (candidate->IsA())->GetName());
 916    Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
 917    if(isselected==1 || isselected==3) {
 918
 919       if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[0],ptD);
 920
 921       FillMassHistograms(masses[0], ptD);
 922       hsDphiz->Fill(point,1.);
 923    }
 924    if(isselected>=2) {
 925       if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[1],ptD);
 926
 927       FillMassHistograms(masses[1], ptD);
 928       point[4]=masses[1];
 929       hsDphiz->Fill(point,1.);
 930    }
 931
 932 }
 933
 934 //_______________________________________________________________________________
 935
 936 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi,  Double_t phiD, Double_t z, Double_t ptD, Double_t ptj){
 937   //dPhi and z not used at the moment,but will be (re)added
 938
 939    AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
 940    Double_t deltamass= dstar->DeltaInvMass();
 941    //Double_t massD0= dstar->InvMassD0();
 942
 943
 944    TH1F* hPtPion=(TH1F*)fOutput->FindObject("hPtPion");
 945    hPtPion->Fill(softpi->Pt());
 946
 947    TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
 948    THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
 949    Int_t isSB=0;//IsDzeroSideBand(dstar);
 950
 951    Double_t point[8]={z,dPhi,ptj,ptD,deltamass,isSB, fIsDInJet ? 1 : 0,phiD};
 952
 953    if(fIsDInJet) hPtJetWithD->Fill(ptj,deltamass,ptD);
 954
 955    FillMassHistograms(deltamass, ptD);
 956    hsDphiz->Fill(point,1.);
 957
 958 }
 959
 960 //_______________________________________________________________________________
 961
 962 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t phiD,Double_t z,Double_t ptD,Double_t ptjet){
 963
 964    Double_t pdgmass=0;
 965    TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
 966    THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
 967    Double_t point[8]={z,dPhi,ptjet,ptD,pdgmass,0, fIsDInJet ? 1 : 0,phiD};
 968
 969    if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 970    if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
 971    point[4]=pdgmass;
 972    hsDphiz->Fill(point,1.);
 973    if(fIsDInJet) {
 974      hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet
 975    }
 976
 977 }
 978
 979 //_______________________________________________________________________________
 980
 981 void AliAnalysisTaskFlavourJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD){
 982
 983    if(fIsDInJet) {
 984       TH2F* hInvMassptD=(TH2F*)fOutput->FindObject("hInvMassptD");
 985       hInvMassptD->Fill(mass,ptD);
 986    }
 987 }
 988
 989 //________________________________________________________________________________
 990
 991 void AliAnalysisTaskFlavourJetCorrelations::FlagFlavour(AliEmcalJet *jet){
 992
 993    AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
 994    if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;
 995    if (fIsDInJet) jet->AddFlavourTag(tag);
 996
 997    return;
 998
 999 }
1000
1001 //_______________________________________________________________________________
1002
1003 void AliAnalysisTaskFlavourJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){
1004
1005    //  D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
1006    // (expected detector resolution) on the left and right frm the D0 mass. Each band
1007    //  has a width of ~5 sigmas. Two band needed  for opening angle considerations
1008    TH2F* hDiffSideBand=(TH2F*)fOutput->FindObject("hDiffSideBand");
1009    TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
1010    THnSparseF* hsDphiz=(THnSparseF*)fOutput->FindObject("hsDphiz");
1011
1012    Double_t deltaM=candDstar->DeltaInvMass();
1013    //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);
1014    Double_t z=Z(candDstar,jet);
1015    Double_t ptD=candDstar->Pt();
1016
1017    Double_t dPhi=jet->Phi()-candDstar->Phi();
1018
1019    if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1020    if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1021
1022    Int_t isSideBand=1;//IsDzeroSideBand(candDstar);
1023    Double_t point[7]={z,dPhi,fPtJet,ptD,deltaM,isSideBand, fIsDInJet ? 1 : 0};
1024    //fill the background histogram with the side bands or when looking at MC with the real background
1025    if(isSideBand==1){
1026       hDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background
1027       //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);
1028       hsDphiz->Fill(point,1.);
1029       if(fIsDInJet){  // evaluate in the near side
1030          hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
1031       }
1032    }
1033 }
1034
1035 //_______________________________________________________________________________
1036
1037 void AliAnalysisTaskFlavourJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg){
1038
1039    //need mass, deltaR, pt jet, ptD
1040    //two cases: D0 and Dstar
1041
1042    Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());
1043    if(!isselected) return;
1044
1045    Double_t ptD=candbg->Pt();
1046
1047    TH2F* hInvMassptDbg=(TH2F*)fOutput->FindObject("hInvMassptDbg");
1048    TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
1049
1050
1051    if(fCandidateType==kDstartoKpipi){
1052       AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;
1053       Double_t deltaM=dstarbg->DeltaInvMass();
1054       hInvMassptDbg->Fill(deltaM,ptD);
1055       if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
1056    }
1057
1058    if(fCandidateType==kD0toKpi){
1059       Double_t masses[2]={0.,0.};
1060       Int_t pdgdaughtersD0[2]={211,321};//pi,K
1061       Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1062
1063       masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1064       masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1065
1066
1067       if(isselected==1 || isselected==3) {
1068          if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[0],ptD);
1069          hInvMassptDbg->Fill(masses[0],ptD);
1070       }
1071       if(isselected>=2) {
1072          if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[1],ptD);
1073          hInvMassptDbg->Fill(masses[1],ptD);
1074       }
1075
1076
1077    }
1078 }
1079
1080 //_______________________________________________________________________________
1081
1082 Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliVParticle *p1, AliVParticle *p2) const {
1083    //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1084
1085    if(!p1 || !p2) return -1;
1086    Double_t phi1=p1->Phi(),eta1=p1->Eta();
1087    Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
1088
1089    Double_t dPhi=phi1-phi2;
1090    if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1091    if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1092
1093    Double_t dEta=eta1-eta2;
1094    Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
1095    return deltaR;
1096
1097 }
1098
1099 //_______________________________________________________________________________
1100
1101 Int_t AliAnalysisTaskFlavourJetCorrelations::IsDzeroSideBand(AliAODRecoCascadeHF *candDstar){
1102
1103    Double_t ptD=candDstar->Pt();
1104    Int_t bin = fCuts->PtBin(ptD);
1105    if (bin < 0){
1106       // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
1107       bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?
1108       return -1;  // out of bounds
1109    }
1110
1111    Double_t invM=candDstar->InvMassD0();
1112    Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1113
1114    Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
1115    Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
1116
1117    if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
1118    else return 0;
1119
1120 }
1121
1122 //_______________________________________________________________________________
1123
1124 Bool_t AliAnalysisTaskFlavourJetCorrelations::AreDaughtersInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Int_t*& daughOutOfJetID, AliAODTrack**& dtrks, Bool_t fillH){
1125    //returns true/false and the array of particles not found among jet constituents
1126
1127    TH1I* hControlDInJ=(TH1I*)fOutput->FindObject("hControlDInJ");
1128    TH1I* hIDddaugh   =(TH1I*)fOutput->FindObject("hIDddaugh");
1129    TH1I* hIDddaughOut=(TH1I*)fOutput->FindObject("hIDddaughOut");
1130    TH1I* hIDjetTracks=(TH1I*)fOutput->FindObject("hIDjetTracks");
1131
1132    Int_t daughtersID[3];
1133    Int_t ndaugh=0;
1134    Int_t dmatchedID[3]={0,0,0};
1135    Int_t countmatches=0;
1136    if (fCandidateType==kDstartoKpipi){
1137       AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)charm;
1138       AliAODRecoDecayHF2Prong* dzero=dstar->Get2Prong();
1139       daughtersID[0]=dzero->GetProngID(0);
1140       daughtersID[1]=dzero->GetProngID(1);
1141       daughtersID[2]=dstar->GetBachelor()->GetID();
1142       ndaugh=3;
1143
1144       dtrks=new AliAODTrack*[3];
1145       dtrks[0]=(AliAODTrack*)dzero->GetDaughter(0);
1146       dtrks[1]=(AliAODTrack*)dzero->GetDaughter(1);
1147       dtrks[2]=(AliAODTrack*)dstar->GetBachelor();
1148
1149       //check
1150       if(fillH){
1151          if(daughtersID[0]!=((AliAODTrack*)dtrks[0])->GetID() || daughtersID[1]!=((AliAODTrack*)dtrks[1])->GetID())  hControlDInJ->Fill(6);
1152
1153          hIDddaugh->Fill(daughtersID[0]);
1154          hIDddaugh->Fill(daughtersID[1]);
1155          hIDddaugh->Fill(daughtersID[2]);
1156
1157       }
1158       //Printf("ID daughters %d, %d, %d",daughtersID[0], daughtersID[1], daughtersID[2]);
1159    }
1160
1161    if (fCandidateType==kD0toKpi){
1162       daughtersID[0]=charm->GetProngID(0);
1163       daughtersID[1]=charm->GetProngID(1);
1164       ndaugh=2;
1165       if(fillH){
1166          hIDddaugh->Fill(daughtersID[0]);
1167          hIDddaugh->Fill(daughtersID[1]);
1168       }
1169       dtrks=new AliAODTrack*[2];
1170       dtrks[0]=(AliAODTrack*)charm->GetDaughter(0);
1171       dtrks[1]=(AliAODTrack*)charm->GetDaughter(1);
1172
1173    }
1174
1175    const Int_t cndaugh=ndaugh;
1176    daughOutOfJetID=new Int_t[cndaugh];
1177
1178    Int_t nchtrk=thejet->GetNumberOfTracks();
1179    for(Int_t itk=0;itk<nchtrk;itk++){
1180       AliVTrack* tkjet=((AliPicoTrack*)thejet->TrackAt(itk,fTrackArr))->GetTrack();
1181       Int_t idtkjet=tkjet->GetID();
1182       if(fillH) hIDjetTracks->Fill(idtkjet);
1183       for(Int_t id=0;id<ndaugh;id++){
1184          if(idtkjet==daughtersID[id]) {
1185             dmatchedID[id]=daughtersID[id]; //daughter which matches a track in the jet
1186             countmatches++;
1187
1188          }
1189
1190          if(countmatches==ndaugh) break;
1191       }
1192    }
1193    //reverse: include in the array the ID of daughters not matching
1194
1195    for(Int_t id=0;id<ndaugh;id++){
1196       if(dmatchedID[id]==0) {
1197          daughOutOfJetID[id]=daughtersID[id];
1198          if(fillH) hIDddaughOut->Fill(daughOutOfJetID[id]);
1199       }
1200       else daughOutOfJetID[id]=0;
1201    }
1202    if(fillH){
1203       if((ndaugh-countmatches) == 1) hControlDInJ->Fill(0);
1204       if((ndaugh-countmatches) == 2) hControlDInJ->Fill(1);
1205       if((ndaugh-countmatches) == 3) hControlDInJ->Fill(2);
1206    }
1207    if(ndaugh!=countmatches){
1208       return kFALSE;
1209    }
1210
1211    return kTRUE;
1212 }
1213
1214 //_______________________________________________________________________________
1215
1216 Bool_t AliAnalysisTaskFlavourJetCorrelations::IsDInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Bool_t fillH){
1217
1218    //check the conditions type:
1219    //type 0 : DeltaR < jet radius, don't check daughters (and don't correct pt jet)
1220    //type 1 : DeltaR < jet radius and check for all daughters among jet tracks, don't correct ptjet
1221    //type 2 (default) : DeltaR < jet radius and check for all daughters among jet tracks, if not present, correct the ptjet
1222
1223    TH1I* hControlDInJ=(TH1I*)fOutput->FindObject("hControlDInJ");
1224
1225    Bool_t testDeltaR=kFALSE;
1226    if(DeltaR(thejet,charm)<fJetRadius) testDeltaR=kTRUE;
1227
1228    Int_t* daughOutOfJet;
1229    AliAODTrack** charmDaugh;
1230    Bool_t testDaugh=AreDaughtersInJet(thejet, charm, daughOutOfJet,charmDaugh,fillH);
1231    if(!testDaugh && testDeltaR && fTypeDInJet==2){
1232
1233       Int_t ndaugh=3;
1234       if(fCandidateType==kD0toKpi) ndaugh=2;
1235       Int_t nOut=ndaugh;
1236
1237       fPmissing[0]=0;
1238       fPmissing[1]=0;
1239       fPmissing[2]=0;
1240       for(Int_t id=0;id<ndaugh;id++){
1241          if(daughOutOfJet[id]!=0){ //non-matched daughter
1242             //get track and its momentum
1243             nOut--;
1244             if(fillH) {
1245                hControlDInJ->Fill(3);
1246                if(id<2) hControlDInJ->Fill(4);
1247                if(id==2)hControlDInJ->Fill(5);
1248             }
1249             fPmissing[0]+=charmDaugh[id]->Px();
1250             fPmissing[1]+=charmDaugh[id]->Py();
1251             fPmissing[2]+=charmDaugh[id]->Pz();
1252          }
1253
1254       }
1255
1256       //now the D in within the jet
1257       testDaugh=kTRUE;
1258
1259    }
1260
1261    delete[] charmDaugh;
1262
1263    Bool_t result=0;
1264    switch(fTypeDInJet){
1265    case 0:
1266       result=testDeltaR;
1267    case 1:
1268       result=testDeltaR && testDaugh;
1269    case 2:
1270       result=testDeltaR && testDaugh;
1271    }
1272  return result;
1273 }