]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGHF/vertexingHF/AliAnalysisTaskSEDStarJets.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Fix in the computetion of <pt>
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAnalysisTaskSEDStarJets.cxx
... / ...
CommitLineData
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// Base class for DStar in Jets Analysis
18//
19//-----------------------------------------------------------------------
20// Author A.Grelli
21// ERC-QGP Utrecht University - a.grelli@uu.nl
22//-----------------------------------------------------------------------
23
24#include <TDatabasePDG.h>
25#include <TParticle.h>
26#include <TVector3.h>
27#include "TROOT.h"
28
29#include "AliAnalysisTaskSEDStarJets.h"
30#include "AliRDHFCutsDStartoKpipi.h"
31#include "AliStack.h"
32#include "AliMCEvent.h"
33#include "AliAODMCHeader.h"
34#include "AliAODHandler.h"
35#include "AliAnalysisManager.h"
36#include "AliLog.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"
44
45ClassImp(AliAnalysisTaskSEDStarJets)
46
47//__________________________________________________________________________
48AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets() :
49 AliAnalysisTaskSE(),
50 fEvents(0),
51 fchargeFrCorr(0),
52 fUseMCInfo(kTRUE),
53 fRequireNormalization(kTRUE),
54 fOutput(0),
55 fCuts(0),
56 ftrigger(0),
57 fPtPion(0),
58 fInvMass(0),
59 fRECOPtDStar(0),
60 fRECOPtBkg(0),
61 fDStar(0),
62 fDiff(0),
63 fDiffSideBand(0),
64 fDStarMass(0),
65 fPhi(0),
66 fPhiBkg(0),
67 fTrueDiff(0),
68 fResZ(0),
69 fResZBkg(0),
70 fEjet(0),
71 fPhijet(0),
72 fEtaJet(0),
73 theMCFF(0),
74 fDphiD0Dstar(0),
75 fPtJet(0)
76{
77 //
78 // Default ctor
79 //
80}
81//___________________________________________________________________________
82AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets(const Char_t* name, AliRDHFCutsDStartoKpipi* cuts) :
83 AliAnalysisTaskSE(name),
84 fEvents(0),
85 fchargeFrCorr(0),
86 fUseMCInfo(kTRUE),
87 fRequireNormalization(kTRUE),
88 fOutput(0),
89 fCuts(0),
90 ftrigger(0),
91 fPtPion(0),
92 fInvMass(0),
93 fRECOPtDStar(0),
94 fRECOPtBkg(0),
95 fDStar(0),
96 fDiff(0),
97 fDiffSideBand(0),
98 fDStarMass(0),
99 fPhi(0),
100 fPhiBkg(0),
101 fTrueDiff(0),
102 fResZ(0),
103 fResZBkg(0),
104 fEjet(0),
105 fPhijet(0),
106 fEtaJet(0),
107 theMCFF(0),
108 fDphiD0Dstar(0),
109 fPtJet(0)
110{
111 //
112 // Constructor. Initialization of Inputs and Outputs
113 //
114 fCuts=cuts;
115 Info("AliAnalysisTaskSEDStarJets","Calling Constructor");
116
117 DefineOutput(1,TList::Class()); // histos
118 DefineOutput(2,AliRDHFCutsDStartoKpipi::Class()); // my cuts
119}
120//___________________________________________________________________________
121AliAnalysisTaskSEDStarJets::~AliAnalysisTaskSEDStarJets() {
122 //
123 // destructor
124 //
125
126 Info("~AliAnalysisTaskSEDStarJets","Calling Destructor");
127 if (fOutput) {
128 delete fOutput;
129 fOutput = 0;
130 }
131
132 if (fCuts) {
133 delete fCuts;
134 fCuts = 0;
135 }
136
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;}
157}
158
159//___________________________________________________________
160void AliAnalysisTaskSEDStarJets::Init(){
161 //
162 // Initialization
163 //
164 if(fDebug > 1) printf("AnalysisTaskSEDStarJets::Init() \n");
165 AliRDHFCutsDStartoKpipi* copyfCuts=new AliRDHFCutsDStartoKpipi(*fCuts);
166 // Post the cuts
167 PostData(2,copyfCuts);
168
169 return;
170}
171
172//_________________________________________________
173void AliAnalysisTaskSEDStarJets::UserExec(Option_t *)
174{
175 // user exec
176 if (!fInputEvent) {
177 Error("UserExec","NO EVENT FOUND!");
178 return;
179 }
180
181 fEvents++;
182 AliInfo(Form("Event %d",fEvents));
183 if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
184
185 // Load the event
186 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
187
188 TClonesArray *arrayDStartoD0pi=0;
189
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");
202 }
203 } else {
204 arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
205 }
206
207 if (!arrayDStartoD0pi){
208 AliInfo("Could not find array of HF vertices, skipping the event");
209 return;
210 }else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
211
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;
215
216 // Simulate a jet triggered sample
217 TClonesArray *arrayofJets = (TClonesArray*)aodEvent->GetJets();
218 if(aodEvent->GetNJets()<=0) return;
219
220 // counters for efficiencies
221 Int_t icountReco = 0;
222
223 // Normalization factor
224 if(fRequireNormalization){
225 ftrigger->Fill(1);
226 }
227
228 //D* and D0 prongs needed to MatchToMC method
229 Int_t pdgDgDStartoD0pi[2]={421,211};
230 Int_t pdgDgD0toKpi[2]={321,211};
231
232 Double_t max =0;
233 Double_t ejet = 0;
234 Double_t phiJet = 0;
235 Double_t etaJet = 0;
236 Double_t ptjet = 0;
237
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);
241
242 //jets variables
243 ejet = jet->E();
244
245 if(ejet>max){
246 max = jet->E();
247 phiJet = jet->Phi();
248 etaJet = jet->Eta();
249 ptjet = jet->Pt();
250 }
251
252 // fill energy, eta and phi of the jet
253 fEjet ->Fill(ejet);
254 fPhijet ->Fill(phiJet);
255 fEtaJet ->Fill(etaJet);
256 fPtJet->Fill(ptjet);
257 }
258
259 //loop over D* candidates
260 for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
261
262 // D* candidates
263 AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
264 AliAODRecoDecayHF2Prong* theD0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
265
266 Double_t finvM =-999;
267 Double_t finvMDStar =-999;
268 Double_t dPhi =-999;
269 Bool_t isDStar =0;
270 Int_t mcLabel = -9999;
271
272 // find associated MC particle for D* ->D0toKpi
273 if(fUseMCInfo){
274 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
275 if(!mcArray) {
276 printf("AliAnalysisTaskSEDStarSpectra::UserExec: MC particles not found!\n");
277 return;
278 }
279 mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
280
281 if(mcLabel>=0) isDStar = 1;
282 if(mcLabel>0){
283 Double_t zMC =-999;
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);
288 }
289 }
290
291 // soft pion
292 AliAODTrack *track2 = (AliAODTrack*)dstarD0pi->GetBachelor();
293 Double_t pt = dstarD0pi->Pt();
294
295 // track quality cuts
296 Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
297 if(!isTkSelected) continue;
298
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;
303
304 // fill histos
305 finvM = dstarD0pi->InvMassD0();
306 fInvMass->Fill(finvM);
307
308 //DStar invariant mass
309 finvMDStar = dstarD0pi->InvMassDstarKpipi();
310
311 Double_t EGjet = 0;
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;
317
318 dPhi = phiJet - phiDStar;
319
320 //plot right dphi
321 if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
322 if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
323
324 Double_t corrFactorCharge = (ejet/100)*20;
325 EGjet = ejet + corrFactorCharge;
326
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)
330
331 if(isDStar == 1) {
332 fDphiD0Dstar->Fill(dPhiD0Dstar);
333 fDStarMass->Fill(finvMDStar);
334 fTrueDiff->Fill(finvMDStar-finvM);
335 }
336 if(isDStar == 0) fDphiD0Dstar->Fill(dPhiD0Dstar); // angle between D0 and D*
337
338 fDStar->Fill(finvMDStar);
339 fDiff ->Fill(finvMDStar-finvM);
340
341 Double_t mPDGDstar=TDatabasePDG::Instance()->GetParticle(413)->Mass();
342 Double_t invmassDelta = dstarD0pi->DeltaInvMass();
343
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
347
348 fRECOPtDStar->Fill(pt);
349 fPtPion->Fill(track2->Pt());
350
351 fPhi ->Fill(dPhi);
352
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));
357 }
358 }
359 }
360 // evaluate side band background
361 SideBandBackground(finvM, finvMDStar, dStarMom, EGjet, dPhi);
362
363 } // D* cand
364
365AliDebug(2, Form("Found %i Reco particles that are D*!!",icountReco));
366
367PostData(1,fOutput);
368}
369
370//________________________________________ terminate ___________________________
371void AliAnalysisTaskSEDStarJets::Terminate(Option_t*)
372{
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.
376
377 Info("Terminate"," terminate");
378 AliAnalysisTaskSE::Terminate();
379
380 fOutput = dynamic_cast<TList*> (GetOutputData(1));
381 if (!fOutput) {
382 printf("ERROR: fOutput not available\n");
383 return;
384 }
385
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"));
406
407}
408//___________________________________________________________________________
409
410void AliAnalysisTaskSEDStarJets::UserCreateOutputObjects() {
411 // output
412 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
413
414 //slot #1
415 OpenFile(1);
416 fOutput = new TList();
417 fOutput->SetOwner();
418 // define histograms
419 DefineHistoFroAnalysis();
420
421 return;
422}
423//___________________________________ hiostograms _______________________________________
424
425Bool_t AliAnalysisTaskSEDStarJets::DefineHistoFroAnalysis(){
426
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);
433
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);
439
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");
444 fOutput->Add(fDiff);
445
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);
451
452 fDStarMass = new TH1F("RECODStar2","RECO DStar invariant mass distribution",750,1.5,2.5);
453 fOutput->Add(fDStarMass);
454
455 fTrueDiff = new TH1F("dstar","True Reco diff",750,0,0.2);
456 fOutput->Add(fTrueDiff);
457
458 // trigger normalization
459 ftrigger = new TH1F("Normalization","Normalization factor for correlations",1,0,10);
460 ftrigger->SetStats(kTRUE);
461 fOutput->Add(ftrigger);
462
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");
468 fOutput->Add(fPhi);
469
470 fDphiD0Dstar = new TH1F("phiD0Dstar","Delta phi between D0 and DStar ",1000,-6.5,6.5);
471 fOutput->Add(fDphiD0Dstar);
472
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);
478
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);
485
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);
492
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);
498
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);
504 fOutput->Add(fEjet);
505 fOutput->Add(fPhijet);
506 fOutput->Add(fEtaJet);
507 fOutput->Add(fPtJet);
508
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);
513 fOutput->Add(fResZ);
514 fOutput->Add(fResZBkg);
515
516 return kTRUE;
517}
518
519//______________________________ side band background for D*___________________________________
520
521void AliAnalysisTaskSEDStarJets::SideBandBackground(Double_t invM, Double_t invMDStar, Double_t dStarMomBkg, Double_t EGjet, Double_t dPhi){
522
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
526
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) {
530 fPhiBkg->Fill(dPhi);
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));
535 }
536 }
537 }
538}
539
540//_____________________________________________________________________________________________-
541double AliAnalysisTaskSEDStarJets::FillMCFF(AliAODMCParticle* mcPart, TClonesArray* mcArray, Int_t mcLabel){
542 //
543 // GS from MC
544 // UA1 jet algorithm reproduced in MC
545 //
546 Double_t zMC2 =-999;
547
548 Double_t leading =0;
549 Double_t PartE = 0;
550 Double_t PhiLeading = -999;
551 Double_t EtaLeading = -999;
552 Double_t PtLeading = -999;
553 Int_t counter =-999;
554
555 if (!mcPart) return zMC2;
556 if (!mcArray) return zMC2;
557
558 //find leading particle
559 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
560 AliAODMCParticle* Part = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
561 if (!Part) {
562 AliWarning("MC Particle not found in tree, skipping");
563 continue;
564 }
565
566 // remove quarks and the leading particle (it will be counted later)
567 if(iPart == mcLabel) continue;
568 if(iPart <= 8) continue;
569
570 //remove resonances not directly detected in detector
571 Int_t PDGCode = Part->GetPdgCode();
572
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();
577
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;
581
582 Int_t daug0 = -999;
583 Double_t xd =-999;
584 Double_t yd =-999;
585 Double_t zd =-999;
586
587 daug0 = Part->GetDaughter(0);
588
589 if(daug0>=0){
590 AliAODMCParticle* tdaug = dynamic_cast<AliAODMCParticle*>(mcArray->At(daug0));
591 if(tdaug){
592 xd = tdaug->Xv();
593 yd = tdaug->Yv();
594 zd = tdaug->Zv();
595 }
596 }
597 if(TMath::Abs(xd)<3 || TMath::Abs(yd)<3) continue;
598
599 Bool_t AliceAcc = (TMath::Abs(Part->Eta()) <= 0.9);
600 if(!AliceAcc) continue;
601
602 PartE = Part->E();
603
604 if(PartE>leading){
605 leading = Part->E();
606 PhiLeading = Part->Phi();
607 EtaLeading = Part->Eta();
608 PtLeading = Part->Pt();
609 counter = iPart;
610 }
611
612 }
613
614 Double_t jetEnergy = 0;
615
616 //reconstruct the jet
617 for (Int_t iiPart=0; iiPart<mcArray->GetEntriesFast(); iiPart++) {
618 AliAODMCParticle* tPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iiPart));
619 if (!tPart) {
620 AliWarning("MC Particle not found in tree, skipping");
621 continue;
622 }
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;
627
628 //remove resonances not directly detected in detector
629 Int_t PDGCode = tPart->GetPdgCode();
630
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();
635
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;
639
640
641 Int_t daug0 = -999;
642 Double_t xd =-999;
643 Double_t yd =-999;
644 Double_t zd =-999;
645
646 daug0 = tPart->GetDaughter(0);
647
648 if(daug0>=0){
649 AliAODMCParticle* tdaug = dynamic_cast<AliAODMCParticle*>(mcArray->At(daug0));
650 if(tdaug){
651 xd = tdaug->Xv();
652 yd = tdaug->Yv();
653 zd = tdaug->Zv();
654 }
655 }
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;
661
662 Double_t EtaMCp = tPart->Eta();
663 Double_t PhiMCp = tPart->Phi();
664
665 Double_t DphiMClead = PhiLeading-PhiMCp;
666
667 if(DphiMClead<=-(TMath::Pi())/2) DphiMClead = DphiMClead+2*(TMath::Pi());
668 if(DphiMClead>(3*(TMath::Pi()))/2) DphiMClead = DphiMClead-2*(TMath::Pi());
669
670 Double_t deta = (EtaLeading-EtaMCp);
671 //cone radius
672 Double_t radius = TMath::Sqrt((DphiMClead*DphiMClead)+(deta*deta));
673
674 if(radius>0.4) continue; // in the jet cone
675 if(tPart->Charge()==0) continue; // only charged fraction
676
677 jetEnergy = jetEnergy+(tPart->E());
678 }
679
680 jetEnergy = jetEnergy + leading;
681
682 // delta phi D*, jet axis
683 Double_t dPhi = PhiLeading - (mcPart->Phi());
684
685 //plot right dphi
686 if(dPhi<=-(TMath::Pi())/2) dPhi = dPhi+2*(TMath::Pi());
687 if(dPhi>(3*(TMath::Pi()))/2) dPhi = dPhi-2*(TMath::Pi());
688
689 zMC2 = (TMath::Cos(dPhi)*(mcPart->P()))/jetEnergy;
690
691 return zMC2;
692}