PWG3/vertexingHF/AliCFTaskVertexingHF.cxx

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