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