PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVarMultiStep.cxx

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