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