]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/AliCFTaskForDStarAnalysis.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Added new class for 3prong decays within new CF framework (Francesco)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliCFTaskForDStarAnalysis.cxx
CommitLineData
2854fd06 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// Class for DStar corrections:
18//
19// The D0 cutting varibles position in the container and container
20// binning method from a C.Zampolli example
21// In this way a simple comparison between D0 and D0 from D* is possible.
22//
23//-----------------------------------------------------------------------
24// Author : A.Grelli, Utrecht University
25//
26// a.grelli@uu.nl
27//-----------------------------------------------------------------------
28
29#include <TCanvas.h>
30#include <TParticle.h>
31#include <TDatabasePDG.h>
32#include <TH1I.h>
33#include <TStyle.h>
34#include <TFile.h>
35
36#include "AliCFTaskForDStarAnalysis.h"
37#include "AliStack.h"
38#include "AliMCEvent.h"
39#include "AliCFManager.h"
40#include "AliCFContainer.h"
41#include "AliLog.h"
42#include "AliAnalysisManager.h"
43#include "AliAODHandler.h"
44#include "AliAODEvent.h"
45#include "AliAODRecoDecay.h"
46#include "AliAODRecoDecayHF.h"
47#include "AliAODRecoCascadeHF.h"
48#include "AliAODRecoDecayHF2Prong.h"
49#include "AliAODMCParticle.h"
50#include "AliAODMCHeader.h"
51#include "AliESDtrack.h"
52#include "TChain.h"
53#include "THnSparse.h"
54#include "TH2D.h"
55
56//__________________________________________________________________________
57AliCFTaskForDStarAnalysis::AliCFTaskForDStarAnalysis() :
58 AliAnalysisTaskSE(),
59 fCFManager(0x0),
60 fHistEventsProcessed(0x0),
61 fCorrelation(0x0),
62 fCountRecoDStarSel(0),
63 fEvents(0),
64 fMinITSClusters(5),
65 fMinITSClustersSoft(4),
66 fAcceptanceUnf(kTRUE)
67{
68 //
69 //Default ctor
70 //
71}
72//___________________________________________________________________________
73AliCFTaskForDStarAnalysis::AliCFTaskForDStarAnalysis(const Char_t* name) :
74 AliAnalysisTaskSE(name),
75 fCFManager(0x0),
76 fHistEventsProcessed(0x0),
77 fCorrelation(0x0),
78 fCountRecoDStarSel(0),
79 fEvents(0),
80 fMinITSClusters(5),
81 fMinITSClustersSoft(4),
82 fAcceptanceUnf(kTRUE)
83{
84 //
85 // Constructor. Initialization of Inputs and Outputs
86 //
87 Info("AliCFTaskForDStarAnalysis","Calling Constructor");
88
89 DefineOutput(1,TH1I::Class());
90 DefineOutput(2,AliCFContainer::Class());
91 DefineOutput(3,THnSparseD::Class());
92}
93
94//___________________________________________________________________________
95AliCFTaskForDStarAnalysis& AliCFTaskForDStarAnalysis::operator=(const AliCFTaskForDStarAnalysis& c)
96{
97 //
98 // Assignment operator
99 //
100 if (this!=&c) {
101 AliAnalysisTaskSE::operator=(c) ;
102 fCFManager = c.fCFManager;
103 fHistEventsProcessed = c.fHistEventsProcessed;
104 }
105 return *this;
106}
107
108//___________________________________________________________________________
109AliCFTaskForDStarAnalysis::AliCFTaskForDStarAnalysis(const AliCFTaskForDStarAnalysis& c) :
110 AliAnalysisTaskSE(c),
111 fCFManager(c.fCFManager),
112 fHistEventsProcessed(c.fHistEventsProcessed),
113 fCorrelation(c.fCorrelation),
114 fCountRecoDStarSel(c.fCountRecoDStarSel),
115 fEvents(c.fEvents),
116 fMinITSClusters(c.fMinITSClusters),
117 fMinITSClustersSoft(c.fMinITSClustersSoft),
118 fAcceptanceUnf(c.fAcceptanceUnf)
119{
120 //
121 // Copy Constructor
122 //
123}
124
125//___________________________________________________________________________
126AliCFTaskForDStarAnalysis::~AliCFTaskForDStarAnalysis() {
127 //
128 //destructor
129 //
130 if (fCFManager) delete fCFManager ;
131 if (fHistEventsProcessed) delete fHistEventsProcessed ;
132 if (fCorrelation) delete fCorrelation ;
133}
134
135//_________________________________________________
136void AliCFTaskForDStarAnalysis::UserExec(Option_t *)
137{
138 //
139 // Main loop function
140 //
141
142 if (!fInputEvent) {
143 Error("UserExec","NO EVENT FOUND!");
144 return;
145 }
146
2854fd06 147 AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
148
149 TClonesArray *arrayDStartoD0pi=0;
150
151 if(!aodEvent && AODEvent() && IsStandardAOD()) {
152 // In case there is an AOD handler writing a standard AOD, use the AOD
153 // event in memory rather than the input (ESD) event.
154 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
155 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
156 // have to taken from the AOD event hold by the AliAODExtension
157 AliAODHandler* aodHandler = (AliAODHandler*)
158 ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
159 if(aodHandler->GetExtensions()) {
160 AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
161 AliAODEvent *aodFromExt = ext->GetAOD();
162 arrayDStartoD0pi=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
163 }
164 } else {
165 arrayDStartoD0pi=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
166 }
167
168 if (!arrayDStartoD0pi) {
169 AliError("Could not find array of HF vertices");
170 return;
171 }
172
7c23877d 173 // fix for temporary bug in ESDfilter
174 // the AODs with null vertex pointer didn't pass the PhysSel
175 if(!aodEvent->GetPrimaryVertex()) return;
176
177 fEvents++;
178 if (fEvents%10000 ==0) AliDebug(2,Form("Event %d",fEvents));
179
2854fd06 180 fCFManager->SetRecEventInfo(aodEvent);
181 fCFManager->SetMCEventInfo(aodEvent);
182
183 // event selection
06aeb2ef 184 Double_t containerInput[15] ;
185 Double_t containerInputMC[15] ;
2854fd06 186
187 //loop on the MC event
188
189 TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
190 if (!mcArray) {
191 AliError("Could not find Monte-Carlo in AOD");
192 return;
193 }
194
195 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
196 if (!mcHeader) {
197 AliError("Could not find MC Header in AOD");
198 return;
199 }
200
201 // AOD primary vertex
202 AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
203 Double_t zPrimVertex = vtx1->GetZ();
204 Double_t zMCVertex = mcHeader->GetVtxZ();
205 Bool_t vtxFlag = kTRUE;
206 TString title=vtx1->GetTitle();
207 if(!title.Contains("VertexerTracks")) vtxFlag=kFALSE;
208
209 for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
210 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
211 if (!mcPart) {
212 AliWarning("MC Particle not found in tree, skipping");
213 continue;
214 }
215
216 // check the MC-level cuts
217 if (!fCFManager->CheckParticleCuts(0, mcPart)) continue; // 0 stands for MC level
218
219 // fill the container for Gen-level selection
220 Double_t vectorMC[9] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.,9999.,9999.};
221
222 if (GetDStarMCParticle(mcPart, mcArray, vectorMC)){
223
224 containerInputMC[0] = vectorMC[0];
225 containerInputMC[1] = vectorMC[1] ;
226 containerInputMC[2] = vectorMC[2] ;
227 containerInputMC[3] = vectorMC[3] ;
228 containerInputMC[4] = vectorMC[4] ;
229 containerInputMC[5] = vectorMC[5] ;
230 containerInputMC[6] = 0.;
231 containerInputMC[7] = 0.;
232 containerInputMC[8] = 0.;
233 containerInputMC[9] = -100000.;
234 containerInputMC[10] = 1.01;
235 containerInputMC[11] = vectorMC[6];
236 containerInputMC[12] = zMCVertex; // z vertex
237 containerInputMC[13] = vectorMC[7]; // pt D0 pion
238 containerInputMC[14] = vectorMC[8]; // pt D0 kaon
239
240 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGenerated);
241
242 // check the MC-Acceptance level cuts
243 // since standard CF functions are not applicable, using Kine Cuts on daughters
244
245 Int_t daughter0 = mcPart->GetDaughter(0);
246 Int_t daughter1 = mcPart->GetDaughter(1);
247
248 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
249
250 //D0
251 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
252 // Soft Pion
253 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
254
255 // Acceptance variables for the soft pion
256 Double_t eta1 = mcPartDaughter1->Eta();
257 Double_t pt1 = mcPartDaughter1->Pt();
258
259 Int_t daughD00 = 0;
260 Int_t daughD01 = 0;
261
262 // Just to be sure to take the right particles
263 if(TMath::Abs(mcPartDaughter0->GetPdgCode())==421){
264 daughD00 = mcPartDaughter0->GetDaughter(0);
265 daughD01 = mcPartDaughter0->GetDaughter(1);
266 }else{
267 daughD00 = mcPartDaughter1->GetDaughter(0);
268 daughD01 = mcPartDaughter1->GetDaughter(1);
269 }
270
271 // the D0 daughters
272 AliAODMCParticle* mcD0PartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughD00));
273 AliAODMCParticle* mcD0PartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughD01));
274
275 if (!mcD0PartDaughter0 || !mcD0PartDaughter1) {
276 AliWarning("At least one D0 Daughter Particle not found in tree, but it should be, this check was already done...");
277 }
278
279 // D0 daughters - needed for acceptance
280 Double_t theD0pt0 = mcD0PartDaughter0->Pt();
281 Double_t theD0pt1 = mcD0PartDaughter1->Pt();
282 Double_t theD0eta0 = mcD0PartDaughter0->Eta();
283 Double_t theD0eta1 = mcD0PartDaughter1->Eta();
284
285 // ACCEPTANCE REQUESTS ---------
286
287 // soft pion
288 Bool_t daught1inAcceptance = (TMath::Abs(eta1) <= 0.9 && pt1 >= 0.05);
289 // Do daughters
290 Bool_t D0daught0inAcceptance = (TMath::Abs(theD0eta0) <= 0.9 && theD0pt0 >= 0.1);
291 Bool_t D0daught1inAcceptance = (TMath::Abs(theD0eta1) <= 0.9 && theD0pt1 >= 0.1);
292
293 if (daught1inAcceptance && D0daught0inAcceptance && D0daught1inAcceptance) {
294
295 AliDebug(2, "D* Daughter particles in acceptance");
296
297 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance);
298
299 // check on the vertex
300 if (vtxFlag){
301 // filling the container if the vertex is ok
302 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex) ;
303
304 Bool_t refitFlag = kTRUE;
305 for (Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
306 AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
307
308 // refit only for D0 daughters
309 if ((track->GetLabel() == daughD00) || (track->GetLabel() == daughD01)) {
310 if(!(track->GetStatus()&AliESDtrack::kTPCrefit) || !(track->GetStatus()&AliESDtrack::kITSrefit)) {
311 refitFlag = kFALSE;
312 }
313 }
314 }
315 if (refitFlag){ // refit
316
317 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit);
318
319 } // end of refit
320 } // end of vertex
321 } //end of acceptance
322 } // end of MC D*
323 else {
324 AliDebug(3,"Problems in filling the container");
325 continue;
326 }
327 } // end of MC loop
328
329 //rec
330 AliDebug(2, Form("Found %d D* candidates",arrayDStartoD0pi->GetEntriesFast()));
331
332 //D* and D0 prongs needed to MatchToMC method
333 Int_t pdgDgDStartoD0pi[2]={421,211};
334 Int_t pdgDgD0toKpi[2]={321,211};
335
336 for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
337
338 // D* candidates
339 AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
340
341 // D0 from the reco cascade
342 AliAODRecoDecayHF2Prong* D0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
343 Bool_t unsetvtx=kFALSE;
344
345 // needed for pointing angle
346 if(!D0particle->GetOwnPrimaryVtx()) {
347 D0particle->SetOwnPrimaryVtx(vtx1);
348 unsetvtx=kTRUE;
349 }
350
351 // find associated MC particle for D* ->D0toKpi
352 Int_t mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray);
353
354 // find D0->Kpi ... needed in the following
355 Int_t mcD0Label = D0particle->MatchToMC(421,mcArray,2,pdgDgD0toKpi);
356
357 if (mcLabel == -1 || mcD0Label ==-1)
358 {
359 AliDebug(2,"No MC particle found");
360 }
361 else {
362
363 // the D* and the D0 in MC
364 AliAODMCParticle* mcVtxHFDStar = (AliAODMCParticle*)mcArray->At(mcLabel);
365 AliAODMCParticle* mcVtxHFD0Kpi = (AliAODMCParticle*)mcArray->At(mcD0Label);
366
367 if (!mcVtxHFD0Kpi || !mcVtxHFDStar) {
368 AliWarning("Could not find associated MC D0 and/or D* in AOD MC tree");
369 continue;
370 }
371
372 // soft pion
373 AliAODTrack *track2 = (AliAODTrack*)dstarD0pi->GetBachelor();
374
375 //D0tokpi
376 AliAODTrack *track0 = (AliAODTrack*)D0particle->GetDaughter(0);
377 AliAODTrack *track1 = (AliAODTrack*)D0particle->GetDaughter(1);
378
379 // check if associated MC v0 passes the cuts
380 if (!fCFManager->CheckParticleCuts(0 ,mcVtxHFDStar)) {
381 AliDebug(2, "Skipping the particles due to cuts");
382 continue;
383 }
384
385 // fill the container...
386 Double_t pt = TMath::Sqrt(TMath::Power(D0particle->Pt(),2)+TMath::Power(track2->Pt(),2));
387 Double_t rapidity = dstarD0pi->YDstar();
388 Double_t cosThetaStar = 9999.;
389 Double_t pTpi = 0.;
390 Double_t pTK = 0.;
391 Double_t dca = (D0particle->GetDCA())*1E4;
392 Double_t d0pi = 0.;
393 Double_t d0K = 0.;
394 Double_t d0xd0 = (D0particle->Prodd0d0())*1E8;
395 Double_t cosPointingAngle = D0particle->CosPointingAngle();
396 Double_t phi = D0particle->Phi();
397
398 // Select D0 cutting variables
399 Int_t pdgCode = mcVtxHFD0Kpi->GetPdgCode();
400
401 // D0 related variables
402 if (pdgCode > 0){
403
404 cosThetaStar = D0particle->CosThetaStarD0();
405 pTpi = D0particle->PtProng(0);
406 pTK = D0particle->PtProng(1);
407 d0pi = (D0particle->Getd0Prong(0))*1E4;
408 d0K = (D0particle->Getd0Prong(1))*1E4;
409
410 }
411 else {
412
413 cosThetaStar = D0particle->CosThetaStarD0bar();
414 pTpi = D0particle->PtProng(1);
415 pTK = D0particle->PtProng(0);
416 d0pi = (D0particle->Getd0Prong(1))*1E4;
417 d0K = (D0particle->Getd0Prong(0))*1E4;
418
419 }
420
421 // ct of the D0 from D*
422 Double_t cT = D0particle->CtD0();
423
424 containerInput[0] = pt;
425 containerInput[1] = rapidity;
426 containerInput[2] = cosThetaStar;
427 containerInput[3] = track2->Pt();
428 containerInput[4] = D0particle->Pt();
429 containerInput[5] = cT*1.E4; // in micron
430 containerInput[6] = dca; // in micron
431 containerInput[7] = d0pi; // in micron
432 containerInput[8] = d0K; // in micron
433 containerInput[9] = d0xd0; // in micron^2
434 containerInput[10] = cosPointingAngle; // in micron
435 containerInput[11] = phi;
436 containerInput[12] = zPrimVertex; // z of reconstructed of primary vertex
437 containerInput[13] = pTpi; // D0 pion
438 containerInput[14] = pTK; // D0 kaon
439
440 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;
441
442 // refit in ITS and TPC for D0 daughters
443 if((!(track0->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track1->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track0->GetStatus()&AliESDtrack::kITSrefit)) || (!(track1->GetStatus()&AliESDtrack::kITSrefit))) {
444 continue;
445 }
446
447 // reft in ITS for soft pion
448 if((!(track2->GetStatus()&AliESDtrack::kITSrefit))) {
449 continue;
450 }
451
452 // cut in acceptance for the soft pion and for the D0 daughters
453 Bool_t acceptanceProng0 = (TMath::Abs(D0particle->EtaProng(0))<= 0.9 && D0particle->PtProng(0) >= 0.1);
454 Bool_t acceptanceProng1 = (TMath::Abs(D0particle->EtaProng(1))<= 0.9 && D0particle->PtProng(1) >= 0.1);
455
456 // soft pion acceptance ... is it fine 0.9?????
457 Bool_t acceptanceProng2 = (TMath::Abs(track2->Eta())<= 0.9 && track2->Pt() >= 0.05);
458
459 if (acceptanceProng0 && acceptanceProng1 && acceptanceProng2) {
460 AliDebug(2,"D* reco daughters in acceptance");
461 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance) ;
462
463 if(fAcceptanceUnf){
464 Double_t fill[4]; //fill response matrix
465
466 // dimensions 0&1 : pt,eta (Rec)
467 fill[0] = pt ;
468 fill[1] = rapidity;
469 // dimensions 2&3 : pt,eta (MC)
470 fill[2] = mcVtxHFDStar->Pt();
471 fill[3] = mcVtxHFDStar->Y();
472 fCorrelation->Fill(fill);
473 }
474
475 // cut on the min n. of clusters in ITS for the D0 and soft pion
476 Int_t ncls0=0,ncls1=0,ncls2=0;
477 for(Int_t l=0;l<6;l++) {
478 if(TESTBIT(track0->GetITSClusterMap(),l)) ncls0++;
479 if(TESTBIT(track1->GetITSClusterMap(),l)) ncls1++;
480 if(TESTBIT(track2->GetITSClusterMap(),l)) ncls2++;
481 }
482 // see AddTask for soft pion ITS clusters request
483 AliDebug(2, Form("n clusters = %d", ncls0));
484
485 if (ncls0 >= fMinITSClusters && ncls1 >= fMinITSClusters && ncls2>= fMinITSClustersSoft) {
486 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters) ;
487
488 // D0 cuts optimized for D* analysis
489 Double_t cuts[7] = {9999999., 1.1, 0., 9999999., 9999999., 0.,0.027};
490
491 // needed for cuts
492 Double_t theD0pt = D0particle->Pt();
493
494 if (theD0pt <= 1){ // first bin not optimized
495 cuts[0] = 400;
496 cuts[1] = 0.8;
497 cuts[2] = 0.21;
498 cuts[3] = 500;
499 cuts[4] = 500;
500 cuts[5] = -20000;
501 cuts[6] = 0.6;
502 }
503 else if (theD0pt > 1 && theD0pt <= 2){
504 cuts[0] = 200;
505 cuts[1] = 0.7;
506 cuts[2] = 0.8;
507 cuts[3] = 210;
508 cuts[4] = 210;
509 cuts[5] = -20000;
510 cuts[6] = 0.9;
511 }
512 else if (theD0pt > 2 && theD0pt <= 3){
513 cuts[0] = 400;
514 cuts[1] = 0.8;
515 cuts[2] = 0.8;
516 cuts[3] = 420;
517 cuts[4] = 350;
518 cuts[5] = -8500;
519 cuts[6] = 0.9;
520 }
521 else if (theD0pt > 3 && theD0pt <= 5){
522 cuts[0] = 160;
523 cuts[1] = 1.0;
524 cuts[2] = 1.2;
525 cuts[3] = 560;
526 cuts[4] = 420;
527 cuts[5] = -8500;
528 cuts[6] = 0.9;
529 }
530 else if (theD0pt > 5){
531 cuts[0] = 800;
532 cuts[1] = 1.0;
533 cuts[2] = 1.2;
534 cuts[3] = 700;
535 cuts[4] = 700;
536 cuts[5] = 10000;
537 cuts[6] = 0.9;
538 }
539 if (dca < cuts[0]
540 && TMath::Abs(cosThetaStar) < cuts[1]
541 && pTpi > cuts[2]
542 && pTK > cuts[2]
543 && TMath::Abs(d0pi) < cuts[3]
544 && TMath::Abs(d0K) < cuts[4]
545 && d0xd0 < cuts[5]
546 && cosPointingAngle > cuts[6]
547 ){
548
549 AliDebug(2,"Particle passed D* selection cuts");
550 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoCuts) ;
551
552 if(!fAcceptanceUnf){ // unfolding
553
554 Double_t fill[4]; //fill response matrix
555
556 // dimensions 0&1 : pt,eta (Rec)
557 fill[0] = pt ;
558 fill[1] = rapidity;
559 // dimensions 2&3 : pt,eta (MC)
560 fill[2] = mcVtxHFDStar->Pt();
561 fill[3] = mcVtxHFDStar->Y();
562
563 fCorrelation->Fill(fill);
564 }
565 }
566 }
567 }
568 }
569 if(unsetvtx) D0particle->UnsetOwnPrimaryVtx();
570 } // end loop on D*->Kpipi
571
572 fHistEventsProcessed->Fill(0);
573
574 PostData(1,fHistEventsProcessed) ;
575 PostData(2,fCFManager->GetParticleContainer()) ;
576 PostData(3,fCorrelation) ;
577}
578
579
580//___________________________________________________________________________
581void AliCFTaskForDStarAnalysis::Terminate(Option_t*)
582{
583 // The Terminate()
584 AliAnalysisTaskSE::Terminate();
585
586 // draw correlation matrix
587 AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
588 if(!cont) {
589 printf("CONTAINER NOT FOUND\n");
590 return;
591 }
592}
593
594//___________________________________________________________________________
595void AliCFTaskForDStarAnalysis::UserCreateOutputObjects() {
596
597 //
598 // useroutput
599 //
600
601 Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
602
603 //slot #1
604 OpenFile(1);
605 fHistEventsProcessed = new TH1I("CFDSchist0","",1,0,1) ;
606}
607
608//________________________________________________________________________________
609Bool_t AliCFTaskForDStarAnalysis::GetDStarMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC)const {
610
611 //
612 // fill the D* and D0 MC container
613 //
614
615 Bool_t isDStar = kFALSE;
616
617 if(TMath::Abs(mcPart->GetPdgCode())!=413) return isDStar;
618
619 // getting the daughters
620 Int_t daughter0 = mcPart->GetDaughter(0);
621 Int_t daughter1 = mcPart->GetDaughter(1);
622
623 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
624 if (daughter0 == 0 || daughter1 == 0) {
625 AliDebug(2, "Error! the D* MC doesn't have correct daughters!!");
626 return isDStar;
627 }
628
629 if (TMath::Abs(daughter1 - daughter0) != 1) { // should be everytime true - see PDGdatabooklet
630 AliDebug(2, "The D* MC doesn't come from a 2-prong decay, skipping!!");
631 return isDStar;
632 }
633
634 AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
635 AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
636 if (!mcPartDaughter0 || !mcPartDaughter1) {
637 AliWarning("D*: At least one Daughter Particle not found in tree, skipping");
638 return isDStar;
639 }
640
641 if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==421 &&
642 TMath::Abs(mcPartDaughter1->GetPdgCode())==211) &&
643 !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
644 TMath::Abs(mcPartDaughter1->GetPdgCode())==421)) {
645 AliDebug(2, "The D* MC doesn't come from a Kpi decay, skipping!!");
646 return isDStar;
647 }
648
649 Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0
650 Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1
651
652 // getting vertex from daughters
653 mcPartDaughter0->XvYvZv(vtx2daughter0);
654 mcPartDaughter1->XvYvZv(vtx2daughter1);
655
656 // check if the secondary vertex is the same for both
657 if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
658 AliError("The D* daughters have different secondary vertex, skipping the track");
659 return isDStar;
660 }
661
662 AliAODMCParticle* neutralDaugh = mcPartDaughter0;
663
664 Double_t VectorD0[2] ={0.,0.};
665
666 if (!EvaluateIfD0toKpi(neutralDaugh,mcArray,VectorD0)) {
667 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
668 return isDStar;
669 }
670 // get the pT of the daughters
671
672 Double_t pTpi = 0.;
673 Double_t pTD0 = 0.;
674
675 if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
676 pTpi = mcPartDaughter0->Pt();
677 pTD0 = mcPartDaughter1->Pt();
678 }
679 else {
680 pTpi = mcPartDaughter1->Pt();
681 pTD0 = mcPartDaughter0->Pt();
682 }
683
684 vectorMC[0] = mcPart->Pt();
685 vectorMC[1] = mcPart->Y() ;
686 vectorMC[2] = 0;
687 vectorMC[3] = pTpi ;
688 vectorMC[4] = pTD0 ;
689 vectorMC[5] = 0;
690 vectorMC[6] = mcPart->Phi() ;
691 vectorMC[7] = VectorD0[0] ;
692 vectorMC[8] = VectorD0[1] ;
693
694 isDStar = kTRUE;
695
696 return isDStar;
697}
698//________________________________________________________________________________________________
699
700Bool_t AliCFTaskForDStarAnalysis::EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh, TClonesArray* mcArray, Double_t* VectorD0)const{
701
702 //
703 // chack wether D0 is decaing into kpi
704 //
705
706 Bool_t isHadronic = kFALSE;
707
708 Int_t daughterD00 = neutralDaugh->GetDaughter(0);
709 Int_t daughterD01 = neutralDaugh->GetDaughter(1);
710
711 AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughterD00,daughterD01));
712 if (daughterD00 == 0 || daughterD01 == 0) {
713 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
714 return isHadronic;
715 }
716
717 if (TMath::Abs(daughterD01 - daughterD00) != 1) { // should be everytime true - see PDGdatabooklet
718 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
719 return isHadronic;
720 }
721
722 AliAODMCParticle* mcPartDaughterD00 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD00));
723 AliAODMCParticle* mcPartDaughterD01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD01));
724 if (!mcPartDaughterD00 || !mcPartDaughterD01) {
725 AliWarning("D0 MC analysis: At least one Daughter Particle not found in tree, skipping");
726 return isHadronic;
727 }
728
729 if (!(TMath::Abs(mcPartDaughterD00->GetPdgCode())==321 &&
730 TMath::Abs(mcPartDaughterD01->GetPdgCode())==211) &&
731 !(TMath::Abs(mcPartDaughterD00->GetPdgCode())==211 &&
732 TMath::Abs(mcPartDaughterD01->GetPdgCode())==321)) {
733 AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
734 return isHadronic;
735 }
736
737 Double_t pTD0pi = 0;
738 Double_t pTD0K = 0;
739
740
741 if (TMath::Abs(mcPartDaughterD00->GetPdgCode()) == 211) {
742 pTD0pi = mcPartDaughterD00->Pt();
743 pTD0K = mcPartDaughterD01->Pt();
744 }
745 else {
746 pTD0pi = mcPartDaughterD01->Pt();
747 pTD0K = mcPartDaughterD00->Pt();
748 }
749
750 isHadronic = kTRUE;
751
752 VectorD0[0] = pTD0pi;
753 VectorD0[1] = pTD0K;
754
755 return isHadronic;
756
757}