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