]> 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
fixed bugs with the rule checker.
[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++;
06f6542f 238 AliDebug(2,Form("Event %d",fEvents));
239 if (fEvents%10000 ==0) AliDebug(2,Form("Event %d",fEvents));
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
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
954ac830 761 AliAnalysisTaskSE::Terminate();
762
954ac830 763 AliInfo(Form("Found %i of that MC D*->D0pi(D0->kpi) are in acceptance, in %d events",fCountDStar,fEvents));
764
765 AliInfo(Form("Found %i RECO particles after cuts that are DStar, in %d events",fCountReco,fEvents));
c088a94c 766 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));
767 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));
768 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 769
770 fOutput = dynamic_cast<TList*> (GetOutputData(1));
771 if (!fOutput) {
772 printf("ERROR: fOutput not available\n");
773 return;
774 }
775
954ac830 776 fcharmpt = dynamic_cast<TH1F*>(fOutput->FindObject("fcharmpt"));
777 fdstarE = dynamic_cast<TH1F*>(fOutput->FindObject("fdstarE"));
778 fdstarpt = dynamic_cast<TH1F*>(fOutput->FindObject("fdstarpt"));
779 fDStarMass = dynamic_cast<TH1F*>(fOutput->FindObject("fDStarMass"));
780 fTrueDiff = dynamic_cast<TH1F*>(fOutput->FindObject("fTrueDiff"));
781 fInvMass = dynamic_cast<TH1F*>(fOutput->FindObject("fInvMass"));
782 fPtPion = dynamic_cast<TH1F*>(fOutput->FindObject("fPtPion "));
783 fDStar = dynamic_cast<TH1F*>(fOutput->FindObject("fDStar"));
784 fDiff = dynamic_cast<TH1F*>(fOutput->FindObject("fDiff"));
785 fDiffSideBand = dynamic_cast<TH1F*>(fOutput->FindObject("fDiffSideBand"));
786 ftrigger = dynamic_cast<TH1F*>(fOutput->FindObject("ftrigger"));
787 fRECOPtDStar = dynamic_cast<TH1F*>(fOutput->FindObject("fRECOPtDStar"));
788 fEjet = dynamic_cast<TH1F*>(fOutput->FindObject("fEjet"));
789 fPhijet = dynamic_cast<TH1F*>(fOutput->FindObject("fPhijet"));
790 fEtaJet = dynamic_cast<TH1F*>(fOutput->FindObject("fEtaJet"));
791 fPhi = dynamic_cast<TH1F*>(fOutput->FindObject("fPhi"));
792 fResZ = dynamic_cast<TH1F*>(fOutput->FindObject("fResZ"));
793 fResZBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fResZBkg"));
794 fPhiBkg = dynamic_cast<TH1F*>(fOutput->FindObject("fPhiBkg"));
795 fD0ptvsSoftPtSignal = dynamic_cast<TH2F*>(fOutput->FindObject("fD0ptvsSoftPtSignal"));
796 fD0ptvsSoftPt = dynamic_cast<TH2F*>(fOutput->FindObject("fD0ptvsSoftPt"));
797
798}
799//___________________________________________________________________________
800
801void AliAnalysisTaskSEDStarJets::UserCreateOutputObjects() {
802 // output
803
804 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
805
806 //slot #1
807 OpenFile(1);
808 fOutput = new TList();
809 fOutput->SetOwner();
810 // define histograms
811 DefineHistoFroAnalysis();
812
813 return;
814}
954ac830 815
c088a94c 816//_______________________________D0 invariant mass________________________________
954ac830 817
c088a94c 818Double_t AliAnalysisTaskSEDStarJets::GetInvariantMass(TLorentzVector LorentzTrack1, TLorentzVector LorentzTrack2){
819
820 return TMath::Abs((LorentzTrack1+LorentzTrack2).M()); // invariant mass of two tracks
821}
822//______________________________D* invariant mass_________________________________
954ac830 823
c088a94c 824Double_t AliAnalysisTaskSEDStarJets::GetInvariantMassDStar(TLorentzVector LorentzTrack4, TLorentzVector LorentzTrack3){
825
826 return TMath::Abs((LorentzTrack4+LorentzTrack3).M()); // invariant mass of two tracks
827}
828//_________________________________D* in MC _______________________________________
954ac830 829
c088a94c 830Bool_t AliAnalysisTaskSEDStarJets::DstarInMC(AliAODMCParticle* const mcPart, TClonesArray* mcArray){
831 // D* in MC
832
833 Bool_t isOk = kFALSE;
834
835 if(TMath::Abs(mcPart->GetPdgCode())!=413) return isOk;
954ac830 836
837 // getting the daughters
838 Int_t daughter0 = mcPart->GetDaughter(0);
839 Int_t daughter1 = mcPart->GetDaughter(1);
840
841 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
842 if (daughter0 == 0 || daughter1 == 0) {
c088a94c 843 AliDebug(2, "Error! the D* MC doesn't have correct daughters!!");
954ac830 844 return isOk;
845 }
c088a94c 846
847 if (TMath::Abs(daughter1 - daughter0) != 1) { // should be everytime true - see PDGdatabooklet
848 AliDebug(2, "The D* MC doesn't come from a 2-prong decay, skipping!!");
954ac830 849 return isOk;
850 }
954ac830 851
c088a94c 852 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
853 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
954ac830 854 if (!mcPartDaughter0 || !mcPartDaughter1) {
855 AliWarning("At least one Daughter Particle not found in tree, skipping");
856 return isOk;
857 }
858
c088a94c 859 if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==421 &&
860 TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
861 !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
862 TMath::Abs(mcPartDaughter1->GetPdgCode())==421)) {
863 AliDebug(2, "The D* MC doesn't come from a Kpi decay, skipping!!");
864 return isOk;
865 }
866
867 Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0
868 Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1
869
870 // getting vertex from daughters
871 mcPartDaughter0->XvYvZv(vtx2daughter0); // cm
872 mcPartDaughter1->XvYvZv(vtx2daughter1); //cm
873
874 // check if the secondary vertex is the same for both
875 if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
876 AliError("Daughters have different secondary vertex, skipping the track");
877 return isOk;
878 }
879
880 AliAODMCParticle* neutralDaugh = mcPartDaughter0;
881
882 if (!EvaluateIfD0toKpi(neutralDaugh,mcArray)) {
883 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
884 return isOk;
885 }
954ac830 886
887 isOk = kTRUE;
888
889 return isOk;
954ac830 890
954ac830 891}
954ac830 892
c088a94c 893Bool_t AliAnalysisTaskSEDStarJets::EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh, TClonesArray* mcArray)const{
954ac830 894
c088a94c 895 //
896 // chack wether D0 is decaing into kpi
897 //
954ac830 898
c088a94c 899 Bool_t isHadronic = kFALSE;
954ac830 900
c088a94c 901 Int_t daughterD00 = neutralDaugh->GetDaughter(0);
902 Int_t daughterD01 = neutralDaugh->GetDaughter(1);
954ac830 903
c088a94c 904 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughterD00,daughterD01));
905 if (daughterD00 == 0 || daughterD01 == 0) {
906 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
907 return isHadronic;
954ac830 908 }
c088a94c 909
910 if (TMath::Abs(daughterD01 - daughterD00) != 1) { // should be everytime true - see PDGdatabooklet
911 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
912 return isHadronic;
954ac830 913 }
c088a94c 914
915 AliAODMCParticle* mcPartDaughterD00 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD00));
916 AliAODMCParticle* mcPartDaughterD01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD01));
917 if (!mcPartDaughterD00 || !mcPartDaughterD01) {
954ac830 918 AliWarning("At least one Daughter Particle not found in tree, skipping");
c088a94c 919 return isHadronic;
920 }
954ac830 921
c088a94c 922 if (!(TMath::Abs(mcPartDaughterD00->GetPdgCode())==321 &&
923 TMath::Abs(mcPartDaughterD01->GetPdgCode())==211) &&
924 !(TMath::Abs(mcPartDaughterD00->GetPdgCode())==211 &&
925 TMath::Abs(mcPartDaughterD01->GetPdgCode())==321)) {
926 AliDebug(2, "The D* MC doesn't come from a Kpi decay, skipping!!");
927 return isHadronic;
928 }
954ac830 929
c088a94c 930 isHadronic = kTRUE;
931
932
933 return isHadronic;
954ac830 934
935}
936
c088a94c 937
954ac830 938//___________________________________ hiostograms _______________________________________
939
940Bool_t AliAnalysisTaskSEDStarJets::DefineHistoFroAnalysis(){
941
942 // Invariant mass related histograms
943 fInvMass = new TH1F("invMass","Kpi invariant mass distribution",1500,.5,3.5);
944 fInvMass->SetStats(kTRUE);
945 fInvMass->GetXaxis()->SetTitle("GeV/c");
946 fInvMass->GetYaxis()->SetTitle("Entries");
947 fOutput->Add(fInvMass);
948
949 fDStar = new TH1F("invMassDStar","DStar invariant mass after D0 cuts ",600,1.8,2.4);
950 fDStar->SetStats(kTRUE);
951 fDStar->GetXaxis()->SetTitle("GeV/c");
952 fDStar->GetYaxis()->SetTitle("Entries");
953 fOutput->Add(fDStar);
954
955 fDiff = new TH1F("Diff","M(kpipi)-M(kpi)",750,0.1,0.2);
956 fDiff->SetStats(kTRUE);
957 fDiff->GetXaxis()->SetTitle("M(kpipi)-M(kpi) GeV/c^2");
958 fDiff->GetYaxis()->SetTitle("Entries");
959 fOutput->Add(fDiff);
960
961 fDiffSideBand = new TH1F("DiffSide","M(kpipi)-M(kpi) Side Band Background",750,0.1,0.2);
962 fDiffSideBand->SetStats(kTRUE);
963 fDiffSideBand->GetXaxis()->SetTitle("M(kpipi)-M(Kpi) GeV/c^2");
964 fDiffSideBand->GetYaxis()->SetTitle("Entries");
965 fOutput->Add(fDiffSideBand);
966
967 fDStarMass = new TH1F("RECODStar2","RECO DStar invariant mass distribution",750,1.5,2.5);
968 fOutput->Add(fDStarMass);
969
970 fTrueDiff = new TH1F("dstar","True Reco diff",750,0,0.2);
971 fOutput->Add(fTrueDiff);
972
973 // trigger normalization
974 ftrigger = new TH1F("Normalization","Normalization factor for correlations",1,0,10);
975 ftrigger->SetStats(kTRUE);
976 fOutput->Add(ftrigger);
977
978 //correlation fistograms
979 fPhi = new TH1F("phi","Delta phi between Jet axis and DStar ",25,-1.57,4.72);
980 fPhi->SetStats(kTRUE);
981 fPhi->GetXaxis()->SetTitle("#Delta #phi (rad)");
982 fPhi->GetYaxis()->SetTitle("Entries");
983 fOutput->Add(fPhi);
984
985 fPhiBkg = new TH1F("phiBkg","Delta phi between Jet axis and DStar background ",25,-1.57,4.72);
986 fPhiBkg->SetStats(kTRUE);
987 fPhiBkg->GetXaxis()->SetTitle("#Delta #phi (rad)");
988 fPhiBkg->GetYaxis()->SetTitle("Entries");
989 fOutput->Add(fPhiBkg);
990
991 fRECOPtDStar = new TH1F("RECODStar1","RECO DStar pt distribution",600,0,15);
992 fRECOPtDStar->SetStats(kTRUE);
993 fRECOPtDStar->SetLineColor(2);
994 fRECOPtDStar->GetXaxis()->SetTitle("GeV/c");
995 fRECOPtDStar->GetYaxis()->SetTitle("Entries");
996 fOutput->Add(fRECOPtDStar);
997
998 fPtPion = new TH1F("pionpt","Primary pions candidates pt ",500,0,5);
999 fPtPion->SetStats(kTRUE);
1000 fPtPion->GetXaxis()->SetTitle("GeV/c");
1001 fPtPion->GetYaxis()->SetTitle("Entries");
1002 fOutput->Add(fPtPion);
1003
1004 fcharmpt = new TH1F("charmpt","MC Charm pt distribution", 10000,0,5000);
1005 fdstarE = new TH1F("dstare", "MC D* energy distribution", 10000,0,5000);
1006 fdstarpt = new TH1F("dstarpt","MC D* pt distribution", 10000,0,5000);
1007 fOutput->Add(fcharmpt);
1008 fOutput->Add(fdstarE);
1009 fOutput->Add(fdstarpt);
1010
1011 // jet related fistograms
1012 fEjet = new TH1F("ejet", "UA1 algorithm jet energy distribution",1000,0,500);
1013 fPhijet = new TH1F("Phijet","UA1 algorithm jet #phi distribution", 200,-7,7);
1014 fEtaJet = new TH1F("Etajet","UA1 algorithm jet #eta distribution", 200,-7,7);
1015 fOutput->Add(fEjet);
1016 fOutput->Add(fPhijet);
1017 fOutput->Add(fEtaJet);
954ac830 1018
1019 fResZ = new TH1F("FragFunc","Fragmentation function ",50,0,1);
1020 fResZBkg = new TH1F("FragFuncBkg","Fragmentation function background",50,0,1);
1021 fOutput->Add(fResZ);
1022 fOutput->Add(fResZBkg);
1023
1024 fD0ptvsSoftPt = new TH2F("D0piRec","Candidates (background + signal)",100,0,3,100,0,15);
1025 fD0ptvsSoftPtSignal = new TH2F("D0PiSignal","Signal",100,0,3,100,0,15);
1026 fOutput->Add(fD0ptvsSoftPt);
1027 fOutput->Add(fD0ptvsSoftPtSignal);
1028
1029 return kTRUE;
1030
1031}
1032
1033//______________________________Phase space cut alternative to PPR _______________________________________
1034
1035Bool_t AliAnalysisTaskSEDStarJets::EvaluateCutOnPiD0pt(AliAODRecoDecayHF2Prong* const vtx, AliAODTrack* const aodTrack) {
1036
1037 // The soft pion pt and DO pt are strongly correlated. It can be shown that ~ 95% of the signal in constrained
1038 // into a narrow band defined by 10 < D0pt/softPt < 18. This cut can be used with a relaxed CosThetaStar cut
1039 // to reconstruct the D*.
1040
1041 Double_t softPt = 0;
1042 Double_t d0ptReco = 0;
1043
1044 softPt = aodTrack->Pt();
1045 d0ptReco = vtx->Pt();
1046 fD0ptvsSoftPt->Fill(softPt,d0ptReco);
1047
1048 if(softPt>0){
1049 Double_t ratio = d0ptReco/softPt;
1050 if( ratio <=10. || ratio>=18. ) return kFALSE;
1051 }
1052
1053 fD0ptvsSoftPtSignal->Fill(softPt,d0ptReco);
1054 return kTRUE;
1055}
1056
1057//______________________________ side band background for D*___________________________________
1058
feaae220 1059void AliAnalysisTaskSEDStarJets::SideBandBackground(Double_t invM, Double_t invMDStar, Double_t ejet, Double_t dPhi, Int_t nJets){
954ac830 1060
1061 // D* side band background method. Two side bands, in M(Kpi) are taken at ~6 sigmas
1062 // (expected detector resolution) on the left and right frm the D0 mass. Each band
1063 // has a width of ~5 sigmas. Two band needed for opening angle considerations
1064
1065 if((invM>=1.763 && invM<=1.811) || (invM>=1.919 && invM<=1.963)){
1066
feaae220 1067 if(nJets==1)fDiffSideBand->Fill(invMDStar-invM); // M(Kpipi)-M(Kpi) side band background
954ac830 1068
feaae220 1069 if ((invMDStar-invM)<=0.148 && (invMDStar-invM)>=0.142) {
954ac830 1070 fPhiBkg->Fill(dPhi);
1071
1072 if(dPhi>=-0.5 && dPhi<=0.5){ // evaluate in the near side
1073
1074 Double_t dStarMomBkg = (fLorentzTrack3 + fLorentzTrack4).P();
1075 Double_t zFragBkg = (TMath::Cos(dPhi)*dStarMomBkg)/ejet;
1076 fResZBkg->Fill(TMath::Abs(zFragBkg));
1077
1078 }
1079 }
1080 }
1081}