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