PWGHF/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 #include <TF1.h>
  39
  40 #include "AliCFTaskVertexingHF.h"
  41 #include "AliStack.h"
  42 #include "AliMCEvent.h"
  43 #include "AliCFManager.h"
  44 #include "AliCFContainer.h"
  45 #include "AliLog.h"
  46 #include "AliAnalysisManager.h"
  47 #include "AliAODHandler.h"
  48 #include "AliAODEvent.h"
  49 #include "AliAODRecoDecay.h"
  50 #include "AliAODRecoDecayHF.h"
  51 #include "AliAODRecoDecayHF2Prong.h"
  52 #include "AliAODRecoDecayHF3Prong.h"
  53 #include "AliAODRecoDecayHF4Prong.h"
  54 #include "AliAODRecoCascadeHF.h"
  55 #include "AliAODMCParticle.h"
  56 #include "AliAODMCHeader.h"
  57 #include "AliESDtrack.h"
  58 #include "TChain.h"
  59 #include "THnSparse.h"
  60 #include "TH2D.h"
  61 #include "AliESDtrackCuts.h"
  62 #include "AliRDHFCuts.h"
  63 #include "AliRDHFCutsD0toKpi.h"
  64 #include "AliRDHFCutsDplustoKpipi.h"
  65 #include "AliRDHFCutsDStartoKpipi.h"
  66 #include "AliRDHFCutsDstoKKpi.h"
  67 #include "AliRDHFCutsLctopKpi.h"
  68 #include "AliRDHFCutsD0toKpipipi.h"
  69 #include "AliCFVertexingHF2Prong.h"
  70 #include "AliCFVertexingHF3Prong.h"
  71 #include "AliCFVertexingHFCascade.h"
  72 #include "AliCFVertexingHF.h"
  73 #include "AliVertexingHFUtils.h"
  74 #include "AliAnalysisDataSlot.h"
  75 #include "AliAnalysisDataContainer.h"
  76
  77 //__________________________________________________________________________
  78 AliCFTaskVertexingHF::AliCFTaskVertexingHF() :
  79         AliAnalysisTaskSE(),
  80         fCFManager(0x0),
  81         fHistEventsProcessed(0x0),
  82         fCorrelation(0x0),
  83         fCountMC(0),
  84         fCountAcc(0),
  85         fCountVertex(0),
  86         fCountRefit(0),
  87         fCountReco(0),
  88         fCountRecoAcc(0),
  89         fCountRecoITSClusters(0),
  90         fCountRecoPPR(0),
  91         fCountRecoPID(0),
  92         fEvents(0),
  93         fDecayChannel(0),
  94         fFillFromGenerated(kFALSE),
  95         fOriginDselection(0),
  96         fAcceptanceUnf(kTRUE),
  97         fCuts(0),
  98         fUseWeight(kFALSE),
  99         fWeight(1.),
 100         fUseFlatPtWeight(kFALSE),
 101         fUseZWeight(kFALSE),
 102         fUseNchWeight(kFALSE),
 103         fNvar(0),
 104         fPartName(""),
 105         fDauNames(""),
 106         fSign(2),
 107         fCentralitySelection(kTRUE),
 108         fFakeSelection(0),
 109         fRejectIfNoQuark(kTRUE),
 110         fUseMCVertex(kFALSE),
 111         fDsOption(1),
 112         fGenDsOption(3),
 113         fConfiguration(kCheetah), // by default, setting the fast configuration
 114         fFuncWeight(0x0),
 115         fHistoMeasNch(0x0),
 116         fHistoMCNch(0x0)
 117 {
 118         //
 119         //Default ctor
 120         //
 121 }
 122 //___________________________________________________________________________
 123 AliCFTaskVertexingHF::AliCFTaskVertexingHF(const Char_t* name, AliRDHFCuts* cuts, TF1* func) :
 124         AliAnalysisTaskSE(name),
 125         fCFManager(0x0),
 126         fHistEventsProcessed(0x0),
 127         fCorrelation(0x0),
 128         fCountMC(0),
 129         fCountAcc(0),
 130         fCountVertex(0),
 131         fCountRefit(0),
 132         fCountReco(0),
 133         fCountRecoAcc(0),
 134         fCountRecoITSClusters(0),
 135         fCountRecoPPR(0),
 136         fCountRecoPID(0),
 137         fEvents(0),
 138         fDecayChannel(0),
 139         fFillFromGenerated(kFALSE),
 140         fOriginDselection(0),
 141         fAcceptanceUnf(kTRUE),
 142         fCuts(cuts),
 143         fUseWeight(kFALSE),
 144         fWeight(1.),
 145         fUseFlatPtWeight(kFALSE),
 146         fUseZWeight(kFALSE),
 147         fUseNchWeight(kFALSE),
 148         fNvar(0),
 149         fPartName(""),
 150         fDauNames(""),
 151         fSign(2),
 152         fCentralitySelection(kTRUE),
 153         fFakeSelection(0),
 154         fRejectIfNoQuark(kTRUE),
 155         fUseMCVertex(kFALSE),
 156         fDsOption(1),
 157         fGenDsOption(3),
 158         fConfiguration(kCheetah),  // by default, setting the fast configuration
 159         fFuncWeight(func),
 160         fHistoMeasNch(0x0),
 161         fHistoMCNch(0x0)
 162 {
 163         //
 164         // Constructor. Initialization of Inputs and Outputs
 165         //
 166         /*
 167           DefineInput(0) and DefineOutput(0)
 168           are taken care of by AliAnalysisTaskSE constructor
 169         */
 170         DefineOutput(1,TH1I::Class());
 171         DefineOutput(2,AliCFContainer::Class());
 172         DefineOutput(3,THnSparseD::Class());
 173         DefineOutput(4,AliRDHFCuts::Class());
 174
 175         fCuts->PrintAll();
 176 }
 177
 178 //___________________________________________________________________________
 179 AliCFTaskVertexingHF& AliCFTaskVertexingHF::operator=(const AliCFTaskVertexingHF& c)
 180 {
 181         //
 182         // Assignment operator
 183         //
 184         if (this!=&c) {
 185                 AliAnalysisTaskSE::operator=(c) ;
 186                 fCFManager  = c.fCFManager;
 187                 fHistEventsProcessed = c.fHistEventsProcessed;
 188                 fCuts = c.fCuts;
 189                 fFuncWeight = c.fFuncWeight;
 190                 fHistoMeasNch = c.fHistoMeasNch;
 191                 fHistoMCNch = c.fHistoMCNch;
 192         }
 193         return *this;
 194 }
 195
 196 //___________________________________________________________________________
 197 AliCFTaskVertexingHF::AliCFTaskVertexingHF(const AliCFTaskVertexingHF& c) :
 198         AliAnalysisTaskSE(c),
 199         fCFManager(c.fCFManager),
 200         fHistEventsProcessed(c.fHistEventsProcessed),
 201         fCorrelation(c.fCorrelation),
 202         fCountMC(c.fCountMC),
 203         fCountAcc(c.fCountAcc),
 204         fCountVertex(c.fCountVertex),
 205         fCountRefit(c.fCountRefit),
 206         fCountReco(c.fCountReco),
 207         fCountRecoAcc(c.fCountRecoAcc),
 208         fCountRecoITSClusters(c.fCountRecoITSClusters),
 209         fCountRecoPPR(c.fCountRecoPPR),
 210         fCountRecoPID(c.fCountRecoPID),
 211         fEvents(c.fEvents),
 212         fDecayChannel(c.fDecayChannel),
 213         fFillFromGenerated(c.fFillFromGenerated),
 214         fOriginDselection(c.fOriginDselection),
 215         fAcceptanceUnf(c.fAcceptanceUnf),
 216         fCuts(c.fCuts),
 217         fUseWeight(c.fUseWeight),
 218         fWeight(c.fWeight),
 219         fUseFlatPtWeight(c.fUseFlatPtWeight),
 220         fUseZWeight(c.fUseZWeight),
 221         fUseNchWeight(c.fUseNchWeight),
 222         fNvar(c.fNvar),
 223         fPartName(c.fPartName),
 224         fDauNames(c.fDauNames),
 225         fSign(c.fSign),
 226         fCentralitySelection(c.fCentralitySelection),
 227         fFakeSelection(c.fFakeSelection),
 228         fRejectIfNoQuark(c.fRejectIfNoQuark),
 229         fUseMCVertex(c.fUseMCVertex),
 230         fDsOption(c.fDsOption),
 231         fGenDsOption(c.fGenDsOption),
 232         fConfiguration(c.fConfiguration),
 233         fFuncWeight(c.fFuncWeight),
 234         fHistoMeasNch(c.fHistoMeasNch),
 235         fHistoMCNch(c.fHistoMCNch)
 236 {
 237         //
 238         // Copy Constructor
 239         //
 240 }
 241
 242 //___________________________________________________________________________
 243 AliCFTaskVertexingHF::~AliCFTaskVertexingHF()
 244 {
 245         //
 246         //destructor
 247         //
 248         if (fCFManager)           delete fCFManager ;
 249         if (fHistEventsProcessed) delete fHistEventsProcessed ;
 250         if (fCorrelation)         delete fCorrelation ;
 251         if (fCuts)                delete fCuts;
 252         if (fFuncWeight)          delete fFuncWeight;
 253         if (fHistoMeasNch)        delete fHistoMeasNch;
 254         if (fHistoMCNch)          delete fHistoMCNch;
 255 }
 256
 257 //_________________________________________________________________________-
 258 void AliCFTaskVertexingHF::Init()
 259 {
 260         //
 261         // Initialization
 262         //
 263
 264         if (fDebug>1) printf("AliCFTaskVertexingHF::Init()");
 265         if(fUseWeight && fUseZWeight) { AliFatal("Can not use at the same time pt and z-vtx weights, please choose"); return; }
 266         if(fUseWeight && fUseNchWeight) { AliFatal("Can not use at the same time pt and Nch weights, please choose"); return; }
 267         if(fUseNchWeight && !fHistoMCNch) { AliFatal("Need to pass the MC Nch distribution to use Nch weights"); return; }
 268         if(fUseNchWeight) CreateMeasuredNchHisto();
 269
 270         AliRDHFCuts *copyfCuts = 0x0;
 271         if (!fCuts){
 272                 AliFatal("No cuts defined - Exiting...");
 273                 return;
 274         }
 275
 276         switch (fDecayChannel){
 277         case 2:{
 278                 copyfCuts = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
 279                 switch (fConfiguration) {
 280                 case kSnail:  // slow configuration: all variables in
 281                         fNvar = 16;
 282                         break;
 283                 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
 284                         fNvar = 8;
 285                         break;
 286                 }
 287                 fPartName="D0";
 288                 fDauNames="K+pi";
 289                 break;
 290         }
 291         case 21:{
 292                 copyfCuts = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
 293                 switch (fConfiguration) {
 294                 case kSnail:  // slow configuration: all variables in
 295                         fNvar = 16;
 296                         break;
 297                 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
 298                         fNvar = 8;
 299                         break;
 300                 }
 301                 fPartName="Dstar";
 302                 fDauNames="K+pi+pi";
 303                 break;
 304         }
 305         case 31:{
 306                 copyfCuts = new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fCuts)));
 307                 switch (fConfiguration) {
 308                 case kSnail:  // slow configuration: all variables in
 309                         fNvar = 14;
 310                         break;
 311                 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
 312                         fNvar = 8;
 313                         break;
 314                 }
 315                 fPartName="Dplus";
 316                 fDauNames="K+pi+pi";
 317                 break;
 318         }
 319         case 32:{
 320                 copyfCuts = new AliRDHFCutsLctopKpi(*(static_cast<AliRDHFCutsLctopKpi*>(fCuts)));
 321                 switch (fConfiguration) {
 322                 case kSnail:  // slow configuration: all variables in
 323                         fNvar = 18;
 324                         break;
 325                 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
 326                         fNvar = 8;
 327                         break;
 328                 }
 329                 fPartName="Lambdac";
 330                 fDauNames="p+K+pi";
 331                 break;
 332         }
 333         case 33:{
 334                 copyfCuts = new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fCuts)));
 335                 switch (fConfiguration) {
 336                 case kSnail:  // slow configuration: all variables in
 337                         fNvar = 14;
 338                         break;
 339                 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
 340                         fNvar = 8;
 341                         break;
 342                 }
 343                 fPartName="Ds";
 344                 fDauNames="K+K+pi";
 345                 break;
 346         }
 347         case 4:{
 348                 copyfCuts = new AliRDHFCutsD0toKpipipi(*(static_cast<AliRDHFCutsD0toKpipipi*>(fCuts)));
 349                 switch (fConfiguration) {
 350                 case kSnail:  // slow configuration: all variables in
 351                         fNvar = 16;
 352                         break;
 353                 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
 354                         fNvar = 8;
 355                         break;
 356                 }
 357                 fPartName="D0";
 358                 fDauNames="K+pi+pi+pi";
 359                 break;
 360         }
 361         default:
 362                 AliFatal("The decay channel MUST be defined according to AliCFVertexing::DecayChannel - Exiting...");
 363                 break;
 364         }
 365
 366         const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
 367         if (copyfCuts){
 368                 copyfCuts->SetName(nameoutput);
 369
 370                 //Post the data
 371                 PostData(4, copyfCuts);
 372         }
 373         else{
 374                 AliFatal("Failing initializing AliRDHFCuts object - Exiting...");
 375         }
 376
 377         return;
 378 }
 379
 380 //_________________________________________________
 381 void AliCFTaskVertexingHF::UserExec(Option_t *)
 382 {
 383         //
 384         // Main loop function
 385         //
 386
 387         PostData(1,fHistEventsProcessed) ;
 388         PostData(2,fCFManager->GetParticleContainer()) ;
 389         PostData(3,fCorrelation) ;
 390
 391         if (fFillFromGenerated){
 392                 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
 393         }
 394
 395         if (!fInputEvent) {
 396                 Error("UserExec","NO EVENT FOUND!");
 397                 return;
 398         }
 399
 400         AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 401
 402         TClonesArray *arrayBranch=0;
 403
 404         if(!aodEvent && AODEvent() && IsStandardAOD()) {
 405                 // In case there is an AOD handler writing a standard AOD, use the AOD
 406                 // event in memory rather than the input (ESD) event.
 407                 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
 408                 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 409                 // have to taken from the AOD event hold by the AliAODExtension
 410                 AliAODHandler* aodHandler = (AliAODHandler*)
 411                         ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 412                 if(aodHandler->GetExtensions()) {
 413                         AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 414                         AliAODEvent *aodFromExt = ext->GetAOD();
 415
 416                         switch (fDecayChannel){
 417                         case 2:{
 418                                 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
 419                                 break;
 420                         }
 421                         case 21:{
 422                                 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
 423                                 break;
 424                         }
 425                         case 31:
 426                         case 32:
 427                         case 33:{
 428                                 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 429                                 break;
 430                         }
 431                         case 4:{
 432                                 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
 433                                 break;
 434                         }
 435                         default:
 436                                 break;
 437                         }
 438                 }
 439         }
 440         else {
 441                 switch (fDecayChannel){
 442                 case 2:{
 443                         arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
 444                         break;
 445                 }
 446                 case 21:{
 447                         arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
 448                         break;
 449                 }
 450                 case 31:
 451                 case 32:
 452                 case 33:{
 453                         arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
 454                         break;
 455                 }
 456                 case 4:{
 457                         arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm4Prong");
 458                         break;
 459                 }
 460                 default:
 461                         break;
 462                 }
 463         }
 464
 465         AliAODVertex *aodVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
 466         if (!aodVtx) return;
 467
 468         if (!arrayBranch) {
 469                 AliError("Could not find array of HF vertices");
 470                 return;
 471         }
 472
 473         fEvents++;
 474
 475         fCFManager->SetRecEventInfo(aodEvent);
 476         fCFManager->SetMCEventInfo(aodEvent);
 477
 478         //******** DEFINE number of variables of the container***** for now set at 13, in the future in the config macro.
 479
 480         //loop on the MC event
 481
 482         TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 483         if (!mcArray) {
 484                 AliError("Could not find Monte-Carlo in AOD");
 485                 return;
 486         }
 487         Int_t icountMC = 0;
 488         Int_t icountAcc = 0;
 489         Int_t icountReco = 0;
 490         Int_t icountVertex = 0;
 491         Int_t icountRefit = 0;
 492         Int_t icountRecoAcc = 0;
 493         Int_t icountRecoITSClusters = 0;
 494         Int_t icountRecoPPR = 0;
 495         Int_t icountRecoPID = 0;
 496         Int_t cquarks = 0;
 497
 498         AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 499         if (!mcHeader) {
 500                 AliError("Could not find MC Header in AOD");
 501                 return;
 502         }
 503
 504         Double_t* containerInput = new Double_t[fNvar];
 505         Double_t* containerInputMC = new Double_t[fNvar];
 506
 507
 508         AliCFVertexingHF* cfVtxHF=0x0;
 509         switch (fDecayChannel){
 510         case 2:{
 511           cfVtxHF = new AliCFVertexingHF2Prong(mcArray, fOriginDselection);
 512           break;
 513         }
 514         case 21:{
 515           cfVtxHF = new AliCFVertexingHFCascade(mcArray, fOriginDselection);
 516           break;
 517         }
 518         case 31:
 519         case 32:
 520         case 33:{
 521           cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel);
 522           if(fDecayChannel==33){
 523             cfVtxHF->SetGeneratedDsOption(fGenDsOption);
 524           }
 525                 break;
 526         }
 527         case 4:{
 528                 //cfVtxHF = new AliCFVertexingHF4Prong(mcArray, originDselection);  // not there yet
 529                 break;
 530         }
 531         default:
 532                 break;
 533         }
 534         if (!cfVtxHF){
 535                 AliError("No AliCFVertexingHF initialized");
 536                 delete[] containerInput;
 537                 delete[] containerInputMC;
 538                 return;
 539         }
 540
 541         Double_t zPrimVertex = aodVtx ->GetZ();
 542         Double_t zMCVertex = mcHeader->GetVtxZ();
 543         Int_t runnumber = aodEvent->GetRunNumber();
 544         fWeight=1.;
 545         if(fUseZWeight) fWeight *= GetZWeight(zMCVertex,runnumber);
 546         if(fUseNchWeight){
 547           Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(mcArray,-1.0,1.0);
 548           fWeight *= GetNchWeight(nChargedMCPhysicalPrimary);
 549           AliDebug(2,Form("Using Nch weights, Mult=%d Weight=%f\n",nChargedMCPhysicalPrimary,fWeight));
 550         }
 551
 552         if (TMath::Abs(zMCVertex) > fCuts->GetMaxVtxZ()){
 553           AliDebug(3,Form("z coordinate of MC vertex = %f, it was required to be within [-%f, +%f], skipping event", zMCVertex, fCuts->GetMaxVtxZ(), fCuts->GetMaxVtxZ()));
 554           delete[] containerInput;
 555           delete[] containerInputMC;
 556           return;
 557         }
 558
 559         AliESDtrackCuts** trackCuts = new AliESDtrackCuts*[cfVtxHF->GetNProngs()];
 560         if (fDecayChannel == 21){
 561                 // for the D*, setting the third element of the array of the track cuts to those for the soft pion
 562                 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs()-1; iProng++){
 563                         trackCuts[iProng]=fCuts->GetTrackCuts();
 564                 }
 565                 trackCuts[2] = fCuts->GetTrackCutsSoftPi();
 566         }
 567         else {
 568                 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs(); iProng++){
 569                         trackCuts[iProng]=fCuts->GetTrackCuts();
 570                 }
 571         }
 572
 573         //General settings: vertex, feed down and fill reco container with generated values.
 574         cfVtxHF->SetRecoPrimVertex(zPrimVertex);
 575         cfVtxHF->SetMCPrimaryVertex(zMCVertex);
 576         cfVtxHF->SetFillFromGenerated(fFillFromGenerated);
 577         cfVtxHF->SetNVar(fNvar);
 578         cfVtxHF->SetFakeSelection(fFakeSelection);
 579         cfVtxHF->SetRejectCandidateIfNotFromQuark(fRejectIfNoQuark);
 580         cfVtxHF->SetConfiguration(fConfiguration);
 581
 582         // switch-off the trigger class selection (doesn't work for MC)
 583         fCuts->SetTriggerClass("");
 584
 585         // MC vertex, to be used, in case, for pp
 586         if (fUseMCVertex) fCuts->SetUseMCVertex();
 587
 588         if (fCentralitySelection){ // keep only the requested centrality
 589           if(fCuts->IsEventSelectedInCentrality(aodEvent)!=0) {
 590             delete[] containerInput;
 591             delete[] containerInputMC;
 592             delete [] trackCuts;
 593             return;
 594           }
 595         }  else { // keep all centralities
 596           fCuts->SetMinCentrality(0.);
 597           fCuts->SetMaxCentrality(100.);
 598         }
 599
 600         Float_t centValue = fCuts->GetCentrality(aodEvent);
 601         cfVtxHF->SetCentralityValue(centValue);
 602
 603         // number of tracklets - multiplicity estimator
 604         Double_t multiplicity = (Double_t)(aodEvent->GetTracklets()->GetNumberOfTracklets()); // casted to double because the CF is filled with doubles
 605         cfVtxHF->SetMultiplicity(multiplicity);
 606
 607         for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
 608           AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
 609           if (!mcPart){
 610             AliError("Failed casting particle from MC array!, Skipping particle");
 611             continue;
 612           }
 613           // check the MC-level cuts, must be the desidered particle
 614           if (!fCFManager->CheckParticleCuts(0, mcPart)) {
 615             continue;  // 0 stands for MC level
 616           }
 617           cfVtxHF->SetMCCandidateParam(iPart);
 618
 619           //counting c quarks
 620           cquarks += cfVtxHF->MCcquarkCounting(mcPart);
 621
 622           if (!(cfVtxHF->SetLabelArray())){
 623             AliDebug(2,Form("Impossible to set the label array (decaychannel = %d)",fDecayChannel));
 624             continue;
 625           }
 626
 627           //check the candiate family at MC level
 628           if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) {
 629             AliDebug(2,Form("Check on the family wrong!!! (decaychannel = %d)",fDecayChannel));
 630             continue;
 631           }
 632                 else{
 633                         AliDebug(2,Form("Check on the family OK!!! (decaychannel = %d)",fDecayChannel));
 634                 }
 635
 636                 //Fill the MC container
 637                 Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC);
 638                 if (mcContainerFilled) {
 639                         if (fUseWeight){
 640                                 if (fFuncWeight){ // using user-defined function
 641                                         AliDebug(2,"Using function");
 642                                         fWeight = fFuncWeight->Eval(containerInputMC[0]);
 643                                 }
 644                                 else{ // using FONLL
 645                                         AliDebug(2,"Using FONLL");
 646                                         fWeight = GetWeight(containerInputMC[0]);
 647                                 }
 648                                 AliDebug(2,Form("pt = %f, weight = %f",containerInputMC[0], fWeight));
 649                         }
 650                         if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) continue;
 651                         //MC Limited Acceptance
 652                         if (TMath::Abs(containerInputMC[1]) < 0.5) {
 653                                 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc, fWeight);
 654                                 AliDebug(3,"MC Lim Acc container filled\n");
 655                         }
 656
 657                         //MC
 658                         fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated, fWeight);
 659                         icountMC++;
 660                         AliDebug(3,"MC cointainer filled \n");
 661
 662                         // MC in acceptance
 663                         // check the MC-Acceptance level cuts
 664                         // since standard CF functions are not applicable, using Kine Cuts on daughters
 665                         Bool_t mcAccepStep = cfVtxHF-> MCAcceptanceStep();
 666                         if (mcAccepStep){
 667                                 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance, fWeight);
 668                                 AliDebug(3,"MC acceptance cut passed\n");
 669                                 icountAcc++;
 670
 671                                 //MC Vertex step
 672                                 if (fCuts->IsEventSelected(aodEvent)){
 673                                         // filling the container if the vertex is ok
 674                                         fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex, fWeight) ;
 675                                         AliDebug(3,"Vertex cut passed and container filled\n");
 676                                         icountVertex++;
 677
 678                                         //mc Refit requirement
 679                                         Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, &trackCuts[0]);
 680                                         if (mcRefitStep){
 681                                                 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit, fWeight);
 682                                                 AliDebug(3,"MC Refit cut passed and container filled\n");
 683                                                 icountRefit++;
 684                                         }
 685                                         else{
 686                                                 AliDebug(3,"MC Refit cut not passed\n");
 687                                                 continue;
 688                                         }
 689                                 }
 690                                 else{
 691                                   AliDebug (3, "MC vertex step not passed\n");
 692                                   continue;
 693                                 }
 694                         }
 695                         else{
 696                                 AliDebug (3, "MC in acceptance step not passed\n");
 697                                 continue;
 698                         }
 699                 }
 700                 else {
 701                         AliDebug (3, "MC container not filled\n");
 702                 }
 703         }
 704
 705         if (cquarks<2) AliDebug(2,Form("Event with %d c-quarks", cquarks));
 706         AliDebug(2,Form("Found %i MC particles that are %s!!",icountMC,fPartName.Data()));
 707         AliDebug(2,Form("Found %i MC particles that are %s and satisfy Acc cuts!!",icountAcc,fPartName.Data()));
 708         AliDebug(2,Form("Found %i MC particles that are %s and satisfy Vertex cuts!!",icountVertex,fPartName.Data()));
 709         AliDebug(2,Form("Found %i MC particles that are %s and satisfy Refit cuts!!",icountRefit,fPartName.Data()));
 710
 711         // Now go to rec level
 712         fCountMC += icountMC;
 713         fCountAcc += icountAcc;
 714         fCountVertex+= icountVertex;
 715         fCountRefit+= icountRefit;
 716
 717         AliDebug(2,Form("Found %d vertices for decay channel %d",arrayBranch->GetEntriesFast(),fDecayChannel));
 718
 719         for(Int_t iCandid = 0; iCandid<arrayBranch->GetEntriesFast();iCandid++){
 720                 AliAODRecoDecayHF* charmCandidate=0x0;
 721                 switch (fDecayChannel){
 722                 case 2:{
 723                         charmCandidate = (AliAODRecoDecayHF2Prong*)arrayBranch->At(iCandid);
 724                         break;
 725                 }
 726                 case 21:{
 727                         charmCandidate = (AliAODRecoCascadeHF*)arrayBranch->At(iCandid);
 728                         break;
 729                 }
 730                 case 31:
 731                 case 32:
 732                 case 33:{
 733                         charmCandidate = (AliAODRecoDecayHF3Prong*)arrayBranch->At(iCandid);
 734                         break;
 735                 }
 736                 case 4:{
 737                         charmCandidate = (AliAODRecoDecayHF4Prong*)arrayBranch->At(iCandid);
 738                         break;
 739                 }
 740                 default:
 741                         break;
 742                 }
 743
 744                 Bool_t unsetvtx=kFALSE;
 745                 if(!charmCandidate->GetOwnPrimaryVtx()) {
 746                         charmCandidate->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables
 747                         unsetvtx=kTRUE;
 748                 }
 749
 750                 Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)charmCandidate);
 751                 if (!signAssociation){
 752                         charmCandidate = 0x0;
 753                         continue;
 754                 }
 755
 756                 Int_t isPartOrAntipart = cfVtxHF->CheckReflexion(fSign);
 757                 if (isPartOrAntipart == 0){
 758                         AliDebug(2, Form("The candidate pdg code doesn't match the requirement set in the task (fSign = %d)",fSign));
 759                         continue;
 760                 }
 761
 762
 763                 Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
 764                 if (recoContFilled){
 765
 766                         // weight according to pt
 767                         if (fUseWeight){
 768                                 if (fFuncWeight){ // using user-defined function
 769                                         AliDebug(2, "Using function");
 770                                         fWeight = fFuncWeight->Eval(containerInput[0]);
 771                                 }
 772                                 else{ // using FONLL
 773                                         AliDebug(2, "Using FONLL");
 774                                         fWeight = GetWeight(containerInput[0]);
 775                                 }
 776                                 AliDebug(2, Form("pt = %f, weight = %f",containerInput[0], fWeight));
 777                         }
 778
 779                         if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue;
 780
 781                         //Reco Step
 782                         Bool_t recoStep = cfVtxHF->RecoStep();
 783                         Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
 784
 785                         if (recoStep && recoContFilled && vtxCheck){
 786                                 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ;
 787                                 icountReco++;
 788                                 AliDebug(3,"Reco step  passed and container filled\n");
 789
 790                                 //Reco in the acceptance -- take care of UNFOLDING!!!!
 791                                 Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(&trackCuts[0]);
 792                                 if (recoAcceptanceStep) {
 793                                         fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ;
 794                                         icountRecoAcc++;
 795                                         AliDebug(3,"Reco acceptance cut passed and container filled\n");
 796
 797                                         if(fAcceptanceUnf){
 798                                                 Double_t fill[4]; //fill response matrix
 799                                                 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
 800                                                 if (bUnfolding) fCorrelation->Fill(fill);
 801                                         }
 802
 803                                         //Number of ITS cluster requirements
 804                                         Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks);
 805                                         if (recoITSnCluster){
 806                                                 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ;
 807                                                 icountRecoITSClusters++;
 808                                                 AliDebug(3,"Reco n ITS cluster cut passed and container filled\n");
 809
 810                                                 Bool_t iscutsusingpid = fCuts->GetIsUsePID();
 811                                                 Int_t recoAnalysisCuts = -1, recoPidSelection = -1;
 812                                                 fCuts->SetUsePID(kFALSE);
 813                                                 recoAnalysisCuts = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
 814
 815                                                 if (fDecayChannel==33){ // Ds case, where more possibilities are considered
 816                                                   Bool_t keepDs=ProcessDs(recoAnalysisCuts);
 817                                                   if(keepDs) recoAnalysisCuts=3;
 818                                                 }
 819
 820
 821                                                 fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object
 822                                                 if (recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart){
 823                                                         fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR, fWeight);
 824                                                         icountRecoPPR++;
 825                                                         AliDebug(3,"Reco Analysis cuts passed and container filled \n");
 826                                                         //pid selection
 827                                                         //recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kPID);
 828                                                         //if((fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==isPartOrAntipart)||(fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==3)){
 829                                                         recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
 830
 831                                                         if (fDecayChannel==33){ // Ds case, where more possibilities are considered
 832                                                           Bool_t keepDs=ProcessDs(recoPidSelection);
 833                                                           if(keepDs) recoPidSelection=3;
 834                                                         }
 835
 836                                                         if (recoPidSelection == 3 || recoPidSelection == isPartOrAntipart){
 837                                                                 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID, fWeight);
 838                                                                 icountRecoPID++;
 839                                                                 AliDebug(3,"Reco PID cuts passed and container filled \n");
 840                                                                 if(!fAcceptanceUnf){
 841                                                                         Double_t fill[4]; //fill response matrix
 842                                                                         Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
 843                                                                         if (bUnfolding) fCorrelation->Fill(fill);
 844                                                                 }
 845                                                         }
 846                                                         else {
 847                                                                 AliDebug(3, "Analysis Cuts step not passed \n");
 848                                                                 continue;
 849                                                         }
 850                                                 }
 851                                                 else {
 852                                                         AliDebug(3, "PID selection not passed \n");
 853                                                         continue;
 854                                                 }
 855                                         }
 856                                         else{
 857                                                 AliDebug(3, "Number of ITS cluster step not passed\n");
 858                                                 continue;
 859                                         }
 860                                 }
 861                                 else{
 862                                         AliDebug(3, "Reco acceptance step not passed\n");
 863                                         continue;
 864                                 }
 865                         }
 866                         else {
 867                                 AliDebug(3, "Reco step not passed\n");
 868                                 continue;
 869                         }
 870                 }
 871
 872                 if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
 873         } // end loop on candidate
 874
 875         fCountReco+= icountReco;
 876         fCountRecoAcc+= icountRecoAcc;
 877         fCountRecoITSClusters+= icountRecoITSClusters;
 878         fCountRecoPPR+= icountRecoPPR;
 879         fCountRecoPID+= icountRecoPID;
 880
 881         fHistEventsProcessed->Fill(0);
 882
 883         delete[] containerInput;
 884         delete[] containerInputMC;
 885         delete cfVtxHF;
 886         if (trackCuts){
 887         //      for (Int_t i=0; i<cfVtxHF->GetNProngs(); i++){
 888         //              delete [] trackCuts[i];
 889         //      }
 890         delete [] trackCuts;
 891         }
 892
 893
 894 }
 895
 896 //___________________________________________________________________________
 897 void AliCFTaskVertexingHF::Terminate(Option_t*)
 898 {
 899         // The Terminate() function is the last function to be called during
 900         // a query. It always runs on the client, it can be used to present
 901         // the results graphically or save the results to file.
 902
 903         AliAnalysisTaskSE::Terminate();
 904
 905         AliInfo(Form("Found %i MC particles that are %s in MC, in %d events",fCountMC,fPartName.Data(),fEvents));
 906         AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, in %d events",fCountAcc,fPartName.Data(),fEvents));
 907         AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, and satisfy Vertex requirement in %d events",fCountVertex,fPartName.Data(),fEvents));
 908         AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, and satisfy ITS+TPC refit requirementin %d events",fCountRefit,fPartName.Data(),fEvents));
 909         AliInfo(Form("Found %i reco %s that are decaying in %s, in %d events",fCountReco,fPartName.Data(),fDauNames.Data(),fEvents));
 910         AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and are in the requested acceptance, in %d events",fCountRecoAcc,fPartName.Data(),fDauNames.Data(),fEvents));
 911         AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and have at least 5 clusters in ITS, in %d events",fCountRecoITSClusters,fPartName.Data(),fDauNames.Data(),fEvents));
 912         AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and satisfy PPR cuts, in %d events",fCountRecoPPR,fPartName.Data(),fDauNames.Data(),fEvents));
 913         AliInfo(Form("Among the above, found %i reco %s that are decaying in %s and satisfy PPR+PID cuts, in %d events",fCountRecoPID,fPartName.Data(),fDauNames.Data(),fEvents));
 914
 915         // draw some example plots....
 916         AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
 917         if(!cont) {
 918                 printf("CONTAINER NOT FOUND\n");
 919                 return;
 920         }
 921         // projecting the containers to obtain histograms
 922         // first argument = variable, second argument = step
 923
 924         TH1D** h = new TH1D*[3];
 925         if (fConfiguration == kSnail){
 926                 //h = new TH1D[3][12];
 927                 for (Int_t ih = 0; ih<3; ih++){
 928                         h[ih] = new TH1D[12];
 929                 }
 930                 for(Int_t iC=1;iC<12; iC++){
 931                         // MC-level
 932                         h[0][iC] =   *(cont->ShowProjection(iC,0));
 933                         // MC-Acceptance level
 934                         h[1][iC] =   *(cont->ShowProjection(iC,1));
 935                         // Reco-level
 936                         h[2][iC] =   *(cont->ShowProjection(iC,4));
 937                 }
 938         }
 939         else{
 940                 //h = new TH1D[3][12];
 941                 for (Int_t ih = 0; ih<3; ih++){
 942                         h[ih] = new TH1D[8];
 943                 }
 944                 for(Int_t iC=0;iC<8; iC++){
 945                         // MC-level
 946                         h[0][iC] =   *(cont->ShowProjection(iC,0));
 947                         // MC-Acceptance level
 948                         h[1][iC] =   *(cont->ShowProjection(iC,1));
 949                         // Reco-level
 950                         h[2][iC] =   *(cont->ShowProjection(iC,4));
 951                 }
 952         }
 953         TString* titles;
 954         Int_t nvarToPlot = 0;
 955         if (fConfiguration == kSnail){
 956                 nvarToPlot = 12;
 957                 titles = new TString[nvarToPlot];
 958                 if(fDecayChannel==31){
 959                         titles[0]="pT_Dplus (GeV/c)";
 960                         titles[1]="rapidity";
 961                         titles[2]="phi (rad)";
 962                         titles[3]="cT (#mum)";
 963                         titles[4]="cosPointingAngle";
 964                         titles[5]="pT_1 (GeV/c)";
 965                         titles[6]="pT_2 (GeV/c)";
 966                         titles[7]="pT_3 (GeV/c)";
 967                         titles[8]="d0_1 (#mum)";
 968                         titles[9]="d0_2 (#mum)";
 969                         titles[10]="d0_3 (#mum)";
 970                         titles[11]="zVertex (cm)";
 971                 }else{
 972                         titles[0]="pT_D0 (GeV/c)";
 973                         titles[1]="rapidity";
 974                         titles[2]="cosThetaStar";
 975                         titles[3]="pT_pi (GeV/c)";
 976                         titles[4]="pT_K (Gev/c)";
 977                         titles[5]="cT (#mum)";
 978                         titles[6]="dca (#mum)";
 979                         titles[7]="d0_pi (#mum)";
 980                         titles[8]="d0_K (#mum)";
 981                         titles[9]="d0xd0 (#mum^2)";
 982                         titles[10]="cosPointingAngle";
 983                         titles[11]="phi (rad)";
 984
 985                 }
 986         }
 987         else{
 988                 nvarToPlot = 8;
 989                 titles = new TString[nvarToPlot];
 990                 titles[0]="pT_candidate (GeV/c)";
 991                 titles[1]="rapidity";
 992                 titles[2]="cT (#mum)";
 993                 titles[3]="phi";
 994                 titles[4]="z_{vtx}";
 995                 titles[5]="centrality";
 996                 titles[6]="fake";
 997                 titles[7]="multiplicity";
 998         }
 999
1000         Int_t markers[12]={20,24,21,25,27,28,
1001                            20,24,21,25,27,28};
1002         Int_t colors[3]={2,8,4};
1003         for(Int_t iC=0;iC<nvarToPlot; iC++){
1004                 for(Int_t iStep=0;iStep<3;iStep++){
1005                         h[iStep][iC].SetTitle(titles[iC].Data());
1006                         h[iStep][iC].GetXaxis()->SetTitle(titles[iC].Data());
1007                         Double_t maxh=h[iStep][iC].GetMaximum();
1008                         h[iStep][iC].GetYaxis()->SetRangeUser(0,maxh*1.2);
1009                         h[iStep][iC].SetMarkerStyle(markers[iC]);
1010                         h[iStep][iC].SetMarkerColor(colors[iStep]);
1011                 }
1012         }
1013
1014         gStyle->SetCanvasColor(0);
1015         gStyle->SetFrameFillColor(0);
1016         gStyle->SetTitleFillColor(0);
1017         gStyle->SetStatColor(0);
1018
1019         // drawing in 2 separate canvas for a matter of clearity
1020         TCanvas * c1 =new TCanvas(Form("c1New_%d",fDecayChannel),"Vars 0, 1, 2, 3",1100,1200);
1021         c1->Divide(3,4);
1022         Int_t iPad=1;
1023         for(Int_t iVar=0; iVar<4; iVar++){
1024                 c1->cd(iPad++);
1025                 h[0][iVar].DrawCopy("p");
1026                 c1->cd(iPad++);
1027                 h[1][iVar].DrawCopy("p");
1028                 c1->cd(iPad++);
1029                 h[2][iVar].DrawCopy("p");
1030         }
1031
1032         TCanvas * c2 =new TCanvas(Form("c2New_%d",fDecayChannel),"Vars 4, 5, 6, 7",1100,1200);
1033         c2->Divide(3,4);
1034         iPad=1;
1035         for(Int_t iVar=4; iVar<8; iVar++){
1036                 c2->cd(iPad++);
1037                 h[0][iVar].DrawCopy("p");
1038                 c2->cd(iPad++);
1039                 h[1][iVar].DrawCopy("p");
1040                 c2->cd(iPad++);
1041                 h[2][iVar].DrawCopy("p");
1042         }
1043
1044         if (fConfiguration == kSnail){
1045                 TCanvas * c3 =new TCanvas(Form("c3New_%d",fDecayChannel),"Vars 8, 9, 10, 11",1100,1200);
1046                 c3->Divide(3,4);
1047                 iPad=1;
1048                 for(Int_t iVar=8; iVar<12; iVar++){
1049                         c3->cd(iPad++);
1050                         h[0][iVar].DrawCopy("p");
1051                         c3->cd(iPad++);
1052                         h[1][iVar].DrawCopy("p");
1053                         c3->cd(iPad++);
1054                         h[2][iVar].DrawCopy("p");
1055                 }
1056         }
1057
1058
1059         THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
1060
1061         TH2D* corr1 =hcorr->Projection(0,2);
1062         TH2D* corr2 = hcorr->Projection(1,3);
1063
1064         TCanvas * c7 =new TCanvas(Form("c7New_%d",fDecayChannel),"",800,400);
1065         c7->Divide(2,1);
1066         c7->cd(1);
1067         corr1->DrawCopy("text");
1068         c7->cd(2);
1069         corr2->DrawCopy("text");
1070
1071         TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE");
1072
1073         corr1->Write();
1074         corr2->Write();
1075         for(Int_t iC=0;iC<nvarToPlot; iC++){
1076                 for(Int_t iStep=0;iStep<3;iStep++){
1077                         h[iStep][iC].Write(Form("Step%d_%s",iStep,titles[iC].Data()));
1078                 }
1079         }
1080         file_projection->Close();
1081         for (Int_t ih = 0; ih<3; ih++) delete [] h[ih];
1082         delete [] h;
1083
1084
1085 }
1086
1087 //___________________________________________________________________________
1088 void AliCFTaskVertexingHF::UserCreateOutputObjects()
1089 {
1090         //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1091         //TO BE SET BEFORE THE EXECUTION OF THE TASK
1092         //
1093         Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
1094
1095         //slot #1
1096         OpenFile(1);
1097         fHistEventsProcessed = new TH1I("CFHFchist0","",1,0,1) ;
1098
1099         PostData(1,fHistEventsProcessed) ;
1100         PostData(2,fCFManager->GetParticleContainer()) ;
1101         PostData(3,fCorrelation) ;
1102
1103 }
1104
1105
1106 //_________________________________________________________________________
1107 Double_t AliCFTaskVertexingHF::GetWeight(Float_t pt)
1108 {
1109         //
1110         // calculating the weight to fill the container
1111         //
1112
1113         // FNOLL central:
1114         // p0 = 1.63297e-01 --> 0.322643
1115         // p1 = 2.96275e+00
1116         // p2 = 2.30301e+00
1117         // p3 = 2.50000e+00
1118
1119         // PYTHIA
1120         // p0 = 1.85906e-01 --> 0.36609
1121         // p1 = 1.94635e+00
1122         // p2 = 1.40463e+00
1123         // p3 = 2.50000e+00
1124
1125         Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
1126         Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
1127
1128         Double_t dndpt_func1 = dNdptFit(pt,func1);
1129         if(fUseFlatPtWeight) dndpt_func1 = 1./30.;
1130         Double_t dndpt_func2 = dNdptFit(pt,func2);
1131         AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
1132         return dndpt_func1/dndpt_func2;
1133 }
1134
1135 //__________________________________________________________________________________________________
1136 Double_t AliCFTaskVertexingHF::dNdptFit(Float_t pt, Double_t* par)
1137 {
1138         //
1139         // calculating dNdpt
1140         //
1141
1142         Double_t denom =  TMath::Power((pt/par[1]), par[3] );
1143         Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
1144
1145         return dNdpt;
1146 }
1147
1148 //__________________________________________________________________________________________________
1149 Double_t AliCFTaskVertexingHF::GetZWeight(Float_t z, Int_t runnumber){
1150   //
1151   //  calculating the z-vtx weight for the given run range
1152   //
1153
1154   if(runnumber>146824 || runnumber<146803) return 1.0;
1155
1156   Double_t func1[3] = {1.0, -0.5, 6.5 };
1157   Double_t func2[3] = {1.0, -0.5, 5.5 };
1158
1159   Double_t dzFunc1 = DodzFit(z,func1);
1160   Double_t dzFunc2 = DodzFit(z,func2);
1161
1162   return dzFunc1/dzFunc2;
1163 }
1164
1165 //__________________________________________________________________________________________________
1166 Double_t AliCFTaskVertexingHF::DodzFit(Float_t z, Double_t* par) {
1167
1168   //
1169   // Gaussian z-vtx shape
1170   //
1171   //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
1172
1173   Double_t value =  par[0]/TMath::Sqrt(2.*TMath::Pi())/par[2]*TMath::Exp(-(z-par[1])*(z-par[1])/2./par[2]/par[2]);
1174
1175   return value;
1176 }
1177 //__________________________________________________________________________________________________
1178 Double_t AliCFTaskVertexingHF::GetNchWeight(Int_t nch){
1179   //
1180   //  calculates the Nch weight using the measured and generateed Nch distributions
1181   //
1182   if(nch<=0) return 0.;
1183   Double_t pMeas=fHistoMeasNch->GetBinContent(fHistoMeasNch->FindBin(nch));
1184   Double_t pMC=fHistoMCNch->GetBinContent(fHistoMCNch->FindBin(nch));
1185   return pMeas/pMC;
1186 }
1187 //__________________________________________________________________________________________________
1188 void AliCFTaskVertexingHF::CreateMeasuredNchHisto(){
1189   // creates historgam with measured multiplcity distribution in pp 7 TeV collisions (from Eur. Phys. J. C (2010) 68: 345–354)
1190   Double_t nchbins[66]={0.50,1.50,2.50,3.50,4.50,5.50,6.50,7.50,8.50,9.50,
1191                         10.50,11.50,12.50,13.50,14.50,15.50,16.50,17.50,18.50,19.50,
1192                         20.50,21.50,22.50,23.50,24.50,25.50,26.50,27.50,28.50,29.50,
1193                         30.50,31.50,32.50,33.50,34.50,35.50,36.50,37.50,38.50,39.50,
1194                         40.50,41.50,42.50,43.50,44.50,45.50,46.50,47.50,48.50,49.50,
1195                         50.50,51.50,52.50,53.50,54.50,55.50,56.50,57.50,58.50,59.50,
1196                         60.50,62.50,64.50,66.50,68.50,70.50};
1197   Double_t pch[65]={0.062011,0.072943,0.070771,0.067245,0.062834,0.057383,0.051499,0.04591,0.041109,0.036954,
1198                     0.03359,0.030729,0.028539,0.026575,0.024653,0.0229,0.021325,0.019768,0.018561,0.017187,
1199                     0.01604,0.014836,0.013726,0.012576,0.011481,0.010393,0.009502,0.008776,0.008024,0.007452,
1200                     0.006851,0.006428,0.00594,0.005515,0.005102,0.00469,0.004162,0.003811,0.003389,0.003071,
1201                     0.002708,0.002422,0.002184,0.001968,0.00186,0.00165,0.001577,0.001387,0.001254,0.001118,
1202                     0.001037,0.000942,0.000823,0.000736,0.000654,0.000579,0.000512,0.00049,0.00045,0.000355,
1203                     0.000296,0.000265,0.000193,0.00016,0.000126};
1204
1205   if(fHistoMeasNch) delete fHistoMeasNch;
1206   fHistoMeasNch=new TH1F("hMeaseNch","",65,nchbins);
1207   for(Int_t i=0; i<65; i++){
1208     fHistoMeasNch->SetBinContent(i+1,pch[i]);
1209     fHistoMeasNch->SetBinError(i+1,0.);
1210   }
1211 }
1212
1213 Bool_t AliCFTaskVertexingHF::ProcessDs(Int_t recoAnalysisCuts) const{
1214   // processes output of Ds is selected
1215   Bool_t keep=kFALSE;
1216   if(recoAnalysisCuts > 0){
1217     Int_t isKKpi=recoAnalysisCuts&1;
1218     Int_t ispiKK=recoAnalysisCuts&2;
1219     Int_t isPhiKKpi=recoAnalysisCuts&4;
1220     Int_t isPhipiKK=recoAnalysisCuts&8;
1221     Int_t isK0starKKpi=recoAnalysisCuts&16;
1222     Int_t isK0starpiKK=recoAnalysisCuts&32;
1223     if(fDsOption==1){
1224       if(isKKpi && isPhiKKpi) keep=kTRUE;
1225       if(ispiKK && isPhipiKK) keep=kTRUE;
1226     }
1227     else if(fDsOption==2){
1228       if(isKKpi && isK0starKKpi) keep=kTRUE;
1229       if(ispiKK && isK0starpiKK) keep=kTRUE;
1230     }
1231     else if(fDsOption==3)keep=kTRUE;
1232   }
1233   return keep;
1234 }