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 "AliAnalysisManager.h"
  47 #include "AliAODHandler.h"
  48 #include "AliAODEvent.h"
  49 #include "AliAODRecoDecay.h"
  50 #include "AliAODRecoDecayHF.h"
  51 #include "AliAODRecoDecayHF2Prong.h"
  52 #include "AliAODRecoDecayHF3Prong.h"
  53 #include "AliAODRecoDecayHF4Prong.h"
  54 #include "AliAODRecoCascadeHF.h"
  55 #include "AliAODMCParticle.h"
  56 #include "AliAODMCHeader.h"
  57 #include "AliESDtrack.h"
  58 #include "TChain.h"
  59 #include "THnSparse.h"
  60 #include "TH2D.h"
  61 #include "AliESDtrackCuts.h"
  62 #include "AliRDHFCuts.h"
  63 #include "AliRDHFCutsD0toKpi.h"
  64 #include "AliRDHFCutsDplustoKpipi.h"
  65 #include "AliRDHFCutsDStartoKpipi.h"
  66 #include "AliRDHFCutsDstoKKpi.h"
  67 #include "AliRDHFCutsLctopKpi.h"
  68 #include "AliRDHFCutsD0toKpipipi.h"
  69 #include "AliCFVertexingHF2Prong.h"
  70 #include "AliCFVertexingHF3Prong.h"
  71 #include "AliCFVertexingHFCascade.h"
  72 #include "AliCFVertexingHF.h"
  73 #include "AliAnalysisDataSlot.h"
  74 #include "AliAnalysisDataContainer.h"
  75
  76 //__________________________________________________________________________
  77 AliCFTaskVertexingHF::AliCFTaskVertexingHF() :
  78         AliAnalysisTaskSE(),
  79         fCFManager(0x0),
  80         fHistEventsProcessed(0x0),
  81         fCorrelation(0x0),
  82         fCountMC(0),
  83         fCountAcc(0),
  84         fCountVertex(0),
  85         fCountRefit(0),
  86         fCountReco(0),
  87         fCountRecoAcc(0),
  88         fCountRecoITSClusters(0),
  89         fCountRecoPPR(0),
  90         fCountRecoPID(0),
  91         fEvents(0),
  92         fDecayChannel(0),
  93         fFillFromGenerated(kFALSE),
  94         fOriginDselection(0),
  95         fAcceptanceUnf(kTRUE),
  96         fCuts(0),
  97         fUseWeight(kFALSE),
  98         fWeight(1.),
  99         fUseFlatPtWeight(kFALSE),
 100         fUseZWeight(kFALSE),
 101         fNvar(0),
 102         fPartName(""),
 103         fDauNames(""),
 104         fSign(2),
 105         fCentralitySelection(kTRUE),
 106         fFakeSelection(0),
 107         fRejectIfNoQuark(kTRUE),
 108         fUseMCVertex(kFALSE),
 109         fDsOption(1),
 110         fGenDsOption(3),
 111         fConfiguration(kCheetah), // by default, setting the fast configuration
 112         fFuncWeight(0x0)
 113 {
 114         //
 115         //Default ctor
 116         //
 117 }
 118 //___________________________________________________________________________
 119 AliCFTaskVertexingHF::AliCFTaskVertexingHF(const Char_t* name, AliRDHFCuts* cuts, TF1* func) :
 120         AliAnalysisTaskSE(name),
 121         fCFManager(0x0),
 122         fHistEventsProcessed(0x0),
 123         fCorrelation(0x0),
 124         fCountMC(0),
 125         fCountAcc(0),
 126         fCountVertex(0),
 127         fCountRefit(0),
 128         fCountReco(0),
 129         fCountRecoAcc(0),
 130         fCountRecoITSClusters(0),
 131         fCountRecoPPR(0),
 132         fCountRecoPID(0),
 133         fEvents(0),
 134         fDecayChannel(0),
 135         fFillFromGenerated(kFALSE),
 136         fOriginDselection(0),
 137         fAcceptanceUnf(kTRUE),
 138         fCuts(cuts),
 139         fUseWeight(kFALSE),
 140         fWeight(1.),
 141         fUseFlatPtWeight(kFALSE),
 142         fUseZWeight(kFALSE),
 143         fNvar(0),
 144         fPartName(""),
 145         fDauNames(""),
 146         fSign(2),
 147         fCentralitySelection(kTRUE),
 148         fFakeSelection(0),
 149         fRejectIfNoQuark(kTRUE),
 150         fUseMCVertex(kFALSE),
 151         fDsOption(1),
 152         fGenDsOption(3),
 153         fConfiguration(kCheetah),  // by default, setting the fast configuration
 154         fFuncWeight(func)
 155 {
 156         //
 157         // Constructor. Initialization of Inputs and Outputs
 158         //
 159         /*
 160           DefineInput(0) and DefineOutput(0)
 161           are taken care of by AliAnalysisTaskSE constructor
 162         */
 163         DefineOutput(1,TH1I::Class());
 164         DefineOutput(2,AliCFContainer::Class());
 165         DefineOutput(3,THnSparseD::Class());
 166         DefineOutput(4,AliRDHFCuts::Class());
 167
 168         fCuts->PrintAll();
 169 }
 170
 171 //___________________________________________________________________________
 172 AliCFTaskVertexingHF& AliCFTaskVertexingHF::operator=(const AliCFTaskVertexingHF& c)
 173 {
 174         //
 175         // Assignment operator
 176         //
 177         if (this!=&c) {
 178                 AliAnalysisTaskSE::operator=(c) ;
 179                 fCFManager  = c.fCFManager;
 180                 fHistEventsProcessed = c.fHistEventsProcessed;
 181                 fCuts = c.fCuts;
 182                 fFuncWeight = c.fFuncWeight;
 183         }
 184         return *this;
 185 }
 186
 187 //___________________________________________________________________________
 188 AliCFTaskVertexingHF::AliCFTaskVertexingHF(const AliCFTaskVertexingHF& c) :
 189         AliAnalysisTaskSE(c),
 190         fCFManager(c.fCFManager),
 191         fHistEventsProcessed(c.fHistEventsProcessed),
 192         fCorrelation(c.fCorrelation),
 193         fCountMC(c.fCountMC),
 194         fCountAcc(c.fCountAcc),
 195         fCountVertex(c.fCountVertex),
 196         fCountRefit(c.fCountRefit),
 197         fCountReco(c.fCountReco),
 198         fCountRecoAcc(c.fCountRecoAcc),
 199         fCountRecoITSClusters(c.fCountRecoITSClusters),
 200         fCountRecoPPR(c.fCountRecoPPR),
 201         fCountRecoPID(c.fCountRecoPID),
 202         fEvents(c.fEvents),
 203         fDecayChannel(c.fDecayChannel),
 204         fFillFromGenerated(c.fFillFromGenerated),
 205         fOriginDselection(c.fOriginDselection),
 206         fAcceptanceUnf(c.fAcceptanceUnf),
 207         fCuts(c.fCuts),
 208         fUseWeight(c.fUseWeight),
 209         fWeight(c.fWeight),
 210         fUseFlatPtWeight(c.fUseFlatPtWeight),
 211         fUseZWeight(c.fUseZWeight),
 212         fNvar(c.fNvar),
 213         fPartName(c.fPartName),
 214         fDauNames(c.fDauNames),
 215         fSign(c.fSign),
 216         fCentralitySelection(c.fCentralitySelection),
 217         fFakeSelection(c.fFakeSelection),
 218         fRejectIfNoQuark(c.fRejectIfNoQuark),
 219         fUseMCVertex(c.fUseMCVertex),
 220         fDsOption(c.fDsOption),
 221         fGenDsOption(c.fGenDsOption),
 222         fConfiguration(c.fConfiguration),
 223         fFuncWeight(c.fFuncWeight)
 224 {
 225         //
 226         // Copy Constructor
 227         //
 228 }
 229
 230 //___________________________________________________________________________
 231 AliCFTaskVertexingHF::~AliCFTaskVertexingHF()
 232 {
 233         //
 234         //destructor
 235         //
 236         if (fCFManager)           delete fCFManager ;
 237         if (fHistEventsProcessed) delete fHistEventsProcessed ;
 238         if (fCorrelation)         delete fCorrelation ;
 239         if (fCuts)                delete fCuts;
 240         if (fFuncWeight)                delete fFuncWeight;
 241 }
 242
 243 //_________________________________________________________________________-
 244 void AliCFTaskVertexingHF::Init()
 245 {
 246         //
 247         // Initialization
 248         //
 249
 250         if (fDebug>1) printf("AliCFTaskVertexingHF::Init()");
 251         if(fUseWeight && fUseZWeight) { AliFatal("Can not use at the same time pt and z-vtx weights, please choose"); return; }
 252         AliRDHFCuts *copyfCuts = 0x0;
 253         if (!fCuts){
 254                 AliFatal("No cuts defined - Exiting...");
 255                 return;
 256         }
 257
 258         switch (fDecayChannel){
 259         case 2:{
 260                 copyfCuts = new AliRDHFCutsD0toKpi(*(static_cast<AliRDHFCutsD0toKpi*>(fCuts)));
 261                 switch (fConfiguration) {
 262                 case kSnail:  // slow configuration: all variables in
 263                         fNvar = 16;
 264                         break;
 265                 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
 266                         fNvar = 8;
 267                         break;
 268                 }
 269                 fPartName="D0";
 270                 fDauNames="K+pi";
 271                 break;
 272         }
 273         case 21:{
 274                 copyfCuts = new AliRDHFCutsDStartoKpipi(*(static_cast<AliRDHFCutsDStartoKpipi*>(fCuts)));
 275                 switch (fConfiguration) {
 276                 case kSnail:  // slow configuration: all variables in
 277                         fNvar = 16;
 278                         break;
 279                 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
 280                         fNvar = 8;
 281                         break;
 282                 }
 283                 fPartName="Dstar";
 284                 fDauNames="K+pi+pi";
 285                 break;
 286         }
 287         case 31:{
 288                 copyfCuts = new AliRDHFCutsDplustoKpipi(*(static_cast<AliRDHFCutsDplustoKpipi*>(fCuts)));
 289                 switch (fConfiguration) {
 290                 case kSnail:  // slow configuration: all variables in
 291                         fNvar = 14;
 292                         break;
 293                 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
 294                         fNvar = 8;
 295                         break;
 296                 }
 297                 fPartName="Dplus";
 298                 fDauNames="K+pi+pi";
 299                 break;
 300         }
 301         case 32:{
 302                 copyfCuts = new AliRDHFCutsLctopKpi(*(static_cast<AliRDHFCutsLctopKpi*>(fCuts)));
 303                 switch (fConfiguration) {
 304                 case kSnail:  // slow configuration: all variables in
 305                         fNvar = 18;
 306                         break;
 307                 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
 308                         fNvar = 8;
 309                         break;
 310                 }
 311                 fPartName="Lambdac";
 312                 fDauNames="p+K+pi";
 313                 break;
 314         }
 315         case 33:{
 316                 copyfCuts = new AliRDHFCutsDstoKKpi(*(static_cast<AliRDHFCutsDstoKKpi*>(fCuts)));
 317                 switch (fConfiguration) {
 318                 case kSnail:  // slow configuration: all variables in
 319                         fNvar = 14;
 320                         break;
 321                 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
 322                         fNvar = 8;
 323                         break;
 324                 }
 325                 fPartName="Ds";
 326                 fDauNames="K+K+pi";
 327                 break;
 328         }
 329         case 4:{
 330                 copyfCuts = new AliRDHFCutsD0toKpipipi(*(static_cast<AliRDHFCutsD0toKpipipi*>(fCuts)));
 331                 switch (fConfiguration) {
 332                 case kSnail:  // slow configuration: all variables in
 333                         fNvar = 16;
 334                         break;
 335                 case kCheetah:// fast configuration: only pt_candidate, y, phi, ct, fake, z_vtx, centrality, multiplicity will be filled
 336                         fNvar = 8;
 337                         break;
 338                 }
 339                 fPartName="D0";
 340                 fDauNames="K+pi+pi+pi";
 341                 break;
 342         }
 343         default:
 344                 AliFatal("The decay channel MUST be defined according to AliCFVertexing::DecayChannel - Exiting...");
 345                 break;
 346         }
 347
 348         const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
 349         if (copyfCuts){
 350                 copyfCuts->SetName(nameoutput);
 351
 352                 //Post the data
 353                 PostData(4, copyfCuts);
 354         }
 355         else{
 356                 AliFatal("Failing initializing AliRDHFCuts object - Exiting...");
 357         }
 358
 359         return;
 360 }
 361
 362 //_________________________________________________
 363 void AliCFTaskVertexingHF::UserExec(Option_t *)
 364 {
 365         //
 366         // Main loop function
 367         //
 368
 369         PostData(1,fHistEventsProcessed) ;
 370         PostData(2,fCFManager->GetParticleContainer()) ;
 371         PostData(3,fCorrelation) ;
 372
 373         if (fFillFromGenerated){
 374                 AliWarning("Flag to fill container with generated value ON ---> dca, d0pi, d0K, d0xd0, cosPointingAngle will be set as dummy!");
 375         }
 376
 377         if (!fInputEvent) {
 378                 Error("UserExec","NO EVENT FOUND!");
 379                 return;
 380         }
 381
 382         AliAODEvent* aodEvent = dynamic_cast<AliAODEvent*>(fInputEvent);
 383
 384         TClonesArray *arrayBranch=0;
 385
 386         if(!aodEvent && AODEvent() && IsStandardAOD()) {
 387                 // In case there is an AOD handler writing a standard AOD, use the AOD
 388                 // event in memory rather than the input (ESD) event.
 389                 aodEvent = dynamic_cast<AliAODEvent*> (AODEvent());
 390                 // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
 391                 // have to taken from the AOD event hold by the AliAODExtension
 392                 AliAODHandler* aodHandler = (AliAODHandler*)
 393                         ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 394                 if(aodHandler->GetExtensions()) {
 395                         AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
 396                         AliAODEvent *aodFromExt = ext->GetAOD();
 397
 398                         switch (fDecayChannel){
 399                         case 2:{
 400                                 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("D0toKpi");
 401                                 break;
 402                         }
 403                         case 21:{
 404                                 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Dstar");
 405                                 break;
 406                         }
 407                         case 31:
 408                         case 32:
 409                         case 33:{
 410                                 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm3Prong");
 411                                 break;
 412                         }
 413                         case 4:{
 414                                 arrayBranch=(TClonesArray*)aodFromExt->GetList()->FindObject("Charm4Prong");
 415                                 break;
 416                         }
 417                         default:
 418                                 break;
 419                         }
 420                 }
 421         }
 422         else {
 423                 switch (fDecayChannel){
 424                 case 2:{
 425                         arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("D0toKpi");
 426                         break;
 427                 }
 428                 case 21:{
 429                         arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Dstar");
 430                         break;
 431                 }
 432                 case 31:
 433                 case 32:
 434                 case 33:{
 435                         arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm3Prong");
 436                         break;
 437                 }
 438                 case 4:{
 439                         arrayBranch=(TClonesArray*)aodEvent->GetList()->FindObject("Charm4Prong");
 440                         break;
 441                 }
 442                 default:
 443                         break;
 444                 }
 445         }
 446
 447         AliAODVertex *aodVtx = (AliAODVertex*)aodEvent->GetPrimaryVertex();
 448         if (!aodVtx) return;
 449
 450         if (!arrayBranch) {
 451                 AliError("Could not find array of HF vertices");
 452                 return;
 453         }
 454
 455         fEvents++;
 456
 457         fCFManager->SetRecEventInfo(aodEvent);
 458         fCFManager->SetMCEventInfo(aodEvent);
 459
 460         //******** DEFINE number of variables of the container***** for now set at 13, in the future in the config macro.
 461
 462         //loop on the MC event
 463
 464         TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 465         if (!mcArray) {
 466                 AliError("Could not find Monte-Carlo in AOD");
 467                 return;
 468         }
 469         Int_t icountMC = 0;
 470         Int_t icountAcc = 0;
 471         Int_t icountReco = 0;
 472         Int_t icountVertex = 0;
 473         Int_t icountRefit = 0;
 474         Int_t icountRecoAcc = 0;
 475         Int_t icountRecoITSClusters = 0;
 476         Int_t icountRecoPPR = 0;
 477         Int_t icountRecoPID = 0;
 478         Int_t cquarks = 0;
 479
 480         AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
 481         if (!mcHeader) {
 482                 AliError("Could not find MC Header in AOD");
 483                 return;
 484         }
 485
 486         Double_t* containerInput = new Double_t[fNvar];
 487         Double_t* containerInputMC = new Double_t[fNvar];
 488
 489
 490         AliCFVertexingHF* cfVtxHF=0x0;
 491         switch (fDecayChannel){
 492         case 2:{
 493           cfVtxHF = new AliCFVertexingHF2Prong(mcArray, fOriginDselection);
 494           break;
 495         }
 496         case 21:{
 497           cfVtxHF = new AliCFVertexingHFCascade(mcArray, fOriginDselection);
 498           break;
 499         }
 500         case 31:
 501         case 32:
 502         case 33:{
 503           cfVtxHF = new AliCFVertexingHF3Prong(mcArray, fOriginDselection, fDecayChannel);
 504           if(fDecayChannel==33){
 505             cfVtxHF->SetGeneratedDsOption(fGenDsOption);
 506           }
 507                 break;
 508         }
 509         case 4:{
 510                 //cfVtxHF = new AliCFVertexingHF4Prong(mcArray, originDselection);  // not there yet
 511                 break;
 512         }
 513         default:
 514                 break;
 515         }
 516         if (!cfVtxHF){
 517                 AliError("No AliCFVertexingHF initialized");
 518                 delete[] containerInput;
 519                 delete[] containerInputMC;
 520                 return;
 521         }
 522
 523         Double_t zPrimVertex = aodVtx ->GetZ();
 524         Double_t zMCVertex = mcHeader->GetVtxZ();
 525         Int_t runnumber = aodEvent->GetRunNumber();
 526         if(fUseZWeight) fWeight = GetZWeight(zMCVertex,runnumber);
 527
 528         if (TMath::Abs(zMCVertex) > fCuts->GetMaxVtxZ()){
 529           AliDebug(3,Form("z coordinate of MC vertex = %f, it was required to be within [-%f, +%f], skipping event", zMCVertex, fCuts->GetMaxVtxZ(), fCuts->GetMaxVtxZ()));
 530           delete[] containerInput;
 531           delete[] containerInputMC;
 532           return;
 533         }
 534
 535         AliESDtrackCuts** trackCuts = new AliESDtrackCuts*[cfVtxHF->GetNProngs()];
 536         if (fDecayChannel == 21){
 537                 // for the D*, setting the third element of the array of the track cuts to those for the soft pion
 538                 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs()-1; iProng++){
 539                         trackCuts[iProng]=fCuts->GetTrackCuts();
 540                 }
 541                 trackCuts[2] = fCuts->GetTrackCutsSoftPi();
 542         }
 543         else {
 544                 for (Int_t iProng = 0; iProng<cfVtxHF->GetNProngs(); iProng++){
 545                         trackCuts[iProng]=fCuts->GetTrackCuts();
 546                 }
 547         }
 548
 549         //General settings: vertex, feed down and fill reco container with generated values.
 550         cfVtxHF->SetRecoPrimVertex(zPrimVertex);
 551         cfVtxHF->SetMCPrimaryVertex(zMCVertex);
 552         cfVtxHF->SetFillFromGenerated(fFillFromGenerated);
 553         cfVtxHF->SetNVar(fNvar);
 554         cfVtxHF->SetFakeSelection(fFakeSelection);
 555         cfVtxHF->SetRejectCandidateIfNotFromQuark(fRejectIfNoQuark);
 556         cfVtxHF->SetConfiguration(fConfiguration);
 557
 558         // switch-off the trigger class selection (doesn't work for MC)
 559         fCuts->SetTriggerClass("");
 560
 561         // MC vertex, to be used, in case, for pp
 562         if (fUseMCVertex) fCuts->SetUseMCVertex();
 563
 564         if (fCentralitySelection){ // keep only the requested centrality
 565           if(fCuts->IsEventSelectedInCentrality(aodEvent)!=0) {
 566             delete[] containerInput;
 567             delete[] containerInputMC;
 568             delete [] trackCuts;
 569             return;
 570           }
 571         }  else { // keep all centralities
 572           fCuts->SetMinCentrality(0.);
 573           fCuts->SetMaxCentrality(100.);
 574         }
 575
 576         Float_t centValue = fCuts->GetCentrality(aodEvent);
 577         cfVtxHF->SetCentralityValue(centValue);
 578
 579         // number of tracklets - multiplicity estimator
 580         Double_t multiplicity = (Double_t)(aodEvent->GetTracklets()->GetNumberOfTracklets()); // casted to double because the CF is filled with doubles
 581         cfVtxHF->SetMultiplicity(multiplicity);
 582
 583         for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
 584           AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(iPart));
 585           if (!mcPart){
 586             AliError("Failed casting particle from MC array!, Skipping particle");
 587             continue;
 588           }
 589           // check the MC-level cuts, must be the desidered particle
 590           if (!fCFManager->CheckParticleCuts(0, mcPart)) {
 591             continue;  // 0 stands for MC level
 592           }
 593           cfVtxHF->SetMCCandidateParam(iPart);
 594
 595           //counting c quarks
 596           cquarks += cfVtxHF->MCcquarkCounting(mcPart);
 597
 598           if (!(cfVtxHF->SetLabelArray())){
 599             AliDebug(2,Form("Impossible to set the label array (decaychannel = %d)",fDecayChannel));
 600             continue;
 601           }
 602
 603           //check the candiate family at MC level
 604           if (!(cfVtxHF->CheckMCPartFamily(mcPart, mcArray))) {
 605             AliDebug(2,Form("Check on the family wrong!!! (decaychannel = %d)",fDecayChannel));
 606             continue;
 607           }
 608                 else{
 609                         AliDebug(2,Form("Check on the family OK!!! (decaychannel = %d)",fDecayChannel));
 610                 }
 611
 612                 //Fill the MC container
 613                 Bool_t mcContainerFilled = cfVtxHF -> FillMCContainer(containerInputMC);
 614                 if (mcContainerFilled) {
 615                         if (fUseWeight){
 616                                 if (fFuncWeight){ // using user-defined function
 617                                         AliDebug(2,"Using function");
 618                                         fWeight = fFuncWeight->Eval(containerInputMC[0]);
 619                                 }
 620                                 else{ // using FONLL
 621                                         AliDebug(2,"Using FONLL");
 622                                         fWeight = GetWeight(containerInputMC[0]);
 623                                 }
 624                                 AliDebug(2,Form("pt = %f, weight = %f",containerInputMC[0], fWeight));
 625                         }
 626                         if (!fCuts->IsInFiducialAcceptance(containerInputMC[0],containerInputMC[1])) continue;
 627                         //MC Limited Acceptance
 628                         if (TMath::Abs(containerInputMC[1]) < 0.5) {
 629                                 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepGeneratedLimAcc, fWeight);
 630                                 AliDebug(3,"MC Lim Acc container filled\n");
 631                         }
 632
 633                         //MC
 634                         fCFManager->GetParticleContainer()->Fill(containerInputMC, kStepGenerated, fWeight);
 635                         icountMC++;
 636                         AliDebug(3,"MC cointainer filled \n");
 637
 638                         // MC in acceptance
 639                         // check the MC-Acceptance level cuts
 640                         // since standard CF functions are not applicable, using Kine Cuts on daughters
 641                         Bool_t mcAccepStep = cfVtxHF-> MCAcceptanceStep();
 642                         if (mcAccepStep){
 643                                 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepAcceptance, fWeight);
 644                                 AliDebug(3,"MC acceptance cut passed\n");
 645                                 icountAcc++;
 646
 647                                 //MC Vertex step
 648                                 if (fCuts->IsEventSelected(aodEvent)){
 649                                         // filling the container if the vertex is ok
 650                                         fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepVertex, fWeight) ;
 651                                         AliDebug(3,"Vertex cut passed and container filled\n");
 652                                         icountVertex++;
 653
 654                                         //mc Refit requirement
 655                                         Bool_t mcRefitStep = cfVtxHF->MCRefitStep(aodEvent, &trackCuts[0]);
 656                                         if (mcRefitStep){
 657                                                 fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepRefit, fWeight);
 658                                                 AliDebug(3,"MC Refit cut passed and container filled\n");
 659                                                 icountRefit++;
 660                                         }
 661                                         else{
 662                                                 AliDebug(3,"MC Refit cut not passed\n");
 663                                                 continue;
 664                                         }
 665                                 }
 666                                 else{
 667                                   AliDebug (3, "MC vertex step not passed\n");
 668                                   continue;
 669                                 }
 670                         }
 671                         else{
 672                                 AliDebug (3, "MC in acceptance step not passed\n");
 673                                 continue;
 674                         }
 675                 }
 676                 else {
 677                         AliDebug (3, "MC container not filled\n");
 678                 }
 679         }
 680
 681         if (cquarks<2) AliDebug(2,Form("Event with %d c-quarks", cquarks));
 682         AliDebug(2,Form("Found %i MC particles that are %s!!",icountMC,fPartName.Data()));
 683         AliDebug(2,Form("Found %i MC particles that are %s and satisfy Acc cuts!!",icountAcc,fPartName.Data()));
 684         AliDebug(2,Form("Found %i MC particles that are %s and satisfy Vertex cuts!!",icountVertex,fPartName.Data()));
 685         AliDebug(2,Form("Found %i MC particles that are %s and satisfy Refit cuts!!",icountRefit,fPartName.Data()));
 686
 687         // Now go to rec level
 688         fCountMC += icountMC;
 689         fCountAcc += icountAcc;
 690         fCountVertex+= icountVertex;
 691         fCountRefit+= icountRefit;
 692
 693         AliDebug(2,Form("Found %d vertices for decay channel %d",arrayBranch->GetEntriesFast(),fDecayChannel));
 694
 695         for(Int_t iCandid = 0; iCandid<arrayBranch->GetEntriesFast();iCandid++){
 696                 AliAODRecoDecayHF* charmCandidate=0x0;
 697                 switch (fDecayChannel){
 698                 case 2:{
 699                         charmCandidate = (AliAODRecoDecayHF2Prong*)arrayBranch->At(iCandid);
 700                         break;
 701                 }
 702                 case 21:{
 703                         charmCandidate = (AliAODRecoCascadeHF*)arrayBranch->At(iCandid);
 704                         break;
 705                 }
 706                 case 31:
 707                 case 32:
 708                 case 33:{
 709                         charmCandidate = (AliAODRecoDecayHF3Prong*)arrayBranch->At(iCandid);
 710                         break;
 711                 }
 712                 case 4:{
 713                         charmCandidate = (AliAODRecoDecayHF4Prong*)arrayBranch->At(iCandid);
 714                         break;
 715                 }
 716                 default:
 717                         break;
 718                 }
 719
 720                 Bool_t unsetvtx=kFALSE;
 721                 if(!charmCandidate->GetOwnPrimaryVtx()) {
 722                         charmCandidate->SetOwnPrimaryVtx(aodVtx); // needed to compute all variables
 723                         unsetvtx=kTRUE;
 724                 }
 725
 726                 Bool_t signAssociation = cfVtxHF->SetRecoCandidateParam((AliAODRecoDecayHF*)charmCandidate);
 727                 if (!signAssociation){
 728                         charmCandidate = 0x0;
 729                         continue;
 730                 }
 731
 732                 Int_t isPartOrAntipart = cfVtxHF->CheckReflexion(fSign);
 733                 if (isPartOrAntipart == 0){
 734                         AliDebug(2, Form("The candidate pdg code doesn't match the requirement set in the task (fSign = %d)",fSign));
 735                         continue;
 736                 }
 737
 738
 739                 Bool_t recoContFilled = cfVtxHF->FillRecoContainer(containerInput);
 740                 if (recoContFilled){
 741
 742                         // weight according to pt
 743                         if (fUseWeight){
 744                                 if (fFuncWeight){ // using user-defined function
 745                                         AliDebug(2, "Using function");
 746                                         fWeight = fFuncWeight->Eval(containerInput[0]);
 747                                 }
 748                                 else{ // using FONLL
 749                                         AliDebug(2, "Using FONLL");
 750                                         fWeight = GetWeight(containerInput[0]);
 751                                 }
 752                                 AliDebug(2, Form("pt = %f, weight = %f",containerInput[0], fWeight));
 753                         }
 754
 755                         if (!fCuts->IsInFiducialAcceptance(containerInput[0],containerInput[1])) continue;
 756
 757                         //Reco Step
 758                         Bool_t recoStep = cfVtxHF->RecoStep();
 759                         Bool_t vtxCheck = fCuts->IsEventSelected(aodEvent);
 760
 761                         if (recoStep && recoContFilled && vtxCheck){
 762                                 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed, fWeight) ;
 763                                 icountReco++;
 764                                 AliDebug(3,"Reco step  passed and container filled\n");
 765
 766                                 //Reco in the acceptance -- take care of UNFOLDING!!!!
 767                                 Bool_t recoAcceptanceStep = cfVtxHF->RecoAcceptStep(&trackCuts[0]);
 768                                 if (recoAcceptanceStep) {
 769                                         fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoAcceptance, fWeight) ;
 770                                         icountRecoAcc++;
 771                                         AliDebug(3,"Reco acceptance cut passed and container filled\n");
 772
 773                                         if(fAcceptanceUnf){
 774                                                 Double_t fill[4]; //fill response matrix
 775                                                 Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
 776                                                 if (bUnfolding) fCorrelation->Fill(fill);
 777                                         }
 778
 779                                         //Number of ITS cluster requirements
 780                                         Int_t recoITSnCluster = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kTracks);
 781                                         if (recoITSnCluster){
 782                                                 fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoITSClusters, fWeight) ;
 783                                                 icountRecoITSClusters++;
 784                                                 AliDebug(3,"Reco n ITS cluster cut passed and container filled\n");
 785
 786                                                 Bool_t iscutsusingpid = fCuts->GetIsUsePID();
 787                                                 Int_t recoAnalysisCuts = -1, recoPidSelection = -1;
 788                                                 fCuts->SetUsePID(kFALSE);
 789                                                 recoAnalysisCuts = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
 790
 791                                                 if (fDecayChannel==33){ // Ds case, where more possibilities are considered
 792                                                   Bool_t keepDs=ProcessDs(recoAnalysisCuts);
 793                                                   if(keepDs) recoAnalysisCuts=3;
 794                                                 }
 795
 796
 797                                                 fCuts->SetUsePID(iscutsusingpid); //restoring usage of the PID from the cuts object
 798                                                 if (recoAnalysisCuts == 3 || recoAnalysisCuts == isPartOrAntipart){
 799                                                         fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPPR, fWeight);
 800                                                         icountRecoPPR++;
 801                                                         AliDebug(3,"Reco Analysis cuts passed and container filled \n");
 802                                                         //pid selection
 803                                                         //recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kPID);
 804                                                         //if((fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==isPartOrAntipart)||(fCuts->CombineSelectionLevels(3,recoAnalysisCuts,recoPidSelection)==3)){
 805                                                         recoPidSelection = fCuts->IsSelected(charmCandidate, AliRDHFCuts::kCandidate, aodEvent);
 806
 807                                                         if (fDecayChannel==33){ // Ds case, where more possibilities are considered
 808                                                           Bool_t keepDs=ProcessDs(recoPidSelection);
 809                                                           if(keepDs) recoPidSelection=3;
 810                                                         }
 811
 812                                                         if (recoPidSelection == 3 || recoPidSelection == isPartOrAntipart){
 813                                                                 fCFManager->GetParticleContainer()->Fill(containerInput, kStepRecoPID, fWeight);
 814                                                                 icountRecoPID++;
 815                                                                 AliDebug(3,"Reco PID cuts passed and container filled \n");
 816                                                                 if(!fAcceptanceUnf){
 817                                                                         Double_t fill[4]; //fill response matrix
 818                                                                         Bool_t bUnfolding = cfVtxHF -> FillUnfoldingMatrix(fill);
 819                                                                         if (bUnfolding) fCorrelation->Fill(fill);
 820                                                                 }
 821                                                         }
 822                                                         else {
 823                                                                 AliDebug(3, "Analysis Cuts step not passed \n");
 824                                                                 continue;
 825                                                         }
 826                                                 }
 827                                                 else {
 828                                                         AliDebug(3, "PID selection not passed \n");
 829                                                         continue;
 830                                                 }
 831                                         }
 832                                         else{
 833                                                 AliDebug(3, "Number of ITS cluster step not passed\n");
 834                                                 continue;
 835                                         }
 836                                 }
 837                                 else{
 838                                         AliDebug(3, "Reco acceptance step not passed\n");
 839                                         continue;
 840                                 }
 841                         }
 842                         else {
 843                                 AliDebug(3, "Reco step not passed\n");
 844                                 continue;
 845                         }
 846                 }
 847
 848                 if(unsetvtx) charmCandidate->UnsetOwnPrimaryVtx();
 849         } // end loop on candidate
 850
 851         fCountReco+= icountReco;
 852         fCountRecoAcc+= icountRecoAcc;
 853         fCountRecoITSClusters+= icountRecoITSClusters;
 854         fCountRecoPPR+= icountRecoPPR;
 855         fCountRecoPID+= icountRecoPID;
 856
 857         fHistEventsProcessed->Fill(0);
 858
 859         delete[] containerInput;
 860         delete[] containerInputMC;
 861         delete cfVtxHF;
 862         if (trackCuts){
 863         //      for (Int_t i=0; i<cfVtxHF->GetNProngs(); i++){
 864         //              delete [] trackCuts[i];
 865         //      }
 866         delete [] trackCuts;
 867         }
 868
 869
 870 }
 871
 872 //___________________________________________________________________________
 873 void AliCFTaskVertexingHF::Terminate(Option_t*)
 874 {
 875         // The Terminate() function is the last function to be called during
 876         // a query. It always runs on the client, it can be used to present
 877         // the results graphically or save the results to file.
 878
 879         AliAnalysisTaskSE::Terminate();
 880
 881         AliInfo(Form("Found %i MC particles that are %s in MC, in %d events",fCountMC,fPartName.Data(),fEvents));
 882         AliInfo(Form("Found %i MC particles that are %s in MC and satisfy Acc cuts, in %d events",fCountAcc,fPartName.Data(),fEvents));
 883         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));
 884         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));
 885         AliInfo(Form("Found %i reco %s that are decaying in %s, in %d events",fCountReco,fPartName.Data(),fDauNames.Data(),fEvents));
 886         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));
 887         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));
 888         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));
 889         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));
 890
 891         // draw some example plots....
 892         AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
 893         if(!cont) {
 894                 printf("CONTAINER NOT FOUND\n");
 895                 return;
 896         }
 897         // projecting the containers to obtain histograms
 898         // first argument = variable, second argument = step
 899
 900         TH1D** h = new TH1D*[3];
 901         if (fConfiguration == kSnail){
 902                 //h = new TH1D[3][12];
 903                 for (Int_t ih = 0; ih<3; ih++){
 904                         h[ih] = new TH1D[12];
 905                 }
 906                 for(Int_t iC=1;iC<12; iC++){
 907                         // MC-level
 908                         h[0][iC] =   *(cont->ShowProjection(iC,0));
 909                         // MC-Acceptance level
 910                         h[1][iC] =   *(cont->ShowProjection(iC,1));
 911                         // Reco-level
 912                         h[2][iC] =   *(cont->ShowProjection(iC,4));
 913                 }
 914         }
 915         else{
 916                 //h = new TH1D[3][12];
 917                 for (Int_t ih = 0; ih<3; ih++){
 918                         h[ih] = new TH1D[8];
 919                 }
 920                 for(Int_t iC=0;iC<8; iC++){
 921                         // MC-level
 922                         h[0][iC] =   *(cont->ShowProjection(iC,0));
 923                         // MC-Acceptance level
 924                         h[1][iC] =   *(cont->ShowProjection(iC,1));
 925                         // Reco-level
 926                         h[2][iC] =   *(cont->ShowProjection(iC,4));
 927                 }
 928         }
 929         TString* titles;
 930         Int_t nvarToPlot = 0;
 931         if (fConfiguration == kSnail){
 932                 nvarToPlot = 12;
 933                 titles = new TString[nvarToPlot];
 934                 if(fDecayChannel==31){
 935                         titles[0]="pT_Dplus (GeV/c)";
 936                         titles[1]="rapidity";
 937                         titles[2]="phi (rad)";
 938                         titles[3]="cT (#mum)";
 939                         titles[4]="cosPointingAngle";
 940                         titles[5]="pT_1 (GeV/c)";
 941                         titles[6]="pT_2 (GeV/c)";
 942                         titles[7]="pT_3 (GeV/c)";
 943                         titles[8]="d0_1 (#mum)";
 944                         titles[9]="d0_2 (#mum)";
 945                         titles[10]="d0_3 (#mum)";
 946                         titles[11]="zVertex (cm)";
 947                 }else{
 948                         titles[0]="pT_D0 (GeV/c)";
 949                         titles[1]="rapidity";
 950                         titles[2]="cosThetaStar";
 951                         titles[3]="pT_pi (GeV/c)";
 952                         titles[4]="pT_K (Gev/c)";
 953                         titles[5]="cT (#mum)";
 954                         titles[6]="dca (#mum)";
 955                         titles[7]="d0_pi (#mum)";
 956                         titles[8]="d0_K (#mum)";
 957                         titles[9]="d0xd0 (#mum^2)";
 958                         titles[10]="cosPointingAngle";
 959                         titles[11]="phi (rad)";
 960
 961                 }
 962         }
 963         else{
 964                 nvarToPlot = 8;
 965                 titles = new TString[nvarToPlot];
 966                 titles[0]="pT_candidate (GeV/c)";
 967                 titles[1]="rapidity";
 968                 titles[2]="cT (#mum)";
 969                 titles[3]="phi";
 970                 titles[4]="z_{vtx}";
 971                 titles[5]="centrality";
 972                 titles[6]="fake";
 973                 titles[7]="multiplicity";
 974         }
 975
 976         Int_t markers[12]={20,24,21,25,27,28,
 977                            20,24,21,25,27,28};
 978         Int_t colors[3]={2,8,4};
 979         for(Int_t iC=0;iC<nvarToPlot; iC++){
 980                 for(Int_t iStep=0;iStep<3;iStep++){
 981                         h[iStep][iC].SetTitle(titles[iC].Data());
 982                         h[iStep][iC].GetXaxis()->SetTitle(titles[iC].Data());
 983                         Double_t maxh=h[iStep][iC].GetMaximum();
 984                         h[iStep][iC].GetYaxis()->SetRangeUser(0,maxh*1.2);
 985                         h[iStep][iC].SetMarkerStyle(markers[iC]);
 986                         h[iStep][iC].SetMarkerColor(colors[iStep]);
 987                 }
 988         }
 989
 990         gStyle->SetCanvasColor(0);
 991         gStyle->SetFrameFillColor(0);
 992         gStyle->SetTitleFillColor(0);
 993         gStyle->SetStatColor(0);
 994
 995         // drawing in 2 separate canvas for a matter of clearity
 996         TCanvas * c1 =new TCanvas(Form("c1New_%d",fDecayChannel),"Vars 0, 1, 2, 3",1100,1200);
 997         c1->Divide(3,4);
 998         Int_t iPad=1;
 999         for(Int_t iVar=0; iVar<4; iVar++){
1000                 c1->cd(iPad++);
1001                 h[0][iVar].DrawCopy("p");
1002                 c1->cd(iPad++);
1003                 h[1][iVar].DrawCopy("p");
1004                 c1->cd(iPad++);
1005                 h[2][iVar].DrawCopy("p");
1006         }
1007
1008         TCanvas * c2 =new TCanvas(Form("c2New_%d",fDecayChannel),"Vars 4, 5, 6, 7",1100,1200);
1009         c2->Divide(3,4);
1010         iPad=1;
1011         for(Int_t iVar=4; iVar<8; iVar++){
1012                 c2->cd(iPad++);
1013                 h[0][iVar].DrawCopy("p");
1014                 c2->cd(iPad++);
1015                 h[1][iVar].DrawCopy("p");
1016                 c2->cd(iPad++);
1017                 h[2][iVar].DrawCopy("p");
1018         }
1019
1020         if (fConfiguration == kSnail){
1021                 TCanvas * c3 =new TCanvas(Form("c3New_%d",fDecayChannel),"Vars 8, 9, 10, 11",1100,1200);
1022                 c3->Divide(3,4);
1023                 iPad=1;
1024                 for(Int_t iVar=8; iVar<12; iVar++){
1025                         c3->cd(iPad++);
1026                         h[0][iVar].DrawCopy("p");
1027                         c3->cd(iPad++);
1028                         h[1][iVar].DrawCopy("p");
1029                         c3->cd(iPad++);
1030                         h[2][iVar].DrawCopy("p");
1031                 }
1032         }
1033
1034
1035         THnSparseD* hcorr = dynamic_cast<THnSparseD*> (GetOutputData(3));
1036
1037         TH2D* corr1 =hcorr->Projection(0,2);
1038         TH2D* corr2 = hcorr->Projection(1,3);
1039
1040         TCanvas * c7 =new TCanvas(Form("c7New_%d",fDecayChannel),"",800,400);
1041         c7->Divide(2,1);
1042         c7->cd(1);
1043         corr1->DrawCopy("text");
1044         c7->cd(2);
1045         corr2->DrawCopy("text");
1046
1047         TFile* file_projection = new TFile("CFtaskHFprojectionNew.root","RECREATE");
1048
1049         corr1->Write();
1050         corr2->Write();
1051         for(Int_t iC=0;iC<nvarToPlot; iC++){
1052                 for(Int_t iStep=0;iStep<3;iStep++){
1053                         h[iStep][iC].Write(Form("Step%d_%s",iStep,titles[iC].Data()));
1054                 }
1055         }
1056         file_projection->Close();
1057         for (Int_t ih = 0; ih<3; ih++) delete [] h[ih];
1058         delete [] h;
1059
1060
1061 }
1062
1063 //___________________________________________________________________________
1064 void AliCFTaskVertexingHF::UserCreateOutputObjects()
1065 {
1066         //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED
1067         //TO BE SET BEFORE THE EXECUTION OF THE TASK
1068         //
1069         Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName());
1070
1071         //slot #1
1072         OpenFile(1);
1073         fHistEventsProcessed = new TH1I("CFHFchist0","",1,0,1) ;
1074
1075         PostData(1,fHistEventsProcessed) ;
1076         PostData(2,fCFManager->GetParticleContainer()) ;
1077         PostData(3,fCorrelation) ;
1078
1079 }
1080
1081
1082 //_________________________________________________________________________
1083 Double_t AliCFTaskVertexingHF::GetWeight(Float_t pt)
1084 {
1085         //
1086         // calculating the weight to fill the container
1087         //
1088
1089         // FNOLL central:
1090         // p0 = 1.63297e-01 --> 0.322643
1091         // p1 = 2.96275e+00
1092         // p2 = 2.30301e+00
1093         // p3 = 2.50000e+00
1094
1095         // PYTHIA
1096         // p0 = 1.85906e-01 --> 0.36609
1097         // p1 = 1.94635e+00
1098         // p2 = 1.40463e+00
1099         // p3 = 2.50000e+00
1100
1101         Double_t func1[4] = {0.322643,2.96275,2.30301,2.5};
1102         Double_t func2[4] = {0.36609,1.94635,1.40463,2.5};
1103
1104         Double_t dndpt_func1 = dNdptFit(pt,func1);
1105         if(fUseFlatPtWeight) dndpt_func1 = 1./30.;
1106         Double_t dndpt_func2 = dNdptFit(pt,func2);
1107         AliDebug(2,Form("pt = %f, FONLL = %f, Pythia = %f, ratio = %f",pt,dndpt_func1,dndpt_func2,dndpt_func1/dndpt_func2));
1108         return dndpt_func1/dndpt_func2;
1109 }
1110
1111 //__________________________________________________________________________________________________
1112 Double_t AliCFTaskVertexingHF::dNdptFit(Float_t pt, Double_t* par)
1113 {
1114         //
1115         // calculating dNdpt
1116         //
1117
1118         Double_t denom =  TMath::Power((pt/par[1]), par[3] );
1119         Double_t dNdpt = par[0]*pt/TMath::Power(1.+denom, par[2]);
1120
1121         return dNdpt;
1122 }
1123
1124 //__________________________________________________________________________________________________
1125 Double_t AliCFTaskVertexingHF::GetZWeight(Float_t z, Int_t runnumber){
1126   //
1127   //  calculating the z-vtx weight for the given run range
1128   //
1129
1130   if(runnumber>146824 || runnumber<146803) return 1.0;
1131
1132   Double_t func1[3] = {1.0, -0.5, 6.5 };
1133   Double_t func2[3] = {1.0, -0.5, 5.5 };
1134
1135   Double_t dzFunc1 = DodzFit(z,func1);
1136   Double_t dzFunc2 = DodzFit(z,func2);
1137
1138   return dzFunc1/dzFunc2;
1139 }
1140
1141 //__________________________________________________________________________________________________
1142 Double_t AliCFTaskVertexingHF::DodzFit(Float_t z, Double_t* par) {
1143
1144   //
1145   // Gaussian z-vtx shape
1146   //
1147   //gaussian = [0]/TMath::Sqrt(2.*TMath::Pi())/[2]*exp[-(x-[1])*(x-[1])/(2*[2]*[2])]
1148
1149   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]);
1150
1151   return value;
1152 }
1153
1154 //__________________________________________________________________________________________________
1155 Bool_t AliCFTaskVertexingHF::ProcessDs(Int_t recoAnalysisCuts) const{
1156   // processes output of Ds is selected
1157   Bool_t keep=kFALSE;
1158   if(recoAnalysisCuts > 0){
1159     Int_t isKKpi=recoAnalysisCuts&1;
1160     Int_t ispiKK=recoAnalysisCuts&2;
1161     Int_t isPhiKKpi=recoAnalysisCuts&4;
1162     Int_t isPhipiKK=recoAnalysisCuts&8;
1163     Int_t isK0starKKpi=recoAnalysisCuts&16;
1164     Int_t isK0starpiKK=recoAnalysisCuts&32;
1165     if(fDsOption==1){
1166       if(isKKpi && isPhiKKpi) keep=kTRUE;
1167       if(ispiKK && isPhipiKK) keep=kTRUE;
1168     }
1169     else if(fDsOption==2){
1170       if(isKKpi && isK0starKKpi) keep=kTRUE;
1171       if(ispiKK && isK0starpiKK) keep=kTRUE;
1172     }
1173     else if(fDsOption==3)keep=kTRUE;
1174   }
1175   return keep;
1176 }