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