PWGJE/FlavourJetTasks/AliAnalysisTaskFlavourJetCorrelations.cxx

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