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