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