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 nbinsptjet=500;
 918    const Int_t nbinsptjetBG=400;
 919    const Int_t nbinsptD=100;
 920    const Int_t nbinsz=100;
 921    const Int_t nbinsphi=200;
 922    const Int_t nbinseta=100;
 923
 924    //binning for THnSparse
 925    const Int_t nbinsSpsmass=50;
 926    const Int_t nbinsSpsptjet=100;
 927    const Int_t nbinsSpsptD=50;
 928    const Int_t nbinsSpsz=100;
 929    const Int_t nbinsSpsphi=100;
 930    const Int_t nbinsSpseta=60;
 931    const Int_t nbinsSpsContrib=100;
 932    const Int_t nbinsSpsA=100;
 933
 934    const Float_t ptjetlims[2]={0.,200.};
 935    const Float_t ptjetlimsBG[2]={-200.,200.};
 936    const Float_t ptDlims[2]={0.,50.};
 937    const Float_t zlims[2]={0.,1.2};
 938    const Float_t philims[2]={0.,6.3};
 939    const Float_t etalims[2]={-1.5,1.5};
 940    const Int_t   nContriblims[2]={0,100};
 941    const Float_t arealims[2]={0.,2};
 942
 943    // jet related fistograms
 944
 945    fhEjetTrks      = new TH1F("hEjetTrks",  "Jet tracks energy distribution;Energy (GeV)",500,0,200);
 946    fhEjetTrks->Sumw2();
 947    fhPhiJetTrks    = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)",  nbinsphi,philims[0],philims[1]);
 948    fhPhiJetTrks->Sumw2();
 949    fhEtaJetTrks    = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta",  nbinseta,etalims[0],etalims[1]);
 950    fhEtaJetTrks->Sumw2();
 951    fhPtJetTrks     = new TH1F("hPtJetTrks",  "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
 952    fhPtJetTrks->Sumw2();
 953
 954    fhEjet      = new TH1F("hEjet",  "Jet energy distribution;Energy (GeV)",500,0,200);
 955    fhEjet->Sumw2();
 956    fhPhiJet    = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)",  nbinsphi,philims[0],philims[1]);
 957    fhPhiJet->Sumw2();
 958    fhEtaJet    = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
 959    fhEtaJet->Sumw2();
 960     if(!fJetCont->GetRhoName().IsNull()){ fhPtJet      = new TH1F("hPtJet",  "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjetBG,ptjetlimsBG[0],ptjetlimsBG[1]);}
 961     else{ fhPtJet      = new TH1F("hPtJet",  "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);}
 962    fhPtJet->Sumw2();
 963
 964    fhdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.);
 965    fhdeltaRJetTracks->Sumw2();
 966
 967    fhNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000);
 968    fhNtrArr->Sumw2();
 969
 970    fhNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);
 971    fhNJetPerEv->Sumw2();
 972
 973    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
 974    //NB at the moment fillet with jet track imp par!!!
 975    fhVtxX = new TH1F("hVtxX", "Vertex X;vtx_x",100,-0.5,0.5);
 976    fhVtxY = new TH1F("hVtxY", "Vertex Y;vtx_y",100,-0.5,0.5);
 977    fhVtxZ = new TH1F("hVtxZ", "Vertex Z;vtx_z",100,-10,10);
 978    /*
 979    if(fUseMCInfo){
 980    //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
 981       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
 982       fOutput->Add(fhImpParB);
 983
 984    }
 985    */
 986
 987    fOutput->Add(fhEjetTrks);
 988    fOutput->Add(fhPhiJetTrks);
 989    fOutput->Add(fhEtaJetTrks);
 990    fOutput->Add(fhPtJetTrks);
 991    fOutput->Add(fhEjet);
 992    fOutput->Add(fhPhiJet);
 993    fOutput->Add(fhEtaJet);
 994    fOutput->Add(fhPtJet);
 995    fOutput->Add(fhdeltaRJetTracks);
 996    fOutput->Add(fhNtrArr);
 997    fOutput->Add(fhNJetPerEv);
 998    fOutput->Add(fhImpPar);
 999    fOutput->Add(fhVtxX);
1000    fOutput->Add(fhVtxY);
1001    fOutput->Add(fhVtxZ);
1002    if(fJetOnlyMode && fSwitchOnSparses){
1003       //thnsparse for jets
1004       const Int_t nAxis=6;
1005        if (!fJetCont->GetRhoName().IsNull()) {
1006            const Int_t nbinsSparse[nAxis]={nbinsSpsphi,nbinsSpseta, nbinsptjetBG, nbinsSpsptjet,nbinsSpsContrib, nbinsSpsA};
1007            const Double_t minSparse[nAxis]={philims[0],etalims[0],ptjetlimsBG[0],ptjetlims[0],static_cast<Double_t>(nContriblims[0]),arealims[0]};
1008            const Double_t maxSparse[nAxis]={philims[1],etalims[1],ptjetlimsBG[1],ptjetlims[1],static_cast<Double_t>(nContriblims[1]),arealims[1]};
1009            fhsJet=new THnSparseF("hsJet","#phi, #eta, p_{T}^{jet}, E^{jet}, N contrib, Area", nAxis, nbinsSparse, minSparse, maxSparse);
1010            fhsJet->Sumw2();
1011        }
1012        else
1013        {
1014            const Int_t nbinsSparse[nAxis]={nbinsSpsphi,nbinsSpseta, nbinsSpsptjet, nbinsSpsptjet,nbinsSpsContrib, nbinsSpsA};
1015            const Double_t minSparse[nAxis]={philims[0],etalims[0],ptjetlims[0],ptjetlims[0],static_cast<Double_t>(nContriblims[0]),arealims[0]};
1016            const Double_t maxSparse[nAxis]={philims[1],etalims[1],ptjetlims[1],ptjetlims[1],static_cast<Double_t>(nContriblims[1]),arealims[1]};
1017            fhsJet=new THnSparseF("hsJet","#phi, #eta, p_{T}^{jet}, E^{jet}, N contrib, Area", nAxis, nbinsSparse, minSparse, maxSparse);
1018            fhsJet->Sumw2();
1019        }
1020
1021       fOutput->Add(fhsJet);
1022
1023    }
1024
1025    if(!fJetOnlyMode){
1026
1027       //debugging histograms
1028       fhControlDInJ=new TH1I("hControlDInJ","Checks D in Jet",8, -0.5,7.5);
1029       fhControlDInJ->GetXaxis()->SetBinLabel(1,"DR In,1 daugh out");
1030       fhControlDInJ->GetXaxis()->SetBinLabel(2,"DR In,2 daugh out");
1031       fhControlDInJ->GetXaxis()->SetBinLabel(3,"DR In,3 daugh out");
1032       fhControlDInJ->GetXaxis()->SetBinLabel(4,"Tot tracks non matched");
1033       fhControlDInJ->GetXaxis()->SetBinLabel(5,"D0 daug missing");
1034       fhControlDInJ->GetXaxis()->SetBinLabel(6,"soft pi missing");
1035       fhControlDInJ->GetXaxis()->SetBinLabel(7,"IDprong diff GetDau");
1036       fhControlDInJ->GetXaxis()->SetBinLabel(8,"still z>1 (cand)");
1037
1038       fhControlDInJ->SetNdivisions(1);
1039       fhControlDInJ->GetXaxis()->SetLabelSize(0.05);
1040       fOutput->Add(fhControlDInJ);
1041
1042       fhmissingp=new TH1F("hmissingp", "Distribution of missing momentum (modulo);|p_{missing}|", 200, 0.,20.);
1043       fOutput->Add(fhmissingp);
1044
1045       fhMissPi=new TH1F*[3];
1046       for(Int_t k=0;k<3;k++) {
1047          fhMissPi[k]=new TH1F(Form("hMissP%d",k), Form("Missing p component {%d};p_{%d}",k,k), 400, -10.,10.);
1048          fOutput->Add(fhMissPi[k]);
1049       }
1050       fhDeltaPtJet=new TH1F("hDeltaPtJet", "Difference between the jet pt before and after correction;p_{T}^{jet,recalc}-p_{T}^{jet,orig}", 200, 0.,20.);
1051
1052       fOutput->Add(fhDeltaPtJet);
1053       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.);
1054       fOutput->Add(fhRelDeltaPtJet);
1055
1056       fhzDinjet=new TH1F("hzDinjet","Z of candidates with daughters in jet;z",nbinsz,zlims[0],zlims[1]);
1057       fOutput->Add(fhzDinjet);
1058       //frag func of tracks in the jet
1059       fhztracksinjet=new TH1F("hztracksinjet","Z of tracks in jet;z",nbinsz,zlims[0],zlims[1]);
1060       fOutput->Add(fhztracksinjet);
1061
1062       //check jet and tracks
1063       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);
1064       fOutput->Add(fhDiffPtTrPtJzok);
1065       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);
1066       fOutput->Add(fhDiffPtTrPtJzNok);
1067       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);
1068       fOutput->Add(fhDiffPzTrPtJzok);
1069       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);
1070       fOutput->Add(fhDiffPzTrPtJzNok);
1071       fhNtrkjzNok=new TH1I("hNtrkjzNok", "Number of tracks in a jet with z>1;N^{tracks} (z>1)",5,0,5);
1072       fOutput->Add(fhNtrkjzNok);
1073
1074       //calculate frag func with pt (simply ptD(or track)\cdot pt jet /ptjet^2)
1075       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]);
1076       fOutput->Add(fhzDT);
1077       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]);
1078       fOutput->Add(fhztracksinjetT);
1079
1080       fhIDddaugh=new TH1I("hIDddaugh", "ID of daughters;ID", 2001,-1000,1000);
1081       fOutput->Add(fhIDddaugh);
1082       fhIDddaughOut=new TH1I("hIDddaughOut", "ID of daughters not found in jet;ID", 2001,-1000,1000);
1083       fOutput->Add(fhIDddaughOut);
1084       fhIDjetTracks=new TH1I("hIDjetTracks", "ID of jet tracks;ID", 2001,-1000,1000);
1085       fOutput->Add(fhIDjetTracks);
1086
1087       fhDRdaughOut=new TH1F("hDRdaughOut", "#Delta R of daughters not belonging to the jet tracks (D in jet);#Delta R",200, 0.,2.);
1088       fOutput->Add(fhDRdaughOut);
1089
1090
1091       if(fCandidateType==kDstartoKpipi)
1092       {
1093          if(fSwitchOnSB){
1094             fhDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
1095             fhDiffSideBand->SetStats(kTRUE);
1096             fhDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
1097             fhDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
1098             fhDiffSideBand->Sumw2();
1099             fOutput->Add(fhDiffSideBand);
1100          }
1101
1102          fhPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
1103          fhPtPion->SetStats(kTRUE);
1104          fhPtPion->GetXaxis()->SetTitle("GeV/c");
1105          fhPtPion->GetYaxis()->SetTitle("Entries");
1106          fhPtPion->Sumw2();
1107          fOutput->Add(fhPtPion);
1108
1109       }
1110       // D related histograms
1111       fhNDPerEvNoJet=new TH1F("hNDPerEvNoJet","Number of candidates per event with no jets; N candidate/ev with no jet", 20, 0., 20.);
1112       fhNDPerEvNoJet->Sumw2();
1113       fOutput->Add(fhNDPerEvNoJet);
1114
1115       fhptDPerEvNoJet=new TH1F("hptDPerEvNoJet","pt distribution of candidates per events with no jets; p_{t}^{D} (GeV/c)",nbinsptD, ptDlims[0],ptDlims[1]);
1116       fhptDPerEvNoJet->Sumw2();
1117       fOutput->Add(fhptDPerEvNoJet);
1118
1119       if(fSwitchOnSparses){
1120          const Int_t    nAxisD=4;
1121          const Int_t    nbinsSparseD[nAxisD]={nbinsSpseta,nbinsSpsphi,nbinsSpsptD,nbinsSpsmass};
1122          const Double_t minSparseD[nAxisD]  ={etalims[0],philims[0],ptDlims[0],fMinMass};
1123          const Double_t maxSparseD[nAxisD]  ={etalims[1],philims[1],ptDlims[1],fMaxMass};
1124          fhsDstandalone=new THnSparseF("hsDstandalone","#phi, #eta, p_{T}^{D}, and mass", nAxisD, nbinsSparseD, minSparseD, maxSparseD);
1125          fhsDstandalone->Sumw2();
1126
1127          fOutput->Add(fhsDstandalone);
1128       }
1129
1130       /*
1131       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]);
1132       hPtJetWithD->Sumw2();
1133       //for the MC this histogram is filled with the real background
1134       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]);
1135       hPtJetWithDsb->Sumw2();
1136
1137       fOutput->Add(hPtJetWithD);
1138       fOutput->Add(hPtJetWithDsb);
1139
1140       */
1141       fhNJetPerEvNoD=new TH1F("hNJetPerEvNoD","Number of jets per event with no D; number of jets/ev with no D",10,-0.5,9.5);
1142       fhNJetPerEvNoD->Sumw2();
1143
1144       fhPtJetPerEvNoD=new TH1F("hPtJetPerEvNoD","pt distribution of jets per event with no D; #it{p}_{T}^{jet} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
1145       fhPtJetPerEvNoD->Sumw2();
1146
1147       fOutput->Add(fhNJetPerEvNoD);
1148       fOutput->Add(fhPtJetPerEvNoD);
1149
1150       fhDeltaRD=new TH1F("hDeltaRD",Form("#Delta R distribution of D candidates %s selected;#Delta R", fUseMCInfo ? "(S+B)" : ""),200, 0.,10.);
1151       fhDeltaRD->Sumw2();
1152       fOutput->Add(fhDeltaRD);
1153
1154       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]);
1155       fhDeltaRptDptj->Sumw2();
1156       fOutput->Add(fhDeltaRptDptj);
1157
1158       if(fUseMCInfo){
1159          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]);
1160          fhDeltaRptDptjB->Sumw2();
1161          fOutput->Add(fhDeltaRptDptjB);
1162       }
1163
1164       //background (side bands for the Dstar and like sign for D0)
1165       AliJetContainer *jetCont=GetJetContainer(0);
1166       if(!jetCont){
1167          Printf("Container 0 not found, try with name %s", fJetArrName.Data());
1168          jetCont=GetJetContainer(fJetArrName);
1169          if(!jetCont) Printf("NOT FOUND AGAIN");
1170          return kFALSE;
1171       }
1172       Printf("CONTAINER NAME IS %s", jetCont->GetArrayName().Data());
1173       //fJetRadius=jetCont->GetJetRadius();
1174       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]);
1175       fhInvMassptD->SetStats(kTRUE);
1176       fhInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
1177       fhInvMassptD->GetYaxis()->SetTitle("#it{p}_{t}^{D} (GeV/c)");
1178       fhInvMassptD->Sumw2();
1179
1180       fOutput->Add(fhInvMassptD);
1181
1182       if(fUseMCInfo){
1183          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]);
1184          fhInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
1185          fhInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
1186          fhInvMassptDbg->Sumw2();
1187          fOutput->Add(fhInvMassptDbg);
1188
1189       }
1190
1191       if(fSwitchOnSparses){
1192          Double_t pi=TMath::Pi(), philims2[2];
1193          philims2[0]=-pi*1./2.;
1194          philims2[1]=pi*3./2.;
1195          fhsDphiz=0x0; //definition below according to the switches
1196
1197          if(fSwitchOnSB && fSwitchOnPhiAxis && fSwitchOnOutOfConeAxis){
1198             AliInfo("Creating a 9 axes container: This might cause large memory usage");
1199             const Int_t nAxis=9;
1200              if(!fJetCont->GetRhoName().IsNull())
1201              {
1202                  const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsphi,nbinsptjetBG,nbinsSpsptD,nbinsSpsmass,2, 2, 2, 2};
1203                  const Double_t minSparse[nAxis]={zlims[0],philims2[0],ptjetlimsBG[0],ptDlims[0],fMinMass,-0.5, -0.5,-0.5,-0.5};
1204                  const Double_t maxSparse[nAxis]={zlims[1],philims2[1],ptjetlimsBG[1],ptDlims[1],fMaxMass, 1.5, 1.5, 1.5 , 1.5};
1205              }
1206              else
1207              {
1208                  const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsphi,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,2, 2, 2, 2};
1209                  const Double_t minSparse[nAxis]={zlims[0],philims2[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5, -0.5,-0.5,-0.5};
1210                  const Double_t maxSparse[nAxis]={zlims[1],philims2[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5, 1.5 , 1.5};
1211              }
1212             fNAxesBigSparse=nAxis;
1213
1214             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);
1215          }
1216
1217          if(!fSwitchOnPhiAxis || !fSwitchOnOutOfConeAxis || !fSwitchOnSB){
1218             fSwitchOnPhiAxis=0;
1219             fSwitchOnOutOfConeAxis=0;
1220             fSwitchOnSB=0;
1221             if(fUseMCInfo){
1222                AliInfo("Creating a 7 axes container (MB background candidates)");
1223                const Int_t nAxis=7;
1224                 if (!fJetCont->GetRhoName().IsNull()) {
1225                     const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsptjetBG,nbinsSpsptD,nbinsSpsmass,2, 2, 2};
1226                     const Double_t minSparse[nAxis]={zlims[0],ptjetlimsBG[0],ptDlims[0],fMinMass, -0.5,-0.5,-0.5};
1227                     const Double_t maxSparse[nAxis]={zlims[1],ptjetlimsBG[1],ptDlims[1],fMaxMass, 1.5, 1.5 , 1.5};
1228                 }
1229                 else
1230                 {
1231                     const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,2, 2, 2};
1232                     const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass, -0.5,-0.5,-0.5};
1233                     const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5 , 1.5};
1234                 }
1235                fNAxesBigSparse=nAxis;
1236                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);
1237
1238             } else{
1239                AliInfo("Creating a 6 axes container");
1240                const Int_t nAxis=6;
1241                 if (!fJetCont->GetRhoName().IsNull()) {
1242                     const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsptjetBG,nbinsSpsptD,nbinsSpsmass, 2, 2};
1243                     const Double_t minSparse[nAxis]={zlims[0],ptjetlimsBG[0],ptDlims[0],fMinMass,-0.5,-0.5};
1244                     const Double_t maxSparse[nAxis]={zlims[1],ptjetlimsBG[1],ptDlims[1],fMaxMass, 1.5, 1.5};
1245                 }
1246                 else
1247                 {
1248                     const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass, 2, 2};
1249                     const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5,-0.5};
1250                     const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5};
1251                 }
1252                fNAxesBigSparse=nAxis;
1253
1254                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);
1255             }
1256          }
1257          if(!fhsDphiz) AliFatal("No THnSparse created");
1258          fhsDphiz->Sumw2();
1259
1260          fOutput->Add(fhsDphiz);
1261       }
1262    }
1263    PostData(1, fOutput);
1264
1265    return kTRUE;
1266 }
1267
1268 //_______________________________________________________________________________
1269
1270 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet,  AliAODEvent* aodEvent){
1271
1272    Double_t ptD=candidate->Pt();
1273    Double_t deltaR=DeltaR(jet,candidate);
1274    Double_t phiD=candidate->Phi();
1275    Double_t deltaphi = jet->Phi()-phiD;
1276    if(deltaphi<=-(TMath::Pi())/2.) deltaphi = deltaphi+2.*(TMath::Pi());
1277    if(deltaphi>(3.*(TMath::Pi()))/2.) deltaphi = deltaphi-2.*(TMath::Pi());
1278    Double_t z=Z(candidate,jet);
1279    /*
1280    if(z>1) {
1281       ((TH1I*)fOutput->FindObject("hControlDInJ"))->Fill(7);
1282       Double_t pmissing=TMath::Sqrt(fPmissing[0]*fPmissing[0]+fPmissing[1]*fPmissing[1]+fPmissing[2]*fPmissing[2]);
1283
1284       for(Int_t k=0;k<3;k++) ((TH1F*)fOutput->FindObject(Form("hMissP%d",k)))->Fill(fPmissing[k]);
1285
1286       ((TH1F*)fOutput->FindObject("hmissingp"))->Fill(pmissing);
1287       Double_t ptdiff=fPtJet-jet->Pt();
1288       ((TH1F*)fOutput->FindObject("hDeltaPtJet"))->Fill(ptdiff);
1289       ((TH1F*)fOutput->FindObject("hRelDeltaPtJet"))->Fill(ptdiff/jet->Pt());
1290
1291
1292    }
1293    */
1294    if(fIsDInJet)fhzDT->Fill(Z(candidate,jet,kTRUE));
1295
1296    fhDeltaRD->Fill(deltaR);
1297    fhDeltaRptDptj->Fill(deltaR,ptD,fPtJet);
1298
1299    Bool_t bDInEMCalAcc=InEMCalAcceptance(candidate);
1300    Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1301    if(fUseReco){
1302       if(fCandidateType==kD0toKpi) {
1303          AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;
1304
1305          FillHistogramsD0JetCorr(dzero,deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc, aodEvent);
1306
1307       }
1308
1309       if(fCandidateType==kDstartoKpipi) {
1310          AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;
1311          FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc);
1312
1313       }
1314    } else{ //generated level
1315       //AliInfo("Non reco");
1316       FillHistogramsMCGenDJetCorr(deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc);
1317    }
1318 }
1319
1320 //_______________________________________________________________________________
1321
1322 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){
1323
1324
1325    Double_t masses[2]={0.,0.};
1326    Int_t pdgdaughtersD0[2]={211,321};//pi,K
1327    Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1328
1329    masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1330    masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1331
1332    //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1333    Double_t *point=0x0;
1334    if(fNAxesBigSparse==9){
1335       point=new Double_t[9];
1336       point[0]=z;
1337       point[1]=dPhi;
1338       point[2]=ptj;
1339       point[3]=ptD;
1340       point[4]=masses[0];
1341       point[5]=0;
1342       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1343       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1344       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1345    }
1346    if(fNAxesBigSparse==6){
1347       point=new Double_t[6];
1348       point[0]=z;
1349       point[1]=ptj;
1350       point[2]=ptD;
1351       point[3]=masses[0];
1352       point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1353       point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1354 }
1355   if(fNAxesBigSparse==7){
1356       point=new Double_t[7];
1357       point[0]=z;
1358       point[1]=ptj;
1359       point[2]=ptD;
1360       point[3]=masses[0];
1361       point[4]=0;
1362       point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1363       point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1364
1365    }
1366
1367
1368    //Printf("Candidate in FillHistogramsD0JetCorr IsA %s", (candidate->IsA())->GetName());
1369    Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
1370    if(isselected==1 || isselected==3) {
1371
1372       //if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[0],ptD);
1373
1374       FillMassHistograms(masses[0], ptD);
1375       if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1376    }
1377    if(isselected>=2) {
1378       //if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[1],ptD);
1379
1380       FillMassHistograms(masses[1], ptD);
1381       if(fNAxesBigSparse==9) point[4]=masses[1];
1382       if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[1];
1383       if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1384    }
1385
1386 }
1387
1388 //_______________________________________________________________________________
1389
1390 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi,  Double_t z, Double_t ptD, Double_t ptj, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc){
1391   //dPhi and z not used at the moment,but will be (re)added
1392
1393    AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
1394    Double_t deltamass= dstar->DeltaInvMass();
1395    //Double_t massD0= dstar->InvMassD0();
1396
1397
1398    fhPtPion->Fill(softpi->Pt());
1399
1400    //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1401    Int_t isSB=0;//IsDzeroSideBand(dstar);
1402
1403    //Double_t point[]={z,dPhi,ptj,ptD,deltamass,static_cast<Double_t>(isSB), static_cast<Double_t>(fIsDInJet ? 1 : 0),bDInEMCalAcc,bJetInEMCalAcc};
1404    Double_t *point=0x0;
1405    if(fNAxesBigSparse==9){
1406       point=new Double_t[9];
1407       point[0]=z;
1408       point[1]=dPhi;
1409       point[2]=ptj;
1410       point[3]=ptD;
1411       point[4]=deltamass;
1412       point[5]=static_cast<Double_t>(isSB);
1413       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1414       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1415       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1416    }
1417    if(fNAxesBigSparse==6){
1418       point=new Double_t[6];
1419       point[0]=z;
1420       point[1]=ptj;
1421       point[2]=ptD;
1422       point[3]=deltamass;
1423       point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1424       point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1425    }
1426    if(fNAxesBigSparse==7){
1427       point=new Double_t[7];
1428       point[0]=z;
1429       point[1]=ptj;
1430       point[2]=ptD;
1431       point[3]=deltamass;
1432       point[4]=0;
1433       point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1434       point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1435    }
1436
1437    //if(fIsDInJet) hPtJetWithD->Fill(ptj,deltamass,ptD);
1438
1439    FillMassHistograms(deltamass, ptD);
1440    if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1441
1442 }
1443
1444 //_______________________________________________________________________________
1445
1446 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t z,Double_t ptD,Double_t ptjet, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc){
1447
1448    Double_t pdgmass=0;
1449    if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1450    if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
1451    //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1452
1453    //Double_t point[9]={z,dPhi,ptjet,ptD,pdgmass,0, static_cast<Double_t>(fIsDInJet ? 1 : 0),bDInEMCalAcc,bJetInEMCalAcc};
1454
1455    Double_t *point=0x0;
1456    if(fNAxesBigSparse==9){
1457       point=new Double_t[9];
1458       point[0]=z;
1459       point[1]=dPhi;
1460       point[2]=ptjet;
1461       point[3]=ptD;
1462       point[4]=pdgmass;
1463       point[5]=0;
1464       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1465       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1466       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1467    }
1468    if(fNAxesBigSparse==6){
1469       point=new Double_t[6];
1470       point[0]=z;
1471       point[1]=ptjet;
1472       point[2]=ptD;
1473       point[3]=pdgmass;
1474       point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1475       point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1476    }
1477       if(fNAxesBigSparse==7){
1478       point=new Double_t[7];
1479       point[0]=z;
1480       point[1]=ptjet;
1481       point[2]=ptD;
1482       point[3]=pdgmass;
1483       point[4]=1;
1484       point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1485       point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1486    }
1487
1488
1489
1490    if(fNAxesBigSparse==9) point[4]=pdgmass;
1491    if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=pdgmass;
1492    if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1493    //if(fIsDInJet) {
1494    //  hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet
1495    //}
1496
1497 }
1498
1499 //_______________________________________________________________________________
1500
1501 void AliAnalysisTaskFlavourJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD){
1502
1503    if(fIsDInJet) {
1504       fhInvMassptD->Fill(mass,ptD);
1505    }
1506 }
1507
1508 //________________________________________________________________________________
1509
1510 void AliAnalysisTaskFlavourJetCorrelations::FlagFlavour(AliEmcalJet *jet){
1511
1512    AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
1513    if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;
1514    if (fIsDInJet) jet->AddFlavourTag(tag);
1515
1516    return;
1517
1518 }
1519
1520 //_______________________________________________________________________________
1521
1522 void AliAnalysisTaskFlavourJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){
1523
1524    //  D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
1525    // (expected detector resolution) on the left and right frm the D0 mass. Each band
1526    //  has a width of ~5 sigmas. Two band needed  for opening angle considerations
1527
1528    //TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
1529
1530    Bool_t bDInEMCalAcc=InEMCalAcceptance(candDstar);
1531    Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1532
1533    Double_t deltaM=candDstar->DeltaInvMass();
1534    //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);
1535    Double_t z=Z(candDstar,jet);
1536    Double_t ptD=candDstar->Pt();
1537
1538    Double_t dPhi=jet->Phi()-candDstar->Phi();
1539
1540    if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1541    if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1542
1543    Int_t isSideBand=1;//IsDzeroSideBand(candDstar);
1544    //if no SB analysis we should not enter here, so no need of checking the number of axes
1545    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)};
1546    //fill the background histogram with the side bands or when looking at MC with the real background
1547    if(isSideBand==1){
1548       fhDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background
1549       //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);
1550       if(fSwitchOnSparses) fhsDphiz->Fill(point,1.);
1551       //if(fIsDInJet){  // evaluate in the near side
1552       //         hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
1553       //}
1554    }
1555 }
1556
1557 //_______________________________________________________________________________
1558
1559 void AliAnalysisTaskFlavourJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg,AliEmcalJet* jet){
1560
1561    //need mass, deltaR, pt jet, ptD
1562    //two cases: D0 and Dstar
1563
1564    Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());
1565    if(!isselected) return;
1566
1567    Double_t ptD=candbg->Pt();
1568    Double_t deltaR=DeltaR(jet,candbg);
1569    Double_t phiD=candbg->Phi();
1570    Double_t deltaphi = jet->Phi()-phiD;
1571    if(deltaphi<=-(TMath::Pi())/2.) deltaphi = deltaphi+2.*(TMath::Pi());
1572    if(deltaphi>(3.*(TMath::Pi()))/2.) deltaphi = deltaphi-2.*(TMath::Pi());
1573    Double_t z=Z(candbg,jet);
1574
1575    if(fIsDInJet) fhzDT->Fill(Z(candbg,jet,kTRUE));
1576
1577
1578
1579
1580    fhDeltaRD->Fill(deltaR);
1581    fhDeltaRptDptjB->Fill(deltaR,ptD,fPtJet);
1582
1583    Bool_t bDInEMCalAcc=InEMCalAcceptance(candbg);
1584    Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1585
1586    //TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
1587
1588    Double_t *point=0x0;
1589    if(fNAxesBigSparse==9){
1590       point=new Double_t[9];
1591       point[0]=z;
1592       point[1]=deltaphi;
1593       point[2]=fPtJet;
1594       point[3]=ptD;
1595       point[4]=-999; //set below
1596       point[5]=1;
1597       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1598       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1599       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1600    }
1601
1602    if(fNAxesBigSparse==7){
1603       point=new Double_t[7];
1604       point[0]=z;
1605       point[1]=fPtJet;
1606       point[2]=ptD;
1607       point[3]=-999; //set below
1608       point[4]=1;
1609       point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1610       point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1611    }
1612
1613    if(fCandidateType==kDstartoKpipi){
1614       AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;
1615       Double_t deltaM=dstarbg->DeltaInvMass();
1616       fhInvMassptDbg->Fill(deltaM,ptD);
1617       //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
1618       if(fNAxesBigSparse==9) point[4]=deltaM;
1619       if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=deltaM;
1620       if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1621    }
1622
1623    if(fCandidateType==kD0toKpi){
1624       Double_t masses[2]={0.,0.};
1625       Int_t pdgdaughtersD0[2]={211,321};//pi,K
1626       Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1627
1628       masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1629       masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1630
1631
1632       if(isselected==1 || isselected==3) {
1633          //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[0],ptD);
1634          fhInvMassptDbg->Fill(masses[0],ptD);
1635          if(fNAxesBigSparse==9) point[4]=masses[0];
1636          if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[0];
1637          if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1638      }
1639       if(isselected>=2) {
1640          //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[1],ptD);
1641          fhInvMassptDbg->Fill(masses[1],ptD);
1642          if(fNAxesBigSparse==9) point[4]=masses[1];
1643          if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[1];
1644          if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1645
1646       }
1647
1648
1649    }
1650 }
1651
1652 //_______________________________________________________________________________
1653
1654 Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliEmcalJet *p1, AliVParticle *p2) const {
1655    //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1656    //It recalculates the eta-phi values if it was asked for background subtraction of the jets
1657    if(!p1 || !p2) return -1;
1658
1659     Double_t phi1=p1->Phi(),eta1=p1->Eta();
1660
1661     //It subtracts the backgroud of jets if it was asked for it.
1662     if (!fJetCont->GetRhoName().IsNull()) {
1663         AliRhoParameter *rho = 0;
1664         rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fJetCont->GetRhoName()));
1665         if(rho){
1666             Double_t pj[3];
1667             Bool_t okpj=p1->PxPyPz(pj);
1668             if(!okpj){
1669                 printf("Problems getting momenta\n");
1670                 return -999;
1671             }
1672             Double_t rhoval = rho->GetVal();
1673             pj[0] = p1->Px() - p1->Area()*(rhoval*TMath::Cos(p1->AreaPhi()));
1674             pj[1] = p1->Py() - p1->Area()*(rhoval*TMath::Sin(p1->AreaPhi()));
1675             pj[2] = p1->Pz() - p1->Area()*(rhoval*TMath::SinH(p1->AreaEta()));
1676             //Image of the function Arccos(px/pt) where pt = sqrt(px*px+py*py) is:
1677             //[0,pi]    if py > 0 and
1678             //[pi,2pi]  if py < 0
1679             phi1 = TMath::ACos(pj[0]/TMath::Sqrt(pj[0]*pj[0]+pj[1]*pj[1]));
1680             if(pj[1]<0) phi1 = 2*TMath::Pi() - phi1;
1681             eta1 = TMath::ASinH(pj[2]/TMath::Sqrt(pj[0]*pj[0]+pj[1]*pj[1]));
1682         }
1683     }
1684
1685    Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
1686
1687    Double_t dPhi=phi1-phi2;
1688    if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1689    if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1690
1691    Double_t dEta=eta1-eta2;
1692    Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
1693    return deltaR;
1694
1695 }
1696
1697 //_______________________________________________________________________________
1698
1699 Int_t AliAnalysisTaskFlavourJetCorrelations::IsDzeroSideBand(AliAODRecoCascadeHF *candDstar){
1700
1701    Double_t ptD=candDstar->Pt();
1702    Int_t bin = fCuts->PtBin(ptD);
1703    if (bin < 0){
1704       // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
1705       bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?
1706       return -1;  // out of bounds
1707    }
1708
1709    Double_t invM=candDstar->InvMassD0();
1710    Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1711
1712    Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
1713    Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
1714
1715    if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
1716    else return 0;
1717
1718 }
1719
1720 //_______________________________________________________________________________
1721
1722 Bool_t AliAnalysisTaskFlavourJetCorrelations::AreDaughtersInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Int_t*& daughOutOfJetID, AliAODTrack**& dtrks, Bool_t fillH){
1723    //returns true/false and the array of particles not found among jet constituents
1724
1725
1726    Int_t daughtersID[3];
1727    Int_t ndaugh=0;
1728    Int_t dmatchedID[3]={0,0,0};
1729    Int_t countmatches=0;
1730    if (fCandidateType==kDstartoKpipi){
1731       AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)charm;
1732       AliAODRecoDecayHF2Prong* dzero=dstar->Get2Prong();
1733       daughtersID[0]=dzero->GetProngID(0);
1734       daughtersID[1]=dzero->GetProngID(1);
1735       daughtersID[2]=dstar->GetBachelor()->GetID();
1736       ndaugh=3;
1737
1738       dtrks=new AliAODTrack*[3];
1739       dtrks[0]=(AliAODTrack*)dzero->GetDaughter(0);
1740       dtrks[1]=(AliAODTrack*)dzero->GetDaughter(1);
1741       dtrks[2]=(AliAODTrack*)dstar->GetBachelor();
1742
1743       //check
1744       if(fillH){
1745          if(daughtersID[0]!=((AliAODTrack*)dtrks[0])->GetID() || daughtersID[1]!=((AliAODTrack*)dtrks[1])->GetID())  fhControlDInJ->Fill(6);
1746
1747          fhIDddaugh->Fill(daughtersID[0]);
1748          fhIDddaugh->Fill(daughtersID[1]);
1749          fhIDddaugh->Fill(daughtersID[2]);
1750
1751       }
1752       //Printf("ID daughters %d, %d, %d",daughtersID[0], daughtersID[1], daughtersID[2]);
1753    }
1754
1755    if (fCandidateType==kD0toKpi){
1756       daughtersID[0]=charm->GetProngID(0);
1757       daughtersID[1]=charm->GetProngID(1);
1758       ndaugh=2;
1759       if(fillH){
1760          fhIDddaugh->Fill(daughtersID[0]);
1761          fhIDddaugh->Fill(daughtersID[1]);
1762       }
1763       dtrks=new AliAODTrack*[2];
1764       dtrks[0]=(AliAODTrack*)charm->GetDaughter(0);
1765       dtrks[1]=(AliAODTrack*)charm->GetDaughter(1);
1766
1767    }
1768
1769    const Int_t cndaugh=ndaugh;
1770    daughOutOfJetID=new Int_t[cndaugh];
1771
1772    Int_t nchtrk=thejet->GetNumberOfTracks();
1773    for(Int_t itk=0;itk<nchtrk;itk++){
1774       AliVTrack* tkjet=((AliPicoTrack*)thejet->TrackAt(itk,fTrackArr))->GetTrack();
1775       if(!tkjet) continue;
1776       Int_t idtkjet=tkjet->GetID();
1777       if(fillH) fhIDjetTracks->Fill(idtkjet);
1778       for(Int_t id=0;id<ndaugh;id++){
1779          if(idtkjet==daughtersID[id]) {
1780             dmatchedID[id]=daughtersID[id]; //daughter which matches a track in the jet
1781             countmatches++;
1782
1783          }
1784
1785          if(countmatches==ndaugh) break;
1786       }
1787    }
1788    //reverse: include in the array the ID of daughters not matching
1789
1790    for(Int_t id=0;id<ndaugh;id++){
1791       if(dmatchedID[id]==0) {
1792          daughOutOfJetID[id]=daughtersID[id];
1793          if(fillH) fhIDddaughOut->Fill(daughOutOfJetID[id]);
1794       }
1795       else daughOutOfJetID[id]=0;
1796    }
1797    if(fillH){
1798       if((ndaugh-countmatches) == 1) fhControlDInJ->Fill(0);
1799       if((ndaugh-countmatches) == 2) fhControlDInJ->Fill(1);
1800       if((ndaugh-countmatches) == 3) fhControlDInJ->Fill(2);
1801    }
1802    if(ndaugh!=countmatches){
1803       return kFALSE;
1804    }
1805
1806    return kTRUE;
1807 }
1808
1809 //_______________________________________________________________________________
1810
1811 Bool_t AliAnalysisTaskFlavourJetCorrelations::IsDInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Bool_t fillH){
1812
1813    //check the conditions type:
1814    //type 0 : DeltaR < jet radius, don't check daughters (and don't correct pt jet)
1815    //type 1 : DeltaR < jet radius and check for all daughters among jet tracks, don't correct ptjet
1816    //type 2 (default) : DeltaR < jet radius and check for all daughters among jet tracks, if not present, correct the ptjet
1817    //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
1818
1819    fPmissing[0]=0;
1820    fPmissing[1]=0;
1821    fPmissing[2]=0;
1822
1823    Bool_t testDeltaR=kFALSE;
1824    if(DeltaR(thejet,charm)<fJetRadius) testDeltaR=kTRUE;
1825    //for type 3 this check should be redone after the modification of the jet axis
1826
1827    Int_t* daughOutOfJet;
1828    AliAODTrack** charmDaugh;
1829    Bool_t testDaugh=AreDaughtersInJet(thejet, charm, daughOutOfJet,charmDaugh,fillH);
1830    if(testDaugh && testDeltaR) {
1831       //Z of candidates with daughters in the jet
1832       fhzDinjet->Fill(Z(charm,thejet));
1833
1834    }
1835
1836    TVector3 thejetv;    //(x=0,y=0,z=0)
1837    thejetv.SetPtEtaPhi(thejet->Pt(),thejet->Eta(),thejet->Phi());
1838    TVector3 newjet(thejetv);
1839    TVector3 correction; //(x=0,y=0,z=0)
1840    TVector3 charmcand;  //(x=0,y=0,z=0)
1841    charmcand.SetPtEtaPhi(charm->Pt(),charm->Eta(),charm->Phi());
1842
1843    if(!testDaugh && testDeltaR && fTypeDInJet>=2){
1844
1845       Int_t ndaugh=3;
1846       if(fCandidateType==kD0toKpi) ndaugh=2;
1847       Int_t nOut=ndaugh;
1848
1849       for(Int_t id=0;id<ndaugh;id++){
1850          if(daughOutOfJet[id]!=0){ //non-matched daughter
1851             //get track and its momentum
1852             nOut--;
1853             if(fillH) {
1854                fhControlDInJ->Fill(3);
1855                if(id<2) fhControlDInJ->Fill(4);
1856                if(id==2)fhControlDInJ->Fill(5);
1857                fhDRdaughOut->Fill(DeltaR(thejet, charmDaugh[id]));
1858             }
1859             if(fTypeDInJet==2){
1860
1861                newjet.SetX(newjet(0)+charmDaugh[id]->Px());
1862                newjet.SetY(newjet(1)+charmDaugh[id]->Py());
1863                newjet.SetZ(newjet(2)+charmDaugh[id]->Pz());
1864             }
1865             if(fTypeDInJet==3){
1866
1867                Double_t ptdaug  = charmDaugh[id]->Pt();
1868                Double_t ptjet   = newjet.Perp();
1869                Double_t ptn     = ptjet+ptdaug;
1870                Double_t phidaug = charmDaugh[id]->Phi();
1871                Double_t phijet  = newjet.Phi();
1872                Double_t phin    = (phijet*ptjet+phidaug*ptdaug)/(ptjet+ptdaug);
1873
1874                Double_t etadaug = charmDaugh[id]->Eta();
1875                Double_t etajet  = newjet.Eta();
1876                Double_t etan    = (etajet*ptjet+etadaug*ptdaug)/(ptjet+ptdaug);
1877
1878                newjet.SetPtEtaPhi(ptn,etan,phin);
1879
1880             }
1881          }
1882       }
1883
1884       correction=newjet-thejetv;
1885       fPmissing[0]=correction(0);
1886       fPmissing[1]=correction(1);
1887       fPmissing[2]=correction(2);
1888
1889       //now the D is within the jet
1890       testDaugh=kTRUE;
1891       if(fTypeDInJet==3){ //recalc R(D,jet)
1892          Double_t dRnew=newjet.DeltaR(charmcand);
1893          if(dRnew<fJetRadius) testDeltaR=kTRUE;
1894          else testDeltaR=kFALSE;
1895       }
1896    }
1897
1898    delete[] charmDaugh;
1899
1900    Bool_t result=0;
1901    switch(fTypeDInJet){
1902    case 0:
1903       result=testDeltaR;
1904       break;
1905    case 1:
1906       result=testDeltaR && testDaugh;
1907       break;
1908    case 2: //this case defines fPmissing
1909       result=testDeltaR && testDaugh;
1910       break;
1911    case 3: //this case defines fPmissing and recalculate R(jet,D) with the updated jet axis
1912       result=testDeltaR && testDaugh;
1913       break;
1914
1915    default:
1916       AliInfo("Selection type not specified, use 1");
1917       result=testDeltaR && testDaugh;
1918       break;
1919    }
1920    return result;
1921 }
1922
1923 Bool_t AliAnalysisTaskFlavourJetCorrelations::InEMCalAcceptance(AliVParticle *vpart){
1924    //check eta phi of a VParticle: return true if it is in the EMCal acceptance, false otherwise
1925
1926    Double_t phiEMCal[2]={1.405,3.135},etaEMCal[2]={-0.7,0.7};
1927    Bool_t binEMCal=kTRUE;
1928    Double_t phi=vpart->Phi(), eta=vpart->Eta();
1929    if(phi<phiEMCal[0] || phi>phiEMCal[1]) binEMCal=kFALSE;
1930    if(eta<etaEMCal[0] || eta>etaEMCal[1]) binEMCal=kFALSE;
1931    return binEMCal;
1932
1933
1934 }