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