From 572b0139d695376c9ce5c1d238562920d63cd853 Mon Sep 17 00:00:00 2001 From: andronic Date: Tue, 4 May 2010 09:29:41 +0000 Subject: [PATCH] Updates and additions: Classes for signal and spectrum extraction; saving of a debug tree (default off); AddTaskJPSI.C; config macros for data --- PWG3/dielectron/AddTaskJPSI.C | 57 +++ .../AliAnalysisTaskMultiDielectron.cxx | 52 ++- .../AliAnalysisTaskMultiDielectron.h | 10 +- PWG3/dielectron/AliDielectron.cxx | 48 ++- PWG3/dielectron/AliDielectron.h | 10 +- PWG3/dielectron/AliDielectronCF.cxx | 174 ++++++--- PWG3/dielectron/AliDielectronCF.h | 11 +- PWG3/dielectron/AliDielectronCFdraw.cxx | 134 +++++-- PWG3/dielectron/AliDielectronCFdraw.h | 5 +- PWG3/dielectron/AliDielectronDebugTree.cxx | 161 ++++++++ PWG3/dielectron/AliDielectronDebugTree.h | 65 ++++ PWG3/dielectron/AliDielectronHistos.cxx | 156 +++++++- PWG3/dielectron/AliDielectronHistos.h | 14 +- PWG3/dielectron/AliDielectronMC.cxx | 8 +- PWG3/dielectron/AliDielectronMC.h | 4 + PWG3/dielectron/AliDielectronPair.cxx | 19 +- PWG3/dielectron/AliDielectronPair.h | 24 +- PWG3/dielectron/AliDielectronSignalBase.cxx | 96 +++++ PWG3/dielectron/AliDielectronSignalBase.h | 157 ++++++++ PWG3/dielectron/AliDielectronSignalExt.cxx | 248 +++++++++++++ PWG3/dielectron/AliDielectronSignalExt.h | 77 ++++ PWG3/dielectron/AliDielectronSignalFunc.cxx | 349 ++++++++++++++++++ PWG3/dielectron/AliDielectronSignalFunc.h | 80 ++++ PWG3/dielectron/AliDielectronSpectrum.cxx | 349 ++++++++++++++++++ PWG3/dielectron/AliDielectronSpectrum.h | 114 ++++++ PWG3/dielectron/AliDielectronVarCuts.cxx | 17 +- PWG3/dielectron/AliDielectronVarCuts.h | 12 +- PWG3/dielectron/AliDielectronVarManager.cxx | 13 +- PWG3/dielectron/AliDielectronVarManager.h | 108 +++++- PWG3/dielectron/macros/ConfigJpsi2eeData.C | 221 +++++++++++ PWG3/dielectron/macros/ConfigJpsi2eeDataAOD.C | 200 ++++++++++ 31 files changed, 2850 insertions(+), 143 deletions(-) create mode 100644 PWG3/dielectron/AddTaskJPSI.C create mode 100644 PWG3/dielectron/AliDielectronDebugTree.cxx create mode 100644 PWG3/dielectron/AliDielectronDebugTree.h create mode 100644 PWG3/dielectron/AliDielectronSignalBase.cxx create mode 100644 PWG3/dielectron/AliDielectronSignalBase.h create mode 100644 PWG3/dielectron/AliDielectronSignalExt.cxx create mode 100644 PWG3/dielectron/AliDielectronSignalExt.h create mode 100644 PWG3/dielectron/AliDielectronSignalFunc.cxx create mode 100644 PWG3/dielectron/AliDielectronSignalFunc.h create mode 100644 PWG3/dielectron/AliDielectronSpectrum.cxx create mode 100644 PWG3/dielectron/AliDielectronSpectrum.h create mode 100644 PWG3/dielectron/macros/ConfigJpsi2eeData.C create mode 100644 PWG3/dielectron/macros/ConfigJpsi2eeDataAOD.C diff --git a/PWG3/dielectron/AddTaskJPSI.C b/PWG3/dielectron/AddTaskJPSI.C new file mode 100644 index 00000000000..9dc535cd968 --- /dev/null +++ b/PWG3/dielectron/AddTaskJPSI.C @@ -0,0 +1,57 @@ +AliAnalysisTask *AddTaskJPSI(const char* config=""){ + //get the current analysis manager + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + ::Error("AddTasJPSI", "No analysis manager found."); + return NULL; + } + if (!mgr->GetInputEventHandler()) { + ::Error("AddTaskJPSI", "This task requires an input event handler"); + return NULL; + } + + TString configFile("$ALICE_ROOT/PWG3/dielectron/macros/ConfigJpsi2eeData.C"); + if (mgr->GetInputEventHandler()->IsA()==AliAODInputHandler::Class()){ + ::Info("AddTaskJPSI", "Using AOD configuration"); + configFile="$ALICE_ROOT/PWG3/dielectron/macros/ConfigJpsi2eeDataAOD.C"; + } + + //create task and add it to the manager + AliAnalysisTaskMultiDielectron *task=new AliAnalysisTaskMultiDielectron("MultiDie"); + mgr->AddTask(task); + + //load dielectron configuration file + gROOT->LoadMacro(configFile.Data()); + + //add dielectron analysis with different cuts to the task + for (Int_t i=0; iAddDielectron(jpsi); + } + + //---------------------- + //create data containers + //---------------------- + + //find input container + AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer(); + + TString containerName = mgr->GetCommonFileName(); + containerName += ":PWG3_dielectron"; + + //create output container + + AliAnalysisDataContainer *cOutputHist1 = + mgr->CreateContainer("QA", TList::Class(), AliAnalysisManager::kOutputContainer, + containerName.Data()); + + AliAnalysisDataContainer *cOutputHist2 = + mgr->CreateContainer("CF", TList::Class(), AliAnalysisManager::kOutputContainer, + containerName.Data()); + + mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer()); + mgr->ConnectOutput(task, 1, cOutputHist1); + mgr->ConnectOutput(task, 2, cOutputHist2); + + return task; +} diff --git a/PWG3/dielectron/AliAnalysisTaskMultiDielectron.cxx b/PWG3/dielectron/AliAnalysisTaskMultiDielectron.cxx index d2fa05fb5ec..da4c12ef056 100644 --- a/PWG3/dielectron/AliAnalysisTaskMultiDielectron.cxx +++ b/PWG3/dielectron/AliAnalysisTaskMultiDielectron.cxx @@ -22,11 +22,15 @@ #include #include +#include +#include +#include #include #include "AliDielectron.h" #include "AliDielectronHistos.h" #include "AliDielectronCF.h" +#include "AliDielectronMC.h" #include "AliAnalysisTaskMultiDielectron.h" ClassImp(AliAnalysisTaskMultiDielectron) @@ -36,7 +40,8 @@ AliAnalysisTaskMultiDielectron::AliAnalysisTaskMultiDielectron() : AliAnalysisTaskSE(), fListDielectron(), fListHistos(), - fListCF() + fListCF(), + fSelectPhysics(kFALSE) { // // Constructor @@ -48,7 +53,8 @@ AliAnalysisTaskMultiDielectron::AliAnalysisTaskMultiDielectron(const char *name) AliAnalysisTaskSE(name), fListDielectron(), fListHistos(), - fListCF() + fListCF(), + fSelectPhysics(kFALSE) { // // Constructor @@ -68,7 +74,7 @@ void AliAnalysisTaskMultiDielectron::UserCreateOutputObjects() // Add all histogram manager histogram lists to the output TList // - if (!fListHistos.IsEmpty()) return; //already initialised + if (!fListHistos.IsEmpty()||!fListCF.IsEmpty()) return; //already initialised TIter nextDie(&fListDielectron); AliDielectron *die=0; @@ -86,7 +92,32 @@ void AliAnalysisTaskMultiDielectron::UserExec(Option_t *) // Main loop. Called for every event // - if (fListHistos.IsEmpty()) return; + if (fListHistos.IsEmpty()&&fListCF.IsEmpty()) return; + + AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager(); + AliESDInputHandler *esdHandler=0x0; + if ( (esdHandler=dynamic_cast(man->GetInputEventHandler())) && esdHandler->GetESDpid() ){ + AliDielectronVarManager::SetESDpid(esdHandler->GetESDpid()); + } else { + //load esd pid bethe bloch parameters depending on the existance of the MC handler + // yes: MC parameters + // no: data parameters + if (!AliDielectronVarManager::GetESDpid()){ + if (AliDielectronMC::Instance()->HasMC()) { + AliDielectronVarManager::InitESDpid(); + } else { + AliDielectronVarManager::InitESDpid(1); + } + } + } + // Was event selected ? + AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler()); + Bool_t isSelected = kTRUE; + if( fSelectPhysics && inputHandler && inputHandler->GetEventSelection() ) { + isSelected = inputHandler->IsEventSelected(); + } + + if (!isSelected) return; //bz for AliKF Double_t bz = InputEvent()->GetMagneticField(); @@ -103,3 +134,16 @@ void AliAnalysisTaskMultiDielectron::UserExec(Option_t *) PostData(2, &fListCF); } +//_________________________________________________________________________________ +void AliAnalysisTaskMultiDielectron::FinishTaskOutput() +{ + // + // Write debug tree + // + TIter nextDie(&fListDielectron); + AliDielectron *die=0; + while ( (die=static_cast(nextDie())) ){ + die->SaveDebugTree(); + } +} + diff --git a/PWG3/dielectron/AliAnalysisTaskMultiDielectron.h b/PWG3/dielectron/AliAnalysisTaskMultiDielectron.h index 7762bf2bb66..a9942d715a7 100644 --- a/PWG3/dielectron/AliAnalysisTaskMultiDielectron.h +++ b/PWG3/dielectron/AliAnalysisTaskMultiDielectron.h @@ -26,11 +26,13 @@ class AliAnalysisTaskMultiDielectron : public AliAnalysisTaskSE { public: AliAnalysisTaskMultiDielectron(); AliAnalysisTaskMultiDielectron(const char *name); - virtual ~AliAnalysisTaskMultiDielectron(){} + virtual ~AliAnalysisTaskMultiDielectron(){ } - virtual void UserExec(Option_t *option); - virtual void UserCreateOutputObjects(); + virtual void UserExec(Option_t *option); + virtual void UserCreateOutputObjects(); + virtual void FinishTaskOutput(); + void UsePhysicsSelection(Bool_t phy=kTRUE) {fSelectPhysics=phy;} void AddDielectron(AliDielectron * const die) { fListDielectron.Add(die); } @@ -39,6 +41,8 @@ private: TList fListDielectron; // List of dielectron framework instances TList fListHistos; //! List of histogram manager lists in the framework classes TList fListCF; //! List with CF Managers + + Bool_t fSelectPhysics; // Whether to use physics selection AliAnalysisTaskMultiDielectron(const AliAnalysisTaskMultiDielectron &c); AliAnalysisTaskMultiDielectron& operator= (const AliAnalysisTaskMultiDielectron &c); diff --git a/PWG3/dielectron/AliDielectron.cxx b/PWG3/dielectron/AliDielectron.cxx index e95e953f664..a76cb2e92e7 100644 --- a/PWG3/dielectron/AliDielectron.cxx +++ b/PWG3/dielectron/AliDielectron.cxx @@ -57,6 +57,8 @@ The names are available via the function PairClassName(Int_t i) #include "AliDielectronCF.h" #include "AliDielectronMC.h" #include "AliDielectronVarManager.h" +#include "AliDielectronDebugTree.h" + #include "AliDielectron.h" ClassImp(AliDielectron) @@ -90,7 +92,8 @@ AliDielectron::AliDielectron() : fPdgMother(443), fHistos(0x0), fPairCandidates(new TObjArray(10)), - fCfManagerPair(0x0) + fCfManagerPair(0x0), + fDebugTree(0x0) { // // Default constructor @@ -107,7 +110,8 @@ AliDielectron::AliDielectron(const char* name, const char* title) : fPdgMother(443), fHistos(0x0), fPairCandidates(new TObjArray(10)), - fCfManagerPair(0x0) + fCfManagerPair(0x0), + fDebugTree(0x0) { // // Named constructor @@ -123,6 +127,7 @@ AliDielectron::~AliDielectron() // if (fHistos) delete fHistos; if (fPairCandidates) delete fPairCandidates; + if (fDebugTree) delete fDebugTree; } //________________________________________________________________ @@ -132,8 +137,8 @@ void AliDielectron::Init() // Initialise objects // if (fCfManagerPair) fCfManagerPair->InitialiseContainer(fPairFilter); - AliDielectronVarManager::InitESDpid(); //currently this will take the values for MC -} + +} //________________________________________________________________ void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2) @@ -142,10 +147,13 @@ void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2) // Process the events // - AliDielectronVarManager::SetEvent(0x0); + AliDielectronVarManager::SetEvent(ev1); //in case we have MC load the MC event and process the MC particles - if (AliDielectronMC::Instance()->ConnectMCEvent()) ProcessMC(); + if (AliDielectronMC::Instance()->HasMC()) { + AliDielectronMC::Instance()->ConnectMCEvent(); + ProcessMC(); + } //if candidate array doesn't exist, create it if (!fPairCandidates->UncheckedAt(0)) { @@ -178,7 +186,9 @@ void AliDielectron::Process(AliVEvent *ev1, AliVEvent *ev2) //in case there is a histogram manager, fill the QA histograms if (fHistos) FillHistograms(ev1); - + + //fill debug tree if a manager is attached + if (fDebugTree) FillDebugTree(); } //________________________________________________________________ @@ -310,6 +320,28 @@ void AliDielectron::FillPairArrays(Int_t arr1, Int_t arr2) { delete candidate; } +//________________________________________________________________ +void AliDielectron::FillDebugTree() +{ + // + // Fill Histogram information for tracks and pairs + // + + //Fill Debug tree + for (Int_t i=0; i<10; ++i){ + Int_t ntracks=PairArray(i)->GetEntriesFast(); + for (Int_t ipair=0; ipairFill(static_cast(PairArray(i)->UncheckedAt(ipair))); + } + } +} - +//________________________________________________________________ +void AliDielectron::SaveDebugTree() +{ + // + // delete the debug tree, this will also write the tree + // + if (fDebugTree) fDebugTree->DeleteStreamer(); +} diff --git a/PWG3/dielectron/AliDielectron.h b/PWG3/dielectron/AliDielectron.h index c3c625038b5..6f6f8267283 100644 --- a/PWG3/dielectron/AliDielectron.h +++ b/PWG3/dielectron/AliDielectron.h @@ -28,6 +28,7 @@ class AliVEvent; class THashList; class AliDielectronCF; +class AliDielectronDebugTree; //________________________________________________________________ class AliDielectron : public TNamed { @@ -68,9 +69,13 @@ public: void SetCFManagerPair(AliDielectronCF * const cf) { fCfManagerPair=cf; } AliDielectronCF* GetCFManagerPair() const { return fCfManagerPair; } + + void SetDebugTree(AliDielectronDebugTree * const tree) { fDebugTree=tree; } static const char* TrackClassName(Int_t i) { return (i>=0&&i<4)?fgkTrackClassNames[i]:""; } static const char* PairClassName(Int_t i) { return (i>=0&&i<10)?fgkPairClassNames[i]:""; } + + void SaveDebugTree(); private: @@ -96,6 +101,8 @@ private: AliDielectronCF *fCfManagerPair;//Correction Framework Manager for the Pair + AliDielectronDebugTree *fDebugTree; // Debug tree output + void FillTrackArrays(AliVEvent * const ev, Int_t eventNr=0); void FillPairArrays(Int_t arr1, Int_t arr2); @@ -112,11 +119,12 @@ private: void ProcessMC(); void FillHistograms(const AliVEvent *ev); + void FillDebugTree(); AliDielectron(const AliDielectron &c); AliDielectron &operator=(const AliDielectron &c); - ClassDef(AliDielectron,1); + ClassDef(AliDielectron,2); }; inline void AliDielectron::InitPairCandidateArrays() diff --git a/PWG3/dielectron/AliDielectronCF.cxx b/PWG3/dielectron/AliDielectronCF.cxx index e377cffa4bf..02e689ae929 100644 --- a/PWG3/dielectron/AliDielectronCF.cxx +++ b/PWG3/dielectron/AliDielectronCF.cxx @@ -56,9 +56,13 @@ AliDielectronCF::AliDielectronCF() : fStepForAfterAllCuts(kTRUE), fStepsForEachCut(kFALSE), fStepsForCutsIncreasing(kFALSE), + fStepsForSignal(kTRUE), + fStepsForBackground(kFALSE), fNStepMasks(0), fPdgMother(-1), - fCfContainer(0x0) + fCfContainer(0x0), + fHasMC(kFALSE), + fNAddSteps(0) { // // Default constructor @@ -91,9 +95,13 @@ AliDielectronCF::AliDielectronCF(const char* name, const char* title) : fStepForAfterAllCuts(kTRUE), fStepsForEachCut(kFALSE), fStepsForCutsIncreasing(kFALSE), + fStepsForSignal(kTRUE), + fStepsForBackground(kFALSE), fNStepMasks(0), fPdgMother(-1), - fCfContainer(0x0) + fCfContainer(0x0), + fHasMC(kFALSE), + fNAddSteps(0) { // // Named constructor @@ -153,15 +161,28 @@ void AliDielectronCF::InitialiseContainer(const AliAnalysisFilter& filter) fNCuts=filter.GetCuts()->GetEntries(); + fHasMC=AliDielectronMC::Instance()->HasMC(); + fNAddSteps=1; + if (fHasMC){ + if (fStepsForSignal) ++fNAddSteps; + if (fStepsForBackground) ++fNAddSteps; + } else { + //if + fStepForMCtruth=kFALSE; + fStepForNoCutsMCmotherPid=kFALSE; + fStepsForSignal=kFALSE; + fStepsForBackground=kFALSE; + } + fNSteps=0; if (fStepForMCtruth) ++fNSteps; if (fStepForNoCutsMCmotherPid) ++fNSteps; - if (fStepForAfterAllCuts) fNSteps+=2; - - if (fStepsForEachCut&&fNCuts>1) fNSteps+=(2*fNCuts); //one step for each cut + Signal (MC) - if (fStepsForCutsIncreasing&&fNCuts>2) fNSteps+=(2*(fNCuts-2)); //one step for the increasing cuts + Signal (MC) + if (fStepForAfterAllCuts) fNSteps+=fNAddSteps; + + if (fStepsForEachCut&&fNCuts>1) fNSteps+=(fNAddSteps*fNCuts); //one step for each cut + Signal (MC) + if (fStepsForCutsIncreasing&&fNCuts>2) fNSteps+=(fNAddSteps*(fNCuts-2)); //one step for the increasing cuts + Signal (MC) // e.g. cut2&cut3, cut2&cut3&cut4 - fNSteps+=(2*fNStepMasks); // cuts for the additional cut masks + fNSteps+=(fNAddSteps*fNStepMasks); // cuts for the additional cut masks // create the container Int_t *nbins=new Int_t[fNVars+2*fNVarsLeg]; for (Int_t i=0;iSetStepTitle(step++,"No cuts (Signal)"); } - //After All cuts TString cutName; - if (fStepForAfterAllCuts){ - cutName="All Cuts"; //TODO: User GetTitle??? - fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut - cutName+=" (Signal)"; - fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut with MC truth - } - //Steps for each of the cuts if (fStepsForEachCut&&fNCuts>1){ for (Int_t iCut=0; iCutAt(iCut)->GetName(); //TODO: User GetTitle??? + fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut - cutName+=" (Signal)"; - fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut with MC truth + + if (fHasMC){ + if (fStepsForSignal) + fCfContainer->SetStepTitle(step++, (cutName+" (Signal)").Data()); //Step for the cut with MC truth + if (fStepsForBackground) + fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data()); //Step for the cut with MC truth + } } } @@ -252,9 +271,15 @@ void AliDielectronCF::InitialiseContainer(const AliAnalysisFilter& filter) for (Int_t iCut=1; iCutAt(iCut)->GetName(); + fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut - cutName+=" (Signal)"; - fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut with MC truth + + if (fHasMC){ + if (fStepsForSignal) + fCfContainer->SetStepTitle(step++, (cutName+" (Signal)").Data()); //Step for the cut with MC truth + if (fStepsForBackground) + fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data()); //Step for the cut with MC truth + } } } @@ -272,9 +297,34 @@ void AliDielectronCF::InitialiseContainer(const AliAnalysisFilter& filter) } } } + fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut - cutName+=" (Signal)"; - fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut with MC truth + + if (fHasMC){ + if (fStepsForSignal) + fCfContainer->SetStepTitle(step++, (cutName+" (Signal)").Data()); //Step for the cut with MC truth + if (fStepsForBackground) + fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data()); //Step for the cut with MC truth + } + } + + //After All cuts + if (fStepForAfterAllCuts){ + cutName="No pair cuts"; + if (filter.GetCuts()->At(0)){ + cutName=filter.GetCuts()->At(0)->GetName(); //TODO: User GetTitle??? + for (Int_t iCut=1; iCutAt(iCut)->GetName(); + } + fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut + } + if (fHasMC){ + if (fStepsForSignal) + fCfContainer->SetStepTitle(step++, (cutName+" (Signal)").Data()); //Step for the cut with MC truth + if (fStepsForBackground) + fCfContainer->SetStepTitle(step++, (cutName+" (Background)").Data()); //Step for the cut with MC truth + } } if (step!=fNSteps) { @@ -290,7 +340,7 @@ void AliDielectronCF::Fill(UInt_t mask, const AliDielectronPair *particle) // Bool_t isMCTruth=kFALSE; - if (fPdgMother>=0) isMCTruth=AliDielectronMC::Instance()->IsMotherPdg(particle,fPdgMother); + if (fHasMC) isMCTruth=AliDielectronMC::Instance()->IsMotherPdg(particle,fPdgMother); Double_t valuesPair[AliDielectronVarManager::kNMaxValues]; AliDielectronVarManager::Fill(particle,valuesPair); @@ -328,28 +378,25 @@ void AliDielectronCF::Fill(UInt_t mask, const AliDielectronPair *particle) ++step; } - //All cuts - if (fStepForAfterAllCuts){ - if (mask == selectedMask){ - fCfContainer->Fill(fValues,step); - ++step; - if (isMCTruth) fCfContainer->Fill(fValues,step); - ++step; - } else { - step+=2; - } - } - //Steps for each of the cuts if (fStepsForEachCut&&fNCuts>1){ for (Int_t iCut=0; iCutFill(fValues,step); ++step; - if (isMCTruth) fCfContainer->Fill(fValues,step); - ++step; + + if (fHasMC){ + if ( fStepsForSignal){ + if (isMCTruth) fCfContainer->Fill(fValues,step); + ++step; + } + if ( fStepsForBackground ){ + if (!isMCTruth) fCfContainer->Fill(fValues,step); + ++step; + } + } } else { - step+=2; + step+=fNAddSteps; } } } @@ -360,10 +407,19 @@ void AliDielectronCF::Fill(UInt_t mask, const AliDielectronPair *particle) if (mask&(1<<((iCut+1)-1))) { fCfContainer->Fill(fValues,step); ++step; - if (isMCTruth) fCfContainer->Fill(fValues,step); - ++step; + + if (fHasMC){ + if ( fStepsForSignal){ + if (isMCTruth) fCfContainer->Fill(fValues,step); + ++step; + } + if ( fStepsForBackground ){ + if (!isMCTruth) fCfContainer->Fill(fValues,step); + ++step; + } + } } else { - step+=2; + step+=fNAddSteps; } } } @@ -374,13 +430,47 @@ void AliDielectronCF::Fill(UInt_t mask, const AliDielectronPair *particle) if (mask&userMask) { fCfContainer->Fill(fValues,step); ++step; - if (isMCTruth) fCfContainer->Fill(fValues,step); + + if (fHasMC){ + if ( fStepsForSignal){ + if (isMCTruth) fCfContainer->Fill(fValues,step); + ++step; + } + if ( fStepsForBackground ){ + if (!isMCTruth) fCfContainer->Fill(fValues,step); + ++step; + } + } + } else { + step+=fNAddSteps; + } + } + + //All cuts + if (fStepForAfterAllCuts){ + if (mask == selectedMask){ + fCfContainer->Fill(fValues,step); ++step; + + if (fHasMC){ + if ( fStepsForSignal){ + if (isMCTruth) fCfContainer->Fill(fValues,step); + ++step; + } + if ( fStepsForBackground ){ + if (!isMCTruth) fCfContainer->Fill(fValues,step); + ++step; + } + } } else { - step+=2; + step+=fNAddSteps; } } + if (step!=fNSteps) { + AliError("Something went wrong in the step filling!!!"); + } + } //________________________________________________________________ diff --git a/PWG3/dielectron/AliDielectronCF.h b/PWG3/dielectron/AliDielectronCF.h index 609cd42834d..c7f9d6f86d8 100644 --- a/PWG3/dielectron/AliDielectronCF.h +++ b/PWG3/dielectron/AliDielectronCF.h @@ -42,6 +42,8 @@ public: void SetStepForAfterAllCuts(Bool_t steps=kTRUE) { fStepForAfterAllCuts=steps; } void SetStepsForEachCut(Bool_t steps=kTRUE) { fStepsForEachCut=steps; } void SetStepsForCutsIncreasing(Bool_t steps=kTRUE) { fStepsForCutsIncreasing=steps; } + void SetStepsForSignal(Bool_t steps=kTRUE) { fStepsForSignal=steps; } + void SetStepsForBackground(Bool_t steps=kTRUE) { fStepsForBackground=steps; } void SetPdgMother(Int_t pdg) { fPdgMother=pdg; } @@ -75,7 +77,7 @@ private: Int_t fNCuts; // Number of cuts in the filter concerned - Double_t *fValues; // Value array for filling the container + Double_t *fValues; //! Value array for filling the container Bool_t fStepForMCtruth; //create a step for the MC truth Bool_t fStepForNoCutsMCmotherPid; //create a step for before cuts, but with MC truth of the mother @@ -83,12 +85,17 @@ private: Bool_t fStepsForEachCut; //create steps for each cut? Bool_t fStepsForCutsIncreasing; //create steps for increasing cut combinatons? //e.g. cut1&cut2, cut1&cut2&cut3 ... - + Bool_t fStepsForSignal; //steps for pure signal + Bool_t fStepsForBackground; //steps for pure background + UInt_t fStepMasks[kNmaxAddSteps]; //steps for additional cut combinatons UInt_t fNStepMasks; //number of configured step masks Int_t fPdgMother; //Pdg code of MCtruth validation AliCFContainer* fCfContainer; //the CF container + + Bool_t fHasMC; //if MC info is available + Int_t fNAddSteps; //number of additional MC related steps per cut step AliDielectronCF(const AliDielectronCF &c); AliDielectronCF &operator=(const AliDielectronCF &c); diff --git a/PWG3/dielectron/AliDielectronCFdraw.cxx b/PWG3/dielectron/AliDielectronCFdraw.cxx index 323ae408a8a..5656f9cb5b0 100644 --- a/PWG3/dielectron/AliDielectronCFdraw.cxx +++ b/PWG3/dielectron/AliDielectronCFdraw.cxx @@ -36,9 +36,11 @@ #include #include #include +#include #include #include #include +#include #include #include #include @@ -57,7 +59,8 @@ ClassImp(AliDielectronCFdraw) AliDielectronCFdraw::AliDielectronCFdraw() : TNamed(), fCfContainer(0x0), - fEffGrid(0x0) + fEffGrid(0x0), + fVdata(0) { // // Ctor @@ -68,7 +71,8 @@ AliDielectronCFdraw::AliDielectronCFdraw() : AliDielectronCFdraw::AliDielectronCFdraw(const char*name, const char* title) : TNamed(name,title), fCfContainer(0x0), - fEffGrid(0x0) + fEffGrid(0x0), + fVdata(0) { // // Named Ctor @@ -80,7 +84,8 @@ AliDielectronCFdraw::AliDielectronCFdraw(const char*name, const char* title) : AliDielectronCFdraw::AliDielectronCFdraw(AliCFContainer *cont) : TNamed(cont->GetName(), cont->GetTitle()), fCfContainer(cont), - fEffGrid(new AliCFEffGrid("eff","eff",*cont)) + fEffGrid(new AliCFEffGrid("eff","eff",*cont)), + fVdata(0) { // // directly set the CF container @@ -92,7 +97,8 @@ AliDielectronCFdraw::AliDielectronCFdraw(AliCFContainer *cont) : AliDielectronCFdraw::AliDielectronCFdraw(const char* filename) : TNamed(), fCfContainer(0x0), - fEffGrid(0x0) + fEffGrid(0x0), + fVdata(0) { // // get CF container(s) from file 'filename' @@ -120,7 +126,7 @@ void AliDielectronCFdraw::SetCFContainers(const TSeqCollection *arr) fCfContainer=new AliCFContainer(GetName(), GetTitle(), nstep, 1, nbins); //delete unneeded steps - for (Int_t istep=0; istepGetGrid(istep); +// for (Int_t istep=0; istepGetGrid(istep); //add step to the new container Int_t istep=0; @@ -199,23 +205,21 @@ void AliDielectronCFdraw::UnsetRangeUser(Int_t ivar, const char* slices) if (arr->GetEntriesFast()==0){ // all slices in case of 0 entries for (Int_t istep=0; istepGetNStep(); ++istep){ - fCfContainer->SetRangeUser(ivar,fCfContainer->GetAxis(ivar,istep)->GetXmin(), - fCfContainer->GetAxis(ivar,istep)->GetXmax(),istep); + fCfContainer->GetAxis(ivar,istep)->SetRange(0,0); } } else { TIter next(arr); TObjString *ostr=0x0; while ( (ostr=static_cast(next())) ) { Int_t istep=ostr->GetString().Atoi(); - fCfContainer->SetRangeUser(ivar,fCfContainer->GetAxis(ivar,istep)->GetXmin(), - fCfContainer->GetAxis(ivar,istep)->GetXmax(),istep); + fCfContainer->GetAxis(ivar,istep)->SetRange(0,0); } } delete arr; } //________________________________________________________________ -void AliDielectronCFdraw::Draw(const Option_t* varnames, const char* slices, const char* opt) +void AliDielectronCFdraw::Draw(const Option_t* varnames, const char* opt, const char* slices) { // // Draw 'variables' of 'slices' @@ -247,13 +251,13 @@ void AliDielectronCFdraw::Draw(const Option_t* varnames, const char* slices, con switch (entries){ case 1: - Draw(ivar[0],slices,opt); + Draw(ivar[0],opt,slices); break; case 2: - Draw(ivar[1],ivar[0],slices,opt); + Draw(ivar[1],ivar[0],opt,slices); break; case 3: - Draw(ivar[2],ivar[1],ivar[0],slices,opt); + Draw(ivar[2],ivar[1],ivar[0],opt,slices); break; } delete arrVars; @@ -318,6 +322,7 @@ TObjArray* AliDielectronCFdraw::CollectHistosProj(Int_t dim, Int_t *vars, const arrHists=new TObjArray(fCfContainer->GetNStep()); for (Int_t istep=0; istepGetNStep(); ++istep){ TH1 *hproj=Project(dim,vars,istep); + if (!hproj) continue; hproj->SetName(Form("proj_%02d",istep)); hproj->SetTitle(fCfContainer->GetStepTitle(istep)); arrHists->Add(hproj); @@ -329,6 +334,7 @@ TObjArray* AliDielectronCFdraw::CollectHistosProj(Int_t dim, Int_t *vars, const while ( (ostr=static_cast(next())) ) { Int_t istep=ostr->GetString().Atoi(); TH1 *hproj=Project(dim,vars,istep); + if (!hproj) continue; hproj->SetName(Form("proj_%02d",istep)); hproj->SetTitle(fCfContainer->GetStepTitle(istep)); arrHists->Add(hproj); @@ -420,7 +426,9 @@ void AliDielectronCFdraw::DrawEfficiency(Int_t var, const char* nominators, Int_ const Int_t ndim=1; Int_t vars[ndim]={var}; TObjArray *arr=CollectHistosEff(ndim,vars,nominators,denominator,type); - Draw(arr,opt); + TString drawOpt=opt; + drawOpt+="eff"; + Draw(arr,drawOpt); delete arr; } @@ -433,7 +441,9 @@ void AliDielectronCFdraw::DrawEfficiency(Int_t var0, Int_t var1, const char* nom const Int_t ndim=2; Int_t vars[ndim]={var0,var1}; TObjArray *arr=CollectHistosEff(ndim,vars,nominators,denominator,type); - Draw(arr,opt); + TString drawOpt=opt; + drawOpt+="eff"; + Draw(arr,drawOpt); delete arr; } @@ -446,7 +456,9 @@ void AliDielectronCFdraw::DrawEfficiency(Int_t var0, Int_t var1, Int_t var2, con const Int_t ndim=3; Int_t vars[ndim]={var0,var1,var2}; TObjArray *arr=CollectHistosEff(ndim,vars,nominators,denominator,type); - Draw(arr,opt); + TString drawOpt=opt; + drawOpt+="eff"; + Draw(arr,drawOpt); delete arr; } @@ -464,25 +476,32 @@ TObjArray* AliDielectronCFdraw::CollectHistosEff(Int_t dim, Int_t *vars, const c if (arr->GetEntriesFast()==0){ // all slices in case of 0 entries arrHists=new TObjArray(fCfContainer->GetNStep()); + fVdata.ResizeTo(arrHists->GetSize()); for (Int_t istep=0; istepGetNStep(); ++istep){ fEffGrid->CalculateEfficiency(istep,denominator); TH1 *hproj=ProjectEff(dim,vars); + if (!hproj) continue; Float_t eff=fEffGrid->GetAverage(); + fVdata(istep)=eff; hproj->SetName(Form("eff_%02d/%02d",istep,denominator)); - hproj->SetTitle(Form("%s (%f)",fCfContainer->GetStepTitle(istep),eff)); + hproj->SetTitle(Form("%s (%.3f)",fCfContainer->GetStepTitle(istep),eff)); arrHists->Add(hproj); } } else { arrHists=new TObjArray(arr->GetEntriesFast()); TIter next(arr); TObjString *ostr=0x0; + fVdata.ResizeTo(arrHists->GetSize()); + Int_t count=0; while ( (ostr=static_cast(next())) ) { Int_t istep=ostr->GetString().Atoi(); fEffGrid->CalculateEfficiency(istep,denominator); TH1 *hproj=ProjectEff(dim,vars); + if (!hproj) continue; Float_t eff=fEffGrid->GetAverage(); + fVdata(count++)=eff; hproj->SetName(Form("eff_%02d/%02d",istep,denominator)); - hproj->SetTitle(Form("%s (%f)",fCfContainer->GetStepTitle(istep),eff)); + hproj->SetTitle(Form("%s (%.3f)",fCfContainer->GetStepTitle(istep),eff)); arrHists->Add(hproj); } } @@ -495,27 +514,34 @@ TObjArray* AliDielectronCFdraw::CollectHistosEff(Int_t dim, Int_t *vars, const c if (arr->GetEntriesFast()==0){ // all slices in case of 0 entries arrHists=new TObjArray(fCfContainer->GetNStep()); + fVdata.ResizeTo(arrHists->GetSize()); for (Int_t istep=0; istepGetNStep(); ++istep){ TH1 *hproj=Project(dim,vars,istep); + if (!hproj) continue; Float_t eff=0; if (entriesDen>0) eff=hproj->GetEffectiveEntries()/entriesDen; + fVdata(istep)=eff; hproj->Divide(hDen); hproj->SetName(Form("eff_%02d/%02d",istep,denominator)); - hproj->SetTitle(Form("%s (%f)",fCfContainer->GetStepTitle(istep),eff)); + hproj->SetTitle(Form("%s (%.3f)",fCfContainer->GetStepTitle(istep),eff)); arrHists->Add(hproj); } } else { arrHists=new TObjArray(arr->GetEntriesFast()); + fVdata.ResizeTo(arrHists->GetSize()); TIter next(arr); TObjString *ostr=0x0; + Int_t count=0; while ( (ostr=static_cast(next())) ) { Int_t istep=ostr->GetString().Atoi(); TH1 *hproj=Project(dim,vars,istep); + if (!hproj) continue; Float_t eff=0; if (entriesDen>0) eff=hproj->GetEffectiveEntries()/entriesDen; + fVdata(count++)=eff; hproj->Divide(hDen); hproj->SetName(Form("eff_%02d/%02d",istep,denominator)); - hproj->SetTitle(Form("%s (%f)",fCfContainer->GetStepTitle(istep),eff)); + hproj->SetTitle(Form("%s (%.3f)",fCfContainer->GetStepTitle(istep),eff)); arrHists->Add(hproj); } } @@ -555,22 +581,31 @@ void AliDielectronCFdraw::Draw(const TObjArray *arr, const char* opt) // TString optStr(opt); optStr.ToLower(); - Bool_t drawSame=optStr.Contains("same"); - Bool_t optLeg=optStr.Contains("leg"); - optStr.ReplaceAll("same",""); + Bool_t drawSame = optStr.Contains("same"); + Bool_t drawSamePlus = optStr.Contains("same+"); + Bool_t drawEff = optStr.Contains("eff"); + Bool_t optLeg = optStr.Contains("leg"); + Bool_t optScaleMax = optStr.Contains("max"); + + if (!drawSamePlus) optStr.ReplaceAll("same",""); + + optStr.ReplaceAll("+",""); + optStr.ReplaceAll("eff",""); + optStr.ReplaceAll("leg",""); + optStr.ReplaceAll("max",""); if (!gPad) new TCanvas; Int_t nPads = arr->GetEntriesFast(); if (nPads==0) return; - if (nPads==1){ - arr->UncheckedAt(0)->Draw(optStr.Data()); - return; - } +// if (nPads==1){ +// arr->UncheckedAt(0)->Draw(optStr.Data()); +// return; +// } TCanvas *c=gPad->GetCanvas(); - c->Clear(); + if (!gVirtualPS&&!drawSamePlus) c->Clear(); if (!drawSame){ @@ -584,24 +619,53 @@ void AliDielectronCFdraw::Draw(const TObjArray *arr, const char* opt) } } else { TLegend *leg=0; - if (optLeg) leg=new TLegend(.8,.3,.99,.9); + if (drawSamePlus){ + //find already existing legend to attach entries to it + TIter nextPrimitive(gPad->GetListOfPrimitives()); + TObject *o=0x0; + while ((o=nextPrimitive())) if (o->IsA()==TLegend::Class()) leg=(TLegend*)o; + } + + if (optLeg&&!leg) leg=new TLegend(.2,.3,.99,.9); + Int_t addColor=0; + if (leg) addColor=leg->GetNRows(); Int_t offset=20; if (nPads<7) offset=24; for (Int_t i=0; i(arr->UncheckedAt(i)); - hist->SetLineColor(i+1); + hist->SetLineColor(i+1+addColor); hist->SetLineWidth(2); - hist->SetMarkerColor(i+1); - hist->SetMarkerStyle(offset+i); + hist->SetMarkerColor(i+1+addColor); + hist->SetMarkerStyle(offset+i+addColor); hist->Draw(optStr.Data()); + + if (drawEff&&i==0&&!drawSamePlus) { + hist->GetYaxis()->SetRangeUser(0,2); + hist->SetYTitle("Rec. Signal / Gen. Signal"); + } + if (leg) leg->AddEntry(hist,hist->GetTitle(),"lp"); } if (leg){ leg->SetFillColor(10); - leg->SetY1(.9-nPads*.05); + leg->SetY1(.9-leg->GetNRows()*.05); + leg->SetY1NDC(.9-leg->GetNRows()*.05); leg->SetMargin(.1); - leg->Draw(); + if (!drawSamePlus) leg->Draw(); + } + if (optScaleMax){ + TIter nextPrimitive(gPad->GetListOfPrimitives()); + TObject *o=0x0; + TH1 *hfirst=0x0; + Float_t max=0; + while ((o=nextPrimitive())) if (o->InheritsFrom(TH1::Class())){ + TH1 *h=(TH1*)o; + if (!hfirst) hfirst=h; + if (h->GetMaximum()>max) max=h->GetMaximum(); + } + max*=1.1; + hfirst->SetMaximum(max); } } diff --git a/PWG3/dielectron/AliDielectronCFdraw.h b/PWG3/dielectron/AliDielectronCFdraw.h index b662f2fce99..142350b9251 100644 --- a/PWG3/dielectron/AliDielectronCFdraw.h +++ b/PWG3/dielectron/AliDielectronCFdraw.h @@ -22,6 +22,7 @@ #include +#include #include class TObjArray; @@ -72,11 +73,13 @@ public: TObjArray* CollectHistosEff(Int_t dim, Int_t *vars, const char* nominators, Int_t denominator, Int_t type=0); TH1* ProjectEff(Int_t ndim, Int_t *vars); - + const TVectorD& GetData() const {return fVdata;} void Draw(const TObjArray *arr, const char* opt=""); private: AliCFContainer *fCfContainer; // CF container AliCFEffGrid *fEffGrid; // Efficiency calculation + + TVectorD fVdata; // vector with data, like mean efficiencies AliDielectronCFdraw(const AliDielectronCFdraw &c); AliDielectronCFdraw &operator=(const AliDielectronCFdraw &c); diff --git a/PWG3/dielectron/AliDielectronDebugTree.cxx b/PWG3/dielectron/AliDielectronDebugTree.cxx new file mode 100644 index 00000000000..550570280ac --- /dev/null +++ b/PWG3/dielectron/AliDielectronDebugTree.cxx @@ -0,0 +1,161 @@ +/************************************************************************* +* 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. * +**************************************************************************/ + +/////////////////////////////////////////////////////////////////////////// +// Dielectron DebugTree // +// // +// // +/* +register variables from the variable manager. The output will be written +to a tree + +NOTE: Please use with extream care! Only for debugging and test purposes!!! + +*/ +// // +/////////////////////////////////////////////////////////////////////////// + +#include +#include +#include + +#include + +#include "AliDielectronPair.h" + +#include "AliDielectronDebugTree.h" + +ClassImp(AliDielectronDebugTree) + +AliDielectronDebugTree::AliDielectronDebugTree() : + TNamed(), + fFileName("jpsi_debug.root"), + fNVars(0), + fNVarsLeg(0), + fStreamer(0x0) +{ + // + // Default Constructor + // + for (Int_t i=0; iGetFile()->Write(); + delete fStreamer; + } +} + +//______________________________________________ +void AliDielectronDebugTree::Fill(AliDielectronPair *pair) +{ + // + // Fill configured variables to the tree + // + + //is there anything to fill + if (fNVars==0&&fNVarsLeg==0) return; + + //only in local mode!!! + AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager(); + if (man && man->GetAnalysisType()!=AliAnalysisManager::kLocalAnalysis) return; + + if (!fStreamer) fStreamer=new TTreeSRedirector(fFileName.Data()); + + Int_t var=0; + Double_t values[AliDielectronVarManager::kNMaxValues]; + Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues]; + Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues]; +// fill pair values + if (fNVars>0){ + AliDielectronVarManager::Fill(pair,values); + + for (Int_t i=0; iGetFile()->Write(); +} diff --git a/PWG3/dielectron/AliDielectronDebugTree.h b/PWG3/dielectron/AliDielectronDebugTree.h new file mode 100644 index 00000000000..bf46792f695 --- /dev/null +++ b/PWG3/dielectron/AliDielectronDebugTree.h @@ -0,0 +1,65 @@ +#ifndef ALIDIELECTRONDEBUGTREE_H +#define ALIDIELECTRONDEBUGTREE_H + +/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//############################################################# +//# # +//# Class AliDielectronDebugTree # +//# # +//# Authors: # +//# Anton Andronic, GSI / A.Andronic@gsi.de # +//# Ionut C. Arsene, GSI / I.C.Arsene@gsi.de # +//# Julian Book, Uni Ffm / Julian.Book@cern.ch # +//# Frederick Kramer, Uni Ffm, / Frederick.Kramer@cern.ch # +//# Magnus Mager, CERN / Magnus.Mager@cern.ch # +//# WooJin J. Park, GSI / W.J.Park@gsi.de # +//# Jens Wiechula, Uni HD / Jens.Wiechula@cern.ch # +//# # +//############################################################# + +#include +#include + +#include "AliDielectronVarManager.h" + +class TTreeSRedirector; +class AliDielectronPair; + +class AliDielectronDebugTree : public TNamed { +public: + AliDielectronDebugTree(); + AliDielectronDebugTree(const char*name, const char* title); + + virtual ~AliDielectronDebugTree(); + + void SetOutputFileName(const char* file) { fFileName=file; } + + void AddPairVariable(AliDielectronVarManager::ValueTypes type) { fVariables[fNVars++]=(Int_t)type; } + void AddLegVariable(AliDielectronVarManager::ValueTypes type) { fVariablesLeg[fNVarsLeg++]=(Int_t)type; } + + void Fill(AliDielectronPair *pair); + + void DeleteStreamer(); + void WriteTree(); +private: + TString fFileName; //output file name + + Int_t fNVars; //number of configured variables + Int_t fVariables[AliDielectronVarManager::kNMaxValues]; //configured variables + Int_t fNVarsLeg; //number of configured variables + Int_t fVariablesLeg[AliDielectronVarManager::kNMaxValues]; //configured variables for the legs + + TTreeSRedirector *fStreamer; //! Tree Redirector + + AliDielectronDebugTree(const AliDielectronDebugTree &c); + AliDielectronDebugTree &operator=(const AliDielectronDebugTree &c); + + + ClassDef(AliDielectronDebugTree,1) // Dielectron DebugTree +}; + + + +#endif diff --git a/PWG3/dielectron/AliDielectronHistos.cxx b/PWG3/dielectron/AliDielectronHistos.cxx index 8fab66bbb02..160d673ac5f 100644 --- a/PWG3/dielectron/AliDielectronHistos.cxx +++ b/PWG3/dielectron/AliDielectronHistos.cxx @@ -36,6 +36,8 @@ #include #include #include +#include +#include // #include #include "AliDielectronHistos.h" @@ -84,14 +86,25 @@ AliDielectronHistos::~AliDielectronHistos() //_____________________________________________________________________________ void AliDielectronHistos::UserHistogram(const char* histClass,const char *name, const char* title, Int_t nbinsX, Double_t xmin, Double_t xmax, - UInt_t valTypeX) + UInt_t valTypeX, Bool_t logBinX) { // // Default histogram creation 1D case // if (!IsHistogramOk(histClass,name)) return; + + Double_t *binLimX=0x0; + + if (logBinX) { + binLimX=MakeLogBinning(nbinsX, xmin, xmax); + } else { + binLimX=MakeLinBinning(nbinsX, xmin, xmax); + } + + TH1* hist=new TH1F(name,title,nbinsX,binLimX); + + delete [] binLimX; - TH1* hist=new TH1F(name,title,nbinsX,xmin,xmax); Bool_t isReserved=fReservedWords->Contains(histClass); if (isReserved) UserHistogramReservedWords(histClass, hist, valTypeX); @@ -103,14 +116,34 @@ void AliDielectronHistos::UserHistogram(const char* histClass,const char *name, void AliDielectronHistos::UserHistogram(const char* histClass,const char *name, const char* title, Int_t nbinsX, Double_t xmin, Double_t xmax, Int_t nbinsY, Double_t ymin, Double_t ymax, - UInt_t valTypeX, UInt_t valTypeY) + UInt_t valTypeX, UInt_t valTypeY, + Bool_t logBinX, Bool_t logBinY) { // // Default histogram creation 2D case // if (!IsHistogramOk(histClass,name)) return; - TH1* hist=new TH2F(name,title,nbinsX,xmin,xmax,nbinsY,ymin,ymax); + Double_t *binLimX=0x0; + Double_t *binLimY=0x0; + + if (logBinX) { + binLimX=MakeLogBinning(nbinsX, xmin, xmax); + } else { + binLimX=MakeLinBinning(nbinsX, xmin, xmax); + } + if (logBinY) { + binLimY=MakeLogBinning(nbinsY, ymin, ymax); + } else { + binLimY=MakeLinBinning(nbinsY, ymin, ymax); + } + + TH1* hist=new TH2F(name,title,nbinsX,binLimX,nbinsY,binLimY); + + delete [] binLimX; + delete [] binLimY; + + Bool_t isReserved=fReservedWords->Contains(histClass); if (isReserved) UserHistogramReservedWords(histClass, hist, valTypeX+100*valTypeY); @@ -124,14 +157,42 @@ void AliDielectronHistos::UserHistogram(const char* histClass,const char *name, Int_t nbinsX, Double_t xmin, Double_t xmax, Int_t nbinsY, Double_t ymin, Double_t ymax, Int_t nbinsZ, Double_t zmin, Double_t zmax, - UInt_t valTypeX, UInt_t valTypeY, UInt_t valTypeZ) + UInt_t valTypeX, UInt_t valTypeY, UInt_t valTypeZ, + Bool_t logBinX, Bool_t logBinY, Bool_t logBinZ) { // // Default histogram creation 3D case // if (!IsHistogramOk(histClass,name)) return; - TH1* hist=new TH3F(name,title,nbinsX,xmin,xmax,nbinsY,ymin,ymax,nbinsZ,zmin,zmax); + Double_t *binLimX=0x0; + Double_t *binLimY=0x0; + Double_t *binLimZ=0x0; + + if (logBinX) { + binLimX=MakeLogBinning(nbinsX, xmin, xmax); + } else { + binLimX=MakeLinBinning(nbinsX, xmin, xmax); + } + + if (logBinY) { + binLimY=MakeLogBinning(nbinsY, ymin, ymax); + } else { + binLimY=MakeLinBinning(nbinsY, ymin, ymax); + } + + if (logBinZ) { + binLimZ=MakeLogBinning(nbinsZ, zmin, zmax); + } else { + binLimZ=MakeLinBinning(nbinsZ, zmin, zmax); + } + + TH1* hist=new TH3F(name,title,nbinsX,binLimX,nbinsY,binLimY,nbinsZ,binLimZ); + + delete [] binLimX; + delete [] binLimY; + delete [] binLimZ; + Bool_t isReserved=fReservedWords->Contains(histClass); if (isReserved) UserHistogramReservedWords(histClass, hist, valTypeX+100*valTypeY+10000*valTypeZ); @@ -337,9 +398,20 @@ void AliDielectronHistos::Draw(const Option_t* option) drawStr.ToLower(); //options // Bool_t same=drawOpt.Contains("same"); //FIXME not yet implemented + + TCanvas *c=0x0; + if (gVirtualPS) { + if (!gPad){ + Error("Draw","When writing to a file you have to create a canvas before opening the file!!!"); + return; + } + c=gPad->GetCanvas(); +// c=new TCanvas; + } TIter nextClass(&fHistoList); THashList *classTable=0; + Bool_t first=kTRUE; while ( (classTable=(THashList*)nextClass()) ){ //optimised division Int_t nPads = classTable->GetEntries(); @@ -347,13 +419,21 @@ void AliDielectronHistos::Draw(const Option_t* option) Int_t nRows = (Int_t)TMath::Ceil( (Double_t)nPads/(Double_t)nCols ); //create canvas - TString canvasName; - canvasName.Form("c%s",classTable->GetName()); - TCanvas *c=(TCanvas*)gROOT->FindObject(canvasName.Data()); - if (!c) c=new TCanvas(canvasName.Data(),classTable->GetName()); - c->Clear(); - c->Divide(nCols,nRows); - + if (!gVirtualPS){ + TString canvasName; + canvasName.Form("c%s_%s",GetName(),classTable->GetName()); + c=(TCanvas*)gROOT->FindObject(canvasName.Data()); + if (!c) c=new TCanvas(canvasName.Data(),Form("%s: %s",GetName(),classTable->GetName())); + c->Clear(); + } else { + if (first){ + first=kFALSE; + } else { + c->Clear(); + } + } + if (nCols>1||nRows>1) c->Divide(nCols,nRows); + //loop over histograms and draw them TIter nextHist(classTable); Int_t iPad=0; @@ -361,10 +441,15 @@ void AliDielectronHistos::Draw(const Option_t* option) while ( (h=(TH1*)nextHist()) ){ TString drawOpt; if ( (h->InheritsFrom(TH2::Class())) ) drawOpt="colz"; - c->cd(++iPad); + if (nCols>1||nRows>1) c->cd(++iPad); + if ( TMath::Abs(h->GetXaxis()->GetBinWidth(1)-h->GetXaxis()->GetBinWidth(2))>1e-10 ) gPad->SetLogx(); + if ( TMath::Abs(h->GetYaxis()->GetBinWidth(1)-h->GetYaxis()->GetBinWidth(2))>1e-10 ) gPad->SetLogy(); + if ( TMath::Abs(h->GetZaxis()->GetBinWidth(1)-h->GetZaxis()->GetBinWidth(2))>1e-10 ) gPad->SetLogz(); h->Draw(drawOpt.Data()); } + if (gVirtualPS) c->Update(); } +// if (gVirtualPS) delete c; } //_____________________________________________________________________________ @@ -405,6 +490,8 @@ void AliDielectronHistos::SetHistogramList(THashList &list) // // set histogram classes and histograms to this instance. It will take onwnership! // + ResetHistogramList(); + SetName(list.GetName()); TIter next(&list); TObject *o; while ( (o=next()) ){ @@ -460,6 +547,7 @@ void AliDielectronHistos::ReadFromFile(const char* file) f.Close(); } +//_____________________________________________________________________________ void AliDielectronHistos::DrawSame(const char* histName, const Option_t *opt) { // @@ -520,3 +608,43 @@ void AliDielectronHistos::SetReservedWords(const char* words) (*fReservedWords)=words; } +//_____________________________________________________________________________ +Double_t* AliDielectronHistos::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) const +{ + // + // Make logarithmic binning + // the user has to delete the array afterwards!!! + // + + //check limits + if (xmin<1e-20 || xmax<1e-20){ + Error("MakeLogBinning","For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!"); + return MakeLinBinning(nbinsX, xmin, xmax); + } + Double_t *binLim=new Double_t[nbinsX+1]; + Double_t first=xmin; + Double_t last=xmax; + Double_t expMax=TMath::Log(last/first); + for (Int_t i=0; i list of histograms @@ -87,6 +90,9 @@ private: Bool_t IsHistogramOk(const char* classTable, const char* name); + Double_t* MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) const; + Double_t* MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) const; + enum {kNoAutoFill=1000000000}; AliDielectronHistos(const AliDielectronHistos &hist); diff --git a/PWG3/dielectron/AliDielectronMC.cxx b/PWG3/dielectron/AliDielectronMC.cxx index 625bd827583..e2b55fefaf5 100644 --- a/PWG3/dielectron/AliDielectronMC.cxx +++ b/PWG3/dielectron/AliDielectronMC.cxx @@ -49,8 +49,11 @@ AliDielectronMC* AliDielectronMC::Instance() AnalysisType type=kUNSET; if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliESDInputHandler::Class()) type=kESD; else if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliAODHandler::Class()) type=kAOD; + + AliMCEventHandler* mcHandler = dynamic_cast (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); fgInstance=new AliDielectronMC(type); + fgInstance->SetHasMC(mcHandler!=0x0); return fgInstance; } @@ -59,7 +62,8 @@ AliDielectronMC* AliDielectronMC::Instance() AliDielectronMC::AliDielectronMC(AnalysisType type): fMCEvent(0x0), fStack(0x0), - fAnaType(type) + fAnaType(type), + fHasMC(kTRUE) { // // default constructor @@ -397,6 +401,7 @@ Int_t AliDielectronMC::GetLabelMotherWithPdgESD(const AliVParticle *particle1, c if (!mcMother1) return -1; if (lblMother1!=lblMother2) return -1; if (TMath::Abs(mcPart1->PdgCode())!=11) return -1; + if (mcPart1->PdgCode()!=-mcPart2->PdgCode()) return -1; if (mcMother1->PdgCode()!=pdgMother) return -1; return lblMother1; @@ -423,6 +428,7 @@ Int_t AliDielectronMC::GetLabelMotherWithPdgAOD(const AliVParticle *particle1, c if (!mcMother1) return -1; if (lblMother1!=lblMother2) return -1; if (TMath::Abs(mcPart1->GetPdgCode())!=11) return -1; + if (mcPart1->GetPdgCode()!=-mcPart2->GetPdgCode()) return -1; if (mcMother1->GetPdgCode()!=pdgMother) return -1; return lblMother1; diff --git a/PWG3/dielectron/AliDielectronMC.h b/PWG3/dielectron/AliDielectronMC.h index fdf5ae930c4..0079dcaa9ad 100644 --- a/PWG3/dielectron/AliDielectronMC.h +++ b/PWG3/dielectron/AliDielectronMC.h @@ -34,6 +34,9 @@ public: AliDielectronMC(AnalysisType type=kUNSET); virtual ~AliDielectronMC(); + void SetHasMC(Bool_t hasMC) { fHasMC=hasMC; } + Bool_t HasMC() const { return fHasMC; } + static AliDielectronMC* Instance(); void Initialize(); // initialization @@ -72,6 +75,7 @@ private: AliStack *fStack; // MC stack AnalysisType fAnaType; // Analysis type + Bool_t fHasMC; // Do we have an MC handler? static AliDielectronMC* fgInstance; //! singleton pointer diff --git a/PWG3/dielectron/AliDielectronPair.cxx b/PWG3/dielectron/AliDielectronPair.cxx index 96d3e228965..ef30d2795b2 100644 --- a/PWG3/dielectron/AliDielectronPair.cxx +++ b/PWG3/dielectron/AliDielectronPair.cxx @@ -26,10 +26,11 @@ ClassImp(AliDielectronPair) AliDielectronPair::AliDielectronPair() : - fOpeningAngle(-1), fType(-1), fLabel(-1), fPair(), + fD1(), + fD2(), fRefD1(), fRefD2() { @@ -42,10 +43,11 @@ AliDielectronPair::AliDielectronPair() : //______________________________________________ AliDielectronPair::AliDielectronPair(AliVTrack * const particle1, Int_t pid1, AliVTrack * const particle2, Int_t pid2, Char_t type) : - fOpeningAngle(-1), fType(type), fLabel(-1), fPair(), + fD1(), + fD2(), fRefD1(), fRefD2() { @@ -69,24 +71,29 @@ void AliDielectronPair::SetTracks(AliVTrack * const particle1, Int_t pid1, AliVTrack * const particle2, Int_t pid2) { // - // Set the tracks to the pair KF particle + // Sort particles by pt, first particle larget Pt + // set AliKF daughters and pair // fPair.Initialize(); + fD1.Initialize(); + fD2.Initialize(); AliKFParticle kf1(*particle1,pid1); AliKFParticle kf2(*particle2,pid2); - + fPair.AddDaughter(kf1); fPair.AddDaughter(kf2); if (particle1->Pt()>particle2->Pt()){ fRefD1 = particle1; fRefD2 = particle2; + fD1+=kf1; + fD2+=kf2; } else { fRefD1 = particle2; fRefD2 = particle1; + fD1+=kf2; + fD2+=kf1; } - - fOpeningAngle=kf1.GetAngle(kf2); } diff --git a/PWG3/dielectron/AliDielectronPair.h b/PWG3/dielectron/AliDielectronPair.h index 64a496150da..fa5addd15eb 100644 --- a/PWG3/dielectron/AliDielectronPair.h +++ b/PWG3/dielectron/AliDielectronPair.h @@ -76,32 +76,42 @@ public: virtual const Double_t *PID() const { return 0;} //TODO: check // Dummy Int_t PdgCode() const {return 0;} - - Double_t OpeningAngle() const { return fOpeningAngle; } + UChar_t GetType() const { return fType; } void SetType(Char_t type) { fType=type; } void SetLabel(Int_t label) {fLabel=label;} + + //inter leg information + Double_t OpeningAngle() const { return fD1.GetAngle(fD2); } + Double_t DistanceDaughters() const { return fD1.GetDistanceFromParticle(fD2); } + Double_t DistanceDaughtersXY() const { return fD1.GetDistanceFromParticleXY(fD2); } + Double_t DeviationDaughters() const { return fD1.GetDeviationFromParticle(fD2); } + Double_t DeviationDaughtersXY() const { return fD1.GetDeviationFromParticleXY(fD2); } + // internal KF particle - const AliKFParticle& GetKFParticle() const { return fPair; } - + const AliKFParticle& GetKFParticle() const { return fPair; } + const AliKFParticle& GetKFFirstDaughter() const { return fD1; } + const AliKFParticle& GetKFSecondDaughter() const { return fD2; } + // daughter references AliVParticle* GetFirstDaughter() const { return dynamic_cast(fRefD1.GetObject()); } AliVParticle* GetSecondDaughter() const { return dynamic_cast(fRefD2.GetObject()); } private: - Double_t fOpeningAngle; // opening angle of the pair Char_t fType; // type of the pair e.g. like sign SE, unlike sign SE, ... see AliDielectron Int_t fLabel; // MC label AliKFParticle fPair; // KF particle internally used for pair calculation - + AliKFParticle fD1; // KF particle first daughter + AliKFParticle fD2; // KF particle1 second daughter + TRef fRefD1; // Reference to first daughter TRef fRefD2; // Reference to second daughter - ClassDef(AliDielectronPair,2) + ClassDef(AliDielectronPair,3) }; #endif diff --git a/PWG3/dielectron/AliDielectronSignalBase.cxx b/PWG3/dielectron/AliDielectronSignalBase.cxx new file mode 100644 index 00000000000..44adbd31f18 --- /dev/null +++ b/PWG3/dielectron/AliDielectronSignalBase.cxx @@ -0,0 +1,96 @@ +/************************************************************************* +* 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. * +**************************************************************************/ + +/////////////////////////////////////////////////////////////////////////// +// Dielectron SignalBase // +// // +// // +/* +Base class for signal extraction from a histogram or an array of histograms +The histogram is assumed to be an inv. mass spectrum, +the array of histograms is assumed to be an array with inv. mass histograms +resulting from single and mixed events, as defined in AliDielectron.cxx + +*/ +// // +/////////////////////////////////////////////////////////////////////////// + + +#include "TPaveText.h" +#include "AliDielectronSignalBase.h" + +ClassImp(AliDielectronSignalBase) + +AliDielectronSignalBase::AliDielectronSignalBase() : + TNamed(), + fValues(4), + fErrors(4), + fIntMin(2.99), + fIntMax(3.19) +{ + // + // Default Constructor + // + +} + +//______________________________________________ +AliDielectronSignalBase::AliDielectronSignalBase(const char* name, const char* title) : + TNamed(name, title), + fValues(4), + fErrors(4), + fIntMin(2.99), + fIntMax(3.19) +{ + // + // Named Constructor + // +} + +//______________________________________________ +AliDielectronSignalBase::~AliDielectronSignalBase() +{ + // + // Default Destructor + // + +} + +//______________________________________________ +TPaveText* AliDielectronSignalBase::DrawStats(Double_t x1/*=0.*/, Double_t y1/*=0.*/, Double_t x2/*=0.*/, Double_t y2/*=0.*/) +{ + // + // Draw extracted values in a TPaveText + // with the corners x1,y2,x2,y2 + // + if (TMath::Abs(x1)<1e-20&&TMath::Abs(x2)<1e-20){ + x1=.6; + x2=.9; + y1=.7; + y2=.9; + } + TPaveText *t=new TPaveText(x1,y1,x2,y2,"brNDC"); + t->SetFillColor(kWhite); + t->SetBorderSize(1); + t->SetTextAlign(12); + t->AddText(Form("Singal : %.2f #pm %.2f",GetSignal(),GetSignalError())); + t->AddText(Form("Backgnd: %.2f #pm %.2f",GetBackground(),GetBackgroundError())); + t->AddText(Form("Signif.: %.2f #pm %.2f",GetSignificance(),GetSignificanceError())); + t->AddText(Form("SoB : %.2f #pm %.2f",GetSignalOverBackground(),GetSignalOverBackgroundError())); + t->Draw(); + + return t; +} + diff --git a/PWG3/dielectron/AliDielectronSignalBase.h b/PWG3/dielectron/AliDielectronSignalBase.h new file mode 100644 index 00000000000..88260809810 --- /dev/null +++ b/PWG3/dielectron/AliDielectronSignalBase.h @@ -0,0 +1,157 @@ +#ifndef ALIDIELECTRONSIGNALBASE_H +#define ALIDIELECTRONSIGNALBASE_H + +/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//############################################################# +//# # +//# Class AliDielectronSignalBase # +//# Manage Cuts on the legs of the pair # +//# # +//# Authors: # +//# Anton Andronic, GSI / A.Andronic@gsi.de # +//# Ionut C. Arsene, GSI / I.C.Arsene@gsi.de # +//# Julian Book, Uni Ffm / Julian.Book@cern.ch # +//# Frederick Kramer, Uni Ffm, / Frederick.Kramer@cern.ch # +//# Magnus Mager, CERN / Magnus.Mager@cern.ch # +//# WooJin J. Park, GSI / W.J.Park@gsi.de # +//# Jens Wiechula, Uni HD / Jens.Wiechula@cern.ch # +//# # +//############################################################# + +#include +#include +#include + +class TObjArray; +class TPaveText; +class TH1; + +class AliDielectronSignalBase : public TNamed { +public: + AliDielectronSignalBase(); + AliDielectronSignalBase(const char*name, const char* title); + + virtual ~AliDielectronSignalBase(); + + + void SetIntegralRange(Double_t min, Double_t max) {fIntMin=min;fIntMax=max;} + + void SetSignal(Double_t val, Double_t valErr) {fValues(0)=val; fErrors(0)=valErr;} + void SetBackground(Double_t val, Double_t valErr) {fValues(1)=val; fErrors(1)=valErr;} + void SetSignificance(Double_t val, Double_t valErr) {fValues(2)=val; fErrors(2)=valErr;} + void SetSignalOverBackground(Double_t val, Double_t valErr) {fValues(3)=val; fErrors(3)=valErr;} + + const TVectorD& GetValues() const {return fValues;} + const TVectorD& GetErrors() const {return fErrors;} + + Double_t GetIntegralMin() const { return fIntMin; } + Double_t GetIntegralMax() const { return fIntMax; } + + Double_t GetSignal() {return fValues(0);} + Double_t GetBackground() {return fValues(1);} + Double_t GetSignificance() {return fValues(2);} + Double_t GetSignalOverBackground() {return fValues(3);} + + Double_t GetSignalError() {return fErrors(0);} + Double_t GetBackgroundError() {return fErrors(1);} + Double_t GetSignificanceError() {return fErrors(2);} + Double_t GetSignalOverBackgroundError() {return fErrors(3);} + + void GetSignal(Double_t &val, Double_t &valErr) {val=fValues(0); valErr=fErrors(0);} + void GetBackground(Double_t &val, Double_t &valErr) {val=fValues(1); valErr=fErrors(1);} + void GetSignificance(Double_t &val, Double_t &valErr) {val=fValues(2); valErr=fErrors(2);} + void GetSignalOverBackground(Double_t &val, Double_t &valErr) {val=fValues(3); valErr=fErrors(3);} + + /** + This function needs to be implemented by the signal extraction classes. + + The signal extraction is done on the mass spectra. + The TObjArray should contain the Inv. Mass spectra of the 10 possible combinations + for single and mixed events defined in AliDielectron.cxx + In the referece TVectorDs the values and value errors of the result should be stored: + Parameter - 0: Signal entries + 1: Background entries + 2: Significance ( S/sqr(S+B) ) + 3: Signal/Background + + It is enough to calculate the signal and background and then call + SetSignificanceAndSOB(TVectorD &values, TVectorD &errors) + Please also read the description there!!! + */ + virtual void Process(TObjArray * const /*arrhist*/) = 0; + +protected: + + void SetSignificanceAndSOB(); + TPaveText* DrawStats(Double_t x1=0., Double_t y1=0., Double_t x2=0., Double_t y2=0.); + void Reset() {fValues.Zero(); fErrors.Zero();} + +private: + TVectorD fValues; // values + TVectorD fErrors; // value errors + + Double_t fIntMin; // signal extraction range min + Double_t fIntMax; // signal extraction range max + + AliDielectronSignalBase(const AliDielectronSignalBase &c); + AliDielectronSignalBase &operator=(const AliDielectronSignalBase &c); + + ClassDef(AliDielectronSignalBase,2) // Dielectron SignalBase +}; + +// +// Inline functions +// + +inline void AliDielectronSignalBase::SetSignificanceAndSOB() +{ + // + // calculate significance and signal over background from values + // it is expected that: + // - 'values' and 'errors' have the size 4 + // - Parameters 0 are given: signal and signal error + // - Parameters 1 are given: background and background error + // - Optionally parameter 2 can be given: signal+background and its error + // + // The calculated significance and S/B and the corresponding errors + // are written in parameters 2 and 3, the first two parameters (signal and background) + // stay untouched + + Double_t signal=fValues(0), signal_err=fErrors(0); + Double_t background=fValues(1), background_err=fErrors(1); + Double_t sigPlusBack=fValues(2), sigPlusBack_err=fErrors(2); + Double_t significance=0, significance_err=0; + Double_t sob=0, sob_err=0; + + if (sigPlusBack<1e-20){ + //calculate signal plus background + sigPlusBack=signal+background; + sigPlusBack_err=TMath::Sqrt( signal_err*signal_err + background_err*background_err ); + } + + if (sigPlusBack>1e-30){ + + significance=signal/TMath::Sqrt(sigPlusBack); + + significance_err=TMath::Sqrt((signal_err/signal)*(signal_err/signal)+ + (0.5*sigPlusBack_err/sigPlusBack)*(0.5*sigPlusBack_err/sigPlusBack) + )*significance; + + } + + if (background>1e-30){ + sob=signal/background; + sob_err=TMath::Sqrt((signal_err/signal)*(signal_err/signal)+ + (background_err/background)*(background_err/background))*sob; + } + + fValues(2)=significance; + fErrors(2)=significance_err; + + fValues(3)=sob; + fErrors(3)=sob_err; +} + +#endif diff --git a/PWG3/dielectron/AliDielectronSignalExt.cxx b/PWG3/dielectron/AliDielectronSignalExt.cxx new file mode 100644 index 00000000000..0d12f9b71d9 --- /dev/null +++ b/PWG3/dielectron/AliDielectronSignalExt.cxx @@ -0,0 +1,248 @@ +/************************************************************************* +* 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. * +**************************************************************************/ + +/////////////////////////////////////////////////////////////////////////// +// // +// Dielectron SignalExt // +// // +/* +Dielectron signal extraction class using functions as input. + +A function to describe the signal as well as one to describe the background +has to be deployed by the user. Alternatively on of the default implementaions +can be used. + +*/ +// // +/////////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include +#include + +#include + +#include "AliDielectronSignalExt.h" + +ClassImp(AliDielectronSignalExt) + +AliDielectronSignalExt::AliDielectronSignalExt() : + AliDielectronSignalBase(), + fSignPM(0x0), + fSignPP(0x0), + fSignMM(0x0), + fBackground(0x0), + fSignal(0x0), + fMethod(1), + fRebin(1), + fBins(0), + fDrawMin(0.), + fDrawMax(0.), + fFitMin(2.5), + fFitMax(4) +{ + // + // Default Constructor + // + +} + +//______________________________________________ +AliDielectronSignalExt::AliDielectronSignalExt(const char* name, const char* title) : + AliDielectronSignalBase(name, title), + fSignPM(0x0), + fSignPP(0x0), + fSignMM(0x0), + fBackground(0x0), + fSignal(0x0), + fMethod(1), + fRebin(1), + fBins(0), + fDrawMin(0.), + fDrawMax(0.), + fFitMin(2.5), + fFitMax(4) +{ + // + // Named Constructor + // +} + +//______________________________________________ +AliDielectronSignalExt::~AliDielectronSignalExt() +{ + // + // Default Destructor + // + if (fSignPM) delete fSignPM; + if (fSignPP) delete fSignPP; + if (fSignMM) delete fSignMM; + if (fBackground) delete fBackground; + if (fSignal) delete fSignal; +} + +//______________________________________________ +void AliDielectronSignalExt::SetHistograms(TH1F* const unlike, TH1F* const backg, TH1F* const signal) +{ + // + // set histograms + // + fSignPM = (TH1F*)unlike->Clone("fSignPM"); + fBackground = backg; + fSignal = signal; +} + +//______________________________________________ +void AliDielectronSignalExt::Process(TObjArray* const arrhist) +{ + // + // signal subtraction. support like-sign subtraction and event mixing method + // + switch ( fMethod ){ + case 1 : + ProcessLS(arrhist); // process like-sign subtraction method + break; + + case 2 : + ProcessEM(arrhist); // process event mixing method + break; + + default : + AliWarning("Subtraction method not supported. Please check SetMethod() function."); + } +} + +//______________________________________________ +void AliDielectronSignalExt::ProcessLS(TObjArray* const arrhist) +{ + // + // signal subtraction + // + fSignPP = (TH1F*)arrhist->At(0); // like sign : plus-plus + fSignPM = (TH1F*)arrhist->At(1); // unlike-sign : plus-minus + fSignMM = (TH1F*)arrhist->At(2); // like sign : minus-minus + fSignPP->Sumw2(); + fSignPM->Sumw2(); + fSignMM->Sumw2(); + + if ( fRebin>1 ){ Rebin(fRebin); } // rebinning of histogram + + fBins = fSignPM->GetNbinsX(); // number of bins + Double_t minX = fSignPM->GetBinLowEdge(1); // minimum X value in axis + Double_t maxX = fSignPM->GetBinLowEdge(fBins+1); // maximum X value in axis + + AliDebug(10,Form("histogram #bin = %d , min = %f , max = %f\n",fBins,minX,maxX)); + TH1F* hBackground = new TH1F("hBackground","Like-sign background",fBins,minX,maxX); + TH1F* hSubtracted = new TH1F("hSubtracted","Background subtracted",fBins,minX,maxX); + + // fill out background and subtracted histogram + for ( Int_t ibin=1; ibinGetBinContent(ibin) * fSignMM->GetBinContent(ibin) ); + Double_t signal = fSignPM->GetBinContent(ibin) - backgr; + + hBackground->Fill(mass, backgr); + hSubtracted->Fill(mass, signal); + } + SetHistograms(fSignPM, hBackground, hSubtracted); + + + Double_t signal=0, signal_err=0; + Double_t background=0, background_err=0; + + signal = fSignal->IntegralAndError(fSignal->FindBin(GetIntegralMin()), + fSignal->FindBin(GetIntegralMax()), signal_err); + background = fBackground->IntegralAndError(fBackground->FindBin(GetIntegralMin()), + fBackground->FindBin(GetIntegralMax()), background_err); + + //reset result arrays + Reset(); + //set values + SetSignal(signal,signal_err); + SetBackground(background,background_err); + SetSignificanceAndSOB(); + + // cleanup + //delete hBackground; + //delete hSubtracted; +} + +//______________________________________________ +void AliDielectronSignalExt::ProcessEM(TObjArray* const arrhist) +{ + // + // event mixing. not implemented yet. + // + printf("event mixing method is not implemented yet. Like-sign method will be used.\n"); + ProcessLS(arrhist); +} + +//______________________________________________ +void AliDielectronSignalExt::Rebin(Int_t rebin) +{ + // + // rebinning of histograms + // + fSignPM->Rebin(rebin); + fSignPP->Rebin(rebin); + fSignMM->Rebin(rebin); +} + +//______________________________________________ +void AliDielectronSignalExt::Draw(const Option_t* option) +{ + // + // Draw the fitted function + // + TString drawOpt(option); + drawOpt.ToLower(); + + TCanvas *cSub = new TCanvas("cSub","signal, background subtracted",1400,1000); + cSub->Divide(2,2); + cSub->Draw(); + + Double_t minX = fSignPM->GetBinLowEdge(1); // minimum X value in axis + Double_t maxX = fSignPM->GetBinLowEdge(fBins+1); // maximum X value in axis + if ( TMath::Abs(fDrawMin)<1.e-30 ) fDrawMin = minX; + if ( TMath::Abs(fDrawMax)<1.e-30 ) fDrawMax = maxX; + + cSub->cd(4); + fSignPM->GetXaxis()->SetRangeUser(fDrawMin, fDrawMax); + fSignPM->SetLineColor(1); + fSignPM->SetLineWidth(2); + fSignPM->SetMarkerStyle(6); + fSignPM->DrawCopy("P"); + + fBackground->SetLineColor(4); + fBackground->SetMarkerColor(4); + fBackground->SetMarkerStyle(6); + fBackground->DrawCopy("Psame"); + + fSignal->SetMarkerStyle(20); + fSignal->SetMarkerColor(2); + fSignal->DrawCopy("Psame"); + + cSub->cd(1); + fSignal->DrawCopy("P"); + + cSub->cd(2); + fSignPM->DrawCopy("P"); + + cSub->cd(3); + fSignPP->DrawCopy("P"); +} + diff --git a/PWG3/dielectron/AliDielectronSignalExt.h b/PWG3/dielectron/AliDielectronSignalExt.h new file mode 100644 index 00000000000..028bb10668e --- /dev/null +++ b/PWG3/dielectron/AliDielectronSignalExt.h @@ -0,0 +1,77 @@ +#ifndef ALIDIELECTRONSIGNALEXT_H +#define ALIDIELECTRONSIGNALEXT_H + +/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//############################################################# +//# # +//# Class AliDielectronSignalExt # +//# # +//# Authors: # +//# Anton Andronic, GSI / A.Andronic@gsi.de # +//# Ionut C. Arsene, GSI / I.C.Arsene@gsi.de # +//# Julian Book, Uni Ffm / Julian.Book@cern.ch # +//# Frederick Kramer, Uni Ffm, / Frederick.Kramer@cern.ch # +//# Magnus Mager, CERN / Magnus.Mager@cern.ch # +//# WooJin J. Park, GSI / W.J.Park@gsi.de # +//# Jens Wiechula, Uni HD / Jens.Wiechula@cern.ch # +//# # +//############################################################# + +#include +#include +#include + +#include "AliDielectronSignalBase.h" + +class AliDielectronSignalExt : public AliDielectronSignalBase { + +public: + AliDielectronSignalExt(); + AliDielectronSignalExt(const char*name, const char* title); + + virtual ~AliDielectronSignalExt(); + + virtual void Process(TObjArray* const arrhist); + void ProcessLS(TObjArray* const arrhist); // like-sign method + void ProcessEM(TObjArray* const arrhist); // event mixing method + + void SetMethod(Int_t method){ fMethod=method; } + void SetHistograms(TH1F* const unlike, TH1F* const backg, TH1F* const signal); + void SetRebin(Int_t rebin) {fRebin = rebin;} + void SetDrawRange(Double_t min, Double_t max){fDrawMin=min; fDrawMax=max;} + void SetFitRange(Double_t min, Double_t max) {fFitMin=min;fFitMax=max;} + + Double_t GetFitMin() const { return fFitMin; } + Double_t GetFitMax() const { return fFitMax; } + void Rebin(Int_t rebin); + + virtual void Draw(const Option_t* option = ""); + + +private: + + TH1F *fSignPM; // histogram of unlike sign (plus-minus) + TH1F *fSignPP; // histogram of like sign (plus-plus) + TH1F *fSignMM; // histogram of like sign (minus-minus) + TH1F *fBackground; // histogram of like-sign background + TH1F *fSignal; // histogram of subtracted signal + + Int_t fMethod; // subtraction method. 1(like-sign), 2(event mixing) + Int_t fRebin; // number of histogram rebin iteration + Int_t fBins; // number of bins in X axis + Double_t fDrawMin; // minimum X when drawing + Double_t fDrawMax; // maximum X when drawing + Double_t fFitMin; // fit range min + Double_t fFitMax; // fit range max + + AliDielectronSignalExt(const AliDielectronSignalExt &c); + AliDielectronSignalExt &operator=(const AliDielectronSignalExt &c); + + ClassDef(AliDielectronSignalExt,1) // Dielectron SignalFunc +}; + + + +#endif diff --git a/PWG3/dielectron/AliDielectronSignalFunc.cxx b/PWG3/dielectron/AliDielectronSignalFunc.cxx new file mode 100644 index 00000000000..3180cb1787a --- /dev/null +++ b/PWG3/dielectron/AliDielectronSignalFunc.cxx @@ -0,0 +1,349 @@ +/************************************************************************* +* 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. * +**************************************************************************/ + +/////////////////////////////////////////////////////////////////////////// +// Dielectron SignalFunc // +// // +// // +/* +Dielectron signal extraction class using functions as input. + +A function to describe the signal as well as one to describe the background +has to be deployed by the user. Alternatively on of the default implementaions +can be used. + +*/ +// // +/////////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "AliDielectronSignalFunc.h" + +ClassImp(AliDielectronSignalFunc) + +AliDielectronSignalFunc::AliDielectronSignalFunc() : + AliDielectronSignalBase(), + fSignal(0x0), + fBackground(0x0), + fSigBack(0x0), + fFitFunc(0x0), + fVInitParam(0), + fFitOpt("MNQE"), + fUseIntegral(kFALSE), + fFitMin(2.5), + fFitMax(4), + fParM(-1), + fParMres(-1) +{ + // + // Default Constructor + // + +} + +//______________________________________________ +AliDielectronSignalFunc::AliDielectronSignalFunc(const char* name, const char* title) : + AliDielectronSignalBase(name, title), + fSignal(0x0), + fBackground(0x0), + fSigBack(0x0), + fFitFunc(0x0), + fVInitParam(0), + fFitOpt("MNQ"), + fUseIntegral(kFALSE), + fFitMin(2.5), + fFitMax(4), + fParM(-1), + fParMres(-1) +{ + // + // Named Constructor + // +} + +//______________________________________________ +AliDielectronSignalFunc::~AliDielectronSignalFunc() +{ + // + // Default Destructor + // + if (fSigBack) delete fSigBack; +} + + +//______________________________________________ +void AliDielectronSignalFunc::Process(TH1 * const hist) +{ + // + // Fit the mass histogram with the function and retrieve the parameters + // + + //reset result arrays + Reset(); + + //sanity check + if (!fSigBack) { + AliError("Use 'SetFunctions(TF1*,TF1*)' or 'SetDefaults(Int_t)' to setup the signal functions first'!"); + return; + } + + //initial the fit function to its default parameters + if (fVInitParam.GetMatrixArray()) fSigBack->SetParameters(fVInitParam.GetMatrixArray()); + + //Perform fit and retrieve values + hist->Fit(fSigBack,fFitOpt.Data(),"",fFitMin,fFitMax); + + //retrieve values and errors + //TODO: perhpas implement different methods to retrieve the valus + fSignal->SetParameters(fSigBack->GetParameters()); + fBackground->SetParameters(fSigBack->GetParameters()+fSignal->GetNpar()); + + //TODO: proper error estimate?!? + // perhaps this is not possible in a general way? + + // signal and background histograms + TH1 *hSignal=(TH1*)hist->Clone("DieleSignalHist"); + TH1 *hBackground=(TH1*)hist->Clone("DieleBackgroundHist"); + + fSignal->SetRange(hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax()); + fBackground->SetRange(hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax()); + + hSignal->Add(fBackground,-1); + hBackground->Add(fSignal,-1); + + + //get values and errors signal and background + Double_t signal=0, signal_err=0; + Double_t background=0, background_err=0; + Double_t sigPlusBack=0, sigPlusBack_err=0; + + if (!fUseIntegral){ + signal=hSignal->IntegralAndError(hSignal->FindBin(GetIntegralMin()), + hSignal->FindBin(GetIntegralMax()), signal_err); + background=hBackground->IntegralAndError(hBackground->FindBin(GetIntegralMin()), + hBackground->FindBin(GetIntegralMax()), background_err); + sigPlusBack=hist->IntegralAndError(hist->FindBin(GetIntegralMin()), + hist->FindBin(GetIntegralMax()), sigPlusBack_err); + } else { + signal=fSignal->Integral(GetIntegralMin(),GetIntegralMax())/hist->GetBinWidth(hist->FindBin((GetIntegralMax()-GetIntegralMin())/2.)); + background=fBackground->Integral(GetIntegralMin(),GetIntegralMax())/hist->GetBinWidth(hist->FindBin((GetIntegralMax()-GetIntegralMin())/2.)); + sigPlusBack=fSigBack->Integral(GetIntegralMin(),GetIntegralMax())/hist->GetBinWidth(hist->FindBin((GetIntegralMax()-GetIntegralMin())/2.)); + //TODO: Error calculation + } + + //set values + SetSignal(signal,signal_err); + SetBackground(background,background_err); + SetSignificanceAndSOB(); + + //cleanup + delete hSignal; + delete hBackground; +} + +//______________________________________________ +void AliDielectronSignalFunc::Process(TObjArray * const arrhist) +{ + // + // Fit the mass histogram with the function and retrieve the parameters + // + + TH1 *hist=(TH1*)arrhist->At(1); + Process(hist); +} +//______________________________________________ +void AliDielectronSignalFunc::SetFunctions(TF1 * const sig, TF1 * const back, TF1 * const combined, + Int_t parM, Int_t parMres) +{ + // + // Set the signal, background functions and combined fit function + // Note: The process method assumes that the first n parameters in the + // combined fit function correspond to the n parameters of the signal function + // and the n+1 to n+m parameters to the m parameters of the background function!!! + // if parM and parMres are set this information can be used in the drawing function + // to also show in the statistics box the mass and mass resolution + + if (!sig||!back||!combined) { + AliError("Both, signal and background function need to be set!"); + } + fSignal=sig; + fBackground=back; + fSigBack=combined; + + InitParams(); + fParM=parM; + fParMres=parMres; +} + +//______________________________________________ +void AliDielectronSignalFunc::InitParams() +{ + // + // initilise common fit parameters + // + fVInitParam.ResizeTo(fSigBack->GetNpar()); + fVInitParam.SetElements(fSigBack->GetParameters()); +} + +//______________________________________________ +void AliDielectronSignalFunc::SetDefaults(Int_t type) +{ + // + // Setup some default functions: + // type = 0: gaus signal + linear background in 2.5 - 4 GeV inv. mass + // type = 1: gaus signal + exponential background in 2.5 - 4 GeV inv. mass + // type = 2: half gaussian, half exponential signal function + // type = 3: Crystal-Ball function + // type = 4: Crystal-Ball signal + exponential background + // + + + if (type==0){ + fSignal=new TF1("DieleSignal","gaus",2.5,4); + fBackground=new TF1("DieleBackground","pol1",2.5,4); + fSigBack=new TF1("DieleCombined","gaus+pol1(3)",2.5,4); + + fSigBack->SetParameters(1,3.1,.05,2.5,1); + fSigBack->SetParLimits(0,0,10000000); + fSigBack->SetParLimits(1,3.05,3.15); + fSigBack->SetParLimits(2,.02,.1); + + SetFunctions(fSignal,fBackground,fSigBack,1,2); + + } else if (type==1){ + fSignal=new TF1("DieleSignal","gaus",2.5,4); + fBackground=new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])",2.5,4); + fSigBack=new TF1("DieleCombined","gaus+[3]*exp(-(x-[4])/[5])",2.5,4); + + fSigBack->SetParameters(1,3.1,.05,1,2.5,1); + fSigBack->SetParLimits(0,0,10000000); + fSigBack->SetParLimits(1,3.05,3.15); + fSigBack->SetParLimits(2,.02,.1); + + SetFunctions(fSignal,fBackground,fSigBack,1,2); + + } else if (type==2){ + // half gaussian, half exponential signal function + // exponential background + fSignal = new TF1("DieleSignal","(x<[1])*([0]*(exp(-0.5*((x-[1])/[2])^2)+exp((x-[1])/[3])*(1-exp(-0.5*((x-[1])/[2])^2))))+(x>=[1])*([0]*exp(-0.5*((x-[1])/[2])^2))",2.5,4); + fBackground=new TF1("DieleBackground","[0]*exp(-(x-[1])/[2])+[3]",2.5,4); + fSigBack=new TF1("DieleCombined","(x<[1])*([0]*(exp(-0.5*((x-[1])/[2])^2)+exp((x-[1])/[3])*(1-exp(-0.5*((x-[1])/[2])^2))))+(x>=[1])*([0]*exp(-0.5*((x-[1])/[2])^2))+[4]*exp(-(x-[5])/[6])+[7]",2.5,4); + fSigBack->SetParameters(1.,3.1,.05,.1,1,2.5,1,0); + + fSigBack->SetParLimits(0,0,10000000); + fSigBack->SetParLimits(1,3.05,3.15); + fSigBack->SetParLimits(2,.02,.1); + fSigBack->FixParameter(6,2.5); + fSigBack->FixParameter(7,0); + SetFunctions(fSignal,fBackground,fSigBack,1,2); + + } else if (type==3){ + // Crystal-Ball function + RooRealVar m("m","Invariant mass",2.5,3.8); + RooRealVar alpha("alpha","alpha",0.5,0,10); + RooRealVar n("n","n",1,0,10); + RooRealVar m0("m0","m0",3.1,1,5); + RooRealVar sigma("sigma","sigma",0.3,0,3); + RooCBShape jpsi("jpsi","jpsi",m,m0,sigma,alpha,n); + + //RooAddPdf model("model","jpsi fitting function", RooArgList(jpsi,expo), fsig); + RooRealVar fsig("fsig","signal fraction",0.5,0.,1.); + RooAddPdf model("model","jpsi fitting function", RooArgList(jpsi),fsig); + + model.Clone("fFitFunc"); + + } else if (type==4){ + // Crystal-Ball function + RooRealVar m("m","Invariant mass",2.5,3.8); + RooRealVar alpha("alpha","alpha",0.5,0,10); + RooRealVar n("n","n",1,0,10); + RooRealVar m0("m0","m0",3.1,1,5); + RooRealVar sigma("sigma","sigma",0.3,0,3); + RooCBShape jpsi("jpsi","jpsi",m,m0,sigma,alpha,n); + + // Expoenetial function + RooRealVar al("al","al",0.5); + RooExponential expo("expo","exponential", m, al); + + // add two functions + RooRealVar fsig("fsig","signal fraction",0.5,0.,1.); + RooAddPdf model("model","jpsi fitting function", RooArgList(jpsi,expo), fsig); + + model.Clone("fFitFunc"); + } +} + + +//______________________________________________ +void AliDielectronSignalFunc::Draw(const Option_t* option) +{ + // + // Draw the fitted function + // + + TString drawOpt(option); + drawOpt.ToLower(); + + Bool_t optStat=drawOpt.Contains("stat"); + + fSigBack->SetNpx(200); + fSigBack->SetRange(GetIntegralMin(),GetIntegralMax()); + fBackground->SetNpx(200); + fBackground->SetRange(GetIntegralMin(),GetIntegralMax()); + + TGraph *grSig=new TGraph(fSigBack); + grSig->SetFillColor(kGreen); + grSig->SetFillStyle(3001); + + TGraph *grBack=new TGraph(fBackground); + grBack->SetFillColor(kRed); + grBack->SetFillStyle(3001); + + grSig->SetPoint(0,grBack->GetX()[0],grBack->GetY()[0]); + grSig->SetPoint(grSig->GetN()-1,grBack->GetX()[grBack->GetN()-1],grBack->GetY()[grBack->GetN()-1]); + + grBack->SetPoint(0,grBack->GetX()[0],0.); + grBack->SetPoint(grBack->GetN()-1,grBack->GetX()[grBack->GetN()-1],0.); + + fSigBack->SetRange(fFitMin,fFitMax); + fBackground->SetRange(fFitMin,fFitMax); + + if (!drawOpt.Contains("same")){ + grSig->Draw("af"); + } else { + grSig->Draw("f"); + } + grBack->Draw("f"); + fSigBack->Draw("same"); + fBackground->Draw("same"); + + if (optStat) { + TPaveText* pave=DrawStats(); + if (fParM>=0) pave->AddText(Form("Mass: %.3f #pm %.3f", fSigBack->GetParameter(fParM), fSigBack->GetParError(fParM))); + if (fParMres>=0) pave->AddText(Form("Mass res.: %.3f #pm %.3f", fSigBack->GetParameter(fParMres), fSigBack->GetParError(fParMres))); + } +} + diff --git a/PWG3/dielectron/AliDielectronSignalFunc.h b/PWG3/dielectron/AliDielectronSignalFunc.h new file mode 100644 index 00000000000..de1b402c5c7 --- /dev/null +++ b/PWG3/dielectron/AliDielectronSignalFunc.h @@ -0,0 +1,80 @@ +#ifndef ALIDIELECTRONSIGNALFUNC_H +#define ALIDIELECTRONSIGNALFUNC_H + +/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//############################################################# +//# # +//# Class AliDielectronSignalFunc # +//# # +//# Authors: # +//# Anton Andronic, GSI / A.Andronic@gsi.de # +//# Ionut C. Arsene, GSI / I.C.Arsene@gsi.de # +//# Julian Book, Uni Ffm / Julian.Book@cern.ch # +//# Frederick Kramer, Uni Ffm, / Frederick.Kramer@cern.ch # +//# Magnus Mager, CERN / Magnus.Mager@cern.ch # +//# WooJin J. Park, GSI / W.J.Park@gsi.de # +//# Jens Wiechula, Uni HD / Jens.Wiechula@cern.ch # +//# # +//############################################################# + +#include +#include +#include +#include + +#include "AliDielectronSignalBase.h" + +class AliDielectronSignalFunc : public AliDielectronSignalBase { +public: + AliDielectronSignalFunc(); + AliDielectronSignalFunc(const char*name, const char* title); + + virtual ~AliDielectronSignalFunc(); + + void Process(TH1 * const hist); + virtual void Process(TObjArray * const arrhist); + + void SetFunctions(TF1 * const sig, TF1 * const back, TF1 * const combined, Int_t parM=-1, Int_t parMres=-1); + void InitParams(); + void SetFitOption(const char* opt) {fFitOpt=opt;} + void SetDefaults(Int_t type); + void SetFitRange(Double_t min, Double_t max) {fFitMin=min;fFitMax=max;} + + void SetUseIntegral(Bool_t integral=kTRUE) { fUseIntegral=integral; } + + const char* GetFitOption() const { return fFitOpt.Data(); } + TF1* GetSignalFunction() const { return fSignal; } + TF1* GetBackgroundFunction() const { return fBackground; } + TF1* GetCombinedFunction() const { return fSigBack; } + + Double_t GetFitMin() const { return fFitMin; } + Double_t GetFitMax() const { return fFitMax; } + + virtual void Draw(const Option_t* option = ""); + +private: + TF1 *fSignal; // Function for the signal description + TF1 *fBackground; // Function for the background description + TF1 *fSigBack; // Combined function signal plus background + RooAddPdf fFitFunc; // RooFit fit function + TVectorD fVInitParam; // initial fit parameters + TString fFitOpt; // fit option used + + Bool_t fUseIntegral; // use integral instead of bin counts + + Double_t fFitMin; // fit range min + Double_t fFitMax; // fit range max + + Int_t fParM; // Paramter which defines the mass + Int_t fParMres; // Paramter which defines the resolution of the mass + + AliDielectronSignalFunc(const AliDielectronSignalFunc &c); + AliDielectronSignalFunc &operator=(const AliDielectronSignalFunc &c); + + ClassDef(AliDielectronSignalFunc,1) // Dielectron SignalFunc +}; + + +#endif diff --git a/PWG3/dielectron/AliDielectronSpectrum.cxx b/PWG3/dielectron/AliDielectronSpectrum.cxx new file mode 100644 index 00000000000..6d1221a2a36 --- /dev/null +++ b/PWG3/dielectron/AliDielectronSpectrum.cxx @@ -0,0 +1,349 @@ + +/************************************************************************* +* 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. * +**************************************************************************/ + +/////////////////////////////////////////////////////////////////////////// +// Add general description // +// // +// // +/* +Detailed description + + +*/ +// // +/////////////////////////////////////////////////////////////////////////// + +#include +#include +#include + +#include +#include + +#include "AliDielectronSignalBase.h" +#include "AliDielectronSignalFunc.h" + +#include "AliDielectronSpectrum.h" + +ClassImp(AliDielectronSpectrum) + +AliDielectronSpectrum::AliDielectronSpectrum() : + TNamed(), + fCFSignal(0x0), + fCFCorrection(0x0), + fCFSpectrum(0x0), + fCFCorrMatrix(0x0), + fStepSignal(kTRUE), + fStepSignificance(kFALSE), + fStepSOB(kFALSE), + fSignalStep(-1), + fCorrNom(-1), + fCorrDenom(-1), + fSignalMethods(0), + fVariables("Pt"), + fOwnerSpectrum(kTRUE), + fNvars(0), + fVars(0x0), + fNbins(0x0), + fCurrentBins(0x0), + fCurrentPositions(0x0) +{ + // + // Default Constructor + // + fSignalMethods.SetOwner(kFALSE); +} + +//______________________________________________ +AliDielectronSpectrum::AliDielectronSpectrum(const char* name, const char* title) : + TNamed(name, title), + fCFSignal(0x0), + fCFCorrection(0x0), + fCFSpectrum(0x0), + fCFCorrMatrix(0x0), + fStepSignal(kTRUE), + fStepSignificance(kFALSE), + fStepSOB(kFALSE), + fSignalStep(-1), + fCorrNom(-1), + fCorrDenom(-1), + fSignalMethods(0), + fVariables("Pt"), + fOwnerSpectrum(kTRUE), + fNvars(0), + fVars(0x0), + fNbins(0x0), + fCurrentBins(0x0), + fCurrentPositions(0x0) +{ + // + // Named Constructor + // + fSignalMethods.SetOwner(kFALSE); +} + +//______________________________________________ +AliDielectronSpectrum::~AliDielectronSpectrum() +{ + // + // Default Destructor + // + + if (fNbins) delete [] fNbins; + if (fVars) delete [] fVars; + if (fCFCorrMatrix) delete fCFCorrMatrix; + if ( fOwnerSpectrum && fCFSpectrum ) delete fCFSpectrum; +} + +//______________________________________________ +void AliDielectronSpectrum::Process() +{ + // + // Extract signal and perform correction in the specified bins + // + + //sanity checks + if (!fCFSignal){ + AliError("Cannot perform signal extraction, no signal container set"); + return; + } + + if (fSignalMethods.GetEntries()==0){ + AliWarning("No Signal extraction method specified, using a default one"); + AliDielectronSignalFunc *func=new AliDielectronSignalFunc("gaus+exp","Gauss for Signal and Exponential background"); + func->SetDefaults(1); + fSignalMethods.Add(func); + fSignalMethods.SetOwner(); + } + + //setup configured variables + if (!SetupVariables()) return; + + //create container for the spectrum + CreateCFSpectrum(); + + if (!fCFSpectrum){ + AliError("Could not create the Spectrum container"); + return; + } + + //get efficiency map if correction container is available + if (fCFCorrection&&!fCFCorrMatrix){ + CreateCorrectionMatrix(); + } + + //loop over all configured bins and extract the signal + fCurrentBins=new Int_t[fNvars]; + fCurrentPositions=new Double_t[fNvars]; + ExtractSignalInBins(); + delete [] fCurrentBins; + delete [] fCurrentPositions; + if (fSignalMethods.IsOwner()) { + fSignalMethods.Delete(); + fSignalMethods.SetOwner(kFALSE); + } + +} + +//______________________________________________ +Bool_t AliDielectronSpectrum::SetupVariables() +{ + // + // Setup the variables arrays + // + + TObjArray *arr=fVariables.Tokenize(":"); + fNvars=arr->GetEntries(); + fVars=new Int_t[fNvars]; + fNbins=new Int_t[fNvars]; + + for (Int_t iVar=0; iVarGetVar(arr->UncheckedAt(iVar)->GetName()); + if (fVars[iVar]==-1){ + AliError(Form("Variable '%s' not found in Signal container!",arr->UncheckedAt(iVar)->GetName())); + delete [] fVars; + fVars=0x0; + delete [] fNbins; + fNbins=0x0; + delete arr; + return kFALSE; + } + + fNbins[iVar]=fCFSignal->GetNBins(fVars[iVar]); + } + delete arr; + return kTRUE; +} + +//______________________________________________ +void AliDielectronSpectrum::CreateCFSpectrum() +{ + // + // Create CF container for the spectrum + // + + Int_t nAddStep=0; + if (fStepSignal) nAddStep+=2; + if (fStepSignificance) ++nAddStep; + if (fStepSOB) ++nAddStep; + + Int_t nStep=nAddStep*(fSignalMethods.GetEntries()); + if (fSignalMethods.GetEntries()>1) nStep+=nAddStep; + + fCFSpectrum = new AliCFContainer(GetName(), GetTitle(), nStep, fNvars, fNbins); + + // initialize the variables and their bin limits + for (Int_t iVar=0; iVarSetBinLimits(iVar, fCFSignal->GetBinLimits(fVars[iVar])); + fCFSpectrum->SetVarTitle(iVar, fCFSignal->GetVarTitle(fVars[iVar])); + } + + // setup step titles + Int_t steps=0; + for (Int_t iMethod=0; iMethodGetName()); + if (fStepSignal){ + fCFSpectrum->SetStepTitle(steps++,(name+" (Siganl)").Data()); + fCFSpectrum->SetStepTitle(steps++,(name+" (Corrected Siganl)").Data()); + } + if (fStepSignificance){ + fCFSpectrum->SetStepTitle(steps++,(name+" (Significance)").Data()); + } + if (fStepSOB){ + fCFSpectrum->SetStepTitle(steps++,(name+" (S/B)").Data()); + } + } + + if (fSignalMethods.GetEntries()>1){ + fCFSpectrum->SetStepTitle(steps++,"Mean of methods"); + fCFSpectrum->SetStepTitle(steps++,"Mean of methods (Corrected)"); + } + + if (nStep!=steps){ + AliError("Something went wrong in the step creation"); + delete fCFSpectrum; + fCFSpectrum=0x0; + } +} + +//______________________________________________ +void AliDielectronSpectrum::CreateCorrectionMatrix() +{ + // + // Get the correction matrix for the corresponding variables + // + + if (!fCFCorrection) return; + + TObjArray *arr=fVariables.Tokenize(":"); + Int_t nvars=arr->GetEntries(); + Int_t *vars=new Int_t[nvars]; + + for (Int_t iVar=0; iVarGetVar(arr->UncheckedAt(iVar)->GetName()); + if (vars[iVar]==-1){ + AliError(Form("Variable '%s' not found in Correction container!",arr->UncheckedAt(iVar)->GetName())); + delete [] vars; + delete arr; + return; + } + } + delete arr; + + fCFCorrMatrix =fCFCorrection->GetGrid(fCorrNom)->Project(nvars,vars,0,0); + AliCFGridSparse *hnDeNom=fCFCorrection->GetGrid(fCorrDenom)->Project(nvars,vars,0,0); + fCFCorrMatrix->Divide(hnDeNom); + delete hnDeNom; +} + +//______________________________________________ +void AliDielectronSpectrum::ExtractSignalInBins(Int_t variable) +{ + // + // recursively loop over bins and extract signal + // + + Int_t varPairType=fCFSignal->GetVar("PairType"); + Int_t varMass=fCFSignal->GetVar("M"); + + for (Int_t ibin=0; ibinGetGrid(fSignalStep)->GetGrid()->GetAxis(fVars[variable])->SetRange(bin,bin); + fCurrentBins[variable]=bin; + fCurrentPositions[variable]=fCFSignal->GetGrid(fSignalStep)->GetBinCenter(fVars[variable],bin); + + if (variable != fNvars-1) ExtractSignalInBins(variable+1); + + TObjArray arrPairTypes(10); + arrPairTypes.SetOwner(); + + for (Int_t itype=0; itype<3; ++itype){ +// Int_t itype=1; + fCFSignal->SetRangeUser(varPairType,itype,itype,fSignalStep); + TH1 *h=fCFSignal->Project(varMass,fSignalStep); + h->SetDirectory(0); + arrPairTypes.AddAt(h,itype); + } + AliInfo(Form("Processing bin: %d (%.2f)",ibin, fCurrentPositions[variable])); + //loop over all signal extraction methods and retrieve signals + for (Int_t iMethod=0; iMethodProcess(&arrPairTypes); + const TVectorD &val=signalMethod->GetValues(); + const TVectorD &err=signalMethod->GetErrors(); + + //Fill sparse containers + Int_t step=0; + if (fStepSignal){ + //signal + fCFSpectrum->GetGrid(step)->SetElement(fCurrentBins,val(0)); + fCFSpectrum->GetGrid(step)->SetElementError(fCurrentBins,err(0)); + ++step; + + //corrected signal + if (fCFCorrMatrix){ + Float_t corrFactor = fCFCorrMatrix->GetElement(fCurrentPositions); + Float_t corrError = fCFCorrMatrix->GetElementError(fCurrentPositions); + + Float_t value=val(0)*corrFactor; + fCFSpectrum->GetGrid(step)->SetElement(fCurrentBins,value); + Float_t error=TMath::Sqrt( (err(0)/val(0))*(err(0)/val(0)) + + (corrError/corrFactor)*(corrError/corrFactor) + )*value; + fCFSpectrum->GetGrid(step)->SetElementError(fCurrentBins,error); +// printf("corrFactor: %f+-%f\n",value,error); + } + ++step; + } + + if (fStepSignificance) { + fCFSpectrum->GetGrid(step)->SetElement(fCurrentBins,val(2)); + fCFSpectrum->GetGrid(step)->SetElementError(fCurrentBins,err(2)); + ++step; + } + + if (fStepSOB) { + fCFSpectrum->GetGrid(step)->SetElement(fCurrentBins,val(3)); + fCFSpectrum->GetGrid(step)->SetElementError(fCurrentBins,err(3)); + ++step; + } + }// end method loop + + arrPairTypes.Delete(); + + }// end bin loop +} + diff --git a/PWG3/dielectron/AliDielectronSpectrum.h b/PWG3/dielectron/AliDielectronSpectrum.h new file mode 100644 index 00000000000..938ca9b80cb --- /dev/null +++ b/PWG3/dielectron/AliDielectronSpectrum.h @@ -0,0 +1,114 @@ +#ifndef ALIDIELECTRONSPECTRUM_H +#define ALIDIELECTRONSPECTRUM_H + +/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//############################################################# +//# # +//# Class AliDielectronSpectrum # +//# Manage Cuts on the legs of the pair # +//# # +//# Authors: # +//# Anton Andronic, GSI / A.Andronic@gsi.de # +//# Ionut C. Arsene, GSI / I.C.Arsene@gsi.de # +//# Julian Book, Uni Ffm / Julian.Book@cern.ch # +//# Frederick Kramer, Uni Ffm, / Frederick.Kramer@cern.ch # +//# Magnus Mager, CERN / Magnus.Mager@cern.ch # +//# WooJin J. Park, GSI / W.J.Park@gsi.de # +//# Jens Wiechula, Uni HD / Jens.Wiechula@cern.ch # +//# # +//############################################################# + +#include +#include + +#include + +class AliCFGridSparse; +class AliDielectronSignalBase; +class AliCFGridSparse; +class AliCFContainer; + +class AliDielectronSpectrum : public TNamed { +public: + AliDielectronSpectrum(); + AliDielectronSpectrum(const char*name, const char* title); + + virtual ~AliDielectronSpectrum(); + + void AddMethod(AliDielectronSignalBase * const method) {fSignalMethods.Add(method);} + + void SetCorrectionContainer(AliCFContainer * const container, Int_t nominator, Int_t denominator); + void SetSignalContainer(AliCFContainer * const container, Int_t step); + + void SetVariables(const char* vars) { fVariables=vars; } + + void SetStepForSignal(Bool_t step=kTRUE) { fStepSignal=step; } + void SetStepForSignificance(Bool_t step=kTRUE) { fStepSignificance=step; } + void SetStepForSignalOverBackground(Bool_t step=kTRUE) { fStepSOB=step; } + + void SetNoOwnerSpectrum(Bool_t noOwner=kTRUE) { fOwnerSpectrum=!noOwner; } + + AliCFContainer* GetSpectrumContainer() const { return fCFSpectrum; } + AliCFGridSparse* GetCorrectionMatrix() const { return fCFCorrMatrix; } + + void Process(); + + +private: + AliCFContainer *fCFSignal; // CF container with from which to extract the Signal + AliCFContainer *fCFCorrection; // CF container from which to extract the correction matrix + AliCFContainer *fCFSpectrum; // CF container with extracted signal + AliCFGridSparse *fCFCorrMatrix; // correction matrix + + Bool_t fStepSignal; // if to create a step for the signal + Bool_t fStepSignificance; // if to create a step for the significance + Bool_t fStepSOB; // if to create a step for signal over background + + Int_t fSignalStep; // step to use from the signal container + + Int_t fCorrNom; // Nominator to use from corr matrix container + Int_t fCorrDenom; // Deominator to use from corr matrix container + + TObjArray fSignalMethods; // array with signal extraction methods + TString fVariables; // variable names as a function of which to extract the signal + + Bool_t fOwnerSpectrum; // if we own the creted spectrum + + Int_t fNvars; //! number of variables + Int_t *fVars; //! variable numbers translated from fVariables + Int_t *fNbins; //! number of bins for each variable + Int_t *fCurrentBins; //! bin currently selected for each variable + Double_t *fCurrentPositions; //! variables values currently selected + + void Fill(Int_t *bin, Int_t step, Double_t value, Double_t error); + Bool_t SetupVariables(); + void CreateCorrectionMatrix(); + void CreateCFSpectrum(); + void ExtractSignalInBins(Int_t variable=0); + + AliDielectronSpectrum(const AliDielectronSpectrum &c); + AliDielectronSpectrum &operator=(const AliDielectronSpectrum &c); + + ClassDef(AliDielectronSpectrum,1) //Cut class providing cuts for both legs of a pair + +}; + +// +// inline functions +// +inline void AliDielectronSpectrum::SetCorrectionContainer(AliCFContainer * const container, Int_t nominator, Int_t denominator) +{ + fCFCorrection=container; + fCorrNom=nominator; + fCorrDenom=denominator; +} + +inline void AliDielectronSpectrum::SetSignalContainer(AliCFContainer * const container, Int_t step) +{ + fCFSignal=container; + fSignalStep=step; +} + +#endif diff --git a/PWG3/dielectron/AliDielectronVarCuts.cxx b/PWG3/dielectron/AliDielectronVarCuts.cxx index 43161841698..d1ebc865067 100644 --- a/PWG3/dielectron/AliDielectronVarCuts.cxx +++ b/PWG3/dielectron/AliDielectronVarCuts.cxx @@ -49,6 +49,7 @@ AliDielectronVarCuts::AliDielectronVarCuts() : fActiveCuts[i]=0; fCutMin[i]=0; fCutMax[i]=0; + fCutExclude[i]=kFALSE; } } @@ -106,7 +107,7 @@ Bool_t AliDielectronVarCuts::IsSelected(TObject* track) for (Int_t iCut=0; iCutfCutMax[iCut]) ) CLRBIT(fSelectedCutsMask,iCut); + if ( ((values[cut]fCutMax[iCut]))^fCutExclude[iCut] ) CLRBIT(fSelectedCutsMask,iCut); } Bool_t isSelected=(fSelectedCutsMask==fActiveCutsMask); @@ -116,7 +117,7 @@ Bool_t AliDielectronVarCuts::IsSelected(TObject* track) } //________________________________________________________________________ -void AliDielectronVarCuts::AddCut(AliDielectronVarManager::ValueTypes type, Double_t min, Double_t max) +void AliDielectronVarCuts::AddCut(AliDielectronVarManager::ValueTypes type, Double_t min, Double_t max, Bool_t excludeRange) { // // Set cut range and activate it @@ -128,6 +129,7 @@ void AliDielectronVarCuts::AddCut(AliDielectronVarManager::ValueTypes type, Doub } fCutMin[fNActiveCuts]=min; fCutMax[fNActiveCuts]=max; + fCutExclude[fNActiveCuts]=excludeRange; SETBIT(fActiveCutsMask,fNActiveCuts); fActiveCuts[fNActiveCuts]=(UShort_t)type; ++fNActiveCuts; @@ -147,7 +149,14 @@ void AliDielectronVarCuts::Print(const Option_t* /*option*/) const } for (Int_t iCut=0; iCut + #include -#include +#include "AliDielectronVarManager.h" class AliDielectronVarCuts : public AliAnalysisCuts { public: @@ -34,8 +35,8 @@ public: AliDielectronVarCuts(const char* name, const char* title); virtual ~AliDielectronVarCuts(); //TODO: make copy constructor and assignment operator public - void AddCut(AliDielectronVarManager::ValueTypes type, Double_t min, Double_t max); - void AddCut(AliDielectronVarManager::ValueTypes type, Double_t value); + void AddCut(AliDielectronVarManager::ValueTypes type, Double_t min, Double_t max, Bool_t excludeRange=kFALSE); + void AddCut(AliDielectronVarManager::ValueTypes type, Double_t value, Bool_t excludeRange=kFALSE); // setters void SetCutOnMCtruth(Bool_t mc=kTRUE) { fCutOnMCtruth=mc; } @@ -76,6 +77,7 @@ private: Double_t fCutMin[AliDielectronVarManager::kNMaxValues]; // minimum values for the cuts Double_t fCutMax[AliDielectronVarManager::kNMaxValues]; // maximum values for the cuts + Bool_t fCutExclude[AliDielectronVarManager::kNMaxValues]; // inverse cut logic? AliDielectronVarCuts(const AliDielectronVarCuts &c); AliDielectronVarCuts &operator=(const AliDielectronVarCuts &c); @@ -87,13 +89,13 @@ private: // //Inline functions // -inline void AliDielectronVarCuts::AddCut(AliDielectronVarManager::ValueTypes type, Double_t value) +inline void AliDielectronVarCuts::AddCut(AliDielectronVarManager::ValueTypes type, Double_t value, Bool_t excludeRange) { // // Set cut in a small delta around value // const Double_t kDelta=1e-20; - AddCut(type,value-kDelta,value+kDelta); + AddCut(type,value-kDelta,value+kDelta, excludeRange); } #endif diff --git a/PWG3/dielectron/AliDielectronVarManager.cxx b/PWG3/dielectron/AliDielectronVarManager.cxx index e6dc32cf288..6f7c178fa79 100644 --- a/PWG3/dielectron/AliDielectronVarManager.cxx +++ b/PWG3/dielectron/AliDielectronVarManager.cxx @@ -57,11 +57,21 @@ const char* AliDielectronVarManager::fgkParticleNames[AliDielectronVarManager::k "P_InnerParam", "TPC_signal", "TPC_nSigma_Electrons", + "TPC_nSigma_Pions", + "TPC_nSigma_Muons", + "TPC_nSigma_Kaons", + "TPC_nSigma_Protons", + "TOF_nSigma_Pions", + "TOF_nSigma_Muons", + "TOF_nSigma_Kaons", + "TOF_nSigma_Protons", // "Chi2NDF", "DecayLength", "R", "OpeningAngle", + "LegDistance", + "LegDistanceXY", "Merr", "DCA", "PairType", @@ -77,8 +87,9 @@ const char* AliDielectronVarManager::fgkParticleNames[AliDielectronVarManager::k "Nevents" }; -AliESDpid* AliDielectronVarManager::fgESDpid = new AliESDpid; +AliESDpid* AliDielectronVarManager::fgESDpid = 0x0; AliVEvent* AliDielectronVarManager::fgEvent = 0x0; +AliKFVertex* AliDielectronVarManager::fgKFVertex = 0x0; //________________________________________________________________ AliDielectronVarManager::AliDielectronVarManager() : TNamed("AliDielectronVarManager","AliDielectronVarManager") diff --git a/PWG3/dielectron/AliDielectronVarManager.h b/PWG3/dielectron/AliDielectronVarManager.h index 7e82f2c9b13..3a1eeb3b484 100644 --- a/PWG3/dielectron/AliDielectronVarManager.h +++ b/PWG3/dielectron/AliDielectronVarManager.h @@ -12,6 +12,7 @@ //# Anton Andronic, GSI / A.Andronic@gsi.de # //# Ionut C. Arsene, GSI / I.C.Arsene@gsi.de # //# Julian Book, Uni Ffm / Julian.Book@cern.ch # +//# Markus Köhler, GSI / M.Koehler@gsi.de # //# Frederick Kramer, Uni Ffm / Frederick.Kramer@cern.ch # //# Magnus Mager, CERN / Magnus.Mager@cern.ch # //# WooJin J. Park, GSI / W.J.Park@gsi.de # @@ -32,7 +33,9 @@ #include #include #include +#include #include +#include #include #include #include // ? @@ -80,7 +83,18 @@ public: kPdgCode, // PDG code kPIn, // momentum at inner wall of TPC (if available), used for PID kTPCsignal, // TPC dE/dx signal + kTPCnSigmaEle, // number of sigmas to the dE/dx electron line in the TPC + kTPCnSigmaPio, // number of sigmas to the dE/dx pion line in the TPC + kTPCnSigmaMuo, // number of sigmas to the dE/dx muon line in the TPC + kTPCnSigmaKao, // number of sigmas to the dE/dx kaon line in the TPC + kTPCnSigmaPro, // number of sigmas to the dE/dx proton line in the TPC + + kTOFnSigmaPio, // number of sigmas to the pion line in the TOF + kTOFnSigmaMuo, // number of sigmas to the muon line in the TOF + kTOFnSigmaKao, // number of sigmas to the kaon line in the TOF + kTOFnSigmaPro, // number of sigmas to the proton line in the TOF + kParticleMax, // // TODO: kRNClusters ?? // AliDielectronPair specific variables @@ -88,6 +102,8 @@ public: kDecayLength, // decay length kR, // distance to the origin kOpeningAngle, // opening angle + kLegDist, // distance of the legs + kLegDistXY, // distance of the legs in XY kMerr, // error of mass calculation kDCA, // distance of closest approach TODO: not implemented yet kPairType, // type of the pair, like like sign ++ unlikesign ... @@ -113,7 +129,9 @@ public: static void Fill(const TObject* particle, Double_t * const values); static void InitESDpid(Int_t type=0); - static void SetEvent(AliVEvent * const ev) { fgEvent = ev; } + static void SetESDpid(AliESDpid * const pid) {fgESDpid=pid;} + static AliESDpid* GetESDpid() {return fgESDpid;} + static void SetEvent(AliVEvent * const ev); static const char* GetValueName(Int_t i) { return (i>=0&&iNumberOfSigmasTPC(particle,AliPID::kElectron); + values[AliDielectronVarManager::kTPCnSigmaPio]=fgESDpid->NumberOfSigmasTPC(particle,AliPID::kPion); + values[AliDielectronVarManager::kTPCnSigmaMuo]=fgESDpid->NumberOfSigmasTPC(particle,AliPID::kMuon); + values[AliDielectronVarManager::kTPCnSigmaKao]=fgESDpid->NumberOfSigmasTPC(particle,AliPID::kKaon); + values[AliDielectronVarManager::kTPCnSigmaPro]=fgESDpid->NumberOfSigmasTPC(particle,AliPID::kProton); + + Double_t t0=fgESDpid->GetTOFResponse().GetTimeZero(); + values[AliDielectronVarManager::kTOFnSigmaPio]=fgESDpid->NumberOfSigmasTOF(particle,AliPID::kPion,t0); + values[AliDielectronVarManager::kTOFnSigmaMuo]=fgESDpid->NumberOfSigmasTOF(particle,AliPID::kMuon,t0); + values[AliDielectronVarManager::kTOFnSigmaKao]=fgESDpid->NumberOfSigmasTOF(particle,AliPID::kKaon,t0); + values[AliDielectronVarManager::kTOFnSigmaPro]=fgESDpid->NumberOfSigmasTOF(particle,AliPID::kProton,t0); } inline void AliDielectronVarManager::FillVarAODTrack(const AliAODTrack *particle, Double_t * const values) @@ -235,9 +264,66 @@ inline void AliDielectronVarManager::FillVarAODTrack(const AliAODTrack *particle // Fill common AliVParticle interface information FillVarVParticle(particle, values); - // Fill AliAODTrack interface information + + // Reset AliESDtrack interface specific information + values[AliDielectronVarManager::kNclsITS] = 0; + values[AliDielectronVarManager::kNclsTPC] = 0; + values[AliDielectronVarManager::kNFclsTPC] = 0; + values[AliDielectronVarManager::kNclsTRD] = 0; + values[AliDielectronVarManager::kTRDntracklets] = 0; + values[AliDielectronVarManager::kTRDpidQuality] = 0; + + //TODO: This is only an approximation!!! + values[AliDielectronVarManager::kTPCsignalN] = particle->GetTPCClusterMap().CountBits(); + +// Fill AliAODTrack interface information // ... + values[AliDielectronVarManager::kImpactParXY] = particle->DCA(); + values[AliDielectronVarManager::kImpactParZ] = particle->ZAtDCA(); + + values[AliDielectronVarManager::kPIn]=0; + values[AliDielectronVarManager::kTPCsignal]=0; + + values[AliDielectronVarManager::kTPCnSigmaEle]=0; + values[AliDielectronVarManager::kTPCnSigmaPio]=0; + values[AliDielectronVarManager::kTPCnSigmaMuo]=0; + values[AliDielectronVarManager::kTPCnSigmaKao]=0; + values[AliDielectronVarManager::kTPCnSigmaPro]=0; + + AliAODPid *pid=particle->GetDetPid(); + if (pid){ + Double_t mom =pid->GetTPCmomentum(); + //TODO: kTPCsignalN is only an approximation (see above)!! + + Double_t tpcNsigmaEle=fgESDpid->GetTPCResponse().GetNumberOfSigmas(mom,pid->GetTPCsignal(), + TMath::Nint(values[AliDielectronVarManager::kTPCsignalN]) ,AliPID::kElectron); + + Double_t tpcNsigmaPio=fgESDpid->GetTPCResponse().GetNumberOfSigmas(mom,pid->GetTPCsignal(), + TMath::Nint(values[AliDielectronVarManager::kTPCsignalN]),AliPID::kPion); + + Double_t tpcNsigmaMuo=fgESDpid->GetTPCResponse().GetNumberOfSigmas(mom,pid->GetTPCsignal(), + TMath::Nint(values[AliDielectronVarManager::kTPCsignalN]),AliPID::kMuon); + + Double_t tpcNsigmaKao=fgESDpid->GetTPCResponse().GetNumberOfSigmas(mom,pid->GetTPCsignal(), + TMath::Nint(values[AliDielectronVarManager::kTPCsignalN]),AliPID::kKaon); + + Double_t tpcNsigmaPro=fgESDpid->GetTPCResponse().GetNumberOfSigmas(mom,pid->GetTPCsignal(), + TMath::Nint(values[AliDielectronVarManager::kTPCsignalN]),AliPID::kProton); + + values[AliDielectronVarManager::kPIn]=pid->GetTPCmomentum(); + values[AliDielectronVarManager::kTPCsignal]=pid->GetTPCsignal(); + + values[AliDielectronVarManager::kTPCnSigmaEle]=tpcNsigmaEle; + values[AliDielectronVarManager::kTPCnSigmaPio]=tpcNsigmaPio; + values[AliDielectronVarManager::kTPCnSigmaMuo]=tpcNsigmaMuo; + values[AliDielectronVarManager::kTPCnSigmaKao]=tpcNsigmaKao; + values[AliDielectronVarManager::kTPCnSigmaPro]=tpcNsigmaPro; + + values[AliDielectronVarManager::kTRDntracklets] = 0; + values[AliDielectronVarManager::kTRDpidQuality] = 0; + + } } inline void AliDielectronVarManager::FillVarMCParticle(const AliMCParticle *particle, Double_t * const values) @@ -282,7 +368,9 @@ inline void AliDielectronVarManager::FillVarDielectronPair(const AliDielectronPa values[AliDielectronVarManager::kDecayLength] = kfPair.GetDecayLength(); values[AliDielectronVarManager::kR] = kfPair.GetR(); values[AliDielectronVarManager::kOpeningAngle] = pair->OpeningAngle(); - values[AliDielectronVarManager::kMerr] = kfPair.GetErrMass()>0?kfPair.GetErrMass()/kfPair.GetMass():1000000; + values[AliDielectronVarManager::kLegDist] = pair->DistanceDaughters(); + values[AliDielectronVarManager::kLegDistXY] = pair->DistanceDaughtersXY(); + values[AliDielectronVarManager::kMerr] = kfPair.GetErrMass()>1e-30&&kfPair.GetErrMass()>1e-30?kfPair.GetErrMass()/kfPair.GetMass():1000000; values[AliDielectronVarManager::kPairType] = pair->GetType(); } @@ -349,6 +437,7 @@ inline void AliDielectronVarManager::InitESDpid(Int_t type) // type=0 is simulation // type=1 is data + if (!fgESDpid) fgESDpid=new AliESDpid; Double_t alephParameters[5]; // simulation alephParameters[0] = 2.15898e+00/50.; @@ -372,6 +461,15 @@ inline void AliDielectronVarManager::InitESDpid(Int_t type) } + +inline void AliDielectronVarManager::SetEvent(AliVEvent * const ev) +{ + + fgEvent = ev; + if (fgKFVertex) delete fgKFVertex; + fgKFVertex=0x0; + if (ev) fgKFVertex=new AliKFVertex(*ev->GetPrimaryVertex()); +} /* inline void AliDielectronVarManager::FillValues(const TParticle *particle, Double_t *values) { diff --git a/PWG3/dielectron/macros/ConfigJpsi2eeData.C b/PWG3/dielectron/macros/ConfigJpsi2eeData.C new file mode 100644 index 00000000000..35a5a83dba6 --- /dev/null +++ b/PWG3/dielectron/macros/ConfigJpsi2eeData.C @@ -0,0 +1,221 @@ + +void InitHistograms(AliDielectron *die, Int_t cutDefinition); +void InitCF(AliDielectron* die, Int_t cutDefinition); + +void SetupTrackCuts(AliDielectron *die, Int_t cutDefinition); +void SetupPairCuts(AliDielectron *die, Int_t cutDefinition); + +AliESDtrackCuts *SetupESDtrackCuts(Int_t cutDefinition); + +TString names=("trackQ+Pt>1GeV+loose_PID;trackQ+Pt>1GeV+tight_PID;trackQ+Pt>.5GeV"); +TObjArray *arrNames=names.Tokenize(";"); + +const Int_t nDie=arrNames->GetEntries(); + +AliDielectron* ConfigJpsi2ee(Int_t cutDefinition) +{ + // + // Setup the instance of AliDielectron + // + + // create the actual framework object + TString name=Form("%02d",cutDefinition); + if (cutDefinitionGetEntriesFast()){ + name=arrNames->At(cutDefinition)->GetName(); + } + AliDielectron *die = + new AliDielectron(Form("%s",name.Data()), + Form("Track cuts: %s",name.Data())); + + // cut setup + SetupTrackCuts(die,cutDefinition); + SetupPairCuts(die,cutDefinition); + + // + // histogram setup + // only if an AliDielectronHistos object is attached to the + // dielectron framework the QA histograms will be filled + // + InitHistograms(die,cutDefinition); + + // the last definition uses no cuts and only the QA histograms should be filled! + if (cutDefinitionGetTrackFilter().AddCuts(SetupESDtrackCuts(cutDefinition)); + + + //QA no CF + if (cutDefinition==nDie-1) { + //Pt cut + AliDielectronVarCuts *pt = new AliDielectronVarCuts("Pt>.5","Pt>.5"); + pt->AddCut(AliDielectronVarManager::kPt,.5,1e30); + die->GetTrackFilter().AddCuts(pt); + + return; + } + + + AliDielectronVarCuts *pt = new AliDielectronVarCuts("Pt>1","Pt>1"); + pt->AddCut(AliDielectronVarManager::kPt,1.,1e30); + die->GetTrackFilter().AddCuts(pt); + + //loose PID + if (cutDefinition==0){ + //ESD pid cuts (TPC nSigma Electrons) + AliDielectronVarCuts *pid = new AliDielectronVarCuts("TPCpid","TPC nSigma cut |e|<3 pi>2 |P|>2 |K|>2"); + //include + pid->AddCut(AliDielectronVarManager::kTPCnSigmaEle, -3., 3.); + //exclude + pid->AddCut(AliDielectronVarManager::kTPCnSigmaPio, -20., 2., kTRUE); + pid->AddCut(AliDielectronVarManager::kTPCnSigmaPro, -2., 2., kTRUE); + pid->AddCut(AliDielectronVarManager::kTPCnSigmaKao, -2., 2., kTRUE); + die->GetTrackFilter().AddCuts(pid); + } + + //tight PID + if (cutDefinition==1){ + //ESD pid cuts (TPC nSigma Electrons) + AliDielectronVarCuts *pid = new AliDielectronVarCuts("TPCpid","TPC nSigma cut |e|<2 pi>2 |P|>2 |K|>2; |dXY|<200um"); + //include + pid->AddCut(AliDielectronVarManager::kTPCnSigmaEle, -2., 2.); + //exclude + pid->AddCut(AliDielectronVarManager::kTPCnSigmaPio, -2., 2., kTRUE); + pid->AddCut(AliDielectronVarManager::kTPCnSigmaPro, -2., 2., kTRUE); + pid->AddCut(AliDielectronVarManager::kTPCnSigmaKao, -2., 2., kTRUE); + //tighten Impact parameter to 200um + pid->AddCut(AliDielectronVarManager::kImpactParXY, -0.02, 0.02); + die->GetTrackFilter().AddCuts(pid); + } + +} + +//______________________________________________________________________________________ +void SetupPairCuts(AliDielectron *die, Int_t cutDefinition) +{ + // + // Setup the pair cuts + // + + + // reject conversions + // and select mass region + AliDielectronVarCuts *openingAngleCut=new AliDielectronVarCuts("OpeningAngle","Opening angle > 35mrad"); + openingAngleCut->AddCut(AliDielectronVarManager::kOpeningAngle,.035,4.); + openingAngleCut->AddCut(AliDielectronVarManager::kM,2.,4.); + die->GetPairFilter().AddCuts(openingAngleCut); +} + +//______________________________________________________________________________________ +AliESDtrackCuts *SetupESDtrackCuts(Int_t cutDefinition) +{ + // + // Setup default AliESDtrackCuts + // + AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts; + + esdTrackCuts->SetMaxDCAToVertexZ(3.0); + esdTrackCuts->SetMaxDCAToVertexXY(.07); + esdTrackCuts->SetRequireTPCRefit(kTRUE); + esdTrackCuts->SetRequireITSRefit(kTRUE); + esdTrackCuts->SetAcceptKinkDaughters(kFALSE); + esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst); + + esdTrackCuts->SetMinNClustersTPC(75); + esdTrackCuts->SetMaxChi2PerClusterTPC(4); + + return esdTrackCuts; +} + +//______________________________________________________________________________________ +void InitHistograms(AliDielectron *die, Int_t cutDefinition) +{ + // + // Initialise the QA histograms + // + + //Setup QA histograms + AliDielectronHistos *histos= + new AliDielectronHistos(die->GetName(), + die->GetTitle()); + + //Initialise histogram classes + histos->SetReservedWords("Track;Pair"); + + //Event class (only for last QA) + if (cutDefinition==nDie-1) histos->AddClass("Event"); + + //Track classes, only first event + for (Int_t i=0; i<2; ++i){ + histos->AddClass(Form("Track_%s",AliDielectron::TrackClassName(i))); + } + + //Pair classes, only first event + for (Int_t i=0; i<3; ++i){ + histos->AddClass(Form("Pair_%s",AliDielectron::PairClassName(i))); + } + + //Event histograms + if (cutDefinition==nDie-1){ + //add histograms to event class + histos->UserHistogram("Event","nEvents","Number of processed events after cuts;Number events", + 1,0.,1.,AliDielectronVarManager::kNevents); + } + + //add histograms to Track classes + histos->UserHistogram("Track","Pt","Pt;Pt [GeV];#tracks",200,0,20.,AliDielectronVarManager::kPt); + histos->UserHistogram("Track","dEdx_P","dEdx;P [GeV];TPC signal (arb units);#tracks", + 400,0.2,20.,200,0.,200.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCsignal,kTRUE); + histos->UserHistogram("Track","TPCnSigmaEle_P","TPC number of sigmas Electrons;P [GeV];TPC number of sigmas Electrons;#tracks", + 400,0.2,20.,100,-5.,5.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaEle,kTRUE); + histos->UserHistogram("Track","TPCnSigmaPio_P","TPC number of sigmas Pions;P [GeV];TPC number of sigmas Pions;#tracks", + 400,0.2,20.,100,-5.,5.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaPio,kTRUE); + histos->UserHistogram("Track","dXY","dXY;dXY [cm];#tracks",200,-2.,2.,AliDielectronVarManager::kImpactParXY); + + //add histograms to Pair classes + histos->UserHistogram("Pair","InvMass","Inv.Mass;Inv. Mass [GeV];#pairs", + 500,0.,4.,AliDielectronVarManager::kM); + histos->UserHistogram("Pair","Rapidity","Rapidity;Rapidity;#pairs", + 100,-2.,2.,AliDielectronVarManager::kY); + histos->UserHistogram("Pair","OpeningAngle","Opening angle;angle", + 100,0.,3.15,AliDielectronVarManager::kOpeningAngle); + histos->UserHistogram("Pair","Chi2/NDF","#Chi^{2}/NDF;#Chi^{2}/NDF", + 100, 0., 20., AliDielectronVarManager::kChi2NDF); + + die->SetHistogramManager(histos); +} + + +void InitCF(AliDielectron* die, Int_t cutDefinition) +{ + // + // Setupd the CF Manager if needed + // + AliDielectronCF *cf=new AliDielectronCF(die->GetName(),die->GetTitle()); + + //pair variables + cf->AddVariable(AliDielectronVarManager::kPt,20,0,10); + cf->AddVariable(AliDielectronVarManager::kY,40,-2,2); + cf->AddVariable(AliDielectronVarManager::kM,200,-0.01,3.99); + cf->AddVariable(AliDielectronVarManager::kPairType,10,0,10); + //leg variables + cf->AddVariable(AliDielectronVarManager::kPt,20,0,10,kTRUE); + + //only in this case write MC truth info + cf->SetStepsForCutsIncreasing(); + if (cutDefinition==0){ + cf->SetStepForMCtruth(); + } + + die->SetCFManagerPair(cf); +} diff --git a/PWG3/dielectron/macros/ConfigJpsi2eeDataAOD.C b/PWG3/dielectron/macros/ConfigJpsi2eeDataAOD.C new file mode 100644 index 00000000000..ca430af6bba --- /dev/null +++ b/PWG3/dielectron/macros/ConfigJpsi2eeDataAOD.C @@ -0,0 +1,200 @@ + +void InitHistograms(AliDielectron *die, Int_t cutDefinition); +void InitCF(AliDielectron* die, Int_t cutDefinition); + +void SetupTrackCuts(AliDielectron *die, Int_t cutDefinition); +void SetupPairCuts(AliDielectron *die, Int_t cutDefinition); + +TString names=("trackQ+Pt>1GeV+loose_PID;trackQ+Pt>1GeV+tight_PID;trackQ+Pt>.5GeV"); +TObjArray *arrNames=names.Tokenize(";"); + +const Int_t nDie=arrNames->GetEntries(); + +AliDielectron* ConfigJpsi2ee(Int_t cutDefinition) +{ + // + // Setup the instance of AliDielectron + // + + // create the actual framework object + TString name=Form("%02d",cutDefinition); + if (cutDefinitionGetEntriesFast()){ + name=arrNames->At(cutDefinition)->GetName(); + } + AliDielectron *die = + new AliDielectron(Form("%s",name.Data()), + Form("Track cuts: %s",name.Data())); + + // cut setup + SetupTrackCuts(die,cutDefinition); + SetupPairCuts(die,cutDefinition); + + // + // histogram setup + // only if an AliDielectronHistos object is attached to the + // dielectron framework the QA histograms will be filled + // + InitHistograms(die,cutDefinition); + + // the last definition uses no cuts and only the QA histograms should be filled! + if (cutDefinitionGetTrackFilter().AddCuts(SetupESDtrackCuts(cutDefinition)); + + + //QA no CF + if (cutDefinition==nDie-1) { + //Pt cut + AliDielectronVarCuts *pt = new AliDielectronVarCuts("Pt>.5","Pt>.5"); + pt->AddCut(AliDielectronVarManager::kPt,.5,1e30); + die->GetTrackFilter().AddCuts(pt); + + return; + } + + + AliDielectronVarCuts *pt = new AliDielectronVarCuts("Pt>1","Pt>1"); + pt->AddCut(AliDielectronVarManager::kPt,1.,1e30); + pt->AddCut(AliDielectronVarManager::kImpactParZ,-3.,3.); + pt->AddCut(AliDielectronVarManager::kImpactParXY,-.07,.07); + die->GetTrackFilter().AddCuts(pt); + + //loose PID + if (cutDefinition==0){ + //ESD pid cuts (TPC nSigma Electrons) + AliDielectronVarCuts *pid = new AliDielectronVarCuts("TPCpid","TPC nSigma cut |e|<3 pi>2 |P|>2 |K|>2"); + //include + pid->AddCut(AliDielectronVarManager::kTPCnSigmaEle, -3., 3.); + //exclude + pid->AddCut(AliDielectronVarManager::kTPCnSigmaPio, -20., 2., kTRUE); + pid->AddCut(AliDielectronVarManager::kTPCnSigmaPro, -2., 2., kTRUE); + pid->AddCut(AliDielectronVarManager::kTPCnSigmaKao, -2., 2., kTRUE); + die->GetTrackFilter().AddCuts(pid); + } + + //tight PID + if (cutDefinition==1){ + //ESD pid cuts (TPC nSigma Electrons) + AliDielectronVarCuts *pid = new AliDielectronVarCuts("TPCpid","TPC nSigma cut |e|<2 pi>2 |P|>2 |K|>2; |dXY|<200um"); + //include + pid->AddCut(AliDielectronVarManager::kTPCnSigmaEle, -2., 2.); + //exclude + pid->AddCut(AliDielectronVarManager::kTPCnSigmaPio, -2., 2., kTRUE); + pid->AddCut(AliDielectronVarManager::kTPCnSigmaPro, -2., 2., kTRUE); + pid->AddCut(AliDielectronVarManager::kTPCnSigmaKao, -2., 2., kTRUE); + //tighten Impact parameter to 200um + pid->AddCut(AliDielectronVarManager::kImpactParXY, -0.02, 0.02); + die->GetTrackFilter().AddCuts(pid); + } + +} + +//______________________________________________________________________________________ +void SetupPairCuts(AliDielectron *die, Int_t cutDefinition) +{ + // + // Setup the pair cuts + // + + + // reject conversions + // and select mass region + AliDielectronVarCuts *openingAngleCut=new AliDielectronVarCuts("OpeningAngle","Opening angle > 35mrad"); + openingAngleCut->AddCut(AliDielectronVarManager::kOpeningAngle,.035,4.); + openingAngleCut->AddCut(AliDielectronVarManager::kM,2.,4.); + die->GetPairFilter().AddCuts(openingAngleCut); +} + +//______________________________________________________________________________________ +void InitHistograms(AliDielectron *die, Int_t cutDefinition) +{ + // + // Initialise the QA histograms + // + + //Setup QA histograms + AliDielectronHistos *histos= + new AliDielectronHistos(die->GetName(), + die->GetTitle()); + + //Initialise histogram classes + histos->SetReservedWords("Track;Pair"); + + //Event class (only for last QA) + if (cutDefinition==nDie-1) histos->AddClass("Event"); + + //Track classes, only first event + for (Int_t i=0; i<2; ++i){ + histos->AddClass(Form("Track_%s",AliDielectron::TrackClassName(i))); + } + + //Pair classes, only first event + for (Int_t i=0; i<3; ++i){ + histos->AddClass(Form("Pair_%s",AliDielectron::PairClassName(i))); + } + + //Event histograms + if (cutDefinition==nDie-1){ + //add histograms to event class + histos->UserHistogram("Event","nEvents","Number of processed events after cuts;Number events", + 1,0.,1.,AliDielectronVarManager::kNevents); + } + + //add histograms to Track classes + histos->UserHistogram("Track","Pt","Pt;Pt [GeV];#tracks",200,0,20.,AliDielectronVarManager::kPt); + histos->UserHistogram("Track","dEdx_P","dEdx;P [GeV];TPC signal (arb units);#tracks", + 400,0.2,20.,200,0.,200.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCsignal,kTRUE); + histos->UserHistogram("Track","TPCnSigmaEle_P","TPC number of sigmas Electrons;P [GeV];TPC number of sigmas Electrons;#tracks", + 400,0.2,20.,100,-5.,5.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaEle,kTRUE); + histos->UserHistogram("Track","TPCnSigmaPio_P","TPC number of sigmas Pions;P [GeV];TPC number of sigmas Pions;#tracks", + 400,0.2,20.,100,-5.,5.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaPio,kTRUE); + histos->UserHistogram("Track","dXY","dXY;dXY [cm];#tracks",200,-2.,2.,AliDielectronVarManager::kImpactParXY); + + //add histograms to Pair classes + histos->UserHistogram("Pair","InvMass","Inv.Mass;Inv. Mass [GeV];#pairs", + 500,0.,4.,AliDielectronVarManager::kM); + histos->UserHistogram("Pair","Rapidity","Rapidity;Rapidity;#pairs", + 100,-2.,2.,AliDielectronVarManager::kY); + histos->UserHistogram("Pair","OpeningAngle","Opening angle;angle", + 100,0.,3.15,AliDielectronVarManager::kOpeningAngle); + histos->UserHistogram("Pair","Chi2/NDF","#Chi^{2}/NDF;#Chi^{2}/NDF", + 100, 0., 20., AliDielectronVarManager::kChi2NDF); + + die->SetHistogramManager(histos); +} + + +void InitCF(AliDielectron* die, Int_t cutDefinition) +{ + // + // Setupd the CF Manager if needed + // + AliDielectronCF *cf=new AliDielectronCF(die->GetName(),die->GetTitle()); + + //pair variables + cf->AddVariable(AliDielectronVarManager::kPt,20,0,10); + cf->AddVariable(AliDielectronVarManager::kY,40,-2,2); + cf->AddVariable(AliDielectronVarManager::kM,200,-0.01,3.99); + cf->AddVariable(AliDielectronVarManager::kPairType,10,0,10); + //leg variables + cf->AddVariable(AliDielectronVarManager::kPt,20,0,10,kTRUE); + + //only in this case write MC truth info + cf->SetStepsForCutsIncreasing(); + if (cutDefinition==0){ + cf->SetStepForMCtruth(); + } + + die->SetCFManagerPair(cf); +} -- 2.43.5