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