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