From: dainese Date: Thu, 19 Feb 2009 01:27:30 +0000 (+0000) Subject: New class for computing corrections as a function of several variables (Chiara) X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=becfe9890c4906d378cac808d765ff68c3d8c014;p=u%2Fmrichter%2FAliRoot.git New class for computing corrections as a function of several variables (Chiara) --- diff --git a/PWG3/CMake_libPWG3vertexingHF.txt b/PWG3/CMake_libPWG3vertexingHF.txt index 7564d55929a..adfa9954607 100644 --- a/PWG3/CMake_libPWG3vertexingHF.txt +++ b/PWG3/CMake_libPWG3vertexingHF.txt @@ -10,6 +10,7 @@ set(SRCS vertexingHF/AliAnalysisTaskSESelectHF.cxx vertexingHF/AliAnalysisTaskSECompareHF.cxx vertexingHF/AliCFHeavyFlavourTask.cxx + vertexingHF/AliCFHeavyFlavourTaskMultiVar.cxx vertexingHF/AliMultiDimVector.cxx vertexingHF/AliSignificanceCalculator.cxx ) diff --git a/PWG3/PWG3vertexingHFLinkDef.h b/PWG3/PWG3vertexingHFLinkDef.h index 3f85d2b8292..07bb43764ee 100644 --- a/PWG3/PWG3vertexingHFLinkDef.h +++ b/PWG3/PWG3vertexingHFLinkDef.h @@ -15,6 +15,7 @@ #pragma link C++ class AliAnalysisTaskSESelectHF+; #pragma link C++ class AliAnalysisTaskSECompareHF+; #pragma link C++ class AliCFHeavyFlavourTask+; +#pragma link C++ class AliCFHeavyFlavourTaskMultiVar+; #pragma link C++ class AliMultiDimVector+; #pragma link C++ class AliSignificanceCalculator+; diff --git a/PWG3/libPWG3vertexingHF.pkg b/PWG3/libPWG3vertexingHF.pkg index 5101a10f092..2b7178fe32f 100644 --- a/PWG3/libPWG3vertexingHF.pkg +++ b/PWG3/libPWG3vertexingHF.pkg @@ -9,6 +9,7 @@ SRCS:= vertexingHF/AliAODRecoDecayHF.cxx \ vertexingHF/AliAnalysisTaskSESelectHF.cxx \ vertexingHF/AliAnalysisTaskSECompareHF.cxx \ vertexingHF/AliCFHeavyFlavourTask.cxx \ + vertexingHF/AliCFHeavyFlavourTaskMultiVar.cxx \ vertexingHF/AliMultiDimVector.cxx vertexingHF/AliSignificanceCalculator.cxx HDRS:= $(SRCS:.cxx=.h) diff --git a/PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVar.C b/PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVar.C new file mode 100644 index 00000000000..8545a0053bd --- /dev/null +++ b/PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVar.C @@ -0,0 +1,294 @@ +//DEFINITION OF A FEW CONSTANTS +const Double_t ymin = -2.0 ; +const Double_t ymax = 2.0 ; +const Double_t ptmin_0_4 = 0.0 ; +const Double_t ptmax_0_4 = 4.0 ; +const Double_t ptmin_4_8 = 4.0 ; +const Double_t ptmax_4_8 = 8.0 ; +const Double_t ptmin_8_10 = 8.0 ; +const Double_t ptmax_8_10 = 10.0 ; +const Double_t cosmin = -1.; +const Double_t cosmax = 1.; +const Double_t cTmin = 0; // micron +const Double_t cTmax = 500; // micron +const Int_t mintrackrefsTPC = 2 ; +const Int_t mintrackrefsITS = 3 ; +const Int_t charge = 1 ; +const Int_t PDG = 421; +const Int_t minclustersTPC = 50 ; +//---------------------------------------------------- + +Bool_t AliCFHeavyFlavourTaskMultiVar() +{ + + TBenchmark benchmark; + benchmark.Start("AliCFHeavyFlavourTaskMultiVar"); + + AliLog::SetGlobalDebugLevel(0); + + Load() ; //load the required libraries + + + TChain * analysisChain ; + TChain * analysisChainFriend ; + + //here put your input data path + printf("\n\nRunning on local file, please check the path\n\n"); + analysisChain = new TChain("aodTree"); + analysisChainFriend = new TChain("aodTree"); + for (Int_t i = 1; i<2; i++){ + if (i != 4 && i != 5 && i != 6 && i != 10 && i != 11 && i != 14 && i != 15 && i != 16 && i != 18 && i != 20 && i != 21 && i != 24 && i != 26 && (i <30 || i > 43) && i!= 46 && i!=48 && i!=49 && i!=51 && i!=53 && i!=56 && i!=64 && i!=73 && i!=76 && i!=78 && i!=83 && i!=84 && i!=85 && i!=90 && i!=91 && i!=92){ + + TString fileName = "/home/zampolli/HF/FromRenu/Events/180002/"; + TString aodHFFileName = "/home/zampolli/HF/FromRenu/Events/180002/"; + fileName+=Form("%i/AliAOD.root",i); + printf(Form("file = %s \n", fileName.Data())); + aodHFFileName+=Form("%i/AliAOD.VertexingHF.root",i); + analysisChain->Add(fileName); + analysisChainFriend->Add(aodHFFileName); + } + } + + analysisChain->AddFriend(analysisChainFriend); + Info("AliCFHeavyFlavourTaskMultiVar",Form("CHAIN HAS %d ENTRIES",(Int_t)analysisChain->GetEntries())); + + //CONTAINER DEFINITION + Info("AliCFHeavyFlavourTaskMultiVar","SETUP CONTAINER"); + //the sensitive variables (6 in this example), their indices + UInt_t ipt = 0; + UInt_t iy = 1; + UInt_t icosThetaStar = 2; + UInt_t ipTpi = 3; + UInt_t ipTk = 4; + UInt_t icT = 5; + //Setting up the container grid... + UInt_t nstep = 2 ; //number of selection steps MC + const Int_t nvar = 6 ; //number of variables on the grid:pt, y, cosThetaStar, pTpi, pTk, cT + const Int_t nbin0_0_4 = 8 ; //bins in pt from 0 to 4 GeV + const Int_t nbin0_4_8 = 4 ; //bins in pt from 4 to 8 GeV + const Int_t nbin0_8_10 = 1 ; //bins in pt from 8 to 10 GeV + const Int_t nbin1 = 8 ; //bins in y + const Int_t nbin2 = 8 ; //bins in cosThetaStar + const Int_t nbin3_0_4 = 8 ; //bins in ptPi from 0 to 4 GeV + const Int_t nbin3_4_8 = 4 ; //bins in ptPi from 4 to 8 GeV + const Int_t nbin3_8_10 = 1 ; //bins in ptPi from 8 to 10 GeV + const Int_t nbin4_0_4 = 8 ; //bins in ptKa from 0 to 4 GeV + const Int_t nbin4_4_8 = 4 ; //bins in ptKa from 4 to 8 GeV + const Int_t nbin4_8_10 = 1 ; //bins in ptKa from 8 to 10 GeV + const Int_t nbin5 = 24 ; //bins in cT + + //arrays for the number of bins in each dimension + Int_t iBin[nvar]; + iBin[0]=nbin0_0_4+nbin0_4_8+nbin0_8_10; + iBin[1]=nbin1; + iBin[2]=nbin2; + iBin[3]=nbin3_0_4+nbin3_4_8+nbin3_8_10; + iBin[4]=nbin4_0_4+nbin4_4_8+nbin4_8_10; + iBin[5]=nbin5; + + //arrays for lower bounds : + Double_t *binLim0=new Double_t[iBin[0]+1]; + Double_t *binLim1=new Double_t[iBin[1]+1]; + Double_t *binLim2=new Double_t[iBin[2]+1]; + Double_t *binLim3=new Double_t[iBin[3]+1]; + Double_t *binLim4=new Double_t[iBin[4]+1]; + Double_t *binLim5=new Double_t[iBin[5]+1]; + + // checking limits + if (ptmax_0_4 != ptmin_4_8) { + Error("AliCFHeavyFlavourTaskMultiVar","max lim 1st range != min lim 2nd range, please check!"); + } + if (ptmax_4_8 != ptmin_8_10) { + Error("AliCFHeavyFlavourTaskMultiVar","max lim 2nd range != min lim 3rd range, please check!"); + } + + // values for bin lower bounds + // pt + for(Int_t i=0; i<=nbin0_0_4; i++) binLim0[i]=(Double_t)ptmin_0_4 + (ptmax_0_4-ptmin_0_4)/nbin0_0_4*(Double_t)i ; + if (binLim0[nbin0_0_4] != ptmin_4_8) { + Error("AliCFHeavyFlavourTaskMultiVar","Calculated bin lim for pt - 1st range - differs from expected!\n"); + } + for(Int_t i=0; i<=nbin0_4_8; i++) binLim0[i+nbin0_0_4]=(Double_t)ptmin_4_8 + (ptmax_4_8-ptmin_4_8)/nbin0_4_8*(Double_t)i ; + if (binLim0[nbin0_0_4+nbin0_4_8] != ptmin_8_10) { + Error("AliCFHeavyFlavourTaskMultiVar","Calculated bin lim for pt - 2nd range - differs from expected!\n"); + } + for(Int_t i=0; i<=nbin0_8_10; i++) binLim0[i+nbin0_0_4+nbin0_4_8]=(Double_t)ptmin_8_10 + (ptmax_8_10-ptmin_8_10)/nbin0_8_10*(Double_t)i ; + + // y + for(Int_t i=0; i<=nbin1; i++) binLim1[i]=(Double_t)ymin + (ymax-ymin) /nbin1*(Double_t)i ; + + // cosThetaStar + for(Int_t i=0; i<=nbin2; i++) binLim2[i]=(Double_t)cosmin + (cosmax-cosmin) /nbin2*(Double_t)i ; + + // ptPi + for(Int_t i=0; i<=nbin3_0_4; i++) binLim3[i]=(Double_t)ptmin_0_4 + (ptmax_0_4-ptmin_0_4)/nbin3_0_4*(Double_t)i ; + if (binLim3[nbin3_0_4] != ptmin_4_8) { + Error("AliCFHeavyFlavourTaskMultiVar","Calculated bin lim for ptPi - 1st range - differs from expected!"); + } + for(Int_t i=0; i<=nbin3_4_8; i++) binLim3[i+nbin3_0_4]=(Double_t)ptmin_4_8 + (ptmax_4_8-ptmin_4_8)/nbin3_4_8*(Double_t)i ; + if (binLim3[nbin3_0_4+nbin3_4_8] != ptmin_8_10) { + Error("AliCFHeavyFlavourTaskMultiVar","Calculated bin lim for ptPi - 2nd range - differs from expected!\n"); + } + for(Int_t i=0; i<=nbin3_8_10; i++) binLim3[i+nbin3_0_4+nbin3_4_8]=(Double_t)ptmin_8_10 + (ptmax_8_10-ptmin_8_10)/nbin3_8_10*(Double_t)i ; + + // ptKa + for(Int_t i=0; i<=nbin4_0_4; i++) binLim4[i]=(Double_t)ptmin_0_4 + (ptmax_0_4-ptmin_0_4)/nbin4_0_4*(Double_t)i ; + if (binLim4[nbin4_0_4] != ptmin_4_8) { + Error("AliCFHeavyFlavourTaskMultiVar","Calculated bin lim for ptKa - 1st range - differs from expected!"); + } + for(Int_t i=0; i<=nbin4_4_8; i++) binLim4[i+nbin4_0_4]=(Double_t)ptmin_4_8 + (ptmax_4_8-ptmin_4_8)/nbin4_4_8*(Double_t)i ; + if (binLim4[nbin4_0_4+nbin4_4_8] != ptmin_8_10) { + Error("AliCFHeavyFlavourTaskMultiVar","Calculated bin lim for ptKa - 2nd range - differs from expected!\n"); + } + for(Int_t i=0; i<=nbin4_8_10; i++) binLim4[i+nbin4_0_4+nbin4_4_8]=(Double_t)ptmin_8_10 + (ptmax_8_10-ptmin_8_10)/nbin4_8_10*(Double_t)i ; + + // cT + for(Int_t i=0; i<=nbin5; i++) binLim5[i]=(Double_t)cTmin + (cTmax-cTmin) /nbin5*(Double_t)i ; + + Info("AliCFHeavyFlavourTaskMultiVar","Printing lower limits for bins in pt"); + for (Int_t i =0; i<= iBin[0]; i++){ + Info("AliCFHeavyFlavourTaskMultiVar",Form("i-th bin, lower limit = %f", binLim0[i])); + } + Info("Printing lower limits for bins in ptPi"); + for (Int_t i =0; i<= iBin[3]; i++){ + Info("AliCFHeavyFlavourTaskMultiVar",Form("i-th bin, lower limit = %f", binLim3[i])); + } + Info("Printing lower limits for bins in ptKa"); + for (Int_t i =0; i<= iBin[4]; i++){ + Info("AliCFHeavyFlavourTaskMultiVar",Form("i-th bin, lower limit = %f", binLim4[i])); + } + + //one "container" for MC + AliCFContainer* container = new AliCFContainer("container","container for tracks",nstep,nvar,iBin); + //setting the bin limits + container -> SetBinLimits(ipt,binLim0); + container -> SetBinLimits(iy,binLim1); + container -> SetBinLimits(icosThetaStar,binLim2); + container -> SetBinLimits(ipTpi,binLim3); + container -> SetBinLimits(ipTk,binLim4); + container -> SetBinLimits(icT,binLim5); + + //CREATE THE CUTS ----------------------------------------------- + + // Gen-Level kinematic cuts + AliCFTrackKineCuts *mcKineCuts = new AliCFTrackKineCuts("mcKineCuts","MC-level kinematic cuts"); + // mcKineCuts->SetPtRange(ptmin,ptmax); + // mcKineCuts->SetRapidityRange(ymin,ymax); + // mcKineCuts->SetChargeMC(charge); + + //Particle-Level cuts: + AliCFParticleGenCuts* mcGenCuts = new AliCFParticleGenCuts("mcGenCuts","MC particle generation cuts"); + //mcGenCuts->SetRequireIsPrimary(); + mcGenCuts->SetRequirePdgCode(PDG, kTRUE); // kTRUE set in order to include D0_bar + mcGenCuts->SetAODMC(1); //special flag for reading MC in AOD tree (important) + + // Rec-Level kinematic cuts + AliCFTrackKineCuts *recKineCuts = new AliCFTrackKineCuts("recKineCuts","rec-level kine cuts"); + // recKineCuts->SetPtRange(ptmin,ptmax); + // recKineCuts->SetRapidityRange(ymin,ymax); + // recKineCuts->SetChargeRec(charge); + + AliCFTrackQualityCuts *recQualityCuts = new AliCFTrackQualityCuts("recQualityCuts","rec-level quality cuts"); + //recQualityCuts->SetStatus(AliESDtrack::kITSrefit); + + AliCFTrackIsPrimaryCuts *recIsPrimaryCuts = new AliCFTrackIsPrimaryCuts("recIsPrimaryCuts","rec-level isPrimary cuts"); + //recIsPrimaryCuts->SetAODType(AliAODTrack::kPrimary); + + printf("CREATE MC KINE CUTS\n"); + TObjArray* mcList = new TObjArray(0) ; + mcList->AddLast(mcKineCuts); + mcList->AddLast(mcGenCuts); + + printf("CREATE RECONSTRUCTION CUTS\n"); + TObjArray* recList = new TObjArray(0) ; + recList->AddLast(recKineCuts); + recList->AddLast(recQualityCuts); + recList->AddLast(recIsPrimaryCuts); + + //CREATE THE INTERFACE TO CORRECTION FRAMEWORK USED IN THE TASK + printf("CREATE INTERFACE AND CUTS\n"); + AliCFManager* man = new AliCFManager() ; + man->SetParticleContainer (container); + man->SetParticleCutsList(0 , mcList); // MC + man->SetParticleCutsList(1 , recList); // AOD + + //CREATE THE TASK + printf("CREATE TASK\n"); + // create the task + AliCFHeavyFlavourTaskMultiVar *task = new AliCFHeavyFlavourTaskMultiVar("AliCFHeavyFlavourTaskMultiVar"); + task->SetFillFromGenerated(kFALSE); + task->SetCFManager(man); //here is set the CF manager + + //SETUP THE ANALYSIS MANAGER TO READ INPUT CHAIN AND WRITE DESIRED OUTPUTS + printf("CREATE ANALYSIS MANAGER\n"); + // Make the analysis manager + AliAnalysisManager *mgr = new AliAnalysisManager("TestManager"); + + mgr->SetAnalysisType(AliAnalysisManager::kLocalAnalysis); + + AliAODInputHandler* dataHandler = new AliAODInputHandler(); + mgr->SetInputEventHandler(dataHandler); + + // Create and connect containers for input/output + + // ------ input data ------ + AliAnalysisDataContainer *cinput0 = mgr->GetCommonInputContainer(); + + // ----- output data ----- + + //slot 0 : default output tree (by default handled by AliAnalysisTaskSE) + AliAnalysisDataContainer *coutput0 = mgr->CreateContainer("ctree0", TTree::Class(),AliAnalysisManager::kOutputContainer,"output.root"); + + //now comes user's output objects : + + // output TH1I for event counting + AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("chist0", TH1I::Class(),AliAnalysisManager::kOutputContainer,"output.root"); + // output Correction Framework Container (for acceptance & efficiency calculations) + AliAnalysisDataContainer *coutput2 = mgr->CreateContainer("ccontainer0", AliCFContainer::Class(),AliAnalysisManager::kOutputContainer,"output.root"); + + //cinput0->SetData(analysisChain); + + mgr->AddTask(task); + + mgr->ConnectInput(task,0,mgr->GetCommonInputContainer()); + mgr->ConnectOutput(task,0,coutput0); + mgr->ConnectOutput(task,1,coutput1); + mgr->ConnectOutput(task,2,coutput2); + + /* + mgr->ConnectOutput(task,0,mgr->GetCommonOutputContainer()); + mgr->ConnectOutput(task,1,mgr->GetCommonOutputContainer()); + mgr->ConnectOutput(task,2,mgr->GetCommonOutputContainer()); + */ + + printf("READY TO RUN\n"); + //RUN !!! + if (mgr->InitAnalysis()) { + mgr->PrintStatus(); + mgr->StartAnalysis("local",analysisChain); + } + + benchmark.Stop("AliCFHeavyFlavourTaskMultiVar"); + benchmark.Show("AliCFHeavyFlavourTaskMultiVar"); + + return kTRUE ; +} + +void Load() { + + //load the required aliroot libraries + gSystem->Load("libANALYSIS") ; + gSystem->Load("libANALYSISalice") ; + gSystem->Load("libCORRFW.so") ; + gSystem->Load("libPWG3base.so"); + gSystem->Load("libPWG3vertexingHF.so"); + + gSystem->Load("libTree.so"); + gSystem->Load("libGeom.so"); + gSystem->Load("libPhysics.so"); + gSystem->Load("libVMC.so"); + gSystem->Load("libSTEERBase.so"); + gSystem->Load("libESD.so"); + gSystem->Load("libAOD.so"); + +} diff --git a/PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVar.cxx b/PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVar.cxx new file mode 100644 index 00000000000..a1d01d8d08d --- /dev/null +++ b/PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVar.cxx @@ -0,0 +1,806 @@ +/************************************************************************** + * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +//----------------------------------------------------------------------- +// Class for HF corrections as a function of many variables +// +// Example of task running on AliEn +// which provides standard way of calculating acceptance and efficiency +// between different steps of the procedure. +// The ouptut of the task is a AliCFContainer from which the efficiencies +// can be calculated +//----------------------------------------------------------------------- +// Author : C. Zampolli, CERN, on the basis of R. Vernet's example +//----------------------------------------------------------------------- + +#include +#include +#include +#include + +#include "AliCFHeavyFlavourTaskMultiVar.h" +#include "AliStack.h" +#include "AliMCEvent.h" +#include "AliCFManager.h" +#include "AliCFContainer.h" +#include "AliLog.h" +#include "AliAODEvent.h" +#include "AliAODRecoDecay.h" +#include "AliAODRecoDecayHF.h" +#include "AliAODRecoDecayHF2Prong.h" +#include "AliAODMCParticle.h" + +//__________________________________________________________________________ +AliCFHeavyFlavourTaskMultiVar::AliCFHeavyFlavourTaskMultiVar() : + AliAnalysisTaskSE(), + fPDG(0), + fCFManager(0x0), + fHistEventsProcessed(0x0), + fCountMC(0), + fEvents(0), + fFillFromGenerated(kFALSE), + fSkipped(0) +{ + // + //Default ctor + // +} +//___________________________________________________________________________ +AliCFHeavyFlavourTaskMultiVar::AliCFHeavyFlavourTaskMultiVar(const Char_t* name) : + AliAnalysisTaskSE(name), + fPDG(0), + fCFManager(0x0), + fHistEventsProcessed(0x0), + fCountMC(0), + fEvents(0), + fFillFromGenerated(kFALSE), + fSkipped(0) +{ + // + // Constructor. Initialization of Inputs and Outputs + // + Info("AliCFHeavyFlavourTaskMultiVar","Calling Constructor"); + /* + DefineInput(0) and DefineOutput(0) + are taken care of by AliAnalysisTaskSE constructor + */ + DefineOutput(1,TH1I::Class()); + DefineOutput(2,AliCFContainer::Class()); +} + +//___________________________________________________________________________ +AliCFHeavyFlavourTaskMultiVar& AliCFHeavyFlavourTaskMultiVar::operator=(const AliCFHeavyFlavourTaskMultiVar& c) +{ + // + // Assignment operator + // + if (this!=&c) { + AliAnalysisTaskSE::operator=(c) ; + fPDG = c.fPDG; + fCFManager = c.fCFManager; + fHistEventsProcessed = c.fHistEventsProcessed; + } + return *this; +} + +//___________________________________________________________________________ +AliCFHeavyFlavourTaskMultiVar::AliCFHeavyFlavourTaskMultiVar(const AliCFHeavyFlavourTaskMultiVar& c) : + AliAnalysisTaskSE(c), + fPDG(c.fPDG), + fCFManager(c.fCFManager), + fHistEventsProcessed(c.fHistEventsProcessed), + fCountMC(c.fCountMC), + fEvents(c.fEvents), + fFillFromGenerated(c.fFillFromGenerated), + fSkipped(c.fSkipped) +{ + // + // Copy Constructor + // +} + +//___________________________________________________________________________ +AliCFHeavyFlavourTaskMultiVar::~AliCFHeavyFlavourTaskMultiVar() { + // + //destructor + // + Info("~AliCFHeavyFlavourTaskMultiVar","Calling Destructor"); + if (fCFManager) delete fCFManager ; + if (fHistEventsProcessed) delete fHistEventsProcessed ; +} + +//_________________________________________________ +void AliCFHeavyFlavourTaskMultiVar::UserExec(Option_t *) +{ + // + // Main loop function + // + + if (!fInputEvent) { + Error("UserExec","NO EVENT FOUND!"); + return; + } + + fEvents++; + if (fEvents%10000 ==0) AliInfo(Form("Event %d",fEvents)); + AliAODEvent* aodEvent = dynamic_cast(fInputEvent); + fCFManager->SetEventInfo(aodEvent); + + // MC-event selection + Double_t containerInput[6] ; + + //loop on the MC event + + TClonesArray* mcArray = dynamic_cast(aodEvent->FindListObject(AliAODMCParticle::StdBranchName())); + if (!mcArray) AliError("Could not find Monte-Carlo in AOD"); + Int_t icountMC = 0; + + for (Int_t iPart=0; iPartGetEntriesFast(); iPart++) { + AliAODMCParticle* mcPart = dynamic_cast(mcArray->At(iPart)); + if (!mcPart) { + AliWarning("Particle not found in tree, skipping"); + continue; + } + + // check the MC-level cuts + if (!fCFManager->CheckParticleCuts(0, mcPart)) continue; // 0 stands for MC level + + // fill the container for Gen-level selection + Double_t vectorMC[6] = {9999.,9999.,9999.,9999.,9999.,9999.}; + if (GetGeneratedValuesFromMCParticle(mcPart, mcArray, vectorMC)){ + containerInput[0] = vectorMC[0]; + containerInput[1] = vectorMC[1] ; + containerInput[2] = vectorMC[2] ; + containerInput[3] = vectorMC[3] ; + containerInput[4] = vectorMC[4] ; + containerInput[5] = vectorMC[5] ; // in micron + fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated); + icountMC++; + } + else { + AliDebug(3,"Problems in filling the container"); + continue; + } + } + + AliDebug(2, Form("Found %i MC particles that are D0!!",icountMC)); + AliDebug(2, Form("aodEvent = %p ",aodEvent)); + + // Now go to rec level + fCountMC += icountMC; + // load heavy flavour vertices + + TClonesArray *arrayVerticesHF = (TClonesArray*)((aodEvent->GetList())->FindObject("D0toKpi")); + if (!arrayVerticesHF) AliError("Could not find array of HF vertices"); + AliDebug(2, Form("Found %d vertices",arrayVerticesHF->GetEntriesFast())); + + for (Int_t iVertex = 0; iVertexGetEntriesFast(); iVertex++) { + + AliAODRecoDecayHF2Prong* vtx = (AliAODRecoDecayHF2Prong*)arrayVerticesHF->At(iVertex); + + // cuts can't be applied to RecoDecays particles + // if (!fCFManager->CheckParticleCuts(1 , vtx)) continue; // 1 stands for AOD level + + // find associated MC particle + Int_t mcLabel = IsMcVtx(vtx) ; + if (mcLabel == -1) + { + AliDebug(2,"No MC particle found"); + } + else { + AliAODMCParticle* mcVtxHF = (AliAODMCParticle*)mcArray->At(mcLabel); + if (!mcVtxHF) { + AliWarning("Could not find associated MC in AOD MC tree"); + continue; + } + + // check if associated MC v0 passes the cuts + if (!fCFManager->CheckParticleCuts(0 ,mcVtxHF)) { // 0 stands for MC + AliDebug(2, "Skipping the particles due to cuts"); + continue; + } + + // fill the container + // either with reconstructed values.... + Double_t pt = vtx->Pt(); + Double_t rapidity = vtx->YD0(); + + Double_t cosThetaStar = 9999.; + Double_t pTpi = 0.; + Double_t pTK = 0.; + Int_t pdgCode = mcVtxHF->GetPdgCode(); + if (pdgCode > 0){ + cosThetaStar = vtx->CosThetaStarD0(); + pTpi = vtx->PtProng(0); + pTK = vtx->PtProng(1); + } + else { + cosThetaStar = vtx->CosThetaStarD0bar(); + pTpi = vtx->PtProng(1); + pTK = vtx->PtProng(0); + } + + Double_t cT = vtx->CtD0(); + AliDebug(2, Form("cT from reconstructed vertex = %f (micron)",cT*1E4)); + AliDebug(2, Form("pdg code from MC = %d",TMath::Abs(mcVtxHF->GetPdgCode()))); + + // ... or with generated values + + if (!fFillFromGenerated){ + containerInput[0] = pt; + containerInput[1] = rapidity; + containerInput[2] = cosThetaStar; + containerInput[3] = pTpi; + containerInput[4] = pTK; + containerInput[5] = cT*1.E4; // in micron + } + else { + Double_t vectorMC[6] = {9999.,9999.,9999.,9999.,9999.,9999.}; + if (GetGeneratedValuesFromMCParticle(mcVtxHF, mcArray, vectorMC)){ + containerInput[0] = vectorMC[0]; + containerInput[1] = vectorMC[1] ; + containerInput[2] = vectorMC[2] ; + containerInput[3] = vectorMC[3] ; + containerInput[4] = vectorMC[4] ; + containerInput[5] = vectorMC[5] ; // in micron + } + else { + AliDebug(3,"Problems in filling the container"); + continue; + } + } + AliDebug(2, Form("Filling the container with pt = %f, rapidity = %f, cosThetaStar = %f, pTpi = %f, pTK = %f, cT = %f", containerInput[0], containerInput[1], containerInput[2], containerInput[3], containerInput[4], containerInput[5])); + fCFManager->GetParticleContainer()->Fill(containerInput,1) ; + } + /* + // for debugging only, filling container with dummy values + Double_t pt = 1.5; + Double_t rapidity = 0.5; + Double_t cosThetaStar = 0.7; + Double_t pTpi = 2.5; + Double_t pTK = 1; + Double_t cT = 3; + containerInput[0] = pt; + containerInput[1] = rapidity; + containerInput[2] = cosThetaStar; + containerInput[3] = pTpi; + containerInput[4] = pTK; + containerInput[5] = cT; + fCFManager->GetParticleContainer()->Fill(containerInput,1) ; + */ + + } + + fHistEventsProcessed->Fill(0); + /* PostData(0) is taken care of by AliAnalysisTaskSE */ + PostData(1,fHistEventsProcessed) ; + PostData(2,fCFManager->GetParticleContainer()) ; +} + + +//___________________________________________________________________________ +void AliCFHeavyFlavourTaskMultiVar::Terminate(Option_t*) +{ + // The Terminate() function is the last function to be called during + // a query. It always runs on the client, it can be used to present + // the results graphically or save the results to file. + + Info("Terminate",""); + AliAnalysisTaskSE::Terminate(); + + AliInfo(Form("Found %i MC particles that are D0 in MC in %d events",fCountMC,fEvents)); + AliInfo(Form("Found %i reco D0 that are not decaying in K+pi in %d events",fSkipped,fEvents)); + + // draw some example plots.... + + AliCFContainer *cont= dynamic_cast (GetOutputData(2)); + + // projecting the containers to obtain histograms + // first argument = variable, second argument = step + TH1D* h00 = cont->ShowProjection(0,0) ; // pt + TH1D* h10 = cont->ShowProjection(1,0) ; // rapidity + TH1D* h20 = cont->ShowProjection(2,0) ; // cosThetaStar + TH1D* h30 = cont->ShowProjection(3,0) ; // pTpi + TH1D* h40 = cont->ShowProjection(4,0) ; // pTK + TH1D* h50 = cont->ShowProjection(5,0) ; // cT + + TH1D* h01 = cont->ShowProjection(0,1) ; // pt + TH1D* h11 = cont->ShowProjection(1,1) ; // rapidity + TH1D* h21 = cont->ShowProjection(2,1) ; // cosThetaStar + TH1D* h31 = cont->ShowProjection(3,1) ; // pTpi + TH1D* h41 = cont->ShowProjection(4,1) ; // pTK + TH1D* h51 = cont->ShowProjection(5,1) ; // cT + + h00->SetTitle("pT_D0 (GeV/c)"); + h10->SetTitle("rapidity"); + h20->SetTitle("cosThetaStar"); + h30->SetTitle("pT_pi (GeV/c)"); + h40->SetTitle("pT_K (Gev/c)"); + h50->SetTitle("cT (#mum)"); + + h00->GetXaxis()->SetTitle("pT_D0 (GeV/c)"); + h10->GetXaxis()->SetTitle("rapidity"); + h20->GetXaxis()->SetTitle("cosThetaStar"); + h30->GetXaxis()->SetTitle("pT_pi (GeV/c)"); + h40->GetXaxis()->SetTitle("pT_K (Gev/c)"); + h50->GetXaxis()->SetTitle("cT (#mum)"); + + h01->SetTitle("pT_D0 (GeV/c)"); + h11->SetTitle("rapidity"); + h21->SetTitle("cosThetaStar"); + h31->SetTitle("pT_pi (GeV/c)"); + h41->SetTitle("pT_K (Gev/c)"); + h51->SetTitle("cT (#mum)"); + + h01->GetXaxis()->SetTitle("pT_D0 (GeV/c)"); + h11->GetXaxis()->SetTitle("rapidity"); + h21->GetXaxis()->SetTitle("cosThetaStar"); + h31->GetXaxis()->SetTitle("pT_pi (GeV/c)"); + h41->GetXaxis()->SetTitle("pT_K (Gev/c)"); + h51->GetXaxis()->SetTitle("cT (#mum)"); + + Double_t max0 = h00->GetMaximum(); + Double_t max1 = h10->GetMaximum(); + Double_t max2 = h20->GetMaximum(); + Double_t max3 = h30->GetMaximum(); + Double_t max4 = h40->GetMaximum(); + Double_t max5 = h50->GetMaximum(); + + h00->GetYaxis()->SetRangeUser(0,max0*1.2); + h10->GetYaxis()->SetRangeUser(0,max1*1.2); + h20->GetYaxis()->SetRangeUser(0,max2*1.2); + h30->GetYaxis()->SetRangeUser(0,max3*1.2); + h40->GetYaxis()->SetRangeUser(0,max4*1.2); + h50->GetYaxis()->SetRangeUser(0,max5*1.2); + + h01->GetYaxis()->SetRangeUser(0,max0*1.2); + h11->GetYaxis()->SetRangeUser(0,max1*1.2); + h21->GetYaxis()->SetRangeUser(0,max2*1.2); + h31->GetYaxis()->SetRangeUser(0,max3*1.2); + h41->GetYaxis()->SetRangeUser(0,max4*1.2); + h51->GetYaxis()->SetRangeUser(0,max5*1.2); + + h00->SetMarkerStyle(20); + h10->SetMarkerStyle(24); + h20->SetMarkerStyle(21); + h30->SetMarkerStyle(25); + h40->SetMarkerStyle(27); + h50->SetMarkerStyle(28); + + h00->SetMarkerColor(2); + h10->SetMarkerColor(2); + h20->SetMarkerColor(2); + h30->SetMarkerColor(2); + h40->SetMarkerColor(2); + h50->SetMarkerColor(2); + + h01->SetMarkerStyle(20) ; + h11->SetMarkerStyle(24) ; + h21->SetMarkerStyle(21) ; + h31->SetMarkerStyle(25) ; + h41->SetMarkerStyle(27) ; + h51->SetMarkerStyle(28) ; + + h01->SetMarkerColor(4); + h11->SetMarkerColor(4); + h21->SetMarkerColor(4); + h31->SetMarkerColor(4); + h41->SetMarkerColor(4); + h51->SetMarkerColor(4); + + gStyle->SetCanvasColor(0); + gStyle->SetFrameFillColor(0); + gStyle->SetTitleFillColor(0); + gStyle->SetStatColor(0); + + // drawing in 2 separate canvas for a matter of clearity + TCanvas * c1 =new TCanvas("c1","pT, rapidiy, cosThetaStar",800,1600); + TCanvas * c2 =new TCanvas("c2","pTpi, pTK, cT",800,1600); + c1->Divide(2,3); + c2->Divide(2,3); + + c1->cd(1); + h00->Draw("p"); + c1->cd(2); + h01->Draw("p"); + c1->cd(3); + h10->Draw("p"); + c1->cd(4); + h11->Draw("p"); + c1->cd(5); + h20->Draw("p"); + c1->cd(6); + h21->Draw("p"); + c1->cd(); + + c2->cd(1); + h30->Draw("p"); + c2->cd(2); + h31->Draw("p"); + c2->cd(3); + h40->Draw("p"); + c2->cd(4); + h41->Draw("p"); + c2->cd(5); + h50->Draw("p"); + c2->cd(6); + h51->Draw("p"); + c2->cd(); + + c1->SaveAs("Variables/pT_rapidity_cosThetaStar.eps"); + c2->SaveAs("Variables/pTpi_pTK_cT.eps"); + +} + +//___________________________________________________________________________ +void AliCFHeavyFlavourTaskMultiVar::UserCreateOutputObjects() { + //HERE ONE CAN CREATE OUTPUT OBJECTS, IN PARTICULAR IF THE OBJECT PARAMETERS DON'T NEED + //TO BE SET BEFORE THE EXECUTION OF THE TASK + // + Info("UserCreateOutputObjects","CreateOutputObjects of task %s\n", GetName()); + + //slot #1 + OpenFile(1); + fHistEventsProcessed = new TH1I("fHistEventsProcessed","",1,0,1) ; +} + +//___________________________________________________________________________ +Int_t AliCFHeavyFlavourTaskMultiVar::GetVtxLabel(Int_t* labels) { + // + // returns the label of the vertex, given the labels of the 2 daughter tracks + // returns -1 if the V0 is fake + // + + Int_t labD0daugh0=-1; + Int_t labD0daugh1=-1; + Int_t labD0=-1; + Int_t labMother = -1; + Int_t pdgMother = -1; + + AliAODEvent* aodEvent = dynamic_cast(fInputEvent); + TClonesArray* mcArray = dynamic_cast(aodEvent->FindListObject(AliAODMCParticle::StdBranchName())); + if (!mcArray) AliError("Could not find Monte-Carlo in AOD"); + + AliAODMCParticle *part0 = (AliAODMCParticle*)mcArray->At(labels[0]); + if(!part0) { + printf("no MC particle\n"); + return -1; + } + + while(part0->GetMother()>=0) { + labMother=part0->GetMother(); + part0 = (AliAODMCParticle*)mcArray->At(labMother); + if(!part0) { + printf("no MC mother particle\n"); + } + pdgMother = TMath::Abs(part0->GetPdgCode()); + AliDebug(2,Form("pdgMother from particle 0 = %d",pdgMother)); + if(pdgMother==421) { + AliDebug(2,Form("found one particle (first daughter) coming from a D0")); + labD0daugh0=labMother; + } + } + + AliAODMCParticle *part1 = (AliAODMCParticle*)mcArray->At(labels[1]); + if(!part1) { + printf("no MC particle\n"); + return -1; + } + + while(part1->GetMother()>=0) { + labMother=part1->GetMother(); + part1 = (AliAODMCParticle*)mcArray->At(labMother); + if(!part1) { + printf("no MC mother particle\n"); + } + pdgMother = TMath::Abs(part1->GetPdgCode()); + AliDebug(2,Form("pdgMother from particle 1 = %d",pdgMother)); + if(pdgMother==421) { + AliDebug(2,Form("found one particle (second daughter) coming from a D0")); + labD0daugh1=labMother; + } + } + if ((labD0daugh0 >= 0) && (labD0daugh1 >= 0) && (labD0daugh0 == labD0daugh1)){ + //AliDebug(2, Form("found particles from the same D0 mother")); + labD0 = labD0daugh0; + // check that the D0 decays in 2 prongs + AliAODMCParticle *part = (AliAODMCParticle*)mcArray->At(labD0); + if (TMath::Abs(part->GetDaughter(1)-part->GetDaughter(0))==1) { + return labD0; + } + else { + fSkipped+=1; + } + } + + return -1; + +} + +//___________________________________________________________________________ +Int_t AliCFHeavyFlavourTaskMultiVar::IsMcVtx(AliAODRecoDecayHF2Prong* vtx) { + // + // check if the passed vertex is associated to a MC one, + // and returns the corresponding geant label. + // returns -1 if it is fake (i.e. label<0). + // + + + AliAODTrack *trk0 = (AliAODTrack*)vtx->GetDaughter(0); + AliAODTrack *trk1 = (AliAODTrack*)vtx->GetDaughter(1); + + if (!trk0 || ! trk1) { + AliDebug(2, "problems with the tracks while looking for mother"); + //if (!trk0) AliDebug(2, "track 0 gives problems"); + //if (!trk1) AliDebug(2, "track 1 gives problems"); + return -1; + } + + Int_t labels[2] = {-1,-1}; + labels[0] = trk0->GetLabel(); + labels[1] = trk1->GetLabel(); + + if (labels[0] < 0 || labels[1] < 0) { + AliDebug(2, "problems with the labels of the tracks while looking for mother"); + return -1; + } + + return GetVtxLabel(&labels[0]) ; +} +//___________________________________________________________________________ +Double_t AliCFHeavyFlavourTaskMultiVar::CosThetaStar(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const { + + // + // to calculate cos(ThetaStar) for generated particle + // using the K, since mcPartDaughter0 and mcPartDaughter1 always correspond to K and pi respectively + // (see where the function is called) + // + + Int_t pdgvtx = mcPart->GetPdgCode(); + /* if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar) + Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode()); + Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode()); + AliInfo(Form("D0, with pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1)); + AliDebug(2,"This is a D0"); + AliAODMCParticle* mcPartdummy = mcPartDaughter0; + mcPartDaughter0 = mcPartDaughter1; + mcPartDaughter1 = mcPartdummy; + } + else{ + AliInfo("D0bar"); + } + */ + Int_t pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode()); + Int_t pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode()); + if (pdgvtx > 0) { // setting as the first daughter always the kaon, to be used to calculate cos(ThetaStar) + AliDebug(2,"D0"); + } + else{ + AliDebug(2,"D0bar"); + } + if (pdgprong0 == 211){ + AliDebug(2,Form("pdgprong0 = %d, pdgprong1 = %d, switching...",pdgprong0,pdgprong1)); + AliAODMCParticle* mcPartdummy = mcPartDaughter0; + mcPartDaughter0 = mcPartDaughter1; + mcPartDaughter1 = mcPartdummy; + pdgprong0 = TMath::Abs(mcPartDaughter0->GetPdgCode()); + pdgprong1 = TMath::Abs(mcPartDaughter1->GetPdgCode()); + } + + AliDebug(2,Form("After checking, pdgprong0 = %d, pdgprong1 = %d",pdgprong0,pdgprong1)); + Double_t massvtx = TDatabasePDG::Instance()->GetParticle(TMath::Abs(pdgvtx))->Mass(); + Double_t massp[2]; + massp[0] = TDatabasePDG::Instance()->GetParticle(pdgprong0)->Mass(); + massp[1] = TDatabasePDG::Instance()->GetParticle(pdgprong1)->Mass(); + + Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)-4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx); + + Double_t px = mcPartDaughter0->Px()+mcPartDaughter1->Px(); + Double_t py = mcPartDaughter0->Py()+mcPartDaughter1->Py(); + Double_t pz = mcPartDaughter0->Pz()+mcPartDaughter1->Pz(); + Double_t p = TMath::Sqrt(px*px+py*py+pz*pz); + Double_t e = TMath::Sqrt(massvtx*massvtx+p*p); + Double_t beta = p/e; + Double_t gamma = e/massvtx; + TVector3 mom(mcPartDaughter0->Px(),mcPartDaughter0->Py(),mcPartDaughter0->Pz()); + TVector3 momTot(mcPartDaughter0->Px()+mcPartDaughter1->Px(),mcPartDaughter0->Py()+mcPartDaughter1->Py(),mcPartDaughter0->Pz()+mcPartDaughter1->Pz()); + + Double_t qlprong = mom.Dot(momTot)/momTot.Mag(); // analog to AliAODRecoDecay::QlProng(0) + + AliDebug(2,Form("pStar = %f, beta = %f, gamma = %f, qlprong = %f, massp[0] = %f", pStar, beta, gamma, qlprong, massp[0])); + Double_t cts = (qlprong/gamma-beta*TMath::Sqrt(pStar*pStar+massp[0]*massp[0]))/pStar; + AliDebug(2,Form("cts = %f", cts)); + return cts; +} +//___________________________________________________________________________ +Double_t AliCFHeavyFlavourTaskMultiVar::CT(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const { + + // + // to calculate cT for generated particle + // + + Double_t xmcPart[3] = {0,0,0}; + Double_t xdaughter0[3] = {0,0,0}; + Double_t xdaughter1[3] = {0,0,0}; + mcPart->XvYvZv(xmcPart); // cm + mcPartDaughter0->XvYvZv(xdaughter0); // cm + mcPartDaughter1->XvYvZv(xdaughter1); //cm + Double_t prodVtxD0 = TMath::Sqrt(xmcPart[0]*xmcPart[0]+ + xmcPart[1]*xmcPart[1]+ + xmcPart[2]*xmcPart[2]); + Double_t prodVtxDaughter0 = TMath::Sqrt(xdaughter0[0]*xdaughter0[0]+ + xdaughter0[1]*xdaughter0[1]+ + xdaughter0[2]*xdaughter0[2]); + Double_t prodVtxDaughter1 = TMath::Sqrt(xdaughter1[0]*xdaughter1[0]+ + xdaughter1[1]*xdaughter1[1]+ + xdaughter1[2]*xdaughter1[2]); + + AliDebug(2, Form("D0: x = %f, y = %f, z = %f, production vertex distance = %f (cm), %f (micron)", xmcPart[0], xmcPart[1], xmcPart[2], prodVtxD0, prodVtxD0*1.E4)); + AliDebug(2, Form("Daughter0: x = %f, y = %f, z = %f, production vertex distance = %f (cm) %f (micron)", xdaughter0[0], xdaughter0[1], xdaughter0[2], prodVtxDaughter0, prodVtxDaughter0*1E4)); + AliDebug(2, Form("Daughter1: x = %f, y = %f, z = %f, production vertex distance = %f (cm) %f (micron)", xdaughter1[0], xdaughter1[1], xdaughter1[2], prodVtxDaughter1, prodVtxDaughter1*1.E4)); + + Double_t cT0 = TMath::Sqrt((xdaughter0[0]-xmcPart[0])*(xdaughter0[0]-xmcPart[0])+ + (xdaughter0[1]-xmcPart[1])*(xdaughter0[1]-xmcPart[1])+ + (xdaughter0[2]-xmcPart[2])*(xdaughter0[2]-xmcPart[2])); + + Double_t cT1 = TMath::Sqrt((xdaughter1[0]-xmcPart[0])*(xdaughter1[0]-xmcPart[0])+ + (xdaughter1[1]-xmcPart[1])*(xdaughter1[1]-xmcPart[1])+ + (xdaughter1[2]-xmcPart[2])*(xdaughter1[2]-xmcPart[2])); + + if (cT0 != cT1) { + AliWarning("cT from daughter 0 different from cT from daughter 1! Using cT from daughter 0, but PLEASE, CHECK...."); + } + + // calculating cT from cT0 + + Double_t mass = TDatabasePDG::Instance()->GetParticle(mcPart->GetPdgCode())->Mass(); // mcPart->GetPdgCode() should return +/- 421 for the D0/D0bar + Double_t p = mcPart-> P(); + Double_t cT = cT0*mass/p; + AliDebug(2, Form("cT from daughter 0 = %f (micron)", cT0*1E4)); + AliDebug(2, Form("cT from daughter 1 = %f (micron)", cT1*1E4)); + AliDebug(2, Form("cT (with mass = %f and p = %f) = %f (micron)", mass, p, cT*1E4)); + return cT; +} +//________________________________________________________________________________ +Bool_t AliCFHeavyFlavourTaskMultiVar::GetGeneratedValuesFromMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC) const { + + // + // collecting all the necessary info (pt, y, cosThetaStar, ptPi, ptKa, cT) from MC particle + // + + Bool_t isOk = kFALSE; + + // check whether the D0 decays in pi+K + // to be added!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + // could use a cut in the CF, but not clear how to define a TDecayChannel + // implemented for the time being as a cut on the number of daughters - see later when + // getting the daughters + + // getting the daughters + Int_t daughter0 = mcPart->GetDaughter(0); + Int_t daughter1 = mcPart->GetDaughter(1); + AliDebug(2, Form("daughter0 = %d and daughter1 = %d",daughter0,daughter1)); + if (daughter0 == 0 || daughter1 == 0) { + AliDebug(2, "Error! the D0 MC doesn't have correct daughters!!"); + return isOk; + } + if (TMath::Abs(daughter1 - daughter0 != 1)) { + AliDebug(2, "The D0 MC doesn't come from a 2-prong decay, skipping!!"); + return isOk; + } + AliAODMCParticle* mcPartDaughter0 = dynamic_cast(mcArray->At(daughter0)); + AliAODMCParticle* mcPartDaughter1 = dynamic_cast(mcArray->At(daughter1)); + if (!mcPartDaughter0 || !mcPartDaughter1) { + AliWarning("At least one Daughter Particle not found in tree, skipping"); + return isOk; + } + + Double_t vtx1[3] = {0,0,0}; // primary vertex + Double_t vtx2daughter0[3] = {0,0,0}; // secondary vertex from daughter 0 + Double_t vtx2daughter1[3] = {0,0,0}; // secondary vertex from daughter 1 + mcPart->XvYvZv(vtx1); // cm + // getting vertex from daughters + mcPartDaughter0->XvYvZv(vtx2daughter0); // cm + mcPartDaughter1->XvYvZv(vtx2daughter1); //cm + if (vtx2daughter0[0] != vtx2daughter1[0] && vtx2daughter0[1] != vtx2daughter1[1] && vtx2daughter0[2] != vtx2daughter1[2]) { + AliError("Daughters have different secondary vertex, skipping the track"); + return isOk; + } + Int_t nprongs = 2; + Short_t charge = 0; + // always instantiate the AliAODRecoDecay with the positive daughter first, the negative second + AliAODMCParticle* positiveDaugh = mcPartDaughter0; + AliAODMCParticle* negativeDaugh = mcPartDaughter1; + if (mcPartDaughter0->GetPdgCode()<0 && mcPartDaughter1->GetPdgCode()>0){ + // inverting in case the positive daughter is the second one + positiveDaugh = mcPartDaughter1; + negativeDaugh = mcPartDaughter0; + } + // getting the momentum from the daughters + Double_t px[2] = {positiveDaugh->Px(), negativeDaugh->Px()}; + Double_t py[2] = {positiveDaugh->Py(), negativeDaugh->Py()}; + Double_t pz[2] = {positiveDaugh->Pz(), negativeDaugh->Pz()}; + + Double_t d0[2] = {0.,0.}; + + AliAODRecoDecayHF* decay = new AliAODRecoDecayHF(vtx1,vtx2daughter0,nprongs,charge,px,py,pz,d0); + + Double_t cosThetaStar = 0.; + Double_t cosThetaStarD0 = 0.; + Double_t cosThetaStarD0bar = 0.; + cosThetaStarD0 = decay->CosThetaStar(1,421,211,321); + cosThetaStarD0bar = decay->CosThetaStar(0,421,321,211); + if (mcPart->GetPdgCode() == 421){ // D0 + AliDebug(3, Form("D0, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode())); + cosThetaStar = cosThetaStarD0; + } + else if (mcPart->GetPdgCode() == -421){ // D0bar{ + AliDebug(3, Form("D0bar, with pdgprong0 = %d, pdgprong1 = %d",mcPartDaughter0->GetPdgCode(),mcPartDaughter1->GetPdgCode())); + cosThetaStar = cosThetaStarD0bar; + } + else{ + AliWarning("There are problems!! particle was expected to be either a D0 or a D0bar, check..."); + return vectorMC; + } + if (cosThetaStar < -1 || cosThetaStar > 1) { + AliWarning("Invalid value for cosine Theta star, returning"); + return isOk; + } + + // calculate cos(Theta*) according to the method implemented herein + + Double_t cts = 9999.; + cts = CosThetaStar(mcPart, mcPartDaughter0, mcPartDaughter1); + if (cts < -1 || cts > 1) { + AliWarning("Invalid value for cosine Theta star from AliCFHeavyFlavourTaskMultiVar method"); + } + if (TMath::Abs(cts - cosThetaStar)>0.001) { + AliError(Form("cosThetaStar automatically calculated different from that manually calculated!!! cosThetaStar = %f, cosThetaStar = %f", cosThetaStar,cts)); + } + + Double_t cT = decay->Ct(421); + + // calculate cT from the method implemented herein + Double_t cT1 = 0.; + cT1 = CT(mcPart, mcPartDaughter0, mcPartDaughter1); + + if (TMath::Abs(cT1 - cT)>0.001) { + AliError(Form("cT automatically calculated different from that manually calculated!!! cT = %f, cT1 = %f",cT,cT1)); + } + + // get the pT of the daughters + + Double_t pTpi = 0.; + Double_t pTK = 0.; + + if (TMath::Abs(mcPartDaughter0->GetPdgCode()) == 211) { + pTpi = mcPartDaughter0->Pt(); + pTK = mcPartDaughter1->Pt(); + } + else { + pTpi = mcPartDaughter1->Pt(); + pTK = mcPartDaughter0->Pt(); + } + + vectorMC[0] = mcPart->Pt(); + vectorMC[1] = mcPart->Y() ; + vectorMC[2] = cosThetaStar ; + vectorMC[3] = pTpi ; + vectorMC[4] = pTK ; + vectorMC[5] = cT*1.E4 ; // in micron + isOk = kTRUE; + return isOk; +} + diff --git a/PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVar.h b/PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVar.h new file mode 100644 index 00000000000..9a7449d8d71 --- /dev/null +++ b/PWG3/vertexingHF/AliCFHeavyFlavourTaskMultiVar.h @@ -0,0 +1,81 @@ +#ifndef ALICFHEAVYFLAVOURTASKMULTIVAR_H +#define ALICFHEAVYFLAVOURTASKMULTIVAR_H +/************************************************************************** + * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +//----------------------------------------------------------------------- +// Class for HF corrections as a function of many variables +// Author : C. Zampolli, CERN +//----------------------------------------------------------------------- + + +#include "AliAnalysisTaskSE.h" + +class TH1I; +class TParticle ; +class TFile ; +class TClonesArray ; +class AliCFManager; +class AliAODRecoDecay; +class AliAODRecoDecayHF2Prong; +class AliAODMCParticle; + +class AliCFHeavyFlavourTaskMultiVar : public AliAnalysisTaskSE { + public: + + enum { + kStepGenerated = 0, + kStepReconstructed = 1, + }; + + AliCFHeavyFlavourTaskMultiVar(); + AliCFHeavyFlavourTaskMultiVar(const Char_t* name); + AliCFHeavyFlavourTaskMultiVar& operator= (const AliCFHeavyFlavourTaskMultiVar& c); + AliCFHeavyFlavourTaskMultiVar(const AliCFHeavyFlavourTaskMultiVar& c); + virtual ~AliCFHeavyFlavourTaskMultiVar(); + + // ANALYSIS FRAMEWORK STUFF to loop on data and fill output objects + void UserCreateOutputObjects(); + void UserExec(Option_t *option); + void Terminate(Option_t *); + + // CORRECTION FRAMEWORK RELATED FUNCTIONS + void SetCFManager(AliCFManager* io) {fCFManager = io;} // global correction manager + AliCFManager * GetCFManager() {return fCFManager;} // get corr manager + + void SetPDG(Int_t code) {fPDG = code; } // defines the PDG code of searched HF + Int_t IsMcVtx(AliAODRecoDecayHF2Prong* vtx); // checks if the AliAODv0 can be associated, returns mother label + Int_t GetVtxLabel(Int_t* labels); // returns label of vertex given the daughter labels + Double_t CosThetaStar(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const; // returns cos(ThetaStar) of the D0 decay + Double_t CT(AliAODMCParticle* mcPart, AliAODMCParticle* mcPartDaughter0, AliAODMCParticle* mcPartDaughter1) const; // returns cT of the D0 decay + void SetFillFromGenerated(Bool_t flag) {fFillFromGenerated = flag;} + Bool_t GetFillFromGenerated() const {return fFillFromGenerated;} + Bool_t GetGeneratedValuesFromMCParticle(AliAODMCParticle* mcPart, TClonesArray* mcArray, Double_t* vectorMC) const; + + protected: + Int_t fPDG; // PDG code of searched V0's + AliCFManager *fCFManager ; // pointer to the CF manager + TH1I *fHistEventsProcessed; //! simple histo for monitoring the number of events processed + Int_t fCountMC; // MC particle found + Int_t fEvents; // n. of events + Bool_t fFillFromGenerated; // flag to indicate whether data container should be filled + // with generated values also for reconstructed particles + Int_t fSkipped; // n. of reco decays skipped after checking on the decay (MUST be D0->K+pi) + + + ClassDef(AliCFHeavyFlavourTaskMultiVar,0); // class for HF corrections as a function of many variables +}; + +#endif diff --git a/PWG3/vertexingHF/ReadCFHeavyFlavourOutput.C b/PWG3/vertexingHF/ReadCFHeavyFlavourOutput.C new file mode 100644 index 00000000000..1c8368ad4ba --- /dev/null +++ b/PWG3/vertexingHF/ReadCFHeavyFlavourOutput.C @@ -0,0 +1,766 @@ +#include + +extern TRandom *gRandom; +extern TBenchmark *gBenchmark; +extern TSystem *gSystem; + +void ReadCFHeavyFlavourOutput(){ + + // example macro for reading a container coming from a CF analysis + // the generated and reconstructed distributions are got, + // the efficiencies are calculated, + // the reco distributions are corrected, + // the results plotted + + gROOT->SetStyle("Plain"); + gStyle->SetPalette(1); + gStyle->SetOptStat(1110); + gStyle->SetPalette(1); + gStyle->SetCanvasColor(0); + gStyle->SetFrameFillColor(0); + gStyle->SetOptTitle(0); + + gSystem->SetIncludePath("-I. -I$ALICE_ROOT/include -I$ROOTSYS/include"); + gSystem->Load("libANALYSIS.so"); + gSystem->Load("libANALYSISalice.so"); + gSystem->Load("$ALICE_ROOT/CORRFW/libCORRFW.so") ; + + //Setting up the container grid... + + const Int_t nstep=2; //number of selection steps (just 2 in this ex) + + const Int_t nvar=6; //number of variables on the grid: pt, y, cosThetaStar, ptPi, ptK, cT + + // Flag the sel steps. In this example, we have two, may be any nstep + Int_t stepGen=0; + Int_t stepRec=1; + + // the sensitive variables, their indeces + UInt_t ipt =0; + UInt_t iy =1; + UInt_t icTS = 2; // cTS stands for cosThetaStar + UInt_t iptPi= 3; + UInt_t iptK = 4; + UInt_t icT = 5; + + // Read the container from file + TFile *file = new TFile("output.root"); + AliCFContainer *data = (AliCFContainer*) (file->Get("container")); + + // Make some 1 & 2-D projections.. + // MC level + TCanvas *cmc1 = new TCanvas("cmc1","The MC distributions",0,800,900,1200); + TCanvas *cmc2 = new TCanvas("cmc2","The MC distributions",0,800,900,1200); + TCanvas *cmcpt = new TCanvas("cmcpt","pt distribution from MC",50,50,550,550); + TCanvas *cmcy = new TCanvas("cmcy","y distribution from MC",50,50,550,550); + TCanvas *cmccTS = new TCanvas("cmcTS","cosThetaStar distribution from MC",50,50,550,550); // cTS stands for cosThetaStar + TCanvas *cmcptPi = new TCanvas("cmcptPi","pt_pi distribution from MC",50,50,550,550); + TCanvas *cmcptK = new TCanvas("cmcptK","pt_K distribution from MC",50,50,550,550); + TCanvas *cmccT = new TCanvas("cmccT","cT distribution from MC",50,50,550,550); + + // Reco-aod level + TCanvas *cpt = new TCanvas("cpt","pt distribution from reco aod",50,50,550,550); + TCanvas *cy = new TCanvas("cy","y distribution from reco aod",50,50,550,550); + TCanvas *ccTS = new TCanvas("cTS","cosThetaStar distribution from reco aod",50,50,550,550); + TCanvas *cptPi = new TCanvas("cptPi","pt_pi distribution from reco aod",50,50,550,550); + TCanvas *cptK = new TCanvas("cptK","pt_K distribution from reco aod",50,50,550,550); + TCanvas *ccT = new TCanvas("ccT","cT distribution from reco aod",50,50,550,550); + + // some 2D distributions in pt and y + TCanvas *cGen2d = new TCanvas("cGen2d","2D distribution from MC ",50,50,550,550); + TCanvas *cRec2d = new TCanvas("cRec2d","2D distribution from reco aod",50,50,550,550); + + // 2D plots + TH2D *hMCpty2d = data->ShowProjection(ipt, iy, stepGen); + hMCpty2d->SetLineColor(2); + hMCpty2d->SetLineWidth(3); + hMCpty2d->SetMarkerColor(2); + hMCpty2d->SetMarkerStyle(20); + hMCpty2d->GetXaxis()->SetTitleOffset(1.2); + hMCpty2d->GetYaxis()->SetTitleOffset(1.5); + hMCpty2d->GetXaxis()->SetTitle("p_{T} (GeV/c), MC data"); + hMCpty2d->GetYaxis()->SetTitle("y, MC data"); + hMCpty2d->Draw("text"); + cGen2d->cd(); + cGen2d->SetLeftMargin(0.15); + cGen2d->SetRightMargin(0.05); + cGen2d->Update(); + + TH2D *hRECpty2d = data->ShowProjection(ipt, iy, stepRec); + hRECpty2d->Sumw2(); + hRECpty2d->SetLineColor(4); + hRECpty2d->SetLineWidth(3); + hRECpty2d->SetMarkerColor(4); + hRECpty2d->SetMarkerStyle(20); + hRECpty2d->GetXaxis()->SetTitleOffset(1.2); + hRECpty2d->GetYaxis()->SetTitleOffset(1.5); + hRECpty2d->GetXaxis()->SetTitle("p_{T} (GeV/c), AOD"); + hRECpty2d->GetYaxis()->SetTitle("y, AOD"); + hRECpty2d->Draw("text"); + cRec2d->cd(); + cRec2d->SetLeftMargin(0.15); + cRec2d->SetRightMargin(0.05); + cRec2d->Update(); + + // MC + REC 1D plots + // pt, y, cosThetaStar + cmc1->Divide(2,3); + + cmc1->cd(1); + TH1D *hMCpt1D = data->ShowProjection(ipt, stepGen); + Double_t maxpt = hMCpt1D->GetMaximum(); + hMCpt1D->GetYaxis()->SetRangeUser(0,maxpt*1.2); + hMCpt1D->Sumw2(); + hMCpt1D->SetMinimum(0.01); + hMCpt1D->SetLineColor(2); + hMCpt1D->SetLineWidth(3); + hMCpt1D->SetMarkerColor(2); + hMCpt1D->SetFillColor(2); + hMCpt1D->SetFillStyle(3005); + hMCpt1D->SetMarkerStyle(20); + hMCpt1D->GetXaxis()->SetTitleOffset(1.2); + hMCpt1D->GetXaxis()->SetTitle("p_{T} (GeV/c), MC data"); + hMCpt1D->Draw("hist"); + cmcpt->cd(); + cmcpt->SetLeftMargin(0.15); + cmcpt->SetRightMargin(0.05); + hMCpt1D->Draw("hist"); + hMCpt1D->Draw("err same"); + cmcpt->Update(); + + cmc1->cd(2); + TH1D *hRECpt1D = data->ShowProjection(ipt, stepRec); + hRECpt1D->GetYaxis()->SetRangeUser(0,maxpt*1.2); + hRECpt1D->SetLineColor(4); + hRECpt1D->SetLineWidth(3); + hRECpt1D->SetMarkerColor(4); + hRECpt1D->SetFillColor(4); + hRECpt1D->SetFillStyle(3004); + hRECpt1D->SetMarkerStyle(20); + hRECpt1D->GetXaxis()->SetTitleOffset(1.2); + hRECpt1D->Sumw2(); + hRECpt1D->SetMinimum(0.01); + hRECpt1D->GetXaxis()->SetTitle("p_{T} (GeV/c), AOD"); + hRECpt1D->Draw("hist"); + cpt->cd(); + cpt->SetLeftMargin(0.15); + cpt->SetRightMargin(0.05); + hRECpt1D->Draw("hist"); + hRECpt1D->Draw("err same"); + cpt->Update(); + + cmc1->cd(3); + TH1D *hMCy1D = data->ShowProjection(iy, stepGen); + Double_t maxy = hMCy1D->GetMaximum(); + hMCy1D->GetYaxis()->SetRangeUser(0,maxy*1.2); + hMCy1D->SetLineColor(2); + hMCy1D->SetLineWidth(3); + hMCy1D->SetMarkerColor(2); + hMCy1D->SetFillColor(2); + hMCy1D->SetFillStyle(3005); + hMCy1D->SetMarkerStyle(20); + hMCy1D->GetXaxis()->SetTitleOffset(1.2); + hMCy1D->Sumw2(); + hMCy1D->SetMinimum(0.01); + hMCy1D->GetXaxis()->SetTitle("y, MC data"); + hMCy1D->Draw("hist"); + cmcy->cd(); + cmcy->SetLeftMargin(0.15); + cmcy->SetRightMargin(0.05); + hMCy1D->Draw("hist"); + hMCy1D->Draw("err same"); + cmcy->Update(); + + cmc1->cd(4); + TH1D *hRECy1D = data->ShowProjection(iy, stepRec); + hRECy1D->GetYaxis()->SetRangeUser(0,maxy*1.2); + hRECy1D->SetLineColor(4); + hRECy1D->SetLineWidth(3); + hRECy1D->SetMarkerColor(4); + hRECy1D->SetFillColor(4); + hRECy1D->SetFillStyle(3004); + hRECy1D->SetMarkerStyle(20); + hRECy1D->GetXaxis()->SetTitleOffset(1.2); + hRECy1D->Sumw2(); + hRECy1D->SetMinimum(0.01); + hRECy1D->Draw("hist"); + cy->cd(); + cy->SetLeftMargin(0.15); + cy->SetRightMargin(0.05); + hRECy1D->GetXaxis()->SetTitle("y, AOD"); + hRECy1D->Draw("hist"); + hRECy1D->Draw("err same"); + cy->Update(); + + cmc1->cd(5); + TH1D *hMCcTS1D = data->ShowProjection(icTS, stepGen); + Double_t maxcTS = hMCcTS1D->GetMaximum(); + hMCcTS1D->GetYaxis()->SetRangeUser(0,maxcTS*1.2); + hMCcTS1D->SetLineColor(2); + hMCcTS1D->SetLineWidth(3); + hMCcTS1D->SetMarkerColor(2); + hMCcTS1D->SetFillColor(2); + hMCcTS1D->SetFillStyle(3005); + hMCcTS1D->SetMarkerStyle(20); + hMCcTS1D->GetXaxis()->SetTitleOffset(1.2); + hMCcTS1D->Sumw2(); + hMCcTS1D->SetMinimum(0.01); + hMCcTS1D->GetXaxis()->SetTitle("cosThetaStar, MC data"); + hMCcTS1D->Draw("hist"); + cmccTS->cd(); + cmccTS->SetLeftMargin(0.15); + cmccTS->SetRightMargin(0.05); + hMCcTS1D->Draw("hist"); + hMCcTS1D->Draw("err same"); + cmccTS->Update(); + + cmc1->cd(6); + TH1D *hRECcTS1D = data->ShowProjection(icTS, stepRec); + hRECcTS1D->GetYaxis()->SetRangeUser(0,maxcTS*1.2); + hRECcTS1D->SetLineColor(4); + hRECcTS1D->SetLineWidth(3); + hRECcTS1D->SetMarkerColor(4); + hRECcTS1D->SetFillColor(4); + hRECcTS1D->SetFillStyle(3004); + hRECcTS1D->SetMarkerStyle(20); + hRECcTS1D->GetXaxis()->SetTitleOffset(1.2); + hRECcTS1D->Sumw2(); + hRECcTS1D->SetMinimum(0.01); + hRECcTS1D->Draw("hist"); + ccTS->cd(); + ccTS->SetLeftMargin(0.15); + ccTS->SetRightMargin(0.05); + hRECcTS1D->GetXaxis()->SetTitle("cosThetaStar, AOD"); + hRECcTS1D->Draw("hist"); + hRECcTS1D->Draw("err same"); + ccTS->Update(); + + // ptPi, ptK, cT + cmc2->Divide(2,3); + + cmc2->cd(1); + TH1D *hMCptPi1D = data->ShowProjection(iptPi, stepGen); + Double_t maxptPi = hMCptPi1D->GetMaximum(); + hMCptPi1D->GetYaxis()->SetRangeUser(0,maxptPi*1.2); + hMCptPi1D->Sumw2(); + hMCptPi1D->SetMinimum(0.01); + hMCptPi1D->SetLineColor(2); + hMCptPi1D->SetLineWidth(3); + hMCptPi1D->SetMarkerColor(2); + hMCptPi1D->SetFillColor(2); + hMCptPi1D->SetFillStyle(3005); + hMCptPi1D->SetMarkerStyle(20); + hMCptPi1D->GetXaxis()->SetTitleOffset(1.2); + hMCptPi1D->GetXaxis()->SetTitle("p_{T, #pi} (GeV/c), MC data"); + hMCptPi1D->Draw("hist"); + cmcptPi->cd(); + cmcptPi->SetLeftMargin(0.15); + cmcptPi->SetRightMargin(0.05); + hMCptPi1D->Draw("hist"); + hMCptPi1D->Draw("err same"); + cmcptPi->Update(); + + cmc2->cd(2); + TH1D *hRECptPi1D = data->ShowProjection(iptPi, stepRec); + hRECptPi1D->GetYaxis()->SetRangeUser(0,maxptPi*1.2); + hRECptPi1D->SetLineColor(4); + hRECptPi1D->SetLineWidth(3); + hRECptPi1D->SetMarkerColor(4); + hRECptPi1D->SetFillColor(4); + hRECptPi1D->SetFillStyle(3004); + hRECptPi1D->SetMarkerStyle(20); + hRECptPi1D->GetXaxis()->SetTitleOffset(1.2); + hRECptPi1D->Sumw2(); + hRECptPi1D->SetMinimum(0.01); + hRECptPi1D->GetXaxis()->SetTitle("p_{T, #pi} (GeV/c), AOD"); + hRECptPi1D->Draw("hist"); + cptPi->cd(); + cptPi->SetLeftMargin(0.15); + cptPi->SetRightMargin(0.05); + hRECptPi1D->Draw("hist"); + hRECptPi1D->Draw("err same"); + cptPi->Update(); + + cmc2->cd(3); + TH1D *hMCptK1D = data->ShowProjection(iptK, stepGen); + Double_t maxptK = hMCptK1D->GetMaximum(); + hMCptK1D->GetYaxis()->SetRangeUser(0,maxptK*1.2); + hMCptK1D->SetLineColor(2); + hMCptK1D->SetLineWidth(3); + hMCptK1D->SetMarkerColor(2); + hMCptK1D->SetFillColor(2); + hMCptK1D->SetFillStyle(3005); + hMCptK1D->SetMarkerStyle(20); + hMCptK1D->GetXaxis()->SetTitleOffset(1.2); + hMCptK1D->Sumw2(); + hMCptK1D->SetMinimum(0.01); + hMCptK1D->GetXaxis()->SetTitle("p_{T, K} (GeV/c), MC data"); + hMCptK1D->Draw("hist"); + cmcptK->cd(); + cmcptK->SetLeftMargin(0.15); + cmcptK->SetRightMargin(0.05); + hMCptK1D->Draw("hist"); + hMCptK1D->Draw("err same"); + cmcptK->Update(); + + cmc2->cd(4); + TH1D *hRECptK1D = data->ShowProjection(iptK, stepRec); + hRECptK1D->GetYaxis()->SetRangeUser(0,maxptK*1.2); + hRECptK1D->SetLineColor(4); + hRECptK1D->SetLineWidth(3); + hRECptK1D->SetMarkerColor(4); + hRECptK1D->SetFillColor(4); + hRECptK1D->SetFillStyle(3004); + hRECptK1D->SetMarkerStyle(20); + hRECptK1D->GetXaxis()->SetTitleOffset(1.2); + hRECptK1D->Sumw2(); + hRECptK1D->SetMinimum(0.01); + hRECptK1D->Draw("hist"); + cptK->cd(); + cptK->SetLeftMargin(0.15); + cptK->SetRightMargin(0.05); + hRECptK1D->GetXaxis()->SetTitle("p_{T, K} (GeV/c), AOD"); + hRECptK1D->Draw("hist"); + hRECptK1D->Draw("err same"); + cptK->Update(); + + cmc2->cd(5); + TH1D *hMCcT1D = data->ShowProjection(icT, stepGen); + Double_t maxcT = hMCcT1D->GetMaximum(); + hMCcT1D->GetYaxis()->SetRangeUser(0,maxcT*1.2); + hMCcT1D->SetLineColor(2); + hMCcT1D->SetLineWidth(3); + hMCcT1D->SetMarkerColor(2); + hMCcT1D->SetFillColor(2); + hMCcT1D->SetFillStyle(3005); + hMCcT1D->SetMarkerStyle(20); + hMCcT1D->GetXaxis()->SetTitleOffset(1.2); + hMCcT1D->Sumw2(); + hMCcT1D->SetMinimum(0.01); + hMCcT1D->GetXaxis()->SetTitle("ct (#mum), MC data"); + hMCcT1D->Draw("hist"); + cmccT->cd(); + cmccT->SetLeftMargin(0.15); + cmccT->SetRightMargin(0.05); + hMCcT1D->Draw("hist"); + hMCcT1D->Draw("err same"); + cmccT->Update(); + + cmc2->cd(6); + TH1D *hRECcT1D = data->ShowProjection(icT, stepRec); + hRECcT1D->GetYaxis()->SetRangeUser(0,maxcT*1.2); + hRECcT1D->SetLineColor(4); + hRECcT1D->SetLineWidth(3); + hRECcT1D->SetMarkerColor(4); + hRECcT1D->SetFillColor(4); + hRECcT1D->SetFillStyle(3004); + hRECcT1D->SetMarkerStyle(20); + hRECcT1D->GetXaxis()->SetTitleOffset(1.2); + hRECcT1D->Sumw2(); + hRECcT1D->SetMinimum(0.01); + hRECcT1D->Draw("hist"); + ccT->cd(); + ccT->SetLeftMargin(0.15); + ccT->SetRightMargin(0.05); + hRECcT1D->GetXaxis()->SetTitle("c#t (#mum), AOD"); + hRECcT1D->Draw("hist"); + hRECcT1D->Draw("err same"); + ccT->Update(); + + /* + // printing on eps files + cmc1->Print("Plots/dataMC_pt_y_cTS.gif"); + cmc2->Print("Plots/dataMC_ptPi_ptK_cT.gif"); + cmcpt->Print("Plots/pt_Gen.eps"); + cmcy->Print("Plots/y_Gen.eps"); + cmccTS->Print("Plots/cTS_Gen.eps"); + cmcptPi->Print("Plots/ptPi_Gen.eps"); + cmcptK->Print("Plots/ptK_Gen.eps"); + cmccT->Print("Plots/cT_Gen.eps"); + cpt->Print("Plots/pt_Rec.eps"); + cy->Print("Plots/y_Rec.eps"); + ccTS->Print("Plots/cTS_Rec.eps"); + cptPi->Print("Plots/ptPi_Rec.eps"); + cptK->Print("Plots/ptK_Rec.eps"); + ccT->Print("Plots/cT_Rec.eps"); + cGen2d->Print("Plots/pt_y_Gen_2D.eps"); + cRec2d->Print("Plots/pt_y_Rec_2D.eps"); + + // printing on gif files + cmc1->Print("Plots/dataMC_pt_y_cTS.gif"); + cmc2->Print("Plots/dataMC_ptPi_ptK_cT.gif"); + cmcpt->Print("Plots/pt_Gen.eps"); + cmcy->Print("Plots/y_Gen.eps"); + cmccTS->Print("Plots/cTS_Gen.eps"); + cmcptPi->Print("Plots/ptPi_Gen.eps"); + cmcptK->Print("Plots/ptK_Gen.eps"); + cmccT->Print("Plots/cT_Gen.eps"); + cpt->Print("Plots/pt_Rec.eps"); + cy->Print("Plots/y_Rec.eps"); + ccTS->Print("Plots/cTS_Rec.eps"); + cptPi->Print("Plots/ptPi_Rec.eps"); + cptK->Print("Plots/ptK_Rec.eps"); + ccT->Print("Plots/cT_Rec.eps"); + cGen2d->Print("Plots/pt_y_Gen_2D.eps"); + cRec2d->Print("Plots/pt_y_Rec_2D.eps"); + */ + + //construct the efficiency grid from the data container + AliCFEffGrid *eff = new AliCFEffGrid("eff"," The efficiency",*data); + eff->CalculateEfficiency(stepRec,stepGen); //eff= step1/step0 + + //The efficiency along the variables, and some 2-D projections + TCanvas *ceff =new TCanvas("ceff"," Efficiency",0,0,1600,1200); + ceff->Divide(3,2); + TCanvas *ceff2D =new TCanvas("ceff2D"," Efficiency for pt and y",50,50,550,550); + TCanvas *ceffpt = new TCanvas("ceffpt","Efficiency vs pt",50,50,550,550); + TCanvas *ceffy = new TCanvas("ceffy","Efficiency vs y",50,50,550,550); + TCanvas *ceffcTS = new TCanvas("ceffcTS","Efficiency vs cosThetaStar",50,50,550,550); + TCanvas *ceffptPi = new TCanvas("ceffptPi","Efficiency vs ptPi",50,50,550,550); + TCanvas *ceffptK = new TCanvas("ceffptK","Efficiency vs ptK",50,50,550,550); + TCanvas *ceffcT = new TCanvas("ceffcT","Efficiency vs cT",50,50,550,550); + TCanvas *ceff2Dtext = new TCanvas("ceff2Dtext","Text plot for efficiency in pt and y",50,50,550,550); + + ceff->cd(1); + TH1D *hpteffCF = eff->Project(ipt); //the efficiency vs pt + hpteffCF->Sumw2(); + //hpteffCF->SetMinimum(0.01); + hpteffCF->SetLineColor(8); + hpteffCF->SetLineWidth(3); + hpteffCF->SetMarkerColor(8); + hpteffCF->SetMarkerStyle(20); + hpteffCF->GetXaxis()->SetTitleOffset(1.2); + hpteffCF->GetYaxis()->SetTitleOffset(1.5); + hpteffCF->Draw("hist"); + ceffpt->cd(); + ceffpt->SetLeftMargin(0.15); + ceffpt->SetRightMargin(0.05); + hpteffCF->GetXaxis()->SetTitle("p_{T} (GeV/c)"); + hpteffCF->GetYaxis()->SetTitle("Efficiency"); + hpteffCF->Draw("hist"); + hpteffCF->Draw("err same"); + + ceff->cd(2); + TH1D *hyeffCF = eff->Project(iy); //the efficiency vs y + hyeffCF->Sumw2(); + //hyeffCF->SetMinimum(0.01); + hyeffCF->SetLineColor(8); + hyeffCF->SetLineWidth(3); + hyeffCF->SetMarkerColor(8); + hyeffCF->SetMarkerStyle(20); + hyeffCF->GetXaxis()->SetTitleOffset(1.2); + hyeffCF->GetYaxis()->SetTitleOffset(1.5); + hyeffCF->Draw("hist"); + ceffy->cd(); + ceffy->SetLeftMargin(0.15); + ceffy->SetRightMargin(0.05); + hyeffCF->GetXaxis()->SetTitle("y"); + hyeffCF->GetYaxis()->SetTitle("Efficiency"); + hyeffCF->Draw("hist"); + hyeffCF->Draw("err same"); + + ceff->cd(3); + TH1D *hcTSeffCF = eff->Project(icTS); //the efficiency vs cosThetaStar + hcTSeffCF->Sumw2(); + //hcTSeffCF->SetMinimum(0.01); + hcTSeffCF->SetLineColor(8); + hcTSeffCF->SetLineWidth(3); + hcTSeffCF->SetMarkerColor(8); + hcTSeffCF->SetMarkerStyle(20); + hcTSeffCF->GetXaxis()->SetTitleOffset(1.2); + hcTSeffCF->GetYaxis()->SetTitleOffset(1.5); + hcTSeffCF->Draw("hist"); + ceffcTS->cd(); + ceffcTS->SetLeftMargin(0.15); + ceffcTS->SetRightMargin(0.05); + hcTSeffCF->GetXaxis()->SetTitle("cosThetaStar"); + hcTSeffCF->GetYaxis()->SetTitle("Efficiency"); + hcTSeffCF->Draw("hist"); + hcTSeffCF->Draw("err same"); + + ceff->cd(4); + TH1D *hptPieffCF = eff->Project(iptPi); //the efficiency vs ptPi + hptPieffCF->Sumw2(); + //hptPieffCF->SetMinimum(0.01); + hptPieffCF->SetLineColor(8); + hptPieffCF->SetLineWidth(3); + hptPieffCF->SetMarkerColor(8); + hptPieffCF->SetMarkerStyle(20); + hptPieffCF->GetXaxis()->SetTitleOffset(1.2); + hptPieffCF->GetYaxis()->SetTitleOffset(1.5); + hptPieffCF->Draw("hist"); + ceffptPi->cd(); + ceffptPi->SetLeftMargin(0.15); + ceffptPi->SetRightMargin(0.05); + hptPieffCF->GetXaxis()->SetTitle("p_{T, #pi} (GeV/c)"); + hptPieffCF->GetYaxis()->SetTitle("Efficiency"); + hptPieffCF->Draw("hist"); + hptPieffCF->Draw("err same"); + + ceff->cd(5); + TH1D *hptKeffCF = eff->Project(iptK); //the efficiency vs ptK + hptKeffCF->Sumw2(); + //hptKeffCF->SetMinimum(0.01); + hptKeffCF->SetLineColor(8); + hptKeffCF->SetLineWidth(3); + hptKeffCF->SetMarkerColor(8); + hptKeffCF->SetMarkerStyle(20); + hptKeffCF->GetXaxis()->SetTitleOffset(1.2); + hptKeffCF->GetYaxis()->SetTitleOffset(1.5); + hptKeffCF->Draw("hist"); + ceffptK->cd(); + ceffptK->SetLeftMargin(0.15); + ceffptK->SetRightMargin(0.05); + hptKeffCF->GetXaxis()->SetTitle("p_{T, K} (GeV/c)"); + hptKeffCF->GetYaxis()->SetTitle("Efficiency"); + hptKeffCF->Draw("hist"); + hptKeffCF->Draw("err same"); + + ceff->cd(6); + TH1D *hcTeffCF = eff->Project(icT); //the efficiency vs cT + hcTeffCF->Sumw2(); + //hcTeffCF->SetMinimum(0.01); + hcTeffCF->SetLineColor(8); + hcTeffCF->SetLineWidth(3); + hcTeffCF->SetMarkerColor(8); + hcTeffCF->SetMarkerStyle(20); + hcTeffCF->GetXaxis()->SetTitleOffset(1.2); + hcTeffCF->GetYaxis()->SetTitleOffset(1.5); + hcTeffCF->Draw("hist"); + ceffcT->cd(); + ceffcT->SetLeftMargin(0.15); + ceffcT->SetRightMargin(0.05); + hcTeffCF->GetXaxis()->SetTitle("c#t (#mum)"); + hcTeffCF->GetYaxis()->SetTitle("Efficiency"); + hcTeffCF->Draw("hist"); + hcTeffCF->Draw("err same"); + + ceff2D->cd(); + TH2D *hptyeffCF = eff->Project(ipt,iy); //look at the numerator + //hptyeffCF->SetMinimum(0.01); + hptyeffCF->SetMarkerColor(8); + hptyeffCF->SetLineColor(8); + hptyeffCF->SetMinimum(0.01); + hptyeffCF->Draw("lego"); + ceff2Dtext->cd(); + hptyeffCF->GetXaxis()->SetTitle("p_{T} (GeV/c)"); + hptyeffCF->GetXaxis()->SetTitleOffset(1.2); + hptyeffCF->GetYaxis()->SetTitle("y"); + hptyeffCF->GetYaxis()->SetTitleOffset(1.2); + hptyeffCF->Draw("text"); + + /* + // printing eps files + ceff->Print("Plots/efficiencies.eps"); + ceffpt->Print("Plots/effpt.eps"); + ceffy->Print("Plots/effy.eps"); + ceffcTS->Print("Plots/effcTS.eps"); + ceffptPi->Print("Plots/effptPi.eps"); + ceffptK->Print("Plots/effptK.eps"); + ceffcT->Print("Plots/effcT.eps"); + ceff2D->Print("Plots/eff2D_pt_y.eps"); + ceff2Dtext->Print("Plots/eff2d_text_pt_y.eps"); + + // printing gif files + ceff->Print("Plots/efficiencies.gif"); + ceffpt->Print("Plots/effpt.gif"); + ceffy->Print("Plots/effy.gif"); + ceffcTS->Print("Plots/effcTS.gif"); + ceffptPi->Print("Plots/effptPi.gif"); + ceffptK->Print("Plots/effptK.gif"); + ceffcT->Print("Plots/effcT.gif"); + ceff2D->Print("Plots/eff2D_pt_y.gif"); + ceff2Dtext->Print("Plots/eff2d_text_pt_y.gif"); + */ + + // applying efficiencies - using projections, not the ApplyEffCorrection method + + TCanvas *cmultpiplypt = new TCanvas("cmultiplypt","Reco From Eff in pt distribution",50,50,550,550); + TCanvas *cmultpiplyy = new TCanvas("cmultiplyy","Reco From Eff in y distribution",50,50,550,550); + TCanvas *cmultpiplycTS = new TCanvas("cmultiplycTS","Reco From Eff in cTS distribution",50,50,550,550); + TCanvas *cmultpiplyptPi = new TCanvas("cmultiplyptPi","Reco From Eff in ptPi distribution",50,50,550,550); + TCanvas *cmultpiplyptK = new TCanvas("cmultiplyptK","Reco From Eff in ptK distribution",50,50,550,550); + TCanvas *cmultpiplycT = new TCanvas("cmultiplycT","Reco From Eff in cT distribution",50,50,550,550); + + TH1D *hmultiplypt = new TH1D("hmultiplypt","hmultiplypt",13,0,10); + cout << " bin for histo MC = " << hMCpt1D->GetNbinsX() << " while for RECO histo = " << hRECpt1D->GetNbinsX() << " while for efficiency histo = " << hpteffCF->GetNbinsX() << endl; + + const Double_t ptmin_0_4 = 0.0 ; + const Double_t ptmax_0_4 = 4.0 ; + const Double_t ptmin_4_8 = 4.0 ; + const Double_t ptmax_4_8 = 8.0 ; + const Double_t ptmin_8_10 = 8.0 ; + const Double_t ptmax_8_10 = 10.0 ; + const Int_t nbin0_0_4 = 8 ; //bins in pt + const Int_t nbin0_4_8 = 4 ; //bins in pt + const Int_t nbin0_8_10 = 1 ; //bins in pt + const Int_t nbins = nbin0_0_4 + nbin0_4_8 + nbin0_8_10; + Double_t binLim0[nbins+1]; + for(Int_t i=0; i<=nbin0_0_4; i++) binLim0[i]=(Double_t)ptmin_0_4 + (ptmax_0_4-ptmin_0_4)/nbin0_0_4*(Double_t)i ; + if (binLim0[nbin0_0_4] != ptmin_4_8) { + printf("Calculated bin lim for pt - 1st range - differs from expected!\n"); + } + for(Int_t i=0; i<=nbin0_4_8; i++) binLim0[i+nbin0_0_4]=(Double_t)ptmin_4_8 + (ptmax_4_8-ptmin_4_8)/nbin0_4_8*(Double_t)i ; + if (binLim0[nbin0_0_4+nbin0_4_8] != ptmin_8_10) { + printf("Calculated bin lim for pt - 2nd range - differs from expected!\n"); + } + for(Int_t i=0; i<=nbin0_8_10; i++) binLim0[i+nbin0_0_4+nbin0_4_8]=(Double_t)ptmin_8_10 + (ptmax_8_10-ptmin_8_10)/nbin0_8_10*(Double_t)i ; + + hmultiplypt->SetBins(nbins,binLim0); + hmultiplypt->Divide(hRECpt1D,hpteffCF,1,1,"b"); + hmultiplypt->SetMinimum(0); + hmultiplypt->SetMaximum(hMCpt1D->GetMaximum()); + hmultiplypt->SetLineColor(38); + hmultiplypt->SetLineWidth(3); + hmultiplypt->SetMarkerColor(38); + hmultiplypt->SetMarkerStyle(20); + hmultiplypt->GetXaxis()->SetTitleOffset(1.2); + hmultiplypt->GetYaxis()->SetTitleOffset(1.5); + hmultiplypt->GetXaxis()->SetTitle("p_{T} (GeV/c)"); + cmultiplypt->SetLeftMargin(0.15); + cmultiplypt->SetRightMargin(0.05); + hmultiplypt->SetFillStyle(3004); + hmultiplypt->SetFillColor(38); + cmultiplypt->cd(); + hmultiplypt->Draw("hist err"); + hMCpt1D->Draw("hist same err"); + + TH1D *hmultiplyy = new TH1D("hmultiplyy","hmultiplyy",8,-2,2); + hmultiplyy->Divide(hRECy1D,hyeffCF,1,1,"b"); + hmultiplyy->SetMinimum(0); + hmultiplyy->SetMaximum(hMCy1D->GetMaximum()); + hmultiplyy->SetLineColor(38); + hmultiplyy->SetLineWidth(3); + hmultiplyy->SetMarkerColor(38); + hmultiplyy->SetMarkerStyle(20); + hmultiplyy->GetXaxis()->SetTitleOffset(1.2); + hmultiplyy->GetXaxis()->SetTitle("y"); + hmultiplyy->GetYaxis()->SetTitleOffset(1.5); + cmultiplyy->SetLeftMargin(0.15); + cmultiplyy->SetRightMargin(0.05); + cmultiplyy->cd(); + hmultiplyy->SetFillStyle(3004); + hmultiplyy->SetFillColor(38); + hmultiplyy->Draw("hist err"); + hMCy1D->Draw("hist same err"); + + TH1D *hmultiplycTS = new TH1D("hmultiplycTS","hmultiplycTS",8,-1,1); + hmultiplycTS->Divide(hRECcTS1D,hcTSeffCF,1,1,"b"); + hmultiplycTS->SetMinimum(0); + hmultiplycTS->SetMaximum(hMCcTS1D->GetMaximum()); + hmultiplycTS->SetLineColor(38); + hmultiplycTS->SetLineWidth(3); + hmultiplycTS->SetMarkerColor(38); + hmultiplycTS->SetMarkerStyle(20); + hmultiplycTS->GetXaxis()->SetTitleOffset(1.2); + hmultiplycTS->GetXaxis()->SetTitle("cosThetaStar"); + hmultiplycTS->GetYaxis()->SetTitleOffset(1.5); + cmultiplycTS->SetLeftMargin(0.15); + cmultiplycTS->SetRightMargin(0.05); + cmultiplycTS->cd(); + hmultiplycTS->SetFillStyle(3004); + hmultiplycTS->SetFillColor(38); + hmultiplycTS->Draw("hist err"); + hMCcTS1D->Draw("hist same err"); + + TH1D *hmultiplyptPi = new TH1D("hmultiplyptPi","hmultiplyptPi",13,0,10); + hmultiplyptPi->SetBins(nbins,binLim0); + hmultiplyptPi->Divide(hRECptPi1D,hptPieffCF,1,1,"b"); + hmultiplyptPi->SetMinimum(0); + hmultiplyptPi->SetMaximum(hMCptPi1D->GetMaximum()); + hmultiplyptPi->SetLineColor(38); + hmultiplyptPi->SetLineWidth(3); + hmultiplyptPi->SetMarkerColor(38); + hmultiplyptPi->SetMarkerStyle(20); + hmultiplyptPi->GetXaxis()->SetTitleOffset(1.2); + hmultiplyptPi->GetXaxis()->SetTitle("p_{T, #pi} (GeV/c)"); + hmultiplyptPi->GetYaxis()->SetTitleOffset(1.5); + cmultiplyptPi->SetLeftMargin(0.15); + cmultiplyptPi->SetRightMargin(0.05); + cmultiplyptPi->cd(); + hmultiplyptPi->SetFillStyle(3004); + hmultiplyptPi->SetFillColor(38); + hmultiplyptPi->Draw("hist err"); + hMCptPi1D->Draw("hist same err"); + + TH1D *hmultiplyptK = new TH1D("hmultiplyptK","hmultiplyptK",13,0,10); + hmultiplyptK->SetBins(nbins,binLim0); + hmultiplyptK->Divide(hRECptK1D,hptKeffCF,1,1,"b"); + hmultiplyptK->SetMinimum(0); + hmultiplyptK->SetMaximum(hMCptK1D->GetMaximum()); + hmultiplyptK->SetLineColor(38); + hmultiplyptK->SetLineWidth(3); + hmultiplyptK->SetMarkerColor(38); + hmultiplyptK->SetMarkerStyle(20); + hmultiplyptK->GetXaxis()->SetTitleOffset(1.2); + hmultiplyptK->GetXaxis()->SetTitle("p_{T, K} (GeV/c)"); + hmultiplyptK->GetYaxis()->SetTitleOffset(1.5); + cmultiplyptK->SetLeftMargin(0.15); + cmultiplyptK->SetRightMargin(0.05); + cmultiplyptK->cd(); + hmultiplyptK->SetFillStyle(3004); + hmultiplyptK->SetFillColor(38); + hmultiplyptK->Draw("hist err"); + hMCptK1D->Draw("hist same err"); + + TH1D *hmultiplycT = new TH1D("hmultiplycT","hmultiplycT",24,0,500); + hmultiplycT->Divide(hRECcT1D,hcTeffCF,1,1,"b"); + hmultiplycT->SetMinimum(0); + hmultiplycT->SetMaximum(hMCcT1D->GetMaximum()); + hmultiplycT->SetLineColor(38); + hmultiplycT->SetLineWidth(3); + hmultiplycT->SetMarkerColor(38); + hmultiplycT->SetMarkerStyle(20); + hmultiplycT->GetXaxis()->SetTitleOffset(1.2); + hmultiplycT->GetXaxis()->SetTitle("c#t (#mum)"); + hmultiplycT->GetYaxis()->SetTitleOffset(1.5); + cmultiplycT->SetLeftMargin(0.15); + cmultiplycT->SetRightMargin(0.05); + cmultiplycT->cd(); + hmultiplycT->SetFillStyle(3004); + hmultiplycT->SetFillColor(38); + hmultiplycT->Draw("hist err"); + hMCcT1D->Draw("hist same err"); + + /* + TFile* file_histo = new TFile("fileHisto_180004.root","RECREATE"); + hMCpt1D->Write("hMCpt1D"); + hRECpt1D->Write("hRECpt1D"); + hMCy1D->Write("hMCy1D"); + hRECy1D->Write("hRECy1D"); + hMCcTS1D->Write("hMCcTS1D"); + hRECcTS1D->Write("hRECcTS1D"); + hMCptPi1D->Write("hMCptPi1D"); + hRECptPi1D->Write("hRECptPi1D"); + hMCptK1D->Write("hMCptK1D"); + hRECptK1D->Write("hRECptK1D"); + hMCcT1D->Write("hMCcT1D"); + hRECcT1D->Write("hRECcT1D"); + hpteffCF->Write("hpteffCF"); + hyeffCF->Write("hyeffCF"); + hcTSeffCF->Write("hcTSeffCF"); + hptPieffCF->Write("hptPieffCF"); + hptKeffCF->Write("hptKeffCF"); + hcTeffCF->Write("hcTeffCF"); + hmultiplypt->Write(""); + hmultiplyy->Write(""); + hmultiplycTS->Write(""); + hmultiplyptPi->Write(""); + hmultiplyptK->Write(""); + hmultiplycT->Write(""); + file_histo->Close(); + */ + + Int_t nsliceVars = 5; + Int_t sliceVars[5]; + sliceVars[0] = 1; + sliceVars[1] = 2; + sliceVars[2] = 3; + sliceVars[3] = 4; + sliceVars[4] = 5; + Double_t varMin[6] = {0.5, -2, -1, 0., 0., 0. }; + Double_t varMax[6] = {1.5, 2, 1, 10., 10., 500.}; + AliCFContainer* data_sliced = data->MakeSlice(nsliceVars, sliceVars, varMin, varMax); + cout << " the container now has " << data_sliced->GetNVar() << " dimensions " << endl; + +}