PWGJE/FlavourJetTasks/AliAnalysisTaskFlavourJetCorrelations.cxx

   1 /**************************************************************************\r
   2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
   3 *                                                                        *\r
   4 * Author: The ALICE Off-line Project.                                    *\r
   5 * Contributors are mentioned in the code where appropriate.              *\r
   6 *                                                                        *\r
   7 * Permission to use, copy, modify and distribute this software and its   *\r
   8 * documentation strictly for non-commercial purposes is hereby granted   *\r
   9 * without fee, provided that the above copyright notice appears in all   *\r
  10 * copies and that both the copyright notice and this permission notice   *\r
  11 * appear in the supporting documentation. The authors make no claims     *\r
  12 * about the suitability of this software for any purpose. It is          *\r
  13 * provided "as is" without express or implied warranty.                  *\r
  14 **************************************************************************/\r
  15 //\r
  16 //\r
  17 //  Analysis Taks for reconstructed particle correlation \r
  18 //  (first implementation done for D mesons) with jets \r
  19 //  (use the so called Emcal framework)\r
  20 //\r
  21 //-----------------------------------------------------------------------\r
  22 // Authors:\r
  23 // C. Bianchin (Utrecht University) chiara.bianchin@cern.ch\r
  24 // A. Grelli (Utrecht University) a.grelli@uu.nl\r
  25 // X. Zhang (LBNL)  XMZhang@lbl.gov\r
  26 //-----------------------------------------------------------------------\r
  27 \r
  28 #include <TDatabasePDG.h>\r
  29 #include <TParticle.h>\r
  30 #include "TROOT.h"\r
  31 #include <TH3F.h>\r
  32 #include <THnSparse.h>\r
  33 \r
  34 #include "AliAnalysisTaskFlavourJetCorrelations.h"\r
  35 #include "AliAODMCHeader.h"\r
  36 #include "AliAODHandler.h"\r
  37 #include "AliAnalysisManager.h"\r
  38 #include "AliLog.h"\r
  39 #include "AliEmcalJet.h"\r
  40 #include "AliJetContainer.h"\r
  41 #include "AliAODRecoDecay.h"\r
  42 #include "AliAODRecoCascadeHF.h"\r
  43 #include "AliAODRecoDecayHF2Prong.h"\r
  44 #include "AliESDtrack.h"\r
  45 #include "AliAODMCParticle.h"\r
  46 #include "AliPicoTrack.h"\r
  47 #include "AliRDHFCutsD0toKpi.h"\r
  48 #include "AliRDHFCutsDStartoKpipi.h"\r
  49 \r
  50 ClassImp(AliAnalysisTaskFlavourJetCorrelations)\r
  51 \r
  52 \r
  53 //_______________________________________________________________________________\r
  54 \r
  55 AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations() :\r
  56 AliAnalysisTaskEmcalJet("",kFALSE),\r
  57 fUseMCInfo(kTRUE), \r
  58 fUseReco(kTRUE),\r
  59 fCandidateType(),\r
  60 fPDGmother(),\r
  61 fNProngs(),\r
  62 fPDGdaughters(),\r
  63 fBranchName(),\r
  64 fmyOutput(0),\r
  65 fCuts(0),\r
  66 fMinMass(),\r
  67 fMaxMass(),  \r
  68 fJetArrName(0),\r
  69 fCandArrName(0),\r
  70 fLeadingJetOnly(kFALSE),\r
  71 fJetRadius(0)\r
  72 {\r
  73    //\r
  74    // Default ctor\r
  75    //\r
  76    //SetMakeGeneralHistograms(kTRUE);\r
  77    \r
  78 }\r
  79 \r
  80 //_______________________________________________________________________________\r
  81 \r
  82 AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations(const Char_t* name, AliRDHFCuts* cuts,ECandidateType candtype) :\r
  83 AliAnalysisTaskEmcalJet(name,kFALSE),\r
  84 fUseMCInfo(kTRUE),\r
  85 fUseReco(kTRUE),  \r
  86 fCandidateType(),\r
  87 fPDGmother(),\r
  88 fNProngs(),\r
  89 fPDGdaughters(),\r
  90 fBranchName(),\r
  91 fmyOutput(0),\r
  92 fCuts(0),\r
  93 fMinMass(),\r
  94 fMaxMass(),  \r
  95 fJetArrName(0),\r
  96 fCandArrName(0),\r
  97 fLeadingJetOnly(kFALSE),\r
  98 fJetRadius(0)\r
  99 {\r
 100    //\r
 101    // Constructor. Initialization of Inputs and Outputs\r
 102    //\r
 103    \r
 104    Info("AliAnalysisTaskFlavourJetCorrelations","Calling Constructor");\r
 105    fCuts=cuts;\r
 106    fCandidateType=candtype;\r
 107    const Int_t nptbins=fCuts->GetNPtBins();\r
 108    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};\r
 109    \r
 110    switch(fCandidateType){\r
 111    case 0:\r
 112       fPDGmother=421;\r
 113       fNProngs=2;\r
 114       fPDGdaughters[0]=211;//pi \r
 115       fPDGdaughters[1]=321;//K\r
 116       fPDGdaughters[2]=0; //empty\r
 117       fPDGdaughters[3]=0; //empty\r
 118       fBranchName="D0toKpi";\r
 119       fCandArrName="D0";\r
 120       break;\r
 121    case 1: \r
 122       fPDGmother=413;\r
 123       fNProngs=3;\r
 124       fPDGdaughters[1]=211;//pi soft\r
 125       fPDGdaughters[0]=421;//D0\r
 126       fPDGdaughters[2]=211;//pi fromD0\r
 127       fPDGdaughters[3]=321; // K from D0\r
 128       fBranchName="Dstar";\r
 129       fCandArrName="Dstar";\r
 130       \r
 131       if(nptbins<=13){\r
 132          for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]= defaultSigmaD013[ipt];\r
 133       } else {\r
 134          AliFatal(Form("Default sigma D0 not enough for %d pt bins, use SetSigmaD0ForDStar to set them",nptbins));\r
 135       }\r
 136       break;\r
 137    default:\r
 138       printf("%d not accepted!!\n",fCandidateType);\r
 139       break;\r
 140    }\r
 141    \r
 142    if(fCandidateType==kD0toKpi)SetMassLimits(0.15,fPDGmother);\r
 143    if(fCandidateType==kDstartoKpipi) SetMassLimits(0.015, fPDGmother);\r
 144    \r
 145    DefineOutput(1,TList::Class()); // histos\r
 146    DefineOutput(2,AliRDHFCuts::Class()); // my cuts\r
 147    \r
 148 }\r
 149 \r
 150 //_______________________________________________________________________________\r
 151 \r
 152 AliAnalysisTaskFlavourJetCorrelations::~AliAnalysisTaskFlavourJetCorrelations() {\r
 153    //\r
 154    // destructor\r
 155    //\r
 156    \r
 157    Info("~AliAnalysisTaskFlavourJetCorrelations","Calling Destructor");  \r
 158    \r
 159    delete fmyOutput;\r
 160    delete fCuts;\r
 161    \r
 162 }\r
 163 \r
 164 //_______________________________________________________________________________\r
 165 \r
 166 void AliAnalysisTaskFlavourJetCorrelations::Init(){\r
 167    //\r
 168    // Initialization\r
 169    //\r
 170    \r
 171    if(fDebug > 1) printf("AnalysisTaskRecoJetCorrelations::Init() \n");\r
 172    switch(fCandidateType){\r
 173    case 0:\r
 174       {\r
 175          AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));\r
 176          copyfCuts->SetName("AnalysisCutsDzero");\r
 177          // Post the data\r
 178          PostData(2,copyfCuts);\r
 179       }\r
 180       break;\r
 181    case 1:\r
 182       {\r
 183          AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));\r
 184          copyfCuts->SetName("AnalysisCutsDStar");\r
 185          // Post the cuts\r
 186          PostData(2,copyfCuts);\r
 187       }\r
 188       break;\r
 189    default:\r
 190       return;\r
 191    }\r
 192    \r
 193    return;\r
 194 }\r
 195 \r
 196 //_______________________________________________________________________________\r
 197 \r
 198 void AliAnalysisTaskFlavourJetCorrelations::UserCreateOutputObjects() { \r
 199    // output \r
 200    Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());\r
 201    fmyOutput = new TList();\r
 202    fmyOutput->SetOwner();\r
 203    fmyOutput->SetName("pippo");\r
 204    // define histograms\r
 205    DefineHistoForAnalysis();\r
 206    PostData(1,fmyOutput);\r
 207    \r
 208    return;\r
 209 }\r
 210 \r
 211 //_______________________________________________________________________________\r
 212 \r
 213 void AliAnalysisTaskFlavourJetCorrelations::UserExec(Option_t *)\r
 214 {\r
 215    // user exec\r
 216    if (!fInitialized){\r
 217       AliAnalysisTaskEmcalJet::ExecOnce();\r
 218    }\r
 219    \r
 220    // Load the event\r
 221    AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);\r
 222    \r
 223    TClonesArray *arrayDStartoD0pi=0;\r
 224    TClonesArray *trackArr = 0;\r
 225    TClonesArray *candidatesArr = 0;\r
 226    TClonesArray *sbcandArr = 0;\r
 227    \r
 228    if (!aodEvent && AODEvent() && IsStandardAOD()) {\r
 229       \r
 230       // In case there is an AOD handler writing a standard AOD, use the AOD \r
 231       // event in memory rather than the input (ESD) event.    \r
 232       aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());\r
 233       \r
 234       // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)\r
 235       // have to taken from the AOD event hold by the AliAODExtension\r
 236       AliAODHandler* aodHandler = (AliAODHandler*) \r
 237       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());\r
 238       if(aodHandler->GetExtensions()) {\r
 239          AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");\r
 240          AliAODEvent *aodFromExt = ext->GetAOD();\r
 241          arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject(fBranchName.Data());\r
 242       }\r
 243    } else {\r
 244       arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject(fBranchName.Data());\r
 245    }\r
 246    \r
 247    if (!arrayDStartoD0pi) {\r
 248       AliInfo(Form("Could not find array %s, skipping the event",fBranchName.Data()));\r
 249       //  return;\r
 250    } else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));   \r
 251    \r
 252    TClonesArray* mcArray = 0x0;\r
 253    if (fUseMCInfo) {\r
 254       mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));\r
 255       if (!mcArray) {\r
 256          printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");\r
 257          return;\r
 258       }\r
 259    }\r
 260    \r
 261    //retrieve jets\r
 262    trackArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("PicoTracks"));\r
 263    //clusArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClustersCorr"));\r
 264    //jetArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fJetArrName));\r
 265    fJetRadius=GetJetContainer(0)->GetJetRadius();\r
 266    \r
 267    if(!trackArr){\r
 268       AliInfo(Form("Could not find the track array\n"));\r
 269       return;\r
 270    }\r
 271    \r
 272    \r
 273    TString arrname="fCandidateArray";\r
 274    TString candarrname=Form("%s%s%s",arrname.Data(),fCandArrName.Data(),fUseReco ? "rec" : "gen");\r
 275    candidatesArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(candarrname));\r
 276    if (!candidatesArr) {\r
 277       Printf("%s array not found",candarrname.Data());\r
 278       InputEvent()->GetList()->ls();\r
 279       return;\r
 280    }\r
 281    //Printf("ncandidates found %d",candidatesArr->GetEntriesFast());\r
 282    \r
 283    TString arrSBname="fSideBandArray";\r
 284    TString sbarrname=Form("%s%s%s",arrSBname.Data(),fCandArrName.Data(),fUseReco ? "rec" : "gen");\r
 285    sbcandArr = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(sbarrname));\r
 286    if (fCandidateType==1 && !sbcandArr) {\r
 287       Printf("%s array not found",sbarrname.Data());\r
 288       InputEvent()->GetList()->ls();\r
 289       return;\r
 290    }\r
 291    //Printf("nSBCand or Bkg found %d",sbcandArr->GetEntriesFast());\r
 292    \r
 293    \r
 294    //Histograms\r
 295    TH1I* hstat=(TH1I*)fmyOutput->FindObject("hstat");\r
 296    TH1F* hPtJetTrks=(TH1F*)fmyOutput->FindObject("hPtJetTrks");\r
 297    TH1F* hPhiJetTrks=(TH1F*)fmyOutput->FindObject("hPhiJetTrks");\r
 298    TH1F* hEtaJetTrks=(TH1F*)fmyOutput->FindObject("hEtaJetTrks");\r
 299    TH1F* hEjetTrks=(TH1F*)fmyOutput->FindObject("hEjetTrks");\r
 300    TH1F* hPtJet=(TH1F*)fmyOutput->FindObject("hPtJet");\r
 301    TH1F* hPhiJet=(TH1F*)fmyOutput->FindObject("hPhiJet");\r
 302    TH1F* hEtaJet=(TH1F*)fmyOutput->FindObject("hEtaJet");\r
 303    TH1F* hEjet=(TH1F*)fmyOutput->FindObject("hEjet");\r
 304    TH1F* hNtrArr=(TH1F*)fmyOutput->FindObject("hNtrArr");\r
 305    TH1F* hNJetPerEv=(TH1F*)fmyOutput->FindObject("hNJetPerEv");\r
 306    TH1F* hdeltaRJetTracks=((TH1F*)fmyOutput->FindObject("hdeltaRJetTracks"));\r
 307    \r
 308    hstat->Fill(0);\r
 309    \r
 310    // fix for temporary bug in ESDfilter \r
 311    // the AODs with null vertex pointer didn't pass the PhysSel\r
 312    if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;\r
 313    \r
 314    //Event selection\r
 315    Bool_t iseventselected=fCuts->IsEventSelected(aodEvent);\r
 316    TString firedTriggerClasses=((AliAODEvent*)aodEvent)->GetFiredTriggerClasses();\r
 317    if(!iseventselected) return;\r
 318    \r
 319    hstat->Fill(1);\r
 320    \r
 321    //trigger on jets\r
 322    \r
 323    Int_t njets=GetJetContainer()->GetNJets();\r
 324    if(njets == 0) return;\r
 325 \r
 326    //retrieve charm candidates selected\r
 327    Int_t candidates = candidatesArr->GetEntriesFast();\r
 328    \r
 329    // we start with jets\r
 330    Double_t ejet   = 0;\r
 331    Double_t phiJet = 0;\r
 332    Double_t etaJet = 0;\r
 333    Double_t ptjet = 0;\r
 334    Double_t leadingJet =0;\r
 335    \r
 336    Int_t ntrarr=trackArr->GetEntriesFast();\r
 337    hNtrArr->Fill(ntrarr);\r
 338    \r
 339    for(Int_t i=0;i<ntrarr;i++){\r
 340       AliVTrack *jtrack=static_cast<AliVTrack*>(trackArr->At(i));\r
 341       //reusing histograms\r
 342       hPtJetTrks->Fill(jtrack->Pt());\r
 343       hPhiJetTrks->Fill(jtrack->Phi());\r
 344       hEtaJetTrks->Fill(jtrack->Eta());\r
 345       hEjetTrks->Fill(jtrack->E());\r
 346    }\r
 347    \r
 348    \r
 349    //option to use only the leading jet\r
 350    if(fLeadingJetOnly){\r
 351       for (Int_t iJetsL = 0; iJetsL<njets; iJetsL++) {    \r
 352          AliEmcalJet* jetL = (AliEmcalJet*)GetJetFromArray(iJetsL);\r
 353          ptjet   = jetL->Pt();\r
 354          if(ptjet>leadingJet ) leadingJet = ptjet;\r
 355          \r
 356       }\r
 357    }\r
 358    \r
 359    Int_t cntjet=0;\r
 360    //loop over jets and consider only the leading jet in the event\r
 361    for (Int_t iJets = 0; iJets<njets; iJets++) {    \r
 362       //Printf("Jet N %d",iJets);\r
 363       AliEmcalJet* jet = (AliEmcalJet*)GetJetFromArray(iJets);\r
 364       //Printf("Corr task Accept Jet");\r
 365       if(!AcceptJet(jet)) {\r
 366          hstat->Fill(5);\r
 367          continue;\r
 368       }\r
 369       \r
 370       //jets variables\r
 371       ejet   = jet->E();\r
 372       phiJet = jet->Phi();\r
 373       etaJet = jet->Eta();\r
 374       ptjet = jet->Pt();\r
 375       \r
 376       // choose the leading jet\r
 377       if(fLeadingJetOnly && (ejet<leadingJet)) continue;\r
 378       //used for normalization\r
 379       hstat->Fill(3);\r
 380       cntjet++;\r
 381       // fill energy, eta and phi of the jet\r
 382       hEjet   ->Fill(ejet);\r
 383       hPhiJet ->Fill(phiJet);\r
 384       hEtaJet ->Fill(etaJet);\r
 385       hPtJet  ->Fill(ptjet);\r
 386       \r
 387       //loop on jet particles\r
 388       Int_t ntrjet=  jet->GetNumberOfTracks();    \r
 389       for(Int_t itrk=0;itrk<ntrjet;itrk++){\r
 390          \r
 391          AliPicoTrack* jetTrk=(AliPicoTrack*)jet->TrackAt(itrk,trackArr);     \r
 392          hdeltaRJetTracks->Fill(DeltaR(jet,jetTrk));\r
 393          \r
 394       }//end loop on jet tracks\r
 395       \r
 396       //Printf("N candidates %d ", candidates);\r
 397       for(Int_t ic = 0; ic < candidates; ic++) {\r
 398          \r
 399          // D* candidates\r
 400          AliVParticle* charm=0x0;\r
 401          charm=(AliVParticle*)candidatesArr->At(ic);\r
 402          if(!charm) continue;\r
 403          hstat->Fill(2);\r
 404          \r
 405          FlagFlavour(charm, jet);\r
 406          if (jet->TestFlavourTag(AliEmcalJet::kDStar)) hstat->Fill(4);\r
 407          \r
 408          FillHistogramsRecoJetCorr(charm, jet);\r
 409          \r
 410       }\r
 411       //retrieve side band background candidates for Dstar\r
 412       Int_t nsbcand = 0; if(sbcandArr) nsbcand=sbcandArr->GetEntriesFast();\r
 413       \r
 414       for(Int_t ib=0;ib<nsbcand;ib++){\r
 415          if(fCandidateType==kDstartoKpipi && !fUseMCInfo){\r
 416             AliAODRecoCascadeHF *sbcand=(AliAODRecoCascadeHF*)sbcandArr->At(ib);\r
 417             if(!sbcand) continue;\r
 418             SideBandBackground(sbcand,jet);\r
 419          }\r
 420          if(fUseMCInfo){\r
 421             AliAODRecoDecayHF* charmbg = 0x0;\r
 422             charmbg=(AliAODRecoDecayHF*)candidatesArr->At(ib);\r
 423             if(!charmbg) continue;\r
 424             MCBackground(charmbg,jet);\r
 425          }\r
 426       }\r
 427    } // end of jet loop\r
 428    \r
 429    hNJetPerEv->Fill(cntjet);\r
 430    \r
 431    PostData(1,fmyOutput);\r
 432    \r
 433 }\r
 434 \r
 435 //_______________________________________________________________________________\r
 436 \r
 437 void AliAnalysisTaskFlavourJetCorrelations::Terminate(Option_t*)\r
 438 {    \r
 439    // The Terminate() function is the last function to be called during\r
 440    // a query. It always runs on the client, it can be used to present\r
 441    // the results graphically or save the results to file.\r
 442    \r
 443    Info("Terminate"," terminate");\r
 444    AliAnalysisTaskSE::Terminate();\r
 445    \r
 446    fmyOutput = dynamic_cast<TList*> (GetOutputData(1));\r
 447    if (!fmyOutput) {     \r
 448       printf("ERROR: fmyOutput not available\n");\r
 449       return;\r
 450    }\r
 451 }\r
 452 \r
 453 //_______________________________________________________________________________\r
 454 \r
 455 void  AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t range, Int_t pdg){\r
 456    Float_t mass=0;\r
 457    Int_t abspdg=TMath::Abs(pdg);\r
 458    \r
 459    mass=TDatabasePDG::Instance()->GetParticle(abspdg)->Mass();\r
 460    // compute the Delta mass for the D*\r
 461    if(fCandidateType==kDstartoKpipi){\r
 462       Float_t mass1=0;\r
 463       mass1=TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
 464       mass = mass-mass1;\r
 465    }\r
 466    \r
 467    fMinMass = mass-range;\r
 468    fMaxMass = mass+range;\r
 469    \r
 470    AliInfo(Form("Setting mass limits to %f, %f",fMinMass,fMaxMass));\r
 471    if (fMinMass<0 || fMaxMass<=0 || fMaxMass<fMinMass) AliFatal("Wrong mass limits!\n");\r
 472 }\r
 473 \r
 474 //_______________________________________________________________________________\r
 475 \r
 476 void  AliAnalysisTaskFlavourJetCorrelations::SetMassLimits(Double_t lowlimit, Double_t uplimit){\r
 477    if(uplimit>lowlimit)\r
 478    {\r
 479       fMinMass = lowlimit;\r
 480       fMaxMass = uplimit;\r
 481    }\r
 482    else{\r
 483       printf("Error! Lower limit larger than upper limit!\n");\r
 484       fMinMass = uplimit - uplimit*0.2;\r
 485       fMaxMass = uplimit;\r
 486    }\r
 487 }\r
 488 \r
 489 //_______________________________________________________________________________\r
 490 \r
 491 Bool_t AliAnalysisTaskFlavourJetCorrelations::SetD0WidthForDStar(Int_t nptbins,Float_t *width){\r
 492    if(nptbins>30) {\r
 493       AliInfo("Maximum number of bins allowed is 30!");\r
 494       return kFALSE;\r
 495    }\r
 496    if(!width) return kFALSE;\r
 497    for(Int_t ipt=0;ipt<nptbins;ipt++) fSigmaD0[ipt]=width[ipt];\r
 498    return kTRUE;\r
 499 }\r
 500 \r
 501 //_______________________________________________________________________________\r
 502 \r
 503 Double_t AliAnalysisTaskFlavourJetCorrelations::Z(AliVParticle* part,AliEmcalJet* jet) const{\r
 504    if(!part || !jet){\r
 505       printf("Particle or jet do not exist!\n");\r
 506       return -99;\r
 507    }\r
 508    Double_t p[3],pj[3];\r
 509    Bool_t okpp=part->PxPyPz(p);\r
 510    Bool_t okpj=jet->PxPyPz(pj);\r
 511    if(!okpp || !okpj){\r
 512       printf("Problems getting momenta\n");\r
 513       return -999;\r
 514    }\r
 515    Double_t pjet=jet->P();\r
 516    Double_t z=(p[0]*pj[0]+p[1]*pj[1]+p[2]*pj[2])/(pjet*pjet);\r
 517    return z;\r
 518 }\r
 519 \r
 520 //_______________________________________________________________________________\r
 521 \r
 522 Bool_t  AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){\r
 523    \r
 524    // Statistics \r
 525    TH1I* hstat=new TH1I("hstat","Statistics",6,-0.5,5.5);\r
 526    hstat->GetXaxis()->SetBinLabel(1,"N ev anal");\r
 527    hstat->GetXaxis()->SetBinLabel(2,"N ev sel");\r
 528    hstat->GetXaxis()->SetBinLabel(3,"N cand sel cuts");\r
 529    hstat->GetXaxis()->SetBinLabel(4,"N jets");\r
 530    hstat->GetXaxis()->SetBinLabel(5,"N cand in jet");\r
 531    hstat->GetXaxis()->SetBinLabel(6,"N jet rej");\r
 532    hstat->SetNdivisions(1);\r
 533    fmyOutput->Add(hstat);\r
 534    \r
 535    const Int_t nbinsmass=200;\r
 536    const Int_t nbinsptjet=500;\r
 537    const Int_t nbinsptD=100;\r
 538    const Int_t nbinsz=100;\r
 539    const Int_t nbinsphi=200;\r
 540    \r
 541    const Float_t ptjetlims[2]={0.,200.};\r
 542    const Float_t ptDlims[2]={0.,50.};\r
 543    const Float_t zlims[2]={0.,1.2};\r
 544    const Float_t philims[2]={0.,6.3};\r
 545    \r
 546    if(fCandidateType==kDstartoKpipi) \r
 547    {\r
 548       \r
 549       TH2F* hDiffSideBand = new TH2F("hDiffSideBand","M(kpipi)-M(kpi) Side Band Background",nbinsmass,fMinMass,fMaxMass,nbinsptD, ptDlims[0],ptDlims[1]);\r
 550       hDiffSideBand->SetStats(kTRUE);\r
 551       hDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");\r
 552       hDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");\r
 553       hDiffSideBand->Sumw2();\r
 554       fmyOutput->Add(hDiffSideBand); \r
 555       \r
 556       \r
 557       TH1F* hPtPion = new TH1F("hPtPion","Primary pions candidates pt ",500,0,10);\r
 558       hPtPion->SetStats(kTRUE);\r
 559       hPtPion->GetXaxis()->SetTitle("GeV/c");\r
 560       hPtPion->GetYaxis()->SetTitle("Entries");\r
 561       hPtPion->Sumw2();\r
 562       fmyOutput->Add(hPtPion);\r
 563       \r
 564    }\r
 565    \r
 566    // jet related fistograms\r
 567    \r
 568    TH1F* hEjetTrks      = new TH1F("hEjetTrks",  "Jet tracks energy distribution;Energy (GeV)",500,0,200);\r
 569    hEjetTrks->Sumw2();\r
 570    TH1F* hPhiJetTrks    = new TH1F("hPhiJetTrks","Jet tracks #phi distribution; #phi (rad)",  nbinsphi,philims[0],philims[1]);\r
 571    hPhiJetTrks->Sumw2();\r
 572    TH1F* hEtaJetTrks    = new TH1F("hEtaJetTrks","Jet tracks #eta distribution; #eta",  100,-1.5,1.5);\r
 573    hEtaJetTrks->Sumw2();\r
 574    TH1F* hPtJetTrks     = new TH1F("hPtJetTrks",  "Jet tracks Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);\r
 575    hPtJetTrks->Sumw2();\r
 576    \r
 577    TH1F* hEjet      = new TH1F("hEjet",  "Jet energy distribution;Energy (GeV)",500,0,200);\r
 578    hEjet->Sumw2();\r
 579    TH1F* hPhiJet    = new TH1F("hPhiJet","Jet #phi distribution; #phi (rad)",  nbinsphi,philims[0],philims[1]);\r
 580    hPhiJet->Sumw2();\r
 581    TH1F* hEtaJet    = new TH1F("hEtaJet","Jet #eta distribution; #eta",  100,-1.5,1.5);\r
 582    hEtaJet->Sumw2();\r
 583    TH1F* hPtJet      = new TH1F("hPtJet",  "Jet Pt distribution; p_{T} (GeV/c)",nbinsptjet,ptjetlims[0],ptjetlims[1]);\r
 584    hPtJet->Sumw2();\r
 585    \r
 586    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]);\r
 587    hPtJetWithD->Sumw2();\r
 588    //for the MC this histogram is filled with the real background\r
 589    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]);\r
 590    hPtJetWithDsb->Sumw2();\r
 591    \r
 592    TH1F* hdeltaRJetTracks=new TH1F("hdeltaRJetTracks","Delta R of tracks in the jets",200, 0.,10.);\r
 593    hdeltaRJetTracks->Sumw2();\r
 594    \r
 595    TH1F* hNtrArr= new TH1F("hNtrArr", "Number of tracks in the array of jets; number of tracks",500,0,1000);\r
 596    hNtrArr->Sumw2();\r
 597    TH1F *hNJetPerEv=new TH1F("hNJetPerEv","Number of jets used per event; number of jets/ev",10,-0.5,9.5);\r
 598    hNJetPerEv->Sumw2();\r
 599    \r
 600    fmyOutput->Add(hEjetTrks);\r
 601    fmyOutput->Add(hPhiJetTrks);\r
 602    fmyOutput->Add(hEtaJetTrks);\r
 603    fmyOutput->Add(hPtJetTrks);\r
 604    fmyOutput->Add(hEjet);\r
 605    fmyOutput->Add(hPhiJet);\r
 606    fmyOutput->Add(hEtaJet);\r
 607    fmyOutput->Add(hPtJet);\r
 608    fmyOutput->Add(hPtJetWithD);\r
 609    fmyOutput->Add(hPtJetWithDsb);\r
 610    fmyOutput->Add(hdeltaRJetTracks);\r
 611    fmyOutput->Add(hNtrArr);\r
 612    fmyOutput->Add(hNJetPerEv);\r
 613    \r
 614    TH1F* hDeltaRD=new TH1F("hDeltaRD","#Delta R distribution of D candidates selected;#Delta R",200, 0.,10.);\r
 615    hDeltaRD->Sumw2();\r
 616    fmyOutput->Add(hDeltaRD);\r
 617    \r
 618    //background (side bands for the Dstar and like sign for D0)\r
 619    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]);\r
 620    hInvMassptD->SetStats(kTRUE);\r
 621    hInvMassptD->GetXaxis()->SetTitle("mass (GeV)");\r
 622    hInvMassptD->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");\r
 623    hInvMassptD->Sumw2();\r
 624    \r
 625    fmyOutput->Add(hInvMassptD);\r
 626    \r
 627    if(fUseMCInfo){\r
 628       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]);\r
 629       hInvMassptDbg->GetXaxis()->SetTitle(Form("%s (GeV)",(fCandidateType==kDstartoKpipi) ? "M(Kpipi)" : "M(Kpi)"));\r
 630       hInvMassptDbg->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");\r
 631       hInvMassptDbg->Sumw2();\r
 632       fmyOutput->Add(hInvMassptDbg);\r
 633       \r
 634    }\r
 635    const Int_t nAxis=5;\r
 636    const Int_t nbinsSparse[nAxis]={nbinsz,nbinsphi,nbinsptjet,nbinsptD,nbinsmass};\r
 637    const Double_t minSparse[nAxis]={zlims[0],philims[0],ptjetlims[0],ptDlims[0],fMinMass};\r
 638    const Double_t maxSparse[nAxis]={zlims[1],philims[1],ptjetlims[1],ptDlims[1],fMinMass};\r
 639    THnSparseF *hsDphiz=new THnSparseF("hsDphiz","Z and #Delta#phi vs p_{T}^{jet}, p_{T}^{D}, and mass", nAxis, nbinsSparse, minSparse, maxSparse);\r
 640    hsDphiz->Sumw2();\r
 641    \r
 642    fmyOutput->Add(hsDphiz);\r
 643 \r
 644    PostData(1, fmyOutput);\r
 645    \r
 646    return kTRUE; \r
 647 }\r
 648 \r
 649 //_______________________________________________________________________________\r
 650 \r
 651 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsRecoJetCorr(AliVParticle* candidate, AliEmcalJet *jet){\r
 652    \r
 653    Double_t ptD=candidate->Pt();\r
 654    Double_t ptjet=jet->Pt();\r
 655    Double_t deltaR=DeltaR(candidate,jet);\r
 656    Double_t deltaphi = jet->Phi()-candidate->Phi();\r
 657    if(deltaphi<=-(TMath::Pi())/2) deltaphi = deltaphi+2*(TMath::Pi());\r
 658    if(deltaphi>(3*(TMath::Pi()))/2) deltaphi = deltaphi-2*(TMath::Pi());\r
 659    Double_t z=Z(candidate,jet);\r
 660    \r
 661    TH1F* hDeltaRD=(TH1F*)fmyOutput->FindObject("hDeltaRD");\r
 662    hDeltaRD->Fill(deltaR); \r
 663    if(fUseReco){\r
 664       if(fCandidateType==kD0toKpi) {\r
 665          AliAODRecoDecayHF* dzero=(AliAODRecoDecayHF*)candidate;\r
 666          FillHistogramsD0JetCorr(dzero,deltaphi,z,ptD,ptjet,deltaR, AODEvent());\r
 667          \r
 668       }\r
 669       \r
 670       if(fCandidateType==kDstartoKpipi) {\r
 671          AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;\r
 672          FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,ptjet,deltaR);\r
 673          \r
 674       }\r
 675    } else{ //generated level\r
 676       //AliInfo("Non reco");\r
 677       FillHistogramsMCGenDJetCorr(deltaphi,z,ptD,ptjet,deltaR);\r
 678    }\r
 679 }\r
 680 \r
 681 //_______________________________________________________________________________\r
 682 \r
 683 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj,Double_t deltaR, AliAODEvent* aodEvent){\r
 684 \r
 685   //dPhi and z not used at the moment,but will be (re)added\r
 686 \r
 687    Double_t masses[2]={0.,0.};\r
 688    Int_t pdgdaughtersD0[2]={211,321};//pi,K \r
 689    Int_t pdgdaughtersD0bar[2]={321,211};//K,pi \r
 690    \r
 691    masses[0]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0\r
 692    masses[1]=candidate->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar\r
 693    \r
 694    TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");\r
 695    THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");\r
 696    Double_t point[5]={z,dPhi,ptj,ptD,masses[0]};\r
 697    \r
 698    Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);\r
 699    if(isselected==1 || isselected==3) {\r
 700       \r
 701       if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[0],ptD);\r
 702       \r
 703       FillMassHistograms(masses[0], ptD, deltaR);\r
 704       hsDphiz->Fill(point,1.);\r
 705    }\r
 706    if(isselected>=2) {\r
 707       if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[1],ptD);\r
 708       \r
 709       FillMassHistograms(masses[1], ptD, deltaR);\r
 710       point[4]=masses[1];\r
 711       hsDphiz->Fill(point,1.);\r
 712    }\r
 713    \r
 714 }\r
 715 \r
 716 //_______________________________________________________________________________\r
 717 \r
 718 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj, Double_t deltaR){\r
 719   //dPhi and z not used at the moment,but will be (re)added\r
 720 \r
 721    AliAODTrack *softpi = (AliAODTrack*)dstar->GetBachelor();\r
 722    Double_t deltamass= dstar->DeltaInvMass(); \r
 723    //Double_t massD0= dstar->InvMassD0();\r
 724    \r
 725    \r
 726    TH1F* hPtPion=(TH1F*)fmyOutput->FindObject("hPtPion");\r
 727    hPtPion->Fill(softpi->Pt());\r
 728    \r
 729    TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");\r
 730    THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");\r
 731    Double_t point[5]={z,dPhi,ptj,ptD,deltamass};\r
 732 \r
 733    if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,deltamass,ptD);\r
 734    \r
 735    FillMassHistograms(deltamass, ptD, deltaR);\r
 736    hsDphiz->Fill(point,1.);\r
 737    \r
 738 }\r
 739 \r
 740 //_______________________________________________________________________________\r
 741 \r
 742 void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsMCGenDJetCorr(Double_t dPhi,Double_t z,Double_t ptD,Double_t ptjet,Double_t deltaR){\r
 743    \r
 744    Double_t pdgmass=0;\r
 745    TH3F* hPtJetWithD=(TH3F*)fmyOutput->FindObject("hPtJetWithD");\r
 746    THnSparseF* hsDphiz=(THnSparseF*)fmyOutput->FindObject("hsDphiz");\r
 747    Double_t point[5]={z,dPhi,ptjet,ptD,pdgmass};\r
 748 \r
 749    if(fCandidateType==kD0toKpi) pdgmass=TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
 750    if(fCandidateType==kDstartoKpipi) pdgmass=TDatabasePDG::Instance()->GetParticle(413)->Mass()-TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
 751    point[4]=pdgmass;\r
 752 \r
 753    if(deltaR<fJetRadius) {\r
 754      hPtJetWithD->Fill(ptjet,pdgmass,ptD); // candidates within a jet\r
 755      hsDphiz->Fill(point,1.);\r
 756    }\r
 757    \r
 758 }\r
 759 \r
 760 //_______________________________________________________________________________\r
 761 \r
 762 void AliAnalysisTaskFlavourJetCorrelations::FillMassHistograms(Double_t mass,Double_t ptD, Double_t deltaR){\r
 763    \r
 764    if(deltaR<fJetRadius) {\r
 765       TH2F* hInvMassptD=(TH2F*)fmyOutput->FindObject("hInvMassptD");\r
 766       hInvMassptD->Fill(mass,ptD);\r
 767    }\r
 768 }\r
 769 \r
 770 void AliAnalysisTaskFlavourJetCorrelations::FlagFlavour(AliVParticle *charm, AliEmcalJet *jet){\r
 771    Double_t deltaR=DeltaR(charm, jet);\r
 772    AliEmcalJet::EFlavourTag tag=AliEmcalJet::kDStar;\r
 773    if (fCandidateType==kD0toKpi) tag=AliEmcalJet::kD0;\r
 774    if (deltaR<fJetRadius) jet->AddFlavourTag(tag);\r
 775    \r
 776    return;\r
 777    \r
 778 }\r
 779 \r
 780 //_______________________________________________________________________________\r
 781 \r
 782 void AliAnalysisTaskFlavourJetCorrelations::SideBandBackground(AliAODRecoCascadeHF *candDstar, AliEmcalJet *jet){\r
 783    \r
 784    //  D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas \r
 785    // (expected detector resolution) on the left and right frm the D0 mass. Each band\r
 786    //  has a width of ~5 sigmas. Two band needed  for opening angle considerations   \r
 787    TH2F* hDiffSideBand=(TH2F*)fmyOutput->FindObject("hDiffSideBand");\r
 788    TH3F* hPtJetWithDsb=(TH3F*)fmyOutput->FindObject("hPtJetWithDsb");\r
 789    \r
 790    Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
 791    \r
 792    Int_t bin = fCuts->PtBin(candDstar->Pt());\r
 793    Float_t fourSigmal= mPDGD0-4.*fSigmaD0[bin] , sixSigmal= mPDGD0-8.*fSigmaD0[bin];\r
 794    Float_t fourSigmar= mPDGD0+4.*fSigmaD0[bin] , sixSigmar= mPDGD0+8.*fSigmaD0[bin];\r
 795    \r
 796    Double_t invM=candDstar->InvMassD0(),  deltaM=candDstar->DeltaInvMass(); \r
 797    //Printf("Inv mass = %f between %f and %f or %f and %f?",invM, sixSigmal,fourSigmal,fourSigmar,sixSigmar);\r
 798    Double_t ptD=candDstar->Pt();\r
 799    Double_t ptjet=jet->Pt();\r
 800    Double_t dPhi=jet->Phi()-candDstar->Phi();\r
 801    Double_t deltaR=DeltaR(candDstar,jet);\r
 802    if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();\r
 803    if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();\r
 804    \r
 805    //fill the background histogram with the side bands or when looking at MC with the real background\r
 806    if((invM>sixSigmal && invM<fourSigmal) || (invM>fourSigmar && invM<=sixSigmar)){  \r
 807       hDiffSideBand->Fill(deltaM,ptD); // M(Kpipi)-M(Kpi) side band background    \r
 808       //hdeltaPhiDjaB->Fill(deltaM,ptD,dPhi);\r
 809       \r
 810       if(deltaR<fJetRadius){  // evaluate in the near side      \r
 811          //hzptDB->Fill(Z(candDstar,jet),deltaM,ptD);\r
 812          hPtJetWithDsb->Fill(ptjet,deltaM,ptD);\r
 813       }\r
 814    }\r
 815 }\r
 816 \r
 817 //_______________________________________________________________________________\r
 818 \r
 819 void AliAnalysisTaskFlavourJetCorrelations::MCBackground(AliAODRecoDecayHF *candbg, AliEmcalJet *jet){\r
 820    \r
 821    //need mass, deltaR, pt jet, ptD\r
 822    //two cases: D0 and Dstar\r
 823    \r
 824    Int_t isselected=fCuts->IsSelected(candbg,AliRDHFCuts::kAll, AODEvent());\r
 825    if(!isselected) return;\r
 826    \r
 827    Double_t ptD=candbg->Pt();\r
 828    Double_t ptjet=jet->Pt();\r
 829    Double_t deltaR=DeltaR(candbg,jet);\r
 830    \r
 831    TH2F* hInvMassptDbg=(TH2F*)fmyOutput->FindObject("hInvMassptDbg");\r
 832    TH3F* hPtJetWithDsb=(TH3F*)fmyOutput->FindObject("hPtJetWithDsb");\r
 833    \r
 834    \r
 835    if(fCandidateType==kDstartoKpipi){\r
 836       AliAODRecoCascadeHF* dstarbg = (AliAODRecoCascadeHF*)candbg;\r
 837       Double_t deltaM=dstarbg->DeltaInvMass();\r
 838       hInvMassptDbg->Fill(deltaM,ptD);\r
 839       if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,deltaM,ptD);\r
 840    }\r
 841    \r
 842    if(fCandidateType==kD0toKpi){\r
 843       Double_t masses[2]={0.,0.};\r
 844       Int_t pdgdaughtersD0[2]={211,321};//pi,K \r
 845       Int_t pdgdaughtersD0bar[2]={321,211};//K,pi \r
 846       \r
 847       masses[0]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0); //D0\r
 848       masses[1]=candbg->InvMass(fNProngs,(UInt_t*)pdgdaughtersD0bar); //D0bar\r
 849       \r
 850       \r
 851       if(isselected==1 || isselected==3) {\r
 852          if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[0],ptD);\r
 853          hInvMassptDbg->Fill(masses[0],ptD);\r
 854       }\r
 855       if(isselected>=2) {\r
 856          if(deltaR<fJetRadius) hPtJetWithDsb->Fill(ptjet,masses[1],ptD);\r
 857          hInvMassptDbg->Fill(masses[1],ptD);\r
 858       }\r
 859       \r
 860       \r
 861    }\r
 862 }\r
 863 \r
 864 //_______________________________________________________________________________\r
 865 \r
 866 Float_t AliAnalysisTaskFlavourJetCorrelations::DeltaR(AliVParticle *p1, AliVParticle *p2) const {\r
 867    //Calculate DeltaR between p1 and p2: DeltaR=sqrt(Delataphi^2+DeltaEta^2)\r
 868    \r
 869    if(!p1 || !p2) return -1;\r
 870    Double_t phi1=p1->Phi(),eta1=p1->Eta();\r
 871    Double_t phi2 = p2->Phi(),eta2 = p2->Eta() ;\r
 872    \r
 873    Double_t dPhi=phi1-phi2;\r
 874    if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());\r
 875    if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());\r
 876    \r
 877    Double_t dEta=eta1-eta2;\r
 878    Double_t deltaR=TMath::Sqrt(dEta*dEta + dPhi*dPhi );\r
 879    return deltaR;\r
 880    \r
 881 }\r