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     TString sname("Rho");
 532     if (!sname.IsNull()) {
 533         rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(sname));
 534         if(rho) rhoval = rho->GetVal();
 535     }
 536
 537
 538    // we start with jets
 539    Double_t ejet   = 0;
 540    Double_t phiJet = 0;
 541    Double_t etaJet = 0;
 542    Double_t ptjet = 0;
 543    Double_t leadingJet =0;
 544    Double_t pointJ[6];
 545
 546    Int_t ntrarr=fTrackCont->GetNParticles();
 547    fhNtrArr->Fill(ntrarr);
 548
 549    for(Int_t i=0;i<ntrarr;i++){
 550       AliAODTrack *jtrack=static_cast<AliAODTrack*>(fTrackCont->GetParticle(i));
 551       //reusing histograms
 552       if(!jtrack) continue;
 553       fhPtJetTrks->Fill(jtrack->Pt());
 554       fhPhiJetTrks->Fill(jtrack->Phi());
 555       fhEtaJetTrks->Fill(jtrack->Eta());
 556       fhEjetTrks->Fill(jtrack->E());
 557       Double_t vDCAglobalxy;
 558       //Double_t vDCAglobalz;
 559       if(jtrack->IsGlobalConstrained()){
 560       // retrieve impact parameter
 561       Double_t pos[3];
 562
 563       jtrack->GetXYZ(pos);
 564
 565       vDCAglobalxy  = pos[0] - vtx[0]; //should be impact parameter in transverse plane
 566       //vDCAglobalz  = pos[1] - vtx[1]; //should be impact parameter in z direction
 567       }
 568       else{
 569       vDCAglobalxy=jtrack->DCA();
 570       //vDCAglobalz=jtrack->ZAtDCA();
 571       }
 572       fhImpPar->Fill(vDCAglobalxy,jtrack->Pt());
 573       //Printf("Track position  %.3f,%.3f,%.3f",pos[0],pos[1],pos[2]);
 574    }
 575
 576
 577    //option to use only the leading jet
 578    if(fLeadingJetOnly){
 579       for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) {
 580          AliEmcalJet* jetL = (AliEmcalJet*)fJetCont->GetJet(iJetsL);
 581          ptjet   = jetL->Pt() - jetL->Area()*rhoval; //It takes into account the background subtraction
 582          if(ptjet>leadingJet ) leadingJet = ptjet;
 583
 584       }
 585    }
 586
 587    Int_t cntjet=0;
 588    //loop over jets and consider only the leading jet in the event
 589    for (Int_t iJets = 0; iJets<njets; iJets++) {
 590       fPmissing[0]=0;
 591       fPmissing[1]=0;
 592       fPmissing[2]=0;
 593
 594       //Printf("Jet N %d",iJets);
 595       AliEmcalJet* jet = (AliEmcalJet*)fJetCont->GetJet(iJets);
 596       //Printf("Corr task Accept Jet");
 597       if(!AcceptJet(jet)) {
 598          fhstat->Fill(5);
 599          continue;
 600       }
 601
 602       //jets variables
 603       ejet   = jet->E();
 604       phiJet = jet->Phi();
 605       etaJet = jet->Eta();
 606       fPtJet = jet->Pt() - jet->Area()*rhoval; //It takes into account the background subtraction
 607       Double_t origPtJet=fPtJet;
 608
 609       pointJ[0] = phiJet;
 610       pointJ[1] = etaJet;
 611       pointJ[2] = ptjet - jet->Area()*rhoval; //It takes into account the background subtraction
 612       pointJ[3] = ejet;
 613       pointJ[4] = jet->GetNumberOfConstituents();
 614       pointJ[5] = jet->Area();
 615
 616       // choose the leading jet
 617       if(fLeadingJetOnly && (ejet<leadingJet)) continue;
 618       //used for normalization
 619       fhstat->Fill(3);
 620       cntjet++;
 621       // fill energy, eta and phi of the jet
 622       fhEjet   ->Fill(ejet);
 623       fhPhiJet ->Fill(phiJet);
 624       fhEtaJet ->Fill(etaJet);
 625       fhPtJet  ->Fill(fPtJet);
 626       if(fJetOnlyMode && fSwitchOnSparses) fhsJet->Fill(pointJ,1);
 627       //loop on jet particles
 628       Int_t ntrjet=  jet->GetNumberOfTracks();
 629       Double_t sumPtTracks=0,sumPzTracks=0;
 630       Int_t zg1jtrk=0;
 631       for(Int_t itrk=0;itrk<ntrjet;itrk++){
 632
 633          AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,fTrackArr);
 634          fhdeltaRJetTracks->Fill(DeltaR(jet,jetTrk));
 635          Double_t ztrackj=Z(jetTrk,jet);
 636          fhztracksinjet->Fill(ztrackj);
 637          sumPtTracks+=jetTrk->Pt();
 638          sumPzTracks+=jetTrk->Pz();
 639          if(ztrackj>1){
 640             zg1jtrk++;
 641          }
 642
 643          Double_t ztrackjTr=Z(jetTrk,jet,kTRUE);
 644          fhztracksinjetT->Fill(ztrackjTr);
 645
 646       }//end loop on jet tracks
 647
 648       fhNtrkjzNok->Fill(zg1jtrk);
 649       Double_t diffpt=origPtJet-sumPtTracks;
 650       Double_t diffpz=jet->Pz()-sumPzTracks;
 651       if(zg1jtrk>0){
 652          fhDiffPtTrPtJzNok->Fill(diffpt);
 653          fhDiffPzTrPtJzNok->Fill(diffpz);
 654
 655       }else{
 656          fhDiffPtTrPtJzok->Fill(diffpt);
 657          fhDiffPzTrPtJzok->Fill(diffpz);
 658       }
 659
 660       if(candidates==0){
 661
 662          fhPtJetPerEvNoD->Fill(fPtJet);
 663       }
 664       if(!fJetOnlyMode) {
 665          //Printf("N candidates %d ", candidates);
 666          for(Int_t ic = 0; ic < candidates; ic++) {
 667             fhstat->Fill(7);
 668             // D* candidates
 669             AliVParticle* charm=0x0;
 670             charm=(AliVParticle*)fCandidateArray->At(ic);
 671             if(!charm) continue;
 672             AliAODRecoDecayHF *charmdecay=(AliAODRecoDecayHF*) charm;
 673             fIsDInJet=IsDInJet(jet, charmdecay, kTRUE);
 674             if (fIsDInJet) FlagFlavour(jet);
 675             if (jet->TestFlavourTag(AliEmcalJet::kDStar) || jet->TestFlavourTag(AliEmcalJet::kD0)) fhstat->Fill(4);
 676
 677             //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.
 678
 679             Double_t pjet[3];
 680             jet->PxPyPz(pjet);
 681              //It corrects the jet momentum if it was asked for jet background subtraction
 682              pjet[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
 683              pjet[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
 684              pjet[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
 685
 686             //It corrects the jet momentum due to D daughters out of the jet
 687             RecalculateMomentum(pjet,fPmissing);
 688             fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]); //recalculated jet pt
 689
 690
 691             //debugging histograms
 692             Double_t pmissing=TMath::Sqrt(fPmissing[0]*fPmissing[0]+fPmissing[1]*fPmissing[1]+fPmissing[2]*fPmissing[2]); //recalculated jet total momentum
 693             for(Int_t k=0;k<3;k++) fhMissPi[k]->Fill(fPmissing[k]);
 694
 695             fhmissingp->Fill(pmissing);
 696             Double_t ptdiff=fPtJet-origPtJet;
 697             fhDeltaPtJet->Fill(ptdiff);
 698             fhRelDeltaPtJet->Fill(ptdiff/origPtJet);
 699
 700             FillHistogramsRecoJetCorr(charm, jet, aodEvent);
 701
 702          }//end loop on candidates
 703
 704          //retrieve side band background candidates for Dstar
 705          Int_t nsbcand = 0; if(fSideBandArray) nsbcand=fSideBandArray->GetEntriesFast();
 706
 707          for(Int_t ib=0;ib<nsbcand;ib++){
 708             if(fCandidateType==kDstartoKpipi && !fUseMCInfo){
 709                AliAODRecoCascadeHF *sbcand=(AliAODRecoCascadeHF*)fSideBandArray->At(ib);
 710                if(!sbcand) continue;
 711
 712                fIsDInJet=IsDInJet(jet, sbcand,kFALSE);
 713                Double_t pjet[3];
 714                jet->PxPyPz(pjet);
 715                 //It corrects the jet momentum if it was asked for jet background subtraction
 716                 pjet[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
 717                 pjet[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
 718                 pjet[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
 719
 720                //It corrects the jet momentum due to D daughters out of the jet
 721                RecalculateMomentum(pjet,fPmissing);
 722                fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]); //recalculated jet pt
 723
 724                SideBandBackground(sbcand,jet);
 725
 726             }
 727             if(fUseMCInfo){
 728
 729                AliAODRecoDecayHF* charmbg = 0x0;
 730                charmbg=(AliAODRecoDecayHF*)fSideBandArray->At(ib);
 731                if(!charmbg) continue;
 732                fhstat->Fill(8);
 733                fIsDInJet=IsDInJet(jet, charmbg,kFALSE);
 734                if (fIsDInJet) {
 735                   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
 736                   fhstat->Fill(9);
 737                }
 738                Double_t pjet[3];
 739                jet->PxPyPz(pjet);
 740                //It corrects the jet momentum if it was asked for jet background subtraction
 741                pjet[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
 742                pjet[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
 743                pjet[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
 744
 745                //It corrects the jet momentum due to D daughters out of the jet
 746                RecalculateMomentum(pjet,fPmissing);
 747                fPtJet=TMath::Sqrt(pjet[0]*pjet[0]+pjet[1]*pjet[1]); //recalculated jet pt
 748
 749                MCBackground(charmbg,jet);
 750             }
 751          }
 752       }
 753    } // end of jet loop
 754
 755    fhNJetPerEv->Fill(cntjet);
 756    if(candidates==0) fhNJetPerEvNoD->Fill(cntjet);
 757    PostData(1,fOutput);
 758    return kTRUE;
 759 }
 760
 761 //_______________________________________________________________________________
 762
 763 void AliAnalysisTaskFlavourJetCorrelations::Terminate(Option_t*)
 764 {
 765    // The Terminate() function is the last function to be called during
 766    // a query. It always runs on the client, it can be used to present
 767    // the results graphically or save the results to file.
 768
 769    Info("Terminate"," terminate");
 770    AliAnalysisTaskSE::Terminate();
 771
 772    fOutput = dynamic_cast<TList*> (GetOutputData(1));
 773    if (!fOutput) {
 774       printf("ERROR: fOutput not available\n");
 775       return;
 776    }
 777 }
 778
 779 //_______________________________________________________________________________
 780
 781 void  AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){
 782    Float_t mass=0;
 783    Int_t abspdg=TMath::Abs(pdg);
 784
 785    mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
 786    // compute the Delta mass for the D*
 787    if(fCandidateType==kDstartoKpipi){
 788       Float_t mass1=0;
 789       mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 790       mass = mass-mass1;
 791    }
 792
 793    fMinMass = mass-range;
 794    fMaxMass = mass+range;
 795
 796    AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
 797    if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
 798 }
 799
 800 //_______________________________________________________________________________
 801
 802 void  AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){
 803    if(uplimit>lowlimit)
 804    {
 805       fMinMass = lowlimit;
 806       fMaxMass = uplimit;
 807    }
 808    else{
 809       printf("Error! Lower limit larger than upper limit!\n");
 810       fMinMass = uplimit - uplimit*0.2;
 811       fMaxMass = uplimit;
 812    }
 813 }
 814
 815 //_______________________________________________________________________________
 816
 817 Bool_t AliAnalysisTaskFlavourJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){
 818    if(nptbins>30) {
 819       AliInfo("Maximum number of bins allowed is 30!");
 820       return kFALSE;
 821    }
 822    if(!width) return kFALSE;
 823    for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
 824    return kTRUE;
 825 }
 826
 827 //_______________________________________________________________________________
 828
 829 Double_t AliAnalysisTaskFlavourJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet, Bool_t transverse) const{
 830    if(!part || !jet){
 831       printf("Particle or jet do not exist!\n");
 832       return -99;
 833    }
 834    Double_t p[3],pj[3];
 835    Bool_t okpp=part->PxPyPz(p);
 836    Bool_t okpj=jet->PxPyPz(pj);
 837
 838     //Background Subtraction
 839     //It corrects the each component of the jet momentum for Z calculation
 840     AliRhoParameter *rho = 0;
 841     Double_t rhoval = 0;
 842     TString sname("Rho");
 843     if (!sname.IsNull()) {
 844         rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(sname));
 845         if(rho){
 846             rhoval = rho->GetVal();
 847             pj[0] = jet->Px() - jet->Area()*(rhoval*TMath::Cos(jet->AreaPhi()));
 848             pj[1] = jet->Py() - jet->Area()*(rhoval*TMath::Sin(jet->AreaPhi()));
 849             pj[2] = jet->Pz() - jet->Area()*(rhoval*TMath::SinH(jet->AreaEta()));
 850         }
 851     }
 852
 853
 854
 855    if(!okpp || !okpj){
 856       printf("Problems getting momenta\n");
 857       return -999;
 858    }
 859
 860    RecalculateMomentum(pj, fPmissing);
 861    if(transverse) return ZT(p,pj);
 862    else
 863    return Z(p,pj);
 864 }
 865
 866 //_______________________________________________________________________________
 867 Double_t AliAnalysisTaskFlavourJetCorrelations::Z(Double_t* p, Double_t *pj) const{
 868
 869    Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1]+pj[2]*pj[2];
 870    Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet2);
 871    return z;
 872 }
 873
 874
 875 //_______________________________________________________________________________
 876 Double_t AliAnalysisTaskFlavourJetCorrelations::ZT(Double_t* p, Double_t *pj) const{
 877
 878    Double_t pjet2=pj[0]*pj[0]+pj[1]*pj[1];
 879    Double_t z=(p[0]*pj[0]+p[1]*pj[1])/(pjet2);
 880    return z;
 881 }
 882
 883 //_______________________________________________________________________________
 884
 885 void AliAnalysisTaskFlavourJetCorrelations::RecalculateMomentum(Double_t* pj, const Double_t* pmissing) const {
 886
 887    pj[0]+=pmissing[0];
 888    pj[1]+=pmissing[1];
 889    pj[2]+=pmissing[2];
 890
 891 }
 892
 893 //_______________________________________________________________________________
 894
 895 Bool_t  AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){
 896
 897    // Statistics
 898    Int_t nbins=8;
 899    if(fUseMCInfo) nbins+=2;
 900
 901    fhstat=new TH1I("hstat","Statistics",nbins,-0.5,nbins-0.5);
 902    fhstat->GetXaxis()->SetBinLabel(1,"N ev anal");
 903    fhstat->GetXaxis()->SetBinLabel(2,"N ev sel");
 904    fhstat->GetXaxis()->SetBinLabel(3,"N cand sel");
 905    fhstat->GetXaxis()->SetBinLabel(4,"N jets");
 906    fhstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
 907    fhstat->GetXaxis()->SetBinLabel(6,"N jet rej");
 908    fhstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet");
 909    fhstat->GetXaxis()->SetBinLabel(8,"N jets & cand");
 910    if(fUseMCInfo) {
 911     fhstat->GetXaxis()->SetBinLabel(3,"N Signal sel & jet");
 912     fhstat->GetXaxis()->SetBinLabel(5,"N Signal in jet");
 913     fhstat->GetXaxis()->SetBinLabel(9,"N Bkg sel & jet");
 914     fhstat->GetXaxis()->SetBinLabel(10,"N Bkg in jet");
 915    }
 916    fhstat->SetNdivisions(1);
 917    fOutput->Add(fhstat);
 918
 919    const Int_t nbinsmass=300;
 920    const Int_t nbinsptjet=500;
 921    const Int_t nbinsptD=100;
 922    const Int_t nbinsz=100;
 923    const Int_t nbinsphi=200;
 924    const Int_t nbinseta=100;
 925
 926    //binning for THnSparse
 927    const Int_t nbinsSpsmass=50;
 928    const Int_t nbinsSpsptjet=100;
 929    const Int_t nbinsSpsptD=50;
 930    const Int_t nbinsSpsz=100;
 931    const Int_t nbinsSpsphi=100;
 932    const Int_t nbinsSpseta=60;
 933    const Int_t nbinsSpsContrib=100;
 934    const Int_t nbinsSpsA=100;
 935
 936    const Float_t ptjetlims[2]={0.,200.};
 937    const Float_t ptDlims[2]={0.,50.};
 938    const Float_t zlims[2]={0.,1.2};
 939    const Float_t philims[2]={0.,6.3};
 940    const Float_t etalims[2]={-1.5,1.5};
 941    const Int_t   nContriblims[2]={0,100};
 942    const Float_t arealims[2]={0.,2};
 943
 944    // jet related fistograms
 945
 946    fhEjetTrks      = new TH1F("hEjetTrks",  "Jet tracks energy distribution;Energy (GeV)",500,0,200);
 947    fhEjetTrks->Sumw2();
 948    fhPhiJetTrks    = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)",  nbinsphi,philims[0],philims[1]);
 949    fhPhiJetTrks->Sumw2();
 950    fhEtaJetTrks    = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta",  nbinseta,etalims[0],etalims[1]);
 951    fhEtaJetTrks->Sumw2();
 952    fhPtJetTrks     = new TH1F("hPtJetTrks",  "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
 953    fhPtJetTrks->Sumw2();
 954
 955    fhEjet      = new TH1F("hEjet",  "Jet energy distribution;Energy (GeV)",500,0,200);
 956    fhEjet->Sumw2();
 957    fhPhiJet    = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)",  nbinsphi,philims[0],philims[1]);
 958    fhPhiJet->Sumw2();
 959    fhEtaJet    = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
 960    fhEtaJet->Sumw2();
 961    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       const Int_t nbinsSparse[nAxis]={nbinsSpsphi,nbinsSpseta, nbinsSpsptjet, nbinsSpsptjet,nbinsSpsContrib, nbinsSpsA};
1006       const Double_t minSparse[nAxis]={philims[0],etalims[0],ptjetlims[0],ptjetlims[0],static_cast<Double_t>(nContriblims[0]),arealims[0]};
1007       const Double_t
1008         maxSparse[nAxis]={philims[1],etalims[1],ptjetlims[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       fOutput->Add(fhsJet);
1013
1014    }
1015
1016    if(!fJetOnlyMode){
1017
1018       //debugging histograms
1019       fhControlDInJ=new TH1I("hControlDInJ","Checks D in Jet",8, -0.5,7.5);
1020       fhControlDInJ->GetXaxis()->SetBinLabel(1,"DR In,1 daugh out");
1021       fhControlDInJ->GetXaxis()->SetBinLabel(2,"DR In,2 daugh out");
1022       fhControlDInJ->GetXaxis()->SetBinLabel(3,"DR In,3 daugh out");
1023       fhControlDInJ->GetXaxis()->SetBinLabel(4,"Tot tracks non matched");
1024       fhControlDInJ->GetXaxis()->SetBinLabel(5,"D0 daug missing");
1025       fhControlDInJ->GetXaxis()->SetBinLabel(6,"soft pi missing");
1026       fhControlDInJ->GetXaxis()->SetBinLabel(7,"IDprong diff GetDau");
1027       fhControlDInJ->GetXaxis()->SetBinLabel(8,"still z>1 (cand)");
1028
1029       fhControlDInJ->SetNdivisions(1);
1030       fhControlDInJ->GetXaxis()->SetLabelSize(0.05);
1031       fOutput->Add(fhControlDInJ);
1032
1033       fhmissingp=new TH1F("hmissingp", "Distribution of missing momentum (modulo);|p_{missing}|", 200, 0.,20.);
1034       fOutput->Add(fhmissingp);
1035
1036       fhMissPi=new TH1F*[3];
1037       for(Int_t k=0;k<3;k++) {
1038          fhMissPi[k]=new TH1F(Form("hMissP%d",k), Form("Missing p component {%d};p_{%d}",k,k), 400, -10.,10.);
1039          fOutput->Add(fhMissPi[k]);
1040       }
1041       fhDeltaPtJet=new TH1F("hDeltaPtJet", "Difference between the jet pt before and after correction;p_{T}^{jet,recalc}-p_{T}^{jet,orig}", 200, 0.,20.);
1042
1043       fOutput->Add(fhDeltaPtJet);
1044       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.);
1045       fOutput->Add(fhRelDeltaPtJet);
1046
1047       fhzDinjet=new TH1F("hzDinjet","Z of candidates with daughters in jet;z",nbinsz,zlims[0],zlims[1]);
1048       fOutput->Add(fhzDinjet);
1049       //frag func of tracks in the jet
1050       fhztracksinjet=new TH1F("hztracksinjet","Z of tracks in jet;z",nbinsz,zlims[0],zlims[1]);
1051       fOutput->Add(fhztracksinjet);
1052
1053       //check jet and tracks
1054       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);
1055       fOutput->Add(fhDiffPtTrPtJzok);
1056       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);
1057       fOutput->Add(fhDiffPtTrPtJzNok);
1058       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);
1059       fOutput->Add(fhDiffPzTrPtJzok);
1060       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);
1061       fOutput->Add(fhDiffPzTrPtJzNok);
1062       fhNtrkjzNok=new TH1I("hNtrkjzNok", "Number of tracks in a jet with z>1;N^{tracks} (z>1)",5,0,5);
1063       fOutput->Add(fhNtrkjzNok);
1064
1065       //calculate frag func with pt (simply ptD(or track)\cdot pt jet /ptjet^2)
1066       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]);
1067       fOutput->Add(fhzDT);
1068       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]);
1069       fOutput->Add(fhztracksinjetT);
1070
1071       fhIDddaugh=new TH1I("hIDddaugh", "ID of daughters;ID", 2001,-1000,1000);
1072       fOutput->Add(fhIDddaugh);
1073       fhIDddaughOut=new TH1I("hIDddaughOut", "ID of daughters not found in jet;ID", 2001,-1000,1000);
1074       fOutput->Add(fhIDddaughOut);
1075       fhIDjetTracks=new TH1I("hIDjetTracks", "ID of jet tracks;ID", 2001,-1000,1000);
1076       fOutput->Add(fhIDjetTracks);
1077
1078       fhDRdaughOut=new TH1F("hDRdaughOut", "#Delta R of daughters not belonging to the jet tracks (D in jet);#Delta R",200, 0.,2.);
1079       fOutput->Add(fhDRdaughOut);
1080
1081
1082       if(fCandidateType==kDstartoKpipi)
1083       {
1084          if(fSwitchOnSB){
1085             fhDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
1086             fhDiffSideBand->SetStats(kTRUE);
1087             fhDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
1088             fhDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
1089             fhDiffSideBand->Sumw2();
1090             fOutput->Add(fhDiffSideBand);
1091          }
1092
1093          fhPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
1094          fhPtPion->SetStats(kTRUE);
1095          fhPtPion->GetXaxis()->SetTitle("GeV/c");
1096          fhPtPion->GetYaxis()->SetTitle("Entries");
1097          fhPtPion->Sumw2();
1098          fOutput->Add(fhPtPion);
1099
1100       }
1101       // D related histograms
1102       fhNDPerEvNoJet=new TH1F("hNDPerEvNoJet","Number of candidates per event with no jets; N candidate/ev with no jet", 20, 0., 20.);
1103       fhNDPerEvNoJet->Sumw2();
1104       fOutput->Add(fhNDPerEvNoJet);
1105
1106       fhptDPerEvNoJet=new TH1F("hptDPerEvNoJet","pt distribution of candidates per events with no jets; p_{t}^{D} (GeV/c)",nbinsptD, ptDlims[0],ptDlims[1]);
1107       fhptDPerEvNoJet->Sumw2();
1108       fOutput->Add(fhptDPerEvNoJet);
1109
1110       if(fSwitchOnSparses){
1111          const Int_t    nAxisD=4;
1112          const Int_t    nbinsSparseD[nAxisD]={nbinsSpseta,nbinsSpsphi,nbinsSpsptD,nbinsSpsmass};
1113          const Double_t minSparseD[nAxisD]  ={etalims[0],philims[0],ptDlims[0],fMinMass};
1114          const Double_t maxSparseD[nAxisD]  ={etalims[1],philims[1],ptDlims[1],fMaxMass};
1115          fhsDstandalone=new THnSparseF("hsDstandalone","#phi, #eta, p_{T}^{D}, and mass", nAxisD, nbinsSparseD, minSparseD, maxSparseD);
1116          fhsDstandalone->Sumw2();
1117
1118          fOutput->Add(fhsDstandalone);
1119       }
1120
1121       /*
1122       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]);
1123       hPtJetWithD->Sumw2();
1124       //for the MC this histogram is filled with the real background
1125       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]);
1126       hPtJetWithDsb->Sumw2();
1127
1128       fOutput->Add(hPtJetWithD);
1129       fOutput->Add(hPtJetWithDsb);
1130
1131       */
1132       fhNJetPerEvNoD=new TH1F("hNJetPerEvNoD","Number of jets per event with no D; number of jets/ev with no D",10,-0.5,9.5);
1133       fhNJetPerEvNoD->Sumw2();
1134
1135       fhPtJetPerEvNoD=new TH1F("hPtJetPerEvNoD","pt distribution of jets per event with no D; #it{p}_{T}^{jet} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
1136       fhPtJetPerEvNoD->Sumw2();
1137
1138       fOutput->Add(fhNJetPerEvNoD);
1139       fOutput->Add(fhPtJetPerEvNoD);
1140
1141       fhDeltaRD=new TH1F("hDeltaRD",Form("#Delta R distribution of D candidates %s selected;#Delta R", fUseMCInfo ? "(S+B)" : ""),200, 0.,10.);
1142       fhDeltaRD->Sumw2();
1143       fOutput->Add(fhDeltaRD);
1144
1145       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]);
1146       fhDeltaRptDptj->Sumw2();
1147       fOutput->Add(fhDeltaRptDptj);
1148
1149       if(fUseMCInfo){
1150          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]);
1151          fhDeltaRptDptjB->Sumw2();
1152          fOutput->Add(fhDeltaRptDptjB);
1153       }
1154
1155       //background (side bands for the Dstar and like sign for D0)
1156       AliJetContainer *jetCont=GetJetContainer(0);
1157       if(!jetCont){
1158          Printf("Container 0 not found, try with name %s", fJetArrName.Data());
1159          jetCont=GetJetContainer(fJetArrName);
1160          if(!jetCont) Printf("NOT FOUND AGAIN");
1161          return kFALSE;
1162       }
1163       Printf("CONTAINER NAME IS %s", jetCont->GetArrayName().Data());
1164       //fJetRadius=jetCont->GetJetRadius();
1165       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]);
1166       fhInvMassptD->SetStats(kTRUE);
1167       fhInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
1168       fhInvMassptD->GetYaxis()->SetTitle("#it{p}_{t}^{D} (GeV/c)");
1169       fhInvMassptD->Sumw2();
1170
1171       fOutput->Add(fhInvMassptD);
1172
1173       if(fUseMCInfo){
1174          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]);
1175          fhInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
1176          fhInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
1177          fhInvMassptDbg->Sumw2();
1178          fOutput->Add(fhInvMassptDbg);
1179
1180       }
1181
1182       if(fSwitchOnSparses){
1183          Double_t pi=TMath::Pi(), philims2[2];
1184          philims2[0]=-pi*1./2.;
1185          philims2[1]=pi*3./2.;
1186          fhsDphiz=0x0; //definition below according to the switches
1187
1188          if(fSwitchOnSB && fSwitchOnPhiAxis && fSwitchOnOutOfConeAxis){
1189             AliInfo("Creating a 9 axes container: This might cause large memory usage");
1190             const Int_t nAxis=9;
1191             const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsphi,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,2, 2, 2, 2};
1192             const Double_t minSparse[nAxis]={zlims[0],philims2[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5, -0.5,-0.5,-0.5};
1193             const Double_t maxSparse[nAxis]={zlims[1],philims2[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5, 1.5 , 1.5};
1194             fNAxesBigSparse=nAxis;
1195
1196             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);
1197          }
1198
1199          if(!fSwitchOnPhiAxis || !fSwitchOnOutOfConeAxis || !fSwitchOnSB){
1200             fSwitchOnPhiAxis=0;
1201             fSwitchOnOutOfConeAxis=0;
1202             fSwitchOnSB=0;
1203             if(fUseMCInfo){
1204                AliInfo("Creating a 7 axes container (MB background candidates)");
1205                const Int_t nAxis=7;
1206                const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass,2, 2, 2};
1207                const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass, -0.5,-0.5,-0.5};
1208                const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5 , 1.5};
1209                fNAxesBigSparse=nAxis;
1210                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);
1211
1212             } else{
1213                AliInfo("Creating a 6 axes container");
1214                const Int_t nAxis=6;
1215                const Int_t nbinsSparse[nAxis]={nbinsSpsz,nbinsSpsptjet,nbinsSpsptD,nbinsSpsmass, 2, 2};
1216                const Double_t minSparse[nAxis]={zlims[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5,-0.5};
1217                const Double_t maxSparse[nAxis]={zlims[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5, 1.5};
1218                fNAxesBigSparse=nAxis;
1219
1220                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);
1221             }
1222          }
1223          if(!fhsDphiz) AliFatal("No THnSparse created");
1224          fhsDphiz->Sumw2();
1225
1226          fOutput->Add(fhsDphiz);
1227       }
1228    }
1229    PostData(1, fOutput);
1230
1231    return kTRUE;
1232 }
1233
1234 //_______________________________________________________________________________
1235
1236 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet,  AliAODEvent* aodEvent){
1237
1238    Double_t ptD=candidate->Pt();
1239    Double_t deltaR=DeltaR(jet,candidate);
1240    Double_t phiD=candidate->Phi();
1241    Double_t deltaphi = jet->Phi()-phiD;
1242    if(deltaphi<=-(TMath::Pi())/2.) deltaphi = deltaphi+2.*(TMath::Pi());
1243    if(deltaphi>(3.*(TMath::Pi()))/2.) deltaphi = deltaphi-2.*(TMath::Pi());
1244    Double_t z=Z(candidate,jet);
1245    /*
1246    if(z>1) {
1247       ((TH1I*)fOutput->FindObject("hControlDInJ"))->Fill(7);
1248       Double_t pmissing=TMath::Sqrt(fPmissing[0]*fPmissing[0]+fPmissing[1]*fPmissing[1]+fPmissing[2]*fPmissing[2]);
1249
1250       for(Int_t k=0;k<3;k++) ((TH1F*)fOutput->FindObject(Form("hMissP%d",k)))->Fill(fPmissing[k]);
1251
1252       ((TH1F*)fOutput->FindObject("hmissingp"))->Fill(pmissing);
1253       Double_t ptdiff=fPtJet-jet->Pt();
1254       ((TH1F*)fOutput->FindObject("hDeltaPtJet"))->Fill(ptdiff);
1255       ((TH1F*)fOutput->FindObject("hRelDeltaPtJet"))->Fill(ptdiff/jet->Pt());
1256
1257
1258    }
1259    */
1260    if(fIsDInJet)fhzDT->Fill(Z(candidate,jet,kTRUE));
1261
1262    fhDeltaRD->Fill(deltaR);
1263    fhDeltaRptDptj->Fill(deltaR,ptD,fPtJet);
1264
1265    Bool_t bDInEMCalAcc=InEMCalAcceptance(candidate);
1266    Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1267    if(fUseReco){
1268       if(fCandidateType==kD0toKpi) {
1269          AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;
1270
1271          FillHistogramsD0JetCorr(dzero,deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc, aodEvent);
1272
1273       }
1274
1275       if(fCandidateType==kDstartoKpipi) {
1276          AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;
1277          FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc);
1278
1279       }
1280    } else{ //generated level
1281       //AliInfo("Non reco");
1282       FillHistogramsMCGenDJetCorr(deltaphi,z,ptD,fPtJet,bDInEMCalAcc,bJetInEMCalAcc);
1283    }
1284 }
1285
1286 //_______________________________________________________________________________
1287
1288 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){
1289
1290
1291    Double_t masses[2]={0.,0.};
1292    Int_t pdgdaughtersD0[2]={211,321};//pi,K
1293    Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1294
1295    masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1296    masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1297
1298    //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1299    Double_t *point=0x0;
1300    if(fNAxesBigSparse==9){
1301       point=new Double_t[9];
1302       point[0]=z;
1303       point[1]=dPhi;
1304       point[2]=ptj;
1305       point[3]=ptD;
1306       point[4]=masses[0];
1307       point[5]=0;
1308       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1309       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1310       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1311    }
1312    if(fNAxesBigSparse==6){
1313       point=new Double_t[6];
1314       point[0]=z;
1315       point[1]=ptj;
1316       point[2]=ptD;
1317       point[3]=masses[0];
1318       point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1319       point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1320 }
1321   if(fNAxesBigSparse==7){
1322       point=new Double_t[7];
1323       point[0]=z;
1324       point[1]=ptj;
1325       point[2]=ptD;
1326       point[3]=masses[0];
1327       point[4]=0;
1328       point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1329       point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1330
1331    }
1332
1333
1334    //Printf("Candidate in FillHistogramsD0JetCorr IsA %s", (candidate->IsA())->GetName());
1335    Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
1336    if(isselected==1 || isselected==3) {
1337
1338       //if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[0],ptD);
1339
1340       FillMassHistograms(masses[0], ptD);
1341       if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1342    }
1343    if(isselected>=2) {
1344       //if(fIsDInJet) hPtJetWithD->Fill(ptj,masses[1],ptD);
1345
1346       FillMassHistograms(masses[1], ptD);
1347       if(fNAxesBigSparse==9) point[4]=masses[1];
1348       if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[1];
1349       if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1350    }
1351
1352 }
1353
1354 //_______________________________________________________________________________
1355
1356 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi,  Double_t z, Double_t ptD, Double_t ptj, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc){
1357   //dPhi and z not used at the moment,but will be (re)added
1358
1359    AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
1360    Double_t deltamass= dstar->DeltaInvMass();
1361    //Double_t massD0= dstar->InvMassD0();
1362
1363
1364    fhPtPion->Fill(softpi->Pt());
1365
1366    //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1367    Int_t isSB=0;//IsDzeroSideBand(dstar);
1368
1369    //Double_t point[]={z,dPhi,ptj,ptD,deltamass,static_cast<Double_t>(isSB), static_cast<Double_t>(fIsDInJet ? 1 : 0),bDInEMCalAcc,bJetInEMCalAcc};
1370    Double_t *point=0x0;
1371    if(fNAxesBigSparse==9){
1372       point=new Double_t[9];
1373       point[0]=z;
1374       point[1]=dPhi;
1375       point[2]=ptj;
1376       point[3]=ptD;
1377       point[4]=deltamass;
1378       point[5]=static_cast<Double_t>(isSB);
1379       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1380       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1381       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1382    }
1383    if(fNAxesBigSparse==6){
1384       point=new Double_t[6];
1385       point[0]=z;
1386       point[1]=ptj;
1387       point[2]=ptD;
1388       point[3]=deltamass;
1389       point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1390       point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1391    }
1392    if(fNAxesBigSparse==7){
1393       point=new Double_t[7];
1394       point[0]=z;
1395       point[1]=ptj;
1396       point[2]=ptD;
1397       point[3]=deltamass;
1398       point[4]=0;
1399       point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1400       point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1401    }
1402
1403    //if(fIsDInJet) hPtJetWithD->Fill(ptj,deltamass,ptD);
1404
1405    FillMassHistograms(deltamass, ptD);
1406    if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1407
1408 }
1409
1410 //_______________________________________________________________________________
1411
1412 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t z,Double_t ptD,Double_t ptjet, Bool_t bDInEMCalAcc, Bool_t bJetInEMCalAcc){
1413
1414    Double_t pdgmass=0;
1415    if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1416    if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
1417    //TH3F* hPtJetWithD=(TH3F*)fOutput->FindObject("hPtJetWithD");
1418
1419    //Double_t point[9]={z,dPhi,ptjet,ptD,pdgmass,0, static_cast<Double_t>(fIsDInJet ? 1 : 0),bDInEMCalAcc,bJetInEMCalAcc};
1420
1421    Double_t *point=0x0;
1422    if(fNAxesBigSparse==9){
1423       point=new Double_t[9];
1424       point[0]=z;
1425       point[1]=dPhi;
1426       point[2]=ptjet;
1427       point[3]=ptD;
1428       point[4]=pdgmass;
1429       point[5]=0;
1430       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1431       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1432       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1433    }
1434    if(fNAxesBigSparse==6){
1435       point=new Double_t[6];
1436       point[0]=z;
1437       point[1]=ptjet;
1438       point[2]=ptD;
1439       point[3]=pdgmass;
1440       point[4]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1441       point[5]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1442    }
1443       if(fNAxesBigSparse==7){
1444       point=new Double_t[7];
1445       point[0]=z;
1446       point[1]=ptjet;
1447       point[2]=ptD;
1448       point[3]=pdgmass;
1449       point[4]=1;
1450       point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1451       point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1452    }
1453
1454
1455
1456    if(fNAxesBigSparse==9) point[4]=pdgmass;
1457    if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=pdgmass;
1458    if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1459    //if(fIsDInJet) {
1460    //  hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet
1461    //}
1462
1463 }
1464
1465 //_______________________________________________________________________________
1466
1467 void AliAnalysisTaskFlavourJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD){
1468
1469    if(fIsDInJet) {
1470       fhInvMassptD->Fill(mass,ptD);
1471    }
1472 }
1473
1474 //________________________________________________________________________________
1475
1476 void AliAnalysisTaskFlavourJetCorrelations::FlagFlavour(AliEmcalJet *jet){
1477
1478    AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
1479    if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;
1480    if (fIsDInJet) jet->AddFlavourTag(tag);
1481
1482    return;
1483
1484 }
1485
1486 //_______________________________________________________________________________
1487
1488 void AliAnalysisTaskFlavourJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){
1489
1490    //  D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
1491    // (expected detector resolution) on the left and right frm the D0 mass. Each band
1492    //  has a width of ~5 sigmas. Two band needed  for opening angle considerations
1493
1494    //TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
1495
1496    Bool_t bDInEMCalAcc=InEMCalAcceptance(candDstar);
1497    Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1498
1499    Double_t deltaM=candDstar->DeltaInvMass();
1500    //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);
1501    Double_t z=Z(candDstar,jet);
1502    Double_t ptD=candDstar->Pt();
1503
1504    Double_t dPhi=jet->Phi()-candDstar->Phi();
1505
1506    if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1507    if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1508
1509    Int_t isSideBand=1;//IsDzeroSideBand(candDstar);
1510    //if no SB analysis we should not enter here, so no need of checking the number of axes
1511    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)};
1512    //fill the background histogram with the side bands or when looking at MC with the real background
1513    if(isSideBand==1){
1514       fhDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background
1515       //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);
1516       if(fSwitchOnSparses) fhsDphiz->Fill(point,1.);
1517       //if(fIsDInJet){  // evaluate in the near side
1518       //         hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
1519       //}
1520    }
1521 }
1522
1523 //_______________________________________________________________________________
1524
1525 void AliAnalysisTaskFlavourJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg,AliEmcalJet* jet){
1526
1527    //need mass, deltaR, pt jet, ptD
1528    //two cases: D0 and Dstar
1529
1530    Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());
1531    if(!isselected) return;
1532
1533    Double_t ptD=candbg->Pt();
1534    Double_t deltaR=DeltaR(jet,candbg);
1535    Double_t phiD=candbg->Phi();
1536    Double_t deltaphi = jet->Phi()-phiD;
1537    if(deltaphi<=-(TMath::Pi())/2.) deltaphi = deltaphi+2.*(TMath::Pi());
1538    if(deltaphi>(3.*(TMath::Pi()))/2.) deltaphi = deltaphi-2.*(TMath::Pi());
1539    Double_t z=Z(candbg,jet);
1540
1541    if(fIsDInJet) fhzDT->Fill(Z(candbg,jet,kTRUE));
1542
1543
1544
1545
1546    fhDeltaRD->Fill(deltaR);
1547    fhDeltaRptDptjB->Fill(deltaR,ptD,fPtJet);
1548
1549    Bool_t bDInEMCalAcc=InEMCalAcceptance(candbg);
1550    Bool_t bJetInEMCalAcc=InEMCalAcceptance(jet);
1551
1552    //TH3F* hPtJetWithDsb=(TH3F*)fOutput->FindObject("hPtJetWithDsb");
1553
1554    Double_t *point=0x0;
1555    if(fNAxesBigSparse==9){
1556       point=new Double_t[9];
1557       point[0]=z;
1558       point[1]=deltaphi;
1559       point[2]=fPtJet;
1560       point[3]=ptD;
1561       point[4]=-999; //set below
1562       point[5]=1;
1563       point[6]=static_cast<Double_t>(fIsDInJet ? 1 : 0);
1564       point[7]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1565       point[8]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1566    }
1567
1568    if(fNAxesBigSparse==7){
1569       point=new Double_t[7];
1570       point[0]=z;
1571       point[1]=fPtJet;
1572       point[2]=ptD;
1573       point[3]=-999; //set below
1574       point[4]=1;
1575       point[5]=static_cast<Double_t>(bDInEMCalAcc ? 1 : 0);
1576       point[6]=static_cast<Double_t>(bJetInEMCalAcc ? 1 : 0);
1577    }
1578
1579    if(fCandidateType==kDstartoKpipi){
1580       AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;
1581       Double_t deltaM=dstarbg->DeltaInvMass();
1582       fhInvMassptDbg->Fill(deltaM,ptD);
1583       //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,deltaM,ptD);
1584       if(fNAxesBigSparse==9) point[4]=deltaM;
1585       if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=deltaM;
1586       if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1587    }
1588
1589    if(fCandidateType==kD0toKpi){
1590       Double_t masses[2]={0.,0.};
1591       Int_t pdgdaughtersD0[2]={211,321};//pi,K
1592       Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1593
1594       masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1595       masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1596
1597
1598       if(isselected==1 || isselected==3) {
1599          //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[0],ptD);
1600          fhInvMassptDbg->Fill(masses[0],ptD);
1601          if(fNAxesBigSparse==9) point[4]=masses[0];
1602          if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[0];
1603          if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1604      }
1605       if(isselected>=2) {
1606          //if(fIsDInJet) hPtJetWithDsb->Fill(fPtJet,masses[1],ptD);
1607          fhInvMassptDbg->Fill(masses[1],ptD);
1608          if(fNAxesBigSparse==9) point[4]=masses[1];
1609          if(fNAxesBigSparse==6 || fNAxesBigSparse==7) point[3]=masses[1];
1610          if(fSwitchOnSparses && (fSwitchOnOutOfConeAxis || fIsDInJet)) fhsDphiz->Fill(point,1.);
1611
1612       }
1613
1614
1615    }
1616 }
1617
1618 //_______________________________________________________________________________
1619
1620 Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliEmcalJet *p1, AliVParticle *p2) const {
1621    //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1622    //It recalculates the eta-phi values if it was asked for background subtraction of the jets
1623    if(!p1 || !p2) return -1;
1624
1625     Double_t phi1=p1->Phi(),eta1=p1->Eta();
1626
1627     //It subtracts the backgroud of jets if it was asked for it.
1628     TString sname("Rho");
1629     if (!sname.IsNull()) {
1630         AliRhoParameter *rho = 0;
1631         rho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(sname));
1632         if(rho){
1633             Double_t pj[3];
1634             Bool_t okpj=p1->PxPyPz(pj);
1635             if(!okpj){
1636                 printf("Problems getting momenta\n");
1637                 return -999;
1638             }
1639             Double_t rhoval = rho->GetVal();
1640             pj[0] = p1->Px() - p1->Area()*(rhoval*TMath::Cos(p1->AreaPhi()));
1641             pj[1] = p1->Py() - p1->Area()*(rhoval*TMath::Sin(p1->AreaPhi()));
1642             pj[2] = p1->Pz() - p1->Area()*(rhoval*TMath::SinH(p1->AreaEta()));
1643             //Image of the function Arccos(px/pt) where pt = sqrt(px*px+py*py) is:
1644             //[0,pi]    if py > 0 and
1645             //[pi,2pi]  if py < 0
1646             phi1 = TMath::ACos(pj[0]/TMath::Sqrt(pj[0]*pj[0]+pj[1]*pj[1]));
1647             if(pj[1]<0) phi1 = 2*TMath::Pi() - phi1;
1648             eta1 = TMath::ASinH(pj[2]/TMath::Sqrt(pj[0]*pj[0]+pj[1]*pj[1]));
1649         }
1650     }
1651
1652    Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
1653
1654    Double_t dPhi=phi1-phi2;
1655    if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1656    if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1657
1658    Double_t dEta=eta1-eta2;
1659    Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
1660    return deltaR;
1661
1662 }
1663
1664 //_______________________________________________________________________________
1665
1666 Int_t AliAnalysisTaskFlavourJetCorrelations::IsDzeroSideBand(AliAODRecoCascadeHF *candDstar){
1667
1668    Double_t ptD=candDstar->Pt();
1669    Int_t bin = fCuts->PtBin(ptD);
1670    if (bin < 0){
1671       // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
1672       bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?
1673       return -1;  // out of bounds
1674    }
1675
1676    Double_t invM=candDstar->InvMassD0();
1677    Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1678
1679    Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
1680    Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
1681
1682    if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
1683    else return 0;
1684
1685 }
1686
1687 //_______________________________________________________________________________
1688
1689 Bool_t AliAnalysisTaskFlavourJetCorrelations::AreDaughtersInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Int_t*& daughOutOfJetID, AliAODTrack**& dtrks, Bool_t fillH){
1690    //returns true/false and the array of particles not found among jet constituents
1691
1692
1693    Int_t daughtersID[3];
1694    Int_t ndaugh=0;
1695    Int_t dmatchedID[3]={0,0,0};
1696    Int_t countmatches=0;
1697    if (fCandidateType==kDstartoKpipi){
1698       AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)charm;
1699       AliAODRecoDecayHF2Prong* dzero=dstar->Get2Prong();
1700       daughtersID[0]=dzero->GetProngID(0);
1701       daughtersID[1]=dzero->GetProngID(1);
1702       daughtersID[2]=dstar->GetBachelor()->GetID();
1703       ndaugh=3;
1704
1705       dtrks=new AliAODTrack*[3];
1706       dtrks[0]=(AliAODTrack*)dzero->GetDaughter(0);
1707       dtrks[1]=(AliAODTrack*)dzero->GetDaughter(1);
1708       dtrks[2]=(AliAODTrack*)dstar->GetBachelor();
1709
1710       //check
1711       if(fillH){
1712          if(daughtersID[0]!=((AliAODTrack*)dtrks[0])->GetID() || daughtersID[1]!=((AliAODTrack*)dtrks[1])->GetID())  fhControlDInJ->Fill(6);
1713
1714          fhIDddaugh->Fill(daughtersID[0]);
1715          fhIDddaugh->Fill(daughtersID[1]);
1716          fhIDddaugh->Fill(daughtersID[2]);
1717
1718       }
1719       //Printf("ID daughters %d, %d, %d",daughtersID[0], daughtersID[1], daughtersID[2]);
1720    }
1721
1722    if (fCandidateType==kD0toKpi){
1723       daughtersID[0]=charm->GetProngID(0);
1724       daughtersID[1]=charm->GetProngID(1);
1725       ndaugh=2;
1726       if(fillH){
1727          fhIDddaugh->Fill(daughtersID[0]);
1728          fhIDddaugh->Fill(daughtersID[1]);
1729       }
1730       dtrks=new AliAODTrack*[2];
1731       dtrks[0]=(AliAODTrack*)charm->GetDaughter(0);
1732       dtrks[1]=(AliAODTrack*)charm->GetDaughter(1);
1733
1734    }
1735
1736    const Int_t cndaugh=ndaugh;
1737    daughOutOfJetID=new Int_t[cndaugh];
1738
1739    Int_t nchtrk=thejet->GetNumberOfTracks();
1740    for(Int_t itk=0;itk<nchtrk;itk++){
1741       AliVTrack* tkjet=((AliPicoTrack*)thejet->TrackAt(itk,fTrackArr))->GetTrack();
1742       if(!tkjet) continue;
1743       Int_t idtkjet=tkjet->GetID();
1744       if(fillH) fhIDjetTracks->Fill(idtkjet);
1745       for(Int_t id=0;id<ndaugh;id++){
1746          if(idtkjet==daughtersID[id]) {
1747             dmatchedID[id]=daughtersID[id]; //daughter which matches a track in the jet
1748             countmatches++;
1749
1750          }
1751
1752          if(countmatches==ndaugh) break;
1753       }
1754    }
1755    //reverse: include in the array the ID of daughters not matching
1756
1757    for(Int_t id=0;id<ndaugh;id++){
1758       if(dmatchedID[id]==0) {
1759          daughOutOfJetID[id]=daughtersID[id];
1760          if(fillH) fhIDddaughOut->Fill(daughOutOfJetID[id]);
1761       }
1762       else daughOutOfJetID[id]=0;
1763    }
1764    if(fillH){
1765       if((ndaugh-countmatches) == 1) fhControlDInJ->Fill(0);
1766       if((ndaugh-countmatches) == 2) fhControlDInJ->Fill(1);
1767       if((ndaugh-countmatches) == 3) fhControlDInJ->Fill(2);
1768    }
1769    if(ndaugh!=countmatches){
1770       return kFALSE;
1771    }
1772
1773    return kTRUE;
1774 }
1775
1776 //_______________________________________________________________________________
1777
1778 Bool_t AliAnalysisTaskFlavourJetCorrelations::IsDInJet(AliEmcalJet *thejet, AliAODRecoDecayHF* charm, Bool_t fillH){
1779
1780    //check the conditions type:
1781    //type 0 : DeltaR < jet radius, don't check daughters (and don't correct pt jet)
1782    //type 1 : DeltaR < jet radius and check for all daughters among jet tracks, don't correct ptjet
1783    //type 2 (default) : DeltaR < jet radius and check for all daughters among jet tracks, if not present, correct the ptjet
1784    //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
1785
1786    fPmissing[0]=0;
1787    fPmissing[1]=0;
1788    fPmissing[2]=0;
1789
1790    Bool_t testDeltaR=kFALSE;
1791    if(DeltaR(thejet,charm)<fJetRadius) testDeltaR=kTRUE;
1792
1793    Int_t* daughOutOfJet;
1794    AliAODTrack** charmDaugh;
1795    Bool_t testDaugh=AreDaughtersInJet(thejet, charm, daughOutOfJet,charmDaugh,fillH);
1796    if(testDaugh && testDeltaR) {
1797       //Z of candidates with daughters in the jet
1798       fhzDinjet->Fill(Z(charm,thejet));
1799
1800    }
1801    if(!testDaugh && testDeltaR && fTypeDInJet>=2){
1802
1803       Int_t ndaugh=3;
1804       if(fCandidateType==kD0toKpi) ndaugh=2;
1805       Int_t nOut=ndaugh;
1806
1807       for(Int_t id=0;id<ndaugh;id++){
1808          if(daughOutOfJet[id]!=0){ //non-matched daughter
1809             //get track and its momentum
1810             nOut--;
1811             if(fillH) {
1812                fhControlDInJ->Fill(3);
1813                if(id<2) fhControlDInJ->Fill(4);
1814                if(id==2)fhControlDInJ->Fill(5);
1815                fhDRdaughOut->Fill(DeltaR(thejet, charmDaugh[id]));
1816             }
1817             if(fTypeDInJet==2){
1818                fPmissing[0]+=charmDaugh[id]->Px();
1819                fPmissing[1]+=charmDaugh[id]->Py();
1820                fPmissing[2]+=charmDaugh[id]->Pz();
1821             }
1822             if(fTypeDInJet==3){
1823                Double_t ptdaug  = charmDaugh[id]->Pt();
1824                Double_t ptjet   = thejet->Pt();
1825                Double_t ptn     = ptjet+ptdaug;
1826                Double_t phidaug = charmDaugh[id]->Phi();
1827                Double_t phijet  = thejet->Phi();
1828                Double_t phin    = (phijet/ptjet+phidaug/ptdaug)/(1./ptjet+ 1./ptdaug);
1829                Double_t etadaug = charmDaugh[id]->Eta();
1830                Double_t etajet  = thejet->Eta();
1831                Double_t etan    = (etajet/ptjet+etadaug/ptdaug)/(1./ptjet+ 1./ptdaug);
1832
1833                fPmissing[0]+= ptn*TMath::Cos(phin);
1834                fPmissing[1]+= ptn*TMath::Sin(phin);
1835                fPmissing[2]+= ptn*TMath::SinH(etan);
1836             }
1837          }
1838
1839       }
1840
1841       //now the D in within the jet
1842       testDaugh=kTRUE;
1843
1844    }
1845
1846    delete[] charmDaugh;
1847
1848    Bool_t result=0;
1849    switch(fTypeDInJet){
1850    case 0:
1851       result=testDeltaR;
1852       break;
1853    case 1:
1854       result=testDeltaR && testDaugh;
1855       break;
1856    case 2: //this case defines fPmissing
1857       result=testDeltaR && testDaugh;
1858       break;
1859    default:
1860       AliInfo("Selection type not specified, use 1");
1861       result=testDeltaR && testDaugh;
1862       break;
1863    }
1864  return result;
1865 }
1866
1867 Bool_t AliAnalysisTaskFlavourJetCorrelations::InEMCalAcceptance(AliVParticle *vpart){
1868    //check eta phi of a VParticle: return true if it is in the EMCal acceptance, false otherwise
1869
1870    Double_t phiEMCal[2]={1.405,3.135},etaEMCal[2]={-0.7,0.7};
1871    Bool_t binEMCal=kTRUE;
1872    Double_t phi=vpart->Phi(), eta=vpart->Eta();
1873    if(phi<phiEMCal[0] || phi>phiEMCal[1]) binEMCal=kFALSE;
1874    if(eta<etaEMCal[0] || eta>etaEMCal[1]) binEMCal=kFALSE;
1875    return binEMCal;
1876
1877
1878 }