Coverity
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliCFTaskForDStarAnalysis.cxx
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 //__________________________________________________________________________
57 AliCFTaskForDStarAnalysis::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 //___________________________________________________________________________
73 AliCFTaskForDStarAnalysis::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 //___________________________________________________________________________
95 AliCFTaskForDStarAnalysis& 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 //___________________________________________________________________________
109 AliCFTaskForDStarAnalysis::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 //___________________________________________________________________________
126 AliCFTaskForDStarAnalysis::~AliCFTaskForDStarAnalysis() {
127   //
128   //destructor
129   //
130   if (fCFManager)           delete fCFManager ;
131   if (fHistEventsProcessed) delete fHistEventsProcessed ;
132   if (fCorrelation)         delete fCorrelation ;
133 }
134
135 //_________________________________________________
136 void 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   
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   
173   // fix for temporary bug in ESDfilter 
174   // the AODs with null vertex pointer didn't pass the PhysSel
175   if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
176
177   fEvents++;
178   if (fEvents%10000 ==0) AliDebug(2,Form("Event %d",fEvents));
179
180   fCFManager->SetRecEventInfo(aodEvent);
181   fCFManager->SetMCEventInfo(aodEvent);
182   
183   // event selection
184   Double_t containerInput[15] ;
185   Double_t containerInputMC[15] ;
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       if (!mcPartDaughter0) {
253         AliError("Could not find Monte-Carlo in AOD");
254         return;
255       }
256   
257       // Soft Pion
258       AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
259       if (!mcPartDaughter1) {
260         AliError("Could not find Monte-Carlo in AOD");
261         return;
262       }
263       
264       // Acceptance variables for the soft pion
265       Double_t eta1 = mcPartDaughter1->Eta();
266       Double_t pt1  = mcPartDaughter1->Pt(); 
267       
268       Int_t daughD00 = 0;
269       Int_t daughD01 = 0;
270       
271       // Just to be sure to take the right particles    
272       if(TMath::Abs(mcPartDaughter0->GetPdgCode())==421){
273         daughD00 = mcPartDaughter0->GetDaughter(0);
274         daughD01 = mcPartDaughter0->GetDaughter(1);
275       }else{
276         daughD00 = mcPartDaughter1->GetDaughter(0);
277         daughD01 = mcPartDaughter1->GetDaughter(1);
278       }
279       
280       // the D0 daughters
281       AliAODMCParticle* mcD0PartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughD00));
282       AliAODMCParticle* mcD0PartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughD01));
283       
284       if (!mcD0PartDaughter0 || !mcD0PartDaughter1) {
285         AliWarning("At least one D0 Daughter Particle not found in tree, but it should be, this check was already done...");
286         return; 
287       }
288       
289       // D0 daughters - needed for acceptance
290       Double_t theD0pt0 =  mcD0PartDaughter0->Pt();
291       Double_t theD0pt1 =  mcD0PartDaughter1->Pt();
292       Double_t theD0eta0 = mcD0PartDaughter0->Eta();
293       Double_t theD0eta1 = mcD0PartDaughter1->Eta();
294       
295       // ACCEPTANCE REQUESTS ---------
296       
297       // soft pion 
298       Bool_t daught1inAcceptance = (TMath::Abs(eta1) <= 0.9 && pt1 >= 0.05);
299       // Do daughters
300       Bool_t D0daught0inAcceptance = (TMath::Abs(theD0eta0) <= 0.9 && theD0pt0 >= 0.1); 
301       Bool_t D0daught1inAcceptance = (TMath::Abs(theD0eta1) <= 0.9 && theD0pt1 >= 0.1);
302       
303       if (daught1inAcceptance && D0daught0inAcceptance && D0daught1inAcceptance) {
304         
305         AliDebug(2, "D* Daughter particles in acceptance");
306         
307         fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance);
308         
309         // check on the vertex
310         if (vtxFlag){
311           // filling the container if the vertex is ok
312           fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex) ;
313           
314           Bool_t refitFlag = kTRUE;
315           for (Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
316             AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
317             
318             // refit only for D0 daughters
319             if ((track->GetLabel() == daughD00) || (track->GetLabel() == daughD01)) {
320               if(!(track->GetStatus()&AliESDtrack::kTPCrefit) || !(track->GetStatus()&AliESDtrack::kITSrefit)) {
321                 refitFlag = kFALSE;
322               }
323             }
324           } 
325           if (refitFlag){ // refit
326       
327             fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit);
328
329           } // end of refit       
330         } // end of vertex
331       } //end of acceptance             
332     } // end of MC D*
333     else {
334       AliDebug(3,"Problems in filling the container");
335       continue;
336     }
337   } // end of MC loop
338   
339   //rec
340   AliDebug(2, Form("Found %d D* candidates",arrayDStartoD0pi->GetEntriesFast()));
341   
342   //D* and D0 prongs needed to MatchToMC method
343   Int_t pdgDgDStartoD0pi[2]={421,211};
344   Int_t pdgDgD0toKpi[2]={321,211};
345   
346   for (Int_t iDStartoD0pi = 0; iDStartoD0pi<arrayDStartoD0pi->GetEntriesFast(); iDStartoD0pi++) {
347     
348     // D* candidates
349     AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)arrayDStartoD0pi->At(iDStartoD0pi);
350     
351     // D0 from the reco cascade
352     AliAODRecoDecayHF2Prong* D0particle = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
353     Bool_t unsetvtx=kFALSE;
354     
355     // needed for pointing angle
356     if(!D0particle->GetOwnPrimaryVtx()) {
357       D0particle->SetOwnPrimaryVtx(vtx1);
358       unsetvtx=kTRUE;
359     }
360     
361     // find associated MC particle for D* ->D0toKpi
362     Int_t mcLabel = dstarD0pi->MatchToMC(413,421,pdgDgDStartoD0pi,pdgDgD0toKpi,mcArray); 
363     
364     // find D0->Kpi ... needed in the following
365     Int_t mcD0Label = D0particle->MatchToMC(421,mcArray,2,pdgDgD0toKpi); 
366     
367     if (mcLabel == -1 || mcD0Label ==-1) 
368       {
369         AliDebug(2,"No MC particle found");
370       }
371     else {
372       
373       // the D* and the D0 in MC
374       AliAODMCParticle* mcVtxHFDStar = (AliAODMCParticle*)mcArray->At(mcLabel);
375       AliAODMCParticle* mcVtxHFD0Kpi = (AliAODMCParticle*)mcArray->At(mcD0Label);
376       
377       if (!mcVtxHFD0Kpi || !mcVtxHFDStar) {
378         AliWarning("Could not find associated MC D0 and/or D* in AOD MC tree");
379         continue;
380       }
381       
382       // soft pion
383       AliAODTrack *track2 = (AliAODTrack*)dstarD0pi->GetBachelor(); 
384
385       //D0tokpi
386       AliAODTrack *track0 = (AliAODTrack*)D0particle->GetDaughter(0);
387       AliAODTrack *track1 = (AliAODTrack*)D0particle->GetDaughter(1);
388       
389       // check if associated MC v0 passes the cuts
390       if (!fCFManager->CheckParticleCuts(0 ,mcVtxHFDStar)) { 
391         AliDebug(2, "Skipping the particles due to cuts");
392         continue; 
393       }
394       
395       // fill the container...
396       Double_t pt = TMath::Sqrt(TMath::Power(D0particle->Pt(),2)+TMath::Power(track2->Pt(),2));
397       Double_t rapidity = dstarD0pi->YDstar();
398       Double_t cosThetaStar = 9999.;
399       Double_t pTpi  = 0.;
400       Double_t pTK   = 0.;
401       Double_t dca   = (D0particle->GetDCA())*1E4;
402       Double_t d0pi  = 0.;
403       Double_t d0K   = 0.;
404       Double_t d0xd0 = (D0particle->Prodd0d0())*1E8;
405       Double_t cosPointingAngle = D0particle->CosPointingAngle();
406       Double_t phi = D0particle->Phi();
407            
408       // Select D0 cutting variables
409       Int_t pdgCode = mcVtxHFD0Kpi->GetPdgCode();
410
411       // D0 related variables
412       if (pdgCode > 0){
413
414         cosThetaStar = D0particle->CosThetaStarD0();
415         pTpi    = D0particle->PtProng(0);
416         pTK     = D0particle->PtProng(1);
417         d0pi    = (D0particle->Getd0Prong(0))*1E4;
418         d0K     = (D0particle->Getd0Prong(1))*1E4;
419    
420       }
421       else {
422
423         cosThetaStar = D0particle->CosThetaStarD0bar();
424         pTpi    = D0particle->PtProng(1);
425         pTK     = D0particle->PtProng(0);
426         d0pi    = (D0particle->Getd0Prong(1))*1E4;
427         d0K     = (D0particle->Getd0Prong(0))*1E4;
428    
429       }
430       
431       // ct of the D0 from D*
432       Double_t cT = D0particle->CtD0();
433       
434       containerInput[0]  = pt;
435       containerInput[1]  = rapidity;
436       containerInput[2]  = cosThetaStar;
437       containerInput[3]  = track2->Pt();
438       containerInput[4]  = D0particle->Pt();
439       containerInput[5]  = cT*1.E4;  // in micron
440       containerInput[6]  = dca;  // in micron
441       containerInput[7]  = d0pi;  // in micron
442       containerInput[8]  = d0K;  // in micron
443       containerInput[9]  = d0xd0;  // in micron^2
444       containerInput[10] = cosPointingAngle;  // in micron
445       containerInput[11] = phi;  
446       containerInput[12] = zPrimVertex;    // z of reconstructed of primary vertex
447       containerInput[13] = pTpi;  // D0 pion
448       containerInput[14] = pTK;   // D0 kaon  
449
450       fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;   
451       
452       // refit in ITS and TPC for D0 daughters
453       if((!(track0->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track1->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track0->GetStatus()&AliESDtrack::kITSrefit)) || (!(track1->GetStatus()&AliESDtrack::kITSrefit))) {
454         continue;
455       }
456
457       // reft in ITS for soft pion
458       if((!(track2->GetStatus()&AliESDtrack::kITSrefit))) {
459         continue;
460       }
461
462       // cut in acceptance for the soft pion and for the D0 daughters      
463       Bool_t acceptanceProng0 = (TMath::Abs(D0particle->EtaProng(0))<= 0.9 && D0particle->PtProng(0) >= 0.1);
464       Bool_t acceptanceProng1 = (TMath::Abs(D0particle->EtaProng(1))<= 0.9 && D0particle->PtProng(1) >= 0.1);
465       
466       // soft pion acceptance ... is it fine 0.9?????
467       Bool_t acceptanceProng2 = (TMath::Abs(track2->Eta())<= 0.9 && track2->Pt() >= 0.05);
468       
469       if (acceptanceProng0 && acceptanceProng1 && acceptanceProng2) {
470         AliDebug(2,"D* reco daughters in acceptance");
471         fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance) ;
472         
473         if(fAcceptanceUnf){       
474           Double_t fill[4]; //fill response matrix
475           
476           // dimensions 0&1 : pt,eta (Rec)              
477           fill[0] = pt ;
478           fill[1] = rapidity;           
479           // dimensions 2&3 : pt,eta (MC)                                       
480           fill[2] =  mcVtxHFDStar->Pt();
481           fill[3] =  mcVtxHFDStar->Y();   
482           fCorrelation->Fill(fill);             
483         }  
484         
485         // cut on the min n. of clusters in ITS for the D0 and soft pion
486         Int_t ncls0=0,ncls1=0,ncls2=0;
487         for(Int_t l=0;l<6;l++) {
488           if(TESTBIT(track0->GetITSClusterMap(),l)) ncls0++;
489           if(TESTBIT(track1->GetITSClusterMap(),l)) ncls1++;
490           if(TESTBIT(track2->GetITSClusterMap(),l)) ncls2++;
491         }       
492         // see AddTask for soft pion ITS clusters request
493         AliDebug(2, Form("n clusters = %d", ncls0));
494
495         if (ncls0 >= fMinITSClusters && ncls1 >= fMinITSClusters && ncls2>= fMinITSClustersSoft) {
496           fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters) ;
497           
498           // D0 cuts optimized for D* analysis
499           Double_t cuts[7] = {9999999., 1.1, 0., 9999999., 9999999., 0.,0.027}; 
500
501           // needed for cuts      
502           Double_t theD0pt = D0particle->Pt();
503
504           if (theD0pt <= 1){ // first bin not optimized
505             cuts[0] = 400;
506             cuts[1] = 0.8;
507             cuts[2] = 0.21;
508             cuts[3] = 500;
509             cuts[4] = 500;
510             cuts[5] = -20000;
511             cuts[6] = 0.6;  
512           }
513           else if (theD0pt > 1 && theD0pt <= 2){
514             cuts[0] = 200; 
515             cuts[1] = 0.7; 
516             cuts[2] = 0.8; 
517             cuts[3] = 210;
518             cuts[4] = 210;
519             cuts[5] = -20000;
520             cuts[6] = 0.9;  
521           }
522           else if (theD0pt > 2 && theD0pt <= 3){
523             cuts[0] = 400;
524             cuts[1] = 0.8; 
525             cuts[2] = 0.8;
526             cuts[3] = 420;
527             cuts[4] = 350; 
528             cuts[5] = -8500;
529             cuts[6] = 0.9;   
530           }
531           else if (theD0pt > 3 && theD0pt <= 5){
532             cuts[0] = 160;  
533             cuts[1] = 1.0; 
534             cuts[2] = 1.2;  
535             cuts[3] = 560; 
536             cuts[4] = 420; 
537             cuts[5] = -8500;
538             cuts[6] = 0.9;  
539           }
540           else if (theD0pt > 5){
541             cuts[0] = 800;
542             cuts[1] = 1.0;
543             cuts[2] = 1.2; 
544             cuts[3] = 700;
545             cuts[4] = 700; 
546             cuts[5] = 10000;
547             cuts[6] = 0.9;  
548           }
549           if (dca < cuts[0] 
550               && TMath::Abs(cosThetaStar) < cuts[1]  
551               && pTpi > cuts[2] 
552               && pTK > cuts[2]  
553               && TMath::Abs(d0pi) < cuts[3] 
554               && TMath::Abs(d0K) < cuts[4]  
555               && d0xd0 < cuts[5] 
556               && cosPointingAngle > cuts[6]
557             ){
558             
559             AliDebug(2,"Particle passed D* selection cuts");
560             fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoCuts) ;   
561             
562             if(!fAcceptanceUnf){ // unfolding
563               
564               Double_t fill[4]; //fill response matrix
565               
566               // dimensions 0&1 : pt,eta (Rec)              
567               fill[0] = pt ;
568               fill[1] = rapidity;                   
569               // dimensions 2&3 : pt,eta (MC)               
570               fill[2] =  mcVtxHFDStar->Pt();
571               fill[3] =  mcVtxHFDStar->Y();
572               
573               fCorrelation->Fill(fill);             
574             }
575           }             
576         }
577       }
578     }
579     if(unsetvtx) D0particle->UnsetOwnPrimaryVtx();
580   } // end loop on D*->Kpipi
581   
582   fHistEventsProcessed->Fill(0);
583   
584   PostData(1,fHistEventsProcessed) ;
585   PostData(2,fCFManager->GetParticleContainer()) ;
586   PostData(3,fCorrelation) ;
587 }
588
589
590 //___________________________________________________________________________
591 void AliCFTaskForDStarAnalysis::Terminate(Option_t*)
592 {
593   // The Terminate()    
594   AliAnalysisTaskSE::Terminate();
595   
596   // draw correlation matrix
597   AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
598   if(!cont) {
599     printf("CONTAINER NOT FOUND\n");
600     return;
601   } 
602 }
603
604 //___________________________________________________________________________
605 void AliCFTaskForDStarAnalysis::UserCreateOutputObjects() {
606
607   //
608   // useroutput
609   //
610   
611   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
612   
613   //slot #1
614   OpenFile(1);
615   fHistEventsProcessed = new TH1I("CFDSchist0","",1,0,1) ;
616 }
617
618 //________________________________________________________________________________
619 Bool_t AliCFTaskForDStarAnalysis::GetDStarMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC)const {
620   
621   // 
622   // fill the D* and D0 MC container
623   //
624   
625   Bool_t isDStar = kFALSE;
626   
627   if(TMath::Abs(mcPart->GetPdgCode())!=413) return isDStar;
628   
629   // getting the daughters
630   Int_t daughter0 = mcPart->GetDaughter(0);
631   Int_t daughter1 = mcPart->GetDaughter(1);
632   
633   AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
634   if (daughter0 == 0 || daughter1 == 0) {
635     AliDebug(2, "Error! the D* MC doesn't have correct daughters!!");
636     return isDStar;  
637   }
638   
639   if (TMath::Abs(daughter1 - daughter0) != 1) { // should be everytime true - see PDGdatabooklet
640     AliDebug(2, "The D* MC doesn't come from a 2-prong decay, skipping!!");
641     return isDStar;  
642   }
643   
644   AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
645   AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
646   if (!mcPartDaughter0 || !mcPartDaughter1) {
647     AliWarning("D*: At least one Daughter Particle not found in tree, skipping"); 
648     return isDStar;  
649   }
650   
651   if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==421 &&
652         TMath::Abs(mcPartDaughter1->GetPdgCode())==211) && 
653       !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
654         TMath::Abs(mcPartDaughter1->GetPdgCode())==421)) {
655     AliDebug(2, "The D* MC doesn't come from a Kpi decay, skipping!!");
656     return isDStar;  
657   }
658   
659   Double_t vtx2daughter0[3] = {0,0,0};   // secondary vertex from daughter 0
660   Double_t vtx2daughter1[3] = {0,0,0};   // secondary vertex from daughter 1
661
662   // getting vertex from daughters
663   mcPartDaughter0->XvYvZv(vtx2daughter0);  
664   mcPartDaughter1->XvYvZv(vtx2daughter1);  
665   
666   // check if the secondary vertex is the same for both
667   if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
668     AliError("The D* daughters have different secondary vertex, skipping the track");
669     return isDStar;
670   }
671   
672   AliAODMCParticle* neutralDaugh = mcPartDaughter0;
673   
674   Double_t VectorD0[2] ={0.,0.};
675   
676   if (!EvaluateIfD0toKpi(neutralDaugh,mcArray,VectorD0)) {
677     AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
678     return isDStar;  
679   }
680   // get the pT of the daughters
681   
682   Double_t pTpi = 0.;
683   Double_t pTD0 = 0.;
684   
685   if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
686     pTpi = mcPartDaughter0->Pt();
687     pTD0 = mcPartDaughter1->Pt();
688   }
689   else {
690     pTpi = mcPartDaughter1->Pt();
691     pTD0 = mcPartDaughter0->Pt();
692   }
693   
694   vectorMC[0] = mcPart->Pt();
695   vectorMC[1] = mcPart->Y() ;
696   vectorMC[2] = 0;
697   vectorMC[3] = pTpi ;
698   vectorMC[4] = pTD0 ;
699   vectorMC[5] = 0;
700   vectorMC[6] = mcPart->Phi() ;
701   vectorMC[7] = VectorD0[0] ;
702   vectorMC[8] = VectorD0[1] ;
703   
704   isDStar = kTRUE;
705   
706   return isDStar;
707 }
708 //________________________________________________________________________________________________
709
710 Bool_t AliCFTaskForDStarAnalysis::EvaluateIfD0toKpi(AliAODMCParticle* neutralDaugh, TClonesArray* mcArray, Double_t* VectorD0)const{
711   
712   //
713   // chack wether D0 is decaing into kpi
714   //
715   
716   Bool_t isHadronic = kFALSE;
717   
718   Int_t daughterD00 = neutralDaugh->GetDaughter(0);
719   Int_t daughterD01 = neutralDaugh->GetDaughter(1);
720   
721   AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughterD00,daughterD01));
722   if (daughterD00 == 0 || daughterD01 == 0) {
723     AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
724     return isHadronic;  
725   }
726   
727   if (TMath::Abs(daughterD01 - daughterD00) != 1) { // should be everytime true - see PDGdatabooklet
728     AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
729     return isHadronic;  
730   }
731   
732   AliAODMCParticle* mcPartDaughterD00 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD00));
733   AliAODMCParticle* mcPartDaughterD01 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughterD01));
734   if (!mcPartDaughterD00 || !mcPartDaughterD01) {
735     AliWarning("D0 MC analysis: At least one Daughter Particle not found in tree, skipping"); 
736     return isHadronic;  
737   }
738   
739   if (!(TMath::Abs(mcPartDaughterD00->GetPdgCode())==321 &&
740         TMath::Abs(mcPartDaughterD01->GetPdgCode())==211) && 
741       !(TMath::Abs(mcPartDaughterD00->GetPdgCode())==211 &&
742         TMath::Abs(mcPartDaughterD01->GetPdgCode())==321)) {
743     AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
744     return isHadronic;  
745   }
746
747   Double_t pTD0pi = 0;
748   Double_t pTD0K = 0;
749
750  
751   if (TMath::Abs(mcPartDaughterD00->GetPdgCode()) == 211) {
752     pTD0pi = mcPartDaughterD00->Pt();
753     pTD0K = mcPartDaughterD01->Pt();
754   }
755   else {
756     pTD0pi = mcPartDaughterD01->Pt();
757     pTD0K  = mcPartDaughterD00->Pt();
758   }
759  
760   isHadronic = kTRUE;
761
762   VectorD0[0] = pTD0pi;
763   VectorD0[1] = pTD0K;
764
765   return isHadronic;
766
767 }