PWGJE/FlavourJetTasks/AliAnalysisTaskFlavourJetCorrelations.cxx

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