PWGHF/vertexingHF/AliCFTaskVertexingHF.cxx

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