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