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