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