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