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