]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/AliAnalysisTaskSEDStarJets.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added possiblity to run without MC info (Alessandro)
[u/mrichter/AliRoot.git] / PWG3 / 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 **************************************************************************/
15//
16//
17// Base class for DStar in Jets Analysis
18//
19// The D* (+ and -) is reconstructed inside jets. Two different cuts are
20// implemented:
21//
22// 1) C.Ivan D* cuts adapted for correlation analysis
23// 2) Topological cut enforcing the correlation D0-softPion pt + relaxed
24// CosThetaStar. This second should be better for correlations.
25//
26// USAGE:
27//
28// The analysis is performed separately for D*+ and D*-. A flag in the .C is
29// taking care to decide which analysis.
30//
31// The cut number 2 can be activaded with a flag in the .C (not active in this version)
32// Cuts 1) are the default.
33//
34// The task requires reconstructed jets in the AODs
35//
36//-----------------------------------------------------------------------
37// Author A.Grelli
38// ERC-QGP Utrecht University - a.grelli@uu.nl
39//-----------------------------------------------------------------------
40
41#include <TDatabasePDG.h>
42#include <TParticle.h>
43#include <TVector3.h>
44
45#include "AliAnalysisTaskSEDStarJets.h"
46#include "AliStack.h"
47#include "AliMCEvent.h"
48#include "AliAODMCHeader.h"
49#include "AliAnalysisManager.h"
50#include "AliAODHandler.h"
51#include "AliLog.h"
52#include "AliAODVertex.h"
53#include "AliAODJet.h"
54#include "AliAODRecoDecay.h"
55#include "AliAODRecoDecayHF.h"
56#include "AliAODRecoDecayHF2Prong.h"
57#include "AliESDtrack.h"
58#include "AliAODMCParticle.h"
59
60ClassImp(AliAnalysisTaskSEDStarJets)
61
62//__________________________________________________________________________
63AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets() :
64 AliAnalysisTaskSE(),
954ac830 65 fCountReco(0),
66 fCountRecoAcc(0),
67 fCountRecoITSClusters(0),
68 fCountRecoPPR(0),
69 fCountDStar(0),
954ac830 70 fEvents(0),
71 fMinITSClusters(0),
72 fComputeD0(kTRUE),
feaae220 73 fUseMCInfo(kTRUE),
954ac830 74 ftopologicalCut(kFALSE),
75 fRequireNormalization(kTRUE),
76 fLorentzTrack1(0,0,0,0),
77 fLorentzTrack2(0,0,0,0),
78 fLorentzTrack3(0,0,0,0),
79 fLorentzTrack4(0,0,0,0),
80 fOutput(0),
81 fD0ptvsSoftPtSignal(0),
82 fD0ptvsSoftPt(0),
83 ftrigger(0),
84 fPtPion(0),
85 fInvMass(0),
86 fRECOPtDStar(0),
87 fDStar(0),
88 fDiff(0),
89 fDiffSideBand(0),
90 fDStarMass(0),
91 fPhi(0),
92 fPhiBkg(0),
93 fTrueDiff(0),
94 fResZ(0),
95 fResZBkg(0),
96 fcharmpt(0),
97 fdstarE(0),
98 fEjet(0),
99 fPhijet(0),
100 fEtaJet(0),
c088a94c 101 fdstarpt(0)
954ac830 102{
103 //
104 // Default ctor
105 //
106}
107//___________________________________________________________________________
108AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets(const Char_t* name) :
109 AliAnalysisTaskSE(name),
954ac830 110 fCountReco(0),
111 fCountRecoAcc(0),
112 fCountRecoITSClusters(0),
113 fCountRecoPPR(0),
114 fCountDStar(0),
954ac830 115 fEvents(0),
116 fMinITSClusters(0),
117 fComputeD0(kTRUE),
feaae220 118 fUseMCInfo(kTRUE),
954ac830 119 ftopologicalCut(kFALSE),
120 fRequireNormalization(kTRUE),
121 fLorentzTrack1(0,0,0,0),
122 fLorentzTrack2(0,0,0,0),
123 fLorentzTrack3(0,0,0,0),
124 fLorentzTrack4(0,0,0,0),
125 fOutput(0),
126 fD0ptvsSoftPtSignal(0),
127 fD0ptvsSoftPt(0),
128 ftrigger(0),
129 fPtPion(0),
130 fInvMass(0),
131 fRECOPtDStar(0),
132 fDStar(0),
133 fDiff(0),
134 fDiffSideBand(0),
135 fDStarMass(0),
136 fPhi(0),
137 fPhiBkg(0),
138 fTrueDiff(0),
139 fResZ(0),
140 fResZBkg(0),
141 fcharmpt(0),
142 fdstarE(0),
143 fEjet(0),
144 fPhijet(0),
145 fEtaJet(0),
c088a94c 146 fdstarpt(0)
954ac830 147{
148 //
149 // Constructor. Initialization of Inputs and Outputs
150 //
151 Info("AliAnalysisTaskSEDStarJets","Calling Constructor");
152
153 DefineOutput(1,TList::Class());
154}
155
156//___________________________________________________________________________
157AliAnalysisTaskSEDStarJets& AliAnalysisTaskSEDStarJets::operator=(const AliAnalysisTaskSEDStarJets& c)
158{
159 //
160 // Assignment operator
161 //
162 if (this!=&c) {
163 AliAnalysisTaskSE::operator=(c) ;
164 }
165
166 return *this;
167}
168
169//___________________________________________________________________________
170AliAnalysisTaskSEDStarJets::AliAnalysisTaskSEDStarJets(const AliAnalysisTaskSEDStarJets& c) :
171 AliAnalysisTaskSE(c),
954ac830 172 fCountReco(c.fCountReco),
173 fCountRecoAcc(c.fCountRecoAcc),
174 fCountRecoITSClusters(c.fCountRecoITSClusters),
175 fCountRecoPPR(c.fCountRecoPPR),
176 fCountDStar(c.fCountDStar),
954ac830 177 fEvents(c.fEvents),
178 fMinITSClusters(c.fMinITSClusters),
179 fComputeD0(c.fComputeD0),
feaae220 180 fUseMCInfo(c.fUseMCInfo),
954ac830 181 ftopologicalCut(c.ftopologicalCut),
182 fRequireNormalization(c.fRequireNormalization),
183 fLorentzTrack1(c.fLorentzTrack1),
184 fLorentzTrack2(c.fLorentzTrack2),
185 fLorentzTrack3(c.fLorentzTrack3),
186 fLorentzTrack4(c.fLorentzTrack4),
187 fOutput(c.fOutput),
188 fD0ptvsSoftPtSignal(c.fD0ptvsSoftPtSignal),
189 fD0ptvsSoftPt(c.fD0ptvsSoftPt),
190 ftrigger(c.ftrigger),
191 fPtPion(c.fPtPion),
192 fInvMass(c.fInvMass),
193 fRECOPtDStar(c.fRECOPtDStar),
194 fDStar(c.fDStar),
195 fDiff(c.fDiff),
196 fDiffSideBand(c.fDiffSideBand),
197 fDStarMass(c.fDStarMass),
198 fPhi(c.fPhi),
199 fPhiBkg(c.fPhiBkg),
200 fTrueDiff(c.fTrueDiff),
201 fResZ(c.fResZ),
202 fResZBkg(c.fResZBkg),
203 fcharmpt(c.fcharmpt),
204 fdstarE(c.fdstarE),
205 fEjet(c.fEjet),
206 fPhijet(c.fPhijet),
207 fEtaJet(c.fEtaJet),
c088a94c 208 fdstarpt(c.fdstarpt)
954ac830 209
210{
211 //
212 // Copy Constructor
213 //
214}
215
216//___________________________________________________________________________
217AliAnalysisTaskSEDStarJets::~AliAnalysisTaskSEDStarJets() {
218 // destructor
219
220 Info("~AliAnalysisTaskSEDStarJets","Calling Destructor");
221 if (fOutput) {
222 delete fOutput;
223 fOutput = 0;
224 }
225}
226
227//_________________________________________________
228void AliAnalysisTaskSEDStarJets::UserExec(Option_t *)
229{
230 // user exec
231 if (!fInputEvent) {
232 Error("UserExec","NO EVENT FOUND!");
233 return;
234 }
235
236 // Load the event
237 fEvents++;
238 AliInfo(Form("Event %d",fEvents));
239 if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
240 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
241
242 TClonesArray *arrayVerticesHF=0;
243
244 if(!aodEvent && AODEvent() && IsStandardAOD()) {
245 // In case there is an AOD handler writing a standard AOD, use the AOD
246 // event in memory rather than the input (ESD) event.
247 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
248 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
249 // have to taken from the AOD event hold by the AliAODExtension
250 AliAODHandler* aodHandler = (AliAODHandler*)
251 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
252 if(aodHandler->GetExtensions()) {
253 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
254 AliAODEvent *aodFromExt = ext->GetAOD();
255 arrayVerticesHF=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
256 }
257 } else {
258 arrayVerticesHF=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
259 }
260
261 if (!arrayVerticesHF){
262 AliInfo("Could not find array of HF vertices, skipping the event");
263 return;
264 }else AliDebug(2, Form("Found %d vertices",arrayVerticesHF->GetEntriesFast()));
265
266 // Simulate a jet triggered sample
267 TClonesArray *arrayofJets = (TClonesArray*)aodEvent->GetJets();
268 if(aodEvent->GetNJets()<=0) return;
269 AliInfo("found a jet: processing D* in jet analysis");
954ac830 270
c088a94c 271 // AOD primary vertex
272 AliAODVertex *prVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
273 Double_t primaryPos[3];
274 prVtx->GetXYZ(primaryPos);
275 Bool_t vtxFlag = kTRUE;
276 TString title=prVtx->GetTitle();
277 if(!title.Contains("VertexerTracks")) vtxFlag=kFALSE;
278
954ac830 279 // counters for efficiencies
954ac830 280 Int_t icountReco = 0;
281 Int_t icountRecoAcc = 0;
282 Int_t icountRecoITSClusters = 0;
283 Int_t icountRecoPPR = 0;
284 Int_t fiDstar = 0;
285 Int_t fDStarD0 = 0;
feaae220 286
287 // TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
288
289 if(fUseMCInfo){
290 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
291 //loop on the MC event - some basic MC info on D*, D0 and soft pion
292 if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
954ac830 293
feaae220 294 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
295 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
296 if (!mcPart) {
297 AliWarning("Particle not found in tree, skipping");
298 continue;
299 }
c088a94c 300
feaae220 301 // charm pt
302 if(TMath::Abs(mcPart->GetPdgCode())==4){
303 fcharmpt->Fill(mcPart->Pt());
c088a94c 304 }
305
feaae220 306 // fill energy and pt for D* in acceptance with correct prongs
307 Bool_t isOk = DstarInMC(mcPart,mcArray);
c088a94c 308
feaae220 309 if (isOk){ //D*
310 AliDebug(2, "Found a DStar in MC with correct prongs and in acceptance");
311 fdstarE ->Fill(mcPart->E());
312 fdstarpt->Fill(mcPart->Pt());
954ac830 313
feaae220 314 // check the MC-Acceptance level cuts
315 // since standard CF functions are not applicable, using Kine Cuts on daughters
316
317 Int_t daughter0 = mcPart->GetDaughter(0);
318 Int_t daughter1 = mcPart->GetDaughter(1);
319
320 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
321
322 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
323 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
324
325 Double_t eta0 = mcPartDaughter0->Eta();
326 Double_t eta1 = mcPartDaughter1->Eta();
327 Double_t y0 = mcPartDaughter0->Y();
328 Double_t y1 = mcPartDaughter1->Y();
329 Double_t pt0 = mcPartDaughter0->Pt();
330 Double_t pt1 = mcPartDaughter1->Pt();
331
332 AliDebug(2, Form("D* Daughter 0: eta = %f, y = %f, pt = %f", eta0, y0, pt0));
333 AliDebug(2, Form("D* Daughter 1: eta = %f, y = %f, pt = %f", eta1, y1, pt1));
334
335 Int_t daughD00 = 0;
336 Int_t daughD01 = 0;
337
338 // D0 daughters - do not need to check D0-kpi, already done
339
340 if(TMath::Abs(mcPartDaughter0->GetPdgCode())==421){
341 daughD00 = mcPartDaughter0->GetDaughter(0);
342 daughD01 = mcPartDaughter0->GetDaughter(1);
343 }else{
344 daughD00 = mcPartDaughter1->GetDaughter(0);
345 daughD01 = mcPartDaughter1->GetDaughter(1);
346 }
347
348 AliAODMCParticle* mcD0PartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughD00));
349 AliAODMCParticle* mcD0PartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughD01));
350
351 if (!mcD0PartDaughter0 || !mcD0PartDaughter1) {
352 AliWarning("At least one Daughter Particle not found in tree, but it should be, this check was already done...");
353 }
354
355 // D0 daughters - needed for acceptance
356 Double_t pD0pt0 = mcD0PartDaughter0->Pt();
357 Double_t pD0pt1 = mcD0PartDaughter1->Pt();
358 Double_t pD0eta0 = mcD0PartDaughter0->Eta();
359 Double_t pD0eta1 = mcD0PartDaughter1->Eta();
360
361 // ACCEPTANCE REQUESTS ---------
362
363 // soft pion
364 Bool_t daught1inAcceptance = (TMath::Abs(eta1) <= 0.9 && pt1 > 0.05);
365 // Do daughters
366 Bool_t theD0daught0inAcceptance = (TMath::Abs(pD0eta0) <= 0.9 && pD0pt0 >= 0.1);
367 Bool_t theD0daught1inAcceptance = (TMath::Abs(pD0eta1) <= 0.9 && pD0pt1 >= 0.1);
368
369 if (daught1inAcceptance && theD0daught0inAcceptance && theD0daught1inAcceptance) {
370
371 AliDebug(2, "Daughter particles in acceptance");
372
373 // check on the vertex
374 if (vtxFlag){
375 printf("Vertex cut passed 2\n");
376 fDStarD0++;
377 Bool_t refitFlag = kTRUE;
378 for (Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
379 AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
380
381 // refit only for D0 daughters
382 if ((track->GetLabel() == daughD00) || (track->GetLabel() == daughD01)) {
383 if(!(track->GetStatus()&AliESDtrack::kTPCrefit) || !(track->GetStatus()&AliESDtrack::kITSrefit)) {
384 refitFlag = kFALSE;
385 }
c088a94c 386 }
387 }
feaae220 388 if (refitFlag){
389 printf("Refit cut passed\n");
390 }
391 else{
392 AliDebug(3,"Refit cut not passed\n");
393 }
c088a94c 394 }
395 else{
feaae220 396 AliDebug(3,"Vertex cut not passed\n");
397 }
954ac830 398 }
c088a94c 399 else{
feaae220 400 AliDebug(3,"Acceptance cut not passed\n");
401 }
954ac830 402 }
954ac830 403 }
feaae220 404
405 AliDebug(2, Form("Found %i MC particles that are D* in Kpipi and satisfy Acc cuts!!",fDStarD0));
406 fCountDStar += fDStarD0;
407
408 } // End of MC
954ac830 409
410 // Now perform the D* in jet reconstruction
feaae220 411
954ac830 412 // Normalization factor
413 if(fRequireNormalization){
414 ftrigger->Fill(1);
415 }
416
417 Int_t nJets = 0; // for reco D0 countings
418
419 for (Int_t iJets = 0; iJets<arrayofJets->GetEntriesFast(); iJets++) { // main loop on jets
420
421 AliAODJet* jet = (AliAODJet*)arrayofJets->At(iJets);
422
423 //jets variables
424 Double_t ejet = jet->E();
425 Double_t phiJet = jet->Phi();
426 Double_t etaJet = jet->Eta();
427
428 // fill energy, eta and phi of the jet
429 fEjet ->Fill(ejet);
430 fPhijet ->Fill(phiJet);
431 fEtaJet ->Fill(etaJet);
432
433 nJets++;
434
435 // loop over the tracks to search for candidates soft pions
436 for (Int_t i=0; i<aodEvent->GetNTracks(); i++){
437
438 AliAODTrack* aodTrack = aodEvent->GetTrack(i);
954ac830 439 Double_t vD0[4] = {0.,0.,0.,0.};
440 Double_t vD0bar[4] ={0.,0.,0.,0.};
441
442 Int_t pCharge = aodTrack->Charge();
443
444 // few selections on soft pion
445 Bool_t tPrimCand = aodTrack->IsPrimaryCandidate(); // is it primary?
446
447 if(aodTrack->Pt()>= 5 || aodTrack->Pt()< 0.08 ) continue; //cut on soft pion pt VERY large
448 //~ D*s of pt >80GeV with a soft pion of 5GeV!
449 if(TMath::Abs(aodTrack->Eta())>0.9) continue;
450
c088a94c 451 // if D*+ analysis tha D0 and pi+ otherwise pi-
954ac830 452 if(fComputeD0 && pCharge!= 1 ) continue;
954ac830 453 if(!fComputeD0 && pCharge!= -1 ) continue;
454
455 if(tPrimCand && arrayVerticesHF->GetEntriesFast()>0){ // isPion and is Primary, no PID for now
456
457 // label to the candidate soft pion
feaae220 458 Int_t pLabel = -1;
459 if(fUseMCInfo) pLabel = aodTrack->GetLabel();
954ac830 460
461 // prepare the TLorentz vector for the pion
462 Float_t pionMass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
463 fLorentzTrack3.SetPxPyPzE(aodTrack->Px(),aodTrack->Py(),aodTrack->Pz(),aodTrack->E(pionMass));
464
465 // search for the D0
466 for (Int_t iVertex = 0; iVertex<arrayVerticesHF->GetEntriesFast(); iVertex++) {
467
468 Double_t invM = 0;
469 Double_t invMDStar = 0;
c088a94c 470 Double_t dPhi = 0;
954ac830 471
472 AliAODRecoDecayHF2Prong* vtx = (AliAODRecoDecayHF2Prong*)arrayVerticesHF->At(iVertex);
473
474 Double_t pt = vtx->Pt();
c088a94c 475
476 // acceptance variables
477 Bool_t acceptanceProng0 = (TMath::Abs(vtx->EtaProng(0))<= 0.9 && vtx->PtProng(0) >= 0.1);
478 Bool_t acceptanceProng1 = (TMath::Abs(vtx->EtaProng(1))<= 0.9 && vtx->PtProng(1) >= 0.1);
479 Bool_t acceptanceSoftPi = (TMath::Abs(aodTrack->Eta())<= 0.9 && aodTrack->Pt() >= 0.05);
954ac830 480
954ac830 481 Int_t pdgDgD0toKpi[2]={321,211};
c088a94c 482
483 Int_t mcLabel =-1;
c088a94c 484
485 Bool_t isDStar = kFALSE; // just to count
486
feaae220 487 //switch of MC for real data
488 if(fUseMCInfo){
489 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
490 mcLabel = vtx->MatchToMC(421, mcArray,2,pdgDgD0toKpi) ; //MC D0
491 // matching to MC D*
492 if(mcLabel !=-1 && pLabel!=-1 && nJets ==1) { // count only once in case of multijets
493
494 // search for a D0 and a pi with mother and check it is a D*
495 AliAODMCParticle* theMCpion = (AliAODMCParticle*)mcArray->At(pLabel);
496 Int_t motherMCPion = theMCpion->GetMother();
497 AliAODMCParticle* theMCD0 = (AliAODMCParticle*)mcArray->At(mcLabel);
498 Int_t motherMCD0 = theMCD0->GetMother();
499
500 if(motherMCPion!=-1 && (motherMCD0 == motherMCPion)){
501 AliAODMCParticle* mcMother = (AliAODMCParticle*)mcArray->At(motherMCPion);
502 if(TMath::Abs(mcMother->GetPdgCode()) == 413 && TMath::Abs(theMCpion->GetPdgCode())==211){
503 isDStar = kTRUE;
504 fiDstar++;
505 }
c088a94c 506 }
507 }
508 }
954ac830 509
feaae220 510
c088a94c 511 if (acceptanceProng0 && acceptanceProng1 && acceptanceSoftPi) {
954ac830 512
513 AliDebug(2,"D0 reco daughters in acceptance");
c088a94c 514 if(isDStar && nJets==1) icountRecoAcc++;
954ac830 515
c088a94c 516 // D0 daughters
954ac830 517 AliAODTrack *track0 = (AliAODTrack*)vtx->GetDaughter(0);
518 AliAODTrack *track1 = (AliAODTrack*)vtx->GetDaughter(1);
519
520 // check for ITS refit (already required at the ESD filter level )
c088a94c 521 Bool_t kRefitITS = kTRUE;
954ac830 522
523 if((!(track0->GetStatus()&AliESDtrack::kITSrefit)|| (!(track1->GetStatus()&AliESDtrack::kITSrefit)))) {
524 kRefitITS = kFALSE;
c088a94c 525 }
954ac830 526
c088a94c 527 Int_t ncls0=0,ncls1=0,ncls2=0;
528
954ac830 529 for(Int_t l=0;l<6;l++) {
530 if(TESTBIT(track0->GetITSClusterMap(),l)) ncls0++;
531 if(TESTBIT(track1->GetITSClusterMap(),l)) ncls1++;
c088a94c 532 if(TESTBIT(aodTrack->GetITSClusterMap(),l)) ncls2++;
533 }
534
535 // clusters in ITS for D0 daugthers and soft pion
feaae220 536 if (ncls0 >= fMinITSClusters && ncls1 >= fMinITSClusters && ncls2>=4 && kRefitITS) {
954ac830 537
c088a94c 538 if(isDStar && nJets==1) icountRecoITSClusters++;
954ac830 539
540 // D0 cutting varibles
541 Double_t cosThetaStar = 9999.;
542 Double_t pTpi = 0.;
543 Double_t pTK = 0.;
c088a94c 544 Double_t d01 = 0;
545 Double_t d00 = 0;
954ac830 546
547 // D0, D0bar
c088a94c 548 if (fComputeD0){
549 cosThetaStar = vtx->CosThetaStarD0();
954ac830 550 pTpi = vtx->PtProng(0);
c088a94c 551 pTK = vtx->PtProng(1);
552 d01 = vtx->Getd0Prong(0)*1E4;
553 d00 = vtx->Getd0Prong(1)*1E4;
554 }else{
555 cosThetaStar = vtx->CosThetaStarD0bar();
954ac830 556 pTpi = vtx->PtProng(1);
c088a94c 557 pTK = vtx->PtProng(0);
558 d01 = vtx->Getd0Prong(1)*1E4;
559 d00 = vtx->Getd0Prong(0)*1E4;
954ac830 560 }
c088a94c 561
954ac830 562 AliDebug(2,"D0 reco daughters in acceptance");
c088a94c 563
954ac830 564 Double_t dca = vtx->GetDCA()*1E4;
954ac830 565 Double_t d0d0 = vtx->Prodd0d0()*1E8;
c088a94c 566
954ac830 567 TVector3 mom(vtx->Px(),vtx->Py(),vtx->Pz());
568 TVector3 flight((vtx->Xv())- primaryPos[0],(vtx->Yv())- primaryPos[1],(vtx->Zv())- primaryPos[2]);
569 Double_t pta = mom.Angle(flight);
570 Double_t cosPointingAngle = TMath::Cos(pta);
c088a94c 571
572 // Crstian Ivan D* cuts. Multidimensional optimization
573 Double_t cuts[7] = {9999999., 1.1, 0., 9999999.,999999., 9999999., 0.};
954ac830 574
575 if (pt <= 1){ // first bin not optimized
576 cuts[0] = 400;
577 cuts[1] = 0.8;
578 cuts[2] = 0.21;
579 cuts[3] = 500;
c088a94c 580 cuts[4] = 500;
581 cuts[5] = -20000;
582 cuts[6] = 0.6;
954ac830 583 }
584 else if (pt > 1 && pt <= 2){
c088a94c 585 cuts[0] = 200;
954ac830 586 cuts[1] = 0.7;
587 cuts[2] = 0.8;
c088a94c 588 cuts[3] = 210;
589 cuts[4] = 210;
590 cuts[5] = -20000;
591 cuts[6] = 0.9;
954ac830 592 }
593 else if (pt > 2 && pt <= 3){
594 cuts[0] = 400;
595 cuts[1] = 0.8;
596 cuts[2] = 0.8;
c088a94c 597 cuts[3] = 420;
598 cuts[4] = 350;
599 cuts[5] = -8500;
600 cuts[6] = 0.9;
954ac830 601 }
602 else if (pt > 3 && pt <= 5){
c088a94c 603 cuts[0] = 160;
604 cuts[1] = 1.0;
954ac830 605 cuts[2] = 1.2;
c088a94c 606 cuts[3] = 420;
607 cuts[4] = 560;
608 cuts[5] = -8500;
609 cuts[6] = 0.9;
954ac830 610 }
611 else if (pt > 5){
612 cuts[0] = 800;
613 cuts[1] = 1.0;
614 cuts[2] = 1.2;
c088a94c 615 cuts[3] = 700;
616 cuts[4] = 700;
617 cuts[5] = 10000;
618 cuts[6] = 0.9;
954ac830 619 }
620 // apply D0 cuts
621
622 if (dca < cuts[0]
c088a94c 623 && TMath::Abs(cosThetaStar) < cuts[1]
624 && pTpi > cuts[2]
625 && pTK > cuts[2]
626 && TMath::Abs(d01) < cuts[3]
627 && TMath::Abs(d00) < cuts[4]
628 && d0d0 < cuts[5]
629 && cosPointingAngle > cuts[6]
630 ){
954ac830 631
632 if(fComputeD0){ // D0 from default
633
634 if(vtx->ChargeProng(0)==-1){
635 fLorentzTrack1.SetPxPyPzE( vtx->PxProng(0),vtx->PyProng(0), vtx->PzProng(0),vtx->EProng(0,321) );
636 fLorentzTrack2.SetPxPyPzE( vtx->PxProng(1),vtx->PyProng(1), vtx->PzProng(1),vtx->EProng(1,211) );
637 }else{
638 fLorentzTrack1.SetPxPyPzE( vtx->PxProng(0),vtx->PyProng(0), vtx->PzProng(0),vtx->EProng(0,211) );
639 fLorentzTrack2.SetPxPyPzE( vtx->PxProng(1),vtx->PyProng(1), vtx->PzProng(1),vtx->EProng(1,321) );
640 }
641
954ac830 642 vD0[0] = (fLorentzTrack1+fLorentzTrack2).Px();
643 vD0[1] = (fLorentzTrack1+fLorentzTrack2).Py();
644 vD0[2] = (fLorentzTrack1+fLorentzTrack2).Pz();
645 vD0[3] = (fLorentzTrack1+fLorentzTrack2).E();
646
647 fLorentzTrack4.SetPxPyPzE(vD0[0],vD0[1],vD0[2],vD0[3]);
648
649 }else{ // D0bar analysis
650
651 if(vtx->ChargeProng(0)==-1){
652 fLorentzTrack1.SetPxPyPzE( vtx->PxProng(0),vtx->PyProng(0), vtx->PzProng(0),vtx->EProng(0,211) );
653 fLorentzTrack2.SetPxPyPzE( vtx->PxProng(1),vtx->PyProng(1), vtx->PzProng(1),vtx->EProng(1,321) );
654 }else{
655 fLorentzTrack1.SetPxPyPzE( vtx->PxProng(0),vtx->PyProng(0), vtx->PzProng(0),vtx->EProng(0,321) );
656 fLorentzTrack2.SetPxPyPzE( vtx->PxProng(1),vtx->PyProng(1), vtx->PzProng(1),vtx->EProng(1,211) );
657 }
658
954ac830 659 vD0bar[0] = (fLorentzTrack1+fLorentzTrack2).Px();
660 vD0bar[1] = (fLorentzTrack1+fLorentzTrack2).Py();
661 vD0bar[2] = (fLorentzTrack1+fLorentzTrack2).Pz();
662 vD0bar[3] = (fLorentzTrack1+fLorentzTrack2).E();
663
664 fLorentzTrack4.SetPxPyPzE(vD0bar[0],vD0bar[1],vD0bar[2],vD0bar[3]);
665 }
666
667 // D0-D0bar related variables
668
669 invM = GetInvariantMass(fLorentzTrack1,fLorentzTrack2);
c088a94c 670 if(nJets==1) fInvMass->Fill(invM);
954ac830 671
c088a94c 672 if(isDStar && nJets==1) icountRecoPPR++;
673
954ac830 674 //DStar invariant mass
675 invMDStar = GetInvariantMassDStar(fLorentzTrack3,fLorentzTrack4);
676
677 //conversion from phi TLorentzVerctor to phi aliroot
678 Double_t kconvert = 0;
679 Double_t phiDStar = (fLorentzTrack3 + fLorentzTrack4).Phi();
680 kconvert = phiDStar;
681 if(phiDStar<0) kconvert = phiDStar + 2*(TMath::Pi());
682 phiDStar = kconvert;
683
684 // dphi between jet and D*
685 dPhi = phiJet - phiDStar;
686
687 //Just for plotting pourposal
688 if(dPhi<=-1.58) dPhi = TMath::Abs(dPhi);
689 if(dPhi>4.72){
690 dPhi = dPhi-2*(TMath::Pi());
691 }
692
693 //Alternative cutting procedure
694 //the cut on cosThetaStar is to reduce near collinear
695 //background from jet fragmentation
696
697 Bool_t tCut = EvaluateCutOnPiD0pt(vtx,aodTrack);
698
699 if(ftopologicalCut && tCut) continue;
700 if(ftopologicalCut && TMath::Abs(cosThetaStar)>cuts[1]) continue;
701
702 if(invM >= 1.829 && invM <= 1.901){ // ~4 sigma cut on D0 mass
703
c088a94c 704 if(isDStar && nJets==1) {
705 fDStarMass->Fill(invMDStar);
706 fTrueDiff->Fill(invMDStar-invM);
954ac830 707 }
708
c088a94c 709 if(nJets==1) { // avoid double counting
710 fDStar->Fill(invMDStar);
711 fDiff->Fill(invMDStar-invM); // M(Kpipi)-M(Kpi)
712 }
713
954ac830 714 // now the dphi signal and the fragmentation function
feaae220 715 if((invMDStar-invM)<=0.148 && (invMDStar-invM)>=0.142) {
954ac830 716
717 //fill candidates D* and soft pion reco pt
c088a94c 718 if(nJets==1) fPtPion->Fill(aodTrack->Pt());
719 if(nJets==1) fRECOPtDStar->Fill((fLorentzTrack3 + fLorentzTrack4).Pt());
954ac830 720 fPhi ->Fill(dPhi);
721
722 if(dPhi>=-0.5 && dPhi<=0.5){ // evaluate in the near side
723 Double_t dStarMom = (fLorentzTrack3 + fLorentzTrack4).P();
724 Double_t zFrag = (TMath::Cos(dPhi)*dStarMom)/ejet;
725 fResZ->Fill(TMath::Abs(zFrag));
726 }
727 }
728 }
729
730 // evaluate side band background
feaae220 731 SideBandBackground(invM, invMDStar, ejet, dPhi, nJets);
954ac830 732
733 invM = 0;
734 invMDStar = 0;
735
c088a94c 736 } // end PPR cuts
737 } // end ITS cluster
738 } // end acceptance
739 } // D0 cand
740 } // softpion
741 } // tracks
954ac830 742 } // jets
743
744 AliDebug(2, Form("Found %i Reco particles that are D0!!",icountReco));
745
746 fCountReco+= fiDstar;
747 fCountRecoAcc+= icountRecoAcc;
748 fCountRecoITSClusters+= icountRecoITSClusters;
749 fCountRecoPPR+= icountRecoPPR;
750
751 PostData(1,fOutput);
752}
753
754//________________________________________ terminate ___________________________
755void AliAnalysisTaskSEDStarJets::Terminate(Option_t*)
756{
757 // The Terminate() function is the last function to be called during
758 // a query. It always runs on the client, it can be used to present
759 // the results graphically or save the results to file.
760
761 Info("Terminate","");
762 AliAnalysisTaskSE::Terminate();
763
954ac830 764 AliInfo(Form("Found %i of that MC D*->D0pi(D0->kpi) are in acceptance, in %d events",fCountDStar,fEvents));
765
766 AliInfo(Form("Found %i RECO particles after cuts that are DStar, in %d events",fCountReco,fEvents));
c088a94c 767 AliInfo(Form("Among the above, found %i reco D* that are decaying in K+pi+pi and are in the requested acceptance, in %d events",fCountRecoAcc,fEvents));
768 AliInfo(Form("Among the above, found %i reco D* that are decaying in K+pi+pi and have at least %d clusters in ITS, in %d events",fCountRecoITSClusters,fMinITSClusters,fEvents));
769 AliInfo(Form("Among the above, found %i reco D* that are decaying in K+pi+pi and satisfy D* cuts, in %d events",fCountRecoPPR,fEvents));
954ac830 770
771 fOutput = dynamic_cast<TList*> (GetOutputData(1));
772 if (!fOutput) {
773 printf("ERROR: fOutput not available\n");
774 return;
775 }
776
954ac830 777 fcharmpt = dynamic_cast<TH1F*>(fOutput->FindObject("fcharmpt"));
778 fdstarE = dynamic_cast<TH1F*>(fOutput->FindObject("fdstarE"));
779 fdstarpt = dynamic_cast<TH1F*>(fOutput->FindObject("fdstarpt"));
780 fDStarMass = dynamic_cast<TH1F*>(fOutput->FindObject("fDStarMass"));
781 fTrueDiff = dynamic_cast<TH1F*>(fOutput->FindObject("fTrueDiff"));
782 fInvMass = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass"));
783 fPtPion = dynamic_cast<TH1F*>(fOutput->FindObject("fPtPion "));
784 fDStar = dynamic_cast<TH1F*>(fOutput->FindObject("fDStar"));
785 fDiff = dynamic_cast<TH1F*>(fOutput->FindObject("fDiff"));
786 fDiffSideBand = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand"));
787 ftrigger = dynamic_cast<TH1F*>(fOutput->FindObject("ftrigger"));
788 fRECOPtDStar = dynamic_cast<TH1F*>(fOutput->FindObject("fRECOPtDStar"));
789 fEjet = dynamic_cast<TH1F*>(fOutput->FindObject("fEjet"));
790 fPhijet = dynamic_cast<TH1F*>(fOutput->FindObject("fPhijet"));
791 fEtaJet = dynamic_cast<TH1F*>(fOutput->FindObject("fEtaJet"));
792 fPhi = dynamic_cast<TH1F*>(fOutput->FindObject("fPhi"));
793 fResZ = dynamic_cast<TH1F*>(fOutput->FindObject("fResZ"));
794 fResZBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fResZBkg"));
795 fPhiBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fPhiBkg"));
796 fD0ptvsSoftPtSignal = dynamic_cast<TH2F*>(fOutput->FindObject("fD0ptvsSoftPtSignal"));
797 fD0ptvsSoftPt = dynamic_cast<TH2F*>(fOutput->FindObject("fD0ptvsSoftPt"));
798
799}
800//___________________________________________________________________________
801
802void AliAnalysisTaskSEDStarJets::UserCreateOutputObjects() {
803 // output
804
805 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
806
807 //slot #1
808 OpenFile(1);
809 fOutput = new TList();
810 fOutput->SetOwner();
811 // define histograms
812 DefineHistoFroAnalysis();
813
814 return;
815}
954ac830 816
c088a94c 817//_______________________________D0 invariant mass________________________________
954ac830 818
c088a94c 819Double_t AliAnalysisTaskSEDStarJets::GetInvariantMass(TLorentzVector LorentzTrack1, TLorentzVector LorentzTrack2){
820
821 return TMath::Abs((LorentzTrack1+LorentzTrack2).M()); // invariant mass of two tracks
822}
823//______________________________D* invariant mass_________________________________
954ac830 824
c088a94c 825Double_t AliAnalysisTaskSEDStarJets::GetInvariantMassDStar(TLorentzVector LorentzTrack4, TLorentzVector LorentzTrack3){
826
827 return TMath::Abs((LorentzTrack4+LorentzTrack3).M()); // invariant mass of two tracks
828}
829//_________________________________D* in MC _______________________________________
954ac830 830
c088a94c 831Bool_t AliAnalysisTaskSEDStarJets::DstarInMC(AliAODMCParticle* const mcPart, TClonesArray* mcArray){
832 // D* in MC
833
834 Bool_t isOk = kFALSE;
835
836 if(TMath::Abs(mcPart->GetPdgCode())!=413) return isOk;
954ac830 837
838 // getting the daughters
839 Int_t daughter0 = mcPart->GetDaughter(0);
840 Int_t daughter1 = mcPart->GetDaughter(1);
841
842 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
843 if (daughter0 == 0 || daughter1 == 0) {
c088a94c 844 AliDebug(2, "Error! the D* MC doesn't have correct daughters!!");
954ac830 845 return isOk;
846 }
c088a94c 847
848 if (TMath::Abs(daughter1 - daughter0) != 1) { // should be everytime true - see PDGdatabooklet
849 AliDebug(2, "The D* MC doesn't come from a 2-prong decay, skipping!!");
954ac830 850 return isOk;
851 }
954ac830 852
c088a94c 853 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
854 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
954ac830 855 if (!mcPartDaughter0 || !mcPartDaughter1) {
856 AliWarning("At least one Daughter Particle not found in tree, skipping");
857 return isOk;
858 }
859
c088a94c 860 if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==421 &&
861 TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
862 !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
863 TMath::Abs(mcPartDaughter1->GetPdgCode())==421)) {
864 AliDebug(2, "The D* MC doesn't come from a Kpi decay, skipping!!");
865 return isOk;
866 }
867
868 Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0
869 Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1
870
871 // getting vertex from daughters
872 mcPartDaughter0->XvYvZv(vtx2daughter0); // cm
873 mcPartDaughter1->XvYvZv(vtx2daughter1); //cm
874
875 // check if the secondary vertex is the same for both
876 if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
877 AliError("Daughters have different secondary vertex, skipping the track");
878 return isOk;
879 }
880
881 AliAODMCParticle* neutralDaugh = mcPartDaughter0;
882
883 if (!EvaluateIfD0toKpi(neutralDaugh,mcArray)) {
884 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
885 return isOk;
886 }
954ac830 887
888 isOk = kTRUE;
889
890 return isOk;
954ac830 891
954ac830 892}
954ac830 893
c088a94c 894Bool_t AliAnalysisTaskSEDStarJets::EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh, TClonesArray* mcArray)const{
954ac830 895
c088a94c 896 //
897 // chack wether D0 is decaing into kpi
898 //
954ac830 899
c088a94c 900 Bool_t isHadronic = kFALSE;
954ac830 901
c088a94c 902 Int_t daughterD00 = neutralDaugh->GetDaughter(0);
903 Int_t daughterD01 = neutralDaugh->GetDaughter(1);
954ac830 904
c088a94c 905 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughterD00,daughterD01));
906 if (daughterD00 == 0 || daughterD01 == 0) {
907 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
908 return isHadronic;
954ac830 909 }
c088a94c 910
911 if (TMath::Abs(daughterD01 - daughterD00) != 1) { // should be everytime true - see PDGdatabooklet
912 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
913 return isHadronic;
954ac830 914 }
c088a94c 915
916 AliAODMCParticle* mcPartDaughterD00 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD00));
917 AliAODMCParticle* mcPartDaughterD01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD01));
918 if (!mcPartDaughterD00 || !mcPartDaughterD01) {
954ac830 919 AliWarning("At least one Daughter Particle not found in tree, skipping");
c088a94c 920 return isHadronic;
921 }
954ac830 922
c088a94c 923 if (!(TMath::Abs(mcPartDaughterD00->GetPdgCode())==321 &&
924 TMath::Abs(mcPartDaughterD01->GetPdgCode())==211) &&
925 !(TMath::Abs(mcPartDaughterD00->GetPdgCode())==211 &&
926 TMath::Abs(mcPartDaughterD01->GetPdgCode())==321)) {
927 AliDebug(2, "The D* MC doesn't come from a Kpi decay, skipping!!");
928 return isHadronic;
929 }
954ac830 930
c088a94c 931 isHadronic = kTRUE;
932
933
934 return isHadronic;
954ac830 935
936}
937
c088a94c 938
954ac830 939//___________________________________ hiostograms _______________________________________
940
941Bool_t AliAnalysisTaskSEDStarJets::DefineHistoFroAnalysis(){
942
943 // Invariant mass related histograms
944 fInvMass = new TH1F("invMass","Kpi invariant mass distribution",1500,.5,3.5);
945 fInvMass->SetStats(kTRUE);
946 fInvMass->GetXaxis()->SetTitle("GeV/c");
947 fInvMass->GetYaxis()->SetTitle("Entries");
948 fOutput->Add(fInvMass);
949
950 fDStar = new TH1F("invMassDStar","DStar invariant mass after D0 cuts ",600,1.8,2.4);
951 fDStar->SetStats(kTRUE);
952 fDStar->GetXaxis()->SetTitle("GeV/c");
953 fDStar->GetYaxis()->SetTitle("Entries");
954 fOutput->Add(fDStar);
955
956 fDiff = new TH1F("Diff","M(kpipi)-M(kpi)",750,0.1,0.2);
957 fDiff->SetStats(kTRUE);
958 fDiff->GetXaxis()->SetTitle("M(kpipi)-M(kpi) GeV/c^2");
959 fDiff->GetYaxis()->SetTitle("Entries");
960 fOutput->Add(fDiff);
961
962 fDiffSideBand = new TH1F("DiffSide","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
963 fDiffSideBand->SetStats(kTRUE);
964 fDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
965 fDiffSideBand->GetYaxis()->SetTitle("Entries");
966 fOutput->Add(fDiffSideBand);
967
968 fDStarMass = new TH1F("RECODStar2","RECO DStar invariant mass distribution",750,1.5,2.5);
969 fOutput->Add(fDStarMass);
970
971 fTrueDiff = new TH1F("dstar","True Reco diff",750,0,0.2);
972 fOutput->Add(fTrueDiff);
973
974 // trigger normalization
975 ftrigger = new TH1F("Normalization","Normalization factor for correlations",1,0,10);
976 ftrigger->SetStats(kTRUE);
977 fOutput->Add(ftrigger);
978
979 //correlation fistograms
980 fPhi = new TH1F("phi","Delta phi between Jet axis and DStar ",25,-1.57,4.72);
981 fPhi->SetStats(kTRUE);
982 fPhi->GetXaxis()->SetTitle("#Delta #phi (rad)");
983 fPhi->GetYaxis()->SetTitle("Entries");
984 fOutput->Add(fPhi);
985
986 fPhiBkg = new TH1F("phiBkg","Delta phi between Jet axis and DStar background ",25,-1.57,4.72);
987 fPhiBkg->SetStats(kTRUE);
988 fPhiBkg->GetXaxis()->SetTitle("#Delta #phi (rad)");
989 fPhiBkg->GetYaxis()->SetTitle("Entries");
990 fOutput->Add(fPhiBkg);
991
992 fRECOPtDStar = new TH1F("RECODStar1","RECO DStar pt distribution",600,0,15);
993 fRECOPtDStar->SetStats(kTRUE);
994 fRECOPtDStar->SetLineColor(2);
995 fRECOPtDStar->GetXaxis()->SetTitle("GeV/c");
996 fRECOPtDStar->GetYaxis()->SetTitle("Entries");
997 fOutput->Add(fRECOPtDStar);
998
999 fPtPion = new TH1F("pionpt","Primary pions candidates pt ",500,0,5);
1000 fPtPion->SetStats(kTRUE);
1001 fPtPion->GetXaxis()->SetTitle("GeV/c");
1002 fPtPion->GetYaxis()->SetTitle("Entries");
1003 fOutput->Add(fPtPion);
1004
1005 fcharmpt = new TH1F("charmpt","MC Charm pt distribution", 10000,0,5000);
1006 fdstarE = new TH1F("dstare", "MC D* energy distribution", 10000,0,5000);
1007 fdstarpt = new TH1F("dstarpt","MC D* pt distribution", 10000,0,5000);
1008 fOutput->Add(fcharmpt);
1009 fOutput->Add(fdstarE);
1010 fOutput->Add(fdstarpt);
1011
1012 // jet related fistograms
1013 fEjet = new TH1F("ejet", "UA1 algorithm jet energy distribution",1000,0,500);
1014 fPhijet = new TH1F("Phijet","UA1 algorithm jet #phi distribution", 200,-7,7);
1015 fEtaJet = new TH1F("Etajet","UA1 algorithm jet #eta distribution", 200,-7,7);
1016 fOutput->Add(fEjet);
1017 fOutput->Add(fPhijet);
1018 fOutput->Add(fEtaJet);
954ac830 1019
1020 fResZ = new TH1F("FragFunc","Fragmentation function ",50,0,1);
1021 fResZBkg = new TH1F("FragFuncBkg","Fragmentation function background",50,0,1);
1022 fOutput->Add(fResZ);
1023 fOutput->Add(fResZBkg);
1024
1025 fD0ptvsSoftPt = new TH2F("D0piRec","Candidates (background + signal)",100,0,3,100,0,15);
1026 fD0ptvsSoftPtSignal = new TH2F("D0PiSignal","Signal",100,0,3,100,0,15);
1027 fOutput->Add(fD0ptvsSoftPt);
1028 fOutput->Add(fD0ptvsSoftPtSignal);
1029
1030 return kTRUE;
1031
1032}
1033
1034//______________________________Phase space cut alternative to PPR _______________________________________
1035
1036Bool_t AliAnalysisTaskSEDStarJets::EvaluateCutOnPiD0pt(AliAODRecoDecayHF2Prong* const vtx, AliAODTrack* const aodTrack) {
1037
1038 // The soft pion pt and DO pt are strongly correlated. It can be shown that ~ 95% of the signal in constrained
1039 // into a narrow band defined by 10 < D0pt/softPt < 18. This cut can be used with a relaxed CosThetaStar cut
1040 // to reconstruct the D*.
1041
1042 Double_t softPt = 0;
1043 Double_t d0ptReco = 0;
1044
1045 softPt = aodTrack->Pt();
1046 d0ptReco = vtx->Pt();
1047 fD0ptvsSoftPt->Fill(softPt,d0ptReco);
1048
1049 if(softPt>0){
1050 Double_t ratio = d0ptReco/softPt;
1051 if( ratio <=10. || ratio>=18. ) return kFALSE;
1052 }
1053
1054 fD0ptvsSoftPtSignal->Fill(softPt,d0ptReco);
1055 return kTRUE;
1056}
1057
1058//______________________________ side band background for D*___________________________________
1059
feaae220 1060void AliAnalysisTaskSEDStarJets::SideBandBackground(Double_t invM, Double_t invMDStar, Double_t ejet, Double_t dPhi, Int_t nJets){
954ac830 1061
1062 // D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
1063 // (expected detector resolution) on the left and right frm the D0 mass. Each band
1064 // has a width of ~5 sigmas. Two band needed for opening angle considerations
1065
1066 if((invM>=1.763 && invM<=1.811) || (invM>=1.919 && invM<=1.963)){
1067
feaae220 1068 if(nJets==1)fDiffSideBand->Fill(invMDStar-invM); // M(Kpipi)-M(Kpi) side band background
954ac830 1069
feaae220 1070 if ((invMDStar-invM)<=0.148 && (invMDStar-invM)>=0.142) {
954ac830 1071 fPhiBkg->Fill(dPhi);
1072
1073 if(dPhi>=-0.5 && dPhi<=0.5){ // evaluate in the near side
1074
1075 Double_t dStarMomBkg = (fLorentzTrack3 + fLorentzTrack4).P();
1076 Double_t zFragBkg = (TMath::Cos(dPhi)*dStarMomBkg)/ejet;
1077 fResZBkg->Fill(TMath::Abs(zFragBkg));
1078
1079 }
1080 }
1081 }
1082}