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