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