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