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