]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/FlavourJetTasks/AliAnalysisTaskFlavourJetCorrelations.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updates from Chiara
[u/mrichter/AliRoot.git] / PWGJE / FlavourJetTasks / AliAnalysisTaskFlavourJetCorrelations.cxx
CommitLineData
3150c932 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
3150c932 49ClassImp(AliAnalysisTaskFlavourJetCorrelations)
50
51//__________________________________________________________________________
52AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations() :
a429c767 53 AliAnalysisTaskEmcalJet("",kFALSE),
3150c932 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//___________________________________________________________________________
73AliAnalysisTaskFlavourJetCorrelations::AliAnalysisTaskFlavourJetCorrelations(const Char_t* name, AliRDHFCuts* cuts,ECandidateType candtype, TString jetArrName) :
a429c767 74 AliAnalysisTaskEmcalJet(name,kFALSE),
3150c932 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());
a429c767 135
3150c932 136 DefineOutput(1,TList::Class()); // histos
137 DefineOutput(2,AliRDHFCuts::Class()); // my cuts
a429c767 138
3150c932 139}
140//___________________________________________________________________________
141AliAnalysisTaskFlavourJetCorrelations::~AliAnalysisTaskFlavourJetCorrelations() {
142 //
143 // destructor
144 //
145
146 Info("~AliAnalysisTaskFlavourJetCorrelations","Calling Destructor");
147
148 delete fOutput;
149 delete fCuts;
150
151}
152
153//___________________________________________________________
154void AliAnalysisTaskFlavourJetCorrelations::Init(){
155 //
156 // Initialization
157 //
a429c767 158
3150c932 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
a429c767 184//___________________________________________________________________________
185
186void 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
3150c932 199//_________________________________________________
200void 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);
a429c767 375 if(!AcceptJet(jet)) {
376 hstat->Fill(5);
377 continue;
378 }
3150c932 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));
3150c932 426 PostData(1,fOutput);
427
428}
429
430//________________________________________ terminate ___________________________
431void 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}
3150c932 446
3150c932 447//_________________________________________________________________
448void 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//_________________________________________________________________
467void 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//__________________________________________________________________
481Bool_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//__________________________________________________________________
492Double_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
510Bool_t AliAnalysisTaskFlavourJetCorrelations::DefineHistoForAnalysis(){
511
512 // Statistics
a429c767 513 TH1I* hstat=new TH1I("hstat","Statistics",6,-0.5,5.5);
3150c932 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");
a429c767 519 hstat->GetXaxis()->SetBinLabel(6,"N jet rej");
3150c932 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);
a429c767 536 // hDiff->GetXaxis()->SetTitle("M(kpipi)-M(kpi) GeV");
537 // hDiff->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
3150c932 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);
a429c767 542 hDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV");
543 hDiffSideBand->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
3150c932 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);
a429c767 589 TH1F* hPtJetTrks = new TH1F("hPtJetTrks", "Jet tracks Pt distribution; p_{T} (GeV/c)",500,0,200);
3150c932 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);
a429c767 594 TH1F* hPtJet = new TH1F("hPtJet", "Jet Pt distribution; p_{T} (GeV/c)",500,0,200);
3150c932 595
a429c767 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);
3150c932 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);
a429c767 612 fOutput->Add(hPtJetWithDsb);
3150c932 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());
a429c767 636 hdeltaPhiDja->GetXaxis()->SetTitle("mass (GeV)");
3150c932 637 hdeltaPhiDja->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
638 hdeltaPhiDja->GetYaxis()->SetTitle("#Delta#phi (rad)");
a429c767 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)");
3150c932 642 // hdeltaPhiDjl->GetYaxis()->SetTitle("#Delta#phi (rad)");
a429c767 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());
3150c932 644 // hdeltaPhiDjh->GetXaxis()->SetTitle("mass (GeV/c)");
a429c767 645 // hdeltaPhiDjh->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
3150c932 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());
a429c767 654 hdeltaPhiDjaB->GetXaxis()->SetTitle("mass (GeV)");
3150c932 655 hdeltaPhiDjaB->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
656 hdeltaPhiDjaB->GetYaxis()->SetTitle("#Delta#phi (rad)");
a429c767 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)");
3150c932 660 // hdeltaPhiDjlB->GetYaxis()->SetTitle("#Delta#phi (rad)");
a429c767 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)");
3150c932 664 // hdeltaPhiDjhB->GetYaxis()->SetTitle("#Delta#phi (rad)");
665 fOutput->Add(hdeltaPhiDjaB);
666 // fOutput->Add(hdeltaPhiDjlB);
667 // fOutput->Add(hdeltaPhiDjhB);
668
a429c767 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.);
3150c932 670 hInvMassptD->SetStats(kTRUE);
a429c767 671 hInvMassptD->GetXaxis()->SetTitle("mass (GeV)");
3150c932 672 hInvMassptD->GetYaxis()->SetTitle("p_{t}^{D} (GeV/c)");
673
674 fOutput->Add(hInvMassptD);
a429c767 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.);
3150c932 676 // fMasspjDeltaR->SetStats(kFALSE);
677 // fOutput->Add(fMasspjDeltaR);
678
a429c767 679 TH3F* hzptD=new TH3F("hzptD",Form("Fragmentation function (DeltaR < %f)",fJetRadius),100,0.,1.2,nbinsmass,fMinMass,fMaxMass,100,0,50);
3150c932 680 hzptD->SetStats(kTRUE);
681 hzptD->GetXaxis()->SetTitle("z=p_{D} #cdot p_{j}/|p_{j}|^{2}");
a429c767 682 hzptD->GetYaxis()->SetTitle("mass (GeV)");
3150c932 683 hzptD->GetZaxis()->SetTitle("p_{t}^{D} (GeV/c)");
684 fOutput->Add(hzptD);
685
a429c767 686 TH3F* hzptDB=new TH3F("hzptDB",Form("Fragmentation function (DeltaR < %f) - Side Bands",fJetRadius),100,0.,1.2,nbinsmass,fMinMass,fMaxMass,100,0.,50.);
3150c932 687 hzptDB->SetStats(kTRUE);
688 hzptDB->GetXaxis()->SetTitle("z=p_{D} #cdot p_{j}/|p_{j}|^{2}");
a429c767 689 hzptDB->GetYaxis()->SetTitle("mass (GeV)");
3150c932 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);
a429c767 697 PostData(1, fOutput);
698
3150c932 699 return kTRUE;
700}
701
702void 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);
3150c932 716
717 if(fCandidateType==kD0toKpi) {
718
a429c767 719 FillHistogramsD0JetCorr(candidate,deltaphi,z,ptD,ptjet,deltaR, AODEvent());
3150c932 720
721 }
722
723 if(fCandidateType==kDstartoKpipi) {
724 AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)candidate;
a429c767 725 FillHistogramsDstarJetCorr(dstar,deltaphi,z,ptD,ptjet,deltaR);
3150c932 726
727 }
728
729}
730
a429c767 731void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsD0JetCorr(AliAODRecoDecayHF* candidate, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj,Double_t deltaR, AliAODEvent* aodEvent){
3150c932 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
a429c767 739 TH2F* hPtJetWithD=(TH2F*)fOutput->FindObject("hPtJetWithD");
740
3150c932 741 Int_t isselected=fCuts->IsSelected(candidate,AliRDHFCuts::kAll,aodEvent);
742 if(isselected==1 || isselected==3) {
a429c767 743
744 if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[0]);
3150c932 745
a429c767 746 FillHistograms(masses[0],dPhi,z, ptD, deltaR);
3150c932 747 }
748 if(isselected>=2) {
a429c767 749 if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,masses[1]);
3150c932 750
a429c767 751 FillHistograms(masses[1],dPhi,z, ptD, deltaR);
3150c932 752 }
753
754}
755
a429c767 756void AliAnalysisTaskFlavourJetCorrelations::FillHistogramsDstarJetCorr(AliAODRecoCascadeHF* dstar, Double_t dPhi, Double_t z, Double_t ptD, Double_t ptj, Double_t deltaR){
3150c932 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());
a429c767 763
764 TH2F* hPtJetWithD=(TH2F*)fOutput->FindObject("hPtJetWithD");
765 if(deltaR<fJetRadius) hPtJetWithD->Fill(ptj,deltamass);
766
767 FillHistograms(deltamass,dPhi,z, ptD, deltaR);
3150c932 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
a429c767 775 TH1F* hPtJetWithDsb=(TH1F*)fOutput->FindObject("hPtJetWithDsb");
776
3150c932 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
a429c767 783 if((massD0>sixSigmal && massD0<=fourSigmal) || (massD0>fourSigmar && massD0<=sixSigmar)){
3150c932 784 hDiffSideBand->Fill(deltamass,ptD); // M(Kpipi)-M(Kpi) side band background
785 hdeltaPhiDjaB->Fill(deltamass,ptD,dPhi);
786
a429c767 787 if(deltaR<fJetRadius){ // evaluate in the near side
788 hzptDB->Fill(z,deltamass,ptD);
789 hPtJetWithDsb->Fill(ptj);
3150c932 790 }
a429c767 791
3150c932 792 } // SideBandBackground(dstar, dPhi, z, ptD, deltaR);
793
794}
795
a429c767 796void AliAnalysisTaskFlavourJetCorrelations::FillHistograms(Double_t mass,Double_t dphi, Double_t z,Double_t ptD, Double_t deltaR){
3150c932 797 TH3F* hdeltaPhiDja=((TH3F*)fOutput->FindObject("hdeltaPhiDja"));
798 hdeltaPhiDja->Fill(mass,ptD,dphi);
799
a429c767 800 if(deltaR<fJetRadius) {
3150c932 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
a429c767 840// if(deltaR<fJetRadius){ // evaluate in the near side
3150c932 841// hzptDB->Fill(Z(candDstar,jet),deltaM,ptD);
842// }
843
844// }
845// }
846
847Float_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
867Bool_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
875Bool_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
883Bool_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
889Int_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
908Int_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
932Bool_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
a429c767 985 if(deltaRD<fJetRadius) {
3150c932 986 hstat->Fill(7);
a429c767 987 //comment if you prefer to ask for DeltaR of the daughters < fJetRadius and uncomment below
3150c932 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*/