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