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 <TProfile.h>
  36 #include <TH1I.h>
  37 #include <TStyle.h>
  38 #include <TFile.h>
  39 #include <TF1.h>
  40
  41 #include "AliCFTaskVertexingHF.h"
  42 #include "AliStack.h"
  43 #include "AliMCEvent.h"
  44 #include "AliCFManager.h"
  45 #include "AliCFContainer.h"
  46 #include "AliLog.h"
  47 #include "AliInputEventHandler.h"
  48 #include "AliAnalysisManager.h"
  49 #include "AliAODHandler.h"
  50 #include "AliAODEvent.h"
  51 #include "AliAODRecoDecay.h"
  52 #include "AliAODRecoDecayHF.h"
  53 #include "AliAODRecoDecayHF2Prong.h"
  54 #include "AliAODRecoDecayHF3Prong.h"
  55 #include "AliAODRecoDecayHF4Prong.h"
  56 #include "AliAODRecoCascadeHF.h"
  57 #include "AliAODMCParticle.h"
  58 #include "AliAODMCHeader.h"
  59 #include "AliESDtrack.h"
  60 #include "TChain.h"
  61 #include "THnSparse.h"
  62 #include "TH2D.h"
  63 #include "AliESDtrackCuts.h"
  64 #include "AliRDHFCuts.h"
  65 #include "AliRDHFCutsD0toKpi.h"
  66 #include "AliRDHFCutsDplustoKpipi.h"
  67 #include "AliRDHFCutsDStartoKpipi.h"
  68 #include "AliRDHFCutsDstoKKpi.h"
  69 #include "AliRDHFCutsLctopKpi.h"
  70 #include "AliRDHFCutsD0toKpipipi.h"
  71 #include "AliRDHFCutsLctoV0.h"
  72 #include "AliCFVertexingHF2Prong.h"
  73 #include "AliCFVertexingHF3Prong.h"
  74 #include "AliCFVertexingHFCascade.h"
  75 #include "AliCFVertexingHFLctoV0bachelor.h"
  76 #include "AliCFVertexingHF.h"
  77 #include "AliVertexingHFUtils.h"
  78 #include "AliAnalysisDataSlot.h"
  79 #include "AliAnalysisDataContainer.h"
  80 #include "AliPIDResponse.h"
  81
  82 //__________________________________________________________________________
  83 AliCFTaskVertexingHF::AliCFTaskVertexingHF() :
  84   AliAnalysisTaskSE(),
  85   fCFManager(0x0),
  86   fHistEventsProcessed(0x0),
  87   fCorrelation(0x0),
  88   fListProfiles(0),
  89   fCountMC(0),
  90   fCountAcc(0),
  91   fCountVertex(0),
  92   fCountRefit(0),
  93   fCountReco(0),
  94   fCountRecoAcc(0),
  95   fCountRecoITSClusters(0),
  96   fCountRecoPPR(0),
  97   fCountRecoPID(0),
  98   fEvents(0),
  99   fDecayChannel(0),
 100   fFillFromGenerated(kFALSE),
 101   fOriginDselection(0),
 102   fAcceptanceUnf(kTRUE),
 103   fCuts(0),
 104   fUseWeight(kFALSE),
 105   fWeight(1.),
 106   fUseFlatPtWeight(kFALSE),
 107   fUseZWeight(kFALSE),
 108   fUseNchWeight(kFALSE),
 109   fUseTrackletsWeight(kFALSE),
 110   fUseMultRatioAsWeight(kFALSE),
 111   fNvar(0),
 112   fPartName(""),
 113   fDauNames(""),
 114   fSign(2),
 115   fCentralitySelection(kTRUE),
 116   fFakeSelection(0),
 117   fRejectIfNoQuark(kTRUE),
 118   fUseMCVertex(kFALSE),
 119   fDsOption(1),
 120   fGenDsOption(3),
 121   fConfiguration(kCheetah), // by default, setting the fast configuration
 122   fFuncWeight(0x0),
 123   fHistoMeasNch(0x0),
 124   fHistoMCNch(0x0),
 125   fResonantDecay(0),
 126   fLctoV0bachelorOption(1),
 127   fGenLctoV0bachelorOption(0),
 128   fUseSelectionBit(kTRUE),
 129   fPDGcode(0),
 130   fMultiplicityEstimator(kNtrk10),
 131   fRefMult(9.26),
 132   fZvtxCorrectedNtrkEstimator(kFALSE),
 133   fIsPPData(kFALSE),
 134   fIsPPbData(kFALSE),
 135   fUseAdditionalCuts(kFALSE),
 136   fUseCutsForTMVA(kFALSE),
 137   fUseCascadeTaskForLctoV0bachelor(kFALSE),
 138   fCutOnMomConservation(0.00001)
 139 {
 140   //
 141   //Default ctor
 142   //
 143   for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
 144 }
 145 //___________________________________________________________________________
 146 AliCFTaskVertexingHF::AliCFTaskVertexingHF(const Char_t* name, AliRDHFCuts* cuts, TF1* func) :
 147   AliAnalysisTaskSE(name),
 148   fCFManager(0x0),
 149   fHistEventsProcessed(0x0),
 150   fCorrelation(0x0),
 151   fListProfiles(0),
 152   fCountMC(0),
 153   fCountAcc(0),
 154   fCountVertex(0),
 155   fCountRefit(0),
 156   fCountReco(0),
 157   fCountRecoAcc(0),
 158   fCountRecoITSClusters(0),
 159   fCountRecoPPR(0),
 160   fCountRecoPID(0),
 161   fEvents(0),
 162   fDecayChannel(0),
 163   fFillFromGenerated(kFALSE),
 164   fOriginDselection(0),
 165   fAcceptanceUnf(kTRUE),
 166   fCuts(cuts),
 167   fUseWeight(kFALSE),
 168   fWeight(1.),
 169   fUseFlatPtWeight(kFALSE),
 170   fUseZWeight(kFALSE),
 171   fUseNchWeight(kFALSE),
 172   fUseTrackletsWeight(kFALSE),
 173   fUseMultRatioAsWeight(kFALSE),
 174   fNvar(0),
 175   fPartName(""),
 176   fDauNames(""),
 177   fSign(2),
 178   fCentralitySelection(kTRUE),
 179   fFakeSelection(0),
 180   fRejectIfNoQuark(kTRUE),
 181   fUseMCVertex(kFALSE),
 182   fDsOption(1),
 183   fGenDsOption(3),
 184   fConfiguration(kCheetah),  // by default, setting the fast configuration
 185   fFuncWeight(func),
 186   fHistoMeasNch(0x0),
 187   fHistoMCNch(0x0),
 188   fResonantDecay(0),
 189   fLctoV0bachelorOption(1),
 190   fGenLctoV0bachelorOption(0),
 191   fUseSelectionBit(kTRUE),
 192   fPDGcode(0),
 193   fMultiplicityEstimator(kNtrk10),
 194   fRefMult(9.26),
 195   fZvtxCorrectedNtrkEstimator(kFALSE),
 196   fIsPPData(kFALSE),
 197   fIsPPbData(kFALSE),
 198   fUseAdditionalCuts(kFALSE),
 199   fUseCutsForTMVA(kFALSE),
 200   fUseCascadeTaskForLctoV0bachelor(kFALSE),
 201   fCutOnMomConservation(0.00001)
 202 {
 203   //
 204   // Constructor. Initialization of Inputs and Outputs
 205   //
 206   /*
 207     DefineInput(0) and DefineOutput(0)
 208     are taken care of by AliAnalysisTaskSE constructor
 209   */
 210   DefineOutput(1,TH1I::Class());
 211   DefineOutput(2,AliCFContainer::Class());
 212   DefineOutput(3,THnSparseD::Class());
 213   DefineOutput(4,AliRDHFCuts::Class());
 214   for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=0;
 215   DefineOutput(5,TList::Class()); // slot #5 keeps the zvtx Ntrakclets correction profiles
 216
 217   fCuts->PrintAll();
 218 }
 219
 220 //___________________________________________________________________________
 221 AliCFTaskVertexingHF& AliCFTaskVertexingHF::operator=(const AliCFTaskVertexingHF& c)
 222 {
 223   //
 224   // Assignment operator
 225   //
 226   if (this!=&c) {
 227     AliAnalysisTaskSE::operator=(c) ;
 228     fCFManager  = c.fCFManager;
 229     fHistEventsProcessed = c.fHistEventsProcessed;
 230     fCuts = c.fCuts;
 231     fFuncWeight = c.fFuncWeight;
 232     fHistoMeasNch = c.fHistoMeasNch;
 233     fHistoMCNch = c.fHistoMCNch;
 234     for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=c.fMultEstimatorAvg[i];
 235   }
 236   return *this;
 237 }
 238
 239 //___________________________________________________________________________
 240 AliCFTaskVertexingHF::AliCFTaskVertexingHF(const AliCFTaskVertexingHF& c) :
 241   AliAnalysisTaskSE(c),
 242   fCFManager(c.fCFManager),
 243   fHistEventsProcessed(c.fHistEventsProcessed),
 244   fCorrelation(c.fCorrelation),
 245   fListProfiles(c.fListProfiles),
 246   fCountMC(c.fCountMC),
 247   fCountAcc(c.fCountAcc),
 248   fCountVertex(c.fCountVertex),
 249   fCountRefit(c.fCountRefit),
 250   fCountReco(c.fCountReco),
 251   fCountRecoAcc(c.fCountRecoAcc),
 252   fCountRecoITSClusters(c.fCountRecoITSClusters),
 253   fCountRecoPPR(c.fCountRecoPPR),
 254   fCountRecoPID(c.fCountRecoPID),
 255   fEvents(c.fEvents),
 256   fDecayChannel(c.fDecayChannel),
 257   fFillFromGenerated(c.fFillFromGenerated),
 258   fOriginDselection(c.fOriginDselection),
 259   fAcceptanceUnf(c.fAcceptanceUnf),
 260   fCuts(c.fCuts),
 261   fUseWeight(c.fUseWeight),
 262   fWeight(c.fWeight),
 263   fUseFlatPtWeight(c.fUseFlatPtWeight),
 264   fUseZWeight(c.fUseZWeight),
 265   fUseNchWeight(c.fUseNchWeight),
 266   fUseTrackletsWeight(c.fUseTrackletsWeight),
 267   fUseMultRatioAsWeight(c.fUseMultRatioAsWeight),
 268   fNvar(c.fNvar),
 269   fPartName(c.fPartName),
 270   fDauNames(c.fDauNames),
 271   fSign(c.fSign),
 272   fCentralitySelection(c.fCentralitySelection),
 273   fFakeSelection(c.fFakeSelection),
 274   fRejectIfNoQuark(c.fRejectIfNoQuark),
 275   fUseMCVertex(c.fUseMCVertex),
 276   fDsOption(c.fDsOption),
 277   fGenDsOption(c.fGenDsOption),
 278   fConfiguration(c.fConfiguration),
 279   fFuncWeight(c.fFuncWeight),
 280   fHistoMeasNch(c.fHistoMeasNch),
 281   fHistoMCNch(c.fHistoMCNch),
 282   fResonantDecay(c.fResonantDecay),
 283   fLctoV0bachelorOption(c.fLctoV0bachelorOption),
 284   fGenLctoV0bachelorOption(c.fGenLctoV0bachelorOption),
 285   fUseSelectionBit(c.fUseSelectionBit),
 286   fPDGcode(c.fPDGcode),
 287   fMultiplicityEstimator(c.fMultiplicityEstimator),
 288   fRefMult(c.fRefMult),
 289   fZvtxCorrectedNtrkEstimator(c.fZvtxCorrectedNtrkEstimator),
 290   fIsPPData(c.fIsPPData),
 291   fIsPPbData(c.fIsPPbData),
 292   fUseAdditionalCuts(c.fUseAdditionalCuts),
 293   fUseCutsForTMVA(c.fUseCutsForTMVA),
 294   fUseCascadeTaskForLctoV0bachelor(c.fUseCascadeTaskForLctoV0bachelor),
 295   fCutOnMomConservation(c.fCutOnMomConservation)
 296 {
 297   //
 298   // Copy Constructor
 299   //
 300   for(Int_t i=0; i<4; i++) fMultEstimatorAvg[i]=c.fMultEstimatorAvg[i];
 301 }
 302
 303 //___________________________________________________________________________
 304 AliCFTaskVertexingHF::~AliCFTaskVertexingHF()
 305 {
 306   //
 307   //destructor
 308   //
 309   if (fCFManager)           delete fCFManager ;
 310   if (fHistEventsProcessed) delete fHistEventsProcessed ;
 311   if (fCorrelation)       delete fCorrelation ;
 312   if (fListProfiles)        delete fListProfiles;
 313   if (fCuts)                delete fCuts;
 314   if (fFuncWeight)          delete fFuncWeight;
 315   if (fHistoMeasNch)        delete fHistoMeasNch;
 316   if (fHistoMCNch)          delete fHistoMCNch;
 317   for(Int_t i=0; i<4; i++) { if(fMultEstimatorAvg[i]) delete fMultEstimatorAvg[i]; }
 318 }
 319
 320 //_________________________________________________________________________-
 321 void AliCFTaskVertexingHF::Init()
 322 {
 323   //
 324   // Initialization
 325   //
 326
 327   if (fDebug>1) printf("AliCFTaskVertexingHF::Init()");
 328   if(fUseWeight && fUseZWeight) { AliFatal("Can not use at the same time pt and z-vtx weights, please choose"); return; }
 329   if(fUseWeight && fUseNchWeight) { AliInfo("Beware, using at the same time pt and Nch weights, please check"); }
 330   if(fUseNchWeight && !fHistoMCNch) { AliFatal("Need to pass the MC Nch distribution to use Nch weights"); return; }
 331   if(fUseNchWeight && !fHistoMeasNch) CreateMeasuredNchHisto();
 332
 333   AliRDHFCuts *copyfCuts = 0x0;
 334   if (!fCuts){
 335     AliFatal("No cuts defined - Exiting...");
 336     return;
 337   }
 338
 339   switch (fDecayChannel){
 340   case 2:{
 341     fPDGcode = 421;
 342     copyfCuts = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
 343     switch (fConfiguration) {
 344     case kSnail:  // slow configuration: all variables in
 345       fNvar = 16;
 346       break;
 347     case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
 348       fNvar = 8;
 349       break;
 350     }
 351     fPartName="D0";
 352     fDauNames="K+pi";
 353     break;
 354   }
 355   case 21:{
 356     fPDGcode = 413;
 357     copyfCuts = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
 358     switch (fConfiguration) {
 359     case kSnail:  // slow configuration: all variables in
 360       fNvar = 16;
 361       break;
 362     case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
 363       fNvar = 8;
 364       break;
 365     }
 366     fPartName="Dstar";
 367     fDauNames="K+pi+pi";
 368     break;
 369   }
 370   case 22:{
 371     fPDGcode = 4122;
 372     copyfCuts = new AliRDHFCutsLctoV0(*(static_cast<AliRDHFCutsLctoV0*>(fCuts)));
 373     switch (fConfiguration) {
 374     case kSnail:  // slow configuration: all variables in
 375       fNvar = 16;
 376       break;
 377     case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
 378       fNvar = 8;
 379       break;
 380     }
 381     fPartName="Lambdac";
 382     fDauNames="V0+bachelor";
 383     break;
 384   }
 385   case 31:{
 386     fPDGcode = 411;
 387     copyfCuts = new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fCuts)));
 388     switch (fConfiguration) {
 389     case kSnail:  // slow configuration: all variables in
 390       fNvar = 14;
 391       break;
 392     case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
 393       fNvar = 8;
 394       break;
 395     }
 396     fPartName="Dplus";
 397     fDauNames="K+pi+pi";
 398     break;
 399   }
 400   case 32:{
 401     fPDGcode = 4122;
 402     copyfCuts = new AliRDHFCutsLctopKpi(*(static_cast<AliRDHFCutsLctopKpi*>(fCuts)));
 403     switch (fConfiguration) {
 404     case kSnail:  // slow configuration: all variables in
 405       fNvar = 18;
 406       break;
 407     case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
 408       fNvar = 8;
 409       break;
 410     }
 411     fPartName="Lambdac";
 412     fDauNames="p+K+pi";
 413     break;
 414   }
 415   case 33:{
 416     fPDGcode = 431;
 417     copyfCuts = new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fCuts)));
 418     switch (fConfiguration) {
 419     case kSnail:  // slow configuration: all variables in
 420       fNvar = 14;
 421       break;
 422     case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
 423       fNvar = 8;
 424       break;
 425     }
 426     fPartName="Ds";
 427     fDauNames="K+K+pi";
 428     break;
 429   }
 430   case 4:{
 431     fPDGcode = 421;
 432     copyfCuts = new AliRDHFCutsD0toKpipipi(*(static_cast<AliRDHFCutsD0toKpipipi*>(fCuts)));
 433     switch (fConfiguration) {
 434     case kSnail:  // slow configuration: all variables in
 435       fNvar = 16;
 436       break;
 437     case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
 438       fNvar = 8;
 439       break;
 440     }
 441     fPartName="D0";
 442     fDauNames="K+pi+pi+pi";
 443     break;
 444   }
 445   default:
 446     AliFatal("The decay channel MUST be defined according to AliCFVertexing::DecayChannel - Exiting...");
 447     break;
 448   }
 449
 450   const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
 451   if (copyfCuts){
 452     copyfCuts->SetName(nameoutput);
 453
 454     //Post the data
 455     PostData(4, copyfCuts);
 456   }
 457   else{
 458     AliFatal("Failing initializing AliRDHFCuts object - Exiting...");
 459   }
 460
 461   fListProfiles = new TList();
 462   fListProfiles->SetOwner();
 463   TString period[4];
 464   Int_t nProfiles=4;
 465
 466   if (fIsPPbData) { //if pPb, use only two estimator histos
 467      period[0] = "LHC13b"; period[1] = "LHC13c";
 468      nProfiles = 2;
 469   } else {        // else assume pp (four histos for LHC10)
 470      period[0] = "LHC10b"; period[1] = "LHC10c"; period[2] = "LHC10d"; period[3] = "LHC10e";
 471      nProfiles = 4;
 472   }
 473
 474   for(Int_t i=0; i<nProfiles; i++){
 475     if(fMultEstimatorAvg[i]){
 476       TProfile* hprof=new TProfile(*fMultEstimatorAvg[i]);
 477       hprof->SetName(Form("ProfileTrkVsZvtx%s\n",period[i].Data()));
 478       fListProfiles->Add(hprof);
 479     }
 480   }
 481   PostData(5,fListProfiles);
 482
 483   return;
 484 }
 485
 486 //_________________________________________________
 487 void AliCFTaskVertexingHF::UserExec(Option_t *)
 488 {
 489   //
 490   // Main loop function
 491   //
 492
 493   PostData(1,fHistEventsProcessed) ;
 494   PostData(2,fCFManager->GetParticleContainer()) ;
 495   PostData(3,fCorrelation) ;
 496
 497   AliDebug(3,Form("*** Processing event %d\n", fEvents));
 498
 499   if (fFillFromGenerated){
 500     AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
 501   }
 502
 503   if (!fInputEvent) {
 504     Error("UserExec","NO EVENT FOUND!");
 505     return;
 506   }
 507
 508   AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 509
 510   TClonesArray *arrayBranch=0;
 511
 512   if(!aodEvent && AODEvent() && IsStandardAOD()) {
 513     // In case there is an AOD handler writing a standard AOD, use the AOD
 514     // event in memory rather than the input (ESD) event.
 515     aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
 516     // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 517     // have to taken from the AOD event hold by the AliAODExtension
 518     AliAODHandler* aodHandler = (AliAODHandler*)
 519       ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 520     if(aodHandler->GetExtensions()) {
 521       AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 522       AliAODEvent *aodFromExt = ext->GetAOD();
 523
 524       switch (fDecayChannel){
 525       case 2:{
 526         arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
 527         break;
 528       }
 529       case 21:{
 530         arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
 531         break;
 532       }
 533       case 22:{
 534         arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("CascadesHF");
 535         break;
 536       }
 537       case 31:
 538       case 32:
 539       case 33:{
 540         arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 541         break;
 542       }
 543       case 4:{
 544         arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
 545         break;
 546       }
 547       default:
 548         break;
 549       }
 550     }
 551   }
 552   else {
 553     switch (fDecayChannel){
 554     case 2:{
 555       arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
 556       break;
 557     }
 558     case 21:{
 559       arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
 560       break;
 561     }
 562     case 22:{
 563       arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("CascadesHF");
 564       break;
 565     }
 566     case 31:
 567     case 32:
 568     case 33:{
 569       arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
 570       break;
 571     }
 572     case 4:{
 573       arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm4Prong");
 574       break;
 575     }
 576     default:
 577       break;
 578     }
 579   }
 580
 581   AliAODVertex *aodVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
 582   if (!aodVtx) {
 583     AliDebug(3, "The event was skipped due to missing vertex");
 584     return;
 585   }
 586
 587   if (!arrayBranch) {
 588     AliError("Could not find array of HF vertices");
 589     return;
 590   }
 591
 592   fEvents++;
 593
 594   fCFManager->SetRecEventInfo(aodEvent);
 595   fCFManager->SetMCEventInfo(aodEvent);
 596
 597   //******** DEFINE number of variables of the container***** for now set at 13, in the future in the config macro.
 598
 599   //loop on the MC event
 600
 601   TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 602   if (!mcArray) {
 603     AliError("Could not find Monte-Carlo in AOD");
 604     return;
 605   }
 606   Int_t icountMC = 0;
 607   Int_t icountAcc = 0;
 608   Int_t icountReco = 0;
 609   Int_t icountVertex = 0;
 610   Int_t icountRefit = 0;
 611   Int_t icountRecoAcc = 0;
 612   Int_t icountRecoITSClusters = 0;
 613   Int_t icountRecoPPR = 0;
 614   Int_t icountRecoPID = 0;
 615   Int_t cquarks = 0;
 616
 617   AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 618   if (!mcHeader) {
 619     AliError("Could not find MC Header in AOD");
 620     return;
 621   }
 622
 623   fHistEventsProcessed->Fill(0.5);
 624
 625   Double_t* containerInput = new Double_t[fNvar];
 626   Double_t* containerInputMC = new Double_t[fNvar];
 627
 628
 629   AliCFVertexingHF* cfVtxHF=0x0;
 630   switch (fDecayChannel){
 631   case 2:{
 632     cfVtxHF = new AliCFVertexingHF2Prong(mcArray, fOriginDselection);
 633     break;
 634   }
 635   case 21:{
 636     cfVtxHF = new AliCFVertexingHFCascade(mcArray, fOriginDselection);
 637     ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGcascade(413);
 638     ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGbachelor(211);
 639     ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaugh(421);
 640     ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaughForMC(421);
 641     ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaughPositive(211);
 642     ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaughNegative(321);
 643     ((AliCFVertexingHFCascade*)cfVtxHF)->SetPrimaryVertex(aodVtx);
 644     break;
 645   }
 646   case 22:{
 647     // Lc ->  K0S+proton
 648     if (fUseCascadeTaskForLctoV0bachelor){
 649       cfVtxHF = new AliCFVertexingHFCascade(mcArray, fOriginDselection);
 650       ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGcascade(4122);
 651       ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGbachelor(2212);
 652       ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaugh(310);
 653       ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaughForMC(311);
 654       ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaughPositive(211);
 655       ((AliCFVertexingHFCascade*)cfVtxHF)->SetPDGneutrDaughNegative(211);
 656       ((AliCFVertexingHFCascade*)cfVtxHF)->SetPrimaryVertex(aodVtx);
 657       ((AliCFVertexingHFCascade*)cfVtxHF)->SetCutOnMomConservation(fCutOnMomConservation);
 658       if (fUseAdditionalCuts) ((AliCFVertexingHFCascade*)cfVtxHF)->SetUseCutsForTMVA(fUseCutsForTMVA);
 659     }
 660     else {
 661       cfVtxHF = new AliCFVertexingHFLctoV0bachelor(mcArray, fOriginDselection,fGenLctoV0bachelorOption);
 662     }
 663     break;
 664   }
 665   case 31:
 666     //  case 32:
 667   case 33:{
 668     cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel);
 669     if(fDecayChannel==33){
 670       ((AliCFVertexingHF3Prong*)cfVtxHF)->SetGeneratedDsOption(fGenDsOption);
 671     }
 672     break;
 673   }
 674   case 32:{
 675     cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel,fResonantDecay);
 676   }
 677   case 4:{
 678     //cfVtxHF = new AliCFVertexingHF4Prong(mcArray, originDselection);  // not there yet
 679     break;
 680   }
 681   default:
 682     break;
 683   }
 684   if (!cfVtxHF){
 685     AliError("No AliCFVertexingHF initialized");
 686     delete[] containerInput;
 687     delete[] containerInputMC;
 688     return;
 689   }
 690
 691   Double_t zPrimVertex = aodVtx ->GetZ();
 692   Double_t zMCVertex = mcHeader->GetVtxZ();
 693   Int_t runnumber = aodEvent->GetRunNumber();
 694
 695   // Multiplicity definition with tracklets
 696   Double_t nTracklets = 0;
 697   Int_t nTrackletsEta10 = static_cast<Int_t>(AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent,-1.,1.));
 698   Int_t nTrackletsEta16 = static_cast<Int_t>(AliVertexingHFUtils::GetNumberOfTrackletsInEtaRange(aodEvent,-1.6,1.6));
 699   nTracklets = (Double_t)nTrackletsEta10;
 700   if(fMultiplicityEstimator==kNtrk10to16) { nTracklets = (Double_t)(nTrackletsEta16 - nTrackletsEta10); }
 701
 702   // Apply the Ntracklets z-vtx data driven correction
 703   if(fZvtxCorrectedNtrkEstimator) {
 704     TProfile* estimatorAvg = GetEstimatorHistogram(aodEvent);
 705     if(estimatorAvg) {
 706       Int_t nTrackletsEta10Corr = static_cast<Int_t>(AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,nTrackletsEta10,zPrimVertex,fRefMult));
 707       Int_t nTrackletsEta16Corr = static_cast<Int_t>(AliVertexingHFUtils::GetCorrectedNtracklets(estimatorAvg,nTrackletsEta16,zPrimVertex,fRefMult));
 708       nTracklets = (Double_t)nTrackletsEta10Corr;
 709       if(fMultiplicityEstimator==kNtrk10to16) { nTracklets = (Double_t)(nTrackletsEta16Corr - nTrackletsEta10Corr); }
 710     }
 711   }
 712
 713
 714   fWeight=1.;
 715   if(fUseZWeight) fWeight *= GetZWeight(zMCVertex,runnumber);
 716   if(fUseNchWeight){
 717     Int_t nChargedMCPhysicalPrimary=AliVertexingHFUtils::GetGeneratedPhysicalPrimariesInEtaRange(mcArray,-1.0,1.0);
 718     if(!fUseTrackletsWeight) fWeight *= GetNchWeight(nChargedMCPhysicalPrimary);
 719     else fWeight *= GetNchWeight(static_cast<Int_t>(nTracklets));
 720     AliDebug(2,Form("Using Nch weights, Mult=%d Weight=%f\n",nChargedMCPhysicalPrimary,fWeight));
 721   }
 722   Double_t eventWeight=fWeight;
 723
 724   if (TMath::Abs(zMCVertex) > fCuts->GetMaxVtxZ()){
 725     AliDebug(3,Form("z coordinate of MC vertex = %f, it was required to be within [-%f, +%f], skipping event", zMCVertex, fCuts->GetMaxVtxZ(), fCuts->GetMaxVtxZ()));
 726     delete[] containerInput;
 727     delete[] containerInputMC;
 728     delete cfVtxHF;
 729     return;
 730   }
 731
 732   if(aodEvent->GetTriggerMask()==0 &&
 733      (runnumber>=195344 && runnumber<=195677)){
 734     AliDebug(3,"Event rejected because of null trigger mask");
 735     delete[] containerInput;
 736     delete[] containerInputMC;
 737     delete cfVtxHF;
 738     return;
 739   }
 740
 741   AliESDtrackCuts** trackCuts = new AliESDtrackCuts*[cfVtxHF->GetNProngs()];
 742   if (fDecayChannel == 21){
 743     // for the D*, setting the third element of the array of the track cuts to those for the soft pion
 744     for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs()-1; iProng++){
 745       trackCuts[iProng]=fCuts->GetTrackCuts();
 746     }
 747     trackCuts[2] = fCuts->GetTrackCutsSoftPi();
 748   }
 749   else if (fDecayChannel == 22) {
 750     // for the Lc->V0+bachelor, setting the second and third elements of the array of the track cuts to those for the V0 daughters
 751     trackCuts[0]=fCuts->GetTrackCuts();
 752     trackCuts[1]=fCuts->GetTrackCutsV0daughters();
 753     trackCuts[2]=fCuts->GetTrackCutsV0daughters();
 754   }
 755   else {
 756     for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs(); iProng++){
 757       trackCuts[iProng]=fCuts->GetTrackCuts();
 758     }
 759   }
 760
 761   //General settings: vertex, feed down and fill reco container with generated values.
 762   cfVtxHF->SetRecoPrimVertex(zPrimVertex);
 763   cfVtxHF->SetMCPrimaryVertex(zMCVertex);
 764   cfVtxHF->SetFillFromGenerated(fFillFromGenerated);
 765   cfVtxHF->SetNVar(fNvar);
 766   cfVtxHF->SetFakeSelection(fFakeSelection);
 767   cfVtxHF->SetRejectCandidateIfNotFromQuark(fRejectIfNoQuark);
 768   cfVtxHF->SetConfiguration(fConfiguration);
 769
 770   // switch-off the trigger class selection (doesn't work for MC)
 771   fCuts->SetTriggerClass("");
 772
 773   // MC vertex, to be used, in case, for pp
 774   if (fUseMCVertex) fCuts->SetUseMCVertex();
 775
 776   if (fCentralitySelection){ // keep only the requested centrality
 777
 778     if(fCuts->IsEventSelectedInCentrality(aodEvent)!=0) {
 779       delete[] containerInput;
 780       delete[] containerInputMC;
 781       delete [] trackCuts;
 782       delete cfVtxHF;
 783       return;
 784     }
 785   }  else { // keep all centralities
 786     fCuts->SetMinCentrality(0.);
 787     fCuts->SetMaxCentrality(100.);
 788   }
 789
 790   Float_t centValue = 0.;
 791   if(!fIsPPData) centValue = fCuts->GetCentrality(aodEvent);
 792   cfVtxHF->SetCentralityValue(centValue);
 793
 794   // multiplicity estimator with VZERO
 795   Double_t vzeroMult=0;
 796   AliAODVZERO *vzeroAOD = (AliAODVZERO*)aodEvent->GetVZEROData();
 797   if(vzeroAOD) vzeroMult = vzeroAOD->GetMTotV0A() +  vzeroAOD->GetMTotV0C();
 798
 799   Double_t multiplicity = nTracklets; // set to the Ntracklet estimator
 800   if(fMultiplicityEstimator==kVZERO) { multiplicity = vzeroMult; }
 801
 802   cfVtxHF->SetMultiplicity(multiplicity);
 803
 804   //  printf("Multiplicity estimator %d, value %2.2f\n",fMultiplicityEstimator,multiplicity);
 805
 806   for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
 807     AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
 808     if (!mcPart){
 809       AliError("Failed casting particle from MC array!, Skipping particle");
 810       continue;
 811     }
 812
 813     //counting c quarks
 814     cquarks += cfVtxHF->MCcquarkCounting(mcPart);
 815
 816     // check the MC-level cuts, must be the desidered particle
 817     if (!fCFManager->CheckParticleCuts(0, mcPart)) {
 818       AliDebug(2,"Check the MC-level cuts - not desidered particle");
 819       continue;  // 0 stands for MC level
 820     }
 821     else {
 822       AliDebug(3, Form("\n\n---> COOL! we found a particle (particle %d)!!! with PDG code = %d \n\n", iPart, mcPart->GetPdgCode()));
 823     }
 824     cfVtxHF->SetMCCandidateParam(iPart);
 825
 826
 827     if (!(cfVtxHF->SetLabelArray())){
 828       AliDebug(2,Form("Impossible to set the label array for particle %d (decaychannel = %d)", iPart, fDecayChannel));
 829       continue;
 830     }
 831
 832     //check the candiate family at MC level
 833     if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) {
 834       AliDebug(2,Form("Check on the family wrong for particle %d!!! (decaychannel = %d)", iPart, fDecayChannel));
 835       continue;
 836     }
 837     else{
 838       AliDebug(2,Form("Check on the family OK for particle %d!!! (decaychannel = %d)", iPart, fDecayChannel));
 839     }
 840
 841     //Fill the MC container
 842     Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC);
 843     AliDebug(2, Form("particle = %d mcContainerFilled = %d",iPart, mcContainerFilled));
 844     if (mcContainerFilled) {
 845       if (fUseWeight){
 846         if (fFuncWeight){ // using user-defined function
 847           AliDebug(2,"Using function");
 848           fWeight = eventWeight*fFuncWeight->Eval(containerInputMC[0]);
 849         }
 850         else{ // using FONLL
 851           AliDebug(2,"Using FONLL");
 852           fWeight = eventWeight*GetWeight(containerInputMC[0]);
 853         }
 854         AliDebug(2,Form("pt = %f, weight = %f",containerInputMC[0], fWeight));
 855       }
 856       if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) {
 857         AliDebug(3, Form("Not in limited acceptance, containerInputMC[0] = %f, containerInputMC[1] = %f", containerInputMC[0], containerInputMC[1]));
 858         continue;
 859       }
 860       else{
 861         AliDebug(3, Form("YES!! in limited acceptance, containerInputMC[0] = %f, containerInputMC[1] = %f", containerInputMC[0],containerInputMC[1]));
 862       }
 863
 864       //MC Limited Acceptance
 865       if (TMath::Abs(containerInputMC[1]) < 0.5) {
 866         fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc, fWeight);
 867         AliDebug(3,"MC Lim Acc container filled\n");
 868       }
 869
 870       //MC
 871       fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated, fWeight);
 872       icountMC++;
 873       AliDebug(3,"MC container filled \n");
 874
 875       // MC in acceptance
 876       // check the MC-Acceptance level cuts
 877       // since standard CF functions are not applicable, using Kine Cuts on daughters
 878       Bool_t mcAccepStep = cfVtxHF-> MCAcceptanceStep();
 879       if (mcAccepStep){
 880         fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance, fWeight);
 881         AliDebug(3,"MC acceptance cut passed\n");
 882         icountAcc++;
 883
 884         //MC Vertex step
 885         if (fCuts->IsEventSelected(aodEvent)){
 886           // filling the container if the vertex is ok
 887           fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex, fWeight) ;
 888           AliDebug(3,"Vertex cut passed and container filled\n");
 889           icountVertex++;
 890
 891           //mc Refit requirement
 892           Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, &trackCuts[0]);
 893           if (mcRefitStep){
 894             fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit, fWeight);
 895             AliDebug(3,"MC Refit cut passed and container filled\n");
 896             icountRefit++;
 897           }
 898           else{
 899             AliDebug(3,"MC Refit cut not passed\n");
 900             continue;
 901           }
 902         }
 903         else{
 904           AliDebug (3, "MC vertex step not passed\n");
 905           continue;
 906         }
 907       }
 908       else{
 909         AliDebug (3, "MC in acceptance step not passed\n");
 910         continue;
 911       }
 912     }
 913     else {
 914       AliDebug (3, "MC container not filled\n");
 915     }
 916   }
 917
 918   if (cquarks<2) AliDebug(2,Form("Event with %d c-quarks", cquarks));
 919   AliDebug(2,Form("Found %i MC particles that are %s!!",icountMC,fPartName.Data()));
 920   AliDebug(2,Form("Found %i MC particles that are %s and satisfy Acc cuts!!",icountAcc,fPartName.Data()));
 921   AliDebug(2,Form("Found %i MC particles that are %s and satisfy Vertex cuts!!",icountVertex,fPartName.Data()));
 922   AliDebug(2,Form("Found %i MC particles that are %s and satisfy Refit cuts!!",icountRefit,fPartName.Data()));
 923
 924   // Now go to rec level
 925   fCountMC += icountMC;
 926   fCountAcc += icountAcc;
 927   fCountVertex+= icountVertex;
 928   fCountRefit+= icountRefit;
 929
 930   AliDebug(2,Form("Found %d vertices for decay channel %d",arrayBranch->GetEntriesFast(),fDecayChannel));
 931
 932   for(Int_t iCandid = 0; iCandid<arrayBranch->GetEntriesFast();iCandid++){
 933     AliAODRecoDecayHF* charmCandidate=0x0;
 934     switch (fDecayChannel){
 935     case 2:{
 936       charmCandidate = (AliAODRecoDecayHF2Prong*)arrayBranch->At(iCandid);
 937       break;
 938     }
 939     case 21:
 940     case 22:{
 941       charmCandidate = (AliAODRecoCascadeHF*)arrayBranch->At(iCandid);
 942       break;
 943     }
 944     case 31:
 945     case 32:
 946     case 33:{
 947       charmCandidate = (AliAODRecoDecayHF3Prong*)arrayBranch->At(iCandid);
 948       break;
 949     }
 950     case 4:{
 951       charmCandidate = (AliAODRecoDecayHF4Prong*)arrayBranch->At(iCandid);
 952       break;
 953     }
 954     default:
 955       break;
 956     }
 957
 958     Bool_t unsetvtx=kFALSE;
 959     if(!charmCandidate->GetOwnPrimaryVtx()) {
 960       charmCandidate->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables
 961       unsetvtx=kTRUE;
 962     }
 963
 964     Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)charmCandidate);
 965     if (!signAssociation){
 966       if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
 967       continue;
 968     }
 969
 970     Int_t isPartOrAntipart = cfVtxHF->CheckReflexion(fSign);
 971     if (isPartOrAntipart == 0){
 972       AliDebug(2, Form("The candidate pdg code doesn't match the requirement set in the task (fSign = %d)",fSign));
 973       if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
 974       continue;
 975     }
 976
 977     AliDebug(3,Form("iCandid=%d - signAssociation=%d, isPartOrAntipart=%d",iCandid, signAssociation, isPartOrAntipart));
 978
 979     Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
 980     AliDebug(3, Form("CF task: RecoContFilled for candidate %d is %d", iCandid, (Int_t)recoContFilled));
 981     if (recoContFilled){
 982
 983       // weight according to pt
 984       if (fUseWeight){
 985         if (fFuncWeight){ // using user-defined function
 986           AliDebug(2, "Using function");
 987           fWeight = eventWeight*fFuncWeight->Eval(containerInput[0]);
 988         }
 989         else{ // using FONLL
 990           AliDebug(2, "Using FONLL");
 991           fWeight = eventWeight*GetWeight(containerInput[0]);
 992         }
 993         AliDebug(2, Form("pt = %f, weight = %f",containerInput[0], fWeight));
 994       }
 995
 996       if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])){
 997         if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
 998         continue;
 999       }
1000       //Reco Step
1001       Bool_t recoStep = cfVtxHF->RecoStep();
1002       if (recoStep) AliDebug(2, Form("CF task: Reco step for candidate %d is %d", iCandid, (Int_t)recoStep));
1003       Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
1004
1005
1006       // Selection on the filtering bit
1007       Bool_t isBitSelected = kTRUE;
1008       if(fDecayChannel==2) {
1009         if(fUseSelectionBit && !charmCandidate->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)) isBitSelected = kFALSE;
1010       }else if(fDecayChannel==31){
1011         if(fUseSelectionBit && !charmCandidate->HasSelectionBit(AliRDHFCuts::kDplusCuts)) isBitSelected = kFALSE;
1012       }else if(fDecayChannel==33){
1013         if(fUseSelectionBit && !charmCandidate->HasSelectionBit(AliRDHFCuts::kDsCuts)) isBitSelected = kFALSE;
1014       }
1015       if(!isBitSelected){
1016         if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1017         continue;
1018       }
1019
1020
1021       if (recoStep && recoContFilled && vtxCheck){
1022         fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ;
1023         icountReco++;
1024         AliDebug(3,"Reco step  passed and container filled\n");
1025
1026         //Reco in the acceptance -- take care of UNFOLDING!!!!
1027         Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(&trackCuts[0]);
1028         if (recoAcceptanceStep) {
1029           fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ;
1030           icountRecoAcc++;
1031           AliDebug(3,"Reco acceptance cut passed and container filled\n");
1032
1033           if(fAcceptanceUnf){
1034             Double_t fill[4]; //fill response matrix
1035             Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fPDGcode,fill);
1036             if (bUnfolding) fCorrelation->Fill(fill);
1037           }
1038
1039           //Number of ITS cluster requirements
1040           Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks);
1041           if (recoITSnCluster){
1042             fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ;
1043             icountRecoITSClusters++;
1044             AliDebug(3,"Reco n ITS cluster cut passed and container filled\n");
1045
1046             Bool_t iscutsusingpid = fCuts->GetIsUsePID();
1047             Int_t recoAnalysisCuts = -1, recoPidSelection = -1;
1048             fCuts->SetUsePID(kFALSE);
1049             recoAnalysisCuts = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
1050
1051             if (fDecayChannel==33){ // Ds case, where more possibilities are considered
1052               Bool_t keepDs=ProcessDs(recoAnalysisCuts);
1053               if(keepDs) recoAnalysisCuts=3;
1054             }
1055             else if (fDecayChannel==22){ // Lc->V0+bachelor case, where more possibilities are considered
1056               Bool_t keepLctoV0bachelor=ProcessLctoV0Bachelor(recoAnalysisCuts);
1057               if (keepLctoV0bachelor) recoAnalysisCuts=3;
1058               AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
1059               AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
1060               AliPIDResponse* pidResponse = inputHandler->GetPIDResponse();
1061               if (fUseAdditionalCuts){
1062                 if (!((AliCFVertexingHFCascade*)cfVtxHF)->CheckAdditionalCuts(pidResponse)) recoAnalysisCuts = -1;
1063               }
1064             }
1065
1066             fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object
1067             Bool_t tempAn=(recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart);
1068             if (fDecayChannel == 32) tempAn=(recoAnalysisCuts >0 || recoAnalysisCuts == isPartOrAntipart);
1069
1070             if (tempAn){
1071               fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR, fWeight);
1072               icountRecoPPR++;
1073               AliDebug(3,"Reco Analysis cuts passed and container filled \n");
1074               //pid selection
1075               //recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kPID);
1076               //if((fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==isPartOrAntipart)||(fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==3)){
1077               recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
1078
1079               if (fDecayChannel==33){ // Ds case, where more possibilities are considered
1080                 Bool_t keepDs=ProcessDs(recoPidSelection);
1081                 if(keepDs) recoPidSelection=3;
1082               } else if (fDecayChannel==22){ // Lc->V0+bachelor case, where more possibilities are considered
1083                 Bool_t keepLctoV0bachelor=ProcessLctoV0Bachelor(recoPidSelection);
1084                 if (keepLctoV0bachelor) recoPidSelection=3;
1085               }
1086
1087               Bool_t tempPid=(recoPidSelection == 3 || recoPidSelection == isPartOrAntipart);
1088               if (fDecayChannel == 32) tempPid=(recoPidSelection >0 || recoPidSelection == isPartOrAntipart);
1089
1090               if (tempPid){
1091                 Double_t weigPID = 1.;
1092                 if (fDecayChannel == 2){ // D0 with Bayesian PID using weights
1093                   if(((AliRDHFCutsD0toKpi*)fCuts)->GetCombPID() && (((AliRDHFCutsD0toKpi*)fCuts)->GetBayesianStrategy() == AliRDHFCutsD0toKpi::kBayesWeight || ((AliRDHFCutsD0toKpi*)fCuts)->GetBayesianStrategy() == AliRDHFCutsD0toKpi::kBayesWeightNoFilter)){
1094                     if (isPartOrAntipart == 1){
1095                       weigPID = ((AliRDHFCutsD0toKpi*)fCuts)->GetWeightsNegative()[AliPID::kKaon] * ((AliRDHFCutsD0toKpi*)fCuts)->GetWeightsPositive()[AliPID::kPion];
1096                     }else if (isPartOrAntipart == 2){
1097                       weigPID = ((AliRDHFCutsD0toKpi*)fCuts)->GetWeightsPositive()[AliPID::kKaon] * ((AliRDHFCutsD0toKpi*)fCuts)->GetWeightsNegative()[AliPID::kPion];
1098                     }
1099                     if ((weigPID  < 0) || (weigPID > 1)) weigPID = 0.;
1100                   }
1101                 }else if (fDecayChannel == 33){ // Ds with Bayesian PID using weights
1102                   if(((AliRDHFCutsDstoKKpi*)fCuts)->GetPidOption()==AliRDHFCutsDstoKKpi::kBayesianWeights){
1103                     Int_t labDau0=((AliAODTrack*)charmCandidate->GetDaughter(0))->GetLabel();
1104                     AliAODMCParticle* firstDau=(AliAODMCParticle*)mcArray->UncheckedAt(TMath::Abs(labDau0));
1105                     if(firstDau){
1106                       Int_t pdgCode0=TMath::Abs(firstDau->GetPdgCode());
1107                       if(pdgCode0==321){
1108                         weigPID=((AliRDHFCutsDstoKKpi*)fCuts)->GetWeightForKKpi();
1109                       }else if(pdgCode0==211){
1110                         weigPID=((AliRDHFCutsDstoKKpi*)fCuts)->GetWeightForpiKK();
1111                       }
1112                       if ((weigPID  < 0) || (weigPID > 1)) weigPID = 0.;
1113                     }else{
1114                       weigPID=0.;
1115                     }
1116                   }
1117                 }
1118                 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID, fWeight*weigPID);
1119                 icountRecoPID++;
1120                 AliDebug(3,"Reco PID cuts passed and container filled \n");
1121                 if(!fAcceptanceUnf){
1122                   Double_t fill[4]; //fill response matrix
1123                   Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fPDGcode,fill);
1124                   if (bUnfolding) fCorrelation->Fill(fill);
1125                 }
1126               }
1127               else {
1128                 AliDebug(3, "Analysis Cuts step not passed \n");
1129                 if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1130                 continue;
1131               }
1132             }
1133             else {
1134               AliDebug(3, "PID selection not passed \n");
1135               if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1136               continue;
1137             }
1138           }
1139           else{
1140             AliDebug(3, "Number of ITS cluster step not passed\n");
1141             if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1142             continue;
1143           }
1144         }
1145         else{
1146           AliDebug(3, "Reco acceptance step not passed\n");
1147           if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1148           continue;
1149         }
1150       }
1151       else {
1152         AliDebug(3, "Reco step not passed\n");
1153         if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1154         continue;
1155       }
1156     }
1157
1158     if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1159   } // end loop on candidate
1160
1161   fCountReco+= icountReco;
1162   fCountRecoAcc+= icountRecoAcc;
1163   fCountRecoITSClusters+= icountRecoITSClusters;
1164   fCountRecoPPR+= icountRecoPPR;
1165   fCountRecoPID+= icountRecoPID;
1166
1167   fHistEventsProcessed->Fill(1.5);
1168
1169   delete[] containerInput;
1170   delete[] containerInputMC;
1171   delete cfVtxHF;
1172   if (trackCuts){
1173     //  for (Int_t i=0; i<cfVtxHF->GetNProngs(); i++){
1174     //          delete [] trackCuts[i];
1175     //  }
1176     delete [] trackCuts;
1177   }
1178
1179
1180 }
1181
1182 //___________________________________________________________________________
1183 void AliCFTaskVertexingHF::Terminate(Option_t*)
1184 {
1185   // The Terminate() function is the last function to be called during
1186   // a query. It always runs on the client, it can be used to present
1187   // the results graphically or save the results to file.
1188
1189   AliAnalysisTaskSE::Terminate();
1190
1191   AliInfo(Form("Found %i MC particles that are %s in MC, in %d events",fCountMC,fPartName.Data(),fEvents));
1192   AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, in %d events",fCountAcc,fPartName.Data(),fEvents));
1193   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));
1194   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));
1195   AliInfo(Form("Found %i reco %s that are decaying in %s, in %d events",fCountReco,fPartName.Data(),fDauNames.Data(),fEvents));
1196   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));
1197   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));
1198   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));
1199   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));
1200
1201   // draw some example plots....
1202   AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
1203   if(!cont) {
1204     printf("CONTAINER NOT FOUND\n");
1205     return;
1206   }
1207   // projecting the containers to obtain histograms
1208   // first argument = variable, second argument = step
1209
1210   TH1D** h = new TH1D*[3];
1211   Int_t nvarToPlot = 0;
1212   if (fConfiguration == kSnail){
1213     //h = new TH1D[3][12];
1214     for (Int_t ih = 0; ih<3; ih++){
1215       if(fDecayChannel==22){
1216         nvarToPlot = 16;
1217       } else {
1218         nvarToPlot = 12;
1219       }
1220       h[ih] = new TH1D[nvarToPlot];
1221     }
1222     for(Int_t iC=1;iC<nvarToPlot; iC++){
1223       // MC-level
1224       h[0][iC] =   *(cont->ShowProjection(iC,0));
1225       // MC-Acceptance level
1226       h[1][iC] =   *(cont->ShowProjection(iC,1));
1227       // Reco-level
1228       h[2][iC] =   *(cont->ShowProjection(iC,4));
1229     }
1230   }
1231   else {
1232     //h = new TH1D[3][12];
1233     nvarToPlot = 8;
1234     for (Int_t ih = 0; ih<3; ih++){
1235       h[ih] = new TH1D[nvarToPlot];
1236     }
1237     for(Int_t iC=0;iC<nvarToPlot; iC++){
1238       // MC-level
1239       h[0][iC] =   *(cont->ShowProjection(iC,0));
1240       // MC-Acceptance level
1241       h[1][iC] =   *(cont->ShowProjection(iC,1));
1242       // Reco-level
1243       h[2][iC] =   *(cont->ShowProjection(iC,4));
1244     }
1245   }
1246   TString* titles;
1247   //Int_t nvarToPlot = 0;
1248   if (fConfiguration == kSnail){
1249     if(fDecayChannel==31){
1250       //nvarToPlot = 12;
1251       titles = new TString[nvarToPlot];
1252       titles[0]="pT_Dplus (GeV/c)";
1253       titles[1]="rapidity";
1254       titles[2]="phi (rad)";
1255       titles[3]="cT (#mum)";
1256       titles[4]="cosPointingAngle";
1257       titles[5]="pT_1 (GeV/c)";
1258       titles[6]="pT_2 (GeV/c)";
1259       titles[7]="pT_3 (GeV/c)";
1260       titles[8]="d0_1 (#mum)";
1261       titles[9]="d0_2 (#mum)";
1262       titles[10]="d0_3 (#mum)";
1263       titles[11]="zVertex (cm)";
1264     } else if (fDecayChannel==22) {
1265       //nvarToPlot = 16;
1266       titles = new TString[nvarToPlot];
1267       titles[0]="p_{T}(#Lambda_{c}) [GeV/c]";
1268       titles[1]="y(#Lambda_{c})";
1269       titles[2]="#varphi(#Lambda_{c}) [rad]";
1270       titles[3]="onTheFlyStatusV0";
1271       titles[4]="z_{vtx} [cm]";
1272       titles[5]="centrality";
1273       titles[6]="fake";
1274       titles[7]="multiplicity";
1275       //titles[8]="pT(bachelor) [GeV/c]";
1276       titles[8]="p(bachelor) [GeV/c]";
1277       titles[9]="p_{T}(V0) [GeV/c]";
1278       titles[10]="y(V0)";
1279       titles[11]="#varphi(V0) [rad]";
1280       titles[12]="m_{inv}(#pi^{+}#pi^{+}) [GeV/c^{2}]";
1281       titles[13]="dcaV0 (nSigma)";
1282       titles[14]="cosine pointing angle (V0)";
1283       titles[15]="cosine pointing angle (#Lambda_{c})";
1284       //titles[16]="c#tauV0 (#mum)";
1285       //titles[17]="c#tau (#mum)";
1286     } else {
1287       //nvarToPlot = 12;
1288       titles = new TString[nvarToPlot];
1289       titles[0]="pT_D0 (GeV/c)";
1290       titles[1]="rapidity";
1291       titles[2]="cosThetaStar";
1292       titles[3]="pT_pi (GeV/c)";
1293       titles[4]="pT_K (Gev/c)";
1294       titles[5]="cT (#mum)";
1295       titles[6]="dca (#mum)";
1296       titles[7]="d0_pi (#mum)";
1297       titles[8]="d0_K (#mum)";
1298       titles[9]="d0xd0 (#mum^2)";
1299       titles[10]="cosPointingAngle";
1300       titles[11]="phi (rad)";
1301     }
1302   }
1303   else {
1304     //nvarToPlot = 8;
1305     titles = new TString[nvarToPlot];
1306     if (fDecayChannel==22) {
1307       titles[0]="p_{T}(#Lambda_{c}) [GeV/c]";
1308       titles[1]="y(#Lambda_{c})";
1309       titles[2]="#varphi(#Lambda_{c}) [rad]";
1310       titles[3]="onTheFlyStatusV0";
1311       titles[4]="z_{vtx} [cm]";
1312       titles[5]="centrality";
1313       titles[6]="fake";
1314       titles[7]="multiplicity";
1315     } else {
1316       titles[0]="pT_candidate (GeV/c)";
1317       titles[1]="rapidity";
1318       titles[2]="cT (#mum)";
1319       titles[3]="phi";
1320       titles[4]="z_{vtx}";
1321       titles[5]="centrality";
1322       titles[6]="fake";
1323       titles[7]="multiplicity";
1324     }
1325   }
1326
1327   Int_t markers[16]={20,24,21,25,27,28,
1328                      20,24,21,25,27,28,
1329                      20,24,21,25};
1330   Int_t colors[3]={2,8,4};
1331   for(Int_t iC=0;iC<nvarToPlot; iC++){
1332     for(Int_t iStep=0;iStep<3;iStep++){
1333       h[iStep][iC].SetTitle(titles[iC].Data());
1334       h[iStep][iC].GetXaxis()->SetTitle(titles[iC].Data());
1335       Double_t maxh=h[iStep][iC].GetMaximum();
1336       h[iStep][iC].GetYaxis()->SetRangeUser(0,maxh*1.2);
1337       h[iStep][iC].SetMarkerStyle(markers[iC]);
1338       h[iStep][iC].SetMarkerColor(colors[iStep]);
1339     }
1340   }
1341
1342   gStyle->SetCanvasColor(0);
1343   gStyle->SetFrameFillColor(0);
1344   gStyle->SetTitleFillColor(0);
1345   gStyle->SetStatColor(0);
1346
1347   // drawing in 2 separate canvas for a matter of clearity
1348   TCanvas * c1 =new TCanvas(Form("c1New_%d",fDecayChannel),"Vars 0, 1, 2, 3",1100,1200);
1349   c1->Divide(3,4);
1350   Int_t iPad=1;
1351   for(Int_t iVar=0; iVar<4; iVar++){
1352     c1->cd(iPad++);
1353     h[0][iVar].DrawCopy("p");
1354     c1->cd(iPad++);
1355     h[1][iVar].DrawCopy("p");
1356     c1->cd(iPad++);
1357     h[2][iVar].DrawCopy("p");
1358   }
1359
1360   TCanvas * c2 =new TCanvas(Form("c2New_%d",fDecayChannel),"Vars 4, 5, 6, 7",1100,1200);
1361   c2->Divide(3,4);
1362   iPad=1;
1363   for(Int_t iVar=4; iVar<8; iVar++){
1364     c2->cd(iPad++);
1365     h[0][iVar].DrawCopy("p");
1366     c2->cd(iPad++);
1367     h[1][iVar].DrawCopy("p");
1368     c2->cd(iPad++);
1369     h[2][iVar].DrawCopy("p");
1370   }
1371
1372   if (fConfiguration == kSnail){
1373     TCanvas * c3 =new TCanvas(Form("c3New_%d",fDecayChannel),"Vars 8, 9, 10, 11",1100,1200);
1374     c3->Divide(3,4);
1375     iPad=1;
1376     for(Int_t iVar=8; iVar<12; iVar++){
1377       c3->cd(iPad++);
1378       h[0][iVar].DrawCopy("p");
1379       c3->cd(iPad++);
1380       h[1][iVar].DrawCopy("p");
1381       c3->cd(iPad++);
1382       h[2][iVar].DrawCopy("p");
1383     }
1384     if (fDecayChannel==22) {
1385       TCanvas * c4 =new TCanvas(Form("c4New_%d",fDecayChannel),"Vars 12, 13, 14, 15",1100,1200);
1386       c4->Divide(3,4);
1387       iPad=1;
1388       for(Int_t iVar=12; iVar<16; iVar++){
1389         c4->cd(iPad++);
1390         h[0][iVar].DrawCopy("p");
1391         c4->cd(iPad++);
1392         h[1][iVar].DrawCopy("p");
1393         c4->cd(iPad++);
1394         h[2][iVar].DrawCopy("p");
1395       }
1396     }
1397   }
1398
1399   /*
1400     THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
1401
1402     TH2D* corr1 = hcorr->Projection(0,2);
1403     TH2D* corr2 = hcorr->Projection(1,3);
1404
1405     TCanvas * c7 =new TCanvas(Form("c7New_%d",fDecayChannel),"",800,400);
1406     c7->Divide(2,1);
1407     c7->cd(1);
1408     corr1->DrawCopy("text");
1409     c7->cd(2);
1410     corr2->DrawCopy("text");
1411   */
1412   TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE");
1413
1414   //    corr1->Write();
1415   //    corr2->Write();
1416
1417   for(Int_t iC=0;iC<nvarToPlot; iC++){
1418     for(Int_t iStep=0;iStep<3;iStep++){
1419       h[iStep][iC].Write(Form("Step%d_%s",iStep,titles[iC].Data()));
1420     }
1421   }
1422   file_projection->Close();
1423   for (Int_t ih = 0; ih<3; ih++) delete [] h[ih];
1424   delete [] h;
1425   delete [] titles;
1426
1427 }
1428
1429 //___________________________________________________________________________
1430 void AliCFTaskVertexingHF::UserCreateOutputObjects()
1431 {
1432   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1433   //TO BE SET BEFORE THE EXECUTION OF THE TASK
1434   //
1435   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
1436
1437   //slot #1
1438   OpenFile(1);
1439   const char* nameoutput=GetOutputSlot(1)->GetContainer()->GetName();
1440   fHistEventsProcessed = new TH1I(nameoutput,"",2,0,2) ;
1441   fHistEventsProcessed->GetXaxis()->SetBinLabel(1,"Events processed (all)");
1442   fHistEventsProcessed->GetXaxis()->SetBinLabel(2,"Events analyzed (after selection)");
1443
1444   PostData(1,fHistEventsProcessed) ;
1445   PostData(2,fCFManager->GetParticleContainer()) ;
1446   PostData(3,fCorrelation) ;
1447
1448 }
1449
1450
1451 //_________________________________________________________________________
1452 void AliCFTaskVertexingHF::SetPtWeightsFromFONLL276overLHC12a17a(){
1453   // ad-hoc weight function from ratio of
1454   // pt spectra from FONLL 2.76 TeV and
1455   // pt spectra from MC production LHC12a17a (PYTHIA Perugia0 with pthard bins)
1456   if(fFuncWeight) delete fFuncWeight;
1457   fFuncWeight=new TF1("funcWeight","[0]+[1]*TMath::Exp(-[2]*x)",0.,50.);
1458   fFuncWeight->SetParameter(0,4.63891e-02);
1459   fFuncWeight->SetParameter(1,1.51674e+01);
1460   fFuncWeight->SetParameter(2,4.09941e-01);
1461   fUseWeight=kTRUE;
1462 }
1463 //_________________________________________________________________________
1464 void AliCFTaskVertexingHF::SetPtWeightsFromDataPbPb276overLHC12a17a(){
1465   // ad-hoc weight function from ratio of
1466   // D0 pt spectra in PbPb 2011 0-10% centrality and
1467  // pt spectra from MC production LHC12a17a (PYTHIA Perugia0 with pthard bins)
1468   if(fFuncWeight) delete fFuncWeight;
1469   fFuncWeight=new TF1("funcWeight","[0]+[1]/TMath::Power(x,[2])",0.05,50.);
1470   fFuncWeight->SetParameter(0,1.43116e-02);
1471   fFuncWeight->SetParameter(1,4.37758e+02);
1472   fFuncWeight->SetParameter(2,3.08583);
1473   fUseWeight=kTRUE;
1474 }
1475
1476 //_________________________________________________________________________
1477 void AliCFTaskVertexingHF::SetPtWeightsFromFONLL276overLHC12a17b(){
1478   // weight function from the ratio of the LHC12a17b MC
1479   // and FONLL calculations for pp data
1480   if(fFuncWeight) delete fFuncWeight;
1481   fFuncWeight=new TF1("funcWeight","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x)",0.15,50.);
1482   fFuncWeight->SetParameters(1.92381e+01, 5.05055e+00, 1.05314e+01, 2.5, 1.88214e-03, 3.44871e+00, -9.74325e-02, 1.97671e+00, -3.21278e-01);
1483   fUseWeight=kTRUE;
1484 }
1485
1486 //_________________________________________________________________________
1487 void AliCFTaskVertexingHF::SetPtWeightsFromFONLL276andBAMPSoverLHC12a17b(){
1488   // weight function from the ratio of the LHC12a17b MC
1489   // and FONLL calculations for pp data
1490   // corrected by the BAMPS Raa calculation for 30-50% CC
1491   if(fFuncWeight) delete fFuncWeight;
1492   fFuncWeight=new TF1("funcWeight","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x)",0.15,50.);
1493   fFuncWeight->SetParameters(6.10443e+00, 1.53487e+00, 1.99474e+00, 2.5, 5.51172e-03, 5.86590e+00, -5.46963e-01, 9.41201e-02, -1.64323e-01);
1494   fUseWeight=kTRUE;
1495 }
1496
1497 //_________________________________________________________________________
1498 void AliCFTaskVertexingHF::SetPtWeightsFromFONLL5overLHC13d3(){
1499   // weight function from the ratio of the LHC13d3 MC
1500   // and FONLL calculations for pp data
1501   if(fFuncWeight) delete fFuncWeight;
1502   fFuncWeight=new TF1("funcWeight","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x)",0.15,30.);
1503   fFuncWeight->SetParameters(2.94999e+00,3.47032e+00,2.81278e+00,2.5,1.93370e-02,3.86865e+00,-1.54113e-01,8.86944e-02,2.56267e-02);
1504   fUseWeight=kTRUE;
1505 }
1506
1507 //_________________________________________________________________________
1508 void AliCFTaskVertexingHF::SetPtWeightsFromFONLL7overLHC10f6a(){
1509   // weight function from the ratio of the LHC10f6a MC
1510   // and FONLL calculations for pp data
1511   if(fFuncWeight) delete fFuncWeight;
1512   fFuncWeight=new TF1("funcWeight","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x)",0.15,40.);
1513   fFuncWeight->SetParameters(2.41522e+01,4.92146e+00,6.72495e+00,2.5,6.15361e-03,4.78995e+00,-4.29135e-01,3.99421e-01,-1.57220e-02);
1514   fUseWeight=kTRUE;
1515 }
1516
1517 //_________________________________________________________________________
1518 void AliCFTaskVertexingHF::SetPtWeightsFromFONLL7overLHC12a12(){
1519   // weight function from the ratio of the LHC12a12 MC
1520   // and FONLL calculations for pp data
1521   if(fFuncWeight) delete fFuncWeight;
1522   fFuncWeight=new TF1("funcWeight","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x+[9])",0.15,50.);
1523   fFuncWeight->SetParameters(1.31497e+01,3.30503e+00,3.45594e+00,2.5,2.28642e-02,1.42372e+00,2.32892e-04,5.21986e-02,-2.14236e-01,3.86200e+00);
1524   fUseWeight=kTRUE;
1525 }
1526
1527 //_________________________________________________________________________
1528 void AliCFTaskVertexingHF::SetPtWeightsFromFONLL7overLHC12a12bis(){
1529   // weight function from the ratio of the LHC12a12bis MC
1530   // and FONLL calculations for pp data
1531   if(fFuncWeight) delete fFuncWeight;
1532   fFuncWeight=new TF1("funcWeight","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x+[9])",0.15,50.);
1533   fFuncWeight->SetParameters(6.54519e+00,2.74007e+00,2.48325e+00,2.5,1.61113e-01,-5.32546e-01,-3.75916e-04,2.38189e-01,-2.17561e-01,2.35975e+00);
1534   fUseWeight=kTRUE;
1535 }
1536
1537 //_________________________________________________________________________
1538 void AliCFTaskVertexingHF::SetPtWeightsFromFONLL7overLHC13e2fix(){
1539   // weight function from the ratio of the LHC13e2fix MC
1540   // and FONLL calculations for pp data
1541   if(fFuncWeight) delete fFuncWeight;
1542   fFuncWeight=new TF1("funcWeight","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x+[9])",0.15,50.);
1543   fFuncWeight->SetParameters(1.85862e+01,2.48171e+00,3.39356e+00,2.5,1.70426e-02,2.28453e+00,-4.57749e-02,5.84585e-02,-3.19719e-01,4.16789e+00);
1544   fUseWeight=kTRUE;
1545 }
1546
1547 //________________________________________________________________
1548 void AliCFTaskVertexingHF::SetPtWeightsFromFONLL5overLHC10f6a(){
1549   // weight function from the ratio of the LHC10f6a MC
1550   // and FONLL calculations for pp data
1551   if(fFuncWeight) delete fFuncWeight;
1552   fFuncWeight=new TF1("funcWeight","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x)",0.15,40.);
1553   fFuncWeight->SetParameters(2.77730e+01,4.78942e+00,7.45378e+00,2.5,9.86255e-02,2.30120e+00,-4.16435e-01,3.43770e-01,-2.29380e-02);
1554   fUseWeight=kTRUE;
1555 }
1556
1557 //________________________________________________________________
1558 void AliCFTaskVertexingHF::SetPtWeightsFromFONLL276overLHC10f6a(){
1559   // weight function from the ratio of the LHC10f6a MC
1560   // and FONLL calculations for pp data
1561   if(fFuncWeight) delete fFuncWeight;
1562   fFuncWeight=new TF1("funcWeight","([0]*x)/TMath::Power([2],(1+TMath::Power([3],x/[1])))+[4]*TMath::Exp([5]+[6]*x)+[7]*TMath::Exp([8]*x+[9])",0.15,40.);
1563   fFuncWeight->SetParameters(1.34412e+01,3.20068e+00,5.14481e+00,2.5,7.59405e-04,7.51821e+00,-3.93811e-01,2.16849e-02,-3.37768e-02,2.40308e+00);
1564   fUseWeight=kTRUE;
1565 }
1566
1567 //_________________________________________________________________________
1568 Double_t AliCFTaskVertexingHF::GetWeight(Float_t pt)
1569 {
1570   //
1571   // calculating the weight to fill the container
1572   //
1573
1574   // FNOLL central:
1575   // p0 = 1.63297e-01 --> 0.322643
1576   // p1 = 2.96275e+00
1577   // p2 = 2.30301e+00
1578   // p3 = 2.50000e+00
1579
1580   // PYTHIA
1581   // p0 = 1.85906e-01 --> 0.36609
1582   // p1 = 1.94635e+00
1583   // p2 = 1.40463e+00
1584   // p3 = 2.50000e+00
1585
1586   Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
1587   Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
1588
1589   Double_t dndpt_func1 = dNdptFit(pt,func1);
1590   if(fUseFlatPtWeight) dndpt_func1 = 1./30.;
1591   Double_t dndpt_func2 = dNdptFit(pt,func2);
1592   AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
1593   return dndpt_func1/dndpt_func2;
1594 }
1595
1596 //__________________________________________________________________________________________________
1597 Double_t AliCFTaskVertexingHF::dNdptFit(Float_t pt, Double_t* par)
1598 {
1599   //
1600   // calculating dNdpt
1601   //
1602
1603   Double_t denom =  TMath::Power((pt/par[1]), par[3] );
1604   Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
1605
1606   return dNdpt;
1607 }
1608
1609 //__________________________________________________________________________________________________
1610 Double_t AliCFTaskVertexingHF::GetZWeight(Float_t z, Int_t runnumber){
1611   //
1612   //  calculating the z-vtx weight for the given run range
1613   //
1614
1615   if(runnumber>146824 || runnumber<146803) return 1.0;
1616
1617   Double_t func1[3] = {1.0, -0.5, 6.5 };
1618   Double_t func2[3] = {1.0, -0.5, 5.5 };
1619
1620   Double_t dzFunc1 = DodzFit(z,func1);
1621   Double_t dzFunc2 = DodzFit(z,func2);
1622
1623   return dzFunc1/dzFunc2;
1624 }
1625
1626 //__________________________________________________________________________________________________
1627 Double_t AliCFTaskVertexingHF::DodzFit(Float_t z, Double_t* par) {
1628
1629   //
1630   // Gaussian z-vtx shape
1631   //
1632   //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
1633
1634   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]);
1635
1636   return value;
1637 }
1638 //__________________________________________________________________________________________________
1639 Double_t AliCFTaskVertexingHF::GetNchWeight(Int_t nch){
1640   //
1641   //  calculates the Nch weight using the measured and generateed Nch distributions
1642   //
1643   if(nch<=0) return 0.;
1644   if(!fHistoMeasNch || !fHistoMCNch) { AliError("Input histos to evaluate Nch weights missing"); return 0.; }
1645   Double_t pMeas=fHistoMeasNch->GetBinContent(fHistoMeasNch->FindBin(nch));
1646   Double_t pMC=fHistoMCNch->GetBinContent(fHistoMCNch->FindBin(nch));
1647   Double_t weight = pMC>0 ? pMeas/pMC : 0.;
1648   if(fUseMultRatioAsWeight)  weight = pMC;
1649   return weight;
1650 }
1651 //__________________________________________________________________________________________________
1652 void AliCFTaskVertexingHF::CreateMeasuredNchHisto(){
1653   // creates historgam with measured multiplcity distribution in pp 7 TeV collisions (from Eur. Phys. J. C (2010) 68: 345–354)
1654   //
1655   // for Nch  > 70 the points were obtained with a double NBD distribution fit
1656   // TF1 *fit1 = new TF1("fit1","[0]*(TMath::Gamma(x+[1])/(TMath::Gamma(x+1)*TMath::Gamma([1])))*(TMath::Power(([2]/[1]),x))*(TMath::Power((1+([2]/[1])),-x-[1]))"); fit1->SetParameter(0,1.);// normalization constant
1657   // fit1->SetParameter(1,1.63); // k parameter
1658   // fit1->SetParameter(2,12.8); // mean multiplicity
1659   //
1660   Double_t nchbins[82]={0.50,1.50,2.50,3.50,4.50,5.50,6.50,7.50,8.50,9.50,
1661                         10.50,11.50,12.50,13.50,14.50,15.50,16.50,17.50,18.50,19.50,
1662                         20.50,21.50,22.50,23.50,24.50,25.50,26.50,27.50,28.50,29.50,
1663                         30.50,31.50,32.50,33.50,34.50,35.50,36.50,37.50,38.50,39.50,
1664                         40.50,41.50,42.50,43.50,44.50,45.50,46.50,47.50,48.50,49.50,
1665                         50.50,51.50,52.50,53.50,54.50,55.50,56.50,57.50,58.50,59.50,
1666                         60.50,62.50,64.50,66.50,68.50,70.50,72.50,74.50,76.50,78.50,
1667                         80.50,82.50,84.50,86.50,88.50,90.50,92.50,94.50,96.50,98.50,
1668                         100.50,102.50};
1669   Double_t pch[81]={0.062011,0.072943,0.070771,0.067245,0.062834,0.057383,0.051499,0.04591,0.041109,0.036954,
1670                     0.03359,0.030729,0.028539,0.026575,0.024653,0.0229,0.021325,0.019768,0.018561,0.017187,
1671                     0.01604,0.014836,0.013726,0.012576,0.011481,0.010393,0.009502,0.008776,0.008024,0.007452,
1672                     0.006851,0.006428,0.00594,0.005515,0.005102,0.00469,0.004162,0.003811,0.003389,0.003071,
1673                     0.002708,0.002422,0.002184,0.001968,0.00186,0.00165,0.001577,0.001387,0.001254,0.001118,
1674                     0.001037,0.000942,0.000823,0.000736,0.000654,0.000579,0.000512,0.00049,0.00045,0.000355,
1675                     0.000296,0.000265,0.000193,0.00016,0.000126,0.0000851, 0.0000676,0.0000537,0.0000426, 0.0000338,
1676                     0.0000268,0.0000213,0.0000166,0.0000133,0.0000106,0.00000837,0.00000662, 0.00000524,0.00000414, 0.00000327,
1677                     0.00000258};
1678
1679   if(fHistoMeasNch) delete fHistoMeasNch;
1680   fHistoMeasNch=new TH1F("hMeaseNch","",81,nchbins);
1681   for(Int_t i=0; i<81; i++){
1682     fHistoMeasNch->SetBinContent(i+1,pch[i]);
1683     fHistoMeasNch->SetBinError(i+1,0.);
1684   }
1685 }
1686
1687 //__________________________________________________________________________________________________
1688 Bool_t AliCFTaskVertexingHF::ProcessDs(Int_t recoAnalysisCuts) const{
1689   // processes output of Ds is selected
1690   Bool_t keep=kFALSE;
1691   if(recoAnalysisCuts > 0){
1692     Int_t isKKpi=recoAnalysisCuts&1;
1693     Int_t ispiKK=recoAnalysisCuts&2;
1694     Int_t isPhiKKpi=recoAnalysisCuts&4;
1695     Int_t isPhipiKK=recoAnalysisCuts&8;
1696     Int_t isK0starKKpi=recoAnalysisCuts&16;
1697     Int_t isK0starpiKK=recoAnalysisCuts&32;
1698     if(fDsOption==1){
1699       if(isKKpi && isPhiKKpi) keep=kTRUE;
1700       if(ispiKK && isPhipiKK) keep=kTRUE;
1701     }
1702     else if(fDsOption==2){
1703       if(isKKpi && isK0starKKpi) keep=kTRUE;
1704       if(ispiKK && isK0starpiKK) keep=kTRUE;
1705     }
1706     else if(fDsOption==3)keep=kTRUE;
1707   }
1708   return keep;
1709 }
1710 //__________________________________________________________________________________________________
1711 Bool_t AliCFTaskVertexingHF::ProcessLctoV0Bachelor(Int_t recoAnalysisCuts) const{
1712   // processes output of Lc->V0+bnachelor is selected
1713   Bool_t keep=kFALSE;
1714   if(recoAnalysisCuts > 0){
1715
1716     Int_t isK0Sp = recoAnalysisCuts&1;
1717     Int_t isLambdaBarpi = recoAnalysisCuts&2;
1718     Int_t isLambdapi = recoAnalysisCuts&4;
1719
1720     if(fLctoV0bachelorOption==1){
1721       if(isK0Sp) keep=kTRUE;
1722     }
1723     else if(fLctoV0bachelorOption==2){
1724       if(isLambdaBarpi) keep=kTRUE;
1725     }
1726     else if(fLctoV0bachelorOption==4){
1727       if(isLambdapi) keep=kTRUE;
1728     }
1729     else if(fLctoV0bachelorOption==7) {
1730       if (isK0Sp || isLambdaBarpi || isLambdapi) keep=kTRUE;
1731     }
1732   }
1733   return keep;
1734 }
1735
1736 //____________________________________________________________________________
1737 TProfile* AliCFTaskVertexingHF::GetEstimatorHistogram(const AliVEvent* event){
1738   // Get Estimator Histogram from period event->GetRunNumber();
1739   //
1740   // If you select SPD tracklets in |eta|<1 you should use type == 1
1741   //
1742
1743   Int_t runNo  = event->GetRunNumber();
1744   Int_t period = -1;   // pp:  0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
1745                        // pPb: 0-LHC13b, 1-LHC13c
1746
1747   if (fIsPPbData) {    // setting run numbers for LHC13 if pPb
1748       if (runNo>195343 && runNo<195484) period = 0;
1749       if (runNo>195528 && runNo<195678) period = 1;
1750       if (period<0 || period>1) return 0;
1751   } else {             //else assume pp 2010
1752       if(runNo>114930 && runNo<117223) period = 0;
1753       if(runNo>119158 && runNo<120830) period = 1;
1754       if(runNo>122373 && runNo<126438) period = 2;
1755       if(runNo>127711 && runNo<130841) period = 3;
1756       if(period<0 || period>3) return 0;
1757   }
1758
1759   return fMultEstimatorAvg[period];
1760 }