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 "AliCFVertexingHF.h"
  71 #include "AliAnalysisDataSlot.h"
  72 #include "AliAnalysisDataContainer.h"
  73
  74 //__________________________________________________________________________
  75 AliCFTaskVertexingHF::AliCFTaskVertexingHF() :
  76         AliAnalysisTaskSE(),
  77         fPDG(0),
  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 {
 100         //
 101         //Default ctor
 102         //
 103 }
 104 //___________________________________________________________________________
 105 AliCFTaskVertexingHF::AliCFTaskVertexingHF(const Char_t* name, AliRDHFCuts* cuts) :
 106         AliAnalysisTaskSE(name),
 107         fPDG(0),
 108         fCFManager(0x0),
 109         fHistEventsProcessed(0x0),
 110         fCorrelation(0x0),
 111         fCountMC(0),
 112         fCountAcc(0),
 113         fCountVertex(0),
 114         fCountRefit(0),
 115         fCountReco(0),
 116         fCountRecoAcc(0),
 117         fCountRecoITSClusters(0),
 118         fCountRecoPPR(0),
 119         fCountRecoPID(0),
 120         fEvents(0),
 121         fDecayChannel(0),
 122         fFillFromGenerated(kFALSE),
 123         fOriginDselection(0),
 124         fAcceptanceUnf(kTRUE),
 125         fCuts(cuts),
 126         fUseWeight(kFALSE),
 127         fWeight(1.),
 128         fNvar(0)
 129 {
 130         //
 131         // Constructor. Initialization of Inputs and Outputs
 132         //
 133         /*
 134           DefineInput(0) and DefineOutput(0)
 135           are taken care of by AliAnalysisTaskSE constructor
 136         */
 137         DefineOutput(1,TH1I::Class());
 138         DefineOutput(2,AliCFContainer::Class());
 139         DefineOutput(3,THnSparseD::Class());
 140         DefineOutput(4,AliRDHFCuts::Class());
 141
 142         fCuts->PrintAll();
 143 }
 144
 145 //___________________________________________________________________________
 146 AliCFTaskVertexingHF& AliCFTaskVertexingHF::operator=(const AliCFTaskVertexingHF& c)
 147 {
 148         //
 149         // Assignment operator
 150         //
 151         if (this!=&c) {
 152                 AliAnalysisTaskSE::operator=(c) ;
 153                 fCFManager  = c.fCFManager;
 154                 fHistEventsProcessed = c.fHistEventsProcessed;
 155                 fCuts = c.fCuts;
 156         }
 157         return *this;
 158 }
 159
 160 //___________________________________________________________________________
 161 AliCFTaskVertexingHF::AliCFTaskVertexingHF(const AliCFTaskVertexingHF& c) :
 162         AliAnalysisTaskSE(c),
 163         fPDG(c.fPDG),
 164         fCFManager(c.fCFManager),
 165         fHistEventsProcessed(c.fHistEventsProcessed),
 166         fCorrelation(c.fCorrelation),
 167         fCountMC(c.fCountMC),
 168         fCountAcc(c.fCountAcc),
 169         fCountVertex(c.fCountVertex),
 170         fCountRefit(c.fCountRefit),
 171         fCountReco(c.fCountReco),
 172         fCountRecoAcc(c.fCountRecoAcc),
 173         fCountRecoITSClusters(c.fCountRecoITSClusters),
 174         fCountRecoPPR(c.fCountRecoPPR),
 175         fCountRecoPID(c.fCountRecoPID),
 176         fEvents(c.fEvents),
 177         fDecayChannel(c.fDecayChannel),
 178         fFillFromGenerated(c.fFillFromGenerated),
 179         fOriginDselection(c.fOriginDselection),
 180         fAcceptanceUnf(c.fAcceptanceUnf),
 181         fCuts(c.fCuts),
 182         fUseWeight(c.fUseWeight),
 183         fWeight(c.fWeight),
 184         fNvar(c.fNvar)
 185 {
 186         //
 187         // Copy Constructor
 188         //
 189 }
 190
 191 //___________________________________________________________________________
 192 AliCFTaskVertexingHF::~AliCFTaskVertexingHF()
 193 {
 194         //
 195         //destructor
 196         //
 197         if (fCFManager)           delete fCFManager ;
 198         if (fHistEventsProcessed) delete fHistEventsProcessed ;
 199         if (fCorrelation)         delete fCorrelation ;
 200         if (fCuts)                delete fCuts;
 201 }
 202
 203 //_________________________________________________________________________-
 204 void AliCFTaskVertexingHF::Init()
 205 {
 206         //
 207         // Initialization
 208         //
 209
 210         if (fDebug>1) printf("AliCFTaskVertexingHF::Init()");
 211         AliRDHFCuts *copyfCuts = 0x0;
 212
 213
 214         switch (fDecayChannel){
 215         case 2:{
 216                 copyfCuts = new AliRDHFCutsD0toKpi(*(dynamic_cast<AliRDHFCutsD0toKpi*>(fCuts)));
 217                 fNvar = 13;
 218                 break;
 219         }
 220         case 21:{
 221                 copyfCuts = new AliRDHFCutsDStartoKpipi(*(dynamic_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
 222                 fNvar = 13;
 223                 break;
 224         }
 225         case 31:{
 226                 copyfCuts = new AliRDHFCutsDplustoKpipi(*(dynamic_cast<AliRDHFCutsDplustoKpipi*>(fCuts)));
 227                 fNvar = 12;
 228                 break;
 229         }
 230         case 32:{
 231                 copyfCuts = new AliRDHFCutsLctopKpi(*(dynamic_cast<AliRDHFCutsLctopKpi*>(fCuts)));
 232                 fNvar = 13;
 233                 break;
 234         }
 235         case 33:{
 236                 copyfCuts = new AliRDHFCutsDstoKKpi(*(dynamic_cast<AliRDHFCutsDstoKKpi*>(fCuts)));
 237                 fNvar = 13;
 238                 break;
 239         }
 240         case 4:{
 241                 copyfCuts = new AliRDHFCutsD0toKpipipi(*(dynamic_cast<AliRDHFCutsD0toKpipipi*>(fCuts)));
 242                 fNvar = 13;
 243                 break;
 244         }
 245         default:
 246                 AliFatal("The decay channel MUST be defined according to AliCFVertexing::DecayChannel - Exiting...");
 247                 break;
 248         }
 249
 250         const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
 251         copyfCuts->SetName(nameoutput);
 252
 253         //Post the data
 254         PostData(4, copyfCuts);
 255
 256         return;
 257 }
 258
 259 //_________________________________________________
 260 void AliCFTaskVertexingHF::UserExec(Option_t *)
 261 {
 262         //
 263         // Main loop function
 264         //
 265
 266         PostData(1,fHistEventsProcessed) ;
 267         PostData(2,fCFManager->GetParticleContainer()) ;
 268         PostData(3,fCorrelation) ;
 269
 270         AliESDtrackCuts* trackCuts = fCuts->GetTrackCuts();
 271
 272         if (fFillFromGenerated){
 273                 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
 274         }
 275
 276         if (!fInputEvent) {
 277                 Error("UserExec","NO EVENT FOUND!");
 278                 return;
 279         }
 280
 281         AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 282
 283         TClonesArray *arrayBranch=0;
 284
 285         if(!aodEvent && AODEvent() && IsStandardAOD()) {
 286                 // In case there is an AOD handler writing a standard AOD, use the AOD
 287                 // event in memory rather than the input (ESD) event.
 288                 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
 289                 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 290                 // have to taken from the AOD event hold by the AliAODExtension
 291                 AliAODHandler* aodHandler = (AliAODHandler*)
 292                         ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 293                 if(aodHandler->GetExtensions()) {
 294                         AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 295                         AliAODEvent *aodFromExt = ext->GetAOD();
 296
 297                         switch (fDecayChannel){
 298                         case 2:{
 299                                 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
 300                                 break;
 301                         }
 302                         case 21:{
 303                                 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
 304                                 break;
 305                         }
 306                         case 31:
 307                         case 32:
 308                         case 33:{
 309                                 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 310                                 break;
 311                         }
 312                         case 4:{
 313                                 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
 314                                 break;
 315                         }
 316                         default:
 317                                 break;
 318                         }
 319                 }
 320         }
 321         else {
 322                 switch (fDecayChannel){
 323                 case 2:{
 324                         arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
 325                         break;
 326                 }
 327                 case 21:{
 328                         arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
 329                         break;
 330                 }
 331                 case 31:
 332                 case 32:
 333                 case 33:{
 334                         arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
 335                         break;
 336                 }
 337                 case 4:{
 338                         arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm4Prong");
 339                         break;
 340                 }
 341                 default:
 342                         break;
 343                 }
 344         }
 345
 346         AliAODVertex *aodVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
 347         if (!aodVtx) return;
 348
 349         if (!arrayBranch) {
 350                 AliError("Could not find array of HF vertices");
 351                 return;
 352         }
 353
 354         fEvents++;
 355         if (fEvents%10000 == 0) AliDebug(2,Form("Event %d",fEvents));
 356
 357         fCFManager->SetRecEventInfo(aodEvent);
 358         fCFManager->SetMCEventInfo(aodEvent);
 359
 360         //******** DEFINE number of variables of the container***** for now set at 13, in the future in the config macro.
 361         //      Int_t nVar = 13;
 362
 363         Double_t* containerInput = new Double_t[fNvar];
 364         Double_t* containerInputMC = new Double_t[fNvar];
 365
 366         //loop on the MC event
 367
 368         TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 369         if (!mcArray) {
 370                 AliError("Could not find Monte-Carlo in AOD");
 371                 return;
 372         }
 373         Int_t icountMC = 0;
 374         Int_t icountAcc = 0;
 375         Int_t icountReco = 0;
 376         Int_t icountVertex = 0;
 377         Int_t icountRefit = 0;
 378         Int_t icountRecoAcc = 0;
 379         Int_t icountRecoITSClusters = 0;
 380         Int_t icountRecoPPR = 0;
 381         Int_t icountRecoPID = 0;
 382         Int_t cquarks = 0;
 383         UShort_t originDselection = 0;
 384
 385         AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 386         if (!mcHeader) {
 387                 AliError("Could not find MC Header in AOD");
 388                 return;
 389         }
 390
 391         AliCFVertexingHF* cfVtxHF=0x0;
 392         switch (fDecayChannel){
 393         case 2:{
 394                 cfVtxHF = new AliCFVertexingHF2Prong(mcArray, originDselection);
 395                 break;
 396         }
 397         case 21:{
 398                 //        cfVtxHF = new AliCFVertexingHFCascade(mcArray, originDselection);  // not there yet
 399                 break;
 400         }
 401         case 31:
 402         case 32:
 403         case 33:{
 404                //cfVtxHF = new AliCFVertexingHF3Prong(mcArray, originDselection, fDecayChannel);
 405                 break;
 406         }
 407         case 4:{
 408                 //cfVtxHF = new AliCFVertexingHF4Prong(mcArray, originDselection);  // not there yet
 409                 break;
 410         }
 411         default:
 412                 break;
 413         }
 414
 415         Double_t zPrimVertex = aodVtx ->GetZ();
 416         Double_t zMCVertex = mcHeader->GetVtxZ();
 417
 418         //General settings: vertex, feed down and fill reco container with generated values.
 419         cfVtxHF->SetRecoPrimVertex(zPrimVertex);
 420         cfVtxHF->SetMCPrimaryVertex(zMCVertex);
 421         cfVtxHF->SetFillFromGenerated(fFillFromGenerated);
 422
 423         for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
 424
 425                 AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
 426
 427                 // check the MC-level cuts, must be the desidered particle
 428                 if (!fCFManager->CheckParticleCuts(0, mcPart)) continue;  // 0 stands for MC level
 429
 430                 cfVtxHF->SetMCCandidateParam(iPart);
 431                 cfVtxHF->SetNVar(fNvar);
 432                 //counting c quarks
 433                 cquarks += cfVtxHF->MCcquarkCounting(mcPart);
 434
 435                 //check the candiate family at MC level
 436                 if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) {
 437                         Printf("Check on the family wrong!!! (decaychannel = %d)",fDecayChannel);
 438                         continue;
 439                 }
 440                 else{
 441                         Printf("Check on the family OK!!! (decaychannel = %d)",fDecayChannel);
 442                 }
 443
 444                 //Fill the MC container
 445                 Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC);
 446
 447                 if (mcContainerFilled) {
 448                         if (fUseWeight)fWeight = GetWeight(containerInputMC[0]);
 449                         if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) continue;
 450                         //MC Limited Acceptance
 451                         if (TMath::Abs(containerInputMC[1]) < 0.5) {
 452                                 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc, fWeight);
 453                                 AliDebug(3,"MC Lim Acc container filled\n");
 454                         }
 455
 456                         //MC
 457                         fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated, fWeight);
 458                         icountMC++;
 459                         AliDebug(3,"MC cointainer filled \n");
 460
 461                         // MC in acceptance
 462                         // check the MC-Acceptance level cuts
 463                         // since standard CF functions are not applicable, using Kine Cuts on daughters
 464                         Bool_t mcAccepStep = cfVtxHF-> MCAcceptanceStep();
 465                         if (mcAccepStep){
 466                                 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance, fWeight);
 467                                 AliDebug(3,"MC acceptance cut passed\n");
 468                                 icountAcc++;
 469
 470                                 //MC Vertex step
 471                                 if (fCuts->IsEventSelected(aodEvent)){
 472                                         // filling the container if the vertex is ok
 473                                         fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex, fWeight) ;
 474                                         AliDebug(3,"Vertex cut passed and container filled\n");
 475                                         icountVertex++;
 476
 477                                         //mc Refit requirement
 478                                         Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, trackCuts);
 479                                         if (mcRefitStep){
 480                                                 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit, fWeight);
 481                                                 AliDebug(3,"MC Refit cut passed and container filled\n");
 482                                                 icountRefit++;
 483                                         }
 484                                         else{
 485                                                 AliDebug(3,"MC Refit cut not passed\n");
 486                                                 continue;
 487                                         }
 488                                 }
 489                                 else{
 490                                   AliDebug (3, "MC vertex step not passed\n");
 491                                   continue;
 492                                 }
 493                         }
 494                         else{
 495                                 AliDebug (3, "MC in acceptance step not passed\n");
 496                                 continue;
 497                         }
 498                 }
 499                 else {
 500                         AliDebug (3, "MC container not filled\n");
 501                 }
 502         }
 503
 504         if (cquarks<2) AliDebug(2,Form("Event with %d c-quarks", cquarks));
 505         AliDebug(2,Form("Found %i MC particles that are D0!!",icountMC));
 506         AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Acc cuts!!",icountAcc));
 507         AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Vertex cuts!!",icountVertex));
 508         AliDebug(2,Form("Found %i MC particles that are D0 and satisfy Refit cuts!!",icountRefit));
 509
 510         // Now go to rec level
 511         fCountMC += icountMC;
 512         fCountAcc += icountAcc;
 513         fCountVertex+= icountVertex;
 514         fCountRefit+= icountRefit;
 515
 516         AliDebug(2,Form("Found %d vertices",arrayBranch->GetEntriesFast()));
 517         AliInfo(Form("Found %d vertices for decay channel %d",arrayBranch->GetEntriesFast(),fDecayChannel));
 518
 519         for(Int_t iCandid = 0; iCandid<arrayBranch->GetEntriesFast();iCandid++){
 520                 Printf("iCandid = %d", iCandid);
 521
 522                 AliAODRecoDecayHF* charmCandidate=0x0;
 523                 switch (fDecayChannel){
 524                 case 2:{
 525                         charmCandidate = (AliAODRecoDecayHF2Prong*)arrayBranch->At(iCandid);
 526                         break;
 527                 }
 528                 case 21:{
 529                         charmCandidate = (AliAODRecoCascadeHF*)arrayBranch->At(iCandid);
 530                         break;
 531                 }
 532                 case 31:
 533                 case 32:
 534                 case 33:{
 535                         charmCandidate = (AliAODRecoDecayHF3Prong*)arrayBranch->At(iCandid);
 536                         break;
 537                 }
 538                 case 4:{
 539                         charmCandidate = (AliAODRecoDecayHF4Prong*)arrayBranch->At(iCandid);
 540                         break;
 541                 }
 542                 default:
 543                         break;
 544                 }
 545
 546                 Bool_t unsetvtx=kFALSE;
 547                 if(!charmCandidate->GetOwnPrimaryVtx()) {
 548                         charmCandidate->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables
 549                         unsetvtx=kTRUE;
 550                 }
 551
 552                 Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)charmCandidate);
 553                 if (!signAssociation){
 554                         charmCandidate = 0x0;
 555                         continue;
 556                 }
 557
 558                 Int_t isPartOrAntipart = cfVtxHF->CheckReflexion();
 559
 560                 Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
 561                 if (recoContFilled){
 562
 563                         if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue;
 564
 565                         //Reco Step
 566                         Bool_t recoStep = cfVtxHF->RecoStep();
 567                         Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
 568
 569                         if (recoStep && recoContFilled && vtxCheck){
 570                                 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ;
 571                                 icountReco++;
 572                                 AliDebug(3,"Reco step  passed and container filled\n");
 573
 574                                 //Reco in the acceptance -- take care of UNFOLDING!!!!
 575                                 Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(trackCuts);
 576                                 if (recoAcceptanceStep) {
 577                                         fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ;
 578                                         icountRecoAcc++;
 579                                         AliDebug(3,"Reco acceptance cut passed and container filled\n");
 580
 581                                         if(fAcceptanceUnf){
 582                                                 Double_t fill[4]; //fill response matrix
 583                                                 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
 584                                                 if (bUnfolding) fCorrelation->Fill(fill);
 585                                         }
 586
 587                                         //Number of ITS cluster requirements
 588                                         Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks);
 589                                         if (recoITSnCluster){
 590                                                 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ;
 591                                                 icountRecoITSClusters++;
 592                                                 AliDebug(3,"Reco n ITS cluster cut passed and container filled\n");
 593
 594                                                 Bool_t iscutsusingpid = fCuts->GetIsUsePID();
 595                                                 Int_t recoAnalysisCuts = -1, recoPidSelection = -1;
 596                                                 fCuts->SetUsePID(kFALSE);
 597                                                 recoAnalysisCuts = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
 598
 599                                                 if (recoAnalysisCuts > 3){ // Ds case, where more possibilities are considered
 600                                                         if (recoAnalysisCuts >= 8){
 601                                                                 recoAnalysisCuts -= 8; // removing K0star mass
 602                                                         }
 603                                                         if (recoAnalysisCuts >= 4){
 604                                                                 recoAnalysisCuts -= 4; // removing Phi mass
 605                                                         }
 606                                                 }
 607
 608                                                 fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object
 609                                                 if (recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart){
 610                                                         fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR, fWeight);
 611                                                         icountRecoPPR++;
 612                                                         AliDebug(3,"Reco Analysis cuts passed and container filled \n");
 613                                                         //pid selection
 614                                                         //recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kPID);
 615                                                         //if((fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==isPartOrAntipart)||(fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==3)){
 616                                                         recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
 617                                                         if (recoPidSelection == 3 || recoPidSelection == isPartOrAntipart){
 618                                                                 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID, fWeight);
 619                                                                 icountRecoPID++;
 620                                                                 AliDebug(3,"Reco PID cuts passed and container filled \n");
 621                                                                 if(!fAcceptanceUnf){
 622                                                                         Double_t fill[4]; //fill response matrix
 623                                                                         Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
 624                                                                         if (bUnfolding) fCorrelation->Fill(fill);
 625                                                                 }
 626                                                         }
 627                                                         else {
 628                                                                 AliDebug(3, "Analysis Cuts step not passed \n");
 629                                                                 continue;
 630                                                         }
 631                                                 }
 632                                                 else {
 633                                                         AliDebug(3, "PID selection not passed \n");
 634                                                         continue;
 635                                                 }
 636                                         }
 637                                         else{
 638                                                 AliDebug(3, "Number of ITS cluster step not passed\n");
 639                                                 continue;
 640                                         }
 641                                 }
 642                                 else{
 643                                         AliDebug(3, "Reco acceptance step not passed\n");
 644                                         continue;
 645                                 }
 646                         }
 647                         else {
 648                                 AliDebug(3, "Reco step not passed\n");
 649                                 continue;
 650                         }
 651                 }
 652
 653                 if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
 654         } // end loop on candidate
 655
 656         fCountReco+= icountReco;
 657         fCountRecoAcc+= icountRecoAcc;
 658         fCountRecoITSClusters+= icountRecoITSClusters;
 659         fCountRecoPPR+= icountRecoPPR;
 660         fCountRecoPID+= icountRecoPID;
 661
 662         fHistEventsProcessed->Fill(0);
 663
 664         delete[] containerInput;
 665         containerInput = 0x0;
 666         delete[] containerInputMC;
 667         containerInputMC = 0x0;
 668         delete cfVtxHF;
 669
 670 }
 671
 672 //___________________________________________________________________________
 673 void AliCFTaskVertexingHF::Terminate(Option_t*)
 674 {
 675         // The Terminate() function is the last function to be called during
 676         // a query. It always runs on the client, it can be used to present
 677         // the results graphically or save the results to file.
 678
 679         AliAnalysisTaskSE::Terminate();
 680
 681         AliInfo(Form("Found %i MC particles that are D0 in MC, in %d events",fCountMC,fEvents));
 682         AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, in %d events",fCountAcc,fEvents));
 683         AliInfo(Form("Found %i MC particles that are D0 in MC and satisfy Acc cuts, and satisfy Vertex requirement in %d events",fCountVertex,fEvents));
 684         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));
 685         AliInfo(Form("Found %i reco D0 that are decaying in K+pi, in %d events",fCountReco,fEvents));
 686         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));
 687         AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and have at least 5 clusters in ITS, in %d events",fCountRecoITSClusters,fEvents));
 688         AliInfo(Form("Among the above, found %i reco D0 that are decaying in K+pi and satisfy PPR cuts, in %d events",fCountRecoPPR,fEvents));
 689
 690         // draw some example plots....
 691
 692         AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
 693         if(!cont) {
 694                 printf("CONTAINER NOT FOUND\n");
 695                 return;
 696         }
 697         // projecting the containers to obtain histograms
 698         // first argument = variable, second argument = step
 699
 700         // MC-level
 701         TH1D* h00 =   cont->ShowProjection(0,0) ;   // pt
 702         TH1D* h10 =   cont->ShowProjection(1,0) ;   // rapidity
 703         TH1D* h20 =   cont->ShowProjection(2,0) ;   // cosThetaStar
 704         TH1D* h30 =   cont->ShowProjection(3,0) ;   // pTpi
 705         TH1D* h40 =   cont->ShowProjection(4,0) ;   // pTK
 706         TH1D* h50 =   cont->ShowProjection(5,0) ;   // cT
 707         TH1D* h60 =   cont->ShowProjection(6,0) ;   // dca
 708         TH1D* h70 =   cont->ShowProjection(7,0) ;   // d0pi
 709         TH1D* h80 =   cont->ShowProjection(8,0) ;   // d0K
 710         TH1D* h90 =   cont->ShowProjection(9,0) ;   // d0xd0
 711         TH1D* h100 =   cont->ShowProjection(10,0) ;   // cosPointingAngle
 712         TH1D* h110 =   cont->ShowProjection(11,0) ;   // phi
 713
 714         // MC-Acceptance level
 715         TH1D* h01 =   cont->ShowProjection(0,1) ;   // pt
 716         TH1D* h11 =   cont->ShowProjection(1,1) ;   // rapidity
 717         TH1D* h21 =   cont->ShowProjection(2,1) ;   // cosThetaStar
 718         TH1D* h31 =   cont->ShowProjection(3,1) ;   // pTpi
 719         TH1D* h41 =   cont->ShowProjection(4,1) ;   // pTK
 720         TH1D* h51 =   cont->ShowProjection(5,1) ;   // cT
 721         TH1D* h61 =   cont->ShowProjection(6,1) ;   // dca
 722         TH1D* h71 =   cont->ShowProjection(7,1) ;   // d0pi
 723         TH1D* h81 =   cont->ShowProjection(8,1) ;   // d0K
 724         TH1D* h91 =   cont->ShowProjection(9,1) ;   // d0xd0
 725         TH1D* h101 =   cont->ShowProjection(10,1) ;   // cosPointingAngle
 726         TH1D* h111 =   cont->ShowProjection(11,1) ;   // phi
 727
 728         // Reco-level
 729         TH1D* h04 =   cont->ShowProjection(0,4) ;   // pt
 730         TH1D* h14 =   cont->ShowProjection(1,4) ;   // rapidity
 731         TH1D* h24 =   cont->ShowProjection(2,4) ;   // cosThetaStar
 732         TH1D* h34 =   cont->ShowProjection(3,4) ;   // pTpi
 733         TH1D* h44 =   cont->ShowProjection(4,4) ;   // pTK
 734         TH1D* h54 =   cont->ShowProjection(5,4) ;   // cT
 735         TH1D* h64 =   cont->ShowProjection(6,4) ;   // dca
 736         TH1D* h74 =   cont->ShowProjection(7,4) ;   // d0pi
 737         TH1D* h84 =   cont->ShowProjection(8,4) ;   // d0K
 738         TH1D* h94 =   cont->ShowProjection(9,4) ;   // d0xd0
 739         TH1D* h104 =   cont->ShowProjection(10,4) ;   // cosPointingAngle
 740         TH1D* h114 =   cont->ShowProjection(11,4) ;   // phi
 741
 742         h00->SetTitle("pT_D0 (GeV/c)");
 743         h10->SetTitle("rapidity");
 744         h20->SetTitle("cosThetaStar");
 745         h30->SetTitle("pT_pi (GeV/c)");
 746         h40->SetTitle("pT_K (Gev/c)");
 747         h50->SetTitle("cT (#mum)");
 748         h60->SetTitle("dca (#mum)");
 749         h70->SetTitle("d0_pi (#mum)");
 750         h80->SetTitle("d0_K (#mum)");
 751         h90->SetTitle("d0xd0 (#mum^2)");
 752         h100->SetTitle("cosPointingAngle");
 753         h100->SetTitle("phi (rad)");
 754
 755         h00->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
 756         h10->GetXaxis()->SetTitle("rapidity");
 757         h20->GetXaxis()->SetTitle("cosThetaStar");
 758         h30->GetXaxis()->SetTitle("pT_pi (GeV/c)");
 759         h40->GetXaxis()->SetTitle("pT_K (Gev/c)");
 760         h50->GetXaxis()->SetTitle("cT (#mum)");
 761         h60->GetXaxis()->SetTitle("dca (#mum)");
 762         h70->GetXaxis()->SetTitle("d0_pi (#mum)");
 763         h80->GetXaxis()->SetTitle("d0_K (#mum)");
 764         h90->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
 765         h100->GetXaxis()->SetTitle("cosPointingAngle");
 766         h110->GetXaxis()->SetTitle("phi (rad)");
 767
 768         h01->SetTitle("pT_D0 (GeV/c)");
 769         h11->SetTitle("rapidity");
 770         h21->SetTitle("cosThetaStar");
 771         h31->SetTitle("pT_pi (GeV/c)");
 772         h41->SetTitle("pT_K (Gev/c)");
 773         h51->SetTitle("cT (#mum)");
 774         h61->SetTitle("dca (#mum)");
 775         h71->SetTitle("d0_pi (#mum)");
 776         h81->SetTitle("d0_K (#mum)");
 777         h91->SetTitle("d0xd0 (#mum^2)");
 778         h101->SetTitle("cosPointingAngle");
 779         h111->GetXaxis()->SetTitle("phi (rad)");
 780
 781         h01->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
 782         h11->GetXaxis()->SetTitle("rapidity");
 783         h21->GetXaxis()->SetTitle("cosThetaStar");
 784         h31->GetXaxis()->SetTitle("pT_pi (GeV/c)");
 785         h41->GetXaxis()->SetTitle("pT_K (Gev/c)");
 786         h51->GetXaxis()->SetTitle("cT (#mum)");
 787         h61->GetXaxis()->SetTitle("dca (#mum)");
 788         h71->GetXaxis()->SetTitle("d0_pi (#mum)");
 789         h81->GetXaxis()->SetTitle("d0_K (#mum)");
 790         h91->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
 791         h101->GetXaxis()->SetTitle("cosPointingAngle");
 792         h111->GetXaxis()->SetTitle("phi (rad)");
 793
 794         h04->SetTitle("pT_D0 (GeV/c)");
 795         h14->SetTitle("rapidity");
 796         h24->SetTitle("cosThetaStar");
 797         h34->SetTitle("pT_pi (GeV/c)");
 798         h44->SetTitle("pT_K (Gev/c)");
 799         h54->SetTitle("cT (#mum)");
 800         h64->SetTitle("dca (#mum)");
 801         h74->SetTitle("d0_pi (#mum)");
 802         h84->SetTitle("d0_K (#mum)");
 803         h94->SetTitle("d0xd0 (#mum^2)");
 804         h104->SetTitle("cosPointingAngle");
 805         h114->GetXaxis()->SetTitle("phi (rad)");
 806
 807         h04->GetXaxis()->SetTitle("pT_D0 (GeV/c)");
 808         h14->GetXaxis()->SetTitle("rapidity");
 809         h24->GetXaxis()->SetTitle("cosThetaStar");
 810         h34->GetXaxis()->SetTitle("pT_pi (GeV/c)");
 811         h44->GetXaxis()->SetTitle("pT_K (Gev/c)");
 812         h54->GetXaxis()->SetTitle("cT (#mum)");
 813         h64->GetXaxis()->SetTitle("dca (#mum)");
 814         h74->GetXaxis()->SetTitle("d0_pi (#mum)");
 815         h84->GetXaxis()->SetTitle("d0_K (#mum)");
 816         h94->GetXaxis()->SetTitle("d0xd0 (#mum^2)");
 817         h104->GetXaxis()->SetTitle("cosPointingAngle");
 818         h114->GetXaxis()->SetTitle("phi (rad)");
 819
 820         Double_t max0 = h00->GetMaximum();
 821         Double_t max1 = h10->GetMaximum();
 822         Double_t max2 = h20->GetMaximum();
 823         Double_t max3 = h30->GetMaximum();
 824         Double_t max4 = h40->GetMaximum();
 825         Double_t max5 = h50->GetMaximum();
 826         Double_t max6 = h60->GetMaximum();
 827         Double_t max7 = h70->GetMaximum();
 828         Double_t max8 = h80->GetMaximum();
 829         Double_t max9 = h90->GetMaximum();
 830         Double_t max10 = h100->GetMaximum();
 831         Double_t max11 = h110->GetMaximum();
 832
 833         h00->GetYaxis()->SetRangeUser(0,max0*1.2);
 834         h10->GetYaxis()->SetRangeUser(0,max1*1.2);
 835         h20->GetYaxis()->SetRangeUser(0,max2*1.2);
 836         h30->GetYaxis()->SetRangeUser(0,max3*1.2);
 837         h40->GetYaxis()->SetRangeUser(0,max4*1.2);
 838         h50->GetYaxis()->SetRangeUser(0,max5*1.2);
 839         h60->GetYaxis()->SetRangeUser(0,max6*1.2);
 840         h70->GetYaxis()->SetRangeUser(0,max7*1.2);
 841         h80->GetYaxis()->SetRangeUser(0,max8*1.2);
 842         h90->GetYaxis()->SetRangeUser(0,max9*1.2);
 843         h100->GetYaxis()->SetRangeUser(0,max10*1.2);
 844         h110->GetYaxis()->SetRangeUser(0,max11*1.2);
 845
 846         h01->GetYaxis()->SetRangeUser(0,max0*1.2);
 847         h11->GetYaxis()->SetRangeUser(0,max1*1.2);
 848         h21->GetYaxis()->SetRangeUser(0,max2*1.2);
 849         h31->GetYaxis()->SetRangeUser(0,max3*1.2);
 850         h41->GetYaxis()->SetRangeUser(0,max4*1.2);
 851         h51->GetYaxis()->SetRangeUser(0,max5*1.2);
 852         h61->GetYaxis()->SetRangeUser(0,max6*1.2);
 853         h71->GetYaxis()->SetRangeUser(0,max7*1.2);
 854         h81->GetYaxis()->SetRangeUser(0,max8*1.2);
 855         h91->GetYaxis()->SetRangeUser(0,max9*1.2);
 856         h101->GetYaxis()->SetRangeUser(0,max10*1.2);
 857         h111->GetYaxis()->SetRangeUser(0,max11*1.2);
 858
 859         h00->SetMarkerStyle(20);
 860         h10->SetMarkerStyle(24);
 861         h20->SetMarkerStyle(21);
 862         h30->SetMarkerStyle(25);
 863         h40->SetMarkerStyle(27);
 864         h50->SetMarkerStyle(28);
 865         h60->SetMarkerStyle(20);
 866         h70->SetMarkerStyle(24);
 867         h80->SetMarkerStyle(21);
 868         h90->SetMarkerStyle(25);
 869         h100->SetMarkerStyle(27);
 870         h110->SetMarkerStyle(28);
 871
 872         h00->SetMarkerColor(2);
 873         h10->SetMarkerColor(2);
 874         h20->SetMarkerColor(2);
 875         h30->SetMarkerColor(2);
 876         h40->SetMarkerColor(2);
 877         h50->SetMarkerColor(2);
 878         h60->SetMarkerColor(2);
 879         h70->SetMarkerColor(2);
 880         h80->SetMarkerColor(2);
 881         h90->SetMarkerColor(2);
 882         h100->SetMarkerColor(2);
 883         h110->SetMarkerColor(2);
 884
 885         h01->SetMarkerStyle(20) ;
 886         h11->SetMarkerStyle(24) ;
 887         h21->SetMarkerStyle(21) ;
 888         h31->SetMarkerStyle(25) ;
 889         h41->SetMarkerStyle(27) ;
 890         h51->SetMarkerStyle(28) ;
 891         h61->SetMarkerStyle(20);
 892         h71->SetMarkerStyle(24);
 893         h81->SetMarkerStyle(21);
 894         h91->SetMarkerStyle(25);
 895         h101->SetMarkerStyle(27);
 896         h111->SetMarkerStyle(28);
 897
 898         h01->SetMarkerColor(8);
 899         h11->SetMarkerColor(8);
 900         h21->SetMarkerColor(8);
 901         h31->SetMarkerColor(8);
 902         h41->SetMarkerColor(8);
 903         h51->SetMarkerColor(8);
 904         h61->SetMarkerColor(8);
 905         h71->SetMarkerColor(8);
 906         h81->SetMarkerColor(8);
 907         h91->SetMarkerColor(8);
 908         h101->SetMarkerColor(8);
 909         h111->SetMarkerColor(8);
 910
 911         h04->SetMarkerStyle(20) ;
 912         h14->SetMarkerStyle(24) ;
 913         h24->SetMarkerStyle(21) ;
 914         h34->SetMarkerStyle(25) ;
 915         h44->SetMarkerStyle(27) ;
 916         h54->SetMarkerStyle(28) ;
 917         h64->SetMarkerStyle(20);
 918         h74->SetMarkerStyle(24);
 919         h84->SetMarkerStyle(21);
 920         h94->SetMarkerStyle(25);
 921         h104->SetMarkerStyle(27);
 922         h114->SetMarkerStyle(28);
 923
 924         h04->SetMarkerColor(4);
 925         h14->SetMarkerColor(4);
 926         h24->SetMarkerColor(4);
 927         h34->SetMarkerColor(4);
 928         h44->SetMarkerColor(4);
 929         h54->SetMarkerColor(4);
 930         h64->SetMarkerColor(4);
 931         h74->SetMarkerColor(4);
 932         h84->SetMarkerColor(4);
 933         h94->SetMarkerColor(4);
 934         h104->SetMarkerColor(4);
 935         h114->SetMarkerColor(4);
 936
 937         gStyle->SetCanvasColor(0);
 938         gStyle->SetFrameFillColor(0);
 939         gStyle->SetTitleFillColor(0);
 940         gStyle->SetStatColor(0);
 941
 942         // drawing in 2 separate canvas for a matter of clearity
 943         TCanvas * c1 =new TCanvas("c1New","pT, rapidiy, cosThetaStar",1100,1600);
 944         c1->Divide(3,3);
 945
 946         c1->cd(1);
 947         h00->Draw("p");
 948         c1->cd(1);
 949         c1->cd(2);
 950         h01->Draw("p");
 951         c1->cd(2);
 952         c1->cd(3);
 953         h04->Draw("p");
 954         c1->cd(3);
 955         c1->cd(4);
 956         h10->Draw("p");
 957         c1->cd(4);
 958         c1->cd(5);
 959         h11->Draw("p");
 960         c1->cd(5);
 961         c1->cd(6);
 962         h14->Draw("p");
 963         c1->cd(6);
 964         c1->cd(7);
 965         h20->Draw("p");
 966         c1->cd(7);
 967         c1->cd(8);
 968         h21->Draw("p");
 969         c1->cd(8);
 970         c1->cd(9);
 971         h24->Draw("p");
 972         c1->cd(9);
 973         c1->cd();
 974
 975         TCanvas * c2 =new TCanvas("c2New","pTpi, pTK, cT",1100,1600);
 976         c2->Divide(3,3);
 977         c2->cd(1);
 978         h30->Draw("p");
 979         c2->cd(1);
 980         c2->cd(2);
 981         h31->Draw("p");
 982         c2->cd(2);
 983         c2->cd(3);
 984         h34->Draw("p");
 985         c2->cd(3);
 986         c2->cd(4);
 987         h40->Draw("p");
 988         c2->cd(4);
 989         c2->cd(5);
 990         h41->Draw("p");
 991         c2->cd(5);
 992         c2->cd(6);
 993         h44->Draw("p");
 994         c2->cd(6);
 995         c2->cd(7);
 996         h50->Draw("p");
 997         c2->cd(7);
 998         c2->cd(8);
 999         h51->Draw("p");
1000         c2->cd(8);
1001         c2->cd(9);
1002         h54->Draw("p");
1003         c2->cd(9);
1004         c2->cd();
1005
1006         TCanvas * c3 =new TCanvas("c3New","dca, d0pi, d0K",1100,1600);
1007         c3->Divide(3,3);
1008         c3->cd(1);
1009         h60->Draw("p");
1010         c3->cd(1);
1011         c3->cd(2);
1012         h61->Draw("p");
1013         c3->cd(2);
1014         c3->cd(3);
1015         h64->Draw("p");
1016         c3->cd(3);
1017         c3->cd(4);
1018         h70->Draw("p");
1019         c3->cd(4);
1020         c3->cd(5);
1021         h71->Draw("p");
1022         c3->cd(5);
1023         c3->cd(6);
1024         h74->Draw("p");
1025         c3->cd(6);
1026         c3->cd(7);
1027         h80->Draw("p");
1028         c3->cd(7);
1029         c3->cd(8);
1030         h81->Draw("p");
1031         c3->cd(8);
1032         c3->cd(9);
1033         h84->Draw("p");
1034         c3->cd(9);
1035         c3->cd();
1036
1037         TCanvas * c4 =new TCanvas("c4New","d0xd0, cosPointingAngle, phi",1100,1600);
1038         c4->Divide(3,3);
1039         c4->cd(1);
1040         h90->Draw("p");
1041         c4->cd(1);
1042         c4->cd(2);
1043         h91->Draw("p");
1044         c4->cd(2);
1045         c4->cd(3);
1046         h94->Draw("p");
1047         c4->cd(3);
1048         c4->cd(4);
1049         h100->Draw("p");
1050         c4->cd(4);
1051         c4->cd(5);
1052         h101->Draw("p");
1053         c4->cd(5);
1054         c4->cd(6);
1055         h104->Draw("p");
1056         c4->cd(6);
1057         c4->cd(7);
1058         h110->Draw("p");
1059         c4->cd(7);
1060         c4->cd(8);
1061         h111->Draw("p");
1062         c4->cd(8);
1063         c4->cd(9);
1064         h114->Draw("p");
1065         c4->cd(9);
1066         c4->cd();
1067
1068         THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
1069
1070         TH2D* corr1 =hcorr->Projection(0,2);
1071         TH2D* corr2 = hcorr->Projection(1,3);
1072
1073         TCanvas * c7 =new TCanvas("c7New","",800,400);
1074         c7->Divide(2,1);
1075         c7->cd(1);
1076         corr1->Draw("text");
1077         c7->cd(2);
1078         corr2->Draw("text");
1079
1080
1081         TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE");
1082
1083         corr1->Write();
1084         corr2->Write();
1085         h00->Write("pT_D0_step0");
1086         h10->Write("rapidity_step0");
1087         h20->Write("cosThetaStar_step0");
1088         h30->Write("pT_pi_step0");
1089         h40->Write("pT_K_step0");
1090         h50->Write("cT_step0");
1091         h60->Write("dca_step0");
1092         h70->Write("d0_pi_step0");
1093         h80->Write("d0_K_step0");
1094         h90->Write("d0xd0_step0");
1095         h100->Write("cosPointingAngle_step0");
1096         h110->Write("phi_step0");
1097
1098         h01->Write("pT_D0_step1");
1099         h11->Write("rapidity_step1");
1100         h21->Write("cosThetaStar_step1");
1101         h31->Write("pT_pi_step1");
1102         h41->Write("pT_K_step1");
1103         h51->Write("cT_step1");
1104         h61->Write("dca_step1");
1105         h71->Write("d0_pi_step1");
1106         h81->Write("d0_K_step1");
1107         h91->Write("d0xd0_step1");
1108         h101->Write("cosPointingAngle_step1");
1109         h111->Write("phi_step1");
1110
1111         h04->Write("pT_D0_step2");
1112         h14->Write("rapidity_step2");
1113         h24->Write("cosThetaStar_step2");
1114         h34->Write("pT_pi_step2");
1115         h44->Write("pT_K_step2");
1116         h54->Write("cT_step2");
1117         h64->Write("dca_step2");
1118         h74->Write("d0_pi_step2");
1119         h80->Write("d0_K_step2");
1120         h94->Write("d0xd0_step2");
1121         h104->Write("cosPointingAngle_step2");
1122         h114->Write("phi_step2");
1123
1124         file_projection->Close();
1125
1126
1127
1128         /*
1129          c1->SaveAs("Plots/pT_rapidity_cosThetaStar.eps");
1130          c2->SaveAs("Plots/pTpi_pTK_cT.eps");
1131          c3->SaveAs("Plots/dca_d0pi_d0TK.eps");
1132          c4->SaveAs("Plots/d0xd0_cosPointingAngle.eps");
1133
1134          c1->SaveAs("Plots/pT_rapidity_cosThetaStar.gif");
1135          c2->SaveAs("Plots/pTpi_pTK_cT.gif");
1136          c3->SaveAs("Plots/dca_d0pi_d0TK.gif");
1137          c4->SaveAs("Plots/d0xd0_cosPointingAngle.gif");
1138          */
1139 }
1140
1141 //___________________________________________________________________________
1142 void AliCFTaskVertexingHF::UserCreateOutputObjects()
1143 {
1144         //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1145         //TO BE SET BEFORE THE EXECUTION OF THE TASK
1146         //
1147         Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
1148
1149         //slot #1
1150         OpenFile(1);
1151         fHistEventsProcessed = new TH1I("CFHFchist0","",1,0,1) ;
1152 }
1153
1154
1155 //_________________________________________________________________________
1156 Double_t AliCFTaskVertexingHF::GetWeight(Float_t pt)
1157 {
1158         //
1159         // calculating the weight to fill the container
1160         //
1161
1162         // FNOLL central:
1163         // p0 = 1.63297e-01 --> 0.322643
1164         // p1 = 2.96275e+00
1165         // p2 = 2.30301e+00
1166         // p3 = 2.50000e+00
1167
1168         // PYTHIA
1169         // p0 = 1.85906e-01 --> 0.36609
1170         // p1 = 1.94635e+00
1171         // p2 = 1.40463e+00
1172         // p3 = 2.50000e+00
1173
1174         Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
1175         Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
1176
1177         Double_t dndpt_func1 = dNdptFit(pt,func1);
1178         Double_t dndpt_func2 = dNdptFit(pt,func2);
1179         AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
1180         return dndpt_func1/dndpt_func2;
1181 }
1182
1183 //__________________________________________________________________________________________________
1184 Double_t AliCFTaskVertexingHF::dNdptFit(Float_t pt, Double_t* par)
1185 {
1186         //
1187         // calculating dNdpt
1188         //
1189
1190         Double_t denom =  TMath::Power((pt/par[1]), par[3] );
1191         Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
1192
1193         return dNdpt;
1194 }