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