PWGJE/FlavourJetTasks/AliAnalysisTaskFlavourFilterAndJetCorrelations.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
  34 #include "AliAnalysisTaskFlavourFilterAndJetCorrelations.h"
  35 #include "AliAODMCHeader.h"
  36 #include "AliAODHandler.h"
  37 #include "AliAnalysisManager.h"
  38 #include "AliLog.h"
  39 #include "AliEmcalJet.h"
  40 #include "AliJetContainer.h"
  41 #include "AliAODRecoDecay.h"
  42 #include "AliAODRecoCascadeHF.h"
  43 #include "AliAODRecoDecayHF2Prong.h"
  44 #include "AliESDtrack.h"
  45 #include "AliAODMCParticle.h"
  46 #include "AliPicoTrack.h"
  47 #include "AliRDHFCutsD0toKpi.h"
  48 #include "AliRDHFCutsDStartoKpipi.h"
  49
  50 ClassImp(AliAnalysisTaskFlavourFilterAndJetCorrelations)
  51
  52
  53 //_______________________________________________________________________________
  54
  55 AliAnalysisTaskFlavourFilterAndJetCorrelations::AliAnalysisTaskFlavourFilterAndJetCorrelations() :
  56 AliAnalysisTaskEmcalJet("",kFALSE),
  57 fUseMCInfo(kTRUE),
  58 fUseReco(kTRUE),
  59 fCandidateType(),
  60 fPDGmother(),
  61 fNProngs(),
  62 fPDGdaughters(),
  63 fBranchName(),
  64 fmyOutput(0),
  65 fCuts(0),
  66 fMinMass(),
  67 fMaxMass(),
  68 fJetArrName(0),
  69 fCandArrName(0),
  70 fLeadingJetOnly(kFALSE),
  71 fJetRadius(0)
  72 {
  73    //
  74    // Default ctor
  75    //
  76    //SetMakeGeneralHistograms(kTRUE);
  77
  78 }
  79
  80 //_______________________________________________________________________________
  81
  82 AliAnalysisTaskFlavourFilterAndJetCorrelations::AliAnalysisTaskFlavourFilterAndJetCorrelations(const Char_t* name, AliRDHFCuts* cuts,ECandidateType candtype) :
  83 AliAnalysisTaskEmcalJet(name,kFALSE),
  84 fUseMCInfo(kTRUE),
  85 fUseReco(kTRUE),
  86 fCandidateType(),
  87 fPDGmother(),
  88 fNProngs(),
  89 fPDGdaughters(),
  90 fBranchName(),
  91 fmyOutput(0),
  92 fmyOutputF(0),
  93 fCuts(0),
  94 fMinMass(),
  95 fMaxMass(),
  96 fJetArrName(0),
  97 fCandArrName(0),
  98 fLeadingJetOnly(kFALSE),
  99 fJetRadius(0)
 100 {
 101    //
 102    // Constructor. Initialization of Inputs and Outputs
 103    //
 104
 105    Info("AliAnalysisTaskFlavourFilterAndJetCorrelations","Calling Constructor");
 106    fCuts=cuts;
 107    fCandidateType=candtype;
 108    const Int_t nptbins=fCuts->GetNPtBins();
 109    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};
 110
 111    switch(fCandidateType){
 112    case 0:
 113       fPDGmother=421;
 114       fNProngs=2;
 115       fPDGdaughters[0]=211;//pi
 116       fPDGdaughters[1]=321;//K
 117       fPDGdaughters[2]=0; //empty
 118       fPDGdaughters[3]=0; //empty
 119       fBranchName="D0toKpi";
 120       fCandArrName="D0";
 121       break;
 122    case 1:
 123       fPDGmother=413;
 124       fNProngs=3;
 125       fPDGdaughters[1]=211;//pi soft
 126       fPDGdaughters[0]=421;//D0
 127       fPDGdaughters[2]=211;//pi fromD0
 128       fPDGdaughters[3]=321; // K from D0
 129       fBranchName="Dstar";
 130       fCandArrName="Dstar";
 131
 132       if(nptbins<=13){
 133          for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]= defaultSigmaD013[ipt];
 134       } else {
 135          AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));
 136       }
 137       break;
 138    default:
 139       printf("%d not accepted!!\n",fCandidateType);
 140       break;
 141    }
 142
 143    if(fCandidateType==kD0toKpi)SetMassLimits(0.15,fPDGmother);
 144    if(fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);
 145
 146    DefineOutput(1,TList::Class()); // histos
 147    DefineOutput(2,AliRDHFCuts::Class()); // my cuts
 148    DefineOutput(3,TList::Class()); //histos filter
 149
 150 }
 151
 152 //_______________________________________________________________________________
 153
 154 AliAnalysisTaskFlavourFilterAndJetCorrelations::~AliAnalysisTaskFlavourFilterAndJetCorrelations() {
 155    //
 156    // destructor
 157    //
 158
 159    Info("~AliAnalysisTaskFlavourFilterAndJetCorrelations","Calling Destructor");
 160
 161    delete fmyOutput;
 162    delete fmyOutputF;
 163    delete fCuts;
 164
 165 }
 166
 167 //_______________________________________________________________________________
 168
 169 void AliAnalysisTaskFlavourFilterAndJetCorrelations::Init(){
 170    //
 171    // Initialization
 172    //
 173
 174    if(fDebug > 1) printf("AnalysisTaskRecoJetCorrelations::Init() \n");
 175    switch(fCandidateType){
 176    case 0:
 177       {
 178          AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
 179          copyfCuts->SetName("AnalysisCutsDzero");
 180          // Post the data
 181          PostData(2,copyfCuts);
 182       }
 183       break;
 184    case 1:
 185       {
 186          AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
 187          copyfCuts->SetName("AnalysisCutsDStar");
 188          // Post the cuts
 189          PostData(2,copyfCuts);
 190       }
 191       break;
 192    default:
 193       return;
 194    }
 195
 196    return;
 197 }
 198
 199 //_______________________________________________________________________________
 200
 201 void AliAnalysisTaskFlavourFilterAndJetCorrelations::UserCreateOutputObjects() {
 202    // output
 203    Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 204    fmyOutput = new TList();
 205    fmyOutput->SetOwner();
 206    fmyOutput->SetName("pippo");
 207    // define histograms
 208    DefineHistoForAnalysis();
 209    PostData(1,fmyOutput);
 210
 211    fmyOutputF = new TList();
 212    fmyOutputF->SetOwner();
 213    fmyOutputF->SetName("pluto");
 214    // define histograms
 215    DefineHistoForFiltering();
 216    PostData(3,fmyOutputF);
 217
 218    return;
 219 }
 220
 221 //_______________________________________________________________________________
 222
 223 void AliAnalysisTaskFlavourFilterAndJetCorrelations::UserExec(Option_t *)
 224 {
 225    // user exec
 226    if (!fInitialized){
 227       AliAnalysisTaskEmcalJet::ExecOnce();
 228    }
 229
 230    // Load the event
 231    AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 232
 233    TClonesArray *arrayDStartoD0pi=0;
 234    TClonesArray *trackArr = 0;
 235
 236    if (!aodEvent && AODEvent() && IsStandardAOD()) {
 237
 238       // In case there is an AOD handler writing a standard AOD, use the AOD
 239       // event in memory rather than the input (ESD) event.
 240       aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
 241
 242       // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 243       // have to taken from the AOD event hold by the AliAODExtension
 244       AliAODHandler* aodHandler = (AliAODHandler*)
 245       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 246       if(aodHandler->GetExtensions()) {
 247          AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 248          AliAODEvent *aodFromExt = ext->GetAOD();
 249          arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());
 250       }
 251    } else {
 252       arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());
 253    }
 254
 255    if (!arrayDStartoD0pi) {
 256       AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));
 257       //  return;
 258    } else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
 259
 260    TClonesArray* mcArray = 0x0;
 261    if (fUseMCInfo) {
 262       mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 263       if (!mcArray) {
 264          printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
 265          return;
 266       }
 267    }
 268
 269    //retrieve jets
 270    trackArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("PicoTracks"));
 271    //clusArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClustersCorr"));
 272    //jetArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrName));
 273    fJetRadius=GetJetContainer(0)->GetJetRadius();
 274
 275    if(!trackArr){
 276       AliInfo(Form("Could not find the track array\n"));
 277       return;
 278    }
 279
 280
 281    //FILTER D MESONS
 282    fCandidateArray = new TClonesArray("AliAODRecoDecayHF",0);
 283    fCandidateArray->SetName(Form("fCandidateArray%s%s",fCandArrName.Data(),fUseReco ? "rec" : "gen"));
 284
 285    if (fCandidateType==kDstartoKpipi){
 286       fSideBandArray = new TClonesArray("AliAODRecoCascadeHF",0); //this is for the DStar only!
 287       fSideBandArray->SetName(Form("fSideBandArray%s%s",fCandArrName.Data(),fUseReco ? "rec" : "gen"));
 288    }
 289
 290    //Histograms
 291    TH1I* hstatF = (TH1I*)fmyOutputF->FindObject("hstatF");
 292    TH1F* hnSBCandEv=(TH1F*)fmyOutputF->FindObject("hnSBCandEv");
 293    TH1F* hnCandEv=(TH1F*)fmyOutputF->FindObject("hnCandEv");
 294    TH2F* hInvMassptDF = (TH2F*)fmyOutputF->FindObject("hInvMassptDF");
 295
 296    TH1F* hPtPion=0x0;
 297    if (fCandidateType==kDstartoKpipi) hPtPion = (TH1F*)fmyOutputF->FindObject("hPtPion");
 298    hstatF->Fill(0);
 299
 300    // fix for temporary bug in ESDfilter
 301    // the AODs with null vertex pointer didn't pass the PhysSel
 302    if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
 303
 304    //Event selection
 305    Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
 306    //TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
 307    if (!iseventselected) return;
 308    hstatF->Fill(1);
 309
 310    const Int_t nD = arrayDStartoD0pi->GetEntriesFast();
 311    if(fUseReco) hstatF->Fill(2, nD);
 312
 313    //D* and D0 prongs needed to MatchToMC method
 314    Int_t pdgDgDStartoD0pi[2] = { 421, 211 };  // D0,pi
 315    Int_t pdgDgD0toKpi[2] = { 321, 211 };      // K, pi
 316
 317    //D0 from D0 bar
 318    Int_t pdgdaughtersD0[2] = { 211, 321 };     // pi,K
 319    Int_t pdgdaughtersD0bar[2] = { 321, 211 };  // K,pi
 320
 321    Int_t iCand =0;
 322    Int_t iSBCand=0;
 323    Int_t isSelected = 0;
 324    AliAODRecoDecayHF *charmCand = 0;
 325    AliAODMCParticle *charmPart = 0;
 326    Bool_t isMCBkg=kFALSE;
 327
 328    Double_t mPDGD0 = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 329
 330    Int_t mcLabel = -9999;
 331    Int_t pdgMeson = 413;
 332    if (fCandidateType==kD0toKpi) pdgMeson = 421;
 333
 334    for (Int_t icharm=0; icharm<nD; icharm++) {   //loop over D candidates
 335       charmCand = (AliAODRecoDecayHF*)arrayDStartoD0pi->At(icharm); // D candidates
 336       if (!charmCand) continue;
 337
 338
 339       if (fUseMCInfo) { // Look in MC, try to simulate the z
 340          if (fCandidateType==kDstartoKpipi) {
 341             AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
 342             mcLabel = temp->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
 343          }
 344
 345          if (fCandidateType==kD0toKpi)
 346             mcLabel = charmCand->MatchToMC(421,mcArray,fNProngs,fPDGdaughters);
 347
 348          if (mcLabel<=0) isMCBkg=kTRUE;
 349          else hstatF->Fill(2);
 350          if (!isMCBkg) charmPart=(AliAODMCParticle*)mcArray->At(mcLabel);
 351       }
 352
 353       Double_t ptD = charmCand->Pt();
 354
 355       // region of interest + cuts
 356       if (!fCuts->IsInFiducialAcceptance(ptD,charmCand->Y(pdgMeson))) continue;
 357
 358       if(!fUseMCInfo && fCandidateType==kDstartoKpipi){
 359          //select by track cuts the side band candidates (don't want mass cut)
 360          isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kTracks,aodEvent);
 361          if (!isSelected) continue;
 362          //add a reasonable cut on the invariant mass (e.g. (+-2\sigma, +-10 \sigma), with \sigma = fSigmaD0[bin])
 363          Int_t bin = fCuts->PtBin(ptD);
 364          if(bin<0 || bin>=fCuts->GetNPtBins()) {
 365             AliError(Form("Pt %.3f out of bounds",ptD));
 366             continue;
 367          }
 368          AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
 369          //if data and Dstar from D0 side band
 370          if (((temp->InvMassD0()<=(mPDGD0-3.*fSigmaD0[bin])) && (temp->InvMassD0()>(mPDGD0-10.*fSigmaD0[bin]))) /*left side band*/||  ((temp->InvMassD0()>=(mPDGD0+3.*fSigmaD0[bin])) && (temp->InvMassD0()<(mPDGD0+10.*fSigmaD0[bin])))/*right side band*/){
 371
 372             new ((*fSideBandArray)[iSBCand]) AliAODRecoCascadeHF(*temp);
 373             iSBCand++;
 374          }
 375       }
 376       //candidate selected by cuts and PID
 377       isSelected = fCuts->IsSelected(charmCand,AliRDHFCuts::kAll,aodEvent); //selected
 378       if (!isSelected) continue;
 379
 380       //for data and MC signal fill fCandidateArray
 381       if(!fUseMCInfo || (fUseMCInfo && !isMCBkg)){
 382          // for data or MC with the requirement fUseReco fill with candidates
 383          if(fUseReco) {
 384             new ((*fCandidateArray)[iCand]) AliAODRecoDecayHF(*charmCand);
 385             //Printf("Filling reco");
 386             hstatF->Fill(3);
 387          }
 388          // for MC with requirement particle level fill with AliAODMCParticle
 389          else if (fUseMCInfo) {
 390             new ((*fCandidateArray)[iCand]) AliAODMCParticle(*charmPart);
 391             //Printf("Filling gen");
 392             hstatF->Fill(3);
 393          }
 394
 395          iCand++;
 396       }
 397       //for MC background fill fSideBandArray (which is instead filled above for DStar in case of data for the side bands candidates)
 398       else if(fUseReco){
 399          new ((*fSideBandArray)[iSBCand]) AliAODRecoDecayHF(*charmCand);
 400          iSBCand++;
 401       }
 402
 403
 404       Double_t masses[2];
 405       if (fCandidateType==kDstartoKpipi) { //D*->D0pi->Kpipi
 406
 407          //softpion from D* decay
 408          AliAODRecoCascadeHF *temp = (AliAODRecoCascadeHF*)charmCand;
 409          AliAODTrack *track2 = (AliAODTrack*)temp->GetBachelor();
 410
 411          // select D* in the D0 window.
 412          // In the cut object window is loose to allow for side bands
 413
 414
 415          // retrieve the corresponding pt bin for the candidate
 416          // and set the expected D0 width (x3)
 417          //    static const Int_t n = fCuts->GetNPtBins();
 418          Int_t bin = fCuts->PtBin(ptD);
 419          if(bin<0 || bin>=fCuts->GetNPtBins()) {
 420             AliError(Form("Pt %.3f out of bounds",ptD));
 421             continue;
 422          }
 423
 424          AliInfo(Form("Pt bin %d and sigma D0 %.4f",bin,fSigmaD0[bin]));
 425          //consider the Dstar candidates only if the mass of the D0 is in 3 sigma wrt the PDG value
 426          if ((temp->InvMassD0()>=(mPDGD0-3.*fSigmaD0[bin])) && (temp->InvMassD0()<=(mPDGD0+3.*fSigmaD0[bin]))) {
 427             masses[0] = temp->DeltaInvMass(); //D*
 428             masses[1] = 0.; //dummy for D*
 429
 430             //D*  delta mass
 431             hInvMassptDF->Fill(masses[0], ptD); // 2 D slice for pt bins
 432
 433             // D* pt and soft pion pt for good candidates
 434             Double_t mPDGDstar = TDatabasePDG::Instance()->GetParticle(413)->Mass();
 435             Double_t invmassDelta = temp->DeltaInvMass();
 436             if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0021) hPtPion->Fill(track2->Pt());
 437          }
 438       } //Dstar specific
 439
 440       if (fCandidateType==kD0toKpi) { //D0->Kpi
 441
 442          //needed quantities
 443          masses[0] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0); //D0
 444          masses[1] = charmCand->InvMass(fNProngs, (UInt_t*)pdgdaughtersD0bar); //D0bar
 445          hstatF->Fill(3);
 446
 447          // mass vs pt
 448          if (isSelected==1 || isSelected==3) hInvMassptDF->Fill(masses[0],ptD);
 449          if (isSelected>=2) hInvMassptDF->Fill(masses[1],ptD);
 450       } //D0 specific
 451
 452       charmCand = 0;
 453       isMCBkg=kFALSE;
 454    } // end of D cand loop
 455
 456    hnCandEv->Fill(fCandidateArray->GetEntriesFast());
 457    if (fCandidateType==kDstartoKpipi || fUseMCInfo) {
 458       Int_t nsbcand=fSideBandArray->GetEntriesFast();
 459       hstatF->Fill(4,nsbcand);
 460       hnSBCandEv->Fill(nsbcand);
 461    }
 462
 463    //CORRELATION WITH JETS
 464    /*
 465    TString arrname="fCandidateArray";
 466    TString candarrname=Form("%s%s%s",arrname.Data(),fCandArrName.Data(),fUseReco ? "rec" : "gen");
 467    fCandidateArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(candarrname));
 468    if (!fCandidateArray) {
 469       Printf("%s array not found",candarrname.Data());
 470       InputEvent()->GetList()->ls();
 471       return;
 472    }
 473    //Printf("ncandidates found %d",fCandidateArray->GetEntriesFast());
 474
 475    TString arrSBname="fSideBandArray";
 476    TString sbarrname=Form("%s%s%s",arrSBname.Data(),fCandArrName.Data(),fUseReco ? "rec" : "gen");
 477    fSideBandArray = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sbarrname));
 478    if (fCandidateType==1 && !fSideBandArray) {
 479       Printf("%s array not found",sbarrname.Data());
 480       InputEvent()->GetList()->ls();
 481       return;
 482    }
 483    //Printf("nSBCand or Bkg found %d",fSideBandArray->GetEntriesFast());
 484    */
 485
 486    //Histograms
 487    TH1I* hstat=(TH1I*)fmyOutput->FindObject("hstat");
 488    TH1F* hPtJetTrks=(TH1F*)fmyOutput->FindObject("hPtJetTrks");
 489    TH1F* hPhiJetTrks=(TH1F*)fmyOutput->FindObject("hPhiJetTrks");
 490    TH1F* hEtaJetTrks=(TH1F*)fmyOutput->FindObject("hEtaJetTrks");
 491    TH1F* hEjetTrks=(TH1F*)fmyOutput->FindObject("hEjetTrks");
 492    TH1F* hPtJet=(TH1F*)fmyOutput->FindObject("hPtJet");
 493    TH1F* hPhiJet=(TH1F*)fmyOutput->FindObject("hPhiJet");
 494    TH1F* hEtaJet=(TH1F*)fmyOutput->FindObject("hEtaJet");
 495    TH1F* hEjet=(TH1F*)fmyOutput->FindObject("hEjet");
 496    TH1F* hNtrArr=(TH1F*)fmyOutput->FindObject("hNtrArr");
 497    TH1F* hNJetPerEv=(TH1F*)fmyOutput->FindObject("hNJetPerEv");
 498    TH1F* hdeltaRJetTracks=(TH1F*)fmyOutput->FindObject("hdeltaRJetTracks");
 499    TH1F* hNDPerEvNoJet=(TH1F*)fmyOutput->FindObject("hNDPerEvNoJet");
 500    TH1F* hptDPerEvNoJet=(TH1F*)fmyOutput->FindObject("hptDPerEvNoJet");
 501    TH1F* hNJetPerEvNoD=(TH1F*)fmyOutput->FindObject("hNJetPerEvNoD");
 502    TH1F* hPtJetPerEvNoD=(TH1F*)fmyOutput->FindObject("hPtJetPerEvNoD");
 503    THnSparseF* hnspDstandalone=(THnSparseF*)fmyOutput->FindObject("hsDstandalone");
 504
 505    hstat->Fill(0);
 506
 507    // fix for temporary bug in ESDfilter
 508    // the AODs with null vertex pointer didn't pass the PhysSel
 509    if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
 510
 511    //Event selection
 512    //Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);
 513    TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();
 514    if(!iseventselected) return;
 515
 516    hstat->Fill(1);
 517
 518    //retrieve charm candidates selected
 519    Int_t candidates = fCandidateArray->GetEntriesFast();
 520
 521    //trigger on jets
 522
 523    Int_t njets=GetJetContainer()->GetNJets();
 524    if(njets == 0) {
 525       hstat->Fill(6, candidates);
 526       hNDPerEvNoJet->Fill(candidates);
 527       for(Int_t iD=0;iD<candidates;iD++){
 528          AliVParticle* cand=(AliVParticle*)fCandidateArray->At(iD);
 529          if(!cand) continue;
 530          hptDPerEvNoJet->Fill(cand->Pt());
 531
 532       }
 533       return;
 534
 535    }
 536
 537    //loop on candidates standalone (checking the candidates are there and their phi-eta distributions)
 538
 539    for(Int_t ic = 0; ic < candidates; ic++) {
 540
 541       // D* candidates
 542       AliVParticle* charm=0x0;
 543       charm=(AliVParticle*)fCandidateArray->At(ic);
 544       if(!charm) continue;
 545       hstat->Fill(2);
 546
 547       Double_t candsparse[4]={charm->Eta(), charm->Phi(), charm->Pt(), 0};
 548
 549       if(fCandidateType==kDstartoKpipi) {
 550          AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)charm;
 551          Double_t deltamass= dstar->DeltaInvMass();
 552          candsparse[3]=deltamass;
 553          hnspDstandalone->Fill(candsparse);
 554       }
 555       if(fCandidateType==kD0toKpi){
 556          AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)charm;
 557          Int_t isselected=fCuts->IsSelected(dzero,AliRDHFCuts::kAll,aodEvent);
 558
 559          Double_t masses[2];
 560          //Int_t pdgdaughtersD0[2]={211,321};//pi,K
 561          //Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
 562
 563          masses[0]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
 564          masses[1]=dzero->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
 565          if(isselected==1 || isselected==3) {
 566             candsparse[3]=masses[0];
 567             hnspDstandalone->Fill(candsparse);
 568          }
 569          if(isselected>=2){
 570             candsparse[3]=masses[1];
 571             hnspDstandalone->Fill(candsparse);
 572
 573          }
 574       }
 575    }
 576
 577    // we start with jets
 578    Double_t ejet   = 0;
 579    Double_t phiJet = 0;
 580    Double_t etaJet = 0;
 581    Double_t ptjet = 0;
 582    Double_t leadingJet =0;
 583
 584    Int_t ntrarr=trackArr->GetEntriesFast();
 585    hNtrArr->Fill(ntrarr);
 586
 587    for(Int_t i=0;i<ntrarr;i++){
 588       AliVTrack *jtrack=static_cast<AliVTrack*>(trackArr->At(i));
 589       //reusing histograms
 590       hPtJetTrks->Fill(jtrack->Pt());
 591       hPhiJetTrks->Fill(jtrack->Phi());
 592       hEtaJetTrks->Fill(jtrack->Eta());
 593       hEjetTrks->Fill(jtrack->E());
 594    }
 595
 596
 597    //option to use only the leading jet
 598    if(fLeadingJetOnly){
 599       for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) {
 600          AliEmcalJet* jetL = (AliEmcalJet*)GetJetFromArray(iJetsL);
 601          ptjet   = jetL->Pt();
 602          if(ptjet>leadingJet ) leadingJet = ptjet;
 603
 604       }
 605    }
 606
 607    Int_t cntjet=0;
 608    //loop over jets and consider only the leading jet in the event
 609    for (Int_t iJets = 0; iJets<njets; iJets++) {
 610       //Printf("Jet N %d",iJets);
 611       AliEmcalJet* jet = (AliEmcalJet*)GetJetFromArray(iJets);
 612       //Printf("Corr task Accept Jet");
 613       if(!AcceptJet(jet)) {
 614          hstat->Fill(5);
 615          continue;
 616       }
 617
 618       //jets variables
 619       ejet   = jet->E();
 620       phiJet = jet->Phi();
 621       etaJet = jet->Eta();
 622       ptjet = jet->Pt();
 623
 624       // choose the leading jet
 625       if(fLeadingJetOnly && (ejet<leadingJet)) continue;
 626       //used for normalization
 627       hstat->Fill(3);
 628       cntjet++;
 629       // fill energy, eta and phi of the jet
 630       hEjet   ->Fill(ejet);
 631       hPhiJet ->Fill(phiJet);
 632       hEtaJet ->Fill(etaJet);
 633       hPtJet  ->Fill(ptjet);
 634
 635       //loop on jet particles
 636       Int_t ntrjet=  jet->GetNumberOfTracks();
 637       for(Int_t itrk=0;itrk<ntrjet;itrk++){
 638
 639          AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,trackArr);
 640          hdeltaRJetTracks->Fill(DeltaR(jet,jetTrk));
 641
 642       }//end loop on jet tracks
 643
 644       if(candidates==0){
 645          hstat->Fill(7);
 646          hPtJetPerEvNoD->Fill(jet->Pt());
 647       }
 648
 649       //Printf("N candidates %d ", candidates);
 650       for(Int_t ic = 0; ic < candidates; ic++) {
 651
 652          // D* candidates
 653          AliVParticle* charm=0x0;
 654          charm=(AliVParticle*)fCandidateArray->At(ic);
 655          if(!charm) continue;
 656
 657
 658          FlagFlavour(charm, jet);
 659          if (jet->TestFlavourTag(AliEmcalJet::kDStar)) hstat->Fill(4);
 660
 661          FillHistogramsRecoJetCorr(charm, jet);
 662
 663       }
 664       //retrieve side band background candidates for Dstar
 665       Int_t nsbcand = 0; if(fSideBandArray) nsbcand=fSideBandArray->GetEntriesFast();
 666
 667       for(Int_t ib=0;ib<nsbcand;ib++){
 668          if(fCandidateType==kDstartoKpipi && !fUseMCInfo){
 669             AliAODRecoCascadeHF *sbcand=(AliAODRecoCascadeHF*)fSideBandArray->At(ib);
 670             if(!sbcand) continue;
 671             SideBandBackground(sbcand,jet);
 672          }
 673          if(fUseMCInfo){
 674             AliAODRecoDecayHF* charmbg = 0x0;
 675             charmbg=(AliAODRecoDecayHF*)fCandidateArray->At(ib);
 676             if(!charmbg) continue;
 677             MCBackground(charmbg,jet);
 678          }
 679       }
 680    } // end of jet loop
 681
 682    hNJetPerEv->Fill(cntjet);
 683    if(candidates==0) hNJetPerEvNoD->Fill(cntjet);
 684    PostData(1,fmyOutput);
 685    PostData(3,fmyOutputF);
 686 }
 687
 688 //_______________________________________________________________________________
 689
 690 void AliAnalysisTaskFlavourFilterAndJetCorrelations::Terminate(Option_t*)
 691 {
 692    // The Terminate() function is the last function to be called during
 693    // a query. It always runs on the client, it can be used to present
 694    // the results graphically or save the results to file.
 695
 696    Info("Terminate"," terminate");
 697    AliAnalysisTaskSE::Terminate();
 698
 699    fmyOutput = dynamic_cast<TList*> (GetOutputData(1));
 700    if (!fmyOutput) {
 701       printf("ERROR: fmyOutput not available\n");
 702       return;
 703    }
 704 }
 705
 706 //_______________________________________________________________________________
 707
 708 void  AliAnalysisTaskFlavourFilterAndJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){
 709    Float_t mass=0;
 710    Int_t abspdg=TMath::Abs(pdg);
 711
 712    mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();
 713    // compute the Delta mass for the D*
 714    if(fCandidateType==kDstartoKpipi){
 715       Float_t mass1=0;
 716       mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 717       mass = mass-mass1;
 718    }
 719
 720    fMinMass = mass-range;
 721    fMaxMass = mass+range;
 722
 723    AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));
 724    if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");
 725 }
 726
 727 //_______________________________________________________________________________
 728
 729 void  AliAnalysisTaskFlavourFilterAndJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){
 730    if(uplimit>lowlimit)
 731    {
 732       fMinMass = lowlimit;
 733       fMaxMass = uplimit;
 734    }
 735    else{
 736       printf("Error! Lower limit larger than upper limit!\n");
 737       fMinMass = uplimit - uplimit*0.2;
 738       fMaxMass = uplimit;
 739    }
 740 }
 741
 742 //_______________________________________________________________________________
 743
 744 Bool_t AliAnalysisTaskFlavourFilterAndJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){
 745    if(nptbins>30) {
 746       AliInfo("Maximum number of bins allowed is 30!");
 747       return kFALSE;
 748    }
 749    if(!width) return kFALSE;
 750    for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];
 751    return kTRUE;
 752 }
 753
 754 //_______________________________________________________________________________
 755
 756 Double_t AliAnalysisTaskFlavourFilterAndJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet) const{
 757    if(!part || !jet){
 758       printf("Particle or jet do not exist!\n");
 759       return -99;
 760    }
 761    Double_t p[3],pj[3];
 762    Bool_t okpp=part->PxPyPz(p);
 763    Bool_t okpj=jet->PxPyPz(pj);
 764    if(!okpp || !okpj){
 765       printf("Problems getting momenta\n");
 766       return -999;
 767    }
 768    Double_t pjet=jet->P();
 769    Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet*pjet);
 770    return z;
 771 }
 772
 773 //_______________________________________________________________________________
 774
 775 Bool_t  AliAnalysisTaskFlavourFilterAndJetCorrelations::DefineHistoForAnalysis(){
 776
 777    // Statistics
 778    TH1I* hstat=new TH1I("hstat","Statistics",8,-0.5,7.5);
 779    hstat->GetXaxis()->SetBinLabel(1,"N ev anal");
 780    hstat->GetXaxis()->SetBinLabel(2,"N ev sel");
 781    hstat->GetXaxis()->SetBinLabel(3,"N cand sel & jet");
 782    hstat->GetXaxis()->SetBinLabel(4,"N jets");
 783    hstat->GetXaxis()->SetBinLabel(5,"N cand in jet");
 784    hstat->GetXaxis()->SetBinLabel(6,"N jet rej");
 785    hstat->GetXaxis()->SetBinLabel(7,"N cand sel & !jet");
 786    hstat->GetXaxis()->SetBinLabel(8,"N jets & !D");
 787    hstat->SetNdivisions(1);
 788    fmyOutput->Add(hstat);
 789
 790    const Int_t nbinsmass=200;
 791    const Int_t nbinsptjet=500;
 792    const Int_t nbinsptD=100;
 793    const Int_t nbinsz=100;
 794    const Int_t nbinsphi=200;
 795    const Int_t nbinseta=100;
 796
 797    const Float_t ptjetlims[2]={0.,200.};
 798    const Float_t ptDlims[2]={0.,50.};
 799    const Float_t zlims[2]={0.,1.2};
 800    const Float_t philims[2]={0.,6.3};
 801    const Float_t etalims[2]={-1.5,1.5};
 802
 803    if(fCandidateType==kDstartoKpipi)
 804    {
 805
 806       TH2F* hDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);
 807       hDiffSideBand->SetStats(kTRUE);
 808       hDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
 809       hDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
 810       hDiffSideBand->Sumw2();
 811       fmyOutput->Add(hDiffSideBand);
 812
 813
 814       TH1F* hPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);
 815       hPtPion->SetStats(kTRUE);
 816       hPtPion->GetXaxis()->SetTitle("GeV/c");
 817       hPtPion->GetYaxis()->SetTitle("Entries");
 818       hPtPion->Sumw2();
 819       fmyOutput->Add(hPtPion);
 820
 821    }
 822    // D related histograms
 823    TH1F *hNDPerEvNoJet=new TH1F("hNDPerEvNoJet","Number of candidates per event with no jets; N candidate/ev with no jet", 20, 0., 20.);
 824    hNDPerEvNoJet->Sumw2();
 825    fmyOutput->Add(hNDPerEvNoJet);
 826
 827    TH1F *hptDPerEvNoJet=new TH1F("hptDPerEvNoJet","pt distribution of candidates per events with no jets; p_{t}^{D} (GeV/c)",nbinsptD, ptDlims[0],ptDlims[1]);
 828    hptDPerEvNoJet->Sumw2();
 829    fmyOutput->Add(hptDPerEvNoJet);
 830
 831    const Int_t    nAxisD=4;
 832    const Int_t    nbinsSparseD[nAxisD]={nbinsphi,nbinseta,nbinsptD,nbinsmass};
 833    const Double_t minSparseD[nAxisD]  ={philims[0],etalims[0],ptDlims[0],fMinMass};
 834    const Double_t maxSparseD[nAxisD]  ={philims[1],etalims[1],ptDlims[1],fMaxMass};
 835    THnSparseF *hsDstandalone=new THnSparseF("hsDstandalone","#phi, #eta, p_{T}^{D}, and mass", nAxisD, nbinsSparseD, minSparseD, maxSparseD);
 836    hsDstandalone->Sumw2();
 837
 838    fmyOutput->Add(hsDstandalone);
 839
 840    // jet related fistograms
 841
 842    TH1F* hEjetTrks      = new TH1F("hEjetTrks",  "Jet tracks energy distribution;Energy (GeV)",500,0,200);
 843    hEjetTrks->Sumw2();
 844    TH1F* hPhiJetTrks    = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)",  nbinsphi,philims[0],philims[1]);
 845    hPhiJetTrks->Sumw2();
 846    TH1F* hEtaJetTrks    = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta",  nbinseta,etalims[0],etalims[1]);
 847    hEtaJetTrks->Sumw2();
 848    TH1F* hPtJetTrks     = new TH1F("hPtJetTrks",  "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
 849    hPtJetTrks->Sumw2();
 850
 851    TH1F* hEjet      = new TH1F("hEjet",  "Jet energy distribution;Energy (GeV)",500,0,200);
 852    hEjet->Sumw2();
 853    TH1F* hPhiJet    = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)",  nbinsphi,philims[0],philims[1]);
 854    hPhiJet->Sumw2();
 855    TH1F* hEtaJet    = new TH1F("hEtaJet","Jet #eta distribution; #eta", nbinseta,etalims[0],etalims[1]);
 856    hEtaJet->Sumw2();
 857    TH1F* hPtJet      = new TH1F("hPtJet",  "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
 858    hPtJet->Sumw2();
 859
 860    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]);
 861    hPtJetWithD->Sumw2();
 862    //for the MC this histogram is filled with the real background
 863    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]);
 864    hPtJetWithDsb->Sumw2();
 865
 866    TH1F* hdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.);
 867    hdeltaRJetTracks->Sumw2();
 868
 869    TH1F* hNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000);
 870    hNtrArr->Sumw2();
 871
 872    TH1F *hNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);
 873    hNJetPerEv->Sumw2();
 874
 875    TH1F *hNJetPerEvNoD=new TH1F("hNJetPerEvNoD","Number of jets per event with no D; number of jets/ev with no D",10,-0.5,9.5);
 876    hNJetPerEvNoD->Sumw2();
 877
 878    TH1F *hPtJetPerEvNoD=new TH1F("hPtJetPerEvNoD","pt distribution of jets per event with no D; p_{T}^{jet} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);
 879    hPtJetPerEvNoD->Sumw2();
 880
 881    fmyOutput->Add(hEjetTrks);
 882    fmyOutput->Add(hPhiJetTrks);
 883    fmyOutput->Add(hEtaJetTrks);
 884    fmyOutput->Add(hPtJetTrks);
 885    fmyOutput->Add(hEjet);
 886    fmyOutput->Add(hPhiJet);
 887    fmyOutput->Add(hEtaJet);
 888    fmyOutput->Add(hPtJet);
 889    fmyOutput->Add(hPtJetWithD);
 890    fmyOutput->Add(hPtJetWithDsb);
 891    fmyOutput->Add(hdeltaRJetTracks);
 892    fmyOutput->Add(hNtrArr);
 893    fmyOutput->Add(hNJetPerEv);
 894    fmyOutput->Add(hNJetPerEvNoD);
 895    fmyOutput->Add(hPtJetPerEvNoD);
 896
 897    TH1F* hDeltaRD=new TH1F("hDeltaRD","#Delta R distribution of D candidates selected;#Delta R",200, 0.,10.);
 898    hDeltaRD->Sumw2();
 899    fmyOutput->Add(hDeltaRD);
 900
 901    //background (side bands for the Dstar and like sign for D0)
 902    fJetRadius=GetJetContainer(0)->GetJetRadius();
 903    TH2F* hInvMassptD = new TH2F("hInvMassptD",Form("D (Delta R < %.1f) invariant mass distribution p_{T}^{j} > threshold",fJetRadius),nbinsmass,fMinMass,fMaxMass,nbinsptD,ptDlims[0],ptDlims[1]);
 904    hInvMassptD->SetStats(kTRUE);
 905    hInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
 906    hInvMassptD->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
 907    hInvMassptD->Sumw2();
 908
 909    fmyOutput->Add(hInvMassptD);
 910
 911    if(fUseMCInfo){
 912       TH2F* hInvMassptDbg = 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]);
 913       hInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));
 914       hInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
 915       hInvMassptDbg->Sumw2();
 916       fmyOutput->Add(hInvMassptDbg);
 917
 918    }
 919    Double_t pi=TMath::Pi(), philims2[2];
 920    philims2[0]=-pi*1./2.;
 921    philims2[1]=pi*3./2.;
 922    const Int_t nAxis=6;
 923    const Int_t nbinsSparse[nAxis]={nbinsz,nbinsphi,nbinsptjet,nbinsptD,nbinsmass,2};
 924    const Double_t minSparse[nAxis]={zlims[0],philims2[0],ptjetlims[0],ptDlims[0],fMinMass,-0.5};
 925    const Double_t maxSparse[nAxis]={zlims[1],philims2[1],ptjetlims[1],ptDlims[1],fMaxMass, 1.5};
 926    THnSparseF *hsDphiz=new THnSparseF("hsDphiz","Z and #Delta#phi vs p_{T}^{jet}, p_{T}^{D}, and mass", nAxis, nbinsSparse, minSparse, maxSparse);
 927    hsDphiz->Sumw2();
 928
 929    fmyOutput->Add(hsDphiz);
 930
 931    PostData(1, fmyOutput);
 932
 933    return kTRUE;
 934 }
 935
 936 //_______________________________________________________________________________
 937
 938 Bool_t AliAnalysisTaskFlavourFilterAndJetCorrelations::DefineHistoForFiltering()
 939 {
 940    //
 941    // AliAnalysisTaskSEDmesonsFilterCJ::DefineHistoForAnalysis
 942    //
 943
 944    // Statistics
 945    TH1I* hstat = new TH1I("hstatF","Statistics",5,-0.5,4.5);
 946    hstat->GetXaxis()->SetBinLabel(1, "N ev anal");
 947    hstat->GetXaxis()->SetBinLabel(2, "N ev sel");
 948    if(fUseReco) hstat->GetXaxis()->SetBinLabel(3, "N cand");
 949    else hstat->GetXaxis()->SetBinLabel(3, "N Gen D");
 950    hstat->GetXaxis()->SetBinLabel(4, "N cand sel cuts");
 951    if (fCandidateType==kDstartoKpipi) hstat->GetXaxis()->SetBinLabel(5, "N side band cand");
 952    hstat->SetNdivisions(1);
 953    fmyOutputF->Add(hstat);
 954
 955    TH1F* hnCandEv=new TH1F("hnCandEv", "Number of candidates per event (after cuts);# cand/ev", 100, 0.,100.);
 956    fmyOutputF->Add(hnCandEv);
 957
 958    // Invariant mass related histograms
 959    const Int_t nbinsmass = 200;
 960    TH2F *hInvMass = new TH2F("hInvMassptDF", "D invariant mass distribution", nbinsmass, fMinMass, fMaxMass, 100, 0., 50.);
 961    hInvMass->SetStats(kTRUE);
 962    hInvMass->GetXaxis()->SetTitle("mass (GeV/c)");
 963    hInvMass->GetYaxis()->SetTitle("p_{T} (GeV/c)");
 964    fmyOutputF->Add(hInvMass);
 965
 966    if (fCandidateType==kDstartoKpipi) {
 967       TH1F* hnSBCandEv=new TH1F("hnSBCandEv", "Number of side bands candidates per event (after cuts);# cand/ev", 100, 0.,100.);
 968       fmyOutputF->Add(hnSBCandEv);
 969
 970       TH1F* hPtPion = new TH1F("hPtPion", "Primary pions candidates pt", 500, 0., 10.);
 971       hPtPion->SetStats(kTRUE);
 972       hPtPion->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 973       hPtPion->GetYaxis()->SetTitle("entries");
 974       fmyOutputF->Add(hPtPion);
 975    }
 976    PostData(3, fmyOutputF);
 977    return kTRUE;
 978 }
 979
 980 //_______________________________________________________________________________
 981
 982 void AliAnalysisTaskFlavourFilterAndJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet){
 983
 984    Double_t ptD=candidate->Pt();
 985    Double_t ptjet=jet->Pt();
 986    Double_t deltaR=DeltaR(candidate,jet);
 987    Double_t deltaphi = jet->Phi()-candidate->Phi();
 988    if(deltaphi<=-(TMath::Pi())/2) deltaphi = deltaphi+2*(TMath::Pi());
 989    if(deltaphi>(3*(TMath::Pi()))/2) deltaphi = deltaphi-2*(TMath::Pi());
 990    Double_t z=Z(candidate,jet);
 991
 992    TH1F* hDeltaRD=(TH1F*)fmyOutput->FindObject("hDeltaRD");
 993    hDeltaRD->Fill(deltaR);
 994    if(fUseReco){
 995       if(fCandidateType==kD0toKpi) {
 996          AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;
 997          FillHistogramsD0JetCorr(dzero,deltaphi,z,ptD,ptjet,deltaR, AODEvent());
 998
 999       }
1000
1001       if(fCandidateType==kDstartoKpipi) {
1002          AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;
1003          FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,ptjet,deltaR);
1004
1005       }
1006    } else{ //generated level
1007       //AliInfo("Non reco");
1008       FillHistogramsMCGenDJetCorr(deltaphi,z,ptD,ptjet,deltaR);
1009    }
1010 }
1011
1012 //_______________________________________________________________________________
1013
1014 void AliAnalysisTaskFlavourFilterAndJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj,Double_t deltaR, AliAODEvent* aodEvent){
1015
1016
1017    Double_t masses[2]={0.,0.};
1018    Int_t pdgdaughtersD0[2]={211,321};//pi,K
1019    Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1020
1021    masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1022    masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1023
1024    TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");
1025    THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");
1026    Double_t point[5]={z,dPhi,ptj,ptD,masses[0]};
1027
1028    Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
1029    if(isselected==1 || isselected==3) {
1030
1031       if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[0],ptD);
1032
1033       FillMassHistograms(masses[0], ptD, deltaR);
1034       hsDphiz->Fill(point,1.);
1035    }
1036    if(isselected>=2) {
1037       if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[1],ptD);
1038
1039       FillMassHistograms(masses[1], ptD, deltaR);
1040       point[4]=masses[1];
1041       hsDphiz->Fill(point,1.);
1042    }
1043
1044 }
1045
1046 //_______________________________________________________________________________
1047
1048 void AliAnalysisTaskFlavourFilterAndJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj, Double_t deltaR){
1049   //dPhi and z not used at the moment,but will be (re)added
1050
1051    AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();
1052    Double_t deltamass= dstar->DeltaInvMass();
1053    //Double_t massD0= dstar->InvMassD0();
1054
1055
1056    TH1F* hPtPion=(TH1F*)fmyOutput->FindObject("hPtPion");
1057    hPtPion->Fill(softpi->Pt());
1058
1059    TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");
1060    THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");
1061    Int_t isSB=IsDzeroSideBand(dstar);
1062
1063    Double_t point[6]={z,dPhi,ptj,ptD,deltamass,isSB};
1064
1065    if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,deltamass,ptD);
1066
1067    FillMassHistograms(deltamass, ptD, deltaR);
1068    hsDphiz->Fill(point,1.);
1069
1070 }
1071
1072 //_______________________________________________________________________________
1073
1074 void AliAnalysisTaskFlavourFilterAndJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t z,Double_t ptD,Double_t ptjet,Double_t deltaR){
1075
1076    Double_t pdgmass=0;
1077    TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");
1078    THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");
1079    Double_t point[6]={z,dPhi,ptjet,ptD,pdgmass,0};
1080
1081    if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1082    if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();
1083    point[4]=pdgmass;
1084
1085    if(deltaR<fJetRadius) {
1086      hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet
1087      hsDphiz->Fill(point,1.);
1088    }
1089
1090 }
1091
1092 //_______________________________________________________________________________
1093
1094 void AliAnalysisTaskFlavourFilterAndJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD, Double_t deltaR){
1095
1096    if(deltaR<fJetRadius) {
1097       TH2F* hInvMassptD=(TH2F*)fmyOutput->FindObject("hInvMassptD");
1098       hInvMassptD->Fill(mass,ptD);
1099    }
1100 }
1101
1102 //________________________________________________________________________________
1103
1104 void AliAnalysisTaskFlavourFilterAndJetCorrelations::FlagFlavour(AliVParticle *charm, AliEmcalJet *jet){
1105    Double_t deltaR=DeltaR(charm, jet);
1106    AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;
1107    if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;
1108    if (deltaR<fJetRadius) jet->AddFlavourTag(tag);
1109
1110    return;
1111
1112 }
1113
1114 //_______________________________________________________________________________
1115
1116 void AliAnalysisTaskFlavourFilterAndJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){
1117
1118    //  D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
1119    // (expected detector resolution) on the left and right frm the D0 mass. Each band
1120    //  has a width of ~5 sigmas. Two band needed  for opening angle considerations
1121    TH2F* hDiffSideBand=(TH2F*)fmyOutput->FindObject("hDiffSideBand");
1122    TH3F* hPtJetWithDsb=(TH3F*)fmyOutput->FindObject("hPtJetWithDsb");
1123
1124    Double_t deltaM=candDstar->DeltaInvMass();
1125    //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);
1126    Double_t ptD=candDstar->Pt();
1127    Double_t ptjet=jet->Pt();
1128    Double_t dPhi=jet->Phi()-candDstar->Phi();
1129    Double_t deltaR=DeltaR(candDstar,jet);
1130    if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1131    if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1132
1133    Int_t isSideBand=IsDzeroSideBand(candDstar);
1134    //fill the background histogram with the side bands or when looking at MC with the real background
1135    if(isSideBand==1){
1136       hDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background
1137       //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);
1138
1139       if(deltaR<fJetRadius){  // evaluate in the near side
1140          //hzptDB->Fill(Z(candDstar,jet),deltaM,ptD);
1141          hPtJetWithDsb->Fill(ptjet,deltaM,ptD);
1142       }
1143    }
1144 }
1145
1146 //_______________________________________________________________________________
1147
1148 void AliAnalysisTaskFlavourFilterAndJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg, AliEmcalJet *jet){
1149
1150    //need mass, deltaR, pt jet, ptD
1151    //two cases: D0 and Dstar
1152
1153    Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());
1154    if(!isselected) return;
1155
1156    Double_t ptD=candbg->Pt();
1157    Double_t ptjet=jet->Pt();
1158    Double_t deltaR=DeltaR(candbg,jet);
1159
1160    TH2F* hInvMassptDbg=(TH2F*)fmyOutput->FindObject("hInvMassptDbg");
1161    TH3F* hPtJetWithDsb=(TH3F*)fmyOutput->FindObject("hPtJetWithDsb");
1162
1163
1164    if(fCandidateType==kDstartoKpipi){
1165       AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;
1166       Double_t deltaM=dstarbg->DeltaInvMass();
1167       hInvMassptDbg->Fill(deltaM,ptD);
1168       if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,deltaM,ptD);
1169    }
1170
1171    if(fCandidateType==kD0toKpi){
1172       Double_t masses[2]={0.,0.};
1173       Int_t pdgdaughtersD0[2]={211,321};//pi,K
1174       Int_t pdgdaughtersD0bar[2]={321,211};//K,pi
1175
1176       masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0
1177       masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar
1178
1179
1180       if(isselected==1 || isselected==3) {
1181          if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[0],ptD);
1182          hInvMassptDbg->Fill(masses[0],ptD);
1183       }
1184       if(isselected>=2) {
1185          if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[1],ptD);
1186          hInvMassptDbg->Fill(masses[1],ptD);
1187       }
1188
1189
1190    }
1191 }
1192
1193 //_______________________________________________________________________________
1194
1195 Float_t AliAnalysisTaskFlavourFilterAndJetCorrelations::DeltaR(AliVParticle *p1, AliVParticle *p2) const {
1196    //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)
1197
1198    if(!p1 || !p2) return -1;
1199    Double_t phi1=p1->Phi(),eta1=p1->Eta();
1200    Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;
1201
1202    Double_t dPhi=phi1-phi2;
1203    if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
1204    if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
1205
1206    Double_t dEta=eta1-eta2;
1207    Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );
1208    return deltaR;
1209
1210 }
1211
1212 //_______________________________________________________________________________
1213
1214 Int_t AliAnalysisTaskFlavourFilterAndJetCorrelations::IsDzeroSideBand(AliAODRecoCascadeHF *candDstar){
1215
1216    Double_t ptD=candDstar->Pt();
1217    Int_t bin = fCuts->PtBin(ptD);
1218    if (bin < 0){
1219       // /PWGHF/vertexingHF/AliRDHFCuts::PtBin(Double_t) const may return values below zero depending on config.
1220       bin = 9999; // void result code for coverity (bin later in the code non-zero) - will coverity pick up on this?
1221       return -1;  // out of bounds
1222    }
1223
1224    Double_t invM=candDstar->InvMassD0();
1225    Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1226
1227    Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];
1228    Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];
1229
1230    if((invM>=sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)) return 1;
1231    else return 0;
1232
1233 }