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