]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Updates and additions: Classes for signal and spectrum extraction; saving of
authorandronic <andronic@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 4 May 2010 09:29:41 +0000 (09:29 +0000)
committerandronic <andronic@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 4 May 2010 09:29:41 +0000 (09:29 +0000)
a debug tree (default off); AddTaskJPSI.C; config macros for data

31 files changed:
PWG3/dielectron/AddTaskJPSI.C [new file with mode: 0644]
PWG3/dielectron/AliAnalysisTaskMultiDielectron.cxx
PWG3/dielectron/AliAnalysisTaskMultiDielectron.h
PWG3/dielectron/AliDielectron.cxx
PWG3/dielectron/AliDielectron.h
PWG3/dielectron/AliDielectronCF.cxx
PWG3/dielectron/AliDielectronCF.h
PWG3/dielectron/AliDielectronCFdraw.cxx
PWG3/dielectron/AliDielectronCFdraw.h
PWG3/dielectron/AliDielectronDebugTree.cxx [new file with mode: 0644]
PWG3/dielectron/AliDielectronDebugTree.h [new file with mode: 0644]
PWG3/dielectron/AliDielectronHistos.cxx
PWG3/dielectron/AliDielectronHistos.h
PWG3/dielectron/AliDielectronMC.cxx
PWG3/dielectron/AliDielectronMC.h
PWG3/dielectron/AliDielectronPair.cxx
PWG3/dielectron/AliDielectronPair.h
PWG3/dielectron/AliDielectronSignalBase.cxx [new file with mode: 0644]
PWG3/dielectron/AliDielectronSignalBase.h [new file with mode: 0644]
PWG3/dielectron/AliDielectronSignalExt.cxx [new file with mode: 0644]
PWG3/dielectron/AliDielectronSignalExt.h [new file with mode: 0644]
PWG3/dielectron/AliDielectronSignalFunc.cxx [new file with mode: 0644]
PWG3/dielectron/AliDielectronSignalFunc.h [new file with mode: 0644]
PWG3/dielectron/AliDielectronSpectrum.cxx [new file with mode: 0644]
PWG3/dielectron/AliDielectronSpectrum.h [new file with mode: 0644]
PWG3/dielectron/AliDielectronVarCuts.cxx
PWG3/dielectron/AliDielectronVarCuts.h
PWG3/dielectron/AliDielectronVarManager.cxx
PWG3/dielectron/AliDielectronVarManager.h
PWG3/dielectron/macros/ConfigJpsi2eeData.C [new file with mode: 0644]
PWG3/dielectron/macros/ConfigJpsi2eeDataAOD.C [new file with mode: 0644]

diff --git a/PWG3/dielectron/AddTaskJPSI.C b/PWG3/dielectron/AddTaskJPSI.C
new file mode 100644 (file)
index 0000000..9dc535c
--- /dev/null
@@ -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; i<nDie; ++i){ //nDie defined in config file
+    AliDielectron *jpsi=ConfigJpsi2ee(i);
+    task->AddDielectron(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;
+}
index d2fa05fb5ec95480946fd141188cff9e52a656a1..da4c12ef056ae3cd2a639c36d8af43a67d973541 100644 (file)
 #include <TChain.h>
 
 #include <AliCFContainer.h>
+#include <AliInputEventHandler.h>
+#include <AliESDInputHandler.h>
+#include <AliAnalysisManager.h>
 #include <AliVEvent.h>
 
 #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<AliESDInputHandler*>(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<AliDielectron*>(nextDie())) ){
+    die->SaveDebugTree();
+  }
+}
+
index 7762bf2bb662c2f8fedad07d671c02cb0fc47ddf..a9942d715a7662489383b995ad7db917bdf06c3e 100644 (file)
@@ -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);
index e95e953f6648f5f1cbd3150003c9a7291e353e37..a76cb2e92e72d4aabf55156ba61fc9e6972a7c65 100644 (file)
@@ -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; ipair<ntracks; ++ipair){
+      fDebugTree->Fill(static_cast<AliDielectronPair*>(PairArray(i)->UncheckedAt(ipair)));
+    }
+  }
+}
 
-
+//________________________________________________________________
+void AliDielectron::SaveDebugTree()
+{
+  //
+  // delete the debug tree, this will also write the tree
+  //
+  if (fDebugTree) fDebugTree->DeleteStreamer();
+}
 
index c3c625038b5dc0c7a621ce8a4aaa87b275b3e65d..6f6f8267283bcaf889bc2801ee34422a95149183 100644 (file)
@@ -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()
index e377cffa4bf6c7adda2037b8bea8fe7015ec1284..02e689ae92923dd5b748ee713b90c811954cd4e9 100644 (file)
@@ -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;i<fNVars;++i)    nbins[i]=fNBins[i];
@@ -227,22 +248,20 @@ void AliDielectronCF::InitialiseContainer(const AliAnalysisFilter& filter)
     fCfContainer->SetStepTitle(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; iCut<fNCuts;++iCut) {
       cutName=filter.GetCuts()->At(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; iCut<fNCuts-1;++iCut) {
       cutName+="&";
       cutName+=filter.GetCuts()->At(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; iCut<fNCuts;++iCut) {
+        cutName+="&";
+        cutName+=filter.GetCuts()->At(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; iCut<fNCuts;++iCut) {
       if (mask&(1<<iCut)) {
         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;
       }
     }
   }
@@ -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!!!");
+  }
+  
 }
 
 //________________________________________________________________
index 609cd42834d9cf694562764ffb205a5eff667cb7..c7f9d6f86d88a615796b74ec46916ea5bf18b493 100644 (file)
@@ -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);
index 323ae408a8a4f1f3a1565cc3d260a27897551a56..5656f9cb5b01c937ff8f413de2f0d9fc4ad5821e 100644 (file)
 #include <TList.h>
 #include <TClass.h>
 #include <TObject.h>
+#include <TVirtualPS.h>
 #include <TFile.h>
 #include <TString.h>
 #include <TObjString.h>
+#include <TVectorD.h>
 #include <TMath.h>
 #include <TH1.h>
 #include <TH2.h>
@@ -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; istep<nstep; ++istep) delete fCfContainer->GetGrid(istep);
+//   for (Int_t istep=0; istep<nstep; ++istep) delete fCfContainer->GetGrid(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; istep<fCfContainer->GetNStep(); ++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<TObjString*>(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; istep<fCfContainer->GetNStep(); ++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<TObjString*>(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; istep<fCfContainer->GetNStep(); ++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<TObjString*>(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; istep<fCfContainer->GetNStep(); ++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<TObjString*>(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<nPads; ++i){
-      if (i==1) optStr+="same";
+      if (i==1&&!drawSamePlus) optStr+="same";
       TH1 *hist=static_cast<TH1*>(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);
     }
   }
   
index b662f2fce991a4b80982e5c6c437d6a2af52df2b..142350b92518244e7bb7dbe84f811d266e6236e0 100644 (file)
@@ -22,6 +22,7 @@
 
 
 #include <TNamed.h>
+#include <TVectorD.h>
 #include <AliCFContainer.h>
 
 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 (file)
index 0000000..5505702
--- /dev/null
@@ -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 <TFile.h>
+#include <TTreeStream.h>
+#include <TString.h>
+
+#include <AliAnalysisManager.h>
+
+#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; i<AliDielectronVarManager::kNMaxValues;++i){
+    fVariables[i]=0;
+    fVariablesLeg[i]=0;
+  }
+}
+
+//______________________________________________
+AliDielectronDebugTree::AliDielectronDebugTree(const char* name, const char* title) :
+  TNamed(name, title),
+  fFileName("jpsi_debug.root"),
+  fNVars(0),
+  fNVarsLeg(0),
+  fStreamer(0x0)
+{
+  //
+  // Named Constructor
+  //
+  for (Int_t i=0; i<AliDielectronVarManager::kNMaxValues;++i){
+    fVariables[i]=0;
+    fVariablesLeg[i]=0;
+  }
+}
+
+//______________________________________________
+AliDielectronDebugTree::~AliDielectronDebugTree()
+{
+  //
+  // Default Destructor
+  //
+  if (fStreamer){
+    fStreamer->GetFile()->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; i<fNVars; ++i){
+      var=fVariables[i];
+      (*fStreamer) << "Pair"
+                   << Form("%s=",AliDielectronVarManager::GetValueName(var))
+                   << values[var];
+    }
+  }
+
+  if (fNVarsLeg>0){
+    //leg1
+    AliDielectronVarManager::Fill(pair->GetFirstDaughter(),valuesLeg1);
+    //leg2
+    AliDielectronVarManager::Fill(pair->GetSecondDaughter(),valuesLeg2);
+    
+    for (Int_t i=0; i<fNVarsLeg; ++i){
+      var=fVariablesLeg[i];
+      (*fStreamer) << "Pair"
+                   << Form("Leg1_%s=",AliDielectronVarManager::GetValueName(var))
+                   << valuesLeg1[var]
+                   << Form("Leg2_%s=",AliDielectronVarManager::GetValueName(var))
+                   << valuesLeg2[var];
+    }
+    
+  }
+  
+  (*fStreamer) << "Pair" << "\n";
+    
+  
+}
+
+//______________________________________________
+void AliDielectronDebugTree::DeleteStreamer()
+{
+  //
+  // delete the streamer
+  //
+  if (!fStreamer) return;
+  delete fStreamer;
+  fStreamer=0x0;
+
+}
+
+//______________________________________________
+void AliDielectronDebugTree::WriteTree()
+{
+  //
+  // Write file
+  //
+  if (!fStreamer) return;
+  fStreamer->GetFile()->Write();
+}
diff --git a/PWG3/dielectron/AliDielectronDebugTree.h b/PWG3/dielectron/AliDielectronDebugTree.h
new file mode 100644 (file)
index 0000000..bf46792
--- /dev/null
@@ -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 <TNamed.h>
+#include <TString.h>
+
+#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
index 8fab66bbb0250a88fa02934b699c8f7e4f4bfaa4..160d673ac5fc747b5d9a1896ef569d31353d42f5 100644 (file)
@@ -36,6 +36,8 @@
 #include <TROOT.h>
 #include <TLegend.h>
 #include <TKey.h>
+#include <TAxis.h>
+#include <TVirtualPS.h>
 // #include <TVectorD.h>
 
 #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<nbinsX+1; ++i){
+    binLim[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
+  }
+  return binLim;
+}
+
+//_____________________________________________________________________________
+Double_t* AliDielectronHistos::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) const
+{
+  //
+  // Make logarithmic binning
+  // the user has to delete the array afterwards!!!
+  //
+  Double_t *binLim=new Double_t[nbinsX+1];
+  Double_t first=xmin;
+  Double_t last=xmax;
+  Double_t binWidth=(last-first)/nbinsX;
+  for (Int_t i=0; i<nbinsX+1; ++i){
+    binLim[i]=first+binWidth*(Double_t)i;
+  }
+  return binLim;
+}
+
index 1be96bef7e175187cb828348fdbf90b26f53bc88..3e47ca9f409811f019e18199a2fdcd093305d0f7 100644 (file)
@@ -31,16 +31,18 @@ public:
   
   void UserHistogram(const char* histClass,const char *name, const char* title,
                      Int_t nbinsX, Double_t xmin, Double_t xmax,
-                     UInt_t valTypeX=kNoAutoFill);
+                     UInt_t valTypeX=kNoAutoFill, Bool_t logBinX=kFALSE);
   void 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=kNoAutoFill, UInt_t valTypeY=0);
+                     UInt_t valTypeX=kNoAutoFill, UInt_t valTypeY=0,
+                     Bool_t logBinX=kFALSE, Bool_t logBinY=kFALSE);
   void 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,
                      Int_t nbinsZ, Double_t zmin, Double_t zmax,
-                     UInt_t valTypeX=kNoAutoFill, UInt_t valTypeY=0, UInt_t valTypeZ=0);
+                     UInt_t valTypeX=kNoAutoFill, UInt_t valTypeY=0, UInt_t valTypeZ=0,
+                     Bool_t logBinX=kFALSE, Bool_t logBinY=kFALSE, Bool_t logBinZ=kFALSE);
   
   void UserHistogram(const char* histClass, TH1* hist, UInt_t valTypes=kNoAutoFill);
 
@@ -56,6 +58,7 @@ public:
   TH1* GetHistogram(const char* histClass, const char* name) const;
 
   void SetHistogramList(THashList &list);
+  void ResetHistogramList(){fHistoList.Clear();}
   const THashList* GetHistogramList() const {return &fHistoList;}
   
   void AddClass(const char* histClass);
@@ -74,7 +77,7 @@ public:
 //   virtual TObject  **GetObjectRef(const TObject *obj) const { return 0; }
 //   virtual TIterator *MakeIterator(Bool_t dir = kIterForward) const ;
 //   virtual TObject   *Remove(TObject *obj) { return 0; }
-  
+
 private:
   THashList fHistoList;             //-> 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);
index 625bd827583e5105c146aa92bc7bf7523a61569d..e2b55fefaf5350daab51a3e777d61e6a70d2f0f9 100644 (file)
@@ -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<AliMCEventHandler*> (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;
index fdf5ae930c423395ce71f7b6e357def850fc2dcc..0079dcaa9ada597eae048542bbbd467adf3d744a 100644 (file)
@@ -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
   
index 96d3e228965fc9bc4aa2c97ee3cf5a8bf693c8b5..ef30d2795b2cb18921c9edd895491eeac0e7eb83 100644 (file)
 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);
 }
 
index 64a496150da10f7d32b940698e8f2544201cacd0..fa5addd15eb4b3ee01c248f07e59db83ba691bcd 100644 (file)
@@ -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<AliVParticle*>(fRefD1.GetObject()); }
   AliVParticle* GetSecondDaughter()  const { return dynamic_cast<AliVParticle*>(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 (file)
index 0000000..44adbd3
--- /dev/null
@@ -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 (file)
index 0000000..8826080
--- /dev/null
@@ -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 <TNamed.h>
+#include <TVectorT.h>
+#include <TMath.h>
+
+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 (file)
index 0000000..0d12f9b
--- /dev/null
@@ -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 <TF1.h>
+#include <TH1.h>
+#include <TCanvas.h>
+#include <TMath.h>
+#include <TString.h>
+
+#include <AliLog.h>
+
+#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; ibin<fBins+1; ibin++ ){
+    Double_t mass  = ibin*(maxX-minX)/(Double_t)fBins;
+    Double_t backgr = 2.*sqrt( fSignPP->GetBinContent(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 (file)
index 0000000..028bb10
--- /dev/null
@@ -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 <TVectorT.h>
+#include <TString.h>
+#include <TH1.h>
+
+#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 (file)
index 0000000..3180cb1
--- /dev/null
@@ -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 <TF1.h>
+#include <TH1.h>
+#include <TGraph.h>
+#include <TMath.h>
+#include <TString.h>
+#include <TPaveText.h>
+#include <RooRealVar.h>
+#include <RooCBShape.h>
+#include <RooExponential.h>
+
+#include <AliLog.h>
+
+#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 (file)
index 0000000..de1b402
--- /dev/null
@@ -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 <TVectorT.h>
+#include <TString.h>
+#include <RooFit.h>
+#include <RooAddPdf.h>
+
+#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 (file)
index 0000000..6d1221a
--- /dev/null
@@ -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 <TMath.h>
+#include <TVectorT.h>
+#include <TH1.h>
+
+#include <AliCFContainer.h>
+#include <AliCFGridSparse.h>
+
+#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; iVar<fNvars; ++iVar){
+    fVars[iVar]=fCFSignal->GetVar(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; iVar<fNvars; iVar++) {
+    fCFSpectrum->SetBinLimits(iVar, fCFSignal->GetBinLimits(fVars[iVar]));
+    fCFSpectrum->SetVarTitle(iVar, fCFSignal->GetVarTitle(fVars[iVar]));
+  }
+
+  // setup step titles
+  Int_t steps=0;
+  for (Int_t iMethod=0; iMethod<fSignalMethods.GetEntries(); ++iMethod){
+    TString name(fSignalMethods.UncheckedAt(iMethod)->GetName());
+    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; iVar<fNvars; ++iVar){
+    vars[iVar]=fCFCorrection->GetVar(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; ibin<fNbins[variable]; ++ibin){
+    Int_t bin=ibin+1;
+    fCFSignal->GetGrid(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; iMethod<fSignalMethods.GetEntries(); ++iMethod){
+      AliDielectronSignalBase *signalMethod=(AliDielectronSignalBase*)fSignalMethods.At(iMethod);
+      signalMethod->Process(&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 (file)
index 0000000..938ca9b
--- /dev/null
@@ -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 <TObjArray.h>
+#include <TString.h>
+
+#include <TNamed.h>
+
+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
index 4316184169821133076066950e3c925701fda78f..d1ebc8650671a8e1d42830abb25d3691a31c4f87 100644 (file)
@@ -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; iCut<fNActiveCuts; ++iCut){
     Int_t cut=fActiveCuts[iCut];
     SETBIT(fSelectedCutsMask,iCut);
-    if ( (values[cut]<fCutMin[iCut]) || (values[cut]>fCutMax[iCut]) ) CLRBIT(fSelectedCutsMask,iCut);
+    if ( ((values[cut]<fCutMin[iCut]) || (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<fNActiveCuts; ++iCut){
     Int_t cut=(Int_t)fActiveCuts[iCut];
-    printf("Cut %02d: %f < %s < %f\n", iCut,
-           fCutMin[iCut], AliDielectronVarManager::GetValueName((Int_t)cut), fCutMax[iCut]);
+    Bool_t inverse=fCutExclude[iCut];
+
+    if (!inverse){
+      printf("Cut %02d: %f < %s < %f\n", iCut,
+             fCutMin[iCut], AliDielectronVarManager::GetValueName((Int_t)cut), fCutMax[iCut]);
+    } else {
+      printf("Cut %02d: !(%f < %s < %f)\n", iCut,
+             fCutMin[iCut], AliDielectronVarManager::GetValueName((Int_t)cut), fCutMax[iCut]);
+    }
   }
 }
index 61e46dcef24981138f09987df81d0d88c0553747..a5c089bfeda246bdc3dcbd65c81ca8b9513c7905 100644 (file)
@@ -22,8 +22,9 @@
 //#############################################################
 
 #include <Rtypes.h>
+
 #include <AliAnalysisCuts.h>
-#include <AliDielectronVarManager.h>
+#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
index e6dc32cf288b0de5e953bb5e9e5a5be09b00c06c..6f7c178fa79dd4a9f04ba689d67f749ccc9b9511 100644 (file)
@@ -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")
index 7e82f2c9b13486a528a6c07ef8f346e6fefe60b1..3a1eeb3b484d41a1918be6e82721de2470d6f59c 100644 (file)
@@ -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 <AliVParticle.h>
 #include <AliESDtrack.h>
 #include <AliAODTrack.h>
+#include <AliAODPid.h>
 #include <AliKFParticle.h>
+#include <AliKFVertex.h>
 #include <AliMCParticle.h>
 #include <AliAODMCParticle.h>
 #include <AliVTrack.h>  // ?
@@ -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&&i<kNMaxValues)?fgkParticleNames[i]:""; }
 private:
@@ -123,7 +141,7 @@ private:
   static void FillVarVParticle(const AliVParticle *particle,         Double_t * const values);
   static void FillVarESDtrack(const AliESDtrack *particle,           Double_t * const values);
   static void FillVarAODTrack(const AliAODTrack *particle,           Double_t * const values);
-  static  void FillVarMCParticle(const AliMCParticle *particle,      Double_t * const values);
+  static void FillVarMCParticle(const AliMCParticle *particle,      Double_t * const values);
   static void FillVarAODMCParticle(const AliAODMCParticle *particle, Double_t * const values);
   static void FillVarDielectronPair(const AliDielectronPair *pair,   Double_t * const values);
   static void FillVarVEvent(const AliVEvent *event,                  Double_t * const values);
@@ -133,7 +151,8 @@ private:
 
   static AliESDpid* fgESDpid;                 // ESD pid object
   static AliVEvent* fgEvent;                  // current event pointer
-
+  static AliKFVertex *fgKFVertex;          // kf vertex
+  
   AliDielectronVarManager(const AliDielectronVarManager &c);
   AliDielectronVarManager &operator=(const AliDielectronVarManager &c);
   
@@ -225,6 +244,16 @@ inline void AliDielectronVarManager::FillVarESDtrack(const AliESDtrack *particle
   // TODO: for the moment we set the bethe bloch parameters manually
   //       this should be changed in future!
   values[AliDielectronVarManager::kTPCnSigmaEle]=fgESDpid->NumberOfSigmasTPC(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 (file)
index 0000000..35a5a83
--- /dev/null
@@ -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 (cutDefinition<arrNames->GetEntriesFast()){
+    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 (cutDefinition<nDie-1) InitCF(die,cutDefinition);
+  
+  return die;
+}
+
+//______________________________________________________________________________________
+void SetupTrackCuts(AliDielectron *die, Int_t cutDefinition)
+{
+  //
+  // Setup the track cuts
+  //
+  
+  //ESD quality cuts
+  die->GetTrackFilter().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 (file)
index 0000000..ca430af
--- /dev/null
@@ -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 (cutDefinition<arrNames->GetEntriesFast()){
+    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 (cutDefinition<nDie-1) InitCF(die,cutDefinition);
+  
+  return die;
+}
+
+//______________________________________________________________________________________
+void SetupTrackCuts(AliDielectron *die, Int_t cutDefinition)
+{
+  //
+  // Setup the track cuts
+  //
+  
+  //ESD quality cuts
+  die->GetTrackFilter().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);
+}