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