PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVarMultiStep.cxx

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