Use SpecialOutput for the CF container (Chiara Z)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliCFHeavyFlavourTaskMultiVarMultiStep.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 HF corrections as a function of many variables
18 // 6 Steps introduced: MC, MC Acc, Reco, Reco Acc, Reco Acc + ITS Cl, 
19 // Reco Acc + ITS Cl + PPR cuts
20 // 12 variables used: pt, y, cosThetaStar, ptPi, ptK, ct,
21 // dca, d0Pi, d0K, d0Pixd0K, cosPointingAngle, phi
22 //
23 //-----------------------------------------------------------------------
24 // Author : C. Zampolli, CERN
25 //-----------------------------------------------------------------------
26 //-----------------------------------------------------------------------
27 // Base class for HF Unfolding (pt and eta)
28 // correlation matrix filled at Acceptance and PPR level
29 // Author: A.Grelli ,  Utrecht - agrelli@uu.nl
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 "AliCFHeavyFlavourTaskMultiVarMultiStep.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 "AliAODRecoDecayHF2Prong.h"
50 #include "AliAODMCParticle.h"
51 #include "AliAODMCHeader.h"
52 #include "AliESDtrack.h"
53 #include "TChain.h"
54 #include "THnSparse.h"
55 #include "TH2D.h"
56
57 //__________________________________________________________________________
58 AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep() :
59         AliAnalysisTaskSE(),
60         fPDG(0),
61         fCFManager(0x0),
62         fHistEventsProcessed(0x0),
63         fCorrelation(0x0),
64         fCountMC(0),
65         fCountAcc(0),
66         fCountVertex(0),
67         fCountRefit(0),
68         fCountReco(0),
69         fCountRecoAcc(0),
70         fCountRecoITSClusters(0),
71         fCountRecoPPR(0),
72         fEvents(0),
73         fFillFromGenerated(kFALSE),
74         fMinITSClusters(5),
75         fAcceptanceUnf(kTRUE)
76 {
77         //
78         //Default ctor
79         //
80 }
81 //___________________________________________________________________________
82 AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(const Char_t* name) :
83         AliAnalysisTaskSE(name),
84         fPDG(0),
85         fCFManager(0x0),
86         fHistEventsProcessed(0x0),
87         fCorrelation(0x0),
88         fCountMC(0),
89         fCountAcc(0),
90         fCountVertex(0),
91         fCountRefit(0),
92         fCountReco(0),
93         fCountRecoAcc(0),
94         fCountRecoITSClusters(0),
95         fCountRecoPPR(0),
96         fEvents(0),
97         fFillFromGenerated(kFALSE),
98         fMinITSClusters(5),
99         fAcceptanceUnf(kTRUE)
100 {
101         //
102         // Constructor. Initialization of Inputs and Outputs
103         //
104         Info("AliCFHeavyFlavourTaskMultiVarMultiStep","Calling Constructor");
105         /*
106           DefineInput(0) and DefineOutput(0)
107           are taken care of by AliAnalysisTaskSE constructor
108         */
109         DefineOutput(1,TH1I::Class());
110         DefineOutput(2,AliCFContainer::Class());
111         DefineOutput(3,THnSparseD::Class());
112 }
113
114 //___________________________________________________________________________
115 AliCFHeavyFlavourTaskMultiVarMultiStep& AliCFHeavyFlavourTaskMultiVarMultiStep::operator=(const AliCFHeavyFlavourTaskMultiVarMultiStep& c) 
116 {
117         //
118         // Assignment operator
119         //
120         if (this!=&c) {
121                 AliAnalysisTaskSE::operator=(c) ;
122                 fPDG      = c.fPDG;
123                 fCFManager  = c.fCFManager;
124                 fHistEventsProcessed = c.fHistEventsProcessed;
125         }
126         return *this;
127 }
128
129 //___________________________________________________________________________
130 AliCFHeavyFlavourTaskMultiVarMultiStep::AliCFHeavyFlavourTaskMultiVarMultiStep(const AliCFHeavyFlavourTaskMultiVarMultiStep& c) :
131         AliAnalysisTaskSE(c),
132         fPDG(c.fPDG),
133         fCFManager(c.fCFManager),
134         fHistEventsProcessed(c.fHistEventsProcessed),
135         fCorrelation(c.fCorrelation),
136         fCountMC(c.fCountMC),
137         fCountAcc(c.fCountAcc),
138         fCountVertex(c.fCountVertex),
139         fCountRefit(c.fCountRefit),
140         fCountReco(c.fCountReco),
141         fCountRecoAcc(c.fCountRecoAcc),
142         fCountRecoITSClusters(c.fCountRecoITSClusters),
143         fCountRecoPPR(c.fCountRecoPPR),
144         fEvents(c.fEvents),
145         fFillFromGenerated(c.fFillFromGenerated),
146         fMinITSClusters(c.fMinITSClusters),
147         fAcceptanceUnf(c.fAcceptanceUnf)
148 {
149         //
150         // Copy Constructor
151         //
152 }
153
154 //___________________________________________________________________________
155 AliCFHeavyFlavourTaskMultiVarMultiStep::~AliCFHeavyFlavourTaskMultiVarMultiStep() {
156         //
157         //destructor
158         //
159         if (fCFManager)           delete fCFManager ;
160         if (fHistEventsProcessed) delete fHistEventsProcessed ;
161         if (fCorrelation) delete fCorrelation ;
162 }
163
164 //_________________________________________________
165 void AliCFHeavyFlavourTaskMultiVarMultiStep::UserExec(Option_t *)
166 {
167         //
168         // Main loop function
169         //
170         
171         if (fFillFromGenerated){
172                 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
173         }
174
175         if (!fInputEvent) {
176                 Error("UserExec","NO EVENT FOUND!");
177                 return;
178         }
179         
180         fEvents++;
181         if (fEvents%10000 ==0) AliDebug(2,Form("Event %d",fEvents));
182         AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
183
184         TClonesArray *arrayD0toKpi=0;
185
186         if(!aodEvent && AODEvent() && IsStandardAOD()) {
187           // In case there is an AOD handler writing a standard AOD, use the AOD 
188           // event in memory rather than the input (ESD) event.    
189           aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
190           // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
191           // have to taken from the AOD event hold by the AliAODExtension
192           AliAODHandler* aodHandler = (AliAODHandler*) 
193             ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
194           if(aodHandler->GetExtensions()) {
195             AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
196             AliAODEvent *aodFromExt = ext->GetAOD();
197             arrayD0toKpi=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
198           }
199         } else {
200           arrayD0toKpi=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
201         }
202
203
204         if (!arrayD0toKpi) {
205           AliError("Could not find array of HF vertices");
206           return;
207         }
208
209
210         fCFManager->SetRecEventInfo(aodEvent);
211         fCFManager->SetMCEventInfo(aodEvent);
212         
213         // MC-event selection
214         Double_t containerInput[12] ;
215         Double_t containerInputMC[12] ;
216         
217         //loop on the MC event
218         
219         TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
220         if (!mcArray) {
221                 AliError("Could not find Monte-Carlo in AOD");
222                 return;
223         }
224         Int_t icountMC = 0;
225         Int_t icountAcc = 0;
226         Int_t icountReco = 0;
227         Int_t icountVertex = 0;
228         Int_t icountRefit = 0;
229         Int_t icountRecoAcc = 0;
230         Int_t icountRecoITSClusters = 0;
231         Int_t icountRecoPPR = 0;
232         
233         AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
234         if (!mcHeader) {
235                 AliError("Could not find MC Header in AOD");
236                 return;
237         }
238
239         Int_t cquarks = 0;
240                 
241         // AOD primary vertex
242         AliAODVertex *vtx1 = (AliAODVertex*)aodEvent->GetPrimaryVertex();
243         Double_t zPrimVertex = vtx1->GetZ();
244         Double_t zMCVertex = mcHeader->GetVtxZ();
245         Bool_t vtxFlag = kTRUE;
246         TString title=vtx1->GetTitle();
247         if(!title.Contains("VertexerTracks")) vtxFlag=kFALSE;
248
249         for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) { 
250                 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
251                 if (mcPart->GetPdgCode() == 4) cquarks++; 
252                 if (mcPart->GetPdgCode() == -4) cquarks++; 
253                 if (!mcPart) {
254                         AliWarning("Particle not found in tree, skipping"); 
255                         continue;
256                 } 
257                 
258                 // check the MC-level cuts
259
260                 if (!fCFManager->CheckParticleCuts(0, mcPart)) continue;  // 0 stands for MC level
261                 Int_t pdgGranma = CheckOrigin(mcPart, mcArray);
262                 Int_t abspdgGranma = TMath::Abs(pdgGranma);
263                 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
264                         AliDebug(2,Form("Particle has a b-meson, or b-baryon mother (pdg code mother = %d )--> not coming from a c-quark, skipping...", pdgGranma));
265                         continue;  // skipping particles that don't come from c quark
266                 }
267                 
268                 //              if (TMath::Abs(pdgGranma)!=4) {
269
270                 // fill the container for Gen-level selection
271                 Double_t vectorMC[7] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.};
272                 if (GetGeneratedValuesFromMCParticle(mcPart, mcArray, vectorMC)){
273                         containerInputMC[0] = vectorMC[0];
274                         containerInputMC[1] = vectorMC[1] ;
275                         containerInputMC[2] = vectorMC[2] ;
276                         containerInputMC[3] = vectorMC[3] ;
277                         containerInputMC[4] = vectorMC[4] ;
278                         containerInputMC[5] = vectorMC[5] ;  // in micron
279                         containerInputMC[6] = 0.;    // dummy value, meaningless in MC, in micron
280                         containerInputMC[7] = 0.;   // dummy value, meaningless in MC, in micron
281                         containerInputMC[8] = 0.;   // dummy value, meaningless in MC, in micron
282                         containerInputMC[9] = -100000.; // dummy value, meaningless in MC, in micron^2
283                         containerInputMC[10] = 1.01;    // dummy value, meaningless in MC
284                         containerInputMC[11] = vectorMC[6];    // dummy value, meaningless in MC
285                         containerInputMC[12] = zMCVertex;    // z of reconstructed of primary vertex
286                         fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGenerated);
287                         icountMC++;
288
289                         // check the MC-Acceptance level cuts
290                         // since standard CF functions are not applicable, using Kine Cuts on daughters
291                         
292                         Int_t daughter0 = mcPart->GetDaughter(0);
293                         Int_t daughter1 = mcPart->GetDaughter(1);
294                         AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
295                         if (daughter0 == 0 || daughter1 == 0) {
296                                 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!! But it should have, this check was already done...");
297                         }
298                         if (TMath::Abs(daughter1 - daughter0) != 1) {
299                                 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, but it should be, this check was already done...");
300                         }
301                         AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
302                         AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
303                         if (!mcPartDaughter0 || !mcPartDaughter1) {
304                                 AliWarning("At least one Daughter Particle not found in tree, but it should be, this check was already done..."); 
305                         }
306                         Double_t eta0 = mcPartDaughter0->Eta();
307                         Double_t eta1 = mcPartDaughter1->Eta();
308                         Double_t y0 = mcPartDaughter0->Y();
309                         Double_t y1 = mcPartDaughter1->Y();
310                         Double_t pt0 = mcPartDaughter0->Pt();
311                         Double_t pt1 = mcPartDaughter1->Pt();
312                         AliDebug(2, Form("Daughter 0: eta = %f, y = %f, pt = %f", eta0, y0, pt0));
313                         AliDebug(2, Form("Daughter 1: eta = %f, y = %f, pt = %f", eta1, y1, pt1));
314                         Bool_t daught0inAcceptance = (TMath::Abs(eta0) < 0.9 && pt0 > 0.1); 
315                         Bool_t daught1inAcceptance = (TMath::Abs(eta1) < 0.9 && pt1 > 0.1); 
316                         if (daught0inAcceptance && daught1inAcceptance) {
317                                 // checking whether the cuts implemented in the CF are equivalent - simply a cross-check
318                                 AliDebug(2, "Daughter particles in acceptance");
319                                 if (!fCFManager->CheckParticleCuts(1, mcPartDaughter0)) {
320                                         AliDebug(2,"Inconsistency with CF cut for daughter 0!");
321                                 } 
322                                 if (!fCFManager->CheckParticleCuts(1, mcPartDaughter1)) {
323                                         AliDebug(2,"Inconsistency with CF cut for daughter 1!");
324                                 } 
325                                 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance);
326                                 icountAcc++;
327
328                                 // check on the vertex
329                                 if (vtxFlag){
330                                         printf("Vertex cut passed\n");
331                                         // filling the container if the vertex is ok
332                                         fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex) ;
333                                         printf("Vertex cut passed 2\n");
334                                         icountVertex++;
335                                         // check on the kTPCrefit and kITSrefit conditions of the daughters
336                                         Bool_t refitFlag = kTRUE;
337                                         for (Int_t iaod =0; iaod<aodEvent->GetNumberOfTracks(); iaod++){
338                                                 AliAODTrack *track = (AliAODTrack*)aodEvent->GetTrack(iaod);
339                                                 if ((track->GetLabel() == daughter0) || (track->GetLabel() == daughter1)) {
340                                                   if(!(track->GetStatus()&AliESDtrack::kTPCrefit) || !(track->GetStatus()&AliESDtrack::kITSrefit)) {
341                                                     refitFlag = kFALSE;
342                                                   }
343                                                 }
344                                         }
345                                         if (refitFlag){
346                                                 printf("Refit cut passed\n");
347                                                 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit);
348                                                 icountRefit++;
349                                         }
350                                         else{
351                                                 AliDebug(3,"Refit cut not passed\n");
352                                         }
353                                 }
354                                 else{
355                                         AliDebug(3,"Vertex cut not passed\n");
356                                 }                       
357                         }
358                         else{
359                                 AliDebug(3,"Acceptance cut not passed\n");
360                         }
361                 }
362                 else {
363                         AliDebug(3,"Problems in filling the container");
364                         continue;
365                 }
366         }
367
368         if (cquarks<2) AliDebug(2,Form("Event found with %d c-quarks", cquarks));
369         
370         AliDebug(2,Form("Found %i MC particles that are D0!!",icountMC));
371         AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Acc cuts!!",icountAcc));
372         AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Vertex cuts!!",icountVertex));
373         AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Refit cuts!!",icountRefit));
374
375         // Now go to rec level
376         fCountMC += icountMC;
377         fCountAcc += icountAcc;
378
379         AliDebug(2, Form("Found %d vertices",arrayD0toKpi->GetEntriesFast()));
380
381         Int_t pdgDgD0toKpi[2]={321,211};
382
383         for (Int_t iD0toKpi = 0; iD0toKpi<arrayD0toKpi->GetEntriesFast(); iD0toKpi++) {
384                 
385                 AliAODRecoDecayHF2Prong* d0tokpi = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->At(iD0toKpi);
386                 Bool_t unsetvtx=kFALSE;
387                 if(!d0tokpi->GetOwnPrimaryVtx()) {
388                   d0tokpi->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
389                   unsetvtx=kTRUE;
390                 }
391
392                 // find associated MC particle
393                 Int_t mcLabel = d0tokpi->MatchToMC(421,mcArray,2,pdgDgD0toKpi) ;
394                 if (mcLabel == -1) 
395                         {
396                                 AliDebug(2,"No MC particle found");
397                         }
398                 else {
399                         AliAODMCParticle* mcVtxHF = (AliAODMCParticle*)mcArray->At(mcLabel);
400                         if (!mcVtxHF) {
401                                 AliWarning("Could not find associated MC in AOD MC tree");
402                                 continue;
403                         }
404                         // check whether the daughters have kTPCrefit and kITSrefit set
405                         AliAODTrack *track0 = (AliAODTrack*)d0tokpi->GetDaughter(0);
406                         AliAODTrack *track1 = (AliAODTrack*)d0tokpi->GetDaughter(1);
407                         if((!(track0->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track1->GetStatus()&AliESDtrack::kTPCrefit)) || (!(track0->GetStatus()&AliESDtrack::kITSrefit)) || (!(track1->GetStatus()&AliESDtrack::kITSrefit))) {
408                                 // skipping if at least one daughter does not have kTPCrefit or kITSrefit
409                                 continue;
410                         }
411
412                         // check on the vertex -- could be moved outside the loop on the reconstructed D0... 
413                         if(!vtxFlag) {
414                                 // skipping cause vertex is not ok
415                                 continue;
416                         }
417
418                         const Double_t d0tokpiCuts[9] = {0.3,999999.,1.1,0.,0.,999999.,999999.,999999.,0.};
419                         Int_t okD0, okD0bar;
420                         if (!(d0tokpi->SelectD0(&d0tokpiCuts[0],okD0,okD0bar))){
421                                 // skipping candidate
422                                 continue;
423                         } 
424
425                         // check if associated MC v0 passes the cuts
426                         if (!fCFManager->CheckParticleCuts(0 ,mcVtxHF)) {  // 0 stands for MC
427                                 AliDebug(2, "Skipping the particles due to cuts");
428                                 continue; 
429                         }
430                         Int_t pdgGranma = CheckOrigin(mcVtxHF, mcArray);
431                         Int_t abspdgGranma = TMath::Abs(pdgGranma);
432                         if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
433                                 AliDebug(2,Form("At Reco level, from MC info: Particle has a b-meson, or b-baryon mother (pdg code mother = %d )--> not coming from a c-quark, skipping...", pdgGranma));
434                                 continue;  // skipping particles that don't come from c quark
435                         }
436
437                         // fill the container...
438                         Double_t pt = d0tokpi->Pt();
439                         Double_t rapidity = d0tokpi->YD0();
440                         Double_t invMass=0.;
441                         Double_t cosThetaStar = 9999.;
442                         Double_t pTpi = 0.;
443                         Double_t pTK = 0.;
444                         Double_t dca = d0tokpi->GetDCA();
445                         Double_t d0pi = 0.;
446                         Double_t d0K = 0.;
447                         Double_t d0xd0 = d0tokpi->Prodd0d0();
448                         Double_t cosPointingAngle = d0tokpi->CosPointingAngle();
449                         Double_t phi = d0tokpi->Phi();
450                         Int_t pdgCode = mcVtxHF->GetPdgCode();
451                         if (pdgCode > 0){
452                                 cosThetaStar = d0tokpi->CosThetaStarD0();
453                                 pTpi = d0tokpi->PtProng(0);
454                                 pTK = d0tokpi->PtProng(1);
455                                 d0pi = d0tokpi->Getd0Prong(0);
456                                 d0K = d0tokpi->Getd0Prong(1);
457                                 invMass=d0tokpi->InvMassD0();
458                         }
459                         else {
460                                 cosThetaStar = d0tokpi->CosThetaStarD0bar();
461                                 pTpi = d0tokpi->PtProng(1);
462                                 pTK = d0tokpi->PtProng(0);
463                                 d0pi = d0tokpi->Getd0Prong(1);
464                                 d0K = d0tokpi->Getd0Prong(0);
465                                 invMass=d0tokpi->InvMassD0bar();
466                         }
467
468                         Double_t cT = d0tokpi->CtD0();
469
470                         if (!fFillFromGenerated){
471                                 // ...either with reconstructed values....
472                                 containerInput[0] = pt;
473                                 containerInput[1] = rapidity;
474                                 containerInput[2] = cosThetaStar;
475                                 containerInput[3] = pTpi;
476                                 containerInput[4] = pTK;
477                                 containerInput[5] = cT*1.E4;  // in micron
478                                 containerInput[6] = dca*1.E4;  // in micron
479                                 containerInput[7] = d0pi*1.E4;  // in micron
480                                 containerInput[8] = d0K*1.E4;  // in micron
481                                 containerInput[9] = d0xd0*1.E8;  // in micron^2
482                                 containerInput[10] = cosPointingAngle;  // in micron
483                                 containerInput[11] = phi;  
484                                 containerInputMC[12] = zPrimVertex;    // z of reconstructed of primary vertex
485                         }
486                         else {
487                                 // ... or with generated values                         
488                                 Double_t vectorMC[7] = {9999.,9999.,9999.,9999.,9999.,9999.,9999.};
489                                 if (GetGeneratedValuesFromMCParticle(mcVtxHF, mcArray, vectorMC)){
490                                         containerInput[0] = vectorMC[0];
491                                         containerInput[1] = vectorMC[1] ;
492                                         containerInput[2] = vectorMC[2] ;
493                                         containerInput[3] = vectorMC[3] ;
494                                         containerInput[4] = vectorMC[4] ;
495                                         containerInput[5] = vectorMC[5] ;  // in micron
496                                         containerInput[6] = 0.;    // dummy value, meaningless in MC, in micron
497                                         containerInput[7] = 0.;   // dummy value, meaningless in MC, in micron
498                                         containerInput[8] = 0.;   // dummy value, meaningless in MC, in micron
499                                         containerInput[9] = 100000.; // dummy value, meaningless in MC, in micron^2
500                                         containerInput[10] = 1.01;    // dummy value, meaningless in MC
501                                         containerInput[11] = vectorMC[6];   
502                                         containerInputMC[12] = zMCVertex;    // z of reconstructed of primary vertex
503                                 }
504                                 else {
505                                         AliDebug(3,"Problems in filling the container");
506                                         continue;
507                                 }
508                         }
509                         AliDebug(2, Form("Filling the container with pt = %f, rapidity = %f, cosThetaStar = %f, pTpi = %f, pTK = %f, cT = %f", containerInput[0], containerInput[1], containerInput[2], containerInput[3], containerInput[4], containerInput[5]));
510                         icountReco++;
511                         fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed) ;   
512
513                         // cut in acceptance
514                         Bool_t acceptanceProng0 = (TMath::Abs(d0tokpi->EtaProng(0))< 0.9 && d0tokpi->PtProng(0) > 0.1);
515                         Bool_t acceptanceProng1 = (TMath::Abs(d0tokpi->EtaProng(1))< 0.9 && d0tokpi->PtProng(1) > 0.1);
516                         if (acceptanceProng0 && acceptanceProng1) {
517                                 AliDebug(2,"D0 reco daughters in acceptance");
518                                 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance) ;
519                                 icountRecoAcc++; 
520                                 
521                                 if(fAcceptanceUnf){
522                                         
523                                         Double_t fill[4]; //fill response matrix
524                                         
525                                         // dimensions 0&1 : pt,eta (Rec)
526                                         
527                                         fill[0] = pt ;
528                                         fill[1] = rapidity;
529                                         
530                                         // dimensions 2&3 : pt,eta (MC)
531                                         
532                                         fill[2] = mcVtxHF->Pt();
533                                         fill[3] = mcVtxHF->Y();
534                                         
535                                         fCorrelation->Fill(fill);
536                                         
537                                 }  
538                                 
539                                 // cut on the min n. of clusters in ITS
540                                 Int_t ncls0=0,ncls1=0;
541                                 for(Int_t l=0;l<6;l++) {
542                                   if(TESTBIT(track0->GetITSClusterMap(),l)) ncls0++;
543                                   if(TESTBIT(track1->GetITSClusterMap(),l)) ncls1++;
544                                 }
545                                 AliDebug(2, Form("n clusters = %d", ncls0));
546                                 if (ncls0 >= fMinITSClusters && ncls1 >= fMinITSClusters) {
547                                         fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters) ;
548                                         icountRecoITSClusters++;   
549                                         AliDebug(2,Form("pT = %f, dca = %f, cosThetaStar = %f, pTpi = %f, pTK = %f, d0pi = %f, d0K = %f, d0xd0 = %f, cosPointingAngle = %f", pt, dca, cosThetaStar,pTpi, pTK, d0pi*1E4, d0K*1E4, d0xd0*1E8, cosPointingAngle));
550                                         
551                                         // PPR cuts 
552                                         //added mass cut for D0 analysis
553                                         Double_t cuts[7] = {9999999., 1.1, 0., 9999999., 9999999., 0.,0.027}; 
554                                         //cuts of the D0 analysis (looser) see AliAnalysisTaskSED0Mass.cxx
555                                         if (pt <= 1){
556                                           
557                                           //coutputmassD02
558                                           cuts[0] = 400; //dca (um)
559                                           cuts[1] = 0.8; //cosTstar
560                                           cuts[2] = 0.5; //pt  (GeV/c)
561                                           cuts[3] = 500; //d0  (um)
562                                           cuts[4] = -25000; //d0xd0 (um^2)
563                                           cuts[5] = 0.7; //cosTpointing
564                                           /*
565                                           //PPR
566                                           cuts[0] = 400; //dca (um)
567                                           cuts[1] = 0.8; //cosTstar
568                                           cuts[2] = 0.5; //pt  (GeV/c)
569                                           cuts[3] = 500; //d0  (um)
570                                           cuts[4] = -20000; //d0xd0 (um^2)
571                                           cuts[5] = 0.5; //cosTpointing
572                                           */
573                                         }
574                                         /* //not same cuts for pt = 1 to 2 and 1 to 3 GeV in case must match them because of poor stat
575                                         else if (pt > 1 && pt <= 2){
576                                          
577                                           //coutputmassD02
578                                           cuts[0] = 300;
579                                           cuts[1] = 0.8;
580                                           cuts[2] = 0.6;
581                                           cuts[3] = 1000;
582                                           cuts[4] = -25000;
583                                           cuts[5] = 0.7;
584                                           
585                                           //PPR
586                                           cuts[0] = 300;
587                                           cuts[1] = 0.8;
588                                           cuts[2] = 0.6;
589                                           cuts[3] = 500;
590                                           cuts[4] = -20000;
591                                           cuts[5] = 0.6;
592                                           
593                                           }
594                                         */
595                                         else if (pt > 1 && pt <= 3){
596                                           
597                                           //coutputmassD02
598                                           cuts[0] = 200;
599                                           cuts[1] = 0.8;
600                                           cuts[2] = 0.7;
601                                           cuts[3] = 1000;
602                                           cuts[4] = -25000;
603                                           cuts[5] = 0.8;
604                                           /*
605                                           //PPR
606                                           cuts[0] = 300;
607                                           cuts[1] = 0.8;
608                                           cuts[2] = 0.6;
609                                           cuts[3] = 500;
610                                           cuts[4] = -20000;
611                                           cuts[5] = 0.6;
612                                           */
613                                         }
614                                         else if (pt > 3 && pt <= 5){
615                                           
616                                           //coutputmassD02
617                                           cuts[0] = 200;
618                                           cuts[1] = 0.8;
619                                           cuts[2] = 0.7;
620                                           cuts[3] = 500;
621                                           cuts[4] = -15000;
622                                           cuts[5] = 0.8;
623                                           /*
624                                           //PPR
625                                           cuts[0] = 200;
626                                           cuts[1] = 0.8;
627                                           cuts[2] = 0.7;
628                                           cuts[3] = 500;
629                                           cuts[4] = -10000;
630                                           cuts[5] = 0.8;
631                                           */
632                                         }
633                                         else if (pt > 5){
634                                           
635                                           //coutputmassD02
636                                           cuts[0] = 200;
637                                           cuts[1] = 0.8;
638                                           cuts[2] = 0.7;
639                                           cuts[3] = 500;
640                                           cuts[4] = -15000;
641                                           cuts[5] = 0.9;
642                                           /*
643                                           //PPR
644                                           cuts[0] = 200;
645                                           cuts[1] = 0.8;
646                                           cuts[2] = 0.7;
647                                           cuts[3] = 500;
648                                           cuts[4] = -5000;
649                                           cuts[5] = 0.8;
650                                           */
651                                         }
652                                 
653                                         Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
654                                         if (TMath::Abs(invMass-mD0PDG) < cuts[6]
655                                             && dca*1E4 < cuts[0] 
656                                             && TMath::Abs(cosThetaStar) < cuts[1]  
657                                             && pTpi > cuts[2] 
658                                             && pTK > cuts[2]  
659                                             && TMath::Abs(d0pi*1E4) < cuts[3] 
660                                             && TMath::Abs(d0K*1E4) < cuts[3]  
661                                             && d0xd0*1E8 < cuts[4] 
662                                             && cosPointingAngle > cuts[5]
663                                             ){
664                                                 
665                                                 AliDebug(2,"Particle passed PPR cuts (actually cuts for D0 analysis!)");
666                                                 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPPR) ;   
667                                                 icountRecoPPR++;
668                                                 
669                                                 if(!fAcceptanceUnf){ // unfolding
670                                                         
671                                                         Double_t fill[4]; //fill response matrix
672                                                         
673                                                         // dimensions 0&1 : pt,eta (Rec)
674                                                         
675                                                         fill[0] = pt ;
676                                                         fill[1] = rapidity;
677                                                         
678                                                         // dimensions 2&3 : pt,eta (MC)
679                                                         
680                                                         fill[2] = mcVtxHF->Pt();
681                                                         fill[3] = mcVtxHF->Y();
682                                                         
683                                                         fCorrelation->Fill(fill);
684                                                         
685                                                 }
686                                         }
687                                         else{
688                                                 AliDebug(2,"Particle skipped due to PPR cuts");
689                                                 if (dca*1E4 > cuts[0]){
690                                                         AliDebug(2,"Problems with dca");
691                                                 }
692                                                 if ( TMath::Abs(cosThetaStar) > cuts[1]){
693                                                         AliDebug(2,"Problems with cosThetaStar");
694                                                 }
695                                                 if (pTpi < cuts[2]){
696                                                         AliDebug(2,"Problems with pTpi");
697                                                 }
698                                                 if (pTK < cuts[2]){
699                                                         AliDebug(2,"Problems with pTK");
700                                                 }
701                                                 if (TMath::Abs(d0pi*1E4) > cuts[3]){
702                                                         AliDebug(2,"Problems with d0pi");
703                                                 }
704                                                 if (TMath::Abs(d0K*1E4) > cuts[3]){
705                                                         AliDebug(2,"Problems with d0K");
706                                                 }
707                                                 if (d0xd0*1E8 > cuts[4]){
708                                                         AliDebug(2,"Problems with d0xd0");
709                                                 }
710                                                 if (cosPointingAngle < cuts[5]){
711                                                         AliDebug(2,"Problems with cosPointingAngle");
712                                                 }
713                                         }
714                                 }
715                         }
716                 }
717                 if(unsetvtx) d0tokpi->UnsetOwnPrimaryVtx();
718         } // end loop on D0->Kpi
719
720         AliDebug(2, Form("Found %i Reco particles that are D0!!",icountReco));
721
722         fCountReco+= icountReco;
723         fCountVertex+= icountVertex;
724         fCountRefit+= icountRefit;
725         fCountRecoAcc+= icountRecoAcc;
726         fCountRecoITSClusters+= icountRecoITSClusters;
727         fCountRecoPPR+= icountRecoPPR;
728         
729         fHistEventsProcessed->Fill(0);
730         /* PostData(0) is taken care of by AliAnalysisTaskSE */
731         PostData(1,fHistEventsProcessed) ;
732         PostData(2,fCFManager->GetParticleContainer()) ;
733         PostData(3,fCorrelation) ;
734 }
735
736
737 //___________________________________________________________________________
738 void AliCFHeavyFlavourTaskMultiVarMultiStep::Terminate(Option_t*)
739 {
740         // The Terminate() function is the last function to be called during
741         // a query. It always runs on the client, it can be used to present
742         // the results graphically or save the results to file.
743         
744         AliAnalysisTaskSE::Terminate();
745         
746         AliInfo(Form("Found %i MC particles that are D0 in MC, in %d events",fCountMC,fEvents));
747         AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, in %d events",fCountAcc,fEvents));
748         AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, and satisfy Vertex requirement in %d events",fCountVertex,fEvents));
749         AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, and satisfy ITS+TPC refit requirementin %d events",fCountRefit,fEvents));
750         AliInfo(Form("Found %i reco D0 that are decaying in K+pi, in %d events",fCountReco,fEvents));
751         AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and are in the requested acceptance, in %d events",fCountRecoAcc,fEvents));
752         AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and have at least %d clusters in ITS, in %d events",fCountRecoITSClusters,fMinITSClusters,fEvents));
753         AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and satisfy PPR cuts, in %d events",fCountRecoPPR,fEvents));
754         
755         // draw some example plots....
756         
757         AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
758         if(!cont) {
759           printf("CONTAINER NOT FOUND\n");
760           return;
761         }
762         // projecting the containers to obtain histograms
763         // first argument = variable, second argument = step
764
765         // MC-level
766         TH1D* h00 =   cont->ShowProjection(0,0) ;   // pt
767         TH1D* h10 =   cont->ShowProjection(1,0) ;   // rapidity
768         TH1D* h20 =   cont->ShowProjection(2,0) ;   // cosThetaStar
769         TH1D* h30 =   cont->ShowProjection(3,0) ;   // pTpi
770         TH1D* h40 =   cont->ShowProjection(4,0) ;   // pTK
771         TH1D* h50 =   cont->ShowProjection(5,0) ;   // cT
772         TH1D* h60 =   cont->ShowProjection(6,0) ;   // dca
773         TH1D* h70 =   cont->ShowProjection(7,0) ;   // d0pi
774         TH1D* h80 =   cont->ShowProjection(8,0) ;   // d0K
775         TH1D* h90 =   cont->ShowProjection(9,0) ;   // d0xd0
776         TH1D* h100 =   cont->ShowProjection(10,0) ;   // cosPointingAngle
777         TH1D* h110 =   cont->ShowProjection(11,0) ;   // phi
778         
779         // MC-Acceptance level
780         TH1D* h01 =   cont->ShowProjection(0,1) ;   // pt
781         TH1D* h11 =   cont->ShowProjection(1,1) ;   // rapidity
782         TH1D* h21 =   cont->ShowProjection(2,1) ;   // cosThetaStar
783         TH1D* h31 =   cont->ShowProjection(3,1) ;   // pTpi
784         TH1D* h41 =   cont->ShowProjection(4,1) ;   // pTK
785         TH1D* h51 =   cont->ShowProjection(5,1) ;   // cT
786         TH1D* h61 =   cont->ShowProjection(6,1) ;   // dca
787         TH1D* h71 =   cont->ShowProjection(7,1) ;   // d0pi
788         TH1D* h81 =   cont->ShowProjection(8,1) ;   // d0K
789         TH1D* h91 =   cont->ShowProjection(9,1) ;   // d0xd0
790         TH1D* h101 =   cont->ShowProjection(10,1) ;   // cosPointingAngle
791         TH1D* h111 =   cont->ShowProjection(11,1) ;   // phi
792
793         // Reco-level
794         TH1D* h04 =   cont->ShowProjection(0,4) ;   // pt
795         TH1D* h14 =   cont->ShowProjection(1,4) ;   // rapidity
796         TH1D* h24 =   cont->ShowProjection(2,4) ;   // cosThetaStar
797         TH1D* h34 =   cont->ShowProjection(3,4) ;   // pTpi
798         TH1D* h44 =   cont->ShowProjection(4,4) ;   // pTK
799         TH1D* h54 =   cont->ShowProjection(5,4) ;   // cT
800         TH1D* h64 =   cont->ShowProjection(6,4) ;   // dca
801         TH1D* h74 =   cont->ShowProjection(7,4) ;   // d0pi
802         TH1D* h84 =   cont->ShowProjection(8,4) ;   // d0K
803         TH1D* h94 =   cont->ShowProjection(9,4) ;   // d0xd0
804         TH1D* h104 =   cont->ShowProjection(10,4) ;   // cosPointingAngle
805         TH1D* h114 =   cont->ShowProjection(11,4) ;   // phi
806         
807         h00->SetTitle("pT_D0 (GeV/c)");
808         h10->SetTitle("rapidity");
809         h20->SetTitle("cosThetaStar");
810         h30->SetTitle("pT_pi (GeV/c)");
811         h40->SetTitle("pT_K (Gev/c)");
812         h50->SetTitle("cT (#mum)");
813         h60->SetTitle("dca (#mum)");
814         h70->SetTitle("d0_pi (#mum)");
815         h80->SetTitle("d0_K (#mum)");
816         h90->SetTitle("d0xd0 (#mum^2)");
817         h100->SetTitle("cosPointingAngle");
818         h100->SetTitle("phi (rad)");
819
820         h00->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
821         h10->GetXaxis()->SetTitle("rapidity");
822         h20->GetXaxis()->SetTitle("cosThetaStar");
823         h30->GetXaxis()->SetTitle("pT_pi (GeV/c)");
824         h40->GetXaxis()->SetTitle("pT_K (Gev/c)");
825         h50->GetXaxis()->SetTitle("cT (#mum)");
826         h60->GetXaxis()->SetTitle("dca (#mum)");
827         h70->GetXaxis()->SetTitle("d0_pi (#mum)");
828         h80->GetXaxis()->SetTitle("d0_K (#mum)");
829         h90->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
830         h100->GetXaxis()->SetTitle("cosPointingAngle");
831         h110->GetXaxis()->SetTitle("phi (rad)");
832
833         h01->SetTitle("pT_D0 (GeV/c)");
834         h11->SetTitle("rapidity");
835         h21->SetTitle("cosThetaStar");
836         h31->SetTitle("pT_pi (GeV/c)");
837         h41->SetTitle("pT_K (Gev/c)");
838         h51->SetTitle("cT (#mum)");
839         h61->SetTitle("dca (#mum)");
840         h71->SetTitle("d0_pi (#mum)");
841         h81->SetTitle("d0_K (#mum)");
842         h91->SetTitle("d0xd0 (#mum^2)");
843         h101->SetTitle("cosPointingAngle");
844         h111->GetXaxis()->SetTitle("phi (rad)");
845
846         h01->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
847         h11->GetXaxis()->SetTitle("rapidity");
848         h21->GetXaxis()->SetTitle("cosThetaStar");
849         h31->GetXaxis()->SetTitle("pT_pi (GeV/c)");
850         h41->GetXaxis()->SetTitle("pT_K (Gev/c)");
851         h51->GetXaxis()->SetTitle("cT (#mum)");
852         h61->GetXaxis()->SetTitle("dca (#mum)");
853         h71->GetXaxis()->SetTitle("d0_pi (#mum)");
854         h81->GetXaxis()->SetTitle("d0_K (#mum)");
855         h91->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
856         h101->GetXaxis()->SetTitle("cosPointingAngle");
857         h111->GetXaxis()->SetTitle("phi (rad)");
858
859         h04->SetTitle("pT_D0 (GeV/c)");
860         h14->SetTitle("rapidity");
861         h24->SetTitle("cosThetaStar");
862         h34->SetTitle("pT_pi (GeV/c)");
863         h44->SetTitle("pT_K (Gev/c)");
864         h54->SetTitle("cT (#mum)");
865         h64->SetTitle("dca (#mum)");
866         h74->SetTitle("d0_pi (#mum)");
867         h84->SetTitle("d0_K (#mum)");
868         h94->SetTitle("d0xd0 (#mum^2)");
869         h104->SetTitle("cosPointingAngle");
870         h114->GetXaxis()->SetTitle("phi (rad)");
871
872         h04->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
873         h14->GetXaxis()->SetTitle("rapidity");
874         h24->GetXaxis()->SetTitle("cosThetaStar");
875         h34->GetXaxis()->SetTitle("pT_pi (GeV/c)");
876         h44->GetXaxis()->SetTitle("pT_K (Gev/c)");
877         h54->GetXaxis()->SetTitle("cT (#mum)");
878         h64->GetXaxis()->SetTitle("dca (#mum)");
879         h74->GetXaxis()->SetTitle("d0_pi (#mum)");
880         h84->GetXaxis()->SetTitle("d0_K (#mum)");
881         h94->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
882         h104->GetXaxis()->SetTitle("cosPointingAngle");
883         h114->GetXaxis()->SetTitle("phi (rad)");
884
885         Double_t max0 = h00->GetMaximum();
886         Double_t max1 = h10->GetMaximum();
887         Double_t max2 = h20->GetMaximum();
888         Double_t max3 = h30->GetMaximum();
889         Double_t max4 = h40->GetMaximum();
890         Double_t max5 = h50->GetMaximum();
891         Double_t max6 = h60->GetMaximum();
892         Double_t max7 = h70->GetMaximum();
893         Double_t max8 = h80->GetMaximum();
894         Double_t max9 = h90->GetMaximum();
895         Double_t max10 = h100->GetMaximum();
896         Double_t max11 = h110->GetMaximum();
897         
898         h00->GetYaxis()->SetRangeUser(0,max0*1.2);
899         h10->GetYaxis()->SetRangeUser(0,max1*1.2);
900         h20->GetYaxis()->SetRangeUser(0,max2*1.2);
901         h30->GetYaxis()->SetRangeUser(0,max3*1.2);
902         h40->GetYaxis()->SetRangeUser(0,max4*1.2);
903         h50->GetYaxis()->SetRangeUser(0,max5*1.2);
904         h60->GetYaxis()->SetRangeUser(0,max6*1.2);
905         h70->GetYaxis()->SetRangeUser(0,max7*1.2);
906         h80->GetYaxis()->SetRangeUser(0,max8*1.2);
907         h90->GetYaxis()->SetRangeUser(0,max9*1.2);
908         h100->GetYaxis()->SetRangeUser(0,max10*1.2);
909         h110->GetYaxis()->SetRangeUser(0,max11*1.2);
910         
911         h01->GetYaxis()->SetRangeUser(0,max0*1.2);
912         h11->GetYaxis()->SetRangeUser(0,max1*1.2);
913         h21->GetYaxis()->SetRangeUser(0,max2*1.2);
914         h31->GetYaxis()->SetRangeUser(0,max3*1.2);
915         h41->GetYaxis()->SetRangeUser(0,max4*1.2);
916         h51->GetYaxis()->SetRangeUser(0,max5*1.2);
917         h61->GetYaxis()->SetRangeUser(0,max6*1.2);
918         h71->GetYaxis()->SetRangeUser(0,max7*1.2);
919         h81->GetYaxis()->SetRangeUser(0,max8*1.2);
920         h91->GetYaxis()->SetRangeUser(0,max9*1.2);
921         h101->GetYaxis()->SetRangeUser(0,max10*1.2);
922         h111->GetYaxis()->SetRangeUser(0,max11*1.2);
923         
924         h00->SetMarkerStyle(20);
925         h10->SetMarkerStyle(24);
926         h20->SetMarkerStyle(21);
927         h30->SetMarkerStyle(25);
928         h40->SetMarkerStyle(27);
929         h50->SetMarkerStyle(28);
930         h60->SetMarkerStyle(20);
931         h70->SetMarkerStyle(24);
932         h80->SetMarkerStyle(21);
933         h90->SetMarkerStyle(25);
934         h100->SetMarkerStyle(27);
935         h110->SetMarkerStyle(28);
936
937         h00->SetMarkerColor(2);
938         h10->SetMarkerColor(2);
939         h20->SetMarkerColor(2);
940         h30->SetMarkerColor(2);
941         h40->SetMarkerColor(2);
942         h50->SetMarkerColor(2);
943         h60->SetMarkerColor(2);
944         h70->SetMarkerColor(2);
945         h80->SetMarkerColor(2);
946         h90->SetMarkerColor(2);
947         h100->SetMarkerColor(2);
948         h110->SetMarkerColor(2);
949
950         h01->SetMarkerStyle(20) ;
951         h11->SetMarkerStyle(24) ;
952         h21->SetMarkerStyle(21) ;
953         h31->SetMarkerStyle(25) ;
954         h41->SetMarkerStyle(27) ;
955         h51->SetMarkerStyle(28) ;
956         h61->SetMarkerStyle(20);
957         h71->SetMarkerStyle(24);
958         h81->SetMarkerStyle(21);
959         h91->SetMarkerStyle(25);
960         h101->SetMarkerStyle(27);
961         h111->SetMarkerStyle(28);
962
963         h01->SetMarkerColor(8);
964         h11->SetMarkerColor(8);
965         h21->SetMarkerColor(8);
966         h31->SetMarkerColor(8);
967         h41->SetMarkerColor(8);
968         h51->SetMarkerColor(8);
969         h61->SetMarkerColor(8);
970         h71->SetMarkerColor(8);
971         h81->SetMarkerColor(8);
972         h91->SetMarkerColor(8);
973         h101->SetMarkerColor(8);
974         h111->SetMarkerColor(8);
975
976         h04->SetMarkerStyle(20) ;
977         h14->SetMarkerStyle(24) ;
978         h24->SetMarkerStyle(21) ;
979         h34->SetMarkerStyle(25) ;
980         h44->SetMarkerStyle(27) ;
981         h54->SetMarkerStyle(28) ;
982         h64->SetMarkerStyle(20);
983         h74->SetMarkerStyle(24);
984         h84->SetMarkerStyle(21);
985         h94->SetMarkerStyle(25);
986         h104->SetMarkerStyle(27);
987         h114->SetMarkerStyle(28);
988
989         h04->SetMarkerColor(4);
990         h14->SetMarkerColor(4);
991         h24->SetMarkerColor(4);
992         h34->SetMarkerColor(4);
993         h44->SetMarkerColor(4);
994         h54->SetMarkerColor(4);
995         h64->SetMarkerColor(4);
996         h74->SetMarkerColor(4);
997         h84->SetMarkerColor(4);
998         h94->SetMarkerColor(4);
999         h104->SetMarkerColor(4);
1000         h114->SetMarkerColor(4);
1001         
1002         gStyle->SetCanvasColor(0);
1003         gStyle->SetFrameFillColor(0);
1004         gStyle->SetTitleFillColor(0);
1005         gStyle->SetStatColor(0);
1006
1007         // drawing in 2 separate canvas for a matter of clearity
1008         TCanvas * c1 =new TCanvas("c1","pT, rapidiy, cosThetaStar",1100,1600);
1009         c1->Divide(3,3);
1010         
1011         c1->cd(1);
1012         h00->Draw("p");
1013         c1->cd(1);
1014         c1->cd(2);
1015         h01->Draw("p");
1016         c1->cd(2);
1017         c1->cd(3);
1018         h04->Draw("p");
1019         c1->cd(3);
1020         c1->cd(4);
1021         h10->Draw("p");
1022         c1->cd(4);
1023         c1->cd(5);
1024         h11->Draw("p");
1025         c1->cd(5);
1026         c1->cd(6);
1027         h14->Draw("p");
1028         c1->cd(6);
1029         c1->cd(7);
1030         h20->Draw("p");
1031         c1->cd(7);
1032         c1->cd(8);
1033         h21->Draw("p");
1034         c1->cd(8);
1035         c1->cd(9);
1036         h24->Draw("p");
1037         c1->cd(9);
1038         c1->cd();
1039
1040         TCanvas * c2 =new TCanvas("c2","pTpi, pTK, cT",1100,1600);
1041         c2->Divide(3,3);
1042         c2->cd(1);
1043         h30->Draw("p");
1044         c2->cd(1);
1045         c2->cd(2);
1046         h31->Draw("p");
1047         c2->cd(2);
1048         c2->cd(3);
1049         h34->Draw("p");
1050         c2->cd(3);
1051         c2->cd(4);
1052         h40->Draw("p");
1053         c2->cd(4);
1054         c2->cd(5);
1055         h41->Draw("p");
1056         c2->cd(5);
1057         c2->cd(6);
1058         h44->Draw("p");
1059         c2->cd(6);
1060         c2->cd(7);
1061         h50->Draw("p");
1062         c2->cd(7);
1063         c2->cd(8);
1064         h51->Draw("p");
1065         c2->cd(8);
1066         c2->cd(9);
1067         h54->Draw("p");
1068         c2->cd(9);
1069         c2->cd();
1070
1071         TCanvas * c3 =new TCanvas("c3","dca, d0pi, d0K",1100,1600);
1072         c3->Divide(3,3);
1073         c3->cd(1);
1074         h60->Draw("p");
1075         c3->cd(1);
1076         c3->cd(2);
1077         h61->Draw("p");
1078         c3->cd(2);
1079         c3->cd(3);
1080         h64->Draw("p");
1081         c3->cd(3);
1082         c3->cd(4);
1083         h70->Draw("p");
1084         c3->cd(4);
1085         c3->cd(5);
1086         h71->Draw("p");
1087         c3->cd(5);
1088         c3->cd(6);
1089         h74->Draw("p");
1090         c3->cd(6);
1091         c3->cd(7);
1092         h80->Draw("p");
1093         c3->cd(7);
1094         c3->cd(8);
1095         h81->Draw("p");
1096         c3->cd(8);
1097         c3->cd(9);
1098         h84->Draw("p");
1099         c3->cd(9);
1100         c3->cd();
1101
1102         TCanvas * c4 =new TCanvas("c4","d0xd0, cosPointingAngle, phi",1100,1600);
1103         c4->Divide(3,3);
1104         c4->cd(1);
1105         h90->Draw("p");
1106         c4->cd(1);
1107         c4->cd(2);
1108         h91->Draw("p");
1109         c4->cd(2);
1110         c4->cd(3);
1111         h94->Draw("p");
1112         c4->cd(3);
1113         c4->cd(4);
1114         h100->Draw("p");
1115         c4->cd(4);
1116         c4->cd(5);
1117         h101->Draw("p");
1118         c4->cd(5);
1119         c4->cd(6);
1120         h104->Draw("p");
1121         c4->cd(6);
1122         c4->cd(7);
1123         h110->Draw("p");
1124         c4->cd(7);
1125         c4->cd(8);
1126         h111->Draw("p");
1127         c4->cd(8);
1128         c4->cd(9);
1129         h114->Draw("p");
1130         c4->cd(9);
1131         c4->cd();
1132
1133         THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
1134
1135         TH2D* corr1 =hcorr->Projection(0,2);
1136         TH2D* corr2 = hcorr->Projection(1,3);
1137
1138         TCanvas * c7 =new TCanvas("c7","",800,400);
1139         c7->Divide(2,1);
1140         c7->cd(1);
1141         corr1->Draw("text");
1142         c7->cd(2);
1143         corr2->Draw("text");
1144       
1145
1146         TFile* file_projection = new TFile("CFtaskHFprojection.root","RECREATE");
1147
1148         corr1->Write();
1149         corr2->Write();
1150         h00->Write("pT_D0_step0");
1151         h10->Write("rapidity_step0");
1152         h20->Write("cosThetaStar_step0");
1153         h30->Write("pT_pi_step0");
1154         h40->Write("pT_K_step0");
1155         h50->Write("cT_step0");
1156         h60->Write("dca_step0");
1157         h70->Write("d0_pi_step0");
1158         h80->Write("d0_K_step0");
1159         h90->Write("d0xd0_step0");
1160         h100->Write("cosPointingAngle_step0");
1161         h110->Write("phi_step0");
1162
1163         h01->Write("pT_D0_step1");
1164         h11->Write("rapidity_step1");
1165         h21->Write("cosThetaStar_step1");
1166         h31->Write("pT_pi_step1");
1167         h41->Write("pT_K_step1");
1168         h51->Write("cT_step1");
1169         h61->Write("dca_step1");
1170         h71->Write("d0_pi_step1");
1171         h81->Write("d0_K_step1");
1172         h91->Write("d0xd0_step1");
1173         h101->Write("cosPointingAngle_step1");
1174         h111->Write("phi_step1");
1175
1176         h04->Write("pT_D0_step2");
1177         h14->Write("rapidity_step2");
1178         h24->Write("cosThetaStar_step2");
1179         h34->Write("pT_pi_step2");
1180         h44->Write("pT_K_step2");
1181         h54->Write("cT_step2");
1182         h64->Write("dca_step2");
1183         h74->Write("d0_pi_step2");
1184         h80->Write("d0_K_step2");
1185         h94->Write("d0xd0_step2");
1186         h104->Write("cosPointingAngle_step2");
1187         h114->Write("phi_step2");
1188
1189         file_projection->Close();
1190
1191     
1192
1193         /*
1194         c1->SaveAs("Plots/pT_rapidity_cosThetaStar.eps");
1195         c2->SaveAs("Plots/pTpi_pTK_cT.eps");
1196         c3->SaveAs("Plots/dca_d0pi_d0TK.eps");
1197         c4->SaveAs("Plots/d0xd0_cosPointingAngle.eps");
1198
1199         c1->SaveAs("Plots/pT_rapidity_cosThetaStar.gif");
1200         c2->SaveAs("Plots/pTpi_pTK_cT.gif");
1201         c3->SaveAs("Plots/dca_d0pi_d0TK.gif");
1202         c4->SaveAs("Plots/d0xd0_cosPointingAngle.gif");
1203         */      
1204 }
1205
1206 //___________________________________________________________________________
1207 void AliCFHeavyFlavourTaskMultiVarMultiStep::UserCreateOutputObjects() {
1208         //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1209         //TO BE SET BEFORE THE EXECUTION OF THE TASK
1210         //
1211         Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
1212         
1213         //slot #1
1214         OpenFile(1);
1215         fHistEventsProcessed = new TH1I("CFHFchist0","",1,0,1) ;
1216 }
1217
1218 //___________________________________________________________________________
1219 Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CosThetaStar(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
1220
1221         //
1222         // to calculate cos(ThetaStar) for generated particle
1223         // using the K, since mcPartDaughter0 and mcPartDaughter1 always correspond to K and pi respectively 
1224         // (see where the function is called)
1225         //
1226
1227         Int_t pdgvtx = mcPart->GetPdgCode();
1228         /*      if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
1229                 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1230                 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1231                 AliInfo(Form("D0, with pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
1232                 AliDebug(2,"This is a D0");
1233                 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
1234                 mcPartDaughter0 = mcPartDaughter1;
1235                 mcPartDaughter1 = mcPartdummy;
1236         } 
1237         else{
1238                 AliInfo("D0bar");
1239         }
1240         */
1241         Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1242         Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1243         if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
1244                 AliDebug(2,"D0");
1245         }
1246         else{
1247                 AliDebug(2,"D0bar");
1248         }
1249         if (pdgprong0 == 211){
1250                 AliDebug(2,Form("pdgprong0 = %d, pdgprong1 = %d, switching...",pdgprong0,pdgprong1));
1251                 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
1252                 mcPartDaughter0 = mcPartDaughter1;
1253                 mcPartDaughter1 = mcPartdummy;
1254                 pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
1255                 pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
1256         } 
1257
1258         AliDebug(2,Form("After checking, pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
1259         Double_t massvtx = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdgvtx))->Mass();
1260         Double_t massp[2];
1261         massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
1262         massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
1263
1264         Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)-4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx);
1265
1266         Double_t px = mcPartDaughter0->Px()+mcPartDaughter1->Px();
1267         Double_t py = mcPartDaughter0->Py()+mcPartDaughter1->Py();
1268         Double_t pz = mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
1269         Double_t p = TMath::Sqrt(px*px+py*py+pz*pz);
1270         Double_t e =  TMath::Sqrt(massvtx*massvtx+p*p);
1271         Double_t beta = p/e;
1272         Double_t gamma = e/massvtx;
1273         TVector3 mom(mcPartDaughter0->Px(),mcPartDaughter0->Py(),mcPartDaughter0->Pz());
1274         TVector3 momTot(mcPartDaughter0->Px()+mcPartDaughter1->Px(),mcPartDaughter0->Py()+mcPartDaughter1->Py(),mcPartDaughter0->Pz()+mcPartDaughter1->Pz());
1275
1276         Double_t qlprong = mom.Dot(momTot)/momTot.Mag();  // analog to AliAODRecoDecay::QlProng(0)
1277         
1278         AliDebug(2,Form("pStar = %f, beta = %f, gamma = %f, qlprong = %f, massp[0] = %f", pStar, beta, gamma, qlprong, massp[0]));
1279         Double_t cts = (qlprong/gamma-beta*TMath::Sqrt(pStar*pStar+massp[0]*massp[0]))/pStar;
1280         AliDebug(2,Form("cts = %f", cts));
1281         return cts;
1282 }
1283 //___________________________________________________________________________
1284 Double_t AliCFHeavyFlavourTaskMultiVarMultiStep::CT(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
1285
1286         //
1287         // to calculate cT for generated particle
1288         //
1289
1290         Double_t xmcPart[3] = {0,0,0};
1291         Double_t xdaughter0[3] = {0,0,0};
1292         Double_t xdaughter1[3] = {0,0,0};
1293         mcPart->XvYvZv(xmcPart);  // cm
1294         mcPartDaughter0->XvYvZv(xdaughter0);  // cm
1295         mcPartDaughter1->XvYvZv(xdaughter1);  //cm
1296         Double_t prodVtxD0 = TMath::Sqrt(xmcPart[0]*xmcPart[0]+
1297                                          xmcPart[1]*xmcPart[1]+
1298                                          xmcPart[2]*xmcPart[2]);
1299         Double_t prodVtxDaughter0 = TMath::Sqrt(xdaughter0[0]*xdaughter0[0]+
1300                                                 xdaughter0[1]*xdaughter0[1]+
1301                                                 xdaughter0[2]*xdaughter0[2]);
1302         Double_t prodVtxDaughter1 = TMath::Sqrt(xdaughter1[0]*xdaughter1[0]+
1303                                                 xdaughter1[1]*xdaughter1[1]+
1304                                                 xdaughter1[2]*xdaughter1[2]);
1305  
1306         AliDebug(2, Form("D0:        x = %f, y = %f, z = %f, production vertex distance = %f (cm), %f (micron)", xmcPart[0], xmcPart[1], xmcPart[2], prodVtxD0, prodVtxD0*1.E4));
1307         AliDebug(2, Form("Daughter0: x = %f, y = %f, z = %f, production vertex distance = %f (cm) %f (micron)", xdaughter0[0], xdaughter0[1], xdaughter0[2], prodVtxDaughter0, prodVtxDaughter0*1E4));
1308         AliDebug(2, Form("Daughter1: x = %f, y = %f, z = %f, production vertex distance = %f (cm) %f (micron)", xdaughter1[0], xdaughter1[1], xdaughter1[2], prodVtxDaughter1, prodVtxDaughter1*1.E4));
1309
1310         Double_t cT0 = TMath::Sqrt((xdaughter0[0]-xmcPart[0])*(xdaughter0[0]-xmcPart[0])+
1311                                      (xdaughter0[1]-xmcPart[1])*(xdaughter0[1]-xmcPart[1])+
1312                                      (xdaughter0[2]-xmcPart[2])*(xdaughter0[2]-xmcPart[2]));
1313
1314         Double_t cT1 = TMath::Sqrt((xdaughter1[0]-xmcPart[0])*(xdaughter1[0]-xmcPart[0])+
1315                                      (xdaughter1[1]-xmcPart[1])*(xdaughter1[1]-xmcPart[1])+
1316                                      (xdaughter1[2]-xmcPart[2])*(xdaughter1[2]-xmcPart[2]));
1317
1318         if ((cT0 - cT1)>1E-5) {
1319                 AliWarning(Form("cT from daughter 0 (%f) different from cT from daughter 1 (%f)! Using cT from daughter 0, but PLEASE, CHECK....",cT0,cT1));
1320         }
1321
1322         // calculating cT from cT0
1323
1324         Double_t mass = TDatabasePDG::Instance()->GetParticle(mcPart->GetPdgCode())->Mass(); // mcPart->GetPdgCode() should return +/- 421 for the D0/D0bar
1325         Double_t p = mcPart-> P();
1326         Double_t cT = cT0*mass/p;
1327         AliDebug(2, Form("cT from daughter 0 = %f (micron)", cT0*1E4)); 
1328         AliDebug(2, Form("cT from daughter 1 = %f (micron)", cT1*1E4)); 
1329         AliDebug(2, Form("cT (with mass = %f and p = %f) = %f (micron)", mass, p, cT*1E4));
1330         return cT;
1331 }
1332 //________________________________________________________________________________
1333 Bool_t AliCFHeavyFlavourTaskMultiVarMultiStep::GetGeneratedValuesFromMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC) const {
1334
1335         // 
1336         // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
1337         //
1338
1339         Bool_t isOk = kFALSE;
1340
1341         // check whether the D0 decays in pi+K
1342         // to be added!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1343         // could use a cut in the CF, but not clear how to define a TDecayChannel
1344         // implemented for the time being as a cut on the number of daughters - see later when 
1345         // getting the daughters
1346
1347         // getting the daughters
1348         Int_t daughter0 = mcPart->GetDaughter(0);
1349         Int_t daughter1 = mcPart->GetDaughter(1);
1350         AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
1351         if (daughter0 == 0 || daughter1 == 0) {
1352                 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
1353                 return isOk;  
1354         }
1355         if (TMath::Abs(daughter1 - daughter0) != 1) {
1356                 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
1357                 return isOk;  
1358         }
1359         AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
1360         AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
1361         if (!mcPartDaughter0 || !mcPartDaughter1) {
1362                 AliWarning("At least one Daughter Particle not found in tree, skipping"); 
1363                 return isOk;  
1364         }
1365         if (!(TMath::Abs(mcPartDaughter0->GetPdgCode())==321 &&
1366               TMath::Abs(mcPartDaughter1->GetPdgCode())==211) && 
1367             !(TMath::Abs(mcPartDaughter0->GetPdgCode())==211 &&
1368               TMath::Abs(mcPartDaughter1->GetPdgCode())==321)) {
1369           AliDebug(2, "The D0 MC doesn't come from a Kpi decay, skipping!!");
1370           return isOk;  
1371         }
1372
1373         Double_t vtx1[3] = {0,0,0};   // primary vertex         
1374         Double_t vtx2daughter0[3] = {0,0,0};   // secondary vertex from daughter 0
1375         Double_t vtx2daughter1[3] = {0,0,0};   // secondary vertex from daughter 1
1376         mcPart->XvYvZv(vtx1);  // cm
1377         // getting vertex from daughters
1378         mcPartDaughter0->XvYvZv(vtx2daughter0);  // cm
1379         mcPartDaughter1->XvYvZv(vtx2daughter1);  //cm
1380         if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
1381                 AliError("Daughters have different secondary vertex, skipping the track");
1382                 return isOk;
1383         }
1384         Int_t nprongs = 2;
1385         Short_t charge = 0;
1386         // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
1387         AliAODMCParticle* positiveDaugh = mcPartDaughter0;
1388         AliAODMCParticle* negativeDaugh = mcPartDaughter1;
1389         if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
1390                 // inverting in case the positive daughter is the second one
1391                 positiveDaugh = mcPartDaughter1;
1392                 negativeDaugh = mcPartDaughter0;
1393         }
1394         // getting the momentum from the daughters
1395         Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};            
1396         Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};            
1397         Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
1398
1399         Double_t d0[2] = {0.,0.};               
1400
1401         AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
1402
1403         Double_t cosThetaStar = 0.;
1404         Double_t cosThetaStarD0 = 0.;
1405         Double_t cosThetaStarD0bar = 0.;
1406         cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
1407         cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
1408         if (mcPart->GetPdgCode() == 421){  // D0
1409                 AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1410                 cosThetaStar = cosThetaStarD0;
1411         }
1412         else if (mcPart->GetPdgCode() == -421){  // D0bar{
1413                 AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
1414                 cosThetaStar = cosThetaStarD0bar;
1415         }
1416         else{
1417                 AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
1418                 return vectorMC;
1419         }
1420         if (cosThetaStar < -1 || cosThetaStar > 1) {
1421                 AliWarning("Invalid value for cosine Theta star, returning");
1422                 return isOk;
1423         }
1424
1425         // calculate cos(Theta*) according to the method implemented herein
1426
1427         Double_t cts = 9999.;
1428         cts = CosThetaStar(mcPart, mcPartDaughter0, mcPartDaughter1);
1429         if (cts < -1 || cts > 1) {
1430                 AliWarning("Invalid value for cosine Theta star from AliCFHeavyFlavourTaskMultiVarMultiStep method");
1431         }
1432         if (TMath::Abs(cts - cosThetaStar)>0.001) {
1433                 AliError(Form("cosThetaStar automatically calculated different from that manually calculated!!! cosThetaStar = %f, cosThetaStar = %f", cosThetaStar,cts));
1434         }
1435         
1436         Double_t cT = decay->Ct(421);
1437
1438         // calculate cT from the method implemented herein
1439         Double_t cT1 = 0.;
1440         cT1 = CT(mcPart, mcPartDaughter0, mcPartDaughter1);
1441
1442         if (TMath::Abs(cT1 - cT)>0.001) {
1443                 AliError(Form("cT automatically calculated different from that manually calculated!!! cT = %f, cT1 = %f",cT,cT1));
1444         }
1445         
1446         // get the pT of the daughters
1447         
1448         Double_t pTpi = 0.;
1449         Double_t pTK = 0.;
1450         
1451         if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
1452                 pTpi = mcPartDaughter0->Pt();
1453                 pTK = mcPartDaughter1->Pt();
1454         }
1455         else {
1456                 pTpi = mcPartDaughter1->Pt();
1457                 pTK = mcPartDaughter0->Pt();
1458         }
1459
1460         vectorMC[0] = mcPart->Pt();
1461         vectorMC[1] = mcPart->Y() ;
1462         vectorMC[2] = cosThetaStar ;
1463         vectorMC[3] = pTpi ;
1464         vectorMC[4] = pTK ;
1465         vectorMC[5] = cT*1.E4 ;  // in micron
1466         vectorMC[6] = mcPart->Phi() ;  
1467         isOk = kTRUE;
1468         return isOk;
1469 }
1470 //_________________________________________________________________________________________________
1471 Int_t AliCFHeavyFlavourTaskMultiVarMultiStep::CheckOrigin(AliAODMCParticle* mcPart, TClonesArray* mcArray)const{
1472
1473         //
1474         // checking whether the very mother of the D0 is a charm or a bottom quark
1475         //
1476
1477         Int_t pdgGranma = 0;
1478         Int_t mother = 0;
1479         mother = mcPart->GetMother();
1480         Int_t istep = 0;
1481         while (mother >0 ){
1482                 istep++;
1483                 AliDebug(2,Form("mother at step %d = %d", istep, mother));
1484                 AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(mcArray->At(mother));
1485                 pdgGranma = mcGranma->GetPdgCode();
1486                 AliDebug(2,Form("Pdg mother at step %d = %d", istep, pdgGranma));
1487                 Int_t abspdgGranma = TMath::Abs(pdgGranma);
1488                 if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)) {
1489                         break;
1490                 }
1491                 mother = mcGranma->GetMother();
1492         }
1493         return pdgGranma;
1494 }
1495