PWG3/vertexingHF/AliCFTaskVertexingHF.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 //          D. Caffarri, Univ & INFN Padova  caffarri@pd.infn.it
  26 //-----------------------------------------------------------------------
  27 //-----------------------------------------------------------------------
  28 // Base class for HF Unfolding (pt and eta)
  29 // correlation matrix filled at Acceptance and PPR level
  30 // Author: A.Grelli ,  Utrecht - agrelli@uu.nl
  31 //-----------------------------------------------------------------------
  32 #include <TCanvas.h>
  33 #include <TParticle.h>
  34 #include <TDatabasePDG.h>
  35 #include <TH1I.h>
  36 #include <TStyle.h>
  37 #include <TFile.h>
  38
  39 #include "AliCFTaskVertexingHF.h"
  40 #include "AliStack.h"
  41 #include "AliMCEvent.h"
  42 #include "AliCFManager.h"
  43 #include "AliCFContainer.h"
  44 #include "AliLog.h"
  45 #include "AliAnalysisManager.h"
  46 #include "AliAODHandler.h"
  47 #include "AliAODEvent.h"
  48 #include "AliAODRecoDecay.h"
  49 #include "AliAODRecoDecayHF.h"
  50 #include "AliAODRecoDecayHF2Prong.h"
  51 #include "AliAODRecoDecayHF3Prong.h"
  52 #include "AliAODRecoDecayHF4Prong.h"
  53 #include "AliAODRecoCascadeHF.h"
  54 #include "AliAODMCParticle.h"
  55 #include "AliAODMCHeader.h"
  56 #include "AliESDtrack.h"
  57 #include "TChain.h"
  58 #include "THnSparse.h"
  59 #include "TH2D.h"
  60 #include "AliESDtrackCuts.h"
  61 #include "AliRDHFCuts.h"
  62 #include "AliRDHFCutsD0toKpi.h"
  63 #include "AliRDHFCutsDplustoKpipi.h"
  64 #include "AliRDHFCutsDStartoKpipi.h"
  65 #include "AliRDHFCutsDstoKKpi.h"
  66 #include "AliRDHFCutsLctopKpi.h"
  67 #include "AliRDHFCutsD0toKpipipi.h"
  68 #include "AliCFVertexingHF2Prong.h"
  69 #include "AliCFVertexingHF3Prong.h"
  70 #include "AliCFVertexingHFCascade.h"
  71 #include "AliCFVertexingHF.h"
  72 #include "AliAnalysisDataSlot.h"
  73 #include "AliAnalysisDataContainer.h"
  74
  75 //__________________________________________________________________________
  76 AliCFTaskVertexingHF::AliCFTaskVertexingHF() :
  77         AliAnalysisTaskSE(),
  78         fCFManager(0x0),
  79         fHistEventsProcessed(0x0),
  80         fCorrelation(0x0),
  81         fCountMC(0),
  82         fCountAcc(0),
  83         fCountVertex(0),
  84         fCountRefit(0),
  85         fCountReco(0),
  86         fCountRecoAcc(0),
  87         fCountRecoITSClusters(0),
  88         fCountRecoPPR(0),
  89         fCountRecoPID(0),
  90         fEvents(0),
  91         fDecayChannel(0),
  92         fFillFromGenerated(kFALSE),
  93         fOriginDselection(0),
  94         fAcceptanceUnf(kTRUE),
  95         fCuts(0),
  96         fUseWeight(kFALSE),
  97         fWeight(1.),
  98         fNvar(0),
  99         fPartName(""),
 100         fDauNames(""),
 101         fSign(2),
 102         fCentralitySelection(kTRUE),
 103         fFakeSelection(0)
 104 {
 105         //
 106         //Default ctor
 107         //
 108 }
 109 //___________________________________________________________________________
 110 AliCFTaskVertexingHF::AliCFTaskVertexingHF(const Char_t* name, AliRDHFCuts* cuts) :
 111         AliAnalysisTaskSE(name),
 112         fCFManager(0x0),
 113         fHistEventsProcessed(0x0),
 114         fCorrelation(0x0),
 115         fCountMC(0),
 116         fCountAcc(0),
 117         fCountVertex(0),
 118         fCountRefit(0),
 119         fCountReco(0),
 120         fCountRecoAcc(0),
 121         fCountRecoITSClusters(0),
 122         fCountRecoPPR(0),
 123         fCountRecoPID(0),
 124         fEvents(0),
 125         fDecayChannel(0),
 126         fFillFromGenerated(kFALSE),
 127         fOriginDselection(0),
 128         fAcceptanceUnf(kTRUE),
 129         fCuts(cuts),
 130         fUseWeight(kFALSE),
 131         fWeight(1.),
 132         fNvar(0),
 133         fPartName(""),
 134         fDauNames(""),
 135         fSign(2),
 136         fCentralitySelection(kTRUE),
 137         fFakeSelection(0)
 138 {
 139         //
 140         // Constructor. Initialization of Inputs and Outputs
 141         //
 142         /*
 143           DefineInput(0) and DefineOutput(0)
 144           are taken care of by AliAnalysisTaskSE constructor
 145         */
 146         DefineOutput(1,TH1I::Class());
 147         DefineOutput(2,AliCFContainer::Class());
 148         DefineOutput(3,THnSparseD::Class());
 149         DefineOutput(4,AliRDHFCuts::Class());
 150
 151         fCuts->PrintAll();
 152 }
 153
 154 //___________________________________________________________________________
 155 AliCFTaskVertexingHF& AliCFTaskVertexingHF::operator=(const AliCFTaskVertexingHF& c)
 156 {
 157         //
 158         // Assignment operator
 159         //
 160         if (this!=&c) {
 161                 AliAnalysisTaskSE::operator=(c) ;
 162                 fCFManager  = c.fCFManager;
 163                 fHistEventsProcessed = c.fHistEventsProcessed;
 164                 fCuts = c.fCuts;
 165         }
 166         return *this;
 167 }
 168
 169 //___________________________________________________________________________
 170 AliCFTaskVertexingHF::AliCFTaskVertexingHF(const AliCFTaskVertexingHF& c) :
 171         AliAnalysisTaskSE(c),
 172         fCFManager(c.fCFManager),
 173         fHistEventsProcessed(c.fHistEventsProcessed),
 174         fCorrelation(c.fCorrelation),
 175         fCountMC(c.fCountMC),
 176         fCountAcc(c.fCountAcc),
 177         fCountVertex(c.fCountVertex),
 178         fCountRefit(c.fCountRefit),
 179         fCountReco(c.fCountReco),
 180         fCountRecoAcc(c.fCountRecoAcc),
 181         fCountRecoITSClusters(c.fCountRecoITSClusters),
 182         fCountRecoPPR(c.fCountRecoPPR),
 183         fCountRecoPID(c.fCountRecoPID),
 184         fEvents(c.fEvents),
 185         fDecayChannel(c.fDecayChannel),
 186         fFillFromGenerated(c.fFillFromGenerated),
 187         fOriginDselection(c.fOriginDselection),
 188         fAcceptanceUnf(c.fAcceptanceUnf),
 189         fCuts(c.fCuts),
 190         fUseWeight(c.fUseWeight),
 191         fWeight(c.fWeight),
 192         fNvar(c.fNvar),
 193         fPartName(c.fPartName),
 194         fDauNames(c.fDauNames),
 195         fSign(c.fSign),
 196         fCentralitySelection(c.fCentralitySelection),
 197         fFakeSelection(c.fFakeSelection)
 198 {
 199         //
 200         // Copy Constructor
 201         //
 202 }
 203
 204 //___________________________________________________________________________
 205 AliCFTaskVertexingHF::~AliCFTaskVertexingHF()
 206 {
 207         //
 208         //destructor
 209         //
 210         if (fCFManager)           delete fCFManager ;
 211         if (fHistEventsProcessed) delete fHistEventsProcessed ;
 212         if (fCorrelation)         delete fCorrelation ;
 213         if (fCuts)                delete fCuts;
 214 }
 215
 216 //_________________________________________________________________________-
 217 void AliCFTaskVertexingHF::Init()
 218 {
 219         //
 220         // Initialization
 221         //
 222
 223         if (fDebug>1) printf("AliCFTaskVertexingHF::Init()");
 224         AliRDHFCuts *copyfCuts = 0x0;
 225         if (!fCuts){
 226                 AliFatal("No cuts defined - Exiting...");
 227                 return;
 228         }
 229
 230         switch (fDecayChannel){
 231         case 2:{
 232                 copyfCuts = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
 233                 fNvar = 14;
 234                 fPartName="D0";
 235                 fDauNames="K+pi";
 236                 break;
 237         }
 238         case 21:{
 239                 copyfCuts = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
 240                 fNvar = 14;
 241                 fPartName="Dstar";
 242                 fDauNames="K+pi+pi";
 243                 break;
 244         }
 245         case 31:{
 246                 copyfCuts = new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fCuts)));
 247                 fNvar = 13;
 248                 fPartName="Dplus";
 249                 fDauNames="K+pi+pi";
 250                 break;
 251         }
 252         case 32:{
 253                 copyfCuts = new AliRDHFCutsLctopKpi(*(static_cast<AliRDHFCutsLctopKpi*>(fCuts)));
 254                 fNvar = 13;
 255                 fPartName="Lambdac";
 256                 fDauNames="p+K+pi";
 257                 break;
 258         }
 259         case 33:{
 260                 copyfCuts = new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fCuts)));
 261                 fNvar = 13;
 262                 fPartName="Ds";
 263                 fDauNames="K+K+pi";
 264                 break;
 265         }
 266         case 4:{
 267                 copyfCuts = new AliRDHFCutsD0toKpipipi(*(static_cast<AliRDHFCutsD0toKpipipi*>(fCuts)));
 268                 fNvar = 14;
 269                 fPartName="D0";
 270                 fDauNames="K+pi+pi+pi";
 271                 break;
 272         }
 273         default:
 274                 AliFatal("The decay channel MUST be defined according to AliCFVertexing::DecayChannel - Exiting...");
 275                 break;
 276         }
 277
 278         const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
 279         if (copyfCuts){
 280                 copyfCuts->SetName(nameoutput);
 281
 282                 //Post the data
 283                 PostData(4, copyfCuts);
 284         }
 285         else{
 286                 AliFatal("Failing initializing AliRDHFCuts object - Exiting...");
 287         }
 288
 289         return;
 290 }
 291
 292 //_________________________________________________
 293 void AliCFTaskVertexingHF::UserExec(Option_t *)
 294 {
 295         //
 296         // Main loop function
 297         //
 298
 299         PostData(1,fHistEventsProcessed) ;
 300         PostData(2,fCFManager->GetParticleContainer()) ;
 301         PostData(3,fCorrelation) ;
 302
 303
 304         if (fFillFromGenerated){
 305                 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
 306         }
 307
 308         if (!fInputEvent) {
 309                 Error("UserExec","NO EVENT FOUND!");
 310                 return;
 311         }
 312
 313         AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 314
 315         TClonesArray *arrayBranch=0;
 316
 317         if(!aodEvent && AODEvent() && IsStandardAOD()) {
 318                 // In case there is an AOD handler writing a standard AOD, use the AOD
 319                 // event in memory rather than the input (ESD) event.
 320                 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
 321                 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 322                 // have to taken from the AOD event hold by the AliAODExtension
 323                 AliAODHandler* aodHandler = (AliAODHandler*)
 324                         ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 325                 if(aodHandler->GetExtensions()) {
 326                         AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 327                         AliAODEvent *aodFromExt = ext->GetAOD();
 328
 329                         switch (fDecayChannel){
 330                         case 2:{
 331                                 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
 332                                 break;
 333                         }
 334                         case 21:{
 335                                 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
 336                                 break;
 337                         }
 338                         case 31:
 339                         case 32:
 340                         case 33:{
 341                                 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 342                                 break;
 343                         }
 344                         case 4:{
 345                                 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
 346                                 break;
 347                         }
 348                         default:
 349                                 break;
 350                         }
 351                 }
 352         }
 353         else {
 354                 switch (fDecayChannel){
 355                 case 2:{
 356                         arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
 357                         break;
 358                 }
 359                 case 21:{
 360                         arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
 361                         break;
 362                 }
 363                 case 31:
 364                 case 32:
 365                 case 33:{
 366                         arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
 367                         break;
 368                 }
 369                 case 4:{
 370                         arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm4Prong");
 371                         break;
 372                 }
 373                 default:
 374                         break;
 375                 }
 376         }
 377
 378         AliAODVertex *aodVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
 379         if (!aodVtx) return;
 380
 381         if (!arrayBranch) {
 382                 AliError("Could not find array of HF vertices");
 383                 return;
 384         }
 385
 386         fEvents++;
 387
 388         fCFManager->SetRecEventInfo(aodEvent);
 389         fCFManager->SetMCEventInfo(aodEvent);
 390
 391         //******** DEFINE number of variables of the container***** for now set at 13, in the future in the config macro.
 392
 393         //loop on the MC event
 394
 395         TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 396         if (!mcArray) {
 397                 AliError("Could not find Monte-Carlo in AOD");
 398                 return;
 399         }
 400         Int_t icountMC = 0;
 401         Int_t icountAcc = 0;
 402         Int_t icountReco = 0;
 403         Int_t icountVertex = 0;
 404         Int_t icountRefit = 0;
 405         Int_t icountRecoAcc = 0;
 406         Int_t icountRecoITSClusters = 0;
 407         Int_t icountRecoPPR = 0;
 408         Int_t icountRecoPID = 0;
 409         Int_t cquarks = 0;
 410
 411         AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 412         if (!mcHeader) {
 413                 AliError("Could not find MC Header in AOD");
 414                 return;
 415         }
 416
 417         Double_t* containerInput = new Double_t[fNvar];
 418         Double_t* containerInputMC = new Double_t[fNvar];
 419
 420
 421         AliCFVertexingHF* cfVtxHF=0x0;
 422         switch (fDecayChannel){
 423         case 2:{
 424           cfVtxHF = new AliCFVertexingHF2Prong(mcArray, fOriginDselection);
 425           break;
 426         }
 427         case 21:{
 428           cfVtxHF = new AliCFVertexingHFCascade(mcArray, fOriginDselection);
 429           break;
 430         }
 431         case 31:
 432         case 32:
 433         case 33:{
 434           cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel);
 435                 break;
 436         }
 437         case 4:{
 438                 //cfVtxHF = new AliCFVertexingHF4Prong(mcArray, originDselection);  // not there yet
 439                 break;
 440         }
 441         default:
 442                 break;
 443         }
 444         if (!cfVtxHF){
 445                 AliError("No AliCFVertexingHF initialized");
 446                 delete[] containerInput;
 447                 delete[] containerInputMC;
 448                 return;
 449         }
 450
 451         Double_t zPrimVertex = aodVtx ->GetZ();
 452         Double_t zMCVertex = mcHeader->GetVtxZ();
 453
 454         AliESDtrackCuts** trackCuts = new AliESDtrackCuts*[cfVtxHF->GetNProngs()];
 455         if (fDecayChannel == 21){
 456                 // for the D*, setting the third element of the array of the track cuts to those for the soft pion
 457                 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs()-1; iProng++){
 458                         trackCuts[iProng]=fCuts->GetTrackCuts();
 459                 }
 460                 trackCuts[2] = fCuts->GetTrackCutsSoftPi();
 461         }
 462         else {
 463                 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs(); iProng++){
 464                         trackCuts[iProng]=fCuts->GetTrackCuts();
 465                 }
 466         }
 467
 468         //General settings: vertex, feed down and fill reco container with generated values.
 469         cfVtxHF->SetRecoPrimVertex(zPrimVertex);
 470         cfVtxHF->SetMCPrimaryVertex(zMCVertex);
 471         cfVtxHF->SetFillFromGenerated(fFillFromGenerated);
 472         cfVtxHF->SetNVar(fNvar);
 473         cfVtxHF->SetFakeSelection(fFakeSelection);
 474
 475         if (fCentralitySelection)
 476           if(fCuts->IsEventSelectedInCentrality(aodEvent)!=0) {
 477             delete[] containerInput;
 478             delete[] containerInputMC;
 479             delete [] trackCuts;
 480             return;
 481           }
 482
 483         Float_t centValue = fCuts->GetCentrality(aodEvent);
 484         cfVtxHF->SetCentralityValue(centValue);
 485
 486         for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
 487           AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
 488           if (!mcPart){
 489             AliError("Failed casting particle from MC array!, Skipping particle");
 490             continue;
 491           }
 492           // check the MC-level cuts, must be the desidered particle
 493           if (!fCFManager->CheckParticleCuts(0, mcPart)) {
 494             continue;  // 0 stands for MC level
 495           }
 496           cfVtxHF->SetMCCandidateParam(iPart);
 497
 498           //counting c quarks
 499           cquarks += cfVtxHF->MCcquarkCounting(mcPart);
 500
 501           if (!(cfVtxHF->SetLabelArray())){
 502             AliDebug(2,Form("Impossible to set the label array (decaychannel = %d)",fDecayChannel));
 503             continue;
 504           }
 505
 506           //check the candiate family at MC level
 507           if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) {
 508             AliDebug(2,Form("Check on the family wrong!!! (decaychannel = %d)",fDecayChannel));
 509             continue;
 510           }
 511                 else{
 512                         AliDebug(2,Form("Check on the family OK!!! (decaychannel = %d)",fDecayChannel));
 513                 }
 514
 515                 //Fill the MC container
 516                 Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC);
 517                 if (mcContainerFilled) {
 518                         if (fUseWeight)fWeight = GetWeight(containerInputMC[0]);
 519                         if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) continue;
 520                         //MC Limited Acceptance
 521                         if (TMath::Abs(containerInputMC[1]) < 0.5) {
 522                                 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc, fWeight);
 523                                 AliDebug(3,"MC Lim Acc container filled\n");
 524                         }
 525
 526                         //MC
 527                         fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated, fWeight);
 528                         icountMC++;
 529                         AliDebug(3,"MC cointainer filled \n");
 530
 531                         // MC in acceptance
 532                         // check the MC-Acceptance level cuts
 533                         // since standard CF functions are not applicable, using Kine Cuts on daughters
 534                         Bool_t mcAccepStep = cfVtxHF-> MCAcceptanceStep();
 535                         if (mcAccepStep){
 536                                 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance, fWeight);
 537                                 AliDebug(3,"MC acceptance cut passed\n");
 538                                 icountAcc++;
 539
 540                                 //MC Vertex step
 541                                 if (fCuts->IsEventSelected(aodEvent)){
 542                                         // filling the container if the vertex is ok
 543                                         fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex, fWeight) ;
 544                                         AliDebug(3,"Vertex cut passed and container filled\n");
 545                                         icountVertex++;
 546
 547                                         //mc Refit requirement
 548                                         Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, &trackCuts[0]);
 549                                         if (mcRefitStep){
 550                                                 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit, fWeight);
 551                                                 AliDebug(3,"MC Refit cut passed and container filled\n");
 552                                                 icountRefit++;
 553
 554
 555                                         }
 556                                         else{
 557                                                 AliDebug(3,"MC Refit cut not passed\n");
 558                                                 continue;
 559                                         }
 560                                 }
 561                                 else{
 562                                   AliDebug (3, "MC vertex step not passed\n");
 563                                   continue;
 564                                 }
 565                         }
 566                         else{
 567                                 AliDebug (3, "MC in acceptance step not passed\n");
 568                                 continue;
 569                         }
 570                 }
 571                 else {
 572                         AliDebug (3, "MC container not filled\n");
 573                 }
 574         }
 575
 576         if (cquarks<2) AliDebug(2,Form("Event with %d c-quarks", cquarks));
 577         AliDebug(2,Form("Found %i MC particles that are %s!!",icountMC,fPartName.Data()));
 578         AliDebug(2,Form("Found %i MC particles that are %s and satisfy Acc cuts!!",icountAcc,fPartName.Data()));
 579         AliDebug(2,Form("Found %i MC particles that are %s and satisfy Vertex cuts!!",icountVertex,fPartName.Data()));
 580         AliDebug(2,Form("Found %i MC particles that are %s and satisfy Refit cuts!!",icountRefit,fPartName.Data()));
 581
 582         // Now go to rec level
 583         fCountMC += icountMC;
 584         fCountAcc += icountAcc;
 585         fCountVertex+= icountVertex;
 586         fCountRefit+= icountRefit;
 587
 588         AliDebug(2,Form("Found %d vertices for decay channel %d",arrayBranch->GetEntriesFast(),fDecayChannel));
 589
 590         for(Int_t iCandid = 0; iCandid<arrayBranch->GetEntriesFast();iCandid++){
 591                 AliAODRecoDecayHF* charmCandidate=0x0;
 592                 switch (fDecayChannel){
 593                 case 2:{
 594                         charmCandidate = (AliAODRecoDecayHF2Prong*)arrayBranch->At(iCandid);
 595                         break;
 596                 }
 597                 case 21:{
 598                         charmCandidate = (AliAODRecoCascadeHF*)arrayBranch->At(iCandid);
 599                         break;
 600                 }
 601                 case 31:
 602                 case 32:
 603                 case 33:{
 604                         charmCandidate = (AliAODRecoDecayHF3Prong*)arrayBranch->At(iCandid);
 605                         break;
 606                 }
 607                 case 4:{
 608                         charmCandidate = (AliAODRecoDecayHF4Prong*)arrayBranch->At(iCandid);
 609                         break;
 610                 }
 611                 default:
 612                         break;
 613                 }
 614
 615                 Bool_t unsetvtx=kFALSE;
 616                 if(!charmCandidate->GetOwnPrimaryVtx()) {
 617                         charmCandidate->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables
 618                         unsetvtx=kTRUE;
 619                 }
 620
 621                 Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)charmCandidate);
 622                 if (!signAssociation){
 623                         charmCandidate = 0x0;
 624                         continue;
 625                 }
 626
 627                 Int_t isPartOrAntipart = cfVtxHF->CheckReflexion(fSign);
 628                 if (isPartOrAntipart == 0){
 629                         AliDebug(2, Form("The candidate pdg code doesn't match the requirement set in the task (fSign = %d)",fSign));
 630                         continue;
 631                 }
 632
 633                 Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
 634                 if (recoContFilled){
 635
 636                         if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue;
 637
 638                         //Reco Step
 639                         Bool_t recoStep = cfVtxHF->RecoStep();
 640                         Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
 641
 642                         if (recoStep && recoContFilled && vtxCheck){
 643                                 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ;
 644                                 icountReco++;
 645                                 AliDebug(3,"Reco step  passed and container filled\n");
 646
 647                                 //Reco in the acceptance -- take care of UNFOLDING!!!!
 648                                 Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(&trackCuts[0]);
 649                                 if (recoAcceptanceStep) {
 650                                         fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ;
 651                                         icountRecoAcc++;
 652                                         AliDebug(3,"Reco acceptance cut passed and container filled\n");
 653
 654                                         if(fAcceptanceUnf){
 655                                                 Double_t fill[4]; //fill response matrix
 656                                                 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
 657                                                 if (bUnfolding) fCorrelation->Fill(fill);
 658                                         }
 659
 660                                         //Number of ITS cluster requirements
 661                                         Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks);
 662                                         if (recoITSnCluster){
 663                                                 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ;
 664                                                 icountRecoITSClusters++;
 665                                                 AliDebug(3,"Reco n ITS cluster cut passed and container filled\n");
 666
 667                                                 Bool_t iscutsusingpid = fCuts->GetIsUsePID();
 668                                                 Int_t recoAnalysisCuts = -1, recoPidSelection = -1;
 669                                                 fCuts->SetUsePID(kFALSE);
 670                                                 recoAnalysisCuts = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
 671
 672                                                 if (recoAnalysisCuts > 3){ // Ds case, where more possibilities are considered
 673                                                         if (recoAnalysisCuts >= 8){
 674                                                                 recoAnalysisCuts -= 8; // removing K0star mass
 675                                                         }
 676                                                         if (recoAnalysisCuts >= 4){
 677                                                                 recoAnalysisCuts -= 4; // removing Phi mass
 678                                                         }
 679                                                 }
 680
 681                                                 fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object
 682                                                 if (recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart){
 683                                                         fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR, fWeight);
 684                                                         icountRecoPPR++;
 685                                                         AliDebug(3,"Reco Analysis cuts passed and container filled \n");
 686                                                         //pid selection
 687                                                         //recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kPID);
 688                                                         //if((fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==isPartOrAntipart)||(fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==3)){
 689                                                         recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
 690                                                         if (recoPidSelection == 3 || recoPidSelection == isPartOrAntipart){
 691                                                                 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID, fWeight);
 692                                                                 icountRecoPID++;
 693                                                                 AliDebug(3,"Reco PID cuts passed and container filled \n");
 694                                                                 if(!fAcceptanceUnf){
 695                                                                         Double_t fill[4]; //fill response matrix
 696                                                                         Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
 697                                                                         if (bUnfolding) fCorrelation->Fill(fill);
 698                                                                 }
 699                                                         }
 700                                                         else {
 701                                                                 AliDebug(3, "Analysis Cuts step not passed \n");
 702                                                                 continue;
 703                                                         }
 704                                                 }
 705                                                 else {
 706                                                         AliDebug(3, "PID selection not passed \n");
 707                                                         continue;
 708                                                 }
 709                                         }
 710                                         else{
 711                                                 AliDebug(3, "Number of ITS cluster step not passed\n");
 712                                                 continue;
 713                                         }
 714                                 }
 715                                 else{
 716                                         AliDebug(3, "Reco acceptance step not passed\n");
 717                                         continue;
 718                                 }
 719                         }
 720                         else {
 721                                 AliDebug(3, "Reco step not passed\n");
 722                                 continue;
 723                         }
 724                 }
 725
 726                 if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
 727         } // end loop on candidate
 728
 729         fCountReco+= icountReco;
 730         fCountRecoAcc+= icountRecoAcc;
 731         fCountRecoITSClusters+= icountRecoITSClusters;
 732         fCountRecoPPR+= icountRecoPPR;
 733         fCountRecoPID+= icountRecoPID;
 734
 735         fHistEventsProcessed->Fill(0);
 736
 737         delete[] containerInput;
 738         delete[] containerInputMC;
 739         delete cfVtxHF;
 740         if (trackCuts){
 741         //      for (Int_t i=0; i<cfVtxHF->GetNProngs(); i++){
 742         //              delete [] trackCuts[i];
 743         //      }
 744         delete [] trackCuts;
 745         }
 746 }
 747
 748 //___________________________________________________________________________
 749 void AliCFTaskVertexingHF::Terminate(Option_t*)
 750 {
 751         // The Terminate() function is the last function to be called during
 752         // a query. It always runs on the client, it can be used to present
 753         // the results graphically or save the results to file.
 754
 755         AliAnalysisTaskSE::Terminate();
 756
 757         AliInfo(Form("Found %i MC particles that are %s in MC, in %d events",fCountMC,fPartName.Data(),fEvents));
 758         AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, in %d events",fCountAcc,fPartName.Data(),fEvents));
 759         AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, and satisfy Vertex requirement in %d events",fCountVertex,fPartName.Data(),fEvents));
 760         AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, and satisfy ITS+TPC refit requirementin %d events",fCountRefit,fPartName.Data(),fEvents));
 761         AliInfo(Form("Found %i reco %s that are decaying in %s, in %d events",fCountReco,fPartName.Data(),fDauNames.Data(),fEvents));
 762         AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and are in the requested acceptance, in %d events",fCountRecoAcc,fPartName.Data(),fDauNames.Data(),fEvents));
 763         AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and have at least 5 clusters in ITS, in %d events",fCountRecoITSClusters,fPartName.Data(),fDauNames.Data(),fEvents));
 764         AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and satisfy PPR cuts, in %d events",fCountRecoPPR,fPartName.Data(),fDauNames.Data(),fEvents));
 765         AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and satisfy PPR+PID cuts, in %d events",fCountRecoPID,fPartName.Data(),fDauNames.Data(),fEvents));
 766
 767         // draw some example plots....
 768
 769         AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
 770         if(!cont) {
 771                 printf("CONTAINER NOT FOUND\n");
 772                 return;
 773         }
 774         // projecting the containers to obtain histograms
 775         // first argument = variable, second argument = step
 776
 777         TH1D* h[3][12];
 778         for(Int_t iC=0;iC<12; iC++){
 779           // MC-level
 780           h[0][iC] =   cont->ShowProjection(iC,0);
 781           // MC-Acceptance level
 782           h[1][iC] =   cont->ShowProjection(iC,1);
 783           // Reco-level
 784           h[2][iC] =   cont->ShowProjection(iC,4);
 785         }
 786         TString titles[12];
 787         if(fDecayChannel==31){
 788           titles[0]="pT_Dplus (GeV/c)";
 789           titles[1]="rapidity";
 790           titles[2]="phi (rad)";
 791           titles[3]="cT (#mum)";
 792           titles[4]="cosPointingAngle";
 793           titles[5]="pT_1 (GeV/c)";
 794           titles[6]="pT_2 (GeV/c)";
 795           titles[7]="pT_3 (GeV/c)";
 796           titles[8]="d0_1 (#mum)";
 797           titles[9]="d0_2 (#mum)";
 798           titles[10]="d0_3 (#mum)";
 799           titles[11]="zVertex (cm)";
 800         }else{
 801           titles[0]="pT_D0 (GeV/c)";
 802           titles[1]="rapidity";
 803           titles[2]="cosThetaStar";
 804           titles[3]="pT_pi (GeV/c)";
 805           titles[4]="pT_K (Gev/c)";
 806           titles[5]="cT (#mum)";
 807           titles[6]="dca (#mum)";
 808           titles[7]="d0_pi (#mum)";
 809           titles[8]="d0_K (#mum)";
 810           titles[9]="d0xd0 (#mum^2)";
 811           titles[10]="cosPointingAngle";
 812           titles[11]="phi (rad)";
 813         }
 814         Int_t markers[12]={20,24,21,25,27,28,
 815                            20,24,21,25,27,28};
 816         Int_t colors[3]={2,8,4};
 817         for(Int_t iC=0;iC<12; iC++){
 818           for(Int_t iStep=0;iStep<3;iStep++){
 819             h[iStep][iC]->SetTitle(titles[iC].Data());
 820             h[iStep][iC]->GetXaxis()->SetTitle(titles[iC].Data());
 821             Double_t maxh=h[iStep][iC]->GetMaximum();
 822             h[iStep][iC]->GetYaxis()->SetRangeUser(0,maxh*1.2);
 823             h[iStep][iC]->SetMarkerStyle(markers[iC]);
 824             h[iStep][iC]->SetMarkerColor(colors[iStep]);
 825           }
 826         }
 827
 828
 829
 830         gStyle->SetCanvasColor(0);
 831         gStyle->SetFrameFillColor(0);
 832         gStyle->SetTitleFillColor(0);
 833         gStyle->SetStatColor(0);
 834
 835         // drawing in 2 separate canvas for a matter of clearity
 836         TCanvas * c1 =new TCanvas("c1New","Vars 0,1,2",1100,1600);
 837         c1->Divide(3,3);
 838         Int_t iPad=1;
 839         for(Int_t iVar=0; iVar<3; iVar++){
 840           c1->cd(iPad++);
 841           h[0][iVar]->Draw("p");
 842           c1->cd(iPad++);
 843           h[1][iVar]->Draw("p");
 844           c1->cd(iPad++);
 845           h[2][iVar]->Draw("p");
 846         }
 847
 848         TCanvas * c2 =new TCanvas("c2New","Vars 3,4,5",1100,1600);
 849         c2->Divide(3,3);
 850         iPad=1;
 851         for(Int_t iVar=3; iVar<6; iVar++){
 852           c2->cd(iPad++);
 853           h[0][iVar]->Draw("p");
 854           c2->cd(iPad++);
 855           h[1][iVar]->Draw("p");
 856           c2->cd(iPad++);
 857           h[2][iVar]->Draw("p");
 858         }
 859
 860
 861         TCanvas * c3 =new TCanvas("c3New","Vars 6,7,8",1100,1600);
 862         c3->Divide(3,3);
 863         iPad=1;
 864         for(Int_t iVar=6; iVar<9; iVar++){
 865           c3->cd(iPad++);
 866           h[0][iVar]->Draw("p");
 867           c3->cd(iPad++);
 868           h[1][iVar]->Draw("p");
 869           c3->cd(iPad++);
 870           h[2][iVar]->Draw("p");
 871         }
 872
 873
 874         TCanvas * c4 =new TCanvas("c4New","Vars 9,10,11",1100,1600);
 875         c4->Divide(3,3);
 876         iPad=1;
 877         for(Int_t iVar=9; iVar<11; iVar++){
 878           c4->cd(iPad++);
 879           h[0][iVar]->Draw("p");
 880           c4->cd(iPad++);
 881           h[1][iVar]->Draw("p");
 882           c4->cd(iPad++);
 883           h[2][iVar]->Draw("p");
 884         }
 885
 886
 887         THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
 888
 889         TH2D* corr1 =hcorr->Projection(0,2);
 890         TH2D* corr2 = hcorr->Projection(1,3);
 891
 892         TCanvas * c7 =new TCanvas("c7New","",800,400);
 893         c7->Divide(2,1);
 894         c7->cd(1);
 895         corr1->Draw("text");
 896         c7->cd(2);
 897         corr2->Draw("text");
 898
 899
 900         TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE");
 901
 902         corr1->Write();
 903         corr2->Write();
 904         for(Int_t iC=0;iC<12; iC++){
 905           for(Int_t iStep=0;iStep<3;iStep++){
 906             h[iStep][iC]->Write(Form("Step%d_%s",iStep,titles[iC].Data()));
 907           }
 908         }
 909         file_projection->Close();
 910
 911
 912
 913         /*
 914          c1->SaveAs("Plots/pT_rapidity_cosThetaStar.eps");
 915          c2->SaveAs("Plots/pTpi_pTK_cT.eps");
 916          c3->SaveAs("Plots/dca_d0pi_d0TK.eps");
 917          c4->SaveAs("Plots/d0xd0_cosPointingAngle.eps");
 918
 919          c1->SaveAs("Plots/pT_rapidity_cosThetaStar.gif");
 920          c2->SaveAs("Plots/pTpi_pTK_cT.gif");
 921          c3->SaveAs("Plots/dca_d0pi_d0TK.gif");
 922          c4->SaveAs("Plots/d0xd0_cosPointingAngle.gif");
 923          */
 924 }
 925
 926 //___________________________________________________________________________
 927 void AliCFTaskVertexingHF::UserCreateOutputObjects()
 928 {
 929         //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
 930         //TO BE SET BEFORE THE EXECUTION OF THE TASK
 931         //
 932         Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
 933
 934         //slot #1
 935         OpenFile(1);
 936         fHistEventsProcessed = new TH1I("CFHFchist0","",1,0,1) ;
 937 }
 938
 939
 940 //_________________________________________________________________________
 941 Double_t AliCFTaskVertexingHF::GetWeight(Float_t pt)
 942 {
 943         //
 944         // calculating the weight to fill the container
 945         //
 946
 947         // FNOLL central:
 948         // p0 = 1.63297e-01 --> 0.322643
 949         // p1 = 2.96275e+00
 950         // p2 = 2.30301e+00
 951         // p3 = 2.50000e+00
 952
 953         // PYTHIA
 954         // p0 = 1.85906e-01 --> 0.36609
 955         // p1 = 1.94635e+00
 956         // p2 = 1.40463e+00
 957         // p3 = 2.50000e+00
 958
 959         Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
 960         Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
 961
 962         Double_t dndpt_func1 = dNdptFit(pt,func1);
 963         Double_t dndpt_func2 = dNdptFit(pt,func2);
 964         AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
 965         return dndpt_func1/dndpt_func2;
 966 }
 967
 968 //__________________________________________________________________________________________________
 969 Double_t AliCFTaskVertexingHF::dNdptFit(Float_t pt, Double_t* par)
 970 {
 971         //
 972         // calculating dNdpt
 973         //
 974
 975         Double_t denom =  TMath::Power((pt/par[1]), par[3] );
 976         Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
 977
 978         return dNdpt;
 979 }