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