PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVar.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 //
  19 // Example of task running on AliEn
  20 // which provides standard way of calculating acceptance and efficiency
  21 // between different steps of the procedure.
  22 // The ouptut of the task is a AliCFContainer from which the efficiencies
  23 // can be calculated
  24 //-----------------------------------------------------------------------
  25 // Author : C. Zampolli, CERN, on the basis of R. Vernet's example
  26 //-----------------------------------------------------------------------
  27
  28 #include <TCanvas.h>
  29 #include <TParticle.h>
  30 #include <TDatabasePDG.h>
  31 #include <TH1I.h>
  32 #include <TStyle.h>
  33
  34 #include "AliCFHeavyFlavourTaskMultiVar.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 AliCFHeavyFlavourTaskMultiVar::AliCFHeavyFlavourTaskMultiVar() :
  48         AliAnalysisTaskSE(),
  49         fPDG(0),
  50         fCFManager(0x0),
  51         fHistEventsProcessed(0x0),
  52         fCountMC(0),
  53         fEvents(0),
  54         fFillFromGenerated(kFALSE)
  55 {
  56         //
  57         //Default ctor
  58         //
  59 }
  60 //___________________________________________________________________________
  61 AliCFHeavyFlavourTaskMultiVar::AliCFHeavyFlavourTaskMultiVar(const Char_t* name) :
  62         AliAnalysisTaskSE(name),
  63         fPDG(0),
  64         fCFManager(0x0),
  65         fHistEventsProcessed(0x0),
  66         fCountMC(0),
  67         fEvents(0),
  68         fFillFromGenerated(kFALSE)
  69 {
  70         //
  71         // Constructor. Initialization of Inputs and Outputs
  72         //
  73         Info("AliCFHeavyFlavourTaskMultiVar","Calling Constructor");
  74         /*
  75           DefineInput(0) and DefineOutput(0)
  76           are taken care of by AliAnalysisTaskSE constructor
  77         */
  78         DefineOutput(1,TH1I::Class());
  79         DefineOutput(2,AliCFContainer::Class());
  80 }
  81
  82 //___________________________________________________________________________
  83 AliCFHeavyFlavourTaskMultiVar& AliCFHeavyFlavourTaskMultiVar::operator=(const AliCFHeavyFlavourTaskMultiVar& c)
  84 {
  85         //
  86         // Assignment operator
  87         //
  88         if (this!=&c) {
  89                 AliAnalysisTaskSE::operator=(c) ;
  90                 fPDG      = c.fPDG;
  91                 fCFManager  = c.fCFManager;
  92                 fHistEventsProcessed = c.fHistEventsProcessed;
  93         }
  94         return *this;
  95 }
  96
  97 //___________________________________________________________________________
  98 AliCFHeavyFlavourTaskMultiVar::AliCFHeavyFlavourTaskMultiVar(const AliCFHeavyFlavourTaskMultiVar& c) :
  99         AliAnalysisTaskSE(c),
 100         fPDG(c.fPDG),
 101         fCFManager(c.fCFManager),
 102         fHistEventsProcessed(c.fHistEventsProcessed),
 103         fCountMC(c.fCountMC),
 104         fEvents(c.fEvents),
 105         fFillFromGenerated(c.fFillFromGenerated)
 106 {
 107         //
 108         // Copy Constructor
 109         //
 110 }
 111
 112 //___________________________________________________________________________
 113 AliCFHeavyFlavourTaskMultiVar::~AliCFHeavyFlavourTaskMultiVar() {
 114         //
 115         //destructor
 116         //
 117         Info("~AliCFHeavyFlavourTaskMultiVar","Calling Destructor");
 118         if (fCFManager)           delete fCFManager ;
 119         if (fHistEventsProcessed) delete fHistEventsProcessed ;
 120 }
 121
 122 //_________________________________________________
 123 void AliCFHeavyFlavourTaskMultiVar::UserExec(Option_t *)
 124 {
 125         //
 126         // Main loop function
 127         //
 128
 129         if (!fInputEvent) {
 130                 Error("UserExec","NO EVENT FOUND!");
 131                 return;
 132         }
 133
 134         fEvents++;
 135         if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents));
 136         AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 137         fCFManager->SetEventInfo(aodEvent);
 138
 139         // MC-event selection
 140         Double_t containerInput[6] ;
 141
 142         //loop on the MC event
 143
 144         TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 145         if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
 146         Int_t icountMC = 0;
 147
 148         for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
 149                 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
 150                 if (!mcPart) {
 151                         AliWarning("Particle not found in tree, skipping");
 152                         continue;
 153                 }
 154
 155                 // check the MC-level cuts
 156                 if (!fCFManager->CheckParticleCuts(0, mcPart)) continue;  // 0 stands for MC level
 157
 158                 // fill the container for Gen-level selection
 159                 Double_t vectorMC[6] = {9999.,9999.,9999.,9999.,9999.,9999.};
 160                 if (GetGeneratedValuesFromMCParticle(mcPart, mcArray, vectorMC)){
 161                         containerInput[0] = vectorMC[0];
 162                         containerInput[1] = vectorMC[1] ;
 163                         containerInput[2] = vectorMC[2] ;
 164                         containerInput[3] = vectorMC[3] ;
 165                         containerInput[4] = vectorMC[4] ;
 166                         containerInput[5] = vectorMC[5] ;  // in micron
 167                         fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated);
 168                         icountMC++;
 169                 }
 170                 else {
 171                         AliDebug(3,"Problems in filling the container");
 172                         continue;
 173                 }
 174         }
 175
 176         AliDebug(2, Form("Found %i MC particles that are D0!!",icountMC));
 177         AliDebug(2, Form("aodEvent = %p ",aodEvent));
 178
 179         // Now go to rec level
 180         fCountMC += icountMC;
 181         // load heavy flavour vertices
 182
 183         TClonesArray *arrayD0toKpi = (TClonesArray*)((aodEvent->GetList())->FindObject("D0toKpi"));
 184         if (!arrayD0toKpi) AliError("Could not find array of HF vertices");
 185         AliDebug(2, Form("Found %d vertices",arrayD0toKpi->GetEntriesFast()));
 186
 187         for (Int_t iD0toKpi = 0; iD0toKpi<arrayD0toKpi->GetEntriesFast(); iD0toKpi++) {
 188
 189                 AliAODRecoDecayHF2Prong* d0tokpi = (AliAODRecoDecayHF2Prong*)arrayD0toKpi->At(iD0toKpi);
 190
 191                 // cuts can't be applied to RecoDecays particles
 192                 // if (!fCFManager->CheckParticleCuts(1  , d0tokpi)) continue;  // 1 stands for AOD level
 193
 194                 // find associated MC particle
 195                 Int_t mcLabel = d0tokpi->MatchToMC(421,mcArray) ;
 196                 if (mcLabel == -1)
 197                         {
 198                                 AliDebug(2,"No MC particle found");
 199                         }
 200                 else {
 201                         AliAODMCParticle* mcVtxHF = (AliAODMCParticle*)mcArray->At(mcLabel);
 202                         if (!mcVtxHF) {
 203                                 AliWarning("Could not find associated MC in AOD MC tree");
 204                                 continue;
 205                         }
 206
 207                         // check if associated MC v0 passes the cuts
 208                         if (!fCFManager->CheckParticleCuts(0 ,mcVtxHF)) {  // 0 stands for MC
 209                                 AliDebug(2, "Skipping the particles due to cuts");
 210                                 continue;
 211                         }
 212
 213                         // fill the container
 214                         // either with reconstructed values....
 215                         Double_t pt = d0tokpi->Pt();
 216                         Double_t rapidity = d0tokpi->YD0();
 217
 218                         Double_t cosThetaStar = 9999.;
 219                         Double_t pTpi = 0.;
 220                         Double_t pTK = 0.;
 221                         Int_t pdgCode = mcVtxHF->GetPdgCode();
 222                         if (pdgCode > 0){
 223                                 cosThetaStar = d0tokpi->CosThetaStarD0();
 224                                 pTpi = d0tokpi->PtProng(0);
 225                                 pTK = d0tokpi->PtProng(1);
 226                         }
 227                         else {
 228                                 cosThetaStar = d0tokpi->CosThetaStarD0bar();
 229                                 pTpi = d0tokpi->PtProng(1);
 230                                 pTK = d0tokpi->PtProng(0);
 231                         }
 232
 233                         Double_t cT = d0tokpi->CtD0();
 234                         AliDebug(2, Form("cT from reconstructed vertex = %f (micron)",cT*1E4));
 235                         AliDebug(2, Form("pdg code from MC = %d",TMath::Abs(mcVtxHF->GetPdgCode())));
 236
 237                         // ... or with generated values
 238
 239                         if (!fFillFromGenerated){
 240                                 containerInput[0] = pt;
 241                                 containerInput[1] = rapidity;
 242                                 containerInput[2] = cosThetaStar;
 243                                 containerInput[3] = pTpi;
 244                                 containerInput[4] = pTK;
 245                                 containerInput[5] = cT*1.E4;  // in micron
 246                         }
 247                         else {
 248                                 Double_t vectorMC[6] = {9999.,9999.,9999.,9999.,9999.,9999.};
 249                                 if (GetGeneratedValuesFromMCParticle(mcVtxHF, mcArray, vectorMC)){
 250                                         containerInput[0] = vectorMC[0];
 251                                         containerInput[1] = vectorMC[1] ;
 252                                         containerInput[2] = vectorMC[2] ;
 253                                         containerInput[3] = vectorMC[3] ;
 254                                         containerInput[4] = vectorMC[4] ;
 255                                         containerInput[5] = vectorMC[5] ;  // in micron
 256                                 }
 257                                 else {
 258                                         AliDebug(3,"Problems in filling the container");
 259                                         continue;
 260                                 }
 261                         }
 262                         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]));
 263                         fCFManager->GetParticleContainer()->Fill(containerInput,1) ;
 264                 }
 265                 /*
 266                 // for debugging only, filling container with dummy values
 267                 Double_t pt = 1.5;
 268                 Double_t rapidity = 0.5;
 269                 Double_t cosThetaStar = 0.7;
 270                 Double_t pTpi = 2.5;
 271                 Double_t pTK = 1;
 272                 Double_t cT = 3;
 273                 containerInput[0] = pt;
 274                 containerInput[1] = rapidity;
 275                 containerInput[2] = cosThetaStar;
 276                 containerInput[3] = pTpi;
 277                 containerInput[4] = pTK;
 278                 containerInput[5] = cT;
 279                 fCFManager->GetParticleContainer()->Fill(containerInput,1) ;
 280                 */
 281
 282         }
 283
 284         fHistEventsProcessed->Fill(0);
 285         /* PostData(0) is taken care of by AliAnalysisTaskSE */
 286         PostData(1,fHistEventsProcessed) ;
 287         PostData(2,fCFManager->GetParticleContainer()) ;
 288 }
 289
 290
 291 //___________________________________________________________________________
 292 void AliCFHeavyFlavourTaskMultiVar::Terminate(Option_t*)
 293 {
 294         // The Terminate() function is the last function to be called during
 295         // a query. It always runs on the client, it can be used to present
 296         // the results graphically or save the results to file.
 297
 298         Info("Terminate","");
 299         AliAnalysisTaskSE::Terminate();
 300
 301         AliInfo(Form("Found %i MC particles that are D0 in MC in %d events",fCountMC,fEvents));
 302
 303         // draw some example plots....
 304
 305         AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
 306
 307         // projecting the containers to obtain histograms
 308         // first argument = variable, second argument = step
 309         TH1D* h00 =   cont->ShowProjection(0,0) ;   // pt
 310         TH1D* h10 =   cont->ShowProjection(1,0) ;   // rapidity
 311         TH1D* h20 =   cont->ShowProjection(2,0) ;   // cosThetaStar
 312         TH1D* h30 =   cont->ShowProjection(3,0) ;   // pTpi
 313         TH1D* h40 =   cont->ShowProjection(4,0) ;   // pTK
 314         TH1D* h50 =   cont->ShowProjection(5,0) ;   // cT
 315
 316         TH1D* h01 =   cont->ShowProjection(0,1) ;   // pt
 317         TH1D* h11 =   cont->ShowProjection(1,1) ;   // rapidity
 318         TH1D* h21 =   cont->ShowProjection(2,1) ;   // cosThetaStar
 319         TH1D* h31 =   cont->ShowProjection(3,1) ;   // pTpi
 320         TH1D* h41 =   cont->ShowProjection(4,1) ;   // pTK
 321         TH1D* h51 =   cont->ShowProjection(5,1) ;   // cT
 322
 323         h00->SetTitle("pT_D0 (GeV/c)");
 324         h10->SetTitle("rapidity");
 325         h20->SetTitle("cosThetaStar");
 326         h30->SetTitle("pT_pi (GeV/c)");
 327         h40->SetTitle("pT_K (Gev/c)");
 328         h50->SetTitle("cT (#mum)");
 329
 330         h00->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
 331         h10->GetXaxis()->SetTitle("rapidity");
 332         h20->GetXaxis()->SetTitle("cosThetaStar");
 333         h30->GetXaxis()->SetTitle("pT_pi (GeV/c)");
 334         h40->GetXaxis()->SetTitle("pT_K (Gev/c)");
 335         h50->GetXaxis()->SetTitle("cT (#mum)");
 336
 337         h01->SetTitle("pT_D0 (GeV/c)");
 338         h11->SetTitle("rapidity");
 339         h21->SetTitle("cosThetaStar");
 340         h31->SetTitle("pT_pi (GeV/c)");
 341         h41->SetTitle("pT_K (Gev/c)");
 342         h51->SetTitle("cT (#mum)");
 343
 344         h01->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
 345         h11->GetXaxis()->SetTitle("rapidity");
 346         h21->GetXaxis()->SetTitle("cosThetaStar");
 347         h31->GetXaxis()->SetTitle("pT_pi (GeV/c)");
 348         h41->GetXaxis()->SetTitle("pT_K (Gev/c)");
 349         h51->GetXaxis()->SetTitle("cT (#mum)");
 350
 351         Double_t max0 = h00->GetMaximum();
 352         Double_t max1 = h10->GetMaximum();
 353         Double_t max2 = h20->GetMaximum();
 354         Double_t max3 = h30->GetMaximum();
 355         Double_t max4 = h40->GetMaximum();
 356         Double_t max5 = h50->GetMaximum();
 357
 358         h00->GetYaxis()->SetRangeUser(0,max0*1.2);
 359         h10->GetYaxis()->SetRangeUser(0,max1*1.2);
 360         h20->GetYaxis()->SetRangeUser(0,max2*1.2);
 361         h30->GetYaxis()->SetRangeUser(0,max3*1.2);
 362         h40->GetYaxis()->SetRangeUser(0,max4*1.2);
 363         h50->GetYaxis()->SetRangeUser(0,max5*1.2);
 364
 365         h01->GetYaxis()->SetRangeUser(0,max0*1.2);
 366         h11->GetYaxis()->SetRangeUser(0,max1*1.2);
 367         h21->GetYaxis()->SetRangeUser(0,max2*1.2);
 368         h31->GetYaxis()->SetRangeUser(0,max3*1.2);
 369         h41->GetYaxis()->SetRangeUser(0,max4*1.2);
 370         h51->GetYaxis()->SetRangeUser(0,max5*1.2);
 371
 372         h00->SetMarkerStyle(20);
 373         h10->SetMarkerStyle(24);
 374         h20->SetMarkerStyle(21);
 375         h30->SetMarkerStyle(25);
 376         h40->SetMarkerStyle(27);
 377         h50->SetMarkerStyle(28);
 378
 379         h00->SetMarkerColor(2);
 380         h10->SetMarkerColor(2);
 381         h20->SetMarkerColor(2);
 382         h30->SetMarkerColor(2);
 383         h40->SetMarkerColor(2);
 384         h50->SetMarkerColor(2);
 385
 386         h01->SetMarkerStyle(20) ;
 387         h11->SetMarkerStyle(24) ;
 388         h21->SetMarkerStyle(21) ;
 389         h31->SetMarkerStyle(25) ;
 390         h41->SetMarkerStyle(27) ;
 391         h51->SetMarkerStyle(28) ;
 392
 393         h01->SetMarkerColor(4);
 394         h11->SetMarkerColor(4);
 395         h21->SetMarkerColor(4);
 396         h31->SetMarkerColor(4);
 397         h41->SetMarkerColor(4);
 398         h51->SetMarkerColor(4);
 399
 400         gStyle->SetCanvasColor(0);
 401         gStyle->SetFrameFillColor(0);
 402         gStyle->SetTitleFillColor(0);
 403         gStyle->SetStatColor(0);
 404
 405         // drawing in 2 separate canvas for a matter of clearity
 406         TCanvas * c1 =new TCanvas("c1","pT, rapidiy, cosThetaStar",800,1600);
 407         TCanvas * c2 =new TCanvas("c2","pTpi, pTK, cT",800,1600);
 408         c1->Divide(2,3);
 409         c2->Divide(2,3);
 410
 411         c1->cd(1);
 412         h00->Draw("p");
 413         c1->cd(2);
 414         h01->Draw("p");
 415         c1->cd(3);
 416         h10->Draw("p");
 417         c1->cd(4);
 418         h11->Draw("p");
 419         c1->cd(5);
 420         h20->Draw("p");
 421         c1->cd(6);
 422         h21->Draw("p");
 423         c1->cd();
 424
 425         c2->cd(1);
 426         h30->Draw("p");
 427         c2->cd(2);
 428         h31->Draw("p");
 429         c2->cd(3);
 430         h40->Draw("p");
 431         c2->cd(4);
 432         h41->Draw("p");
 433         c2->cd(5);
 434         h50->Draw("p");
 435         c2->cd(6);
 436         h51->Draw("p");
 437         c2->cd();
 438
 439         c1->SaveAs("Variables/pT_rapidity_cosThetaStar.eps");
 440         c2->SaveAs("Variables/pTpi_pTK_cT.eps");
 441
 442 }
 443
 444 //___________________________________________________________________________
 445 void AliCFHeavyFlavourTaskMultiVar::UserCreateOutputObjects() {
 446         //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
 447         //TO BE SET BEFORE THE EXECUTION OF THE TASK
 448         //
 449         Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 450
 451         //slot #1
 452         OpenFile(1);
 453         fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ;
 454 }
 455
 456 //___________________________________________________________________________
 457 Double_t AliCFHeavyFlavourTaskMultiVar::CosThetaStar(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
 458
 459         //
 460         // to calculate cos(ThetaStar) for generated particle
 461         // using the K, since mcPartDaughter0 and mcPartDaughter1 always correspond to K and pi respectively
 462         // (see where the function is called)
 463         //
 464
 465         Int_t pdgvtx = mcPart->GetPdgCode();
 466         /*      if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
 467                 Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
 468                 Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
 469                 AliInfo(Form("D0, with pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
 470                 AliDebug(2,"This is a D0");
 471                 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
 472                 mcPartDaughter0 = mcPartDaughter1;
 473                 mcPartDaughter1 = mcPartdummy;
 474         }
 475         else{
 476                 AliInfo("D0bar");
 477         }
 478         */
 479         Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
 480         Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
 481         if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar)
 482                 AliDebug(2,"D0");
 483         }
 484         else{
 485                 AliDebug(2,"D0bar");
 486         }
 487         if (pdgprong0 == 211){
 488                 AliDebug(2,Form("pdgprong0 = %d, pdgprong1 = %d, switching...",pdgprong0,pdgprong1));
 489                 AliAODMCParticle* mcPartdummy = mcPartDaughter0;
 490                 mcPartDaughter0 = mcPartDaughter1;
 491                 mcPartDaughter1 = mcPartdummy;
 492                 pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode());
 493                 pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode());
 494         }
 495
 496         AliDebug(2,Form("After checking, pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1));
 497         Double_t massvtx = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdgvtx))->Mass();
 498         Double_t massp[2];
 499         massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass();
 500         massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass();
 501
 502         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);
 503
 504         Double_t px = mcPartDaughter0->Px()+mcPartDaughter1->Px();
 505         Double_t py = mcPartDaughter0->Py()+mcPartDaughter1->Py();
 506         Double_t pz = mcPartDaughter0->Pz()+mcPartDaughter1->Pz();
 507         Double_t p = TMath::Sqrt(px*px+py*py+pz*pz);
 508         Double_t e =  TMath::Sqrt(massvtx*massvtx+p*p);
 509         Double_t beta = p/e;
 510         Double_t gamma = e/massvtx;
 511         TVector3 mom(mcPartDaughter0->Px(),mcPartDaughter0->Py(),mcPartDaughter0->Pz());
 512         TVector3 momTot(mcPartDaughter0->Px()+mcPartDaughter1->Px(),mcPartDaughter0->Py()+mcPartDaughter1->Py(),mcPartDaughter0->Pz()+mcPartDaughter1->Pz());
 513
 514         Double_t qlprong = mom.Dot(momTot)/momTot.Mag();  // analog to AliAODRecoDecay::QlProng(0)
 515
 516         AliDebug(2,Form("pStar = %f, beta = %f, gamma = %f, qlprong = %f, massp[0] = %f", pStar, beta, gamma, qlprong, massp[0]));
 517         Double_t cts = (qlprong/gamma-beta*TMath::Sqrt(pStar*pStar+massp[0]*massp[0]))/pStar;
 518         AliDebug(2,Form("cts = %f", cts));
 519         return cts;
 520 }
 521 //___________________________________________________________________________
 522 Double_t AliCFHeavyFlavourTaskMultiVar::CT(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const {
 523
 524         //
 525         // to calculate cT for generated particle
 526         //
 527
 528         Double_t xmcPart[3] = {0,0,0};
 529         Double_t xdaughter0[3] = {0,0,0};
 530         Double_t xdaughter1[3] = {0,0,0};
 531         mcPart->XvYvZv(xmcPart);  // cm
 532         mcPartDaughter0->XvYvZv(xdaughter0);  // cm
 533         mcPartDaughter1->XvYvZv(xdaughter1);  //cm
 534         Double_t prodVtxD0 = TMath::Sqrt(xmcPart[0]*xmcPart[0]+
 535                                          xmcPart[1]*xmcPart[1]+
 536                                          xmcPart[2]*xmcPart[2]);
 537         Double_t prodVtxDaughter0 = TMath::Sqrt(xdaughter0[0]*xdaughter0[0]+
 538                                                 xdaughter0[1]*xdaughter0[1]+
 539                                                 xdaughter0[2]*xdaughter0[2]);
 540         Double_t prodVtxDaughter1 = TMath::Sqrt(xdaughter1[0]*xdaughter1[0]+
 541                                                 xdaughter1[1]*xdaughter1[1]+
 542                                                 xdaughter1[2]*xdaughter1[2]);
 543
 544         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));
 545         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));
 546         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));
 547
 548         Double_t cT0 = TMath::Sqrt((xdaughter0[0]-xmcPart[0])*(xdaughter0[0]-xmcPart[0])+
 549                                      (xdaughter0[1]-xmcPart[1])*(xdaughter0[1]-xmcPart[1])+
 550                                      (xdaughter0[2]-xmcPart[2])*(xdaughter0[2]-xmcPart[2]));
 551
 552         Double_t cT1 = TMath::Sqrt((xdaughter1[0]-xmcPart[0])*(xdaughter1[0]-xmcPart[0])+
 553                                      (xdaughter1[1]-xmcPart[1])*(xdaughter1[1]-xmcPart[1])+
 554                                      (xdaughter1[2]-xmcPart[2])*(xdaughter1[2]-xmcPart[2]));
 555
 556         if (cT0 != cT1) {
 557                 AliWarning("cT from daughter 0 different from cT from daughter 1! Using cT from daughter 0, but PLEASE, CHECK....");
 558         }
 559
 560         // calculating cT from cT0
 561
 562         Double_t mass = TDatabasePDG::Instance()->GetParticle(mcPart->GetPdgCode())->Mass(); // mcPart->GetPdgCode() should return +/- 421 for the D0/D0bar
 563         Double_t p = mcPart-> P();
 564         Double_t cT = cT0*mass/p;
 565         AliDebug(2, Form("cT from daughter 0 = %f (micron)", cT0*1E4));
 566         AliDebug(2, Form("cT from daughter 1 = %f (micron)", cT1*1E4));
 567         AliDebug(2, Form("cT (with mass = %f and p = %f) = %f (micron)", mass, p, cT*1E4));
 568         return cT;
 569 }
 570 //________________________________________________________________________________
 571 Bool_t AliCFHeavyFlavourTaskMultiVar::GetGeneratedValuesFromMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC) const {
 572
 573         //
 574         // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle
 575         //
 576
 577         Bool_t isOk = kFALSE;
 578
 579         // check whether the D0 decays in pi+K
 580         // to be added!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 581         // could use a cut in the CF, but not clear how to define a TDecayChannel
 582         // implemented for the time being as a cut on the number of daughters - see later when
 583         // getting the daughters
 584
 585         // getting the daughters
 586         Int_t daughter0 = mcPart->GetDaughter(0);
 587         Int_t daughter1 = mcPart->GetDaughter(1);
 588         AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1));
 589         if (daughter0 == 0 || daughter1 == 0) {
 590                 AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!");
 591                 return isOk;
 592         }
 593         if (TMath::Abs(daughter1 - daughter0 != 1)) {
 594                 AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!");
 595                 return isOk;
 596         }
 597         AliAODMCParticle* mcPartDaughter0 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter0));
 598         AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(mcArray->At(daughter1));
 599         if (!mcPartDaughter0 || !mcPartDaughter1) {
 600                 AliWarning("At least one Daughter Particle not found in tree, skipping");
 601                 return isOk;
 602         }
 603
 604         Double_t vtx1[3] = {0,0,0};   // primary vertex
 605         Double_t vtx2daughter0[3] = {0,0,0};   // secondary vertex from daughter 0
 606         Double_t vtx2daughter1[3] = {0,0,0};   // secondary vertex from daughter 1
 607         mcPart->XvYvZv(vtx1);  // cm
 608         // getting vertex from daughters
 609         mcPartDaughter0->XvYvZv(vtx2daughter0);  // cm
 610         mcPartDaughter1->XvYvZv(vtx2daughter1);  //cm
 611         if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) {
 612                 AliError("Daughters have different secondary vertex, skipping the track");
 613                 return isOk;
 614         }
 615         Int_t nprongs = 2;
 616         Short_t charge = 0;
 617         // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second
 618         AliAODMCParticle* positiveDaugh = mcPartDaughter0;
 619         AliAODMCParticle* negativeDaugh = mcPartDaughter1;
 620         if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){
 621                 // inverting in case the positive daughter is the second one
 622                 positiveDaugh = mcPartDaughter1;
 623                 negativeDaugh = mcPartDaughter0;
 624         }
 625         // getting the momentum from the daughters
 626         Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()};
 627         Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()};
 628         Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()};
 629
 630         Double_t d0[2] = {0.,0.};
 631
 632         AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0);
 633
 634         Double_t cosThetaStar = 0.;
 635         Double_t cosThetaStarD0 = 0.;
 636         Double_t cosThetaStarD0bar = 0.;
 637         cosThetaStarD0 = decay->CosThetaStar(1,421,211,321);
 638         cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211);
 639         if (mcPart->GetPdgCode() == 421){  // D0
 640                 AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
 641                 cosThetaStar = cosThetaStarD0;
 642         }
 643         else if (mcPart->GetPdgCode() == -421){  // D0bar{
 644                 AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode()));
 645                 cosThetaStar = cosThetaStarD0bar;
 646         }
 647         else{
 648                 AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check...");
 649                 return vectorMC;
 650         }
 651         if (cosThetaStar < -1 || cosThetaStar > 1) {
 652                 AliWarning("Invalid value for cosine Theta star, returning");
 653                 return isOk;
 654         }
 655
 656         // calculate cos(Theta*) according to the method implemented herein
 657
 658         Double_t cts = 9999.;
 659         cts = CosThetaStar(mcPart, mcPartDaughter0, mcPartDaughter1);
 660         if (cts < -1 || cts > 1) {
 661                 AliWarning("Invalid value for cosine Theta star from AliCFHeavyFlavourTaskMultiVar method");
 662         }
 663         if (TMath::Abs(cts - cosThetaStar)>0.001) {
 664                 AliError(Form("cosThetaStar automatically calculated different from that manually calculated!!! cosThetaStar = %f, cosThetaStar = %f", cosThetaStar,cts));
 665         }
 666
 667         Double_t cT = decay->Ct(421);
 668
 669         // calculate cT from the method implemented herein
 670         Double_t cT1 = 0.;
 671         cT1 = CT(mcPart, mcPartDaughter0, mcPartDaughter1);
 672
 673         if (TMath::Abs(cT1 - cT)>0.001) {
 674                 AliError(Form("cT automatically calculated different from that manually calculated!!! cT = %f, cT1 = %f",cT,cT1));
 675         }
 676
 677         // get the pT of the daughters
 678
 679         Double_t pTpi = 0.;
 680         Double_t pTK = 0.;
 681
 682         if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) {
 683                 pTpi = mcPartDaughter0->Pt();
 684                 pTK = mcPartDaughter1->Pt();
 685         }
 686         else {
 687                 pTpi = mcPartDaughter1->Pt();
 688                 pTK = mcPartDaughter0->Pt();
 689         }
 690
 691         vectorMC[0] = mcPart->Pt();
 692         vectorMC[1] = mcPart->Y() ;
 693         vectorMC[2] = cosThetaStar ;
 694         vectorMC[3] = pTpi ;
 695         vectorMC[4] = pTK ;
 696         vectorMC[5] = cT*1.E4 ;  // in micron
 697         isOk = kTRUE;
 698         return isOk;
 699 }
 700