Added LHC10h run list for flow analysis (Giacomo)
[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 /* $Id$ */
17
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 //__________________________________________________________________________
59 AliCFTaskForDStarAnalysis::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 //___________________________________________________________________________
75 AliCFTaskForDStarAnalysis::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 //___________________________________________________________________________
97 AliCFTaskForDStarAnalysis& 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 //___________________________________________________________________________
111 AliCFTaskForDStarAnalysis::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 //___________________________________________________________________________
128 AliCFTaskForDStarAnalysis::~AliCFTaskForDStarAnalysis() {
129   //
130   //destructor
131   //
132   if (fCFManager)           delete fCFManager ;
133   if (fHistEventsProcessed) delete fHistEventsProcessed ;
134   if (fCorrelation)         delete fCorrelation ;
135 }
136
137 //_________________________________________________
138 void 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   
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   
175   // fix for temporary bug in ESDfilter 
176   // the AODs with null vertex pointer didn't pass the PhysSel
177   if(!aodEvent->GetPrimaryVertex() || TMath::Abs(aodEvent->GetMagneticField())<0.001) return;
178
179   fEvents++;
180   if (fEvents%10000 ==0) AliDebug(2,Form("Event %d",fEvents));
181
182   fCFManager->SetRecEventInfo(aodEvent);
183   fCFManager->SetMCEventInfo(aodEvent);
184   
185   // event selection
186   Double_t containerInput[15] ;
187   Double_t containerInputMC[15] ;
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));
254       if (!mcPartDaughter0) {
255         AliError("Could not find Monte-Carlo in AOD");
256         return;
257       }
258   
259       // Soft Pion
260       AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
261       if (!mcPartDaughter1) {
262         AliError("Could not find Monte-Carlo in AOD");
263         return;
264       }
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) {
287         AliWarning("At least one D0 Daughter Particle not found in tree, but it should be, this check was already done...");
288         return; 
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 //___________________________________________________________________________
593 void 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 //___________________________________________________________________________
607 void 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 //________________________________________________________________________________
621 Bool_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
712 Bool_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 }