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