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