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