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       if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
 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       if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
 922       continue;
 923     }
 924
 925     AliDebug(3,Form("iCandid=%d - signAssociation=%d, isPartOrAntipart=%d",iCandid, signAssociation, isPartOrAntipart));
 926
 927     Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
 928     if (recoContFilled){
 929
 930       // weight according to pt
 931       if (fUseWeight){
 932         if (fFuncWeight){ // using user-defined function
 933           AliDebug(2, "Using function");
 934           fWeight = eventWeight*fFuncWeight->Eval(containerInput[0]);
 935         }
 936         else{ // using FONLL
 937           AliDebug(2, "Using FONLL");
 938           fWeight = eventWeight*GetWeight(containerInput[0]);
 939         }
 940         AliDebug(2, Form("pt = %f, weight = %f",containerInput[0], fWeight));
 941       }
 942
 943       if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])){
 944         if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
 945         continue;
 946       }
 947       //Reco Step
 948       Bool_t recoStep = cfVtxHF->RecoStep();
 949       Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
 950
 951
 952       // Selection on the filtering bit
 953       Bool_t isBitSelected = kTRUE;
 954       if(fDecayChannel==2) {
 955         if(fUseSelectionBit && !charmCandidate->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)) isBitSelected = kFALSE;
 956       }else if(fDecayChannel==31){
 957         if(fUseSelectionBit && !charmCandidate->HasSelectionBit(AliRDHFCuts::kDplusCuts)) isBitSelected = kFALSE;
 958       }else if(fDecayChannel==33){
 959         if(fUseSelectionBit && !charmCandidate->HasSelectionBit(AliRDHFCuts::kDsCuts)) isBitSelected = kFALSE;
 960       }
 961       if(!isBitSelected){
 962         if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
 963         continue;
 964       }
 965
 966
 967       if (recoStep && recoContFilled && vtxCheck){
 968         fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ;
 969         icountReco++;
 970         AliDebug(3,"Reco step  passed and container filled\n");
 971
 972         //Reco in the acceptance -- take care of UNFOLDING!!!!
 973         Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(&trackCuts[0]);
 974         if (recoAcceptanceStep) {
 975           fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ;
 976           icountRecoAcc++;
 977           AliDebug(3,"Reco acceptance cut passed and container filled\n");
 978
 979           if(fAcceptanceUnf){
 980             Double_t fill[4]; //fill response matrix
 981             Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fPDGcode,fill);
 982             if (bUnfolding) fCorrelation->Fill(fill);
 983           }
 984
 985           //Number of ITS cluster requirements
 986           Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks);
 987           if (recoITSnCluster){
 988             fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ;
 989             icountRecoITSClusters++;
 990             AliDebug(3,"Reco n ITS cluster cut passed and container filled\n");
 991
 992             Bool_t iscutsusingpid = fCuts->GetIsUsePID();
 993             Int_t recoAnalysisCuts = -1, recoPidSelection = -1;
 994             fCuts->SetUsePID(kFALSE);
 995             recoAnalysisCuts = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
 996
 997             if (fDecayChannel==33){ // Ds case, where more possibilities are considered
 998               Bool_t keepDs=ProcessDs(recoAnalysisCuts);
 999               if(keepDs) recoAnalysisCuts=3;
1000             }
1001             else if (fDecayChannel==22){ // Lc->V0+bachelor case, where more possibilities are considered
1002               Bool_t keepLctoV0bachelor=ProcessLctoV0Bachelor(recoAnalysisCuts);
1003               if (keepLctoV0bachelor) recoAnalysisCuts=3;
1004             }
1005
1006
1007
1008             fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object
1009             Bool_t tempAn=(recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart);
1010             if (fDecayChannel == 32) tempAn=(recoAnalysisCuts >0 || recoAnalysisCuts == isPartOrAntipart);
1011
1012             if (tempAn){
1013               fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR, fWeight);
1014               icountRecoPPR++;
1015               AliDebug(3,"Reco Analysis cuts passed and container filled \n");
1016               //pid selection
1017               //recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kPID);
1018               //if((fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==isPartOrAntipart)||(fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==3)){
1019               recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
1020
1021               if (fDecayChannel==33){ // Ds case, where more possibilities are considered
1022                 Bool_t keepDs=ProcessDs(recoPidSelection);
1023                 if(keepDs) recoPidSelection=3;
1024               } else if (fDecayChannel==22){ // Lc->V0+bachelor case, where more possibilities are considered
1025                 Bool_t keepLctoV0bachelor=ProcessLctoV0Bachelor(recoPidSelection);
1026                 if (keepLctoV0bachelor) recoPidSelection=3;
1027               }
1028
1029               Bool_t tempPid=(recoPidSelection == 3 || recoPidSelection == isPartOrAntipart);
1030               if (fDecayChannel == 32) tempPid=(recoPidSelection >0 || recoPidSelection == isPartOrAntipart);
1031
1032               if (tempPid){
1033                 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID, fWeight);
1034                 icountRecoPID++;
1035                 AliDebug(3,"Reco PID cuts passed and container filled \n");
1036                 if(!fAcceptanceUnf){
1037                   Double_t fill[4]; //fill response matrix
1038                   Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fPDGcode,fill);
1039                   if (bUnfolding) fCorrelation->Fill(fill);
1040                 }
1041               }
1042               else {
1043                 AliDebug(3, "Analysis Cuts step not passed \n");
1044                 if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1045                 continue;
1046               }
1047             }
1048             else {
1049               AliDebug(3, "PID selection not passed \n");
1050               if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1051               continue;
1052             }
1053           }
1054           else{
1055             AliDebug(3, "Number of ITS cluster step not passed\n");
1056             if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1057             continue;
1058           }
1059         }
1060         else{
1061           AliDebug(3, "Reco acceptance step not passed\n");
1062           if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1063           continue;
1064         }
1065       }
1066       else {
1067         AliDebug(3, "Reco step not passed\n");
1068         if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1069         continue;
1070       }
1071     }
1072
1073     if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
1074   } // end loop on candidate
1075
1076   fCountReco+= icountReco;
1077   fCountRecoAcc+= icountRecoAcc;
1078   fCountRecoITSClusters+= icountRecoITSClusters;
1079   fCountRecoPPR+= icountRecoPPR;
1080   fCountRecoPID+= icountRecoPID;
1081
1082   fHistEventsProcessed->Fill(1.5);
1083
1084   delete[] containerInput;
1085   delete[] containerInputMC;
1086   delete cfVtxHF;
1087   if (trackCuts){
1088     //  for (Int_t i=0; i<cfVtxHF->GetNProngs(); i++){
1089     //          delete [] trackCuts[i];
1090     //  }
1091     delete [] trackCuts;
1092   }
1093
1094
1095 }
1096
1097 //___________________________________________________________________________
1098 void AliCFTaskVertexingHF::Terminate(Option_t*)
1099 {
1100   // The Terminate() function is the last function to be called during
1101   // a query. It always runs on the client, it can be used to present
1102   // the results graphically or save the results to file.
1103
1104   AliAnalysisTaskSE::Terminate();
1105
1106   AliInfo(Form("Found %i MC particles that are %s in MC, in %d events",fCountMC,fPartName.Data(),fEvents));
1107   AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, in %d events",fCountAcc,fPartName.Data(),fEvents));
1108   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));
1109   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));
1110   AliInfo(Form("Found %i reco %s that are decaying in %s, in %d events",fCountReco,fPartName.Data(),fDauNames.Data(),fEvents));
1111   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));
1112   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));
1113   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));
1114   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));
1115
1116   // draw some example plots....
1117   AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
1118   if(!cont) {
1119     printf("CONTAINER NOT FOUND\n");
1120     return;
1121   }
1122   // projecting the containers to obtain histograms
1123   // first argument = variable, second argument = step
1124
1125   TH1D** h = new TH1D*[3];
1126   Int_t nvarToPlot = 0;
1127   if (fConfiguration == kSnail){
1128     //h = new TH1D[3][12];
1129     for (Int_t ih = 0; ih<3; ih++){
1130       if(fDecayChannel==22){
1131         nvarToPlot = 16;
1132       } else {
1133         nvarToPlot = 12;
1134       }
1135       h[ih] = new TH1D[nvarToPlot];
1136     }
1137     for(Int_t iC=1;iC<nvarToPlot; iC++){
1138       // MC-level
1139       h[0][iC] =   *(cont->ShowProjection(iC,0));
1140       // MC-Acceptance level
1141       h[1][iC] =   *(cont->ShowProjection(iC,1));
1142       // Reco-level
1143       h[2][iC] =   *(cont->ShowProjection(iC,4));
1144     }
1145   }
1146   else {
1147     //h = new TH1D[3][12];
1148     nvarToPlot = 8;
1149     for (Int_t ih = 0; ih<3; ih++){
1150       h[ih] = new TH1D[nvarToPlot];
1151     }
1152     for(Int_t iC=0;iC<nvarToPlot; iC++){
1153       // MC-level
1154       h[0][iC] =   *(cont->ShowProjection(iC,0));
1155       // MC-Acceptance level
1156       h[1][iC] =   *(cont->ShowProjection(iC,1));
1157       // Reco-level
1158       h[2][iC] =   *(cont->ShowProjection(iC,4));
1159     }
1160   }
1161   TString* titles;
1162   //Int_t nvarToPlot = 0;
1163   if (fConfiguration == kSnail){
1164     if(fDecayChannel==31){
1165       //nvarToPlot = 12;
1166       titles = new TString[nvarToPlot];
1167       titles[0]="pT_Dplus (GeV/c)";
1168       titles[1]="rapidity";
1169       titles[2]="phi (rad)";
1170       titles[3]="cT (#mum)";
1171       titles[4]="cosPointingAngle";
1172       titles[5]="pT_1 (GeV/c)";
1173       titles[6]="pT_2 (GeV/c)";
1174       titles[7]="pT_3 (GeV/c)";
1175       titles[8]="d0_1 (#mum)";
1176       titles[9]="d0_2 (#mum)";
1177       titles[10]="d0_3 (#mum)";
1178       titles[11]="zVertex (cm)";
1179     } else if (fDecayChannel==22) {
1180       //nvarToPlot = 16;
1181       titles = new TString[nvarToPlot];
1182       titles[0]="p_{T}(#Lambda_{c}) [GeV/c]";
1183       titles[1]="y(#Lambda_{c})";
1184       titles[2]="#varphi(#Lambda_{c}) [rad]";
1185       titles[3]="onTheFlyStatusV0";
1186       titles[4]="z_{vtx} [cm]";
1187       titles[5]="centrality";
1188       titles[6]="fake";
1189       titles[7]="multiplicity";
1190       //titles[8]="pT(bachelor) [GeV/c]";
1191       titles[8]="p(bachelor) [GeV/c]";
1192       titles[9]="p_{T}(V0) [GeV/c]";
1193       titles[10]="y(V0)";
1194       titles[11]="#varphi(V0) [rad]";
1195       titles[12]="m_{inv}(#pi^{+}#pi^{+}) [GeV/c^{2}]";
1196       titles[13]="dcaV0 (nSigma)";
1197       titles[14]="cosine pointing angle (V0)";
1198       titles[15]="cosine pointing angle (#Lambda_{c})";
1199       //titles[16]="c#tauV0 (#mum)";
1200       //titles[17]="c#tau (#mum)";
1201     } else {
1202       //nvarToPlot = 12;
1203       titles = new TString[nvarToPlot];
1204       titles[0]="pT_D0 (GeV/c)";
1205       titles[1]="rapidity";
1206       titles[2]="cosThetaStar";
1207       titles[3]="pT_pi (GeV/c)";
1208       titles[4]="pT_K (Gev/c)";
1209       titles[5]="cT (#mum)";
1210       titles[6]="dca (#mum)";
1211       titles[7]="d0_pi (#mum)";
1212       titles[8]="d0_K (#mum)";
1213       titles[9]="d0xd0 (#mum^2)";
1214       titles[10]="cosPointingAngle";
1215       titles[11]="phi (rad)";
1216     }
1217   }
1218   else {
1219     //nvarToPlot = 8;
1220     titles = new TString[nvarToPlot];
1221     if (fDecayChannel==22) {
1222       titles[0]="p_{T}(#Lambda_{c}) [GeV/c]";
1223       titles[1]="y(#Lambda_{c})";
1224       titles[2]="#varphi(#Lambda_{c}) [rad]";
1225       titles[3]="onTheFlyStatusV0";
1226       titles[4]="z_{vtx} [cm]";
1227       titles[5]="centrality";
1228       titles[6]="fake";
1229       titles[7]="multiplicity";
1230     } else {
1231       titles[0]="pT_candidate (GeV/c)";
1232       titles[1]="rapidity";
1233       titles[2]="cT (#mum)";
1234       titles[3]="phi";
1235       titles[4]="z_{vtx}";
1236       titles[5]="centrality";
1237       titles[6]="fake";
1238       titles[7]="multiplicity";
1239     }
1240   }
1241
1242   Int_t markers[16]={20,24,21,25,27,28,
1243                      20,24,21,25,27,28,
1244                      20,24,21,25};
1245   Int_t colors[3]={2,8,4};
1246   for(Int_t iC=0;iC<nvarToPlot; iC++){
1247     for(Int_t iStep=0;iStep<3;iStep++){
1248       h[iStep][iC].SetTitle(titles[iC].Data());
1249       h[iStep][iC].GetXaxis()->SetTitle(titles[iC].Data());
1250       Double_t maxh=h[iStep][iC].GetMaximum();
1251       h[iStep][iC].GetYaxis()->SetRangeUser(0,maxh*1.2);
1252       h[iStep][iC].SetMarkerStyle(markers[iC]);
1253       h[iStep][iC].SetMarkerColor(colors[iStep]);
1254     }
1255   }
1256
1257   gStyle->SetCanvasColor(0);
1258   gStyle->SetFrameFillColor(0);
1259   gStyle->SetTitleFillColor(0);
1260   gStyle->SetStatColor(0);
1261
1262   // drawing in 2 separate canvas for a matter of clearity
1263   TCanvas * c1 =new TCanvas(Form("c1New_%d",fDecayChannel),"Vars 0, 1, 2, 3",1100,1200);
1264   c1->Divide(3,4);
1265   Int_t iPad=1;
1266   for(Int_t iVar=0; iVar<4; iVar++){
1267     c1->cd(iPad++);
1268     h[0][iVar].DrawCopy("p");
1269     c1->cd(iPad++);
1270     h[1][iVar].DrawCopy("p");
1271     c1->cd(iPad++);
1272     h[2][iVar].DrawCopy("p");
1273   }
1274
1275   TCanvas * c2 =new TCanvas(Form("c2New_%d",fDecayChannel),"Vars 4, 5, 6, 7",1100,1200);
1276   c2->Divide(3,4);
1277   iPad=1;
1278   for(Int_t iVar=4; iVar<8; iVar++){
1279     c2->cd(iPad++);
1280     h[0][iVar].DrawCopy("p");
1281     c2->cd(iPad++);
1282     h[1][iVar].DrawCopy("p");
1283     c2->cd(iPad++);
1284     h[2][iVar].DrawCopy("p");
1285   }
1286
1287   if (fConfiguration == kSnail){
1288     TCanvas * c3 =new TCanvas(Form("c3New_%d",fDecayChannel),"Vars 8, 9, 10, 11",1100,1200);
1289     c3->Divide(3,4);
1290     iPad=1;
1291     for(Int_t iVar=8; iVar<12; iVar++){
1292       c3->cd(iPad++);
1293       h[0][iVar].DrawCopy("p");
1294       c3->cd(iPad++);
1295       h[1][iVar].DrawCopy("p");
1296       c3->cd(iPad++);
1297       h[2][iVar].DrawCopy("p");
1298     }
1299     if (fDecayChannel==22) {
1300       TCanvas * c4 =new TCanvas(Form("c4New_%d",fDecayChannel),"Vars 12, 13, 14, 15",1100,1200);
1301       c4->Divide(3,4);
1302       iPad=1;
1303       for(Int_t iVar=12; iVar<16; iVar++){
1304         c4->cd(iPad++);
1305         h[0][iVar].DrawCopy("p");
1306         c4->cd(iPad++);
1307         h[1][iVar].DrawCopy("p");
1308         c4->cd(iPad++);
1309         h[2][iVar].DrawCopy("p");
1310       }
1311     }
1312   }
1313
1314   /*
1315     THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
1316
1317     TH2D* corr1 = hcorr->Projection(0,2);
1318     TH2D* corr2 = hcorr->Projection(1,3);
1319
1320     TCanvas * c7 =new TCanvas(Form("c7New_%d",fDecayChannel),"",800,400);
1321     c7->Divide(2,1);
1322     c7->cd(1);
1323     corr1->DrawCopy("text");
1324     c7->cd(2);
1325     corr2->DrawCopy("text");
1326   */
1327   TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE");
1328
1329   //    corr1->Write();
1330   //    corr2->Write();
1331
1332   for(Int_t iC=0;iC<nvarToPlot; iC++){
1333     for(Int_t iStep=0;iStep<3;iStep++){
1334       h[iStep][iC].Write(Form("Step%d_%s",iStep,titles[iC].Data()));
1335     }
1336   }
1337   file_projection->Close();
1338   for (Int_t ih = 0; ih<3; ih++) delete [] h[ih];
1339   delete [] h;
1340   delete [] titles;
1341
1342 }
1343
1344 //___________________________________________________________________________
1345 void AliCFTaskVertexingHF::UserCreateOutputObjects()
1346 {
1347   //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1348   //TO BE SET BEFORE THE EXECUTION OF THE TASK
1349   //
1350   Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
1351
1352   //slot #1
1353   OpenFile(1);
1354   const char* nameoutput=GetOutputSlot(1)->GetContainer()->GetName();
1355   fHistEventsProcessed = new TH1I(nameoutput,"",2,0,2) ;
1356   fHistEventsProcessed->GetXaxis()->SetBinLabel(1,"Events processed (all)");
1357   fHistEventsProcessed->GetXaxis()->SetBinLabel(2,"Events analyzed (after selection)");
1358
1359   PostData(1,fHistEventsProcessed) ;
1360   PostData(2,fCFManager->GetParticleContainer()) ;
1361   PostData(3,fCorrelation) ;
1362
1363 }
1364
1365
1366 //_________________________________________________________________________
1367 void AliCFTaskVertexingHF::SetPtWeightsFromFONLL276overLHC12a17a(){
1368   // ad-hoc weight function from ratio of
1369   // pt spectra from FONLL 2.76 TeV and
1370   // pt spectra from MC production LHC12a17a (PYTHIA Perugia0 with pthard bins)
1371   if(fFuncWeight) delete fFuncWeight;
1372   fFuncWeight=new TF1("funcWeight","[0]+[1]*TMath::Exp(-[2]*x)",0.,50.);
1373   fFuncWeight->SetParameter(0,4.63891e-02);
1374   fFuncWeight->SetParameter(1,1.51674e+01);
1375   fFuncWeight->SetParameter(2,4.09941e-01);
1376   fUseWeight=kTRUE;
1377 }
1378 //_________________________________________________________________________
1379 void AliCFTaskVertexingHF::SetPtWeightsFromDataPbPb276overLHC12a17a(){
1380   // ad-hoc weight function from ratio of
1381   // D0 pt spectra in PbPb 2011 0-10% centrality and
1382  // pt spectra from MC production LHC12a17a (PYTHIA Perugia0 with pthard bins)
1383   if(fFuncWeight) delete fFuncWeight;
1384   fFuncWeight=new TF1("funcWeight","[0]+[1]/TMath::Power(x,[2])",0.05,50.);
1385   fFuncWeight->SetParameter(0,1.43116e-02);
1386   fFuncWeight->SetParameter(1,4.37758e+02);
1387   fFuncWeight->SetParameter(2,3.08583);
1388   fUseWeight=kTRUE;
1389 }
1390
1391 //_________________________________________________________________________
1392 void AliCFTaskVertexingHF::SetPtWeightsFromFONLL276overLHC12a17b(){
1393   // weight function from the ratio of the LHC12a17b MC
1394   // and FONLL calculations for pp data
1395   if(fFuncWeight) delete fFuncWeight;
1396   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.);
1397   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);
1398   fUseWeight=kTRUE;
1399 }
1400
1401 //_________________________________________________________________________
1402 void AliCFTaskVertexingHF::SetPtWeightsFromFONLL276andBAMPSoverLHC12a17b(){
1403   // weight function from the ratio of the LHC12a17b MC
1404   // and FONLL calculations for pp data
1405   // corrected by the BAMPS Raa calculation for 30-50% CC
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,50.);
1408   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);
1409   fUseWeight=kTRUE;
1410 }
1411
1412 //_________________________________________________________________________
1413 void AliCFTaskVertexingHF::SetPtWeightsFromFONLL5overLHC13d3(){
1414   // weight function from the ratio of the LHC13d3 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,30.);
1418   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);
1419   fUseWeight=kTRUE;
1420 }
1421
1422 //_________________________________________________________________________
1423 void AliCFTaskVertexingHF::SetPtWeightsFromFONLL7overLHC10f6a(){
1424   // weight function from the ratio of the LHC10f6a 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)",0.15,40.);
1428   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);
1429   fUseWeight=kTRUE;
1430 }
1431
1432 //_________________________________________________________________________
1433 void AliCFTaskVertexingHF::SetPtWeightsFromFONLL7overLHC12a12(){
1434   // weight function from the ratio of the LHC12a12 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(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);
1439   fUseWeight=kTRUE;
1440 }
1441
1442 //_________________________________________________________________________
1443 void AliCFTaskVertexingHF::SetPtWeightsFromFONLL7overLHC12a12bis(){
1444   // weight function from the ratio of the LHC12a12bis 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(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);
1449   fUseWeight=kTRUE;
1450 }
1451
1452 //_________________________________________________________________________
1453 void AliCFTaskVertexingHF::SetPtWeightsFromFONLL7overLHC13e2fix(){
1454   // weight function from the ratio of the LHC13e2fix 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+[9])",0.15,50.);
1458   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);
1459   fUseWeight=kTRUE;
1460 }
1461
1462 //________________________________________________________________
1463 void AliCFTaskVertexingHF::SetPtWeightsFromFONLL5overLHC10f6a(){
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)",0.15,40.);
1468   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);
1469   fUseWeight=kTRUE;
1470 }
1471
1472 //________________________________________________________________
1473 void AliCFTaskVertexingHF::SetPtWeightsFromFONLL276overLHC10f6a(){
1474   // weight function from the ratio of the LHC10f6a MC
1475   // and FONLL calculations for pp data
1476   if(fFuncWeight) delete fFuncWeight;
1477   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.);
1478   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);
1479   fUseWeight=kTRUE;
1480 }
1481
1482 //_________________________________________________________________________
1483 Double_t AliCFTaskVertexingHF::GetWeight(Float_t pt)
1484 {
1485   //
1486   // calculating the weight to fill the container
1487   //
1488
1489   // FNOLL central:
1490   // p0 = 1.63297e-01 --> 0.322643
1491   // p1 = 2.96275e+00
1492   // p2 = 2.30301e+00
1493   // p3 = 2.50000e+00
1494
1495   // PYTHIA
1496   // p0 = 1.85906e-01 --> 0.36609
1497   // p1 = 1.94635e+00
1498   // p2 = 1.40463e+00
1499   // p3 = 2.50000e+00
1500
1501   Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
1502   Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
1503
1504   Double_t dndpt_func1 = dNdptFit(pt,func1);
1505   if(fUseFlatPtWeight) dndpt_func1 = 1./30.;
1506   Double_t dndpt_func2 = dNdptFit(pt,func2);
1507   AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
1508   return dndpt_func1/dndpt_func2;
1509 }
1510
1511 //__________________________________________________________________________________________________
1512 Double_t AliCFTaskVertexingHF::dNdptFit(Float_t pt, Double_t* par)
1513 {
1514   //
1515   // calculating dNdpt
1516   //
1517
1518   Double_t denom =  TMath::Power((pt/par[1]), par[3] );
1519   Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
1520
1521   return dNdpt;
1522 }
1523
1524 //__________________________________________________________________________________________________
1525 Double_t AliCFTaskVertexingHF::GetZWeight(Float_t z, Int_t runnumber){
1526   //
1527   //  calculating the z-vtx weight for the given run range
1528   //
1529
1530   if(runnumber>146824 || runnumber<146803) return 1.0;
1531
1532   Double_t func1[3] = {1.0, -0.5, 6.5 };
1533   Double_t func2[3] = {1.0, -0.5, 5.5 };
1534
1535   Double_t dzFunc1 = DodzFit(z,func1);
1536   Double_t dzFunc2 = DodzFit(z,func2);
1537
1538   return dzFunc1/dzFunc2;
1539 }
1540
1541 //__________________________________________________________________________________________________
1542 Double_t AliCFTaskVertexingHF::DodzFit(Float_t z, Double_t* par) {
1543
1544   //
1545   // Gaussian z-vtx shape
1546   //
1547   //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
1548
1549   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]);
1550
1551   return value;
1552 }
1553 //__________________________________________________________________________________________________
1554 Double_t AliCFTaskVertexingHF::GetNchWeight(Int_t nch){
1555   //
1556   //  calculates the Nch weight using the measured and generateed Nch distributions
1557   //
1558   if(nch<=0) return 0.;
1559   if(!fHistoMeasNch || !fHistoMCNch) { AliError("Input histos to evaluate Nch weights missing"); return 0.; }
1560   Double_t pMeas=fHistoMeasNch->GetBinContent(fHistoMeasNch->FindBin(nch));
1561   Double_t pMC=fHistoMCNch->GetBinContent(fHistoMCNch->FindBin(nch));
1562   Double_t weight = pMC>0 ? pMeas/pMC : 0.;
1563   if(fUseTrackletsWeight)  weight = pMC;
1564   return weight;
1565 }
1566 //__________________________________________________________________________________________________
1567 void AliCFTaskVertexingHF::CreateMeasuredNchHisto(){
1568   // creates historgam with measured multiplcity distribution in pp 7 TeV collisions (from Eur. Phys. J. C (2010) 68: 345–354)
1569   //
1570   // for Nch  > 70 the points were obtained with a double NBD distribution fit
1571   // 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
1572   // fit1->SetParameter(1,1.63); // k parameter
1573   // fit1->SetParameter(2,12.8); // mean multiplicity
1574   //
1575   Double_t nchbins[82]={0.50,1.50,2.50,3.50,4.50,5.50,6.50,7.50,8.50,9.50,
1576                         10.50,11.50,12.50,13.50,14.50,15.50,16.50,17.50,18.50,19.50,
1577                         20.50,21.50,22.50,23.50,24.50,25.50,26.50,27.50,28.50,29.50,
1578                         30.50,31.50,32.50,33.50,34.50,35.50,36.50,37.50,38.50,39.50,
1579                         40.50,41.50,42.50,43.50,44.50,45.50,46.50,47.50,48.50,49.50,
1580                         50.50,51.50,52.50,53.50,54.50,55.50,56.50,57.50,58.50,59.50,
1581                         60.50,62.50,64.50,66.50,68.50,70.50,72.50,74.50,76.50,78.50,
1582                         80.50,82.50,84.50,86.50,88.50,90.50,92.50,94.50,96.50,98.50,
1583                         100.50,102.50};
1584   Double_t pch[81]={0.062011,0.072943,0.070771,0.067245,0.062834,0.057383,0.051499,0.04591,0.041109,0.036954,
1585                     0.03359,0.030729,0.028539,0.026575,0.024653,0.0229,0.021325,0.019768,0.018561,0.017187,
1586                     0.01604,0.014836,0.013726,0.012576,0.011481,0.010393,0.009502,0.008776,0.008024,0.007452,
1587                     0.006851,0.006428,0.00594,0.005515,0.005102,0.00469,0.004162,0.003811,0.003389,0.003071,
1588                     0.002708,0.002422,0.002184,0.001968,0.00186,0.00165,0.001577,0.001387,0.001254,0.001118,
1589                     0.001037,0.000942,0.000823,0.000736,0.000654,0.000579,0.000512,0.00049,0.00045,0.000355,
1590                     0.000296,0.000265,0.000193,0.00016,0.000126,0.0000851, 0.0000676,0.0000537,0.0000426, 0.0000338,
1591                     0.0000268,0.0000213,0.0000166,0.0000133,0.0000106,0.00000837,0.00000662, 0.00000524,0.00000414, 0.00000327,
1592                     0.00000258};
1593
1594   if(fHistoMeasNch) delete fHistoMeasNch;
1595   fHistoMeasNch=new TH1F("hMeaseNch","",81,nchbins);
1596   for(Int_t i=0; i<81; i++){
1597     fHistoMeasNch->SetBinContent(i+1,pch[i]);
1598     fHistoMeasNch->SetBinError(i+1,0.);
1599   }
1600 }
1601
1602 //__________________________________________________________________________________________________
1603 Bool_t AliCFTaskVertexingHF::ProcessDs(Int_t recoAnalysisCuts) const{
1604   // processes output of Ds is selected
1605   Bool_t keep=kFALSE;
1606   if(recoAnalysisCuts > 0){
1607     Int_t isKKpi=recoAnalysisCuts&1;
1608     Int_t ispiKK=recoAnalysisCuts&2;
1609     Int_t isPhiKKpi=recoAnalysisCuts&4;
1610     Int_t isPhipiKK=recoAnalysisCuts&8;
1611     Int_t isK0starKKpi=recoAnalysisCuts&16;
1612     Int_t isK0starpiKK=recoAnalysisCuts&32;
1613     if(fDsOption==1){
1614       if(isKKpi && isPhiKKpi) keep=kTRUE;
1615       if(ispiKK && isPhipiKK) keep=kTRUE;
1616     }
1617     else if(fDsOption==2){
1618       if(isKKpi && isK0starKKpi) keep=kTRUE;
1619       if(ispiKK && isK0starpiKK) keep=kTRUE;
1620     }
1621     else if(fDsOption==3)keep=kTRUE;
1622   }
1623   return keep;
1624 }
1625 //__________________________________________________________________________________________________
1626 Bool_t AliCFTaskVertexingHF::ProcessLctoV0Bachelor(Int_t recoAnalysisCuts) const{
1627   // processes output of Lc->V0+bnachelor is selected
1628   Bool_t keep=kFALSE;
1629   if(recoAnalysisCuts > 0){
1630
1631     Int_t isK0Sp = recoAnalysisCuts&1;
1632     Int_t isLambdaBarpi = recoAnalysisCuts&2;
1633     Int_t isLambdapi = recoAnalysisCuts&4;
1634
1635     if(fLctoV0bachelorOption==1){
1636       if(isK0Sp) keep=kTRUE;
1637     }
1638     else if(fLctoV0bachelorOption==2){
1639       if(isLambdaBarpi) keep=kTRUE;
1640     }
1641     else if(fLctoV0bachelorOption==4){
1642       if(isLambdapi) keep=kTRUE;
1643     }
1644     else if(fLctoV0bachelorOption==7) {
1645       if (isK0Sp || isLambdaBarpi || isLambdapi) keep=kTRUE;
1646     }
1647   }
1648   return keep;
1649 }
1650
1651 //____________________________________________________________________________
1652 TProfile* AliCFTaskVertexingHF::GetEstimatorHistogram(const AliVEvent* event){
1653   // Get Estimator Histogram from period event->GetRunNumber();
1654   //
1655   // If you select SPD tracklets in |eta|<1 you should use type == 1
1656   //
1657
1658   Int_t runNo  = event->GetRunNumber();
1659   Int_t period = -1;   // pp:  0-LHC10b, 1-LHC10c, 2-LHC10d, 3-LHC10e
1660                        // pPb: 0-LHC13b, 1-LHC13c
1661
1662   if (fIsPPbData) {    // setting run numbers for LHC13 if pPb
1663       if (runNo>195343 && runNo<195484) period = 0;
1664       if (runNo>195528 && runNo<195678) period = 1;
1665       if (period<0 || period>1) return 0;
1666   } else {             //else assume pp 2010
1667       if(runNo>114930 && runNo<117223) period = 0;
1668       if(runNo>119158 && runNo<120830) period = 1;
1669       if(runNo>122373 && runNo<126438) period = 2;
1670       if(runNo>127711 && runNo<130841) period = 3;
1671       if(period<0 || period>3) return 0;
1672   }
1673
1674   return fMultEstimatorAvg[period];
1675 }