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