]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/AliAnalysisTaskSEDStarJets.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix in the computetion of <pt>
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliAnalysisTaskSEDStarJets.cxx
CommitLineData
954ac830 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 **************************************************************************/
954ac830 15//
954ac830 16//
dee71c6f 17// Base class for DStar in Jets Analysis
954ac830 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>
dee71c6f 27#include "TROOT.h"
954ac830 28
29#include "AliAnalysisTaskSEDStarJets.h"
dee71c6f 30#include "AliRDHFCutsDStartoKpipi.h"
954ac830 31#include "AliStack.h"
32#include "AliMCEvent.h"
33#include "AliAODMCHeader.h"
954ac830 34#include "AliAODHandler.h"
dee71c6f 35#include "AliAnalysisManager.h"
954ac830 36#include "AliLog.h"
37#include "AliAODVertex.h"
38#include "AliAODJet.h"
39#include "AliAODRecoDecay.h"
dee71c6f 40#include "AliAODRecoCascadeHF.h"
954ac830 41#include "AliAODRecoDecayHF2Prong.h"
42#include "AliESDtrack.h"
43#include "AliAODMCParticle.h"
44
45ClassImp(AliAnalysisTaskSEDStarJets)
46
47//__________________________________________________________________________
48AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets() :
49 AliAnalysisTaskSE(),
954ac830 50 fEvents(0),
dee71c6f 51 fchargeFrCorr(0),
52 fUseMCInfo(kTRUE),
954ac830 53 fRequireNormalization(kTRUE),
954ac830 54 fOutput(0),
dee71c6f 55 fCuts(0),
954ac830 56 ftrigger(0),
57 fPtPion(0),
58 fInvMass(0),
dee71c6f 59 fRECOPtDStar(0),
60 fRECOPtBkg(0),
954ac830 61 fDStar(0),
62 fDiff(0),
63 fDiffSideBand(0),
64 fDStarMass(0),
65 fPhi(0),
66 fPhiBkg(0),
67 fTrueDiff(0),
68 fResZ(0),
dee71c6f 69 fResZBkg(0),
954ac830 70 fEjet(0),
71 fPhijet(0),
72 fEtaJet(0),
dee71c6f 73 theMCFF(0),
74 fDphiD0Dstar(0),
75 fPtJet(0)
954ac830 76{
77 //
78 // Default ctor
79 //
80}
81//___________________________________________________________________________
dee71c6f 82AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets(const Char_t* name, AliRDHFCutsDStartoKpipi* cuts) :
954ac830 83 AliAnalysisTaskSE(name),
954ac830 84 fEvents(0),
dee71c6f 85 fchargeFrCorr(0),
feaae220 86 fUseMCInfo(kTRUE),
954ac830 87 fRequireNormalization(kTRUE),
954ac830 88 fOutput(0),
dee71c6f 89 fCuts(0),
954ac830 90 ftrigger(0),
91 fPtPion(0),
92 fInvMass(0),
dee71c6f 93 fRECOPtDStar(0),
94 fRECOPtBkg(0),
954ac830 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),
954ac830 104 fEjet(0),
105 fPhijet(0),
106 fEtaJet(0),
dee71c6f 107 theMCFF(0),
108 fDphiD0Dstar(0),
109 fPtJet(0)
954ac830 110{
111 //
112 // Constructor. Initialization of Inputs and Outputs
113 //
dee71c6f 114 fCuts=cuts;
954ac830 115 Info("AliAnalysisTaskSEDStarJets","Calling Constructor");
116
dee71c6f 117 DefineOutput(1,TList::Class()); // histos
118 DefineOutput(2,AliRDHFCutsDStartoKpipi::Class()); // my cuts
954ac830 119}
954ac830 120//___________________________________________________________________________
dee71c6f 121AliAnalysisTaskSEDStarJets::~AliAnalysisTaskSEDStarJets() {
954ac830 122 //
dee71c6f 123 // destructor
954ac830 124 //
954ac830 125
954ac830 126 Info("~AliAnalysisTaskSEDStarJets","Calling Destructor");
127 if (fOutput) {
128 delete fOutput;
129 fOutput = 0;
dee71c6f 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;
954ac830 170}
171
172//_________________________________________________
173void AliAnalysisTaskSEDStarJets::UserExec(Option_t *)
174{
175 // user exec
176 if (!fInputEvent) {
177 Error("UserExec","NO EVENT FOUND!");
178 return;
179 }
dee71c6f 180
181 fEvents++;
182 AliInfo(Form("Event %d",fEvents));
183 if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
184
954ac830 185 // Load the event
954ac830 186 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
dee71c6f 187
188 TClonesArray *arrayDStartoD0pi=0;
189
954ac830 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();
dee71c6f 201 arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
954ac830 202 }
203 } else {
dee71c6f 204 arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
954ac830 205 }
206
dee71c6f 207 if (!arrayDStartoD0pi){
954ac830 208 AliInfo("Could not find array of HF vertices, skipping the event");
209 return;
dee71c6f 210 }else AliDebug(2, Form("Found %d vertices",arrayDStartoD0pi->GetEntriesFast()));
954ac830 211
7c23877d 212 // fix for temporary bug in ESDfilter
213 // the AODs with null vertex pointer didn't pass the PhysSel
5806c290 214 if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
7c23877d 215
954ac830 216 // Simulate a jet triggered sample
217 TClonesArray *arrayofJets = (TClonesArray*)aodEvent->GetJets();
218 if(aodEvent->GetNJets()<=0) return;
c088a94c 219
954ac830 220 // counters for efficiencies
954ac830 221 Int_t icountReco = 0;
48815220 222
954ac830 223 // Normalization factor
224 if(fRequireNormalization){
225 ftrigger->Fill(1);
226 }
dee71c6f 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;
954ac830 237
dee71c6f 238 //loop over jets and consider only the leading jey in the event
239 for (Int_t iJets = 0; iJets<arrayofJets->GetEntriesFast(); iJets++) {
954ac830 240 AliAODJet* jet = (AliAODJet*)arrayofJets->At(iJets);
dee71c6f 241
954ac830 242 //jets variables
dee71c6f 243 ejet = jet->E();
954ac830 244
dee71c6f 245 if(ejet>max){
246 max = jet->E();
247 phiJet = jet->Phi();
248 etaJet = jet->Eta();
249 ptjet = jet->Pt();
250 }
251
954ac830 252 // fill energy, eta and phi of the jet
253 fEjet ->Fill(ejet);
254 fPhijet ->Fill(phiJet);
255 fEtaJet ->Fill(etaJet);
dee71c6f 256 fPtJet->Fill(ptjet);
257 }
258
259 //loop over D* candidates
260 for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
954ac830 261
dee71c6f 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 }
954ac830 290
dee71c6f 291 // soft pion
292 AliAODTrack *track2 = (AliAODTrack*)dstarD0pi->GetBachelor();
293 Double_t pt = dstarD0pi->Pt();
954ac830 294
dee71c6f 295 // track quality cuts
296 Int_t isTkSelected = fCuts->IsSelected(dstarD0pi,AliRDHFCuts::kTracks); // quality cuts on tracks
297 if(!isTkSelected) continue;
954ac830 298
dee71c6f 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);
954ac830 307
dee71c6f 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;
954ac830 319
dee71c6f 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();
954ac830 343
dee71c6f 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
954ac830 347
dee71c6f 348 fRECOPtDStar->Fill(pt);
349 fPtPion->Fill(track2->Pt());
954ac830 350
dee71c6f 351 fPhi ->Fill(dPhi);
954ac830 352
dee71c6f 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);
954ac830 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
dee71c6f 377 Info("Terminate"," terminate");
954ac830 378 AliAnalysisTaskSE::Terminate();
954ac830 379
380 fOutput = dynamic_cast<TList*> (GetOutputData(1));
381 if (!fOutput) {
382 printf("ERROR: fOutput not available\n");
383 return;
384 }
385
954ac830 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"));
dee71c6f 395 fRECOPtBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fRECOPtBkg"));
954ac830 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"));
dee71c6f 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
954ac830 407}
408//___________________________________________________________________________
409
410void AliAnalysisTaskSEDStarJets::UserCreateOutputObjects() {
dee71c6f 411 // output
954ac830 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}
954ac830 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
dee71c6f 470 fDphiD0Dstar = new TH1F("phiD0Dstar","Delta phi between D0 and DStar ",1000,-6.5,6.5);
471 fOutput->Add(fDphiD0Dstar);
472
954ac830 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
dee71c6f 479 fRECOPtDStar = new TH1F("RECODStar1","RECO DStar pt distribution",30,0,30);
954ac830 480 fRECOPtDStar->SetStats(kTRUE);
481 fRECOPtDStar->SetLineColor(2);
482 fRECOPtDStar->GetXaxis()->SetTitle("GeV/c");
483 fRECOPtDStar->GetYaxis()->SetTitle("Entries");
484 fOutput->Add(fRECOPtDStar);
dee71c6f 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);
954ac830 492
dee71c6f 493 fPtPion = new TH1F("pionpt","Primary pions candidates pt ",500,0,10);
954ac830 494 fPtPion->SetStats(kTRUE);
495 fPtPion->GetXaxis()->SetTitle("GeV/c");
496 fPtPion->GetYaxis()->SetTitle("Entries");
497 fOutput->Add(fPtPion);
dee71c6f 498
954ac830 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);
dee71c6f 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);
954ac830 504 fOutput->Add(fEjet);
505 fOutput->Add(fPhijet);
506 fOutput->Add(fEtaJet);
dee71c6f 507 fOutput->Add(fPtJet);
508
509 theMCFF = new TH1F("FragFuncMC","Fragmentation function in MC for FC ",100,0,10);
954ac830 510 fResZ = new TH1F("FragFunc","Fragmentation function ",50,0,1);
511 fResZBkg = new TH1F("FragFuncBkg","Fragmentation function background",50,0,1);
dee71c6f 512 fOutput->Add(theMCFF);
954ac830 513 fOutput->Add(fResZ);
514 fOutput->Add(fResZBkg);
515
dee71c6f 516 return kTRUE;
954ac830 517}
518
519//______________________________ side band background for D*___________________________________
520
dee71c6f 521void AliAnalysisTaskSEDStarJets::SideBandBackground(Double_t invM, Double_t invMDStar, Double_t dStarMomBkg, Double_t EGjet, Double_t dPhi){
954ac830 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
dee71c6f 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
c6d629dd 555 if (!mcPart) return zMC2;
556 if (!mcArray) return zMC2;
557
dee71c6f 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;
954ac830 569
dee71c6f 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();
954ac830 577
dee71c6f 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();
954ac830 595 }
596 }
dee71c6f 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
954ac830 612 }
dee71c6f 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;
954ac830 692}