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