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