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