1 /**************************************************************************
2 * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 // Base class for DStar in Jets Analysis
19 //-----------------------------------------------------------------------
21 // ERC-QGP Utrecht University - a.grelli@uu.nl
22 //-----------------------------------------------------------------------
24 #include <TDatabasePDG.h>
25 #include <TParticle.h>
29 #include "AliAnalysisTaskSEDStarJets.h"
30 #include "AliRDHFCutsDStartoKpipi.h"
32 #include "AliMCEvent.h"
33 #include "AliAODMCHeader.h"
34 #include "AliAODHandler.h"
35 #include "AliAnalysisManager.h"
37 #include "AliAODVertex.h"
38 #include "AliAODJet.h"
39 #include "AliAODRecoDecay.h"
40 #include "AliAODRecoCascadeHF.h"
41 #include "AliAODRecoDecayHF2Prong.h"
42 #include "AliESDtrack.h"
43 #include "AliAODMCParticle.h"
45 ClassImp(AliAnalysisTaskSEDStarJets)
47 //__________________________________________________________________________
48 AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets() :
53 fRequireNormalization(kTRUE),
81 //___________________________________________________________________________
82 AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets(const Char_t* name, AliRDHFCutsDStartoKpipi* cuts) :
83 AliAnalysisTaskSE(name),
87 fRequireNormalization(kTRUE),
112 // Constructor. Initialization of Inputs and Outputs
115 Info("AliAnalysisTaskSEDStarJets","Calling Constructor");
117 DefineOutput(1,TList::Class()); // histos
118 DefineOutput(2,AliRDHFCutsDStartoKpipi::Class()); // my cuts
120 //___________________________________________________________________________
121 AliAnalysisTaskSEDStarJets::~AliAnalysisTaskSEDStarJets() {
126 Info("~AliAnalysisTaskSEDStarJets","Calling Destructor");
137 if (ftrigger) { delete ftrigger; ftrigger = 0;}
138 if (fPtPion) { delete fPtPion; fPtPion = 0;}
139 if (fInvMass) { delete fInvMass; fInvMass = 0;}
140 if (fRECOPtDStar) { delete fRECOPtDStar; fRECOPtDStar = 0;}
141 if (fRECOPtBkg) { delete fRECOPtBkg; fRECOPtBkg = 0;}
142 if (fDStar) { delete fDStar; fDStar = 0;}
143 if (fDiff) { delete fDiff; fDiff = 0;}
144 if (fDiffSideBand) { delete fDiffSideBand; fDiffSideBand = 0;}
145 if (fDStarMass) { delete fDStarMass; fDStarMass = 0;}
146 if (fPhi) { delete fPhi; fPhi = 0;}
147 if (fPhiBkg) { delete fPhiBkg; fPhiBkg = 0;}
148 if (fTrueDiff){ delete fTrueDiff; fTrueDiff = 0;}
149 if (fResZ) { delete fResZ; fResZ = 0;}
150 if (fResZBkg) { delete fResZBkg; fResZBkg = 0;}
151 if (fEjet) { delete fEjet; fEjet = 0;}
152 if (fPhijet) { delete fPhijet; fPhijet = 0;}
153 if (fEtaJet) { delete fEtaJet; fEtaJet = 0;}
154 if (theMCFF) { delete theMCFF; theMCFF = 0;}
155 if (fDphiD0Dstar) { delete fDphiD0Dstar; fDphiD0Dstar = 0;}
156 if (fPtJet) { delete fPtJet; fPtJet = 0;}
159 //___________________________________________________________
160 void AliAnalysisTaskSEDStarJets::Init(){
164 if(fDebug > 1) printf("AnalysisTaskSEDStarJets::Init() \n");
165 AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
167 PostData(2,copyfCuts);
172 //_________________________________________________
173 void AliAnalysisTaskSEDStarJets::UserExec(Option_t *)
177 Error("UserExec","NO EVENT FOUND!");
182 AliInfo(Form("Event %d",fEvents));
183 if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
186 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
188 TClonesArray *arrayDStartoD0pi=0;
190 if(!aodEvent && AODEvent() && IsStandardAOD()) {
191 // In case there is an AOD handler writing a standard AOD, use the AOD
192 // event in memory rather than the input (ESD) event.
193 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
194 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
195 // have to taken from the AOD event hold by the AliAODExtension
196 AliAODHandler* aodHandler = (AliAODHandler*)
197 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
198 if(aodHandler->GetExtensions()) {
199 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
200 AliAODEvent *aodFromExt = ext->GetAOD();
201 arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
204 arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
207 if (!arrayDStartoD0pi){
208 AliInfo("Could not find array of HF vertices, skipping the event");
210 }else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
212 // fix for temporary bug in ESDfilter
213 // the AODs with null vertex pointer didn't pass the PhysSel
214 if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
216 // Simulate a jet triggered sample
217 TClonesArray *arrayofJets = (TClonesArray*)aodEvent->GetJets();
218 if(aodEvent->GetNJets()<=0) return;
220 // counters for efficiencies
221 Int_t icountReco = 0;
223 // Normalization factor
224 if(fRequireNormalization){
228 //D* and D0 prongs needed to MatchToMC method
229 Int_t pdgDgDStartoD0pi[2]={421,211};
230 Int_t pdgDgD0toKpi[2]={321,211};
238 //loop over jets and consider only the leading jey in the event
239 for (Int_t iJets = 0; iJets<arrayofJets->GetEntriesFast(); iJets++) {
240 AliAODJet* jet = (AliAODJet*)arrayofJets->At(iJets);
252 // fill energy, eta and phi of the jet
254 fPhijet ->Fill(phiJet);
255 fEtaJet ->Fill(etaJet);
259 //loop over D* candidates
260 for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
263 AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
264 AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
266 Double_t finvM =-999;
267 Double_t finvMDStar =-999;
270 Int_t mcLabel = -9999;
272 // find associated MC particle for D* ->D0toKpi
274 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
276 printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
279 mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
281 if(mcLabel>=0) isDStar = 1;
284 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcLabel));
285 //fragmentation function in mc
286 zMC= FillMCFF(mcPart,mcArray,mcLabel);
287 if(zMC>0) theMCFF->Fill(zMC);
292 AliAODTrack *track2 = (AliAODTrack*)dstarD0pi->GetBachelor();
293 Double_t pt = dstarD0pi->Pt();
295 // track quality cuts
296 Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
297 if(!isTkSelected) continue;
299 // region of interest + topological cuts + PID
300 if(!fCuts->IsInFiducialAcceptance(dstarD0pi->Pt(),dstarD0pi->YDstar())) continue;
301 Int_t isSelected=fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kCandidate); //selected
302 if (!isSelected) continue;
305 finvM = dstarD0pi->InvMassD0();
306 fInvMass->Fill(finvM);
308 //DStar invariant mass
309 finvMDStar = dstarD0pi->InvMassDstarKpipi();
312 Double_t dStarMom = dstarD0pi->P();
313 Double_t phiDStar = dstarD0pi->Phi();
314 Double_t phiD0 = theD0particle->Phi();
315 //check suggested by Federico
316 Double_t dPhiD0Dstar = phiD0 - phiDStar;
318 dPhi = phiJet - phiDStar;
321 if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
322 if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
324 Double_t corrFactorCharge = (ejet/100)*20;
325 EGjet = ejet + corrFactorCharge;
327 // fill D* candidates
328 Double_t mPDGD0=TDatabasePDG::Instance()->GetParticle(421)->Mass();
329 if(finvM >= (mPDGD0-0.05) && finvM <=(mPDGD0+0.05)){ // ~3 sigma (sigma=17MeV, conservative)
332 fDphiD0Dstar->Fill(dPhiD0Dstar);
333 fDStarMass->Fill(finvMDStar);
334 fTrueDiff->Fill(finvMDStar-finvM);
336 if(isDStar == 0) fDphiD0Dstar->Fill(dPhiD0Dstar); // angle between D0 and D*
338 fDStar->Fill(finvMDStar);
339 fDiff ->Fill(finvMDStar-finvM);
341 Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
342 Double_t invmassDelta = dstarD0pi->DeltaInvMass();
344 // now the dphi signal and the fragmentation function
345 if (TMath::Abs(invmassDelta-(mPDGDstar-mPDGD0))<0.0019){
346 //fill candidates D* and soft pion reco pt
348 fRECOPtDStar->Fill(pt);
349 fPtPion->Fill(track2->Pt());
353 Double_t jetCone = 0.4;
354 if(dPhi>=-jetCone && dPhi<=jetCone){ // evaluate in the near side inside UA1 radius
355 Double_t zFrag = (TMath::Cos(dPhi)*dStarMom)/EGjet;
356 fResZ->Fill(TMath::Abs(zFrag));
360 // evaluate side band background
361 SideBandBackground(finvM, finvMDStar, dStarMom, EGjet, dPhi);
365 AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
370 //________________________________________ terminate ___________________________
371 void AliAnalysisTaskSEDStarJets::Terminate(Option_t*)
373 // The Terminate() function is the last function to be called during
374 // a query. It always runs on the client, it can be used to present
375 // the results graphically or save the results to file.
377 Info("Terminate"," terminate");
378 AliAnalysisTaskSE::Terminate();
380 fOutput = dynamic_cast<TList*> (GetOutputData(1));
382 printf("ERROR: fOutput not available\n");
386 fDStarMass = dynamic_cast<TH1F*>(fOutput->FindObject("fDStarMass"));
387 fTrueDiff = dynamic_cast<TH1F*>(fOutput->FindObject("fTrueDiff"));
388 fInvMass = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass"));
389 fPtPion = dynamic_cast<TH1F*>(fOutput->FindObject("fPtPion "));
390 fDStar = dynamic_cast<TH1F*>(fOutput->FindObject("fDStar"));
391 fDiff = dynamic_cast<TH1F*>(fOutput->FindObject("fDiff"));
392 fDiffSideBand = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand"));
393 ftrigger = dynamic_cast<TH1F*>(fOutput->FindObject("ftrigger"));
394 fRECOPtDStar = dynamic_cast<TH1F*>(fOutput->FindObject("fRECOPtDStar"));
395 fRECOPtBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fRECOPtBkg"));
396 fEjet = dynamic_cast<TH1F*>(fOutput->FindObject("fEjet"));
397 fPhijet = dynamic_cast<TH1F*>(fOutput->FindObject("fPhijet"));
398 fEtaJet = dynamic_cast<TH1F*>(fOutput->FindObject("fEtaJet"));
399 fPhi = dynamic_cast<TH1F*>(fOutput->FindObject("fPhi"));
400 fResZ = dynamic_cast<TH1F*>(fOutput->FindObject("fResZ"));
401 fResZBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fResZBkg"));
402 fPhiBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fPhiBkg"));
403 theMCFF = dynamic_cast<TH1F*>(fOutput->FindObject("theMCFF"));
404 fDphiD0Dstar = dynamic_cast<TH1F*>(fOutput->FindObject("fDphiD0Dstar"));
405 fPtJet = dynamic_cast<TH1F*>(fOutput->FindObject("fPtJet"));
408 //___________________________________________________________________________
410 void AliAnalysisTaskSEDStarJets::UserCreateOutputObjects() {
412 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
416 fOutput = new TList();
419 DefineHistoFroAnalysis();
423 //___________________________________ hiostograms _______________________________________
425 Bool_t AliAnalysisTaskSEDStarJets::DefineHistoFroAnalysis(){
427 // Invariant mass related histograms
428 fInvMass = new TH1F("invMass","Kpi invariant mass distribution",1500,.5,3.5);
429 fInvMass->SetStats(kTRUE);
430 fInvMass->GetXaxis()->SetTitle("GeV/c");
431 fInvMass->GetYaxis()->SetTitle("Entries");
432 fOutput->Add(fInvMass);
434 fDStar = new TH1F("invMassDStar","DStar invariant mass after D0 cuts ",600,1.8,2.4);
435 fDStar->SetStats(kTRUE);
436 fDStar->GetXaxis()->SetTitle("GeV/c");
437 fDStar->GetYaxis()->SetTitle("Entries");
438 fOutput->Add(fDStar);
440 fDiff = new TH1F("Diff","M(kpipi)-M(kpi)",750,0.1,0.2);
441 fDiff->SetStats(kTRUE);
442 fDiff->GetXaxis()->SetTitle("M(kpipi)-M(kpi) GeV/c^2");
443 fDiff->GetYaxis()->SetTitle("Entries");
446 fDiffSideBand = new TH1F("DiffSide","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
447 fDiffSideBand->SetStats(kTRUE);
448 fDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
449 fDiffSideBand->GetYaxis()->SetTitle("Entries");
450 fOutput->Add(fDiffSideBand);
452 fDStarMass = new TH1F("RECODStar2","RECO DStar invariant mass distribution",750,1.5,2.5);
453 fOutput->Add(fDStarMass);
455 fTrueDiff = new TH1F("dstar","True Reco diff",750,0,0.2);
456 fOutput->Add(fTrueDiff);
458 // trigger normalization
459 ftrigger = new TH1F("Normalization","Normalization factor for correlations",1,0,10);
460 ftrigger->SetStats(kTRUE);
461 fOutput->Add(ftrigger);
463 //correlation fistograms
464 fPhi = new TH1F("phi","Delta phi between Jet axis and DStar ",25,-1.57,4.72);
465 fPhi->SetStats(kTRUE);
466 fPhi->GetXaxis()->SetTitle("#Delta #phi (rad)");
467 fPhi->GetYaxis()->SetTitle("Entries");
470 fDphiD0Dstar = new TH1F("phiD0Dstar","Delta phi between D0 and DStar ",1000,-6.5,6.5);
471 fOutput->Add(fDphiD0Dstar);
473 fPhiBkg = new TH1F("phiBkg","Delta phi between Jet axis and DStar background ",25,-1.57,4.72);
474 fPhiBkg->SetStats(kTRUE);
475 fPhiBkg->GetXaxis()->SetTitle("#Delta #phi (rad)");
476 fPhiBkg->GetYaxis()->SetTitle("Entries");
477 fOutput->Add(fPhiBkg);
479 fRECOPtDStar = new TH1F("RECODStar1","RECO DStar pt distribution",30,0,30);
480 fRECOPtDStar->SetStats(kTRUE);
481 fRECOPtDStar->SetLineColor(2);
482 fRECOPtDStar->GetXaxis()->SetTitle("GeV/c");
483 fRECOPtDStar->GetYaxis()->SetTitle("Entries");
484 fOutput->Add(fRECOPtDStar);
486 fRECOPtBkg = new TH1F("RECOptBkg","RECO pt distribution side bands",30,0,30);
487 fRECOPtBkg->SetStats(kTRUE);
488 fRECOPtBkg->SetLineColor(2);
489 fRECOPtBkg->GetXaxis()->SetTitle("GeV/c");
490 fRECOPtBkg->GetYaxis()->SetTitle("Entries");
491 fOutput->Add(fRECOPtBkg);
493 fPtPion = new TH1F("pionpt","Primary pions candidates pt ",500,0,10);
494 fPtPion->SetStats(kTRUE);
495 fPtPion->GetXaxis()->SetTitle("GeV/c");
496 fPtPion->GetYaxis()->SetTitle("Entries");
497 fOutput->Add(fPtPion);
499 // jet related fistograms
500 fEjet = new TH1F("ejet", "UA1 algorithm jet energy distribution",1000,0,500);
501 fPhijet = new TH1F("Phijet","UA1 algorithm jet #phi distribution", 200,-7,7);
502 fEtaJet = new TH1F("Etajet","UA1 algorithm jet #eta distribution", 200,-7,7);
503 fPtJet = new TH1F("PtJet", "UA1 algorithm jet Pt distribution",1000,0,500);
505 fOutput->Add(fPhijet);
506 fOutput->Add(fEtaJet);
507 fOutput->Add(fPtJet);
509 theMCFF = new TH1F("FragFuncMC","Fragmentation function in MC for FC ",100,0,10);
510 fResZ = new TH1F("FragFunc","Fragmentation function ",50,0,1);
511 fResZBkg = new TH1F("FragFuncBkg","Fragmentation function background",50,0,1);
512 fOutput->Add(theMCFF);
514 fOutput->Add(fResZBkg);
519 //______________________________ side band background for D*___________________________________
521 void AliAnalysisTaskSEDStarJets::SideBandBackground(Double_t invM, Double_t invMDStar, Double_t dStarMomBkg, Double_t EGjet, Double_t dPhi){
523 // D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
524 // (expected detector resolution) on the left and right frm the D0 mass. Each band
525 // has a width of ~5 sigmas. Two band needed for opening angle considerations
527 if((invM>=1.7 && invM<=1.8) || (invM>=1.92 && invM<=2.19)){
528 fDiffSideBand->Fill(invMDStar-invM); // M(Kpipi)-M(Kpi) side band background
529 if ((invMDStar-invM)<=0.14732 && (invMDStar-invM)>=0.14352) {
531 fRECOPtBkg->Fill(dStarMomBkg);
532 if(dPhi>=-0.4 && dPhi<=0.4){ // evaluate in the near side
533 Double_t zFragBkg = (TMath::Cos(dPhi)*dStarMomBkg)/EGjet;
534 fResZBkg->Fill(TMath::Abs(zFragBkg));
540 //_____________________________________________________________________________________________-
541 double AliAnalysisTaskSEDStarJets::FillMCFF(AliAODMCParticle* mcPart, TClonesArray* mcArray, Int_t mcLabel){
544 // UA1 jet algorithm reproduced in MC
550 Double_t PhiLeading = -999;
551 Double_t EtaLeading = -999;
552 Double_t PtLeading = -999;
555 if (!mcPart) return zMC2;
556 if (!mcArray) return zMC2;
558 //find leading particle
559 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
560 AliAODMCParticle* Part = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
562 AliWarning("MC Particle not found in tree, skipping");
566 // remove quarks and the leading particle (it will be counted later)
567 if(iPart == mcLabel) continue;
568 if(iPart <= 8) continue;
570 //remove resonances not directly detected in detector
571 Int_t PDGCode = Part->GetPdgCode();
573 // be sure the particle reach the detector
574 Double_t x = Part->Xv();
575 Double_t y = Part->Yv();
576 Double_t z = Part->Zv();
578 if(TMath::Abs(PDGCode)== 2212 && x<3 && y<3) continue;
579 if(TMath::Abs(x)>30 || TMath::Abs(y)>30 || TMath::Abs(z)>30 ) continue;
580 if(TMath::Abs(PDGCode)!=211 && TMath::Abs(PDGCode)!=321 && TMath::Abs(PDGCode)!=11 && TMath::Abs(PDGCode)!=13 && TMath::Abs(PDGCode)!=2212) continue;
587 daug0 = Part->GetDaughter(0);
590 AliAODMCParticle* tdaug = dynamic_cast<AliAODMCParticle*>(mcArray->At(daug0));
597 if(TMath::Abs(xd)<3 || TMath::Abs(yd)<3) continue;
599 Bool_t AliceAcc = (TMath::Abs(Part->Eta()) <= 0.9);
600 if(!AliceAcc) continue;
606 PhiLeading = Part->Phi();
607 EtaLeading = Part->Eta();
608 PtLeading = Part->Pt();
614 Double_t jetEnergy = 0;
616 //reconstruct the jet
617 for (Int_t iiPart=0; iiPart<mcArray->GetEntriesFast(); iiPart++) {
618 AliAODMCParticle* tPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iiPart));
620 AliWarning("MC Particle not found in tree, skipping");
623 // remove quarks and the leading particle (it will be counted later)
624 if(iiPart == counter) continue; // do not count again the leading particle
625 if(iiPart == mcLabel) continue;
626 if(iiPart <= 8) continue;
628 //remove resonances not directly detected in detector
629 Int_t PDGCode = tPart->GetPdgCode();
631 // be sure the particle reach the detector
632 Double_t x = tPart->Xv();
633 Double_t y = tPart->Yv();
634 Double_t z = tPart->Zv();
636 if(TMath::Abs(PDGCode)== 2212 && (x<3 && y<3)) continue;
637 if(TMath::Abs(x)>30 || TMath::Abs(y)>30 || TMath::Abs(z)>30 ) continue; // has to be generated at least in the silicon tracker or beam pipe
638 if(TMath::Abs(PDGCode)!=211 && TMath::Abs(PDGCode)!=321 && TMath::Abs(PDGCode)!=11 && TMath::Abs(PDGCode)!=13 && TMath::Abs(PDGCode)!=2212) continue;
646 daug0 = tPart->GetDaughter(0);
649 AliAODMCParticle* tdaug = dynamic_cast<AliAODMCParticle*>(mcArray->At(daug0));
656 if(TMath::Abs(xd)<3 && TMath::Abs(yd)<3) continue;
657 //remove particles not in ALICE acceptance
658 if(tPart->Pt()<0.07) continue;
659 Bool_t AliceAcc = (TMath::Abs(tPart->Eta()) <= 0.9);
660 if(!AliceAcc) continue;
662 Double_t EtaMCp = tPart->Eta();
663 Double_t PhiMCp = tPart->Phi();
665 Double_t DphiMClead = PhiLeading-PhiMCp;
667 if(DphiMClead<=-(TMath::Pi())/2) DphiMClead = DphiMClead+2*(TMath::Pi());
668 if(DphiMClead>(3*(TMath::Pi()))/2) DphiMClead = DphiMClead-2*(TMath::Pi());
670 Double_t deta = (EtaLeading-EtaMCp);
672 Double_t radius = TMath::Sqrt((DphiMClead*DphiMClead)+(deta*deta));
674 if(radius>0.4) continue; // in the jet cone
675 if(tPart->Charge()==0) continue; // only charged fraction
677 jetEnergy = jetEnergy+(tPart->E());
680 jetEnergy = jetEnergy + leading;
682 // delta phi D*, jet axis
683 Double_t dPhi = PhiLeading - (mcPart->Phi());
686 if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
687 if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
689 zMC2 = (TMath::Cos(dPhi)*(mcPart->P()))/jetEnergy;