From 7b5556ef967c0aa38f411e019dbffc2fa9d3d2ca Mon Sep 17 00:00:00 2001 From: snelling Date: Wed, 16 Jan 2013 08:12:31 +0000 Subject: [PATCH] changes Paul --- PWG/CMakeLists.txt | 10 - PWG/CMakelibPWGTools.pkg | 7 + PWG/CMakelibPWGflowBase.pkg | 3 + PWG/CMakelibPWGflowTasks.pkg | 5 + PWG/CMakelibPWGmuon.pkg | 7 + PWG/FLOW/Base/AliFlowAnalysisWithMSP.cxx | 439 +++++++++++++----- PWG/FLOW/Base/AliFlowAnalysisWithMSP.h | 33 +- PWG/FLOW/Base/AliFlowEventSimple.cxx | 11 + PWG/FLOW/Base/AliFlowEventSimple.h | 1 + PWG/FLOW/Base/AliFlowMSPHistograms.cxx | 51 ++ PWG/FLOW/Base/AliFlowMSPHistograms.h | 3 + PWG/FLOW/Base/AliFlowVector.cxx | 43 +- PWG/FLOW/Base/AliFlowVector.h | 25 +- .../AliAnalysisTwoParticleResonanceFlowTask.h | 2 +- PWG/PWGflowBaseLinkDef.h | 2 + PWG/PWGflowTasksLinkDef.h | 1 + PWG/muon/AliDimuInfoStoreMC.h | 2 +- 17 files changed, 501 insertions(+), 144 deletions(-) diff --git a/PWG/CMakeLists.txt b/PWG/CMakeLists.txt index 90fd61fe57c..e69de29bb2d 100644 --- a/PWG/CMakeLists.txt +++ b/PWG/CMakeLists.txt @@ -1,10 +0,0 @@ -# AliRoot Build System CMakeLists for PWG -# -# Author: Johny Jose m(johny.jose@cern.ch) -# Port of previous Makefile build to cmake - -cmake_minimum_required(VERSION 2.8.4 FATAL_ERROR) - -file(GLOB PACKAGES CMake*.pkg) - -ALICE_BuildModule() diff --git a/PWG/CMakelibPWGTools.pkg b/PWG/CMakelibPWGTools.pkg index 4b111d82dd0..54b84e8733c 100644 --- a/PWG/CMakelibPWGTools.pkg +++ b/PWG/CMakelibPWGTools.pkg @@ -42,3 +42,10 @@ set ( DHDR PWGToolsLinkDef.h) string ( REPLACE ".cxx" ".h" EXPORT "${SRCS}" ) set ( EINCLUDE PWG/Tools ANALYSIS CORRFW STEER/AOD STEER/ESD STEER/STEERBase) + +if( ALICE_TARGET STREQUAL "win32gcc") + set ( PACKSOFLAGS ${SOFLAGS} ) + set ( ELIBS --verbose -L${ALICE_ROOT}/lib/tgt_${ALICE_TARGET} -lANALYSISalice -lANALYSIS -lCORRFW -lAOD -lESD -lSTEER -lSTEERBase ${ROOTCLIBS} -lEG -lGraf -lGraf3d -lASImage ) +endif( ALICE_TARGET STREQUAL "win32gcc") + + diff --git a/PWG/CMakelibPWGflowBase.pkg b/PWG/CMakelibPWGflowBase.pkg index c5e7e2f8f11..878024ce7ea 100644 --- a/PWG/CMakelibPWGflowBase.pkg +++ b/PWG/CMakelibPWGflowBase.pkg @@ -53,6 +53,8 @@ set ( SRCS FLOW/Base/AliFlowAnalysisWithFittingQDistribution.cxx FLOW/Base/AliFlowAnalysisWithMixedHarmonics.cxx FLOW/Base/AliFlowAnalysisWithNestedLoops.cxx + FLOW/Base/AliFlowMSPHistograms.cxx + FLOW/Base/AliFlowAnalysisWithMSP.cxx FLOW/Base/AliFlowOnTheFlyEventGenerator.cxx ) @@ -67,5 +69,6 @@ set ( EINCLUDE PWG/FLOW/Base) if( ALICE_TARGET STREQUAL "win32gcc") set ( PACKSOFLAGS ${SOFLAGS} -L${ROOTLIBDIR} -lEG) + set ( ELIBS ${ROOTCLIBS} -lEG) endif( ALICE_TARGET STREQUAL "win32gcc") diff --git a/PWG/CMakelibPWGflowTasks.pkg b/PWG/CMakelibPWGflowTasks.pkg index d3601ddc240..d562b297625 100644 --- a/PWG/CMakelibPWGflowTasks.pkg +++ b/PWG/CMakelibPWGflowTasks.pkg @@ -34,6 +34,7 @@ set ( SRCS FLOW/Tasks/AliFlowCandidateTrack.cxx FLOW/Tasks/AliFlowTrackCuts.cxx FLOW/Tasks/AliAnalysisTaskScalarProduct.cxx + FLOW/Tasks/AliAnalysisTaskMSP.cxx FLOW/Tasks/AliAnalysisTaskMCEventPlane.cxx FLOW/Tasks/AliAnalysisTaskLYZEventPlane.cxx FLOW/Tasks/AliAnalysisTaskCumulants.cxx @@ -64,4 +65,8 @@ set ( DHDR PWGflowTasksLinkDef.h) set ( EXPORT ${HDRS} ) +if( ALICE_TARGET STREQUAL "win32gcc") + set ( ELIBS -L${ALICE_ROOT}/lib/tgt_${ALICE_TARGET} -lPWGflowBase -lPWGmuon -lANALYSISalice -lANALYSIS -lAOD -lESD -lSTEERBase ${ROOTCLIBS} -lEG ) +endif( ALICE_TARGET STREQUAL "win32gcc") + set ( EINCLUDE PWG/FLOW/Base PWG/FLOW/Tasks STEER/ESD STEER/STEERBase PWG/muon) diff --git a/PWG/CMakelibPWGmuon.pkg b/PWG/CMakelibPWGmuon.pkg index 49ce97379c2..00a67be696b 100644 --- a/PWG/CMakelibPWGmuon.pkg +++ b/PWG/CMakelibPWGmuon.pkg @@ -72,3 +72,10 @@ string ( REPLACE ".cxx" ".h" EXPORT "${SRCS}" ) set ( DHDR PWGmuonLinkDef.h) set ( EINCLUDE PWG/muon ANALYSIS STEER/AOD STEER/ESD STEER/STEERBase) + +if( ALICE_TARGET STREQUAL "win32gcc") + set ( PACKSOFLAGS ${SOFLAGS} ) + set ( ELIBS -L${ALICE_ROOT}/lib/tgt_${ALICE_TARGET} -lANALYSIS -lANALYSISalice -lCORRFW -lAOD -lSTEERBase -lESD ${ROOTCLIBS} -lEG) +endif( ALICE_TARGET STREQUAL "win32gcc") + + diff --git a/PWG/FLOW/Base/AliFlowAnalysisWithMSP.cxx b/PWG/FLOW/Base/AliFlowAnalysisWithMSP.cxx index 84680a35e0e..41e47d6b81f 100644 --- a/PWG/FLOW/Base/AliFlowAnalysisWithMSP.cxx +++ b/PWG/FLOW/Base/AliFlowAnalysisWithMSP.cxx @@ -34,36 +34,87 @@ #include AliFlowAnalysisWithMSP::AliFlowAnalysisWithMSP() - : TNamed(), fHarmonic(2), fNUA(kFALSE), fBookCommonHistograms(kFALSE), fCommonHist(0), fQaComponents(0), + : TNamed(), + fHarmonic(2), fNUA(kFALSE), fUseCommonConstants(kFALSE), fBookCommonHistograms(kFALSE), fCommonHist(0), fQaComponents(0), fQbComponents(0), fQaQb(0), fPtUComponents(0), fEtaUComponents(0), fAllStatistics(0), fPtStatistics(0), fEtaStatistics(0), - fIntegratedFlow(0), fDiffFlowPt(0), fDiffFlowEta(0), fFlags(0) + fIntegratedFlow(0), fDiffFlowPt(0), fDiffFlowEta(0), fFlags(0), fHistList(0) { // Default constructor. Intended for root IO purposes only } -AliFlowAnalysisWithMSP::AliFlowAnalysisWithMSP(TDirectoryFile *file) - : TNamed(), fHarmonic(2), fNUA(kFALSE), fUseCommonConstants(kFALSE), fBookCommonHistograms(kFALSE), fCommonHist(0), fQaComponents(0), +AliFlowAnalysisWithMSP::AliFlowAnalysisWithMSP(TDirectory *file) + : TNamed(), + fHarmonic(2), fNUA(kFALSE), fUseCommonConstants(kFALSE), fBookCommonHistograms(kFALSE), fCommonHist(0), fQaComponents(0), fQbComponents(0), fQaQb(0), fPtUComponents(0), fEtaUComponents(0), fAllStatistics(0), fPtStatistics(0), fEtaStatistics(0), - fIntegratedFlow(0), fDiffFlowPt(0), fDiffFlowEta(0), fFlags(0) + fIntegratedFlow(0), fDiffFlowPt(0), fDiffFlowEta(0), fFlags(0), fHistList(0) { + // Constructor reads the internal state from a root file. + // No check for consistency is done. If flags or histograms are not found they are left at default (0). ReadHistograms(file); } +AliFlowAnalysisWithMSP::AliFlowAnalysisWithMSP(TList *list) + : TNamed(), + fHarmonic(2), fNUA(kFALSE), fUseCommonConstants(kFALSE), fBookCommonHistograms(kFALSE), fCommonHist(0), fQaComponents(0), + fQbComponents(0), fQaQb(0), fPtUComponents(0), fEtaUComponents(0), fAllStatistics(0), + fPtStatistics(0), fEtaStatistics(0), + fIntegratedFlow(0), fDiffFlowPt(0), fDiffFlowEta(0), fFlags(0), fHistList(0) +{ + // Constructor reads the internal state from a root file. + // No check for consistency is done. If flags or histograms are not found they are left at default (0). + ReadHistograms(list); +} + + AliFlowAnalysisWithMSP::AliFlowAnalysisWithMSP(const unsigned int harmonic, const bool commonConst, const bool commonHist) - : TNamed(), fHarmonic(harmonic), fNUA(kFALSE), fUseCommonConstants(commonConst), fBookCommonHistograms(commonHist), fCommonHist(0), fQaComponents(0), + : TNamed(), + fHarmonic(harmonic), fNUA(kFALSE), fUseCommonConstants(commonConst), fBookCommonHistograms(commonHist), fCommonHist(0), fQaComponents(0), fQbComponents(0), fQaQb(0), fPtUComponents(0), fEtaUComponents(0), fAllStatistics(0), fPtStatistics(0), fEtaStatistics(0), - fIntegratedFlow(0), fDiffFlowPt(0), fDiffFlowEta(0), fFlags(0) + fIntegratedFlow(0), fDiffFlowPt(0), fDiffFlowEta(0), fFlags(0), fHistList(0) { // Constructor defining the harmonic, usage of non uniform acceptance corrections and the AliFlowCommonHist() histograms + // Calls Init() to set-up the histograms. This is equivalent to: + // AliFlowAnalysisWithMSP *ana= new AliFlowAnalysisWithMSP); + // ana->SetHarmonic(harmonic); + // ana->UseCommonConstants(commonConst); + // ana->EnableCommonHistograms(commonHist); + // ana->Init(); // This is the constructor intended for the user SetNameTitle("MSP","Flow analysis with the Modified Scalar Product method"); + Init(); +} + +AliFlowAnalysisWithMSP::AliFlowAnalysisWithMSP(const AliFlowAnalysisWithMSP &x) + : TNamed(), + fHarmonic(x.fHarmonic), fNUA(x.fNUA), fUseCommonConstants(x.fUseCommonConstants), fBookCommonHistograms(x.fBookCommonHistograms), fCommonHist(0), + fQaComponents(0), fQbComponents(0), fQaQb(0), fPtUComponents(0), fEtaUComponents(0), fAllStatistics(0), + fPtStatistics(0), fEtaStatistics(0), + fIntegratedFlow(0), fDiffFlowPt(0), fDiffFlowEta(0), fFlags(0), fHistList(0) +{ + // Copy constructor + SetNameTitle("MSP","Flow analysis with the Modified Scalar Product method"); + if( x.fQaComponents ) fQaComponents=(AliFlowMSPHistograms *)(x.fQaComponents)->Clone(); + if( x.fQbComponents ) fQbComponents=(AliFlowMSPHistograms *)(x.fQbComponents)->Clone(); + if( x.fQaQb ) fQaQb=(AliFlowMSPHistograms *)(x.fQaQb)->Clone(); + if( x.fPtUComponents) fPtUComponents=(AliFlowMSPHistograms *)(x.fPtUComponents)->Clone(); + if( x.fEtaUComponents ) fEtaUComponents=(AliFlowMSPHistograms *)(x.fEtaUComponents)->Clone(); + if( x.fAllStatistics ) fAllStatistics=(AliFlowMSPHistograms *)(x.fAllStatistics)->Clone(); + if( x.fPtStatistics ) fPtStatistics=(AliFlowMSPHistograms *)(x.fPtStatistics)->Clone(); + if( x.fEtaStatistics ) fEtaStatistics=(AliFlowMSPHistograms *)(x.fEtaStatistics)->Clone(); + if( x.fIntegratedFlow ) fIntegratedFlow=(TH1D *)(x.fIntegratedFlow)->Clone(); + if( x.fDiffFlowPt ) fDiffFlowPt=(TH1D *)(x.fDiffFlowPt)->Clone(); + if( x.fDiffFlowEta ) fDiffFlowEta=(TH1D *)(x.fDiffFlowEta)->Clone(); + if( x.fFlags ) fFlags=(TProfile *)(x.fFlags)->Clone(); + if( x.fCommonHist ) fCommonHist=new AliFlowCommonHist(*(x.fCommonHist)); + // The histogram list fHistList is not cloned because it is regenerated when requested } AliFlowAnalysisWithMSP::~AliFlowAnalysisWithMSP() { + // Destructor. All internal objects are owned by this class and deleted here. delete fCommonHist; delete fQaComponents; delete fQbComponents; @@ -77,12 +128,16 @@ AliFlowAnalysisWithMSP::~AliFlowAnalysisWithMSP() delete fDiffFlowPt; delete fDiffFlowEta; delete fFlags; + if(fHistList) { + fHistList->SetOwner(kFALSE); // Histograms were already deleted manually + delete fHistList; // Delete the TList itself + } } void AliFlowAnalysisWithMSP::Init() { // Create all output objects. Memory consumption can be reduced by switching off some of the control histograms. - // + // pt and eta binning are set to the values defined in the static class AliFlowCommonConstants if enabled. Otherwise defaults are set. delete fCommonHist; fCommonHist=0; // Delete existing histograms delete fQaComponents; fQaComponents=0; delete fQbComponents; fQbComponents=0; @@ -96,9 +151,14 @@ void AliFlowAnalysisWithMSP::Init() delete fDiffFlowPt; fDiffFlowPt=0; delete fDiffFlowEta; fDiffFlowEta=0; delete fFlags; fFlags=0; + if(fHistList) { + fHistList->SetOwner(kFALSE); // Histograms were already deleted manually + delete fHistList; // Delete the TList itself + fHistList=0; // Clear pointer which is invalid afcter delete + } // Default binning for histograms - // TODO: allow variable binning + // TODO: allow variable binning in a simpler way than by AliFlowCommonConstants() only // Defaults int nBinsPt=100; @@ -182,7 +242,7 @@ void AliFlowAnalysisWithMSP::Make(AliFlowEventSimple *event) if( !event ) return; // Protect against running without events. Can this happen??? AliFlowVector flowVectors[2]; // Calculate the two subevent Q vectors: - event->Get2Qsub(flowVectors, fHarmonic); // No phi, pt or eta weights implemented yet + event->Get2Qsub(flowVectors, fHarmonic); // TODO: Check how do the phi, pt, eta weights work in the flow event? AliFlowVector &Qa=flowVectors[0]; // Define some mnemonics for the subevent flow vectors: AliFlowVector &Qb=flowVectors[1]; @@ -194,60 +254,57 @@ void AliFlowAnalysisWithMSP::Make(AliFlowEventSimple *event) if(fCommonHist)fCommonHist->FillControlHistograms(event); // Standard for all flow analysis - - const double qaxy[]={Qa.X()/QaW,Qa.Y()/QaW}; // Two variables expected + // Fill NUA correction histograms for Qa, Qb and QaQb per event: + const double qaxy[]={Qa.X()/QaW,Qa.Y()/QaW}; // Two variables: Qa components const double wqaxy[]={QaW,QaW}; fQaComponents->Fill(1, qaxy, wqaxy); // only one bin (all pt) - const double qbxy[]={Qb.X()/QbW,Qb.Y()/QbW}; // Two variables expected + const double qbxy[]={Qb.X()/QbW,Qb.Y()/QbW}; // Two variables: Qb components const double wqbxy[]={QbW,QbW}; fQbComponents->Fill(1, qbxy, wqbxy); // only one bin (all pt) const double QaQbW=QaW*QbW; - const double weightedQaQb = (Qa*Qb)/QaQbW; // Scalar product of subevent Q vectors with weight + const double weightedQaQb = (Qa*Qb)/QaQbW; // Scalar product of subevent Q vectors with weight, only 1 variable fQaQb->Fill(1,&weightedQaQb,&QaQbW); // Average of QaQb per event - int iTrack=0; while( AliFlowTrackSimple *track=event->GetTrack(iTrack++) ) {// Loop over the tracks in the event - // Get the track vector - if(! track->InPOISelection() ) continue; - const double trackWeight=track->Weight(); + if(! track->InPOISelection() ) continue; // Ignore if not a POI + const double trackWeight=track->Weight(); // Get the track vector const double phi=track->Phi(); const double pt=track->Pt(); const double eta=track->Eta(); AliFlowVector u; - u.SetMagPhi(1, fHarmonic*phi, trackWeight); - u.SetMult(1); + u.SetMagPhi(trackWeight, fHarmonic*phi, trackWeight); // Length = track weight = multiplicity for a single track // Remove track from subevent a - AliFlowVector mQa(Qa); // Initialize with Qa flow vector - if( track->InSubevent(0) ) { - mQa-=u; // Should introduce phi weights here + AliFlowVector mQa(Qa); // Initialize with Qa flow vector of this event + if( track->InSubevent(0) ) { // Correct for autocorrelation with POI for this track + mQa-=u; } // Remove track from subevent b - AliFlowVector mQb(Qb); // Initialize with Qb flow vector - if( track->InSubevent(1) ) { - mQb-=u; // Should introduce phi weights here + AliFlowVector mQb(Qb); // Initialize with Qb flow vector of this event + if( track->InSubevent(1) ) { // Correct for autocorrelation with POI for this track + mQb-=u; } - const double uQaW = mQa.GetMult(); // Weight is multiplicity of Q vector - const double uQbW = mQb.GetMult(); - const double uQa=u*mQa; // Correlate POI with subevent + const double uQaW = mQa.GetMult()*trackWeight; // Weight is multiplicity of Q vector times weight of POI + const double uQbW = mQb.GetMult()*trackWeight; + const double uQa=u*mQa; // Correlate POI with subevent const double uQb=u*mQb; const double uxy[]={u.X(),u.Y()}; const double wxy[]={1.0,1.0}; - fPtUComponents->Fill(pt, uxy, wxy); // vs pt - fEtaUComponents->Fill(eta, uxy, wxy); // vs eta + fPtUComponents->Fill(pt, uxy, wxy); // NUA for POI vs pt + fEtaUComponents->Fill(eta, uxy, wxy); // NUA for POI vs eta - const double par[]={uQa/uQaW, uQb/uQbW, weightedQaQb}; + const double par[]={uQa/uQaW, uQb/uQbW, weightedQaQb}; // Three variables to correlate: uQa, uQb and QaQb const double wgt[]={uQaW, uQbW, QaQbW}; - fAllStatistics->Fill(1, par, wgt ); - fPtStatistics->Fill(pt, par, wgt ); - fEtaStatistics->Fill(eta, par, wgt ); + fAllStatistics->Fill(1, par, wgt ); // only 1 bin, integrated over pt and eta + fPtStatistics->Fill(pt, par, wgt ); // pt differential correlations + fEtaStatistics->Fill(eta, par, wgt ); // eta differential correlations } } @@ -258,7 +315,6 @@ void AliFlowAnalysisWithMSP::Finish() // If the output histograms already exist then they are replaced by the newly calculated result Print(); // Print a summary of the NUA terms and integrated flow - // TODO: Create result histograms for integrated flow and store the flags to be able to restore the initial state from histograms only // Create result histograms fFlags->Fill("NUA",fNUA); @@ -266,7 +322,11 @@ void AliFlowAnalysisWithMSP::Finish() delete fDiffFlowPt; fDiffFlowPt=0; delete fDiffFlowEta; fDiffFlowEta=0; + bool oldAddStatus=TH1::AddDirectoryStatus(); // Do not store the next histograms automatically + TH1::AddDirectory(kFALSE); // We need full control over the writing of the hostograms + fIntegratedFlow=new TH1D("IntegratedFlow","Integrated flow results",10,0.5,10.5); + fIntegratedFlow->SetDirectory(0); double vn, vnerror; Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 0); fIntegratedFlow->SetBinContent(1,vn); @@ -291,11 +351,11 @@ void AliFlowAnalysisWithMSP::Finish() fDiffFlowPt->SetYTitle(Form("v_{%d}",fHarmonic)); for(int i=1; i<=nbinspt; ++i) { - double vn=0; - double evn=0; - if( Calculate(vn, evn, fPtStatistics, fPtUComponents, i) && evn<1 ) { - fDiffFlowPt->SetBinContent(i,vn); - fDiffFlowPt->SetBinError(i,evn); + double vnpt=0; + double evnpt=0; + if( Calculate(vnpt, evnpt, fPtStatistics, fPtUComponents, i) && evnpt<1 ) { + fDiffFlowPt->SetBinContent(i,vnpt); + fDiffFlowPt->SetBinError(i,evnpt); } } @@ -305,33 +365,42 @@ void AliFlowAnalysisWithMSP::Finish() fDiffFlowEta= new TH1D("DiffFlowEta","flow of POI vs #eta",nbinseta,etalow,etahigh); fDiffFlowEta->SetDirectory(0); // Do not automatically store in file fDiffFlowEta->SetStats(kFALSE); - fDiffFlowEta->SetXTitle("p_{t}"); + fDiffFlowEta->SetXTitle("#eta"); fDiffFlowEta->SetYTitle(Form("v_{%d}",fHarmonic)); for(int i=1; i<=nbinseta; ++i) { - double vn=0; - double evn=0; - if( Calculate(vn, evn, fEtaStatistics, fEtaUComponents, i) && evn<1 ) { - fDiffFlowEta->SetBinContent(i,vn); - fDiffFlowEta->SetBinError(i,evn); + double vneta=0; + double evneta=0; + if( Calculate(vneta, evneta, fEtaStatistics, fEtaUComponents, i) && evneta<1 ) { + fDiffFlowEta->SetBinContent(i,vneta); + fDiffFlowEta->SetBinError(i,evneta); } } + + TH1::AddDirectory(oldAddStatus); } void AliFlowAnalysisWithMSP::WriteHistograms(TDirectoryFile *file) const { - //file->Write(file->GetName(), TObject::kSingleKey); // Make sure the directory itself is written + // Write the internal status to file. If the AliFlowCommonHistograms were enabled they are written also. + // Results are written only if they were calculated by Finish(). + // From these histograms the internal state of *this can be reconstructed with ReadHistograms() + + //file->Write(file->GetName(), TObject::kSingleKey); // Make sure the directory itself is written if(fCommonHist) file->WriteTObject(fCommonHist); - file->WriteTObject(fQaComponents,0,"Overwrite"); // Averages of Qa components per event - file->WriteTObject(fQbComponents,0,"Overwrite"); // Averages of Qb components per event - file->WriteTObject(fQaQb,0,"Overwrite"); // Average of QaQb per event - file->WriteTObject(fPtUComponents,0,"Overwrite"); // u components vs pt - file->WriteTObject(fEtaUComponents,0,"Overwrite"); // u components vs eta - file->WriteTObject(fAllStatistics,0,"Overwrite"); // Integrated uQa, uQb and QaQa - file->WriteTObject(fPtStatistics,0,"Overwrite"); // uQa, uQb and QaQb vs pt - file->WriteTObject(fEtaStatistics,0,"Overwrite"); // uQa, uQb and QaQb vs eta + TList *t=new TList(); // This object is written as a marker for redoFinish() + t->SetName("cobjMSP"); + file->WriteTObject(t,0,"Overwrite"); + file->WriteTObject(fQaComponents,0,"Overwrite"); // Averages of Qa components per event + file->WriteTObject(fQbComponents,0,"Overwrite"); // Averages of Qb components per event + file->WriteTObject(fQaQb,0,"Overwrite"); // Average of QaQb per event + file->WriteTObject(fPtUComponents,0,"Overwrite"); // u components vs pt + file->WriteTObject(fEtaUComponents,0,"Overwrite"); // u components vs eta + file->WriteTObject(fAllStatistics,0,"Overwrite"); // Integrated uQa, uQb and QaQa + file->WriteTObject(fPtStatistics,0,"Overwrite"); // uQa, uQb and QaQb vs pt + file->WriteTObject(fEtaStatistics,0,"Overwrite"); // uQa, uQb and QaQb vs eta if( fIntegratedFlow ) file->WriteTObject(fIntegratedFlow,0,"Overwrite"); // Integrated flow for POI and subevents if( fDiffFlowPt ) file->WriteTObject(fDiffFlowPt,0,"Overwrite"); // Differential flow vs pt if calculated @@ -344,12 +413,10 @@ void AliFlowAnalysisWithMSP::WriteHistograms(TDirectoryFile *file) const void AliFlowAnalysisWithMSP::WriteCommonResults(TDirectoryFile *file) const { - // Copy the results to a AliFlowCommonHistResults() class and write to file - // If the results were not calculated again then Finish() is called to generate them - // The global AliFlowCommonConstants() is modified to adapt to the binning used in this analysis - - AliFlowCommonConstants *c=AliFlowCommonConstants::GetMaster(); - int nBinsPt=fPtStatistics->Nbins(); + // Export the results to a AliFlowCommonHistResults() class and write to file + // If the results were not calculated then Finish() is called to generate them + + int nBinsPt=fPtStatistics->Nbins(); // Get the actually used binning from the Statistics histograms double ptMin=fPtStatistics->XLow(); double ptMax=fPtStatistics->XHigh(); @@ -357,6 +424,15 @@ void AliFlowAnalysisWithMSP::WriteCommonResults(TDirectoryFile *file) const double etaMin=fEtaStatistics->XLow(); double etaMax=fEtaStatistics->XHigh(); + AliFlowCommonConstants *c=AliFlowCommonConstants::GetMaster(); // Get the static common constants object + // Save the old AliFlowCommonConstants status + int oldNbinsPt=c->GetNbinsPt(); + double oldPtMin=c->GetPtMin(); + double oldPtMax=c->GetPtMax(); + int oldNbinsEta=c->GetNbinsEta(); + double oldEtaMin=c->GetEtaMin(); + double oldEtaMax=c->GetEtaMax(); + // Modify AliFlowCommonConstants to make sure that AliFlowCommonResults is generated with the correct binning c->SetNbinsPt(nBinsPt); c->SetPtMin(ptMin); c->SetPtMax(ptMax); @@ -364,8 +440,8 @@ void AliFlowAnalysisWithMSP::WriteCommonResults(TDirectoryFile *file) const c->SetEtaMin(etaMin); c->SetEtaMax(etaMax); - bool oldAddStatus=TH1::AddDirectoryStatus(); - TH1::AddDirectory(kFALSE); + bool oldAddStatus=TH1::AddDirectoryStatus(); // Do not store the next histograms automatically + TH1::AddDirectory(kFALSE); // We need full control over the writing of the hostograms AliFlowCommonHistResults *h=new AliFlowCommonHistResults("AliFlowCommonHistResults_MSP","AliFlowCommonHistResults from the MSP method",fHarmonic); double ivn, ivnerror; @@ -390,15 +466,52 @@ void AliFlowAnalysisWithMSP::WriteCommonResults(TDirectoryFile *file) const file->WriteTObject(h,0,"Overwrite"); - TH1::AddDirectory(oldAddStatus); + TH1::AddDirectory(oldAddStatus); // Restore the automatic storage of histograms to its original status + + // Restore AliFlowCommonConstants to make sure that no other analysis are affected + c->SetNbinsPt(oldNbinsPt); + c->SetPtMin(oldPtMin); + c->SetPtMax(oldPtMax); + c->SetNbinsEta(oldNbinsEta); + c->SetEtaMin(oldEtaMin); + c->SetEtaMax(oldEtaMax); } +TList *AliFlowAnalysisWithMSP::ListHistograms() +{ + if( fHistList ) { + fHistList->SetOwner(kFALSE); + delete fHistList; + } + fHistList = new TList(); + fHistList->SetOwner(kFALSE); + + if(fCommonHist) fHistList->Add(fCommonHist); // Standard control histograms, if enabled + + //Correlations + if(fQaComponents) fHistList->Add(fQaComponents); // Averages of Qa components per event for NUA + if(fQbComponents) fHistList->Add(fQbComponents); // Averages of Qb components per event for NUA + if(fQaQb) fHistList->Add(fQaQb); // Average of QaQb per event + if(fPtUComponents) fHistList->Add(fPtUComponents); // ux and uy per pt bin for NUA + if(fEtaUComponents) fHistList->Add(fEtaUComponents); // ux and uy per eta bin for NUA + if(fAllStatistics) fHistList->Add(fAllStatistics); // Correlations for uQa uQb and QaQb (integrated) + if(fPtStatistics) fHistList->Add(fPtStatistics); // Correlations for uQa uQb and QaQb per pt bin + if(fEtaStatistics) fHistList->Add(fEtaStatistics); // Correlations for uQa uQb and QaQb per eta bin + + // Result histograms (if calculated) + if(fIntegratedFlow) fHistList->Add(fIntegratedFlow); // vn for POI and subevents + if(fDiffFlowPt) fHistList->Add(fDiffFlowPt); // vn as function of pt + if(fDiffFlowEta) fHistList->Add(fDiffFlowEta); // vn as function of eta + if(fFlags) fHistList->Add(fFlags); // Stores fHarmonic and fNUA + + return fHistList; +} void AliFlowAnalysisWithMSP::Print(const Option_t *opt)const { - std::cout << setprecision(3); - + if( opt ) std::cout << std::endl; + std::cout << "****************************************************" << std::endl; std::cout << " Integrated flow from Modified Scalar Product " << std::endl; std::cout << " " << std::endl; @@ -406,16 +519,17 @@ void AliFlowAnalysisWithMSP::Print(const Option_t *opt)const double vn=0; double vnerror=0; + std::cout << setprecision(4); Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 0); - std::cout << "v" << fHarmonic << " for POI : " << setw(10) << vn << " +- " << setw(8) << vnerror << std::endl; + std::cout << "v" << fHarmonic << " for POI : " << setw(11) << vn << " +- " << setw(9) << vnerror << std::endl; Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 1); - std::cout << "v" << fHarmonic << " for subevent A: " << setw(10) << vn << " +- " << setw(8) << vnerror << std::endl; + std::cout << "v" << fHarmonic << " for subevent A: " << setw(11) << vn << " +- " << setw(9) << vnerror << std::endl; Calculate(vn, vnerror, fAllStatistics, fPtUComponents, 0, 2); - std::cout << "v" << fHarmonic << " for subevent B: " << setw(10) << vn << " +- " << setw(8) << vnerror << std::endl; + std::cout << "v" << fHarmonic << " for subevent B: " << setw(11) << vn << " +- " << setw(9) << vnerror << std::endl; std::cout << std::endl; std::cout << "NUA terms: " << (fNUA?"(applied)":"(NOT applied)") << std::endl; - + std::cout << setprecision(3); const double ux=fPtUComponents->Average(0); // Average over all bins const double eux=TMath::Sqrt(fPtUComponents->Variance(0)); std::cout << " " << setw(12) << ux << " +- " << setw(12) << eux << (TMath::Abs(ux)<2*eux?" NOT significant ":" ") << std::endl; @@ -436,11 +550,21 @@ void AliFlowAnalysisWithMSP::Print(const Option_t *opt)const const double euy0pt=TMath::Sqrt(fPtUComponents->Variance(fPtUComponents->FindBin(1.0),1)); std::cout << " pt=1 " << setw(12) << uy0pt << " +- " << setw(12) << euy0pt << (TMath::Abs(uy0pt)<2*euy0pt?" NOT significant ":" ") << std::endl; - std::cout << " " << setw(12) << fQaComponents->Average(0) << " +- " << setw(12) << TMath::Sqrt(fQaComponents->Variance(0)) << std::endl; - std::cout << " " << setw(12) << fQaComponents->Average(1) << " +- " << setw(12) << TMath::Sqrt(fQaComponents->Variance(1)) << std::endl; - std::cout << " " << setw(12) << fQbComponents->Average(0) << " +- " << setw(12) << TMath::Sqrt(fQbComponents->Variance(0)) << std::endl; - std::cout << " " << setw(12) << fQbComponents->Average(1) << " +- " << setw(12) << TMath::Sqrt(fQbComponents->Variance(1)) << std::endl; - std::cout << " " << setw(12) << fQaQb->Average(0) << " +- " << setw(12) << TMath::Sqrt(fQbComponents->Variance(0)) << std::endl; + const double ax=fQaComponents->Average(0); + const double eax=TMath::Sqrt(fQaComponents->Variance(0)); + std::cout << " " << setw(12) << ax << " +- " << setw(12) <Average(1); + const double eay=TMath::Sqrt(fQaComponents->Variance(1)); + std::cout << " " << setw(12) << ay << " +- " << setw(12) << eay << (TMath::Abs(ay)<2*eay?" NOT significant ":" ") << std::endl; + const double bx=fQbComponents->Average(0); + const double ebx=TMath::Sqrt(fQbComponents->Variance(0)); + std::cout << " " << setw(12) << bx << " +- " << setw(12) << ebx << (TMath::Abs(bx)<2*ebx?" NOT significant ":" ") << std::endl; + const double by=fQbComponents->Average(1); + const double eby=TMath::Sqrt(fQbComponents->Variance(1)); + std::cout << " " << setw(12) << by << " +- " << setw(12) << eby << (TMath::Abs(by)<2*eby?" NOT significant ":" ") << std::endl; + const double ab=fQaQb->Average(0); + const double eab=TMath::Sqrt(fQbComponents->Variance(0)); + std::cout << " " << setw(12) << ab << " +- " << setw(12) << eab << (TMath::Abs(ab)<2*eab?" NOT significant ":" ") << std::endl; std::cout << std::endl; std::cout << "Covariance matrix: " << std::endl; @@ -453,6 +577,38 @@ void AliFlowAnalysisWithMSP::Print(const Option_t *opt)const } +AliFlowAnalysisWithMSP &AliFlowAnalysisWithMSP::operator=(const AliFlowAnalysisWithMSP &x) +{ + SetNameTitle("MSP","Flow analysis with the Modified Scalar Product method"); + delete fQaComponents; fQaComponents=0; + if( x.fQaComponents ) fQaComponents=(AliFlowMSPHistograms *)(x.fQaComponents)->Clone(); + delete fQbComponents; fQbComponents=0; + if( x.fQbComponents ) fQbComponents=(AliFlowMSPHistograms *)(x.fQbComponents)->Clone(); + delete fQaQb; fQaQb=0; + if( x.fQaQb ) fQaQb=(AliFlowMSPHistograms *)(x.fQaQb)->Clone(); + delete fPtUComponents; fPtUComponents=0; + if( fPtUComponents) fPtUComponents=(AliFlowMSPHistograms *)(x.fPtUComponents)->Clone(); + delete fEtaUComponents; fEtaUComponents=0; + if( fEtaUComponents ) fEtaUComponents=(AliFlowMSPHistograms *)(x.fEtaUComponents)->Clone(); + delete fAllStatistics; fAllStatistics=0; + if( fAllStatistics ) fAllStatistics=(AliFlowMSPHistograms *)(x.fAllStatistics)->Clone(); + delete fPtStatistics; fPtStatistics=0; + if( fPtStatistics ) fPtStatistics=(AliFlowMSPHistograms *)(x.fPtStatistics)->Clone(); + delete fEtaStatistics; fEtaStatistics=0; + if( fEtaStatistics ) fEtaStatistics=(AliFlowMSPHistograms *)(x.fEtaStatistics)->Clone(); + delete fIntegratedFlow; fIntegratedFlow=0; + if( fIntegratedFlow ) fIntegratedFlow=(TH1D *)(x.fIntegratedFlow)->Clone(); + delete fDiffFlowPt; fDiffFlowPt=0; + if( fDiffFlowPt ) fDiffFlowPt=(TH1D *)(x.fDiffFlowPt)->Clone(); + delete fDiffFlowEta; fDiffFlowEta=0; + if( fDiffFlowEta ) fDiffFlowEta=(TH1D *)(x.fDiffFlowEta)->Clone(); + delete fFlags; fFlags=0; + if( fFlags ) fFlags=(TProfile *)(x.fFlags)->Clone(); + delete fCommonHist; fCommonHist=0; + if( fCommonHist ) fCommonHist=new AliFlowCommonHist(*(x.fCommonHist)); + return *this; +} + // private functions -------------------------------------------------------------------------------------- bool AliFlowAnalysisWithMSP::Calculate(double &vn, double &vnerror, const AliFlowMSPHistograms *hist, const AliFlowMSPHistograms *components, const int bin, const int poi) const { @@ -461,7 +617,7 @@ bool AliFlowAnalysisWithMSP::Calculate(double &vn, double &vnerror, const AliFlo double VuQa=hist->Variance(bin,0); // Var(<>) double uQb=hist->Average(bin,1); // <> double VuQb=hist->Variance(bin,1); // Var(<>) - double QaQb=fQaQb->Average(1,0); // Should not be taken from hist(bin) buit from fQaQa because there QaQb is entered multiple times: <>!! + double QaQb=fQaQb->Average(1,0); // Should not be taken from hist(bin) because there QaQb is entered multiple times: <>!! double VQaQb=fQaQb->Variance(1,0); // V() double CuQauQb=hist->Covariance(bin,0,1); // Cov(<>,<>) double CuQaQaQb=hist->Covariance(bin,0,2); // Cov(<>,) @@ -484,19 +640,19 @@ bool AliFlowAnalysisWithMSP::Calculate(double &vn, double &vnerror, const AliFlo uQb = uQb - uX*QbX - uY*QbY; // Error calculation not fully modified: only above terms, the spread in <> and and are not used // therefore this should only be applied if significant! - // TODO: Check if not fully NUA correcting the error calculation is justified (but this is compatible with the original SP method)! + // Check if not fully NUA correcting the error calculation is justified (but this is compatible with the original SP method)! + // In general it is not justified but this can only be checked by splitting the event sample in many subsamples and looking at + // the variance of the result. This may not be feasible if statistics is low } // Some sanity checks: - if( uQa*uQb*QaQb <= 0 ) { // Catch imaginary results + if( uQa*uQb*QaQb <= 0 ) { // Catch imaginary results vn=-99; vnerror=-99; return false; } - // Decent results, calculate, print and store - // TODO: The three cases are cyclic permutations of u, Qa, Qb so there should be a simpler way to do this. - // However the difficulty is that we deal here with uQa, uQb and QaQb + // Sanity checks passed, calculate, print and store switch (poi) { case 1: // Subevent A reference flow { @@ -505,19 +661,19 @@ bool AliFlowAnalysisWithMSP::Calculate(double &vn, double &vnerror, const AliFlo vnerror=-1; return false; } - double vnB = TMath::Sqrt( uQa*QaQb / uQb ); // vn - double VvnB = QaQb*VuQa/(4*uQa*uQb) // Variance of vn + double vnA = TMath::Sqrt( uQa*QaQb / uQb ); // vn + double VvnA = QaQb*VuQa/(4*uQa*uQb) // Variance of vn + uQa*VQaQb/(4*QaQb*uQb) + uQa*QaQb*VuQb/(4*TMath::Power(uQb,3)) + CuQaQaQb/(2*uQb) - QaQb*CuQauQb/(2*TMath::Power(uQb,2)) - uQa*CuQbQaQb/(2*TMath::Power(uQb,2)); - vn=vnB; - if( VvnB<0 ) { - vnerror=VvnB; + vn=vnA; + if( VvnA<0 ) { + vnerror=VvnA; return false; } - vnerror=TMath::Sqrt(VvnB); + vnerror=TMath::Sqrt(VvnA); } break; case 2: // Subevent B reference flow @@ -527,19 +683,19 @@ bool AliFlowAnalysisWithMSP::Calculate(double &vn, double &vnerror, const AliFlo vnerror=-1; return false; } - double vnA = TMath::Sqrt( uQb*QaQb / uQa ); // vn - double VvnA = uQb*VQaQb/(4*QaQb*uQa) // Variance of vn + double vnB = TMath::Sqrt( uQb*QaQb / uQa ); // vn + double VvnB = uQb*VQaQb/(4*QaQb*uQa) // Variance of vn + QaQb*VuQb/(4*uQb*uQa) + QaQb*uQb*VuQa/(4*TMath::Power(uQa,3)) + CuQbQaQb/(2*uQa) - uQb*CuQaQaQb/(2*TMath::Power(uQa,2)) - - uQa*CuQauQb/(2*TMath::Power(uQa,2)); - vn=vnA; - if( VvnA<0 ) { - vnerror=VvnA; + - QaQb*CuQauQb/(2*TMath::Power(uQa,2)); + vn=vnB; + if( VvnB<0 ) { + vnerror=VvnB; return false; } - vnerror=TMath::Sqrt(VvnA); + vnerror=TMath::Sqrt(VvnB); } break; default: // POI flow @@ -566,7 +722,7 @@ bool AliFlowAnalysisWithMSP::Calculate(double &vn, double &vnerror, const AliFlo + CuQauQb/(2*QaQb) - uQb*CuQaQaQb/(2*TMath::Power(QaQb,2)) - uQa*CuQbQaQb/(2*TMath::Power(QaQb,2)); - vn=vnP; + vn=TMath::Sign(vnP,uQb); if( VvnP<0 ) { vnerror=VvnP; return false; @@ -578,7 +734,7 @@ bool AliFlowAnalysisWithMSP::Calculate(double &vn, double &vnerror, const AliFlo return (vnerror>=0); } -void AliFlowAnalysisWithMSP::ReadHistograms(TDirectoryFile *file) +void AliFlowAnalysisWithMSP::ReadHistograms(TDirectory *file) { delete fCommonHist; fCommonHist=0; // Delete existing histograms delete fQaComponents; fQaComponents=0; @@ -593,6 +749,11 @@ void AliFlowAnalysisWithMSP::ReadHistograms(TDirectoryFile *file) delete fDiffFlowPt; fDiffFlowPt=0; delete fDiffFlowEta; fDiffFlowEta=0; delete fFlags; fFlags=0; + if( fHistList ) { + fHistList->SetOwner(kFALSE); + delete fHistList; + fHistList=0; + } file->GetObject("QaComponents",fQaComponents); file->GetObject("QbComponents",fQbComponents); @@ -606,23 +767,86 @@ void AliFlowAnalysisWithMSP::ReadHistograms(TDirectoryFile *file) std::cerr << "AliFlowAnalysisWithMSP::ReadHistograms() : One or more histograms were not read correctly from " << file->GetPath() << std::endl; } + file->GetObject("Flags",fFlags); // Flags are required + + if( !fFlags ){ + std::cerr << "AliFlowAnalysisWithMSP::ReadHistograms(TDirectoryFile *) : Flags histogram not found, using defaults" << std::endl; + fHarmonic=2; + fNUA=false; + }else{ + fHarmonic=(UInt_t)(fFlags->GetBinContent(1)); + double harmonicSpread=fFlags->GetBinError(1); + if( harmonicSpread!=0 ) { + std::cerr << "AliFlowAnalysisWithMSP::ReadHistograms(TDirectoryFile *) :These histograms seem to be merged from analysis with different Harmonics. Results removed!" << std::endl; + delete fIntegratedFlow; fIntegratedFlow=0; + delete fDiffFlowPt; fDiffFlowPt=0; + delete fDiffFlowEta; fDiffFlowEta=0; + } + fNUA=fFlags->GetBinContent(2); // Mixing NUA does not matter since it needs to be recalculated anyway + double nuaSpread=fFlags->GetBinError(2); + if( nuaSpread!=0 ) { + std::cerr << "AliFlowAnalysisWithMSP::ReadHistograms(TDirectoryFile *) :These histograms seem to be merged from analysis with different NUA corrections. Results removed" << std::endl; + delete fIntegratedFlow; fIntegratedFlow=0; + delete fDiffFlowPt; fDiffFlowPt=0; + delete fDiffFlowEta; fDiffFlowEta=0; + } + } + // Optional histograms, may return a zero pointer - file->GetObject("AliFlowCommonHist_MSP",fCommonHist); // The AliFlowCommonHist is optional - file->GetObject("IntegratedFlow",fIntegratedFlow); + file->GetObject("AliFlowCommonHist_MSP",fCommonHist); // The AliFlowCommonHist is optional + file->GetObject("IntegratedFlow",fIntegratedFlow); // Results are optional file->GetObject("DiffFlowPt",fDiffFlowPt); file->GetObject("DiffFlowEta",fDiffFlowEta); - file->GetObject("Flags",fFlags); + fBookCommonHistograms=(fCommonHist!=0); +} + + +void AliFlowAnalysisWithMSP::ReadHistograms(TList *list) +{ + if( !list ) return; + delete fCommonHist; fCommonHist=0; // Delete existing histograms if any + delete fQaComponents; fQaComponents=0; + delete fQbComponents; fQbComponents=0; + delete fQaQb; fQaQb=0; + delete fPtUComponents; fPtUComponents=0; + delete fEtaUComponents; fEtaUComponents=0; + delete fAllStatistics; fAllStatistics=0; + delete fPtStatistics; fPtStatistics=0; + delete fEtaStatistics; fEtaStatistics=0; + delete fIntegratedFlow; fIntegratedFlow=0; + delete fDiffFlowPt; fDiffFlowPt=0; + delete fDiffFlowEta; fDiffFlowEta=0; + delete fFlags; fFlags=0; + if( fHistList ) { + fHistList->SetOwner(kFALSE); + delete fHistList; + fHistList=0; + } + + fQaComponents = static_cast(list->FindObject("QaComponents")); + fQbComponents = static_cast(list->FindObject("QbComponents")); + fQaQb = static_cast(list->FindObject("QaQb")); + fPtUComponents = static_cast(list->FindObject("PtUComponents")); + fEtaUComponents = static_cast(list->FindObject("EtaUComponents")); + fAllStatistics = static_cast(list->FindObject("AllStatistics")); + fPtStatistics = static_cast(list->FindObject("PtStatistics")); + fEtaStatistics = static_cast(list->FindObject("EtaStatistics")); + if( !fQaComponents || !fQbComponents || !fQaQb || !fPtUComponents || !fEtaUComponents || !fAllStatistics || !fPtStatistics || !fEtaStatistics ) { + std::cerr << "AliFlowAnalysisWithMSP::ReadHistograms(Tlist *) : One or more histograms were not read correctly from TList" << std::endl; + } + + fFlags = static_cast(list->FindObject("Flags")); // Flags are required if( !fFlags ){ - std::cerr << "AliFlowAnalysisWithMSP::ReadHistograms() : Flags histogram not found, using defaults" << std::endl; + std::cerr << "AliFlowAnalysisWithMSP::ReadHistograms(TList *) : Flags histogram not found, using defaults" << std::endl; fHarmonic=2; fNUA=false; }else{ - fHarmonic=fFlags->GetBinContent(1); + fHarmonic=(UInt_t)(fFlags->GetBinContent(1)); double harmonicSpread=fFlags->GetBinError(1); if( harmonicSpread!=0 ) { - std::cerr << "These histograms seem to be merged from analysis with different Harmonics. Results are invalid!" << std::endl; + std::cerr << "These histograms seem to be merged from analysis with different Harmonics. Results removed!" << std::endl; delete fIntegratedFlow; fIntegratedFlow=0; delete fDiffFlowPt; fDiffFlowPt=0; delete fDiffFlowEta; fDiffFlowEta=0; @@ -636,5 +860,12 @@ void AliFlowAnalysisWithMSP::ReadHistograms(TDirectoryFile *file) delete fDiffFlowEta; fDiffFlowEta=0; } } + + // Optional histograms, may return a zero pointer + fCommonHist = static_cast(list->FindObject("AliFlowCommonHist_MSP")); // The AliFlowCommonHist is optional + fIntegratedFlow = static_cast(list->FindObject("IntegratedFlow")); // Results are optional + fDiffFlowPt = static_cast(list->FindObject("DiffFlowPt")); + fDiffFlowEta = static_cast(list->FindObject("DiffFlowEta")); + fBookCommonHistograms=(fCommonHist!=0); -} \ No newline at end of file +} diff --git a/PWG/FLOW/Base/AliFlowAnalysisWithMSP.h b/PWG/FLOW/Base/AliFlowAnalysisWithMSP.h index 4d1a5247d4a..b6859914b6a 100644 --- a/PWG/FLOW/Base/AliFlowAnalysisWithMSP.h +++ b/PWG/FLOW/Base/AliFlowAnalysisWithMSP.h @@ -34,16 +34,18 @@ class AliFlowAnalysisWithMSP : public TNamed public: AliFlowAnalysisWithMSP(); // Default constructor for root I/O AliFlowAnalysisWithMSP(const unsigned int harmonic, const bool commonConst=true, const bool commonHist=false); // Construct and define params - AliFlowAnalysisWithMSP(TDirectoryFile *file); // Construct and read from file - // TODO: derive from TNamed and write a copy constructor so that the entire class can be stored instead of WriteHistograms + AliFlowAnalysisWithMSP(TDirectory *file); // Construct and read from file + AliFlowAnalysisWithMSP(TList *list); // Construct and load histograms from a TList. This object becomes owner. + AliFlowAnalysisWithMSP(const AliFlowAnalysisWithMSP &x); // Copy constructor ~AliFlowAnalysisWithMSP(); // Destructor assumes everything is owned void Init(); // (re)define all output objects, clears results void Make(AliFlowEventSimple *event); // Analyze one event void Finish(); // Calculate final results without redefining NUA usage - void Finish(bool nua){fNUA=nua;Finish();}; // Calculate while redefining if NUA is applied + void Finish(bool nua){fNUA=nua;Finish();}; // Calculate after redefining if NUA is applied void WriteHistograms(TDirectoryFile *file) const; // Write all histograms to a file void WriteCommonResults(TDirectoryFile *file) const; // Write an AliFlowCommonHistResults() to file + TList *ListHistograms(); // Create a TList from all histograms, ownership is kept in this class void Print(const Option_t *opt=0)const; // Summarize results on std::cout @@ -51,14 +53,15 @@ public: void UseCommonConstants(bool b=true){fUseCommonConstants=b;}; // Get parameters from AliFlowCommonConstants. MUST be called before Init() void EnableCommonHistograms(bool c=true){fBookCommonHistograms=c;}; // Create and fill AliFlowCommonHistograms() MUST be called before Init() void EnableNUA(bool c=true){fNUA=c;}; // Switch on/off NUA corrections. MUST be called before Finish() - // TODO: Allow NUA off, NUAu integrated, NUAu per bin instead of On/OFF. Default should be backward compatible - // TODO: Allow setting of pt and eta binning - // TODO: Allow phi-weights + // TODO: Allow NUA off, NUAu integrated, NUAu per bin instead of On(per bin)/OFF. Default should be backward compatible + // TODO: Allow setting of pt and eta binning in a simpler way, i.e. independent of AliFlowCommonConstants. const TH1D* GetIntegratedFlow(){return fIntegratedFlow;}; // Access result const TH1D* GetDiffFlowPt(){return fDiffFlowPt;}; // Access to result const TH1D* GetDiffFlowEta(){return fDiffFlowEta;}; // Access to result + AliFlowAnalysisWithMSP &operator=(const AliFlowAnalysisWithMSP &x); // Assignment + private: UInt_t fHarmonic; // Harmonic (n) in v_n analysis, used in Init() and Make() Bool_t fNUA; // Correct for Non Uniform Acceptance during Finish() @@ -67,6 +70,7 @@ private: AliFlowCommonHist *fCommonHist; // Standard control histograms, if enabled + // Correlations AliFlowMSPHistograms *fQaComponents; // Averages of Qa components per event for NUA AliFlowMSPHistograms *fQbComponents; // Averages of Qb components per event for NUA AliFlowMSPHistograms *fQaQb; // Average of QaQb per event @@ -76,16 +80,21 @@ private: AliFlowMSPHistograms *fPtStatistics; // Correlations for uQa uQb and QaQb per pt bin AliFlowMSPHistograms *fEtaStatistics; // Correlations for uQa uQb and QaQb per eta bin - TH1D *fIntegratedFlow; // vn for POI and subevents - TH1D *fDiffFlowPt; // vn as function of pt - TH1D *fDiffFlowEta; // vn as function of eta - TProfile *fFlags; // Stores fHarmonic and fNUA + // Result histograms (if calculated) + TH1D *fIntegratedFlow; // vn for POI and subevents + TH1D *fDiffFlowPt; // vn as function of pt + TH1D *fDiffFlowEta; // vn as function of eta + TProfile *fFlags; // Stores fHarmonic and fNUA + + TList *fHistList; // List of all histograms if requested by ListHistograms() + // Helper functions bool Calculate(double &vn, double &vnerror, const AliFlowMSPHistograms *hist, const AliFlowMSPHistograms *comp, const int bin=1, const int poi=0)const; // Calculates flow from histograms with NUA corrections bool Calculate(double &vn, double &vnerror, const AliFlowMSPHistograms *hist, const int bin=1, const int poi=0)const{return Calculate(vn,vnerror,hist,0,bin,poi);}; // No NUA version - void ReadHistograms(TDirectoryFile *file); // Restore histograms from a file + void ReadHistograms(TDirectory *file); // Restore histograms from a file + void ReadHistograms(TList *list); // Restore histograms from a TList - ClassDef(AliFlowAnalysisWithMSP,1); // Class version + ClassDef(AliFlowAnalysisWithMSP,0); // Class version }; #endif //ALIFLOWANALYSISWITHMSP_H diff --git a/PWG/FLOW/Base/AliFlowEventSimple.cxx b/PWG/FLOW/Base/AliFlowEventSimple.cxx index fc04f57d00e..32cb654be2f 100644 --- a/PWG/FLOW/Base/AliFlowEventSimple.cxx +++ b/PWG/FLOW/Base/AliFlowEventSimple.cxx @@ -195,6 +195,17 @@ void AliFlowEventSimple::AddTrack( AliFlowTrackSimple* track ) fNumberOfTracks++; } +//----------------------------------------------------------------------- +AliFlowTrackSimple* AliFlowEventSimple::MakeNewTrack() +{ + AliFlowTrackSimple *t=dynamic_cast(fTrackCollection->RemoveAt(fNumberOfTracks)); + if( !t ) { // If there was no track at the end of the list then create a new track + t=new AliFlowTrackSimple(); + } + + return t; +} + //----------------------------------------------------------------------- AliFlowVector AliFlowEventSimple::GetQ(Int_t n, TList *weightsList, Bool_t usePhiWeights, Bool_t usePtWeights, Bool_t useEtaWeights) { diff --git a/PWG/FLOW/Base/AliFlowEventSimple.h b/PWG/FLOW/Base/AliFlowEventSimple.h index dfc3e9df2ee..c19c37b2013 100644 --- a/PWG/FLOW/Base/AliFlowEventSimple.h +++ b/PWG/FLOW/Base/AliFlowEventSimple.h @@ -86,6 +86,7 @@ class AliFlowEventSimple: public TObject { AliFlowTrackSimple* GetTrack(Int_t i); void AddTrack( AliFlowTrackSimple* track ); + AliFlowTrackSimple* MakeNewTrack(); AliFlowVector GetQ(Int_t n=2, TList *weightsList=NULL, Bool_t usePhiWeights=kFALSE, Bool_t usePtWeights=kFALSE, Bool_t useEtaWeights=kFALSE); void Get2Qsub(AliFlowVector* Qarray, Int_t n=2, TList *weightsList=NULL, Bool_t usePhiWeights=kFALSE, Bool_t usePtWeights=kFALSE, Bool_t useEtaWeights=kFALSE); diff --git a/PWG/FLOW/Base/AliFlowMSPHistograms.cxx b/PWG/FLOW/Base/AliFlowMSPHistograms.cxx index eb4f2523b4c..e600ada8a2d 100644 --- a/PWG/FLOW/Base/AliFlowMSPHistograms.cxx +++ b/PWG/FLOW/Base/AliFlowMSPHistograms.cxx @@ -30,8 +30,24 @@ AliFlowMSPHistograms::AliFlowMSPHistograms(const AliFlowMSPHistograms &x) : TNamed(),fVarNames(),fPAvg(x.fPAvg),fPProd(x.fPProd),fXName(x.fXName) { // Copy constructor, clone all the TProfiles from the other object + // The copy constructor made a shallow copy, i.e. all entries point to the originial profiles. + // Therefore we need to Clone each entry. + for( int i=0; iClone(); + } fPAvg.SetOwner(kTRUE); + + for( int i=0; iClone(); + } fPProd.SetOwner(kTRUE); } @@ -42,6 +58,7 @@ AliFlowMSPHistograms::AliFlowMSPHistograms(const int dim, const char *name, cons } AliFlowMSPHistograms::AliFlowMSPHistograms(const int dim, const char *name, const int nBins, const double *xBins) + : TNamed(name,"MSP histograms"),fVarNames(),fPAvg(),fPProd(),fXName() { Init(dim, name, nBins, xBins); } @@ -264,6 +281,40 @@ AliFlowMSPHistograms &AliFlowMSPHistograms::operator+=(const AliFlowMSPHistogram return *this; } +AliFlowMSPHistograms &AliFlowMSPHistograms::operator=(const AliFlowMSPHistograms &x) +{ + // Assignment, clone all the TProfiles from the other object + // We need to Clone each entry and delete existing ones. + for( int i=0; iClone(); + } + fPAvg.SetOwner(kTRUE); + + for( int i=0; iClone(); + } + fPProd.SetOwner(kTRUE); + return *this; +} + +TObject *AliFlowMSPHistograms::Clone(const char *n) const +{ + AliFlowMSPHistograms *x=new AliFlowMSPHistograms(*this); + if(n)x->SetName(n); + return x; +} + // Private functions................................................................................... double AliFlowMSPHistograms::Covariance(double x, double y, double xy, double wx, double wy, double wxy)const diff --git a/PWG/FLOW/Base/AliFlowMSPHistograms.h b/PWG/FLOW/Base/AliFlowMSPHistograms.h index 87db365131d..0bd722aefc1 100644 --- a/PWG/FLOW/Base/AliFlowMSPHistograms.h +++ b/PWG/FLOW/Base/AliFlowMSPHistograms.h @@ -47,6 +47,9 @@ public: Long64_t Merge(TCollection *); // Merge all objects in the collection into this AliFlowMSPHistograms &operator+=(const AliFlowMSPHistograms &x); // Merge this with another AliFlowMSPHistograms object + AliFlowMSPHistograms &operator=(const AliFlowMSPHistograms &x); + + TObject *Clone(const char *n=0)const; // Create a copy of this object private: TObjArray fVarNames; // Names of variables being correlated diff --git a/PWG/FLOW/Base/AliFlowVector.cxx b/PWG/FLOW/Base/AliFlowVector.cxx index f4b15835d87..7a041093580 100644 --- a/PWG/FLOW/Base/AliFlowVector.cxx +++ b/PWG/FLOW/Base/AliFlowVector.cxx @@ -26,9 +26,9 @@ ClassImp(AliFlowVector) //________________________________________________________________________ AliFlowVector::AliFlowVector(): - fMult(0) + TVector2(0.0,0.0),fMult(0.0) { - // default contructor + // default constructor } //________________________________________________________________________ @@ -40,24 +40,43 @@ AliFlowVector::AliFlowVector(const AliFlowVector& aVector): // copy constructor } +//________________________________________________________________________ + +AliFlowVector::AliFlowVector(Double_t *y, const Double_t m): + TVector2(y), + fMult(m) +{ + // Analogue of TVector2 constructor. Sets (x,y) and multiplicity 1. +} + //________________________________________________________________________ AliFlowVector::AliFlowVector(const TVector2 &v, const Double_t m): TVector2(v), fMult(m) { - // custom constructor + // custom constructor, Sets vector and multiplicity +} + + //________________________________________________________________________ + +AliFlowVector::AliFlowVector(const Double_t x, const Double_t y, const Double_t m): + TVector2(x,y), + fMult(m) +{ + // custom constructor analogue of TVector2 constructor } //________________________________________________________________________ AliFlowVector::~AliFlowVector() { - // default constructor + // default destructor } void AliFlowVector::SetMagPhi(double size, double angle, double mult) { + // Analogue to SetMagPhi for a TVector2 but here including a sum of weights TVector2::SetMagPhi(size,angle); SetMult(mult); } @@ -71,7 +90,6 @@ AliFlowVector& AliFlowVector::operator=(const AliFlowVector& aVector) fX = aVector.X(); fY = aVector.Y(); fMult = aVector.GetMult(); - return *this; } @@ -82,17 +100,24 @@ AliFlowVector& AliFlowVector::operator+=(const AliFlowVector& aVector) // addition operator fX += aVector.X(); fY += aVector.Y(); - fMult += aVector.GetMult(); - + fMult += aVector.GetMult(); return *this; } AliFlowVector& AliFlowVector::operator-=(const AliFlowVector& aVector) { - // addition operator + // subtraction operator fX -= aVector.X(); fY -= aVector.Y(); fMult -= aVector.GetMult(); - return *this; } + +AliFlowVector& AliFlowVector::operator*=(const double w) +{ + // multiply by a weight operator + fX*=w; + fY*=w; + fMult*=w; + return *this; +} diff --git a/PWG/FLOW/Base/AliFlowVector.h b/PWG/FLOW/Base/AliFlowVector.h index 929f643bce4..4a64e815164 100644 --- a/PWG/FLOW/Base/AliFlowVector.h +++ b/PWG/FLOW/Base/AliFlowVector.h @@ -17,19 +17,27 @@ class AliFlowVector: public TVector2 { public: AliFlowVector(); AliFlowVector(const AliFlowVector& aVector); - AliFlowVector(const TVector2 &p, const Double_t m); + AliFlowVector(const TVector2 &p, const Double_t m); // constructor: Add a weight to the TVector + AliFlowVector(Double_t *y, const Double_t m=1); // constructor: Analogue to TVector2(y) with multiplicity + AliFlowVector(const Double_t x, const Double_t y, const Double_t m=1); // constructor: Sets the components individually virtual ~AliFlowVector(); - void SetMagPhi(double size, double angle, double mult=1); // Set vector and weighted multiplicity + void SetMagPhi(Double_t size, Double_t angle, Double_t mult=1); // Set vector and weighted multiplicity - AliFlowVector& operator=(const AliFlowVector& aVector); // Assign - AliFlowVector& operator+=(const AliFlowVector& aVector); // Add to self - AliFlowVector& operator-=(const AliFlowVector& aVector); // Subtract from self + AliFlowVector& operator=(const AliFlowVector& aVector); // Assign to self + AliFlowVector& operator+=(const AliFlowVector& aVector); // Add to self + AliFlowVector& operator-=(const AliFlowVector& aVector); // Subtract from self + AliFlowVector& operator*=(const double w); // Multiply by a weight + AliFlowVector& operator/=(const double w){ (*this)*=(1.0/w); return *this;}; // Divide by a weight + const AliFlowVector operator+(const AliFlowVector&a) const { AliFlowVector v(*this); return v+=a; }; // Add and return by value + const AliFlowVector operator-(const AliFlowVector&a) const { AliFlowVector v(*this); return v-=a; }; // Subtract and return by value + const AliFlowVector operator*(const double w) const { AliFlowVector v(*this); return v*=w; }; // Scale and return by value + const AliFlowVector operator/(const double w) const { AliFlowVector v(*this); return v/=w; }; // Scale and return by value Bool_t IsFolder() const {return kTRUE;}; - void SetMult(Double_t const mult) {this->fMult = mult;}; - Double_t GetMult() const {return this->fMult;}; + void SetMult(const Double_t mult) {fMult = mult;}; // Set sum of weights + Double_t GetMult() const {return fMult;}; // Get sum of weights private: Double_t fMult; // multiplicity = sum of weights = w_1 + w_2 + ... + w_n @@ -37,6 +45,8 @@ class AliFlowVector: public TVector2 { ClassDef(AliFlowVector, 1) }; + +/* Old, less efficient code inline AliFlowVector operator+(const AliFlowVector& aVector,const AliFlowVector& bVector) { AliFlowVector cVector; Double_t x = aVector.X() + bVector.X(); @@ -60,5 +70,6 @@ inline AliFlowVector operator-(const AliFlowVector& aVector,const AliFlowVector return cVector; } +*/ #endif diff --git a/PWG/FLOW/Tasks/AliAnalysisTwoParticleResonanceFlowTask.h b/PWG/FLOW/Tasks/AliAnalysisTwoParticleResonanceFlowTask.h index b8b9673c6cd..a31bcb449ad 100644 --- a/PWG/FLOW/Tasks/AliAnalysisTwoParticleResonanceFlowTask.h +++ b/PWG/FLOW/Tasks/AliAnalysisTwoParticleResonanceFlowTask.h @@ -74,7 +74,7 @@ public: TH1F* InitPtSpectraHistograms(Float_t nmin, Float_t nmax); TH1F* BookPtHistogram(const char* name); Bool_t InitializeAnalysis(); - void AddResonanceIdentificationOutputObjects(); + //PK void AddResonanceIdentificationOutputObjects(); virtual void UserCreateOutputObjects(); // setters void SetPtBins(Float_t bin[19], Int_t n) { for(Int_t i = 0; i < n+1; i++) fPtBins[i] = bin[i]; fNPtBins = n; } diff --git a/PWG/PWGflowBaseLinkDef.h b/PWG/PWGflowBaseLinkDef.h index 6c913b53298..6aaacfe0de1 100644 --- a/PWG/PWGflowBaseLinkDef.h +++ b/PWG/PWGflowBaseLinkDef.h @@ -26,11 +26,13 @@ #pragma link C++ class AliFlowCommonHistResults+; #pragma link C++ class AliFlowLYZHist1+; #pragma link C++ class AliFlowLYZHist2+; +#pragma link C++ class AliFlowMSPHistograms+; #pragma link C++ class AliFlowLYZEventPlane+; #pragma link C++ class AliFlowAnalysisWithMCEventPlane+; #pragma link C++ class AliFlowAnalysisWithScalarProduct+; +#pragma link C++ class AliFlowAnalysisWithMSP+; #pragma link C++ class AliFlowAnalysisWithLYZEventPlane+; #pragma link C++ class AliFlowAnalysisWithLeeYangZeros+; #pragma link C++ class AliFlowAnalysisWithCumulants+; diff --git a/PWG/PWGflowTasksLinkDef.h b/PWG/PWGflowTasksLinkDef.h index e5c7de23ca0..9ae59c453e9 100644 --- a/PWG/PWGflowTasksLinkDef.h +++ b/PWG/PWGflowTasksLinkDef.h @@ -12,6 +12,7 @@ #pragma link C++ class AliFlowEventSimpleMaker+; #pragma link C++ class AliAnalysisTaskScalarProduct+; +#pragma link C++ class AliAnalysisTaskMSP+; #pragma link C++ class AliAnalysisTaskMCEventPlane+; #pragma link C++ class AliAnalysisTaskLYZEventPlane+; #pragma link C++ class AliAnalysisTaskLeeYangZeros+; diff --git a/PWG/muon/AliDimuInfoStoreMC.h b/PWG/muon/AliDimuInfoStoreMC.h index 34fb6ba1985..afa86723657 100644 --- a/PWG/muon/AliDimuInfoStoreMC.h +++ b/PWG/muon/AliDimuInfoStoreMC.h @@ -38,7 +38,7 @@ class AliDimuInfoStoreMC : public AliDimuInfoStoreRD { private: void FindDimuonSourceFast(); - void FindDimuonSourceFull(); + void FindDimuonSourceFull(){}; // PK, not implemented in source yet static const TString fgkStdBranchName; // Standard branch name static const Int_t fgkSourcesN; // num. of dimuon sources -- 2.43.0