From: andronic Date: Fri, 10 Dec 2010 22:44:06 +0000 (+0000) Subject: Major dielectron framework update; includes "alignment" to updates in X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=ffbede40cd4ed112a52ccd143d96c07ea4ffb65c Major dielectron framework update; includes "alignment" to updates in CF; config for AOD filtering for PbPb --- diff --git a/PWG3/PWG3dielectronLinkDef.h b/PWG3/PWG3dielectronLinkDef.h index 8c71a55970a..0e51dc718ec 100644 --- a/PWG3/PWG3dielectronLinkDef.h +++ b/PWG3/PWG3dielectronLinkDef.h @@ -14,6 +14,7 @@ #pragma link C++ class AliDielectronMC+; #pragma link C++ class AliDielectronVarManager+; #pragma link C++ class AliAnalysisTaskDielectronSE+; +#pragma link C++ class AliAnalysisTaskDielectronME+; #pragma link C++ class AliAnalysisTaskDielectronFilter+; #pragma link C++ class AliAnalysisTaskDielectronEfficiency+; #pragma link C++ class AliAnalysisTaskMultiDielectron+; @@ -29,4 +30,5 @@ #pragma link C++ class AliDielectronPID+; #pragma link C++ class AliDielectronCutGroup+; #pragma link C++ class AliDielectronEventCuts+; +#pragma link C++ class AliDielectronHelper+; #endif diff --git a/PWG3/dielectron/AliAnalysisTaskDielectronFilter.cxx b/PWG3/dielectron/AliAnalysisTaskDielectronFilter.cxx index c14218ef862..a64cfe09b78 100644 --- a/PWG3/dielectron/AliAnalysisTaskDielectronFilter.cxx +++ b/PWG3/dielectron/AliAnalysisTaskDielectronFilter.cxx @@ -24,6 +24,7 @@ #include #include +#include #include #include #include @@ -134,7 +135,7 @@ void AliAnalysisTaskDielectronFilter::UserExec(Option_t *) } } //AOD case - if (man->GetInputEventHandler()->IsA()==AliESDInputHandler::Class()){ + if (man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class()){ if (!AliDielectronVarManager::GetAODpidUtil()){ if (AliDielectronMC::Instance()->HasMC()) { AliDielectronVarManager::InitAODpidUtil(); diff --git a/PWG3/dielectron/AliAnalysisTaskDielectronME.cxx b/PWG3/dielectron/AliAnalysisTaskDielectronME.cxx new file mode 100644 index 00000000000..3d95e0f8fc4 --- /dev/null +++ b/PWG3/dielectron/AliAnalysisTaskDielectronME.cxx @@ -0,0 +1,193 @@ +/************************************************************************* +* 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. * +**************************************************************************/ + +/////////////////////////////////////////////////////////////////////////// +// // +// Basic Analysis Task // +// // +/////////////////////////////////////////////////////////////////////////// + +#include +#include + +#include +#include +#include +#include +#include + +#include "AliDielectron.h" +#include "AliDielectronHistos.h" +#include "AliDielectronCF.h" +#include "AliDielectronMC.h" +#include "AliAnalysisTaskDielectronME.h" + +ClassImp(AliAnalysisTaskDielectronME) + +//_________________________________________________________________________________ +AliAnalysisTaskDielectronME::AliAnalysisTaskDielectronME() : + AliAnalysisTaskME(), + fListDielectron(), + fListHistos(), + fListCF(), + fPoolDepth(2), + fSelectPhysics(kFALSE), + fTriggerMask(AliVEvent::kMB), + fEventStat(0x0) +{ + // + // Constructor + // +} + +//_________________________________________________________________________________ +AliAnalysisTaskDielectronME::AliAnalysisTaskDielectronME(const char *name) : + AliAnalysisTaskME(name), + fListDielectron(), + fListHistos(), + fListCF(), + fPoolDepth(2), + fSelectPhysics(kFALSE), + fTriggerMask(AliVEvent::kMB), + fEventStat(0x0) +{ + // + // Constructor + // + DefineInput(0,TChain::Class()); + DefineOutput(1, TList::Class()); + DefineOutput(2, TList::Class()); + DefineOutput(3, TH1D::Class()); + fListHistos.SetName("Dielectron_Histos_Multi"); + fListCF.SetName("Dielectron_CF_Multi"); +} + + +//_________________________________________________________________________________ +void AliAnalysisTaskDielectronME::UserCreateOutputObjects() +{ + // + // Add all histogram manager histogram lists to the output TList + // + + if (!fListHistos.IsEmpty()||!fListCF.IsEmpty()) return; //already initialised + + TIter nextDie(&fListDielectron); + AliDielectron *die=0; + while ( (die=static_cast(nextDie())) ){ + die->Init(); + if (die->GetHistogramList()) fListHistos.Add(const_cast(die->GetHistogramList())); + if (die->GetCFManagerPair()) fListCF.Add(const_cast(die->GetCFManagerPair()->GetContainer())); + } + + Int_t cuts=fListDielectron.GetEntries(); + Int_t nbins=2+2*cuts; + if (!fEventStat){ + fEventStat=new TH1D("hEventStat","Event statistics",nbins,0,nbins); + fEventStat->GetXaxis()->SetBinLabel(1,"Before Phys. Sel."); + fEventStat->GetXaxis()->SetBinLabel(2,"After Phys. Sel."); + for (Int_t i=0; iGetXaxis()->SetBinLabel(3+2*i,Form("#splitline{1 candidate}{%s}",fListDielectron.At(i)->GetName())); + fEventStat->GetXaxis()->SetBinLabel(4+2*i,Form("#splitline{With >1 candidate}{%s}",fListDielectron.At(i)->GetName())); + } + } + + PostData(1, &fListHistos); + PostData(2, &fListCF); + PostData(3, fEventStat); +} + +//_________________________________________________________________________________ +void AliAnalysisTaskDielectronME::UserExec(Option_t *) +{ + // + // Main loop. Called for every event + // + + if (fListHistos.IsEmpty()&&fListCF.IsEmpty()) return; + + AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager(); + AliESDInputHandler *esdHandler=0x0; + if ( (esdHandler=dynamic_cast(man->GetInputEventHandler())) && esdHandler->GetESDpid() ){ + AliDielectronVarManager::SetESDpid(esdHandler->GetESDpid()); + } else { + //load esd pid bethe bloch parameters depending on the existance of the MC handler + // yes: MC parameters + // no: data parameters + if (!AliDielectronVarManager::GetESDpid()){ + if (AliDielectronMC::Instance()->HasMC()) { + AliDielectronVarManager::InitESDpid(); + } else { + AliDielectronVarManager::InitESDpid(1); + } + } + } + // Was event selected ? + AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler()); + UInt_t isSelected = AliVEvent::kAny; + if( fSelectPhysics && inputHandler && inputHandler->GetEventSelection() ) { + isSelected = inputHandler->IsEventSelected(); + isSelected&=fTriggerMask; + } + + //Before physics selection + fEventStat->Fill(0.); + if (isSelected==0) { + PostData(3,fEventStat); + return; + } + //after physics selection + fEventStat->Fill(1.); + + //bz for AliKF + Double_t bz = GetEvent(0)->GetMagneticField(); + AliKFParticle::SetField( bz ); + + //Process event in all AliDielectron instances + TIter nextDie(&fListDielectron); + AliDielectron *die=0; + Int_t idie=0; + while ( (die=static_cast(nextDie())) ){ + for (Int_t evt1=0; evt1Process((AliESDEvent*)GetEvent(evt1),(AliESDEvent*)GetEvent(evt2)); + if (die->HasCandidates()){ + Int_t ncandidates=die->GetPairArray(1)->GetEntriesFast(); + if (ncandidates==1) fEventStat->Fill(3+2*idie); + else if (ncandidates>1) fEventStat->Fill(4+2*idie); + } + } + } + ++idie; + } + + PostData(1, &fListHistos); + PostData(2, &fListCF); + PostData(3,fEventStat); +} + +//_________________________________________________________________________________ +void AliAnalysisTaskDielectronME::FinishTaskOutput() +{ + // + // Write debug tree + // + TIter nextDie(&fListDielectron); + AliDielectron *die=0; + while ( (die=static_cast(nextDie())) ){ + die->SaveDebugTree(); + } +} + diff --git a/PWG3/dielectron/AliAnalysisTaskDielectronME.h b/PWG3/dielectron/AliAnalysisTaskDielectronME.h new file mode 100644 index 00000000000..54c42ae10d3 --- /dev/null +++ b/PWG3/dielectron/AliAnalysisTaskDielectronME.h @@ -0,0 +1,65 @@ +#ifndef ALIANALYSISTASKDIELECTRONME_H +#define ALIANALYSISTASKDIELECTRONME_H +/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//##################################################### +//# # +//# Basic Analysis task for Dielectron # +//# single event analysis # +//# # +//# by WooJin J. Park, GSI / W.J.Park@gsi.de # +//# Ionut C. Arsene, GSI / I.C.Arsene@gsi.de # +//# Magnus Mager, CERN / Magnus.Mager@cern.ch # +//# Jens Wiechula, Uni HD / Jens.Wiechula@cern.ch # +//# # +//##################################################### + +#include "TList.h" + +#include "AliAnalysisTaskME.h" + +#include "AliDielectronPID.h" + +class AliDielectron; +class TH1D; + +class AliAnalysisTaskDielectronME : public AliAnalysisTaskME { + +public: + AliAnalysisTaskDielectronME(); + AliAnalysisTaskDielectronME(const char *name); + virtual ~AliAnalysisTaskDielectronME(){ } + + virtual void UserExec(Option_t *option); + virtual void UserCreateOutputObjects(); + virtual void FinishTaskOutput(); + //temporary + //virtual void NotifyRun(){AliDielectronPID::SetCorrVal((Double_t)fCurrentRunNumber);} + virtual void NotifyRun(){AliDielectronPID::SetCorrVal((Double_t)GetEvent(0)->GetRunNumber());} + + void UsePhysicsSelection(Bool_t phy=kTRUE) {fSelectPhysics=phy;} + void SetTriggerMask(UInt_t mask) {fTriggerMask=mask;} + UInt_t GetTriggerMask() const { return fTriggerMask; } + void SetPoolDepth(Int_t depth=2){fPoolDepth=depth;} + + void AddDielectron(AliDielectron * const die) { fListDielectron.Add(die); } + +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 + + Int_t fPoolDepth; // Pool depth for event mixing + Bool_t fSelectPhysics; // Whether to use physics selection + UInt_t fTriggerMask; // Event trigger mask + + TH1D *fEventStat; //! Histogram with event statistics + + AliAnalysisTaskDielectronME(const AliAnalysisTaskDielectronME &c); + AliAnalysisTaskDielectronME& operator= (const AliAnalysisTaskDielectronME &c); + + ClassDef(AliAnalysisTaskDielectronME, 1); //Analysis Task handling multiple instances of AliDielectron +}; +#endif diff --git a/PWG3/dielectron/AliAnalysisTaskDielectronSE.cxx b/PWG3/dielectron/AliAnalysisTaskDielectronSE.cxx index 5581411c54b..83b7bfaabe5 100644 --- a/PWG3/dielectron/AliAnalysisTaskDielectronSE.cxx +++ b/PWG3/dielectron/AliAnalysisTaskDielectronSE.cxx @@ -28,6 +28,7 @@ #include #include #include +#include #include "AliDielectron.h" #include "AliDielectronHistos.h" @@ -124,7 +125,7 @@ void AliAnalysisTaskDielectronSE::UserExec(Option_t *) } } //AOD case - if (man->GetInputEventHandler()->IsA()==AliESDInputHandler::Class()){ + if (man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class()){ if (!AliDielectronVarManager::GetAODpidUtil()){ if (AliDielectronMC::Instance()->HasMC()) { AliDielectronVarManager::InitAODpidUtil(); diff --git a/PWG3/dielectron/AliAnalysisTaskMultiDielectron.cxx b/PWG3/dielectron/AliAnalysisTaskMultiDielectron.cxx index 2638679ca56..2c712fecdb7 100644 --- a/PWG3/dielectron/AliAnalysisTaskMultiDielectron.cxx +++ b/PWG3/dielectron/AliAnalysisTaskMultiDielectron.cxx @@ -25,7 +25,7 @@ #include #include #include -#include +#include #include #include #include @@ -97,7 +97,7 @@ void AliAnalysisTaskMultiDielectron::UserCreateOutputObjects() AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager(); Bool_t isESD=man->GetInputEventHandler()->IsA()==AliESDInputHandler::Class(); -// Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODHandler::Class(); +// Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class(); TIter nextDie(&fListDielectron); AliDielectron *die=0; @@ -144,7 +144,7 @@ void AliAnalysisTaskMultiDielectron::UserExec(Option_t *) AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager(); AliESDInputHandler *esdHandler=0x0; Bool_t isESD=man->GetInputEventHandler()->IsA()==AliESDInputHandler::Class(); - Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODHandler::Class(); + Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class(); if ( (esdHandler=dynamic_cast(man->GetInputEventHandler())) && esdHandler->GetESDpid() ){ AliDielectronVarManager::SetESDpid(esdHandler->GetESDpid()); } else { @@ -223,8 +223,8 @@ void AliAnalysisTaskMultiDielectron::UserExec(Option_t *) die->Process(InputEvent()); if (die->HasCandidates()){ Int_t ncandidates=die->GetPairArray(1)->GetEntriesFast(); - if (ncandidates==1) fEventStat->Fill((kNbinsEvent+1)+2*idie); - else if (ncandidates>1) fEventStat->Fill((kNbinsEvent+2)+2*idie); + if (ncandidates==1) fEventStat->Fill((kNbinsEvent)+2*idie); + else if (ncandidates>1) fEventStat->Fill((kNbinsEvent+1)+2*idie); } ++idie; } diff --git a/PWG3/dielectron/AliAnalysisTaskMultiDielectron.h b/PWG3/dielectron/AliAnalysisTaskMultiDielectron.h index d5cb045919c..1d8b510302f 100644 --- a/PWG3/dielectron/AliAnalysisTaskMultiDielectron.h +++ b/PWG3/dielectron/AliAnalysisTaskMultiDielectron.h @@ -19,7 +19,7 @@ #include "AliAnalysisTaskSE.h" -#include "AliDielectronPID.h" +// #include "AliDielectronPID.h" class AliDielectron; class TH1D; diff --git a/PWG3/dielectron/AliDielectron.cxx b/PWG3/dielectron/AliDielectron.cxx index 7ab83df40de..77fec87a4c2 100644 --- a/PWG3/dielectron/AliDielectron.cxx +++ b/PWG3/dielectron/AliDielectron.cxx @@ -328,22 +328,27 @@ void AliDielectron::FillHistograms(const AliVEvent *ev) } //________________________________________________________________ -void AliDielectron::FillHistogramsPair(AliDielectronPair *pair) +void AliDielectron::FillHistogramsPair(AliDielectronPair *pair,Bool_t fromPreFilter/*=kFALSE*/) { // // Fill Histogram information for pairs and the track in the pair // NOTE: in this funtion the leg information may be filled multiple // times. This funtion is used in the track rotation pairing // and those legs are not saved! - //=fHistos + // TString className,className2; Double_t values[AliDielectronVarManager::kNMaxValues]; //Fill Pair information, separately for all pair candidate arrays and the legs TObjArray arrLegs(100); const Int_t type=pair->GetType(); - className.Form("Pair_%s",fgkPairClassNames[type]); - className2.Form("Track_Legs_%s",fgkPairClassNames[type]); + if (fromPreFilter) { + className.Form("RejPair_%s",fgkPairClassNames[type]); + className2.Form("RejTrack_%s",fgkPairClassNames[type]); + } else { + className.Form("Pair_%s",fgkPairClassNames[type]); + className2.Form("Track_Legs_%s",fgkPairClassNames[type]); + } Bool_t pairClass=fHistos->GetHistogramList()->FindObject(className.Data())!=0x0; Bool_t legClass=fHistos->GetHistogramList()->FindObject(className2.Data())!=0x0; @@ -436,6 +441,7 @@ void AliDielectron::PairPreFilter(Int_t arr1, Int_t arr2, TObjArray &arrTracks1, if (cutMask!=selectedMask) continue; if (fCfManagerPair) fCfManagerPair->Fill(selectedMaskPair+1 ,&candidate); accepted=kTRUE; + FillHistogramsPair(&candidate,kTRUE); //remove the tracks from the Track arrays arrTracks2.AddAt(0x0,itrack2); //in case of like sign remove the track from both arrays! @@ -548,7 +554,7 @@ void AliDielectron::FillPairArrayTR() while ( fTrackRotator->NextCombination() ){ AliDielectronPair candidate; candidate.SetTracks(fTrackRotator->GetTrackP(), fPdgLeg1, fTrackRotator->GetTrackN(), fPdgLeg2); - candidate.SetType(10); + candidate.SetType(kEv1PMRot); //pair cuts UInt_t cutMask=fPairFilter.IsSelected(&candidate); diff --git a/PWG3/dielectron/AliDielectron.h b/PWG3/dielectron/AliDielectron.h index 5be7eb492d3..24e36e3786a 100644 --- a/PWG3/dielectron/AliDielectron.h +++ b/PWG3/dielectron/AliDielectron.h @@ -23,13 +23,13 @@ #include #include "AliDielectronHistos.h" -#include "AliDielectronPair.h" class AliVEvent; class THashList; class AliDielectronCF; class AliDielectronDebugTree; class AliDielectronTrackRotator; +class AliDielectronPair; //________________________________________________________________ class AliDielectron : public TNamed { @@ -145,7 +145,7 @@ private: void ProcessMC(); void FillHistograms(const AliVEvent *ev); - void FillHistogramsPair(AliDielectronPair *pair); + void FillHistogramsPair(AliDielectronPair *pair,Bool_t fromPreFilter=kFALSE); void FillHistogramsTracks(TObjArray **tracks); void FillDebugTree(); diff --git a/PWG3/dielectron/AliDielectronCFdraw.cxx b/PWG3/dielectron/AliDielectronCFdraw.cxx index 7006b5a87d2..a05b4e95cfb 100644 --- a/PWG3/dielectron/AliDielectronCFdraw.cxx +++ b/PWG3/dielectron/AliDielectronCFdraw.cxx @@ -180,14 +180,14 @@ void AliDielectronCFdraw::SetRangeUser(Int_t ivar, Double_t min, Double_t max, c if (arr->GetEntriesFast()==0){ // all slices in case of 0 entries for (Int_t istep=0; istepGetNStep(); ++istep){ - fCfContainer->SetRangeUser(ivar,min,max,istep); + fCfContainer->GetGrid(istep)->SetRangeUser(ivar,min,max); } } else { TIter next(arr); TObjString *ostr=0x0; while ( (ostr=static_cast(next())) ) { Int_t istep=ostr->GetString().Atoi(); - fCfContainer->SetRangeUser(ivar,min,max,istep); + fCfContainer->GetGrid(istep)->SetRangeUser(ivar,min,max); } } delete arr; @@ -240,7 +240,7 @@ void AliDielectronCFdraw::Draw(const Option_t* varnames, const char* opt, const TIter next(arrVars); TObjString *ostr=0x0; Int_t ivar[3]={-1,-1,-1}; - for (Int_t i=0; i=0; --i){ ostr=static_cast(next()); if (ostr->GetString().IsDigit()){ ivar[i]=ostr->GetString().Atoi(); @@ -249,17 +249,7 @@ void AliDielectronCFdraw::Draw(const Option_t* varnames, const char* opt, const } } - switch (entries){ - case 1: - Draw(ivar[0],opt,slices); - break; - case 2: - Draw(ivar[1],ivar[0],opt,slices); - break; - case 3: - Draw(ivar[2],ivar[1],ivar[0],opt,slices); - break; - } + Draw(ivar[0],ivar[1],ivar[2],opt,slices); delete arrVars; } @@ -275,9 +265,8 @@ void AliDielectronCFdraw::Draw(Int_t var, const char* opt, const char* slices) // are drawn in each sub pad // - const Int_t ndim=1; - Int_t vars[ndim]={var}; - TObjArray *arr=CollectHistosProj(ndim,vars,slices); + Int_t vars[3]={var,-1,-1}; + TObjArray *arr=CollectHistosProj(vars,slices); Draw(arr,opt); delete arr; } @@ -288,9 +277,8 @@ void AliDielectronCFdraw::Draw(Int_t var0, Int_t var1, const char* opt, const ch // // Draw 2D case // - const Int_t ndim=2; - Int_t vars[ndim]={var0,var1}; - TObjArray *arr=CollectHistosProj(ndim,vars,slices); + Int_t vars[3]={var0,var1,-1}; + TObjArray *arr=CollectHistosProj(vars,slices); Draw(arr,opt); delete arr; } @@ -301,19 +289,19 @@ void AliDielectronCFdraw::Draw(Int_t var0, Int_t var1, Int_t var2, const char* o // // Draw 3D case // - const Int_t ndim=3; - Int_t vars[ndim]={var0,var1,var2}; - TObjArray *arr=CollectHistosProj(ndim,vars,slices); + Int_t vars[3]={var0,var1,var2}; + TObjArray *arr=CollectHistosProj(vars,slices); Draw(arr,opt); delete arr; } //________________________________________________________________ -TObjArray* AliDielectronCFdraw::CollectHistosProj(Int_t dim, Int_t *vars, const char* slices) +TObjArray* AliDielectronCFdraw::CollectHistosProj(const Int_t vars[3], const char* slices) { // - // Collect histos with 'dim'ension of the 'slices' separated by one of "':;'" + // Collect histos with up to 3 dimension of the 'slices' separated by one of "':;'" // in a TObjArray and return it + // if a dimension is not used it must be set to -1 // TObjArray *arr=TString(slices).Tokenize(",:;"); TObjArray *arrHists=0x0; @@ -321,7 +309,7 @@ TObjArray* AliDielectronCFdraw::CollectHistosProj(Int_t dim, Int_t *vars, const // all slices in case of 0 entries arrHists=new TObjArray(fCfContainer->GetNStep()); for (Int_t istep=0; istepGetNStep(); ++istep){ - TH1 *hproj=Project(dim,vars,istep); + TH1 *hproj=Project(vars,istep); if (!hproj) continue; hproj->SetName(Form("proj_%02d",istep)); hproj->SetTitle(fCfContainer->GetStepTitle(istep)); @@ -333,7 +321,7 @@ TObjArray* AliDielectronCFdraw::CollectHistosProj(Int_t dim, Int_t *vars, const TObjString *ostr=0x0; while ( (ostr=static_cast(next())) ) { Int_t istep=ostr->GetString().Atoi(); - TH1 *hproj=Project(dim,vars,istep); + TH1 *hproj=Project(vars,istep); if (!hproj) continue; hproj->SetName(Form("proj_%02d",istep)); hproj->SetTitle(fCfContainer->GetStepTitle(istep)); @@ -346,23 +334,12 @@ TObjArray* AliDielectronCFdraw::CollectHistosProj(Int_t dim, Int_t *vars, const } //________________________________________________________________ -TH1* AliDielectronCFdraw::Project(Int_t ndim, Int_t *vars, Int_t slice) +TH1* AliDielectronCFdraw::Project(const Int_t *vars, Int_t slice) { // // Do an ndim projection // - switch (ndim){ - case 1: - return fCfContainer->Project(vars[0],slice); - break; - case 2: - return fCfContainer->Project(vars[0],vars[1],slice); - break; - case 3: - return fCfContainer->Project(vars[0],vars[1],vars[2],slice); - break; - } - return 0x0; + return fCfContainer->Project(slice,vars[0],vars[1],vars[2]); } //________________________________________________________________ @@ -379,19 +356,19 @@ TH1* AliDielectronCFdraw::Project(const Option_t* var, Int_t slice) return 0x0; } - TIter next(arrVars); TObjString *ostr=0x0; Int_t ivar[3]={-1,-1,-1}; - for (Int_t i=0; i(next()); + for (Int_t i=entries-1; i>=0; --i){ + ostr=static_cast(arrVars->At(i)); if (ostr->GetString().IsDigit()){ ivar[i]=ostr->GetString().Atoi(); } else { ivar[i]=fCfContainer->GetVar(ostr->GetName()); } } + if (ivar[0]==-1) return 0x0; delete arrVars; - return Project(entries,ivar,slice); + return fCfContainer->Project(slice,ivar[0],ivar[1],ivar[2]); } //________________________________________________________________ @@ -426,17 +403,7 @@ void AliDielectronCFdraw::DrawEfficiency(const char* varnames, const char* numer TString optStr(opt); if (optStr.Contains("2")) type=1; - switch (entries){ - case 1: - DrawEfficiency(ivar[0],numerators, denominator,opt,type); - break; - case 2: - DrawEfficiency(ivar[1],ivar[0], numerators, denominator,opt,type); - break; - case 3: - DrawEfficiency(ivar[2],ivar[1],ivar[0],numerators, denominator,opt,type); - break; - } + DrawEfficiency(ivar[2],ivar[1],ivar[0],numerators, denominator,opt,type); delete arrVars; } @@ -452,9 +419,8 @@ void AliDielectronCFdraw::DrawEfficiency(Int_t var, const char* numerators, Int_ // are drawn in each sub pad // - const Int_t ndim=1; - Int_t vars[ndim]={var}; - TObjArray *arr=CollectHistosEff(ndim,vars,numerators,denominator,type); + Int_t vars[3]={var,-1,-1}; + TObjArray *arr=CollectHistosEff(vars,numerators,denominator,type); TString drawOpt=opt; drawOpt+="eff"; Draw(arr,drawOpt); @@ -467,9 +433,8 @@ void AliDielectronCFdraw::DrawEfficiency(Int_t var0, Int_t var1, const char* num // // Draw 2D case // - const Int_t ndim=2; - Int_t vars[ndim]={var0,var1}; - TObjArray *arr=CollectHistosEff(ndim,vars,numerators,denominator,type); + Int_t vars[3]={var0,var1,-1}; + TObjArray *arr=CollectHistosEff(vars,numerators,denominator,type); TString drawOpt=opt; drawOpt+="eff"; Draw(arr,drawOpt); @@ -482,9 +447,8 @@ void AliDielectronCFdraw::DrawEfficiency(Int_t var0, Int_t var1, Int_t var2, con // // Draw 3D case // - const Int_t ndim=3; - Int_t vars[ndim]={var0,var1,var2}; - TObjArray *arr=CollectHistosEff(ndim,vars,numerators,denominator,type); + Int_t vars[3]={var0,var1,var2}; + TObjArray *arr=CollectHistosEff(vars,numerators,denominator,type); TString drawOpt=opt; drawOpt+="eff"; Draw(arr,drawOpt); @@ -492,7 +456,7 @@ void AliDielectronCFdraw::DrawEfficiency(Int_t var0, Int_t var1, Int_t var2, con } //________________________________________________________________ -TObjArray* AliDielectronCFdraw::CollectHistosEff(Int_t dim, Int_t *vars, const char* numerators, Int_t denominator, Int_t type) +TObjArray* AliDielectronCFdraw::CollectHistosEff(const Int_t vars[3], const char* numerators, Int_t denominator, Int_t type) { // // Collect histos with 'dim'ension of the 'slices' separated by one of "':;'" @@ -508,7 +472,7 @@ TObjArray* AliDielectronCFdraw::CollectHistosEff(Int_t dim, Int_t *vars, const c fVdata.ResizeTo(arrHists->GetSize()); for (Int_t istep=0; istepGetNStep(); ++istep){ fEffGrid->CalculateEfficiency(istep,denominator); - TH1 *hproj=ProjectEff(dim,vars); + TH1 *hproj=ProjectEff(vars); if (!hproj) continue; Float_t eff=fEffGrid->GetAverage(); fVdata(istep)=eff; @@ -525,7 +489,7 @@ TObjArray* AliDielectronCFdraw::CollectHistosEff(Int_t dim, Int_t *vars, const c while ( (ostr=static_cast(next())) ) { Int_t istep=ostr->GetString().Atoi(); fEffGrid->CalculateEfficiency(istep,denominator); - TH1 *hproj=ProjectEff(dim,vars); + TH1 *hproj=ProjectEff(vars); if (!hproj) continue; Float_t eff=fEffGrid->GetAverage(); fVdata(count++)=eff; @@ -538,14 +502,14 @@ TObjArray* AliDielectronCFdraw::CollectHistosEff(Int_t dim, Int_t *vars, const c //second approach if (type==1){ - TH1 *hDen=Project(dim,vars,denominator); + TH1 *hDen=Project(vars,denominator); Double_t entriesDen=hDen->GetEffectiveEntries(); if (arr->GetEntriesFast()==0){ // all slices in case of 0 entries arrHists=new TObjArray(fCfContainer->GetNStep()); fVdata.ResizeTo(arrHists->GetSize()); for (Int_t istep=0; istepGetNStep(); ++istep){ - TH1 *hproj=Project(dim,vars,istep); + TH1 *hproj=Project(vars,istep); if (!hproj) continue; Float_t eff=0; if (entriesDen>0) eff=hproj->GetEffectiveEntries()/entriesDen; @@ -563,7 +527,7 @@ TObjArray* AliDielectronCFdraw::CollectHistosEff(Int_t dim, Int_t *vars, const c Int_t count=0; while ( (ostr=static_cast(next())) ) { Int_t istep=ostr->GetString().Atoi(); - TH1 *hproj=Project(dim,vars,istep); + TH1 *hproj=Project(vars,istep); if (!hproj) continue; Float_t eff=0; if (entriesDen>0) eff=hproj->GetEffectiveEntries()/entriesDen; @@ -583,23 +547,12 @@ TObjArray* AliDielectronCFdraw::CollectHistosEff(Int_t dim, Int_t *vars, const c } //________________________________________________________________ -TH1* AliDielectronCFdraw::ProjectEff(Int_t ndim, Int_t *vars) +TH1* AliDielectronCFdraw::ProjectEff(const Int_t vars[3]) { // // Do an nim projection // - switch (ndim){ - case 1: - return fEffGrid->Project(vars[0]); - break; - case 2: - return fEffGrid->Project(vars[0],vars[1]); - break; - case 3: - return fEffGrid->Project(vars[0],vars[1],vars[2]); - break; - } - return 0x0; + return fEffGrid->Project(vars[0],vars[1],vars[2]); } //________________________________________________________________ @@ -707,9 +660,9 @@ Double_t AliDielectronCFdraw::GetAverageEfficiency(Int_t numerator, Int_t denomi // //use variable 0 as default, since for the average it doesn't matter - TH1 *hDen=fCfContainer->Project(0,denominator); + TH1 *hDen=fCfContainer->Project(denominator,0); Double_t entriesDen=hDen->GetEffectiveEntries(); - TH1 *hproj=fCfContainer->Project(0,numerator); + TH1 *hproj=fCfContainer->Project(numerator,0); if (!hproj) return -1.; Double_t entriesNum=hproj->GetEffectiveEntries(); if (entriesDen<1||entriesNum<1) return -1; diff --git a/PWG3/dielectron/AliDielectronCFdraw.h b/PWG3/dielectron/AliDielectronCFdraw.h index 50bd4aee954..657786aad15 100644 --- a/PWG3/dielectron/AliDielectronCFdraw.h +++ b/PWG3/dielectron/AliDielectronCFdraw.h @@ -61,8 +61,8 @@ public: void Draw(Int_t var0, Int_t var1, const char* opt="", const char* slices=""); void Draw(Int_t var0, Int_t var1, Int_t var2, const char* opt="", const char* slices=""); - TObjArray* CollectHistosProj(Int_t dim, Int_t *vars, const char* slices); - TH1* Project(Int_t ndim, Int_t *vars, Int_t slice); + TObjArray* CollectHistosProj(const Int_t vars[3], const char* slices); + TH1* Project(const Int_t vars[3], Int_t slice); TH1* Project(const Option_t* var, Int_t slice); //Draw efficiencies @@ -71,8 +71,8 @@ public: void DrawEfficiency(Int_t var0, Int_t var1, const char* numerators, Int_t denominator=0, const char* opt="sameleg", Int_t type=0); void DrawEfficiency(Int_t var0, Int_t var1, Int_t var2, const char* numerators, Int_t denominator=0, const char* opt="sameleg", Int_t type=0); - TObjArray* CollectHistosEff(Int_t dim, Int_t *vars, const char* numerators, Int_t denominator, Int_t type=0); - TH1* ProjectEff(Int_t ndim, Int_t *vars); + TObjArray* CollectHistosEff(const Int_t vars[3], const char* numerators, Int_t denominator, Int_t type=0); + TH1* ProjectEff(const Int_t vars[3]); Double_t GetAverageEfficiency(Int_t numerator, Int_t denominator, Double_t &effErr); diff --git a/PWG3/dielectron/AliDielectronEventCuts.cxx b/PWG3/dielectron/AliDielectronEventCuts.cxx index 12de258542b..48049b5bd37 100644 --- a/PWG3/dielectron/AliDielectronEventCuts.cxx +++ b/PWG3/dielectron/AliDielectronEventCuts.cxx @@ -28,6 +28,7 @@ Detailed description #include #include #include +#include #include "AliDielectronEventCuts.h" @@ -39,6 +40,7 @@ AliDielectronEventCuts::AliDielectronEventCuts() : fVtxZmax(0.), fRequireVtx(kFALSE), fMinVtxContributors(0), + fMultITSTPC(kFALSE), fVtxType(kVtxTracks), fRequireV0and(0), fTriggerAnalysis(0x0), @@ -57,6 +59,7 @@ AliDielectronEventCuts::AliDielectronEventCuts(const char* name, const char* tit fVtxZmax(0.), fRequireVtx(kFALSE), fMinVtxContributors(0), + fMultITSTPC(kFALSE), fVtxType(kVtxTracks), fRequireV0and(0), fTriggerAnalysis(0x0), @@ -123,6 +126,13 @@ Bool_t AliDielectronEventCuts::IsSelected(TObject* event) if (!v0AND) return kFALSE; } + + if (fMultITSTPC){ + const AliESDVertex *vtxESDTPC=ev->GetPrimaryVertexTPC(); + const AliMultiplicity *multESD = ev->GetMultiplicity(); + if ( vtxESDTPC && multESD && vtxESDTPC->GetNContributors() < (-10.+0.25*multESD->GetNumberOfITSClusters(0)) ) + return kFALSE; + } return kTRUE; } diff --git a/PWG3/dielectron/AliDielectronEventCuts.h b/PWG3/dielectron/AliDielectronEventCuts.h index 227c4d6cc50..5f7fcaef9c7 100644 --- a/PWG3/dielectron/AliDielectronEventCuts.h +++ b/PWG3/dielectron/AliDielectronEventCuts.h @@ -39,7 +39,7 @@ public: void SetRequireVertex(Bool_t req=kTRUE) { fRequireVtx=req; } void SetRequireV0and(UChar_t type=1) { fRequireV0and=type; } void SetMinVtxContributors(Int_t min=1) { fMinVtxContributors=min; } - + void SetCutOnMultipicityITSTPC(Bool_t mult=kTRUE) { fMultITSTPC=mult; } // //Analysis cuts interface // @@ -53,12 +53,13 @@ private: Double_t fVtxZmax; // maximum z vertex position Bool_t fRequireVtx; // require a vertex Int_t fMinVtxContributors; // min number of vertex contributors + Bool_t fMultITSTPC; // if to cut on the ITS TPC multiplicity correlation (Pb-Pb) EVtxType fVtxType; // vertex type UChar_t fRequireV0and; // use V0and triggered events only AliTriggerAnalysis *fTriggerAnalysis; //! trigger analysis class - const AliESDVertex *fkVertex; //! current vertex + const AliESDVertex *fkVertex; //! current vertex AliDielectronEventCuts(const AliDielectronEventCuts &c); AliDielectronEventCuts &operator=(const AliDielectronEventCuts &c); diff --git a/PWG3/dielectron/AliDielectronHelper.cxx b/PWG3/dielectron/AliDielectronHelper.cxx new file mode 100644 index 00000000000..c10edd518ba --- /dev/null +++ b/PWG3/dielectron/AliDielectronHelper.cxx @@ -0,0 +1,113 @@ +/************************************************************************* +* 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 helper functions wrapped in a namespace +// +// +// Authors: +// Jens Wiechula +// + + + + +#include +#include +#include +#include +#include + +#include "AliDielectronHelper.h" + +//_____________________________________________________________________________ +TVectorD* AliDielectronHelper::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) +{ + // + // Make logarithmic binning + // the user has to delete the array afterwards!!! + // + + //check limits + if (xmin<1e-20 || xmax<1e-20){ + Error("AliDielectronHelper::MakeLogBinning","For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!"); + return AliDielectronHelper::MakeLinBinning(nbinsX, xmin, xmax); + } + if (xmaxGetEntries(); + if (nLimits<2){ + Error("AliDielectronHelper::MakeArbitraryBinning","Need at leas 2 bin limits, cannot add the variable"); + delete arr; + return 0x0; + } + + TVectorD *binLimits=new TVectorD(nLimits); + for (Int_t iLim=0; iLim(arr->At(iLim)))->GetString().Atof(); + } + + delete arr; + return binLimits; +} + diff --git a/PWG3/dielectron/AliDielectronHelper.h b/PWG3/dielectron/AliDielectronHelper.h new file mode 100644 index 00000000000..a3bed4dfe92 --- /dev/null +++ b/PWG3/dielectron/AliDielectronHelper.h @@ -0,0 +1,40 @@ +#ifndef ALIDIELECTRONHELPER_H +#define ALIDIELECTRONHELPER_H +/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +/////////////////////////////////////////////////////////////////////////////////////////// +// // +// Dielectron helpers // +// // +// // +// Authors: // +// Jens Wiechula // +// // +/////////////////////////////////////////////////////////////////////////////////////////// + + +#include + +namespace AliDielectronHelper +{ + + + +TVectorD* MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax); +TVectorD* MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax); +TVectorD* MakeArbitraryBinning(const char* bins); + + + + + + + + + + + +} + +#endif diff --git a/PWG3/dielectron/AliDielectronHistos.cxx b/PWG3/dielectron/AliDielectronHistos.cxx index c941d1fe0b2..24e4b344dd9 100644 --- a/PWG3/dielectron/AliDielectronHistos.cxx +++ b/PWG3/dielectron/AliDielectronHistos.cxx @@ -39,8 +39,9 @@ #include #include #include -// #include +#include +#include "AliDielectronHelper.h" #include "AliDielectronHistos.h" @@ -95,25 +96,16 @@ void AliDielectronHistos::UserHistogram(const char* histClass,const char *name, // // Default histogram creation 1D case // - if (!IsHistogramOk(histClass,name)) return; - Double_t *binLimX=0x0; + TVectorD *binLimX=0x0; if (logBinX) { - binLimX=MakeLogBinning(nbinsX, xmin, xmax); + binLimX=AliDielectronHelper::MakeLogBinning(nbinsX, xmin, xmax); } else { - binLimX=MakeLinBinning(nbinsX, xmin, xmax); + binLimX=AliDielectronHelper::MakeLinBinning(nbinsX, xmin, xmax); } - TH1* hist=new TH1F(name,title,nbinsX,binLimX); - - delete [] binLimX; - - Bool_t isReserved=fReservedWords->Contains(histClass); - if (isReserved) - UserHistogramReservedWords(histClass, hist, valTypeX); - else - UserHistogram(histClass, hist, valTypeX); + UserHistogram(histClass,name,title,binLimX,valTypeX); } //_____________________________________________________________________________ @@ -128,31 +120,21 @@ void AliDielectronHistos::UserHistogram(const char* histClass,const char *name, // if (!IsHistogramOk(histClass,name)) return; - Double_t *binLimX=0x0; - Double_t *binLimY=0x0; + TVectorD *binLimX=0x0; + TVectorD *binLimY=0x0; if (logBinX) { - binLimX=MakeLogBinning(nbinsX, xmin, xmax); + binLimX=AliDielectronHelper::MakeLogBinning(nbinsX, xmin, xmax); } else { - binLimX=MakeLinBinning(nbinsX, xmin, xmax); + binLimX=AliDielectronHelper::MakeLinBinning(nbinsX, xmin, xmax); } if (logBinY) { - binLimY=MakeLogBinning(nbinsY, ymin, ymax); + binLimY=AliDielectronHelper::MakeLogBinning(nbinsY, ymin, ymax); } else { - binLimY=MakeLinBinning(nbinsY, ymin, ymax); + binLimY=AliDielectronHelper::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); - else - UserHistogram(histClass, hist, valTypeX+100*valTypeY); + UserHistogram(histClass,name,title,binLimX,binLimY,valTypeX,valTypeY); } @@ -169,39 +151,135 @@ void AliDielectronHistos::UserHistogram(const char* histClass,const char *name, // if (!IsHistogramOk(histClass,name)) return; - Double_t *binLimX=0x0; - Double_t *binLimY=0x0; - Double_t *binLimZ=0x0; + TVectorD *binLimX=0x0; + TVectorD *binLimY=0x0; + TVectorD *binLimZ=0x0; if (logBinX) { - binLimX=MakeLogBinning(nbinsX, xmin, xmax); + binLimX=AliDielectronHelper::MakeLogBinning(nbinsX, xmin, xmax); } else { - binLimX=MakeLinBinning(nbinsX, xmin, xmax); + binLimX=AliDielectronHelper::MakeLinBinning(nbinsX, xmin, xmax); } if (logBinY) { - binLimY=MakeLogBinning(nbinsY, ymin, ymax); + binLimY=AliDielectronHelper::MakeLogBinning(nbinsY, ymin, ymax); } else { - binLimY=MakeLinBinning(nbinsY, ymin, ymax); + binLimY=AliDielectronHelper::MakeLinBinning(nbinsY, ymin, ymax); } if (logBinZ) { - binLimZ=MakeLogBinning(nbinsZ, zmin, zmax); + binLimZ=AliDielectronHelper::MakeLogBinning(nbinsZ, zmin, zmax); } else { - binLimZ=MakeLinBinning(nbinsZ, zmin, zmax); + binLimZ=AliDielectronHelper::MakeLinBinning(nbinsZ, zmin, zmax); + } + + UserHistogram(histClass,name,title,binLimX,binLimY,binLimZ,valTypeX,valTypeY,valTypeZ); +} + +//_____________________________________________________________________________ +void AliDielectronHistos::UserHistogram(const char* histClass,const char *name, const char* title, + const char* binning, + UInt_t valTypeX) +{ + // + // Histogram creation 1D case with arbitraty binning + // + + TVectorD *binLimX=AliDielectronHelper::MakeArbitraryBinning(binning); + UserHistogram(histClass,name,title,binLimX,valTypeX); +} + +//_____________________________________________________________________________ +void AliDielectronHistos::UserHistogram(const char* histClass,const char *name, const char* title, + const TVectorD * const binsX, + UInt_t valTypeX/*=kNoAutoFill*/) +{ + // + // Histogram creation 1D case with arbitraty binning X + // the TVectorD is assumed to be surplus after the creation and will be deleted!!! + // + + Bool_t isOk=kTRUE; + isOk&=IsHistogramOk(histClass,name); + isOk&=(binsX!=0x0); + + if (isOk){ + TH1* hist=new TH1F(name,title,binsX->GetNrows()-1,binsX->GetMatrixArray()); + + Bool_t isReserved=fReservedWords->Contains(histClass); + if (isReserved) + UserHistogramReservedWords(histClass, hist, valTypeX); + else + UserHistogram(histClass, hist, valTypeX); + } + + delete binsX; +} + +//_____________________________________________________________________________ +void AliDielectronHistos::UserHistogram(const char* histClass,const char *name, const char* title, + const TVectorD * const binsX, const TVectorD * const binsY, + UInt_t valTypeX/*=kNoAutoFill*/, UInt_t valTypeY/*=0*/) +{ + // + // Histogram creation 1D case with arbitraty binning X + // the TVectorD is assumed to be surplus after the creation and will be deleted!!! + // + + Bool_t isOk=kTRUE; + isOk&=IsHistogramOk(histClass,name); + isOk&=(binsX!=0x0); + isOk&=(binsY!=0x0); + + if (isOk){ + TH1* hist=new TH2F(name,title, + binsX->GetNrows()-1,binsX->GetMatrixArray(), + binsY->GetNrows()-1,binsY->GetMatrixArray()); + + Bool_t isReserved=fReservedWords->Contains(histClass); + if (isReserved) + UserHistogramReservedWords(histClass, hist, valTypeX+100*valTypeY); + else + UserHistogram(histClass, hist, valTypeX+100*valTypeY); } - TH1* hist=new TH3F(name,title,nbinsX,binLimX,nbinsY,binLimY,nbinsZ,binLimZ); + delete binsX; + delete binsY; + +} + +//_____________________________________________________________________________ +void AliDielectronHistos::UserHistogram(const char* histClass,const char *name, const char* title, + const TVectorD * const binsX, const TVectorD * const binsY, const TVectorD * const binsZ, + UInt_t valTypeX/*=kNoAutoFill*/, UInt_t valTypeY/*=0*/, UInt_t valTypeZ/*=0*/) +{ + // + // Histogram creation 1D case with arbitraty binning X + // the TVectorD is assumed to be surplus after the creation and will be deleted!!! + // + + Bool_t isOk=kTRUE; + isOk&=IsHistogramOk(histClass,name); + isOk&=(binsX!=0x0); + isOk&=(binsY!=0x0); + isOk&=(binsZ!=0x0); + + if (isOk){ + TH1* hist=new TH3F(name,title, + binsX->GetNrows()-1,binsX->GetMatrixArray(), + binsY->GetNrows()-1,binsY->GetMatrixArray(), + binsZ->GetNrows()-1,binsZ->GetMatrixArray()); - delete [] binLimX; - delete [] binLimY; - delete [] binLimZ; + Bool_t isReserved=fReservedWords->Contains(histClass); + if (isReserved) + UserHistogramReservedWords(histClass, hist, valTypeX+100*valTypeY+10000*valTypeZ); + else + UserHistogram(histClass, hist, valTypeX+100*valTypeY+10000*valTypeZ); + } - Bool_t isReserved=fReservedWords->Contains(histClass); - if (isReserved) - UserHistogramReservedWords(histClass, hist, valTypeX+100*valTypeY+10000*valTypeZ); - else - UserHistogram(histClass, hist, valTypeX+100*valTypeY+10000*valTypeZ); + delete binsX; + delete binsY; + delete binsZ; } //_____________________________________________________________________________ @@ -348,7 +426,7 @@ void AliDielectronHistos::FillClass(const char* histClass, Int_t nValues, const // } //_____________________________________________________________________________ -void AliDielectronHistos::UserHistogramReservedWords(const char* histClass, TH1 *hist, UInt_t valTypes) +void AliDielectronHistos::UserHistogramReservedWords(const char* histClass, const TH1 *hist, UInt_t valTypes) { // // Creation of histogram for all pair types @@ -574,6 +652,8 @@ void AliDielectronHistos::SetHistogramList(THashList &list, Bool_t setOwner/*=kT if (setOwner){ list.SetOwner(kFALSE); fHistoList.SetOwner(kTRUE); + } else { + fHistoList.SetOwner(kFALSE); } } @@ -701,44 +781,3 @@ 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 // #include #include +#include class TH1; class TString; @@ -45,6 +46,20 @@ public: 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,const char *name, const char* title, + const char* binning, + UInt_t valTypeX=kNoAutoFill); + + void UserHistogram(const char* histClass,const char *name, const char* title, + const TVectorD * const binsX, + UInt_t valTypeX=kNoAutoFill); + void UserHistogram(const char* histClass,const char *name, const char* title, + const TVectorD * const binsX, const TVectorD * const binsY, + UInt_t valTypeX=kNoAutoFill, UInt_t valTypeY=0); + void UserHistogram(const char* histClass,const char *name, const char* title, + const TVectorD * const binsX, const TVectorD * const binsY, const TVectorD * const binsZ, + UInt_t valTypeX=kNoAutoFill, UInt_t valTypeY=0, UInt_t valTypeZ=0); + void UserHistogram(const char* histClass, TH1* hist, UInt_t valTypes=kNoAutoFill); @@ -90,7 +105,7 @@ private: TList *fList; //! List of list of histograms TString *fReservedWords; //! list of reserved words - void UserHistogramReservedWords(const char* histClass, TH1 *hist, UInt_t valTypes); + void UserHistogramReservedWords(const char* histClass, const TH1 *hist, UInt_t valTypes); void FillClass(THashTable *classTable, Int_t nValues, Double_t *values); void PrintPDF(Option_t* opt); @@ -98,9 +113,6 @@ private: Bool_t IsHistogramOk(const char* classTable, const char* name); - Double_t* MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) const; - Double_t* MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax) const; - enum {kNoAutoFill=1000000000}; AliDielectronHistos(const AliDielectronHistos &hist); diff --git a/PWG3/dielectron/AliDielectronPID.cxx b/PWG3/dielectron/AliDielectronPID.cxx index b4b1cebd47f..3381b1af4a3 100644 --- a/PWG3/dielectron/AliDielectronPID.cxx +++ b/PWG3/dielectron/AliDielectronPID.cxx @@ -31,9 +31,11 @@ Detailed description #include #include +#include #include #include #include +#include #include "AliDielectronVarManager.h" @@ -204,9 +206,6 @@ Bool_t AliDielectronPID::IsSelected(TObject* track) // test momentum range. In case pMin==pMax use all momenta if ( (TMath::Abs(pMin-pMax)>1e-20) && (mom<=pMin || mom>pMax) ) continue; - // test if we are supposed to use a function for the cut - if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(mom); - if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(mom); switch (fDetType[icut]){ case kITS: @@ -229,7 +228,7 @@ Bool_t AliDielectronPID::IsSelected(TObject* track) } //______________________________________________ -Bool_t AliDielectronPID::IsSelectedITS(AliVTrack * const part, Int_t icut) const +Bool_t AliDielectronPID::IsSelectedITS(AliVTrack * const part, Int_t icut) { // // ITS part of the PID check @@ -239,6 +238,8 @@ Bool_t AliDielectronPID::IsSelectedITS(AliVTrack * const part, Int_t icut) const if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliESDtrack::kITSpid)) return kFALSE; if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliESDtrack::kITSpid)) return kTRUE; + + Double_t mom=part->P(); if (part->IsA()==AliESDtrack::Class()){ // ESD case in case the PID bit is not set, don't use this track! @@ -250,12 +251,17 @@ Bool_t AliDielectronPID::IsSelectedITS(AliVTrack * const part, Int_t icut) const AliAODTrack *track=static_cast(part); numberOfSigmas=fAODpidUtil->NumberOfSigmasITS(track, fPartType[icut]); } + + // test if we are supposed to use a function for the cut + if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(mom); + if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(mom); + Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut]; return selected; } //______________________________________________ -Bool_t AliDielectronPID::IsSelectedTPC(AliVTrack * const part, Int_t icut) const +Bool_t AliDielectronPID::IsSelectedTPC(AliVTrack * const part, Int_t icut) { // // TPC part of the PID check @@ -266,25 +272,36 @@ Bool_t AliDielectronPID::IsSelectedTPC(AliVTrack * const part, Int_t icut) const if (fRequirePIDbit[icut]==AliDielectronPID::kRequire&&!(part->GetStatus()&AliESDtrack::kTPCpid)) return kFALSE; if (fRequirePIDbit[icut]==AliDielectronPID::kIfAvailable&&!(part->GetStatus()&AliESDtrack::kTPCpid)) return kTRUE; + Double_t mom=part->P(); + if (part->IsA()==AliESDtrack::Class()){ // ESD case in case the PID bit is not set, don't use this track! AliESDtrack *track=static_cast(part); + const AliExternalTrackParam *in = track->GetInnerParam(); + if (in) mom = in->GetP(); numberOfSigmas=fESDpid->NumberOfSigmasTPC(track, fPartType[icut]); }else if(part->IsA()==AliAODTrack::Class()){ // AOD case // FIXME: Is there a place to check whether the PID is was set in ESD??? AliAODTrack *track=static_cast(part); + const AliAODPid *pidObj = track->GetDetPid(); + if (pidObj) mom = pidObj->GetTPCmomentum(); numberOfSigmas=fAODpidUtil->NumberOfSigmasTPC(track, fPartType[icut]); } if (fPartType[icut]==AliPID::kElectron){ numberOfSigmas-=fgCorr; } + + // test if we are supposed to use a function for the cut + if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(mom); + if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(mom); + Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut]; return selected; } //______________________________________________ -Bool_t AliDielectronPID::IsSelectedTRD(AliVTrack * const /*part*/, Int_t /*icut*/) const +Bool_t AliDielectronPID::IsSelectedTRD(AliVTrack * const /*part*/, Int_t /*icut*/) { // // TRD part of the pid check @@ -293,7 +310,7 @@ Bool_t AliDielectronPID::IsSelectedTRD(AliVTrack * const /*part*/, Int_t /*icut* } //______________________________________________ -Bool_t AliDielectronPID::IsSelectedTOF(AliVTrack * const part, Int_t icut) const +Bool_t AliDielectronPID::IsSelectedTOF(AliVTrack * const part, Int_t icut) { // // TOF part of the PID check diff --git a/PWG3/dielectron/AliDielectronPID.h b/PWG3/dielectron/AliDielectronPID.h index 3f959e67ca6..8ffe33cec52 100644 --- a/PWG3/dielectron/AliDielectronPID.h +++ b/PWG3/dielectron/AliDielectronPID.h @@ -89,10 +89,10 @@ private: static Double_t fgCorr; //!correction value for current run. Set if fgFitCorr is set and SetCorrVal(run) // was called - Bool_t IsSelectedITS(AliVTrack * const part, Int_t icut) const; - Bool_t IsSelectedTPC(AliVTrack * const part, Int_t icut) const; - Bool_t IsSelectedTRD(AliVTrack * const part, Int_t icut) const; - Bool_t IsSelectedTOF(AliVTrack * const part, Int_t icut) const; + Bool_t IsSelectedITS(AliVTrack * const part, Int_t icut); + Bool_t IsSelectedTPC(AliVTrack * const part, Int_t icut); + Bool_t IsSelectedTRD(AliVTrack * const part, Int_t icut); + Bool_t IsSelectedTOF(AliVTrack * const part, Int_t icut); AliDielectronPID(const AliDielectronPID &c); AliDielectronPID &operator=(const AliDielectronPID &c); diff --git a/PWG3/dielectron/AliDielectronSignalBase.cxx b/PWG3/dielectron/AliDielectronSignalBase.cxx index c5bc8334094..da3cc998f95 100644 --- a/PWG3/dielectron/AliDielectronSignalBase.cxx +++ b/PWG3/dielectron/AliDielectronSignalBase.cxx @@ -52,7 +52,7 @@ AliDielectronSignalBase::AliDielectronSignalBase() : fMethod(kLikeSign), fScaleMin(0.), fScaleMax(0.), - fScaleFactor(0.), + fScaleFactor(1.), fProcessed(kFALSE) { // @@ -78,7 +78,7 @@ AliDielectronSignalBase::AliDielectronSignalBase(const char* name, const char* t fMethod(kLikeSign), fScaleMin(0.), fScaleMax(0.), - fScaleFactor(0.), + fScaleFactor(1.), fProcessed(kFALSE) { // @@ -153,6 +153,15 @@ Double_t AliDielectronSignalBase::ScaleHistograms(TH1* histRaw, TH1* histBackgro // // scale histBackground to match the integral of histRaw in the interval intMin, intMax // + + //protect using over and underflow bins in normalisation calculation + if (intMinGetXaxis()->GetXmin()) intMin=histRaw->GetXaxis()->GetXmin(); + if (intMinGetXaxis()->GetXmin()) intMin=histBackground->GetXaxis()->GetXmin(); + + if (intMax>histRaw->GetXaxis()->GetXmax()) + intMax=histRaw->GetXaxis()->GetXmax()-histRaw->GetBinWidth(histRaw->GetNbinsX())/2.; + if (intMax>histBackground->GetXaxis()->GetXmax()) + intMax=histBackground->GetXaxis()->GetXmax()-histBackground->GetBinWidth(histBackground->GetNbinsX())/2.; Double_t intRaw = histRaw->Integral(histRaw->FindBin(intMin),histRaw->FindBin(intMax)); Double_t intBack = histBackground->Integral(histBackground->FindBin(intMin),histBackground->FindBin(intMax)); diff --git a/PWG3/dielectron/AliDielectronTrackRotator.cxx b/PWG3/dielectron/AliDielectronTrackRotator.cxx index 8fd6ca8844c..544abaa77a1 100644 --- a/PWG3/dielectron/AliDielectronTrackRotator.cxx +++ b/PWG3/dielectron/AliDielectronTrackRotator.cxx @@ -145,9 +145,9 @@ Bool_t AliDielectronTrackRotator::RotateTracks() // Actual track rotation // Find out particle type and perform the rotation // - + const AliVTrack *trackP=dynamic_cast(fkArrTracksP->UncheckedAt(fCurrentTackP)); - const AliVTrack *trackN=dynamic_cast(fkArrTracksP->UncheckedAt(fCurrentTackN)); + const AliVTrack *trackN=dynamic_cast(fkArrTracksN->UncheckedAt(fCurrentTackN)); if (!trackP||!trackN) return kFALSE; diff --git a/PWG3/dielectron/AliDielectronVarManager.cxx b/PWG3/dielectron/AliDielectronVarManager.cxx index 0ef4c84a2d1..a4da04eef55 100644 --- a/PWG3/dielectron/AliDielectronVarManager.cxx +++ b/PWG3/dielectron/AliDielectronVarManager.cxx @@ -90,6 +90,8 @@ const char* AliDielectronVarManager::fgkParticleNames[AliDielectronVarManager::k "TOF_nSigma_Muons", "TOF_nSigma_Kaons", "TOF_nSigma_Protons", + + "KinkIndex0", // "Chi2NDF", "DecayLength", @@ -115,6 +117,7 @@ const char* AliDielectronVarManager::fgkParticleNames[AliDielectronVarManager::k "ZRes", "NTrk", "Tracks", + "Centrality", "Nevents" }; diff --git a/PWG3/dielectron/AliDielectronVarManager.h b/PWG3/dielectron/AliDielectronVarManager.h index 01c097479a5..c86319afa29 100644 --- a/PWG3/dielectron/AliDielectronVarManager.h +++ b/PWG3/dielectron/AliDielectronVarManager.h @@ -42,6 +42,7 @@ #include #include +#include #include #include @@ -124,6 +125,8 @@ public: kTOFnSigmaKao, // number of sigmas to the kaon line in the TOF kTOFnSigmaPro, // number of sigmas to the proton line in the TOF + kKinkIndex0, // kink index 0 + kParticleMax, // // TODO: kRNClusters ?? // AliDielectronPair specific variables @@ -155,6 +158,7 @@ public: kZRes, // primary vertex z-resolution kNTrk, // number of tracks (or tracklets) kTracks, // ESD tracks + kCentrality, // event centrality fraction kNevents, // event counter kNMaxValues // // TODO: (for A+A) ZDCEnergy, impact parameter, Iflag?? @@ -268,7 +272,7 @@ inline void AliDielectronVarManager::FillVarESDtrack(const AliESDtrack *particle Double_t tpcNcls=particle->GetTPCNcls(); values[AliDielectronVarManager::kNclsITS] = particle->GetNcls(0); // TODO: get rid of the plain numbers values[AliDielectronVarManager::kNclsTPC] = tpcNcls; // TODO: get rid of the plain numbers - values[AliDielectronVarManager::kNclsTPC] = particle->GetTPCNclsIter1(); // TODO: get rid of the plain numbers + values[AliDielectronVarManager::kNclsTPCiter1] = particle->GetTPCNclsIter1(); // TODO: get rid of the plain numbers values[AliDielectronVarManager::kNFclsTPC] = particle->GetTPCNclsF(); values[AliDielectronVarManager::kTPCsignalN] = particle->GetTPCsignalN(); values[AliDielectronVarManager::kNclsTRD] = particle->GetNcls(2); // TODO: get rid of the plain numbers @@ -283,6 +287,8 @@ inline void AliDielectronVarManager::FillVarESDtrack(const AliESDtrack *particle particle->GetTRDpid(pidProbs); values[AliDielectronVarManager::kTRDprobEle] = pidProbs[AliPID::kElectron]; values[AliDielectronVarManager::kTRDprobPio] = pidProbs[AliPID::kPion]; + + values[AliDielectronVarManager::kKinkIndex0] = particle->GetKinkIndex(0); Float_t impactParXY, impactParZ; particle->GetImpactParameters(impactParXY, impactParZ); @@ -602,12 +608,17 @@ inline void AliDielectronVarManager::FillVarESDEvent(const AliESDEvent *event, D // Fill common AliVEvent interface information FillVarVEvent(event, values); - + + Double_t centralityF=-1; + AliESDCentrality *esdCentrality = const_cast(event)->GetCentrality(); + if (esdCentrality) centralityF = esdCentrality->GetCentralityPercentile("V0M"); + // Fill AliESDEvent interface specific information const AliESDVertex *primVtx = event->GetPrimaryVertex(); values[AliDielectronVarManager::kXRes] = primVtx->GetXRes(); values[AliDielectronVarManager::kYRes] = primVtx->GetYRes(); values[AliDielectronVarManager::kZRes] = primVtx->GetZRes(); + values[AliDielectronVarManager::kCentrality] = centralityF; } inline void AliDielectronVarManager::FillVarAODEvent(const AliAODEvent *event, Double_t * const values) diff --git a/PWG3/dielectron/macros/AddTaskJPSIFilter.C b/PWG3/dielectron/macros/AddTaskJPSIFilter.C index 9f28d9fe828..ff13f6d608f 100644 --- a/PWG3/dielectron/macros/AddTaskJPSIFilter.C +++ b/PWG3/dielectron/macros/AddTaskJPSIFilter.C @@ -34,6 +34,7 @@ AliAnalysisTask *AddTaskJPSIFilter(){ aodHandler->SetNeedsCaloClustersBranchReplication(); //aodHandler->SetNeedsMCParticlesBranchReplication(); aodHandler->SetNeedsDimuonsBranchReplication(); + if(hasMC) aodHandler->SetNeedsMCParticlesBranchReplication(); } //Create task and add it to the analysis manager diff --git a/PWG3/dielectron/macros/AddTaskJPSIFilterPbPb.C b/PWG3/dielectron/macros/AddTaskJPSIFilterPbPb.C new file mode 100644 index 00000000000..36a26888e6e --- /dev/null +++ b/PWG3/dielectron/macros/AddTaskJPSIFilterPbPb.C @@ -0,0 +1,77 @@ +AliAnalysisTask *AddTaskJPSIFilter(){ + //get the current analysis manager + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + Error("AddTaskJPSIFilter", "No analysis manager found."); + return 0; + } + + //check for output aod handler + if (!mgr->GetOutputEventHandler()||mgr->GetOutputEventHandler()->IsA()!=AliAODHandler::Class()) { + Warning("AddTaskJPSIFilter","No AOD output handler available. Not adding the task!"); + return 0; + } + + //Do we have an MC handler? + Bool_t hasMC=(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()!=0x0); + + //Do we run on AOD? + Bool_t isAOD=mgr->GetInputEventHandler()->IsA()==AliAODInputHandler::Class(); + + if(isAOD) { + //add options to AliAODHandler to duplicate input event + AliAODHandler *aodHandler = (AliAODHandler*)mgr->GetOutputEventHandler(); + aodHandler->SetCreateNonStandardAOD(); + aodHandler->SetNeedsHeaderReplication(); + aodHandler->SetNeedsTracksBranchReplication(); + aodHandler->SetNeedsVerticesBranchReplication(); + aodHandler->SetNeedsV0sBranchReplication(); + aodHandler->SetNeedsCascadesBranchReplication(); + aodHandler->SetNeedsTrackletsBranchReplication(); + aodHandler->SetNeedsPMDClustersBranchReplication(); + aodHandler->SetNeedsJetsBranchReplication(); + aodHandler->SetNeedsFMDClustersBranchReplication(); + aodHandler->SetNeedsCaloClustersBranchReplication(); + //aodHandler->SetNeedsMCParticlesBranchReplication(); + aodHandler->SetNeedsDimuonsBranchReplication(); + if(hasMC) aodHandler->SetNeedsMCParticlesBranchReplication(); + } + + //Create task and add it to the analysis manager + AliAnalysisTaskDielectronFilter *task=new AliAnalysisTaskDielectronFilter("jpsi_DielectronFilter"); + + gROOT->LoadMacro("$ALICE_ROOT/PWG3/dielectron/macros/ConfigJpsi2eeFilterPbPb.C"); + AliDielectron *jpsi=ConfigJpsi2eeFilter(isAOD); + if (!hasMC) task->UsePhysicsSelection(); + task->SetDielectron(jpsi); + mgr->AddTask(task); + + //---------------------- + //create data containers + //---------------------- + + + TString containerName = mgr->GetCommonFileName(); + containerName += ":PWG3_dielectronFilter"; + + //create output container + + AliAnalysisDataContainer *cOutputHist1 = + mgr->CreateContainer("jpsi_FilterQA", + THashList::Class(), + AliAnalysisManager::kOutputContainer, + containerName.Data()); + + AliAnalysisDataContainer *cOutputHist2 = + mgr->CreateContainer("jpsi_FilterEventStat", + TH1D::Class(), + AliAnalysisManager::kOutputContainer, + containerName.Data()); + + + mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer()); + mgr->ConnectOutput(task, 1, cOutputHist1); + mgr->ConnectOutput(task, 2, cOutputHist2); + + return task; +} diff --git a/PWG3/dielectron/macros/ConfigJpsi2eeFilterPbPb.C b/PWG3/dielectron/macros/ConfigJpsi2eeFilterPbPb.C new file mode 100644 index 00000000000..47e7358703d --- /dev/null +++ b/PWG3/dielectron/macros/ConfigJpsi2eeFilterPbPb.C @@ -0,0 +1,177 @@ +void SetupTrackCutsDieleFilter(AliDielectron *diele, Bool_t isAOD); +void SetupPairCutsDieleFilter(AliDielectron *diele, Bool_t isAOD); +void SetupEventCutsDieleFilter(AliDielectron *diele, Int_t cutDefinition); + +void InitHistogramsDieleFilter(AliDielectron *diele); + +AliESDtrackCuts *SetupESDtrackCutsDieleFilter(); + + +AliDielectron* ConfigJpsi2eeFilter(Bool_t isAOD=kFALSE) +{ + // + // Setup the instance of AliDielectron + // + + // create the actual framework object + TString name="trackQ+Pt>1.+652&&|Y|<0.9"; + AliDielectron *diele = new AliDielectron(Form("%s",name.Data()), + Form("Track cuts: %s",name.Data())); + + // cut setup + SetupEventCutsDieleFilter(diele); + + SetupTrackCutsDieleFilter(diele, isAOD); + SetupPairCutsDieleFilter(diele, isAOD); + + // + // QA histogram setup + // + InitHistogramsDieleFilter(diele, isAOD); + + return diele; +} + +//______________________________________________________________________________________ +void SetupEventCutsDieleFilter(AliDielectron *diele) +{ + // + // Setup the event cuts + // + AliDielectronVarCuts *vtxZ = new AliDielectronVarCuts("vtxZ","Vertex z cut"); + vtxZ->AddCut(AliDielectronVarManager::kZvPrim,-15.,15.); + diele->GetEventFilter().AddCuts(vtxZ); +} + +//______________________________________________________________________________________ +void SetupTrackCutsDieleFilter(AliDielectron *diele, Bool_t isAOD) +{ + // + // Setup the track cuts + // + + //ESD quality cuts DielectronTrackCuts + if (!isAOD) { + diele->GetTrackFilter().AddCuts(SetupESDtrackCutsDieleFilter()); + } else { + AliDielectronTrackCuts *trackCuts=new AliDielectronTrackCuts("trackCuts","trackCuts"); + trackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny); + trackCuts->SetRequireTPCRefit(kTRUE); + trackCuts->SetRequireITSRefit(kTRUE); + diele->GetTrackFilter().AddCuts(trackCuts); + } + + //Pt cut + AliDielectronVarCuts *pt = new AliDielectronVarCuts("Pt>.5+601. && 65AddCut(AliDielectronVarManager::kPt,1.,1e30); + pt->AddCut(AliDielectronVarManager::kTPCsignal,65.,90.); + + if (isAOD){ + // TPC #clusteres cut + pt->AddCut(AliDielectronVarManager::kNclsTPC,90.,160.); +// pt->AddCut(AliDielectronVarManager::kEta,-0.88,0.88); + //TODO: DCA cuts to be investigated!!! +// pt->AddCut(AliDielectronVarManager::kImpactParXY,-1.,1.); +// pt->AddCut(AliDielectronVarManager::kImpactParZ,-3.,3.); + } + + diele->GetTrackFilter().AddCuts(pt); +} + +//______________________________________________________________________________________ +void SetupPairCutsDieleFilter(AliDielectron *diele, Bool_t isAOD) +{ + // + // Setup the pair cuts + // + + + //Invarian mass selection + AliDielectronVarCuts *invMassCut=new AliDielectronVarCuts("InvMassY","2AddCut(AliDielectronVarManager::kM,2.,1e30); + invMassCut->AddCut(AliDielectronVarManager::kY,-0.9,0.9); + +// invMassCut->AddCut(AliDielectronVarManager::kPairType,1.); + diele->GetPairFilter().AddCuts(invMassCut); + +} + +//______________________________________________________________________________________ +AliESDtrackCuts *SetupESDtrackCutsDieleFilter() +{ + // + // Setup default AliESDtrackCuts + // + AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts; + + esdTrackCuts->SetEtaRange( -0.9 , 0.9 ); + + esdTrackCuts->SetMaxDCAToVertexZ(3.0); + esdTrackCuts->SetMaxDCAToVertexXY(1.0); + esdTrackCuts->SetRequireTPCRefit(kTRUE); + esdTrackCuts->SetRequireITSRefit(kTRUE); + esdTrackCuts->SetAcceptKinkDaughters(kFALSE); + esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny); + + esdTrackCuts->SetMinNClustersTPC(100); + esdTrackCuts->SetMaxChi2PerClusterTPC(4); + + return esdTrackCuts; +} + + +//______________________________________________________________________________________ +void InitHistogramsDieleFilter(AliDielectron *diele, Bool_t isAOD) +{ + // + // Initialise the histograms + // + +//Setup histogram classes + AliDielectronHistos *histos= + new AliDielectronHistos(diele->GetName(), + diele->GetTitle()); + + //Initialise histogram classes + histos->SetReservedWords("Track;Pair"); + + histos->AddClass("Event"); + //Track classes + //to fill also track info from 2nd event loop until 2 + for (Int_t i=0; i<2; ++i){ + histos->AddClass(Form("Track_%s",AliDielectron::TrackClassName(i))); + } + + //Pair classes + // to fill also mixed event histograms loop until 10 + for (Int_t i=0; i<3; ++i){ + histos->AddClass(Form("Pair_%s",AliDielectron::PairClassName(i))); + } + + //add histograms to event class + histos->UserHistogram("Event","VtxZ","Vertex Z;Z[cm]",300,-15.,15.,AliDielectronVarManager::kZvPrim); + + //add histograms to Track classes + histos->UserHistogram("Track","Pt","Pt;Pt [GeV];#tracks",400,0,20.,AliDielectronVarManager::kPt); + histos->UserHistogram("Track","TPCnCls","Number of Clusters TPC;TPC number clusteres;#tracks",160,-0.5,159.5,AliDielectronVarManager::kNclsTPC); + + histos->UserHistogram("Track","dXY","dXY;dXY [cm];#tracks",500,-1.,1.,AliDielectronVarManager::kImpactParXY); + histos->UserHistogram("Track","dZ","dZ;dZ [cm];#tracks",600,-3.,3.,AliDielectronVarManager::kImpactParZ); + histos->UserHistogram("Track","Eta_Phi","Eta Phi Map; Eta; Phi;#tracks", + 100,-1,1,144,0,6.285,AliDielectronVarManager::kEta,AliDielectronVarManager::kPhi); + + histos->UserHistogram("Track","dEdx_P","dEdx;P [GeV];TPC signal (arb units);#tracks", + 200,0.2,20.,100,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", + 200,0.2,20.,100,-10.,10.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaEle,kTRUE); + + //add histograms to Pair classes + histos->UserHistogram("Pair","InvMass","Inv.Mass;Inv. Mass [GeV];#pairs", + 201,-.01,4.01,AliDielectronVarManager::kM); + histos->UserHistogram("Pair","Rapidity","Rapidity;Rapidity;#pairs", + 100,-1.,1.,AliDielectronVarManager::kY); + histos->UserHistogram("Pair","OpeningAngle","Opening angle;angle", + 100,0.,3.15,AliDielectronVarManager::kOpeningAngle); + + diele->SetHistogramManager(histos); +} diff --git a/PWG3/dielectron/macros/ExtractEfficiencies.C b/PWG3/dielectron/macros/ExtractEfficiencies.C new file mode 100644 index 00000000000..19d992cbf44 --- /dev/null +++ b/PWG3/dielectron/macros/ExtractEfficiencies.C @@ -0,0 +1,1284 @@ +// author: Ionut Cristian Arsene +// Date: 01/09/2010 + +#include +#include +using namespace std; + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "AliCFContainer.h" +#include "AliDielectronCFdraw.h" + +/* +This macro makes projections and saves histograms from a list of CF containers +generated with the dielectron package. These histograms can later be used to +calculate efficiencies. +To use it, the following modifications are needed: +1) Modify the global variables listed below according to your needs. +2) Make projections by applying as many cut sets as needed on the CF containers. + Call the FillHistograms() after each cut set. +3) To extract efficiencies use the ExtractEfficiencies() function where one needs to + specify the indexes for the nominator and denominator histograms. The + indexes are based on the CF step number, cut set number and histogram number + as defined in the global variables below. +*/ + +// The CF variable indexes ------------------------------------------------ +enum Variables { // create an enumeration item for every variable from your CF container + kNothing = -1, // kNothing should be always here + kPt = 0, + kY, + kPairType, + kThetaCS, + kThetaHE, + kM, + kLeg1_Eta, + // kLeg1_TPC_nSigma_Electrons, + kLeg1_Pt, + // kLeg1_P, + kLeg1_NclsTPC, + kLeg2_Eta, + // kLeg2_TPC_nSigma_Electrons, + kLeg2_Pt, + // kLeg2_P, + kLeg2_NclsTPC, + kNVariables // kNVariables should be always here! +}; +const Char_t* gkVarNames[kNVariables] = { // variable names to be put on histograms axes and titles + "p_{T} [GeV/c]", + "y", + "Pair type (0=++; 1=+-; 2=--)", + "cos #theta^{*}_{CS}", + "cos #theta^{*}_{HE}", + "M [GeV/c^{2}]", + "#eta^{leg1}", + // "TPC n #sigma electrons (leg1)", + "p_{T}^{leg1} [GeV/c]", + // "P^{leg1} [GeV/c]", + "# TPC clusters (leg1)", + "#eta^{leg2}", + // "TPC n #sigma electrons (leg2)", + "p_{T}^{leg2} [GeV/c]", + // "P^{leg2} [GeV/c]", + "# TPC clusters (leg2)" +}; +Int_t gNbins[kNVariables]; // number of bins for every variable --> filled automatically +Double_t* gBinLimits[kNVariables]; // bin limits for every variable --> filled automatically +// ------------------------------------------------------------------------ + +// Put here all the CF steps of interest ---------------------------------- +enum Steps { // step indexes in the CF containers to be analyzed + kPureMC = 0, + kESDSPDany = 2, + kESDSPDfirst = 4, + kFullSPDany = 6, + kFullSPDfirst = 8, + kNSteps = 5 // total number of steps (the number of steps above) +}; +const Int_t gkStepNumbers[kNSteps] = { // array with step indexes (all from the enumeration above) + kPureMC, + kESDSPDany, kESDSPDfirst, + kFullSPDany, kFullSPDfirst +}; +const Char_t* gkStepNames[kNSteps][2] = {// names for each CF step + {"PureMC", "Pure MC"}, // NOTE: short names go to histo names, long names go to titles + {"ESDSPDany", "ESD track cuts, SPD any, TPCnclus>90"}, + {"ESDSPDfirst", "ESD track cuts, SPD first, TPCnclus>90"}, + {"FullSPDany", "All track cuts (with SPD any) and TPC-PID"}, + {"FullSPDfirst", "All track cuts (with SPD first) and TPC-PID"} +}; +//------------------------------------------------------------------------ + +// Put here info about the cut sets for which projections will be made --- +const Int_t gkNCutSets = 18; // number of cut sets for which histos will be filled +const Char_t* gkCutSetNames[gkNCutSets][2] = { // short and long names for all the cut sets + // baseline + {"Ycut", "|y_{J/#Psi}|<0.88"}, + {"YcutPt1", "|y_{J/#Psi}|<0.88 & 0.01.0 GeV/c"}, + {"YcutPt1LegsEtaPt1", "|y_{J/#Psi}|<0.88 & 0.01.0 GeV/c"}, + {"YcutPt2LegsEtaPt1", "|y_{J/#Psi}|<0.88 & 0.81.0 GeV/c"}, + {"YcutPt3LegsEtaPt1", "|y_{J/#Psi}|<0.88 & 1.41.0 GeV/c"}, + {"YcutPt4LegsEtaPt1", "|y_{J/#Psi}|<0.88 & 2.81.0 GeV/c"}, + {"YcutPt5LegsEtaPt1", "|y_{J/#Psi}|<0.88 & 5.01.0 GeV/c"}, + // track quality + {"YcutLegsEtaPt1TPC90", "|y_{J/#Psi}|<0.88 & |#eta_{legs}|<0.88 & p_{Tlegs}>1.0 GeV/c TPCnclus>90"}, + {"YcutPt1LegsEtaPt1TPC90", "|y_{J/#Psi}|<0.88 & 0.01.0 GeV/c TPCnclus>90"}, + {"YcutPt2LegsEtaPt1TPC90", "|y_{J/#Psi}|<0.88 & 0.81.0 GeV/c TPCnclus>90"}, + {"YcutPt3LegsEtaPt1TPC90", "|y_{J/#Psi}|<0.88 & 1.41.0 GeV/c TPCnclus>90"}, + {"YcutPt4LegsEtaPt1TPC90", "|y_{J/#Psi}|<0.88 & 2.81.0 GeV/c TPCnclus>90"}, + {"YcutPt5LegsEtaPt1TPC90", "|y_{J/#Psi}|<0.88 & 5.01.0 GeV/c TPCnclus>90"} + // track quality + PID + // {"YcutLegsEtaPt1TPC90PID", "|y_{J/#Psi}|<0.88 & |#eta_{legs}|<0.88 & p_{Tlegs}>1.0 GeV/c TPCnclus>90 & TPC pid"}, + // {"YcutPt1LegsEtaPt1TPC90PID", "|y_{J/#Psi}|<0.88 & 0.01.0 GeV/c TPCnclus>90 & TPC pid"}, + // {"YcutPt2LegsEtaPt1TPC90PID", "|y_{J/#Psi}|<0.88 & 0.81.0 GeV/c TPCnclus>90 & TPC pid"}, + // {"YcutPt3LegsEtaPt1TPC90PID", "|y_{J/#Psi}|<0.88 & 1.41.0 GeV/c TPCnclus>90 & TPC pid"}, + // {"YcutPt4LegsEtaPt1TPC90PID", "|y_{J/#Psi}|<0.88 & 2.81.0 GeV/c TPCnclus>90 & TPC pid"}, + // {"YcutPt5LegsEtaPt1TPC90PID", "|y_{J/#Psi}|<0.88 & 5.01.0 GeV/c TPCnclus>90 & TPC pid"} +}; +// ----------------------------------------------------------------------- + +// Put here info about the histograms to be filled ----------------------- +const Int_t gkNhistos = 6; // how many histograms for every (step,cut set) combination +const Char_t* gkHistoNames[gkNhistos][2] = { // short and long names of the histograms + {"pt","p_{T}(J/#Psi)"}, // NOTE: short names go to the histo name, long name goes to title + {"y","y(J/#Psi)"}, + {"pty","p_{T} vs y(J/#Psi)"}, + // {"m","e^{+}e^{-} invariant mass"}, + {"ThetaCS","cos #theta^{*}_{CS}"}, + {"ThetaHE","cos #theta^{*}_{HE}"}, + {"Minv", "Invariant mass"} +}; +const Int_t gkDims[gkNhistos][4] = { // dimensions and variables for histograms +// ndim xVar yVar zVar + {1, kPt, kNothing, kNothing}, // pt dependence + {1, kY, kNothing, kNothing}, // y dependence + {2, kY, kPt, kNothing}, // pt,y dependence + // {1, kM, kNothing, kNothing}, // inv. mass dependence + {1, kThetaCS, kNothing, kNothing}, // cos theta* CS dependence + {1, kThetaHE, kNothing, kNothing}, // cos theta* HE dependence + {1, kM, kNothing, kNothing} // invariant mass +}; +// ----------------------------------------------------------------------- + +// ******************************************************************************** +// Define here all the efficiencies you want (if any) +// Efficiency = (nominator histogram)/(denominator histogram) +// A histogram is defined by its step,cut set, and number defined according the +// global variables above +// ******************************************************************************** +const Int_t gkNeffs = 210; +const Int_t gkEffs[gkNeffs][6] = { + //nominator: Step Cut Histo | denominator: Step Cut Histo comment + // full corrections, pt dependence + { 3, 12, 0, 0, 0, 0 }, // full correction, SPD any + { 4, 12, 0, 0, 0, 0 }, // full correction, SPD first + { 3, 13, 0, 0, 1, 0 }, // full correction, SPD any, 0.090, pt dependence + { 1, 12, 0, 0, 6, 0 }, // tracking correction, full pt range + { 1, 13, 0, 0, 7, 0 }, // tracking correction, 0.090, y dependence + { 1, 12, 1, 0, 6, 1 }, // tracking correction, full pt range + { 1, 13, 1, 0, 7, 1 }, // tracking correction, 0.090, pt,y dependence + { 1, 12, 2, 0, 6, 2 }, // tracking correction, full pt range + { 1, 13, 2, 0, 7, 2 }, // tracking correction, 0.090, cos Theta* CS dependence + { 1, 12, 3, 0, 6, 3 }, // tracking correction, full pt range + { 1, 13, 3, 0, 7, 3 }, // tracking correction, 0.090, cos Theta* HE dependence + { 1, 12, 4, 0, 6, 4 }, // tracking correction, full pt range + { 1, 13, 4, 0, 7, 4 }, // tracking correction, 0.090, pt dependence + { 2, 12, 0, 0, 6, 0 }, // tracking correction, full pt range + { 2, 13, 0, 0, 7, 0 }, // tracking correction, 0.090, y dependence + { 2, 12, 1, 0, 6, 1 }, // tracking correction, full pt range + { 2, 13, 1, 0, 7, 1 }, // tracking correction, 0.090, pt,y dependence + { 2, 12, 2, 0, 6, 2 }, // tracking correction, full pt range + { 2, 13, 2, 0, 7, 2 }, // tracking correction, 0.090, cos Theta*CS dependence + { 2, 12, 3, 0, 6, 3 }, // tracking correction, full pt range + { 2, 13, 3, 0, 7, 3 }, // tracking correction, 0.090, cos Theta*HE dependence + { 2, 12, 4, 0, 6, 4 }, // tracking correction, full pt range + { 2, 13, 4, 0, 7, 4 }, // tracking correction, 0.01.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt1_pt", "Kinematic correction vs pt, 0.01.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt2_pt", "Kinematic correction vs pt, 0.81.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt3_pt", "Kinematic correction vs pt, 1.41.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt4_pt", "Kinematic correction vs pt, 2.81.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt5_pt", "Kinematic correction vs pt, 5.01.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + // acceptance corrections, y dependence + {"accCorrection_y", "Kinematic correction vs y, (J/#Psi in |y|<0.88 & |#eta_{legs}|<0.88 & p_{Tlegs}>1.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt1_y", "Kinematic correction vs y, 0.01.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt2_y", "Kinematic correction vs y, 0.81.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt3_y", "Kinematic correction vs y, 1.41.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt4_y", "Kinematic correction vs y, 2.81.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt5_y", "Kinematic correction vs y, 5.01.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + // acceptance corrections, pt-y dependence + {"accCorrection_pty", "Kinematic correction vs pt-y, (J/#Psi in |y|<0.88 & |#eta_{legs}|<0.88 & p_{Tlegs}>1.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt1_pty", "Kinematic correction vs pt-y, 0.01.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt2_pty", "Kinematic correction vs pt-y, 0.81.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt3_pty", "Kinematic correction vs pt-y, 1.41.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt4_pty", "Kinematic correction vs pt-y, 2.81.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt5_pty", "Kinematic correction vs pt-y, 5.01.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + // acceptance corrections, cos Theta*CS dependence + {"accCorrection_ThetaCS", "Kinematic correction vs cos #theta^{*}_{CS}, (J/#Psi in |y|<0.88 & |#eta_{legs}|<0.88 & p_{Tlegs}>1.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt1_ThetaCS", "Kinematic correction vs cos #theta^{*}_{CS}, 0.01.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt2_ThetaCS", "Kinematic correction vs cos #theta^{*}_{CS}, 0.81.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt3_ThetaCS", "Kinematic correction vs cos #theta^{*}_{CS}, 1.41.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt4_ThetaCS", "Kinematic correction vs cos #theta^{*}_{CS}, 2.81.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt5_ThetaCS", "Kinematic correction vs cos #theta^{*}_{CS}, 5.01.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + // acceptance corrections, cos Theta*HE dependence + {"accCorrection_ThetaHE", "Kinematic correction vs cos #theta^{*}_{HE}, (J/#Psi in |y|<0.88 & |#eta_{legs}|<0.88 & p_{Tlegs}>1.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt1_ThetaHE", "Kinematic correction vs cos #theta^{*}_{HE}, 0.01.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt2_ThetaHE", "Kinematic correction vs cos #theta^{*}_{HE}, 0.81.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt3_ThetaHE", "Kinematic correction vs cos #theta^{*}_{HE}, 1.41.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt4_ThetaHE", "Kinematic correction vs cos #theta^{*}_{HE}, 2.81.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + {"accCorrectionPt5_ThetaHE", "Kinematic correction vs cos #theta^{*}_{HE}, 5.01.0 GeV/c)/(J/#Psi in |y|<0.88)"}, + // tracking corrections with SPD any and TPCncls>90, pt dependence + {"trackingCorrectionSPDany_pt", "Tracking correction vs pt, 0.090"}, + {"trackingCorrectionSPDanyPt1_pt", "Tracking correction vs pt, 0.090"}, + {"trackingCorrectionSPDanyPt2_pt", "Tracking correction vs pt, 0.890"}, + {"trackingCorrectionSPDanyPt3_pt", "Tracking correction vs pt, 1.490"}, + {"trackingCorrectionSPDanyPt4_pt", "Tracking correction vs pt, 2.890"}, + {"trackingCorrectionSPDanyPt5_pt", "Tracking correction vs pt, 5.090"}, + // tracking corrections with SPD any and TPCncls>90, y dependence + {"trackingCorrectionSPDany_y", "Tracking correction vs y, 0.090"}, + {"trackingCorrectionSPDanyPt1_y", "Tracking correction vs y, 0.090"}, + {"trackingCorrectionSPDanyPt2_y", "Tracking correction vs y, 0.890"}, + {"trackingCorrectionSPDanyPt3_y", "Tracking correction vs y, 1.490"}, + {"trackingCorrectionSPDanyPt4_y", "Tracking correction vs y, 2.890"}, + {"trackingCorrectionSPDanyPt5_y", "Tracking correction vs y, 5.090"}, + // tracking corrections with SPD any and TPCncls>90, pt-y dependence + {"trackingCorrectionSPDany_pty", "Tracking correction vs pt-y, 0.090"}, + {"trackingCorrectionSPDanyPt1_pty", "Tracking correction vs pt-y, 0.090"}, + {"trackingCorrectionSPDanyPt2_pty", "Tracking correction vs pt-y, 0.890"}, + {"trackingCorrectionSPDanyPt3_pty", "Tracking correction vs pt-y, 1.490"}, + {"trackingCorrectionSPDanyPt4_pty", "Tracking correction vs pt-y, 2.890"}, + {"trackingCorrectionSPDanyPt5_pty", "Tracking correction vs pt-y, 5.090"}, + // tracking corrections with SPD any and TPCncls>90, cos Theta*CS dependence + {"trackingCorrectionSPDany_ThetaCS", "Tracking correction vs cos #theta^{*}_{CS}, 0.090"}, + {"trackingCorrectionSPDanyPt1_ThetaCS", "Tracking correction vs cos #theta^{*}_{CS}, 0.090"}, + {"trackingCorrectionSPDanyPt2_ThetaCS", "Tracking correction vs cos #theta^{*}_{CS}, 0.890"}, + {"trackingCorrectionSPDanyPt3_ThetaCS", "Tracking correction vs cos #theta^{*}_{CS}, 1.490"}, + {"trackingCorrectionSPDanyPt4_ThetaCS", "Tracking correction vs cos #theta^{*}_{CS}, 2.890"}, + {"trackingCorrectionSPDanyPt5_ThetaCS", "Tracking correction vs cos #theta^{*}_{CS}, 5.090"}, + // tracking corrections with SPD any and TPCncls>90, cos Theta*HE dependence + {"trackingCorrectionSPDany_ThetaHE", "Tracking correction vs cos #theta^{*}_{HE}, 0.090"}, + {"trackingCorrectionSPDanyPt1_ThetaHE", "Tracking correction vs cos #theta^{*}_{HE}, 0.090"}, + {"trackingCorrectionSPDanyPt2_ThetaHE", "Tracking correction vs cos #theta^{*}_{HE}, 0.890"}, + {"trackingCorrectionSPDanyPt3_ThetaHE", "Tracking correction vs cos #theta^{*}_{HE}, 1.490"}, + {"trackingCorrectionSPDanyPt4_ThetaHE", "Tracking correction vs cos #theta^{*}_{HE}, 2.890"}, + {"trackingCorrectionSPDanyPt5_ThetaHE", "Tracking correction vs cos #theta^{*}_{HE}, 5.090"}, + + // tracking corrections with SPD first and TPCncls>90, pt dependence + {"trackingCorrectionSPDfirst_pt", "Tracking correction vs pt, 0.090"}, + {"trackingCorrectionSPDfirstPt1_pt", "Tracking correction vs pt, 0.090"}, + {"trackingCorrectionSPDfirstPt2_pt", "Tracking correction vs pt, 0.890"}, + {"trackingCorrectionSPDfirstPt3_pt", "Tracking correction vs pt, 1.490"}, + {"trackingCorrectionSPDfirstPt4_pt", "Tracking correction vs pt, 2.890"}, + {"trackingCorrectionSPDfirstPt5_pt", "Tracking correction vs pt, 5.090"}, + // tracking corrections with SPD first and TPCncls>90, y dependence + {"trackingCorrectionSPDfirst_y", "Tracking correction vs y, 0.090"}, + {"trackingCorrectionSPDfirstPt1_y", "Tracking correction vs y, 0.090"}, + {"trackingCorrectionSPDfirstPt2_y", "Tracking correction vs y, 0.890"}, + {"trackingCorrectionSPDfirstPt3_y", "Tracking correction vs y, 1.490"}, + {"trackingCorrectionSPDfirstPt4_y", "Tracking correction vs y, 2.890"}, + {"trackingCorrectionSPDfirstPt5_y", "Tracking correction vs y, 5.090"}, + // tracking corrections with SPD first and TPCncls>90, pt-y dependence + {"trackingCorrectionSPDfirst_pty", "Tracking correction vs pt-y, 0.090"}, + {"trackingCorrectionSPDfirstPt1_pty", "Tracking correction vs pt-y, 0.090"}, + {"trackingCorrectionSPDfirstPt2_pty", "Tracking correction vs pt-y, 0.890"}, + {"trackingCorrectionSPDfirstPt3_pty", "Tracking correction vs pt-y, 1.490"}, + {"trackingCorrectionSPDfirstPt4_pty", "Tracking correction vs pt-y, 2.890"}, + {"trackingCorrectionSPDfirstPt5_pty", "Tracking correction vs pt-y, 5.090"}, + // tracking corrections with SPD first and TPCncls>90, cos Theta*CS dependence + {"trackingCorrectionSPDfirst_ThetaCS", "Tracking correction vs cos #theta^{*}_{CS}, 0.090"}, + {"trackingCorrectionSPDfirstPt1_ThetaCS", "Tracking correction vs cos #theta^{*}_{CS}, 0.090"}, + {"trackingCorrectionSPDfirstPt2_ThetaCS", "Tracking correction vs cos #theta^{*}_{CS}, 0.890"}, + {"trackingCorrectionSPDfirstPt3_ThetaCS", "Tracking correction vs cos #theta^{*}_{CS}, 1.490"}, + {"trackingCorrectionSPDfirstPt4_ThetaCS", "Tracking correction vs cos #theta^{*}_{CS}, 2.890"}, + {"trackingCorrectionSPDfirstPt5_ThetaCS", "Tracking correction vs cos #theta^{*}_{CS}, 5.090"}, + // tracking corrections with SPD first and TPCncls>90, cos Theta*HE dependence + {"trackingCorrectionSPDfirst_ThetaHE", "Tracking correction vs cos #theta^{*}_{HE}, 0.090"}, + {"trackingCorrectionSPDfirstPt1_ThetaHE", "Tracking correction vs cos #theta^{*}_{HE}, 0.090"}, + {"trackingCorrectionSPDfirstPt2_ThetaHE", "Tracking correction vs cos #theta^{*}_{HE}, 0.890"}, + {"trackingCorrectionSPDfirstPt3_ThetaHE", "Tracking correction vs cos #theta^{*}_{HE}, 1.490"}, + {"trackingCorrectionSPDfirstPt4_ThetaHE", "Tracking correction vs cos #theta^{*}_{HE}, 2.890"}, + {"trackingCorrectionSPDfirstPt5_ThetaHE", "Tracking correction vs cos #theta^{*}_{HE}, 5.090"}, + // PID corrections with SPD any and TPCncls>90, pt dependence + {"pidCorrectionSPDany_pt", "PID correction vs pt, 0.090"}, + {"pidCorrectionSPDanyPt1_pt", "PID correction vs pt, 0.090"}, + {"pidCorrectionSPDanyPt2_pt", "PID correction vs pt, 0.890"}, + {"pidCorrectionSPDanyPt3_pt", "PID correction vs pt, 1.490"}, + {"pidCorrectionSPDanyPt4_pt", "PID correction vs pt, 2.890"}, + {"pidCorrectionSPDanyPt5_pt", "PID correction vs pt, 5.090"}, + // PID corrections with SPD any and TPCncls>90, y dependence + {"pidCorrectionSPDany_y", "PID correction vs y, 0.090"}, + {"pidCorrectionSPDanyPt1_y", "PID correction vs y, 0.090"}, + {"pidCorrectionSPDanyPt2_y", "PID correction vs y, 0.890"}, + {"pidCorrectionSPDanyPt3_y", "PID correction vs y, 1.490"}, + {"pidCorrectionSPDanyPt4_y", "PID correction vs y, 2.890"}, + {"pidCorrectionSPDanyPt5_y", "PID correction vs y, 5.090"}, + // PID corrections with SPD any and TPCncls>90, pt-y dependence + {"pidCorrectionSPDany_pty", "PID correction vs pt-y, 0.090"}, + {"pidCorrectionSPDanyPt1_pty", "PID correction vs pt-y, 0.090"}, + {"pidCorrectionSPDanyPt2_pty", "PID correction vs pt-y, 0.890"}, + {"pidCorrectionSPDanyPt3_pty", "PID correction vs pt-y, 1.490"}, + {"pidCorrectionSPDanyPt4_pty", "PID correction vs pt-y, 2.890"}, + {"pidCorrectionSPDanyPt5_pty", "PID correction vs pt-y, 5.090"}, + // PID corrections with SPD any and TPCncls>90, cos Theta*CS dependence + {"pidCorrectionSPDany_ThetaCS", "PID correction vs cos #theta^{*}_{CS}, 0.090"}, + {"pidCorrectionSPDanyPt1_ThetaCS", "PID correction vs cos #theta^{*}_{CS}, 0.090"}, + {"pidCorrectionSPDanyPt2_ThetaCS", "PID correction vs cos #theta^{*}_{CS}, 0.890"}, + {"pidCorrectionSPDanyPt3_ThetaCS", "PID correction vs cos #theta^{*}_{CS}, 1.490"}, + {"pidCorrectionSPDanyPt4_ThetaCS", "PID correction vs cos #theta^{*}_{CS}, 2.890"}, + {"pidCorrectionSPDanyPt5_ThetaCS", "PID correction vs cos #theta^{*}_{CS}, 5.090"}, + // PID corrections with SPD any and TPCncls>90, cos Theta*HE dependence + {"pidCorrectionSPDany_ThetaHE", "PID correction vs cos #theta^{*}_{HE}, 0.090"}, + {"pidCorrectionSPDanyPt1_ThetaHE", "PID correction vs cos #theta^{*}_{HE}, 0.090"}, + {"pidCorrectionSPDanyPt2_ThetaHE", "PID correction vs cos #theta^{*}_{HE}, 0.890"}, + {"pidCorrectionSPDanyPt3_ThetaHE", "PID correction vs cos #theta^{*}_{HE}, 1.490"}, + {"pidCorrectionSPDanyPt4_ThetaHE", "PID correction vs cos #theta^{*}_{HE}, 2.890"}, + {"pidCorrectionSPDanyPt5_ThetaHE", "PID correction vs cos #theta^{*}_{HE}, 5.090"}, + // PID corrections with SPD first and TPCncls>90, pt dependence + {"pidCorrectionSPDfirst_pt", "PID correction vs pt, 0.090"}, + {"pidCorrectionSPDfirstPt1_pt", "PID correction vs pt, 0.090"}, + {"pidCorrectionSPDfirstPt2_pt", "PID correction vs pt, 0.890"}, + {"pidCorrectionSPDfirstPt3_pt", "PID correction vs pt, 1.490"}, + {"pidCorrectionSPDfirstPt4_pt", "PID correction vs pt, 2.890"}, + {"pidCorrectionSPDfirstPt5_pt", "PID correction vs pt, 5.090"}, + // PID corrections with SPD first and TPCncls>90, y dependence + {"pidCorrectionSPDfirst_y", "PID correction vs y, 0.090"}, + {"pidCorrectionSPDfirstPt1_y", "PID correction vs y, 0.090"}, + {"pidCorrectionSPDfirstPt2_y", "PID correction vs y, 0.890"}, + {"pidCorrectionSPDfirstPt3_y", "PID correction vs y, 1.490"}, + {"pidCorrectionSPDfirstPt4_y", "PID correction vs y, 2.890"}, + {"pidCorrectionSPDfirstPt5_y", "PID correction vs y, 5.090"}, + // PID corrections with SPD first and TPCncls>90, pt-y dependence + {"pidCorrectionSPDfirst_pty", "PID correction vs pt-y, 0.090"}, + {"pidCorrectionSPDfirstPt1_pty", "PID correction vs pt-y, 0.090"}, + {"pidCorrectionSPDfirstPt2_pty", "PID correction vs pt-y, 0.890"}, + {"pidCorrectionSPDfirstPt3_pty", "PID correction vs pt-y, 1.490"}, + {"pidCorrectionSPDfirstPt4_pty", "PID correction vs pt-y, 2.890"}, + {"pidCorrectionSPDfirstPt5_pty", "PID correction vs pt-y, 5.090"}, + // PID corrections with SPD first and TPCncls>90, cos Theta*CS dependence + {"pidCorrectionSPDfirst_ThetaCS", "PID correction vs cos #theta^{*}_{CS}, 0.090"}, + {"pidCorrectionSPDfirstPt1_ThetaCS", "PID correction vs cos #theta^{*}_{CS}, 0.090"}, + {"pidCorrectionSPDfirstPt2_ThetaCS", "PID correction vs cos #theta^{*}_{CS}, 0.890"}, + {"pidCorrectionSPDfirstPt3_ThetaCS", "PID correction vs cos #theta^{*}_{CS}, 1.490"}, + {"pidCorrectionSPDfirstPt4_ThetaCS", "PID correction vs cos #theta^{*}_{CS}, 2.890"}, + {"pidCorrectionSPDfirstPt5_ThetaCS", "PID correction vs cos #theta^{*}_{CS}, 5.090"}, + // PID corrections with SPD first and TPCncls>90, cos Theta*HE dependence + {"pidCorrectionSPDfirst_ThetaHE", "PID correction vs cos #theta^{*}_{HE}, 0.090"}, + {"pidCorrectionSPDfirstPt1_ThetaHE", "PID correction vs cos #theta^{*}_{HE}, 0.090"}, + {"pidCorrectionSPDfirstPt2_ThetaHE", "PID correction vs cos #theta^{*}_{HE}, 0.890"}, + {"pidCorrectionSPDfirstPt3_ThetaHE", "PID correction vs cos #theta^{*}_{HE}, 1.490"}, + {"pidCorrectionSPDfirstPt4_ThetaHE", "PID correction vs cos #theta^{*}_{HE}, 2.890"}, + {"pidCorrectionSPDfirstPt5_ThetaHE", "PID correction vs cos #theta^{*}_{HE}, 5.090"} +}; + + +// Function prototypes --------------------------------------------------- +// The user must modify the DefineHistograms() and FillHistograms() functions +// according to need +Double_t* GetBinning(AliCFContainer* cont, Int_t variable, Int_t& nBins); +void FillHistograms(TObjArray* histosArray, AliCFContainer* cont, Int_t currentRangeStep, Bool_t firstTime); +void GetBinLimits(AliCFContainer* cont); +void DefineHistograms(TObjArray* objArray, Int_t iCutSet); +void AddHistogram(TObjArray* objArray, Int_t ndim, + const Char_t* name, const Char_t* title, + Int_t nbinsx, Double_t* binsx, const Char_t* xLabel = "", + Int_t nbinsy=0, Double_t* binsy=0, const Char_t* yLabel = "", + Int_t nbinsz=0, Double_t* binsz=0, const Char_t* zLabel = ""); +void ProjectManyRuns(const Char_t* runList, Int_t howMany=1, Int_t offset = 0); +void ProjectAll(const Char_t* inputList, const Char_t* outfilename="HistosFromCFs.root", + Int_t howMany=1, Int_t offset=0); +void ExtractEfficienciesMany(const Char_t* runList, Int_t howMany=1, Int_t offset=0); +void ExtractEfficiencies(const Char_t* inputFile, const Char_t* outfilename="Efficiencies.root", const Char_t* numbersFile=""); +TH1* DivideHists(TH1* nominator, TH1* denominator); +//------------------------------------------------------------------------- + + +//_______________________________________________________________________________________ +void ProjectManyRuns(const Char_t* runList, Int_t howMany, Int_t offset) { + // + // + // + + // loop over all runs ----------------------- + ifstream input; input.open(runList); + Int_t runCounter = 0; + while(input.good()) { + Char_t readString[256]; + input.getline(readString, 256, '\n'); // get a chunk + TString runStr = readString; + Int_t run = runStr.Atoi(); + if(run<=0) continue; + + if(runCounter=offset+howMany) + break; + + cout << "=================== run " << run << " ============================" << endl; + + ProjectAll(Form("LHC10f7a/invMass_BB1/%s/listCF.txt",readString), + Form("LHC10f7a/invMass_BB1/%s/Projections.root",readString), + 100, 0); + runCounter++; + } +} + + +//_______________________________________________________________________________________ +void ProjectAll(const Char_t* inputList, + const Char_t* outfilename, + Int_t howMany, Int_t offset) { + // + // Main function for making projections from a list of CF containers (inputList). + // The resulting histograms are placed in the ROOT file specified by outfilename + // + // Modify the global variables above to match your requirements + // + + // open the output file + TFile *outFile = new TFile(outfilename,"RECREATE"); + // ----------------------------------------------------------------------------- + + // copy the current ExtractEfficiency macro in the same dir as the output file + TString outStr = ""; + outStr += outfilename; + outStr.ReplaceAll(".root", "_ExtractEfficienciesMacro.C"); + gSystem->Exec(Form("cp ExtractEfficiencies.C %s", outStr.Data())); + // --------------------------------------------------------------------------- + + // create the container for all the histograms --------------------- + TObjArray *histoArray=new TObjArray(); + histoArray->SetOwner(); + //------------------------------------------------------------------ + + + // loop over all CF files, project and merge ----------------------- + ifstream input; input.open(inputList); + Int_t currentFile=0; + Bool_t firstTime = kTRUE; + while(input.good()) { + Char_t readString[256]; + input.getline(readString, 256, '\n'); // get a chunk + TString readStringString = readString; + if(readStringString[0]!='/') continue; + if(!readStringString.Contains(".root")) continue; + + if(currentFile=offset+howMany) + break; + + cout << "file: " << readString << endl; + + AliDielectronCFdraw *cf=new AliDielectronCFdraw(readString); + AliCFContainer* cont=cf->GetCFContainer(); + + // **************************************************************************** + // Below apply all your cut sets then call the FillHistograms() function + // Don't forget to increment the "currentCutSet" variable after every cut set + // **************************************************************************** + + // pair type (0 ++, 1 +-, 2 --) ---------------------------------- + cf->SetRangeUser("PairType", 1, 1); + // Pair rapidity cut + cf->SetRangeUser("Y", -0.899, 0.899); + Int_t currentCutSet = 0; + FillHistograms(histoArray, cont, currentCutSet, firstTime); + + // j/psi 0SetRangeUser("Pt", 0.001, 0.799); + currentCutSet++; + FillHistograms(histoArray,cont,currentCutSet, firstTime); + + // j/psi 0.8SetRangeUser("Pt", 0.801, 1.399); + currentCutSet++; + FillHistograms(histoArray,cont,currentCutSet, firstTime); + + // j/psi 1.4SetRangeUser("Pt", 1.401, 2.799); + currentCutSet++; + FillHistograms(histoArray,cont,currentCutSet, firstTime); + + // j/psi 2.8SetRangeUser("Pt", 2.801, 4.999); + currentCutSet++; + FillHistograms(histoArray,cont,currentCutSet, firstTime); + + // j/psi 5.0SetRangeUser("Pt", 5.001, 9.999); + currentCutSet++; + FillHistograms(histoArray,cont,currentCutSet, firstTime); + + // Leg pseudo-rapidity cut --------------------------------------- + cf->SetRangeUser("Pt", 0.001, 9.999); + cf->SetRangeUser("Leg1_Eta", -0.899, 0.899); + cf->SetRangeUser("Leg2_Eta", -0.899, 0.899); + cf->SetRangeUser("Leg1_Pt", 0.801, 10.0); + cf->SetRangeUser("Leg2_Pt", 0.801, 10.0); + currentCutSet++; + FillHistograms(histoArray,cont,currentCutSet, firstTime); + + // j/psi 0SetRangeUser("Pt", 0.001, 0.799); + currentCutSet++; + FillHistograms(histoArray,cont,currentCutSet, firstTime); + + // j/psi 0.8SetRangeUser("Pt", 0.801, 1.399); + currentCutSet++; + FillHistograms(histoArray,cont,currentCutSet, firstTime); + + // j/psi 1.4SetRangeUser("Pt", 1.401, 2.799); + currentCutSet++; + FillHistograms(histoArray,cont,currentCutSet, firstTime); + + // j/psi 2.8SetRangeUser("Pt", 2.801, 4.999); + currentCutSet++; + FillHistograms(histoArray,cont,currentCutSet, firstTime); + + // j/psi 5.0SetRangeUser("Pt", 5.001, 9.999); + currentCutSet++; + FillHistograms(histoArray,cont,currentCutSet, firstTime); + + // |LegEta|<0.88 & Leg_Pt>1.0 & TPCncls>90 ------------------------------------------ + cf->SetRangeUser("Pt", 0.001, 9.999); + cf->SetRangeUser("Leg1_NclsTPC", 90.1, 160.0); + cf->SetRangeUser("Leg2_NclsTPC", 90.1, 160.0); + currentCutSet++; + FillHistograms(histoArray,cont,currentCutSet, firstTime); + + // j/psi 0SetRangeUser("Pt", 0.001, 0.799); + currentCutSet++; + FillHistograms(histoArray,cont,currentCutSet, firstTime); + + // j/psi 0.8SetRangeUser("Pt", 0.801, 1.399); + currentCutSet++; + FillHistograms(histoArray,cont,currentCutSet, firstTime); + + // j/psi 1.4SetRangeUser("Pt", 1.401, 2.799); + currentCutSet++; + FillHistograms(histoArray,cont,currentCutSet, firstTime); + + // j/psi 2.8SetRangeUser("Pt", 2.801, 4.999); + currentCutSet++; + FillHistograms(histoArray,cont,currentCutSet, firstTime); + + // j/psi 5.0SetRangeUser("Pt", 5.001, 9.999); + currentCutSet++; + FillHistograms(histoArray,cont,currentCutSet, firstTime); + + + currentFile++; + firstTime = kFALSE; + delete cont; + delete cf; + } // end loop over CF files + + outFile->cd(); + histoArray->Write(); + outFile->Close(); + delete histoArray; + return; +} + +//_______________________________________________________________________________________ +void ExtractEfficienciesMany(const Char_t* runList, Int_t howMany, Int_t offset) { + // + // + // + + // loop over all runs ----------------------- + ifstream input; input.open(runList); + Int_t runCounter = 0; + + TGraphErrors* trends[gkNeffs]; + Double_t weightedEffs[gkNeffs]; + Double_t weightedErrs[gkNeffs]; + Double_t nTotalEvents = 0; + for(Int_t iTrend=0; iTrendSetName(gkEffNames[iTrend][0]); + trends[iTrend]->SetTitle(gkEffNames[iTrend][1]); + weightedEffs[iTrend] = 0.0; weightedErrs[iTrend] = 0.0; + } + TFile* file=0x0; + TFile* normalizationFile=0x0; + TNamed* object; + while(input.good()) { + Char_t readString[256]; + input.getline(readString, 256, '\n'); // get a chunk + TString runStr = readString; + Int_t run = runStr.Atoi(); + if(run<=0) continue; + + if(runCounter=offset+howMany) + break; + + cout << "=================== run " << run << " ============================" << endl; + TString periodStr; + if(run<=117222) periodStr = "LHC10b.pass2"; + if(run>117222 && run<=120829) periodStr = "LHC10c.pass2"; + if(run>=122374 && run<=126437) periodStr = "LHC10d.pass1"; + Double_t nPhysicsEvents = 0; + normalizationFile = TFile::Open(Form("/u/iarsene/work/ALICE/normalization/2010-10-03_0628.3293/%s/%d/iarsene_normalization.root", periodStr.Data(), run)); + if(normalizationFile) { + TObjArray *histos=(TObjArray*)normalizationFile->Get("iarsene_normalization"); + TH1I* triggers=(TH1I*)histos->FindObject("TriggersHistogram"); + nPhysicsEvents = triggers->GetBinContent(2); // PHYSICS events + normalizationFile->Close(); + nTotalEvents += nPhysicsEvents; + } + + ExtractEfficiencies(Form("LHC10f7a/iter10_BB2/%s/Projections_iter10.root",readString), + Form("LHC10f7a/iter10_BB2/%s/Efficiencies_iter10.root",readString)); + file = TFile::Open(Form("LHC10f7a/iter10_BB2/%s/Efficiencies_iter10.root",readString)); + if(file && !file->IsZombie()) { + for(Int_t iTrend=0; iTrendGet(Form("%s_value",gkEffNames[iTrend][0])); + if(!object) continue; + Float_t eff = (TString(object->GetTitle())).Atof(); + trends[iTrend]->SetPoint(runCounter, run, eff); + object = (TNamed*)file->Get(Form("%s_error",gkEffNames[iTrend][0])); + if(!object) continue; + Float_t err = (TString(object->GetTitle())).Atof(); + trends[iTrend]->SetPointError(runCounter, 0.0, err); + weightedEffs[iTrend] += nPhysicsEvents*eff; + weightedErrs[iTrend] += nPhysicsEvents*nPhysicsEvents*err*err; + } + file->Close(); + } + if(normalizationFile) + normalizationFile->Close(); + + runCounter++; + } + + TFile *saveTrend = new TFile(Form("%s.trend_iter10.root", runList), "RECREATE"); + TNamed *weightedFactors; + TNamed *weightedErrors; + for(Int_t iTrend=0; iTrendWrite(); + weightedEffs[iTrend] /= nTotalEvents; + weightedErrs[iTrend] = TMath::Sqrt(weightedErrs[gkNeffs]/nTotalEvents/nTotalEvents); + weightedFactors = new TNamed(Form("%s_weighted", gkEffNames[iTrend][0]), + Form("%f", weightedEffs[iTrend])); + weightedErrors = new TNamed(Form("%s_weightedErr", gkEffNames[iTrend][0]), + Form("%f", weightedErrs[iTrend])); + weightedFactors->Write(); + weightedErrors->Write(); + } + weightedFactors = new TNamed("TotalEvents", Form("%f",nTotalEvents)); + weightedFactors->Write(); + + saveTrend->Close(); +} + +//_______________________________________________________________________________________ +void ExtractEfficiencies(const Char_t* inputFilename, + const Char_t* outfilename, + const Char_t* numbersFile) { + // + // Main function to extract efficiencies + // + + // Open the output file + TFile *output = new TFile(outfilename, "RECREATE"); + // --------------------------------------------------------------------------- + + // copy the current ExtractEfficiency macro in the same dir as the output file + TString outStr = ""; + outStr += outfilename; + outStr.ReplaceAll(".root", "_ExtractEfficienciesMacro.C"); + gSystem->Exec(Form("cp ExtractEfficiencies.C %s", outStr.Data())); + // --------------------------------------------------------------------------- + + // open the input file and read out all the histograms + TFile *file = TFile::Open(inputFilename); + if(!file || file->IsZombie()) return; + + TObjArray *effArray = new TObjArray(); + effArray->SetOwner(); + TH1* nominator; + TH1* denominator; + TNamed* effValue; + TNamed* effError; + ofstream asciiOut; + if(numbersFile[0]!='\0') { + asciiOut.open(numbersFile); + asciiOut << "#Format: Name | Value | Abs. Error" << endl; + } + for(Int_t iEff=0; iEffGet(Form("%s_%s_%s",gkStepNames[gkEffs[iEff][0]][0], + gkCutSetNames[gkEffs[iEff][1]][0], + gkHistoNames[gkEffs[iEff][2]][0]))); + denominator = (TH1D*)(file->Get(Form("%s_%s_%s",gkStepNames[gkEffs[iEff][3]][0], + gkCutSetNames[gkEffs[iEff][4]][0], + gkHistoNames[gkEffs[iEff][5]][0]))); + if(!nominator) continue; + if(!denominator) continue; + nominator->GetYaxis()->SetTitle("efficiency"); + } + if(gkDims[gkEffs[iEff][2]][0]==2) { // 2-dim histos + nominator = (TH2D*)(file->Get(Form("%s_%s_%s",gkStepNames[gkEffs[iEff][0]][0], + gkCutSetNames[gkEffs[iEff][1]][0], + gkHistoNames[gkEffs[iEff][2]][0]))); + denominator = (TH2D*)(file->Get(Form("%s_%s_%s",gkStepNames[gkEffs[iEff][3]][0], + gkCutSetNames[gkEffs[iEff][4]][0], + gkHistoNames[gkEffs[iEff][5]][0]))); + if(!nominator) continue; + if(!denominator) continue; + nominator->GetZaxis()->SetTitle("efficiency"); + } + if(gkDims[gkEffs[iEff][2]][0]==3) { // 3-dim histos + nominator = (TH3D*)(file->Get(Form("%s_%s_%s",gkStepNames[gkEffs[iEff][0]][0], + gkCutSetNames[gkEffs[iEff][1]][0], + gkHistoNames[gkEffs[iEff][2]][0]))); + denominator = (TH3D*)(file->Get(Form("%s_%s_%s",gkStepNames[gkEffs[iEff][3]][0], + gkCutSetNames[gkEffs[iEff][4]][0], + gkHistoNames[gkEffs[iEff][5]][0]))); + if(!nominator) continue; + if(!denominator) continue; + } + Double_t nomIntegral = nominator->Integral(); + Double_t denomIntegral = denominator->Integral(); + Double_t eff = (denomIntegral>0 ? nomIntegral/denomIntegral : 0); + Double_t error = (nomIntegral>0 && denomIntegral>0 ? eff*TMath::Sqrt(1.0/nomIntegral + 1.0/denomIntegral) : 0); + //nominator->Divide(denominator); + TH1* ratio = DivideHists(nominator, denominator); + cout << "efficiency = " << nomIntegral << " / " << denomIntegral << " = " + << eff << " +/- " << error << endl; + TString title = gkEffNames[iEff][1]; + title += Form(", integrated eff. = %f #pm %f", eff, error); + ratio->SetTitle(title.Data()); + ratio->SetName(gkEffNames[iEff][0]); + effArray->Add(ratio); + + if(numbersFile[0]!='\0') { + asciiOut << gkEffNames[iEff][0] << "\t" << eff << "\t" << error << endl; + } + effValue = new TNamed(Form("%s_value", gkEffNames[iEff][0]), + Form("%f", eff)); + effArray->Add(effValue); + effError = new TNamed(Form("%s_error", gkEffNames[iEff][0]), + Form("%f", error)); + effArray->Add(effError); + } + + output->cd(); + effArray->Write(); + output->Close(); + file->Close(); + asciiOut.close(); +} + + +//________________________________________________________________________________________ +void DefineHistograms(TObjArray* histoArray, Int_t iCutSet) { + // + // Define the histograms to be filled for every step and a given cut set + // This function is called by the FillHistograms() if the firstTime flag is set + + for(Int_t iStep=0; iStep3) return; + TH1* histo; + if(ndim==1) { + histo = new TH1D(name, title, nbinsx, binsx); + histo->Sumw2(); + histo->GetXaxis()->SetTitle(xLabel); + } + if(ndim==2) { + histo = new TH2D(name, title, nbinsx, binsx, nbinsy, binsy); + histo->Sumw2(); + histo->GetXaxis()->SetTitle(xLabel); + histo->GetYaxis()->SetTitle(yLabel); + } + if(ndim==3) { + histo = new TH3D(name, title, nbinsx, binsx, nbinsy, binsy, nbinsz, binsz); + histo->Sumw2(); + histo->GetXaxis()->SetTitle(xLabel); + histo->GetYaxis()->SetTitle(yLabel); + histo->GetZaxis()->SetTitle(zLabel); + } + histoArray->Add(histo); +} + +//__________________________________________________________________________________________ +void FillHistograms(TObjArray* histosArray, AliCFContainer* cont, Int_t currentCutSet, Bool_t firstTime) { + // + // Fill the user defined histograms for a given cut set + // + // If the firstTime flag is on then update the bin limits and call DefineHistograms() + if(firstTime) { + GetBinLimits(cont); + DefineHistograms(histosArray, currentCutSet); + } + + TH1* histo; + for(Int_t iStep=0; iStepFindObject(Form("%s_%s_%s",gkStepNames[iStep][0], + gkCutSetNames[currentCutSet][0], + gkHistoNames[iHisto][0])); + histo->Add(cont->Project(gkDims[iHisto][1],gkStepNumbers[iStep])); + } + // fill 2-dim histos + if(gkDims[iHisto][0]==2) { + histo = (TH2D*)histosArray->FindObject(Form("%s_%s_%s",gkStepNames[iStep][0], + gkCutSetNames[currentCutSet][0], + gkHistoNames[iHisto][0])); + histo->Add(cont->Project(gkDims[iHisto][1], gkDims[iHisto][2], gkStepNumbers[iStep])); + } + // fill 3-dim histos + if(gkDims[iHisto][0]==3) { + histo = (TH3D*)histosArray->FindObject(Form("%s_%s_%s",gkStepNames[iStep][0], + gkCutSetNames[currentCutSet][0], + gkHistoNames[iHisto][0])); + histo->Add(cont->Project(gkDims[iHisto][1], gkDims[iHisto][2], gkDims[iHisto][3], gkStepNumbers[iStep])); + } + } // end loop over histos + } // end loop over steps +} + +//____________________________________________________________________________________________ +void GetBinLimits(AliCFContainer* cont) { + // + // Extract the bin limits from the CF container + // + cout << "********* New cut set ****************" << endl; + for(Int_t iVar=0; iVarGetVarTitle(iVar) << " : " << gNbins[iVar]; + cout << "; range = " << gBinLimits[iVar][0] << " --> " << gBinLimits[iVar][gNbins[iVar]] << endl; + } +} + +//________________________________________________________________________________________ +Double_t* GetBinning(AliCFContainer* cont, Int_t variable, + Int_t& nBins) { + // + // Get the number of bins and the bin limits for the projection of a given variable + // + TH1D* tempHist = cont->Project(variable, kPureMC); + nBins = tempHist->GetXaxis()->GetNbins(); + Double_t* binLimits = new Double_t[nBins+1]; + for(Int_t i=1; i<=nBins; i++) + binLimits[i-1]=tempHist->GetXaxis()->GetBinLowEdge(i); + binLimits[nBins] = tempHist->GetXaxis()->GetBinLowEdge(nBins) + + tempHist->GetXaxis()->GetBinWidth(nBins); + return binLimits; +} + +//________________________________________________________________________________________ +TH1* DivideHists(TH1* nominator, TH1* denominator) { + // + // divide 2 histograms with error propagation + // + TH1* ratio; + if(nominator->InheritsFrom("TH3")) { + Int_t nBinsXNom = nominator->GetXaxis()->GetNbins(); + Int_t nBinsXDenom = denominator->GetXaxis()->GetNbins(); + Int_t nBinsYNom = nominator->GetYaxis()->GetNbins(); + Int_t nBinsYDenom = denominator->GetYaxis()->GetNbins(); + Int_t nBinsZNom = nominator->GetZaxis()->GetNbins(); + Int_t nBinsZDenom = denominator->GetZaxis()->GetNbins(); + if(nBinsXNom!=nBinsXDenom || nBinsYNom!=nBinsYDenom || nBinsZNom!=nBinsZDenom) { + cout << "Trying to divide histograms with different number of bins" << endl; + return 0x0; + } + ratio = (TH3D*)nominator->Clone("ratio"); + ratio->Reset(); + for(Int_t iXbin=1; iXbin<=nBinsXNom; ++iXbin) { + for(Int_t iYbin=1; iYbin<=nBinsYNom; ++iYbin) { + for(Int_t iZbin=1; iZbin<=nBinsZNom; ++iZbin) { + Double_t countsN = nominator->GetBinContent(iXbin, iYbin, iZbin); + Double_t countsD = denominator->GetBinContent(iXbin, iYbin, iZbin); + if(countsN<1 || countsD<1) continue; // zero entry bins + ratio->SetBinContent(iXbin, iYbin, iZbin, countsN/countsD); + ratio->SetBinError(iXbin, iYbin, iZbin, (countsN/countsD)*TMath::Sqrt(1.0/countsN)+(1.0/countsD)); + } + } + } + return ratio; + } + + if(nominator->InheritsFrom("TH2")) { + Int_t nBinsXNom = nominator->GetXaxis()->GetNbins(); + Int_t nBinsXDenom = denominator->GetXaxis()->GetNbins(); + Int_t nBinsYNom = nominator->GetYaxis()->GetNbins(); + Int_t nBinsYDenom = denominator->GetYaxis()->GetNbins(); + if(nBinsXNom!=nBinsXDenom || nBinsYNom!=nBinsYDenom) { + cout << "Trying to divide histograms with different number of bins" << endl; + return 0x0; + } + ratio = (TH2D*)nominator->Clone("ratio"); + ratio->Reset(); + for(Int_t iXbin=1; iXbin<=nBinsXNom; ++iXbin) { + for(Int_t iYbin=1; iYbin<=nBinsYNom; ++iYbin) { + Double_t countsN = nominator->GetBinContent(iXbin, iYbin); + Double_t countsD = denominator->GetBinContent(iXbin, iYbin); + if(countsN<1 || countsD<1) continue; // zero entry bins + ratio->SetBinContent(iXbin, iYbin, countsN/countsD); + ratio->SetBinError(iXbin, iYbin, (countsN/countsD)*TMath::Sqrt(1.0/countsN)+(1.0/countsD)); + } + } + return ratio; + } + + if(nominator->InheritsFrom("TH1")) { + Int_t nBinsXNom = nominator->GetXaxis()->GetNbins(); + Int_t nBinsXDenom = denominator->GetXaxis()->GetNbins(); + if(nBinsXNom!=nBinsXDenom) { + cout << "Trying to divide histograms with different number of bins" << endl; + return 0x0; + } + ratio = (TH1D*)nominator->Clone("ratio"); + ratio->Reset(); + for(Int_t iXbin=1; iXbin<=nBinsXNom; ++iXbin) { + Double_t countsN = nominator->GetBinContent(iXbin); + Double_t countsD = denominator->GetBinContent(iXbin); + if(countsN<1 || countsD<1) continue; // zero entry bins + ratio->SetBinContent(iXbin, countsN/countsD); + ratio->SetBinError(iXbin, (countsN/countsD)*TMath::Sqrt(1.0/countsN)+(1.0/countsD)); + } + return ratio; + } + + return 0x0; +} diff --git a/PWG3/dielectron/macros/MakeDataReport.C b/PWG3/dielectron/macros/MakeDataReport.C index 0319e4752d1..16778f148f3 100644 --- a/PWG3/dielectron/macros/MakeDataReport.C +++ b/PWG3/dielectron/macros/MakeDataReport.C @@ -343,6 +343,12 @@ Int_t baseColors[5]={kRed, kRed, kRed, kRed, kRed}; alephParameters[2] = 5.04114e-11; alephParameters[3] = 2.12543e+00; alephParameters[4] = 4.88663e+00; + alephParameters[0] = 1.25202/50.; //was 1.79571/55.; + alephParameters[1] = 2.74992e+01; //was 22.0028; + alephParameters[2] = TMath::Exp(-3.31517e+01); //was1.55354e-11; + alephParameters[3] = 2.46246; //was 2.39804; + alephParameters[4] = 6.78938; //was 5.1209; + Double_t mip=50; Color_t color=kRed; @@ -414,9 +420,9 @@ c->SetAlias("nCls","Leg1_NclsTPC>90&&Leg2_NclsTPC>90"); //-- nsigma c->SetAlias("cutE","abs(Leg1_TPC_nSigma_Electrons)<3&&abs(Leg2_TPC_nSigma_Electrons)<3"); -c->SetAlias("cutE","Leg1_TPC_nSigma_Electrons>-1 && Leg2_TPC_nSigma_Electrons>-1"); +// c->SetAlias("cutE","Leg1_TPC_nSigma_Electrons>-1 && Leg2_TPC_nSigma_Electrons>-1"); c->SetAlias("cutPi","abs(Leg1_TPC_nSigma_Pions)>3&&abs(Leg2_TPC_nSigma_Pions)>3"); -c->SetAlias("cutP","(Leg1_TPC_nSigma_Protons)>3.&&(Leg2_TPC_nSigma_Protons)>3."); +c->SetAlias("cutP","(Leg1_TPC_nSigma_Protons)>3&&(Leg2_TPC_nSigma_Protons)>3"); c->SetAlias("pidSig","cutE&&cutPi&&cutP"); //-- Pi param @@ -427,36 +433,43 @@ c->SetAlias("pidParam","eleParam&&cutP"); c->SetAlias("LegEta","abs(Leg1_Eta)<0.9&&abs(Leg2_Eta<0.9)"); +c->SetAlias("LegNcl","Leg1_NclsTPC>90&&Leg2_NclsTPC>90"); c->SetAlias("Rap","abs(Y)<0.9"); c->SetAlias("QA","LegNcl&&LegEta&&Rap"); c->SetAlias("spdFirst","(Leg1_ITS_clusterMap&1)==1 && (Leg2_ITS_clusterMap&1)==1"); +c->SetAlias("LegNclDiffIter1","abs(Leg1_NclsTPC-Leg1_NclsTPCiter1)<10&&abs(Leg2_NclsTPC-Leg2_NclsTPCiter1)<10") +c->SetAlias("LegNclPID","(Leg1_NclsTPC-Leg1_TPCsignalN)<20&&(Leg2_NclsTPC-Leg2_TPCsignalN)<20") - -c->SetAlias("cut","PairType==1&&nCls&&pidSig&&LegEta&&Rap") +c->SetAlias("cut","PairType==1&&QA&&pidSig") c->SetMarkerStyle(20); c->SetMarkerSize(.8); c->SetMarkerColor(kBlack); c->SetLineColor(kBlack); -// c->SetAlias("nCls","Leg1_NclsTPC>120&&Leg2_NclsTPC>120"); +// c->SetAlias("nCls","Leg1_NclsTPC>90&&Leg2_NclsTPC>90"); c->Draw("M>>hM(50,2,4)","cut","e"); c->SetMarkerColor(kBlue); c->SetLineColor(kBlue); -c->SetAlias("cut","PairType==1&&nCls&&pidParam&&LegEta&&Rap") +// c->SetAlias("cut","PairType==1&&nCls&&pidParam&&LegEta&&Rap") +c->SetAlias("cut","PairType==1&&QA&&pidSig&&LegNclPID") // c->SetAlias("nCls","Leg1_NclsTPC>140&&Leg2_NclsTPC>140"); c->Draw("M>>hM2(50,2,4)","cut","esame"); + c->SetMarkerColor(kGreen); c->SetLineColor(kGreen); -c->SetAlias("nCls","Leg1_NclsTPC>150&&Leg2_NclsTPC>150"); +c->SetAlias("cut","PairType==1&&QA&&pidSig") +c->SetAlias("LegNcl","Leg1_NclsTPC>120&&Leg2_NclsTPC>120"); +// c->SetAlias("nCls","Leg1_NclsTPC>150&&Leg2_NclsTPC>150"); c->Draw("M>>hM3(50,2,4)","cut","esame"); -c->SetAlias("cut","PairType==1&&nCls&&LegEta&&Rap") +c->SetAlias("LegNclDiffIter1","(Leg1_NclsTPC-Leg1_NclsTPCiter1)>-1&&(Leg2_NclsTPC-Leg2_NclsTPCiter1)>-1") +c->SetAlias("cut","PairType==1&&LegNclDiffIter1") // histos AliDielectronHistos h("h","h"); h.AddClass("TPCsignal"); @@ -474,6 +487,9 @@ h.GetHistogram("TPCsignal","nSigK")->SetDirectory(gDirectory) h.UserHistogram("TPCsignal","nSigP","TPC n #sigma Protons;P [GeV];TPC n #sigma Protons",400,.3,40.,500,-10,10.,0,0,kTRUE,kFALSE) h.GetHistogram("TPCsignal","nSigP")->SetDirectory(gDirectory) +h.UserHistogram("TPCsignal","nSigDiffP","ncls-nclsXX;P [GeV];ncls-nclsXX",400,.3,40.,200,-40,160.,0,0,kTRUE,kFALSE) +h.GetHistogram("TPCsignal","nSigDiffP")->SetDirectory(gDirectory) + c->Draw("Leg1_TPC_signal:Leg1_P_InnerParam>>sigTPC","cut","colz") c->Draw("Leg2_TPC_signal:Leg2_P_InnerParam>>+sigTPC","cut","colz") @@ -513,6 +529,55 @@ AliDielectronSignalFunc sig; sig.SetDefaults(1); +//---------- +c->Draw("Leg1_NclsTPC>>hNcls(160,-0.5,159.5)","cut","goff"); +c->Draw("Leg2_NclsTPC>>+hNcls","cut","goff"); +hNclsPID->SetLineColor(kBlack) + +c->Draw("Leg1_TPCsignalN>>hNclsPID(160,-0.5,159.5)","cut","goff"); +c->Draw("Leg2_TPCsignalN>>+hNclsPID","cut","goff"); +hNclsPID->SetLineColor(kBlue) + +c->Draw("Leg1_NclsTPCiter1>>hNclsIter1(160,-0.5,159.5)","cut","goff"); +c->Draw("Leg2_NclsTPCiter1>>+hNclsIter1","cut","goff"); +hNclsIter1->SetLineColor(kGreen) + +hNcls->Draw(); +hNclsPID->Draw("same"); +hNclsIter1->Draw("same"); +//----------- + + + +c->Draw("Leg1_TPCsignalN:Leg1_NclsTPC>>hNclsPIDNcls(160,-0.5,159.5,160,-0.5,159.5)","cut","colz"); +c->Draw("Leg2_TPCsignalN:Leg2_NclsTPC>>+hNclsPIDNcls","cut","colz"); + + +c->Draw("Leg1_NclsTPC-Leg1_TPCsignalN:Leg1_P_InnerParam>>nSigDiffP","cut","colz"); +c->Draw("Leg2_NclsTPC-Leg2_TPCsignalN:Leg1_P_InnerParam>>+nSigDiffP","cut","colz"); + + +c->Draw("Leg1_NclsTPC-Leg1_NclsTPCiter1:Leg1_P_InnerParam>>nSigDiffP","cut","colz"); +c->Draw("Leg2_NclsTPC-Leg2_NclsTPCiter1:Leg1_P_InnerParam>>+nSigDiffP","cut","colz"); + + + + + + + + + +c->Draw("Leg1_NclsTPC-Leg1_TPCsignalN:Leg1_TPC_nSigma_Electrons:Leg1_P_InnerParam>>hXX(100,0,10,20,-4,4)","cut","profcolz") +c->Draw("Leg1_NclsTPC-Leg1_TPCsignalN:Leg2_TPC_nSigma_Electrons:Leg2_P_InnerParam>>+hXX","cut","profcolz") + + + + + + + + //WooJins cuts: c->SetAlias("cutE","abs(Leg1_TPC_nSigma_Electrons)<3&&abs(Leg2_TPC_nSigma_Electrons)<3"); diff --git a/PWG3/dielectron/macros/PlotDataResults.C b/PWG3/dielectron/macros/PlotDataResults.C index 2d24597ae66..e24b6677884 100644 --- a/PWG3/dielectron/macros/PlotDataResults.C +++ b/PWG3/dielectron/macros/PlotDataResults.C @@ -41,9 +41,18 @@ void PlotDataResults(const char* filenameData, const char* filenameMC="", Bool_t //Set common Ranges // d.SetRangeUser("Leg1_Pt",0.8,1000.); // d.SetRangeUser("Leg2_Pt",0.8,1000.); -// d.SetRangeUser("Leg1_NclsTPC",120.,170.); -// d.SetRangeUser("Leg2_NclsTPC",120.,170.); - d.SetRangeUser("M",1.,5.); + d.SetRangeUser("Leg1_NclsTPC",140.,170.); + d.SetRangeUser("Leg2_NclsTPC",140.,170.); + d.SetRangeUser("Leg1_Pt",1.01,100000); + d.SetRangeUser("Leg2_Pt",1.01,100000); +// d.SetRangeUser("Leg1_TPC_nSigma_Electrons",-3,3); +// d.SetRangeUser("Leg2_TPC_nSigma_Electrons",-3,3); +// d.SetRangeUser("Leg1_TPC_nSigma_Pions",3,20); +// d.SetRangeUser("Leg2_TPC_nSigma_Pions",3,20); +// d.SetRangeUser("Leg1_TPC_nSigma_Protons",3,20); +// d.SetRangeUser("Leg2_TPC_nSigma_Protons",3,20); + + d.SetRangeUser("M",0.5,5.); //============================ //SPD first // @@ -51,7 +60,7 @@ void PlotDataResults(const char* filenameData, const char* filenameMC="", Bool_t //--- Like sign subtraction AliDielectronSignalBase *sigFirst=GetSignalLS(d,stepFirst); SetStyle(sigFirst,"ITS First - Like Sign subtraction"); - DrawSpectra(sigFirst,"cFirst",hStats,save); +// DrawSpectra(sigFirst,"cFirst",hStats,save); //--- Rotation subtraction AliDielectronSignalBase *sigFirstRot=GetSignalRot(d,stepFirst); SetStyle(sigFirstRot,"ITS First - Track rotation subtraction"); @@ -66,7 +75,7 @@ void PlotDataResults(const char* filenameData, const char* filenameMC="", Bool_t //--- Rotation subtraction AliDielectronSignalBase *sigAnyRot=GetSignalRot(d,stepAny); SetStyle(sigAnyRot,"ITS First - Track rotation subtraction"); -// DrawSpectra(sigAnyRot,"cAnyRot",save); +// DrawSpectra(sigAnyRot,"cAnyRot",hStats,save); //============================= //TOF up to 1.2, parametrisation in TPC ele @@ -77,7 +86,7 @@ void PlotDataResults(const char* filenameData, const char* filenameMC="", Bool_t //--- Rotation subtraction AliDielectronSignalBase *sigTOFmixRot=GetSignalRot(d,stepTOFmix); SetStyle(sigTOFmixRot,"TOF + TPC - Track rotation subtraction"); -// DrawSpectra(sigTOFmixRot,"cTOFTPCrot",save); +// DrawSpectra(sigTOFmixRot,"cTOFTPCrot",hStats,save); if (hStats) delete hStats; } @@ -127,7 +136,7 @@ AliDielectronSignalBase *GetSignalRot(AliDielectronCFdraw &d, Int_t step) arr->AddAt(d.Project("M",step),iType); AliDielectronSignalExt *sig=new AliDielectronSignalExt; - sig->SetScaleRawToBackground(3.2,4.); +// sig->SetScaleRawToBackground(3.2,4.); sig->SetIntegralRange(2.93,3.15); sig->SetMethod(AliDielectronSignalBase::kRotation); sig->Process(arr); @@ -317,22 +326,22 @@ void DrawSpectra(AliDielectronSignalBase *sig, const char* cname, TH1 *hEventSt } if (save){ - c->SaveAs(Form("%s.eps",cname)); +// c->SaveAs(Form("%s.eps",cname)); c->SaveAs(Form("%s.png",cname)); - +/* FILE *out_file; if ( (out_file = fopen(Form("sig_%s.txt",cname), "w")) == NULL ) { fprintf(stderr, "Cannot open file %s\n", Form("sig_%s.txt",cname)); } fprintf(stdout, "Signal file: %s\n", Form("sig_%s.txt",cname)); fprintf(out_file,"%3d %4.1f %3.1f %4.2f %4.1f %4.2f %d\n",(int)sigN,sigEr,sigS2B,sigS2Ber,sigSignif,sigSignifEr,(Int_t)afterPhys); fclose(out_file); - + TFile outMinv(Form("Minv_%s.root",cname), "RECREATE"); hUS->Write(); hBackground->Write(); hSignal->Write(); - hMmc->Write(); - outMinv.Close(); + if (hMmc) hMmc->Write(); + outMinv.Close();*/ } diff --git a/PWG3/dielectron/macros/analyzeJpsiME.C b/PWG3/dielectron/macros/analyzeJpsiME.C new file mode 100644 index 00000000000..469308c8a29 --- /dev/null +++ b/PWG3/dielectron/macros/analyzeJpsiME.C @@ -0,0 +1,69 @@ +void analyzeJpsiME(TString tag="./"){ + + TStopwatch timer; + timer.Start(); + + //_____________Setting up libraries_______________ + gSystem->Load("libTree.so"); + gSystem->Load("libSTEERBase.so"); + gSystem->Load("libVMC.so"); + gSystem->Load("libESD.so"); + gSystem->Load("libANALYSIS.so"); + gSystem->Load("libANALYSISalice"); + gSystem->Load("libCORRFW"); + gSystem->Load("libPWG3dielectron.so"); + gSystem->AddIncludePath("-I$ALICE_ROOT/PWG3/dielectron/ -I$ALICE_ROOT/ANALYSIS -I$ALICE_ROOT/CORRFW"); + + //_____________Load Macro_____________________________ + gROOT->LoadMacro("ConfigJpsi2eeData.C"); // user Config file + + // Make the analysis manager__________________________ + AliAnalysisManager *mgr = new AliAnalysisManager("AnalysisManager"); + + // Event handler setting + Int_t nEventsToLoad = 5; // pool depth + AliMultiEventInputHandler *inputHandler = new AliMultiEventInputHandler(nEventsToLoad,0); // 0 for ESD, 1 for AOD + mgr->SetInputEventHandler(inputHandler); + + // event pool setting + AliEventPoolOTF *myPool = new AliEventPoolOTF("event pool","ESD"); + myPool->SetTagDirectory(tag.Data()); + myPool->SetMultiplicityBin(0,100,100); // (min, max, bin width) + myPool->Init(); + + mgr->SetEventPool(myPool); // link event pool with manager + inputHandler->SetEventPool(myPool); // link event pool with input handler + + // Analysis Task Jpsi->e+e-___________________________ + AliAnalysisTaskDielectronME *task = new AliAnalysisTaskDielectronME("DielectronTaskME"); + task->RequireFreshBuffer(); // refresh mode in the pool + task->SetPoolDepth(nEventsToLoad); // set pool depth in the task + mgr->AddTask(task); + + // add dielectron analysis with different cuts to the task + for (Int_t i=0; iAddDielectron(jpsi); + } + + // Create containers for input/output_________________ + AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("jpsi_treeEff",TTree::Class(),AliAnalysisManager::kExchangeContainer,"jpsi_Effdefault"); + AliAnalysisDataContainer *cOutputHist1 = mgr->CreateContainer("QA", TList::Class(), AliAnalysisManager::kOutputContainer, "jpsi_QA.root"); + AliAnalysisDataContainer *cOutputHist2 = mgr->CreateContainer("CF", TList::Class(), AliAnalysisManager::kOutputContainer, "jpsi_CF.root"); + AliAnalysisDataContainer *cOutputHist3 = mgr->CreateContainer("jpsi_Eff_EventStat",TH1D::Class(),AliAnalysisManager::kOutputContainer,"jpsi_Eff.root"); + + mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer()); + mgr->ConnectOutput(task, 0, coutput1 ); + mgr->ConnectOutput(task, 1, cOutputHist1); + mgr->ConnectOutput(task, 2, cOutputHist2); + mgr->ConnectOutput(task, 3, cOutputHist3); + + if (!mgr->InitAnalysis()) return; + mgr->PrintStatus(); + TChain *chain = NULL; // null pointer to avoid entering grid mode. temporary solution + mgr->StartAnalysis("mix",chain); + //mgr->Terminate(); + + timer.Stop(); + timer.Print(); +} diff --git a/PWG3/libPWG3dielectron.pkg b/PWG3/libPWG3dielectron.pkg index e43a85b46b3..8fc9646f28f 100644 --- a/PWG3/libPWG3dielectron.pkg +++ b/PWG3/libPWG3dielectron.pkg @@ -8,6 +8,7 @@ SRCS= dielectron/AliDielectron.cxx \ dielectron/AliDielectronMC.cxx \ dielectron/AliDielectronVarManager.cxx \ dielectron/AliAnalysisTaskDielectronSE.cxx \ + dielectron/AliAnalysisTaskDielectronME.cxx \ dielectron/AliAnalysisTaskDielectronFilter.cxx \ dielectron/AliAnalysisTaskDielectronEfficiency.cxx \ dielectron/AliAnalysisTaskMultiDielectron.cxx \ @@ -22,7 +23,8 @@ SRCS= dielectron/AliDielectron.cxx \ dielectron/AliDielectronTrackRotator.cxx \ dielectron/AliDielectronPID.cxx \ dielectron/AliDielectronCutGroup.cxx \ - dielectron/AliDielectronEventCuts.cxx + dielectron/AliDielectronEventCuts.cxx \ + dielectron/AliDielectronHelper.cxx HDRS= $(SRCS:.cxx=.h)