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