PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVarMultiStep.cxx

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