vertexingHF/AliAnalysisTaskSESelectHF.cxx
vertexingHF/AliAnalysisTaskSECompareHF.cxx
vertexingHF/AliCFHeavyFlavourTask.cxx
+ vertexingHF/AliCFHeavyFlavourTaskMultiVar.cxx
vertexingHF/AliMultiDimVector.cxx
vertexingHF/AliSignificanceCalculator.cxx
)
#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+;
vertexingHF/AliAnalysisTaskSESelectHF.cxx \
vertexingHF/AliAnalysisTaskSECompareHF.cxx \
vertexingHF/AliCFHeavyFlavourTask.cxx \
+ vertexingHF/AliCFHeavyFlavourTaskMultiVar.cxx \
vertexingHF/AliMultiDimVector.cxx vertexingHF/AliSignificanceCalculator.cxx
HDRS:= $(SRCS:.cxx=.h)
--- /dev/null
+//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");
+
+}
--- /dev/null
+/**************************************************************************
+ * 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 <TCanvas.h>
+#include <TParticle.h>
+#include <TH1I.h>
+#include <TStyle.h>
+
+#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<AliAODEvent*>(fInputEvent);
+ fCFManager->SetEventInfo(aodEvent);
+
+ // MC-event selection
+ Double_t containerInput[6] ;
+
+ //loop on the MC event
+
+ TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(AliAODMCParticle::StdBranchName()));
+ if (!mcArray) AliError("Could not find Monte-Carlo in AOD");
+ Int_t icountMC = 0;
+
+ for (Int_t iPart=0; iPart<mcArray->GetEntriesFast(); iPart++) {
+ AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(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; iVertex<arrayVerticesHF->GetEntriesFast(); 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<AliCFContainer*> (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<AliAODEvent*>(fInputEvent);
+ TClonesArray* mcArray = dynamic_cast<TClonesArray*>(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<AliAODMCParticle*>(mcArray->At(daughter0));
+ AliAODMCParticle* mcPartDaughter1 = dynamic_cast<AliAODMCParticle*>(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;
+}
+
--- /dev/null
+#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
--- /dev/null
+#include <Riostream.h>
+
+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;
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff