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