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