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